Community Sampling


In [87]:
from collections import namedtuple
from itertools import cycle

from os.path import join, exists
from os import makedirs, chdir,getcwd,fchmod
from glob import glob
import shutil

import time

import numpy as np

import pandas as pd

import re

import scipy as sc
from scipy.stats import norm
from scipy.stats import expon
import scipy.interpolate as interpolate
import scipy.integrate
from scipy.fftpack import fft, ifft

from lmfit.models import Model
from lmfit import conf_interval

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

from Bio.Seq import MutableSeq
from Bio.Alphabet import IUPAC

import docker

import configparser

import seaborn as sns
sns.set_context("paper")

Definitions


In [44]:
def create_mount_run(image,mount_dir,cmd,envs):   
    container = cli.create_container(
        image=image, command=cmd, volumes=['/mnt/vol'],
        host_config=cli.create_host_config(binds={
            mount_dir: {
                'bind': '/mnt/vol',
                'mode': 'rw',
            }
        }),
        environment=envs
    )
    ctr=container.get('Id')
    cli.start(ctr)
    cli.wait(ctr,60*60*24*10)
    return cli.logs(ctr)

In [45]:
def Gekv(B,C):
    return 2**(1 - B)/(C*np.log(2))*(1 - 2**(-C))

In [46]:
def Gekv(B,C):
    return 2**(1 - B)/(C*np.log(2))*(1 - 2**(-C))

In [47]:
class Community:
    "A community defined by genome references (Biopython SeqRecords) and corresponding growth parameters."
    def __init__(self,name,acc,growth_param,td,mapper,env):
        self.name = name
        self.conf = local_conf(join(td,name),mapper)
        self.d_conf = local_conf(join('/mnt/vol',name),mapper)
        self.args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join(self.d_conf['node_path'],'project.conf'),'-e','hedani@chalmers.se']
        self.env=env
        
        self.create_dirs(True)
        save_config(self.conf,self.d_conf)
        self.fetch_ref(acc)
        
        self.pop = self.init_population(acc,growth_param)
        self.distribution = self.community_distribution()
        self.samples = []
        
    def init_population(self,acc,growth_param):
        records=open_records(glob(join(self.conf['ref_path'],'Fasta','*.fasta')))
        population = namedtuple("population", "B C D l seq cells")
        
        #refId=[x.id for x in records]
        #acc=keys(comm)
        #refInd=[acc.index(x) for x in refId]
        
        pop = {}
        
        #for i,rec in enumerate(records):
        # add [:-2] to a if . removed 
        for i,a in enumerate(acc): 
            pop[a]=population(B = growth_param[i][0],C = growth_param[i][1],
                             D = growth_param[i][2], l = len(records[a]),
                             seq = records[a], cells=growth_param[i][3])
        return pop
    
    def ptr(self):
        for i,a in enumerate(self.pop.keys()):
            print a+": "+str(2**growth_param[i][1])
    
    def community_distribution(self):
        d=np.array([Gekv(p.B,p.C)*p.l*p.cells for p in self.pop.values()])
        return d/d.sum()

    #def ab_(self):   
    #    return 
    
    def sample(self,nr_samples):
        nr_samples=np.array(nr_samples*self.distribution)
        nr_samples=nr_samples.astype(np.int)
        
        samp_run=[]
        for i,p in enumerate(self.pop.keys()):
            samp_run.append(inverse_transform_sampling(self.pop[p].C,nr_samples[i],self.pop[p].l))
        
        self.samples.append(samp_run)
    
    def write_reads(self):
        for i,samp in enumerate(self.samples):
            write_reads(samp,self,self.conf['data_path'],self.name+str(i))
    
    def compare_fit(self):
        if not self.samples:
            print "The community is not sampled, please run community.sample(nr_samples)"
            return;
        
        err_hfit=[]
        err_pfit=[]
        res_fit=[]
        
        for i,samp in enumerate(self.samples):
            for acc in self.pop.keys():
                try:
                    depth_file=join(self.conf['output_path'],self.name+str(i),'npy',acc+'.depth.npy')
                    best_file=join(self.conf['output_path'],self.name+str(i),'npy',acc+'.depth.best.npy')

                    signal=2**(np.load(depth_file))
                    signal=signal/(signal.sum()/len(signal))#*self.pop[acc].l)
                    from_ptr=np.load(best_file)

                    res=fit_signal(signal,self.pop[acc].l)
                    res_fit.append(res)
                    
                    err_hfit.append((self.pop[acc].C-res.best_values['C'])/self.pop[acc].C)

                    err_pfit.append((self.pop[acc].C-from_ptr[2]+from_ptr[3])/self.pop[acc].C)
                    print "Simulated value: "+str(self.pop[acc].C)
                    print "Error from this fit: "+str(err_hfit[-1])+ ' value: ' + str(res.best_values['C'])
                    print "Error from initial PTR fit "+str(err_pfit[-1])+' value: ' + str(from_ptr[2]-from_ptr[3])
                except Exception as Ex:
                    print "Ex"
                    pass
        return [res_fit,err_hfit,err_pfit]
    
    
    def fetch_ref(self,acc=''):
        acc_path = join(self.conf['node_path'],"acc")
        f = open(acc_path, "w")
        
        if not acc:
            f.write("\n".join(self.pop.keys()))
        else:
            f.write("\n".join(acc))
        f.close()
        
        create_mount_run('centos/hedani',td,self.args0+['fetch-references','-s',join(self.d_conf['node_path'],'acc')],self.env)

    def build_index(self):
        create_mount_run('centos/hedani',td,self.args0+['build-index'],self.env)

    def run_pipeline(self):
        create_mount_run('centos/hedani',td," ".join(['/bin/bash -c "']+self.args0+['make']+[';']+self.args0+['run"']),self.env)
        
    def collect(self):
        create_mount_run('centos/hedani',td,self.args0+['collect'],self.env)
        
    def create_dirs(self,from_init=False):
        if exists(self.conf['node_path']) and from_init:
            shutil.rmtree(self.conf['node_path'])
        
        for d in [self.conf['node_path'],self.conf['data_path'],self.conf['output_path'],self.conf['ref_path']]:
            if not exists(d):
                makedirs(d)

    #def complexity(self):
        #len(self.pop)
        #return complexity

In [48]:
import numpy as np

def cartesian(arrays, out=None):
    """
    Generate a cartesian product of input arrays.

    Parameters
    ----------
    arrays : list of array-like
        1-D arrays to form the cartesian product of.
    out : ndarray
        Array to place the cartesian product in.

    Returns
    -------
    out : ndarray
        2-D array of shape (M, len(arrays)) containing cartesian products
        formed of input arrays.

    Examples
    --------
    >>> cartesian(([1, 2, 3], [4, 5], [6, 7]))
    array([[1, 4, 6],
           [1, 4, 7],
           [1, 5, 6],
           [1, 5, 7],
           [2, 4, 6],
           [2, 4, 7],
           [2, 5, 6],
           [2, 5, 7],
           [3, 4, 6],
           [3, 4, 7],
           [3, 5, 6],
           [3, 5, 7]])

    """

    arrays = [np.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = np.prod([x.size for x in arrays])
    if out is None:
        out = np.zeros([n, len(arrays)], dtype=dtype)

    m = n / arrays[0].size
    out[:,0] = np.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian(arrays[1:], out=out[0:m,1:])
        for j in xrange(1, arrays[0].size):
            out[j*m:(j+1)*m,1:] = out[0:m,1:]
    return out

In [49]:
def pxn(x, C, l):
    return (2**(1 + C - (2*C*x)/float(l))*C*np.log(2))/(float(l)*(-1 + 2**C));

In [50]:
def pxn2(x, C):
    return (2**(1 + C - (2*C*x))*C*np.log(2))/((-1 + 2**C));

In [51]:
def piecewise_prob(x,C,l):
    conds=[(0 < x) & (x <= float(l)/2), x > float(l)/2]
    funcs=[lambda x: pxn(x,C,l)/float(2),lambda x: pxn((l-x),C,l)/float(2)]
    return np.piecewise(x,conds,funcs)

In [52]:
def piecewise_prob2(x,C):
    conds=[(0 < x) & (x <= float(1)/2), x > float(1)/2]
    funcs=[lambda x: pxn(x,C,1)/float(2),lambda x: pxn2((1-x),C,1)/float(2)]
    return np.piecewise(x,conds,funcs)

In [53]:
def inverse_transform_sampling(C,n_samples,l=1):
    x = np.linspace(0, float(1)/1.7, 100)
    y = np.linspace(pxn(.5,C,1),pxn(0,C,1),100)

    pdf=interpolate.interp1d(x,pxn(x,C,1))
    cdf = [scipy.integrate.quad(lambda x: pdf(x),0,i)[0] for i in np.linspace(0,.5,100)]
    inv_cdf = interpolate.interp1d(cdf, np.linspace(0,.5,100))
    
    r = np.random.rand(n_samples)
    v=l*np.round(inv_cdf(r[:int(n_samples/2)]),4)
    v2=l*(1-np.round(inv_cdf(r[int(n_samples/2):]),4))
    v=np.concatenate([v,v2])
    
    return v.astype(np.int)

In [54]:
def read_pair(header,ind,seq,pos,l):
    def circular_yield(x,pos,sub,l):
        if (pos>l):
            return x[pos-l:pos-l+sub]
        elif (pos+sub>l):
            r=pos+sub-l
            x2=x[pos:pos+sub-r]
            return x2+x[0:r]
        else:
            return x[pos:pos+sub]

    #base_error_rate = .02
    #mutation_rate = .001
    #nuc=["A","C","G","T"]
    
    # generate a normally distributued insert size of mean 300 and sd 40 
    #insert_size = int(norm.rvs(size=1,loc=300,scale=30)[0])
    insert_size = int(500)
    
    r1 = circular_yield(seq,pos,100,l)
    r2 = circular_yield(seq,pos+insert_size,100,l)
    
    # flip reads according to seqs error and mutation rate
    #ra=[]
    #for r in [r1,r2]:
    #    rr=np.random.random(100)
    #    m_ind=[i for i,j in enumerate(rr) if j < base_error_rate]
    #    char=np.random.choice(nuc,len(m_ind))
    #    a=list(r)
    #    for i,ind in enumerate(m_ind):
    #        a[ind]=char[i] 
    #    ra.append("".join(a))
        #r_tmp=r_tmp[:m_ind[0]]
        #for i,index in enumerate(m_ind[1:]):
        #    r_tmp = r_tmp[:index] + char[i] + r[index + 1:]
    #[r1,r2]=ra

    #nrs=[str(np.random.randint(low=1000,high=2000)),str(np.random.randint(low=10**4,high=2*10**4)),str(np.random.randint(low=10**5,high=2*10**5))]
    rec1=SeqRecord(r1,id=header+"_"+str(pos)+"."+str(ind),description="RANDXXLAB"+str(pos)+"_"+str(pos+insert_size+100)+":0:0"+"/1")
    rec2=SeqRecord(r2,id=header+"_"+str(pos)+"."+str(ind),description="RANDXXLAB"+str(pos)+"_"+str(pos+insert_size+100)+":0:0"+"/2")
    #rec1=SeqRecord(r1,id=header+"."+str(ind),description="FCB00YLABXX:6:"+nrs[0]+":"+nrs[1]+":"+nrs[2]+"/1")
    #rec2=SeqRecord(r2,id=header+"."+str(ind),description="FCB00YLABXX:6:"+nrs[0]+":"+nrs[1]+":"+nrs[2]+"/2")
    
    #rec1=SeqRecord(r1,id=header+"_"+str(pos)+"_"+str(pos+insert_size)+"/1",description="")
    #rec2=SeqRecord(r2,id=header+"_"+str(pos)+"_"+str(pos+insert_size)+"/2",description="")
    rec2=rec2.reverse_complement(id=header+"_"+str(pos)+"."+str(ind),description="RANDXXLAB"+str(pos)+"_"+str(pos+insert_size+100)+":0:0"+"/2")
    #rec2=SeqRecord(r2,id=header+"."+str(ind)+"/2",description="")
    rec1.letter_annotations['phred_quality']=[17]*100
    rec2.letter_annotations['phred_quality']=[17]*100
    
    #rec1.description=
    
    return(rec1,rec2)

In [55]:
def open_records(fasta):
    records={};
    for f in fasta:
        handle = open(f, "rU")
        tmp=list(SeqIO.parse(handle, "fasta"))
        m=re.match('.*(N[CTZ]_([A-Z]{2})*[0-9]{6}).*',tmp[0].id)
        records[m.group(0)]=tmp[0].seq
        handle.close()
    
    return records

In [56]:
def write_reads(samples,comm,directory,name):
    f1 = open(join(directory,name+"_1.fastq"), "w")
    f2 = open(join(directory,name+"_2.fastq"), "w")
    
    for i,p in enumerate(comm.pop.keys()):
        for j,pos in enumerate( samples[i].tolist()):
            r1,r2 = read_pair(name+"_"+p,str(j+1),comm.pop[p].seq,pos,comm.pop[p].l)

            SeqIO.write(r1, f1, "fastq")
            SeqIO.write(r2, f2, "fastq")
    
    f1.close()
    f2.close()

In [57]:
def plot_samples(samples):
    y=[np.bincount(s) for s in samples]
    for p in y:
        plt.plot(p/float(p.sum()),marker='.',linestyle='None')

In [58]:
def local_conf(td,mapper):   
    return {
            'project': 'ptr_simulation',
            'cluster': '',
            'job_name': 'ptr_simulation',
            'job_nodes': '1',
            'cpu_cores': '4',
            'estimated_time': '',

            'node_path': td,
            'ref_path': join(td,'References'),
            'data_path': join(td,'Data'),
            'output_path': join(td,'Out'),
            'doric_path': join(td,'DoriC'),

            'mapper': mapper,
            'ref_name': 'sim',
            'nr_samples': '1',
            'samples_per_node': '1',
            'email': 'hedani@chalmers.se',
            'data_prefix': '',
            'start_ind': '1',
            'job_range': '1',

            'ftp_url': '',
            'data_url': ''
        }

In [59]:
def run_community_analysis(exec_string,td,mapper,acc,growth_param,community,nr_samples,topname):
    chdir(td)
    comm=[]
    docker_td='/mnt/vol'
    
    for i,c in enumerate(community):
        accs=[acc[x] for x in c]
        grs=[growth_param[x] for x in c]

        wd='comm_'+str(topname)+str(i)
        docker_fd=join(docker_td,wd)
        fd=join(td,wd)
        
        # create dirs
        if 'd' in exec_string:
            if exists(fd):
                shutil.rmtree(fd)
            for d in [wd,join(wd,'Data'),join(wd,'Out')]:
                if not exists(d):
                    makedirs(d)
        
        chdir(wd)
        lconf=local_conf(fd,mapper)
        conf=local_conf(docker_fd,mapper)
        save_config(lconf,conf)
        
        args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join(conf['node_path'],'project.conf'),'-e','hedani@chalmers.se']
        
        if 'f' in exec_string:
            f = open("acc", "w")
            f.write("\n".join(accs))
            f.close()

            create_mount_run('centos/hedani',td,args0+['fetch-references','-s',join(conf['node_path'],'acc')])

        if 's' in exec_string:
            # create community and sample
            # name,acc,growth_param,td,mapper
            comm.append(Community(wd,accs,grs,td,'bowtie2'))
            comm[i].sample(nr_samples)
            comm[i].write_reads()

        if 'b' in exec_string: create_mount_run('centos/hedani',conf['node_path'],args0+['build-index'])

        if 'm' in exec_string: create_mount_run('centos/hedani',conf['node_path'],args0+['make'])

        if 'r' in exec_string: create_mount_run('centos/hedani',conf['node_path'],args0+['run'])

        #chdir('..')
        
    return comm

In [60]:
def fit_signal(signal,l):
    x1=np.linspace(0,1,len(signal))
    
    piecewiseModel=Model(piecewise_prob)
    piecewiseModel.set_param_hint('l',     value=1,vary=False)
    piecewiseModel.set_param_hint('C',   vary=True,  value=.1,min=0,max=1)
    piecewiseModel.make_params()
    
    res=piecewiseModel.fit(signal,x=x1,weights=np.sin(1.5*x1)+1.5)
    return res

In [61]:
def save_config(lconf,conf):
    Config = configparser.ConfigParser()
    Config.optionxform = str
    
    cfgfile = open(join(lconf['node_path'],'project.conf'),'w')

    Config.add_section('Project')
    Config.set('Project','ProjectID',conf['project'])
    Config.set('Project','Cluster',conf['cluster'])
    Config.set('Project','JobName',conf['job_name'])
    Config.set('Project','JobNodes',conf['job_nodes'])
    Config.set('Project','CpuCores',conf['cpu_cores'])
    Config.set('Project','EstimatedTime',conf['estimated_time'])
    
    Config.add_section('Directories')
    Config.set('Directories','Node',conf['node_path'])
    Config.set('Directories','References',conf['ref_path'])
    Config.set('Directories','Data',conf['data_path'])
    Config.set('Directories','Output',conf['output_path'])
    Config.set('Directories','DoriC',conf['doric_path'])
    
    Config.add_section('Other')
    Config.set('Other','Mapper',conf['mapper'])
    Config.set('Other','RefName',conf['ref_name'])
    Config.set('Other','NrSamples',conf['nr_samples'])
    Config.set('Other','SamplesPerNode',conf['samples_per_node'])
    Config.set('Other','Email',conf['email'])
    Config.set('Other','DataPrefix',conf['data_prefix'])
    Config.set('Other','StartInd',conf['start_ind'])
    Config.set('Other','JobRange',conf['job_range'])
    Config.set('Other','FtpURL',conf['ftp_url'])
    Config.set('Other','DataURL',conf['data_url'])
    
    Config.write(cfgfile)
    cfgfile.close()

In [62]:
def mutation_run():
    rec.seq=mutateSeq(mutable_seq,.01)
    
    MutableSeq(str(rec.seq), IUPAC.unambiguous_dna)
    
    rec.id="mut_"+rec.id
    rec.description="mut_"+rec.description
    rec.name="mut_"+rec.name
    with open(sys.argv[2], "w") as output_handle:
        SeqIO.write(rec, output_handle, "fasta")

In [75]:
def compare_fit(c):
    err_hfit=[]
    err_pfit=[]
    res_fit=[]

    for i,samp in enumerate(c.samples):
        res_fit_tmp=[]
        err_hfit_tmp=[]
        err_pfit_tmp=[]
        
        for acc in c.pop.keys():
            try:
                depth_file=join(c.conf['output_path'],c.name+str(i),'npy',acc+'.depth.npy')
                best_file=join(c.conf['output_path'],c.name+str(i),'npy',acc+'.depth.best.npy')

                signal=2**(np.load(depth_file))
                signal=signal/(signal.sum()/len(signal))#*self.pop[acc].l)
                
                # extra deblurring
                signal=smooth_signal(signal,10**3)
                signal=sharpen_filter(signal)
                
                from_ptr=np.load(best_file)

                res=fit_signal(signal,c.pop[acc].l)
                res_fit_tmp.append(res.best_values['C'])

                err_hfit_tmp.append((c.pop[acc].C-res.best_values['C'])/c.pop[acc].C)

                err_pfit_tmp.append((c.pop[acc].C-from_ptr[2]+from_ptr[3])/c.pop[acc].C)
                #print "Simulated value: "+str(c.pop[acc].C)
                #print "Error from this fit: "+str(err_hfit_tmp[-1])+ ' value: ' + str(res.best_values['C'])
                #print "Error from initial PTR fit "+str(err_pfit_tmp[-1])+' value: ' + str(from_ptr[2]-from_ptr[3])
            except Exception as Ex:
                print Ex
                pass
        
        res_fit.append(res_fit_tmp)
        err_hfit.append(err_hfit_tmp)
        err_pfit.append(err_pfit_tmp)
    return [res_fit,np.abs(err_hfit),np.abs(err_pfit)]

In [64]:
def compare_fit_ptrc(c):
    ptrc=[]
    ptre=[]
    for i,samp in enumerate(c.samples):
        ptrc_tmp=[]
        ptre_tmp=[]
        cmd=['/bin/bash -c "']+['cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh']+[c.name+str(i)+'_1.fastq']+[c.name+str(i)+'_2.fastq"']
        create_mount_run('centos/hedani',td," ".join(cmd),envs)
        ptr_file=join(c.conf['output_path'],'ptrc_out')
        ptrc_frame=pd.read_pickle(ptr_file)
        ind_val=ptrc_frame.index.values
        for i in ind_val:#range(ptrc_frame.shape[0]):
            ptrc_tmp.append(ptrc_frame.loc[i]['sm'])#.ix[i]['sm'])
            ptre_tmp.append(np.abs(np.log2(ptrc_frame.loc[i]['sm'])-c.pop[i[:-1]].C)/c.pop[i[:-1]].C)
            
        ptrc.append(ptrc_tmp)
        ptre.append(ptre_tmp)
       
    return [np.log2(ptrc),ptre]

In [65]:
def mv_filt(L,omega):
    return (1/float(L))*(1-np.exp(-omega*1j*L))/(1-np.exp(-omega*1j))

In [66]:
def smooth_signal(signal,length):
    xi=np.linspace(0,1,len(signal))
    xi2=np.linspace(0,1,length)
    iss=scipy.interpolate.UnivariateSpline(xi,signal,k=3)
    iss.set_smoothing_factor(.00001)
    return iss(xi2)

In [67]:
def sharpen_filter(signal):
    x = np.linspace(-0.00001*np.pi/float(1), 1*np.pi/float(1),10**3)
    return np.roll(ifft(fft(signal,10**3)/mv_filt(1.27,x))[0:len(signal)],np.int64(-125/2))

In [68]:
def triangle(length, amplitude, y_offset):
    section = length // 4
    for direction in (1, -1):
        for i in range(section):
            yield y_offset + i * (amplitude / section) * direction
        for i in range(section):
            yield y_offset + (amplitude - (i * (amplitude / section))) * direction

In [ ]:


In [ ]:
df = pd.DataFrame([[1, 2], [3, 4]], columns=list('AB'))
df2 = pd.DataFrame([[5, 6], [7, 8]], columns=list('AB'))
df.append(df2)

In [ ]:
real_c=np.tile([growth_param[sel[0][1]][1],growth_param[sel[0][0]][1]],(5,1))
ptrc_e=np.abs(np.log2(ptrc)-real_c)/real_c

In [31]:
create_mount_run('centos/hedani',td," ".join(cmd),envs)


Out[31]:
'2016-12-11 20:40:13,175-ERROR: The following ids have been listed in the id_species files but not sequence was found: ti|821|org|NC_009614.1\nTraceback (most recent call last):\n  File "/home/Programs/PTRC/PTRC.py", line 138, in <module>\n    create_new_db(args.metadata, args.fasta, args.out)\n  File "cptr.py", line 629, in cptr.create_new_db (cptr.c:23369)\nAssertionError\nTraceback (most recent call last):\n  File "/home/Programs/PTRC/PTRC.py", line 133, in <module>\n    args.outfol, args.qual_format, args.cov_thresh, dbpath = args.db_path_name)\n  File "cptr.py", line 568, in cptr.coverage_analysis (cptr.c:20255)\n  File "cptr.py", line 226, in cptr._runGEM (cptr.c:9054)\n  File "cptr.py", line 185, in cptr._map_pe (cptr.c:7778)\n  File "cptr.py", line 203, in cptr._shell_command (cptr.c:8349)\n  File "/root/miniconda2/lib/python2.7/subprocess.py", line 574, in check_output\n    raise CalledProcessError(retcode, cmd, output=output)\nsubprocess.CalledProcessError: Command \'gem-mapper -I comm0/sim_db.gem -1 comm0/Data/comm00_1.fastq -2 comm0/Data/comm00_1.fastq -o comm0/Data/sm -q offset-33 -d all -p --max-extendable-matches all --max-extensions-per-match 5 -v\' returned non-zero exit status 1\nTraceback (most recent call last):\n  File "/home/Programs/PTRC/PTRC.py", line 136, in <module>\n    dbpath = args.db_path_name, csv_output = args.csv_output)\n  File "cptr.py", line 582, in cptr.calculate_ptrs (cptr.c:21032)\n  File "cptr.py", line 473, in cptr._load (cptr.c:16762)\n  File "/root/miniconda2/lib/python2.7/site-packages/pandas/io/pickle.py", line 65, in read_pickle\n    return try_read(path)\n  File "/root/miniconda2/lib/python2.7/site-packages/pandas/io/pickle.py", line 61, in try_read\n    with open(path, \'rb\') as fh:\nIOError: [Errno 2] No such file or directory: \'comm0/sim_db.pk\'\n'

In [747]:
cmd


Out[747]:
['/bin/bash -c "',
 'cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh',
 'comm04_1.fastq',
 'comm04_1.fastq"']

In [29]:
cmd=['/bin/bash -c "']+['cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh']+[c.name+str(0)+'_1.fastq']+[c.name+str(0)+'_1.fastq"']

In [30]:
" ".join(cmd)


Out[30]:
'/bin/bash -c " cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh comm00_1.fastq comm00_1.fastq"'

In [35]:
ptrc=pd.read_pickle(ptr_file)

In [36]:
ptrc.index.values


Out[36]:
array(['NC_000913.3\n'], dtype=object)

In [34]:
ptr_file=join(c.conf['output_path'],'ptrc_out')

In [716]:
ptrc.loc['Escherichia coli\n']['sm']


Out[716]:
1.3020408163265316

In [662]:
ptrc=compare_fit_ptrc(c)

In [685]:

Run


In [69]:
acc=['NC_000913.3','NC_007779.1','NC_002655.2','NC_009614.1','NC_017218.1']
growth_param=[[.2, .5, .3, 100],[.1, .7, .2, 200],[.2, .8, .0001,250],[0.0001, .8, .2, 250],[.1, .7, .2,150]]
sel=[[0]]
td='/Users/Shared/CommTest'

In [70]:
cli = docker.Client(base_url='unix://var/run/docker.sock')
envs = docker.utils.parse_env_file(join(td,'env'))

In [71]:
totH=[]
totK=[]
for i in range(2):
    c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2',envs)

    c.build_index()

    tot_reads=np.linspace(5*10**3,10**4,2)
    for nr in tot_reads:
        c.sample(nr)

    c.write_reads()

    totK.append(compare_fit_ptrc(c))
    #[k_c,k_err]=

    c.run_pipeline()

    #[res,hfit,pfit]=
    
    
    totH.append(compare_fit(c))

In [85]:
totH[0]


Out[85]:
[[[[<lmfit.model.ModelResult at 0x1160a4ad0>],
   [<lmfit.model.ModelResult at 0x112b1d690>]],
  array([[ 0.05322417],
         [ 0.22688601]]),
  array([[ 0.21196545],
         [ 0.34804447]])],
 [[[<lmfit.model.ModelResult at 0x1160a4110>],
   [<lmfit.model.ModelResult at 0x1160afbd0>]],
  array([[ 0.03857808],
         [ 0.17221465]]),
  array([[ 0.17485312],
         [ 0.32419588]])]]

In [88]:
tips = sns.load_dataset("tips")

In [726]:
[k_c,k_err]=compare_fit_ptrc(c)


---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-726-5d15a506d581> in <module>()
----> 1 [k_c,k_err]=compare_fit_ptrc(c)

<ipython-input-724-5b0971b9bed4> in compare_fit_ptrc(c)
      6         ptre_tmp=[]
      7         cmd=['/bin/bash -c "']+['cd /mnt/vol ; /home/Programs/PTR-Pipeline/extra/ptrc_commands.sh']+[c.name+str(i)+'_1.fastq']+[c.name+str(i)+'_2.fastq"']
----> 8         create_mount_run('centos/hedani',td," ".join(cmd),envs)
      9         ptr_file=join(c.conf['output_path'],'ptrc_out')
     10         ptrc_frame=pd.read_pickle(ptr_file)

<ipython-input-521-80ddf06e139f> in create_mount_run(image, mount_dir, cmd, envs)
     12     ctr=container.get('Id')
     13     cli.start(ctr)
---> 14     cli.wait(ctr,60*60*24*10)
     15     return cli.logs(ctr)

/Library/Python/2.7/site-packages/docker/utils/decorators.pyc in wrapped(self, resource_id, *args, **kwargs)
     19                 'image or container param is undefined'
     20             )
---> 21         return f(self, resource_id, *args, **kwargs)
     22     return wrapped
     23 

/Library/Python/2.7/site-packages/docker/api/container.pyc in wait(self, container, timeout)
    453     def wait(self, container, timeout=None):
    454         url = self._url("/containers/{0}/wait", container)
--> 455         res = self._post(url, timeout=timeout)
    456         self._raise_for_status(res)
    457         json_ = res.json()

/Library/Python/2.7/site-packages/docker/utils/decorators.pyc in inner(self, *args, **kwargs)
     45             else:
     46                 kwargs['headers'].update(self._auth_configs['HttpHeaders'])
---> 47         return f(self, *args, **kwargs)
     48     return inner

/Library/Python/2.7/site-packages/docker/client.pyc in _post(self, url, **kwargs)
    133     @update_headers
    134     def _post(self, url, **kwargs):
--> 135         return self.post(url, **self._set_request_timeout(kwargs))
    136 
    137     @update_headers

/Library/Python/2.7/site-packages/requests/sessions.pyc in post(self, url, data, json, **kwargs)
    533         """
    534 
--> 535         return self.request('POST', url, data=data, json=json, **kwargs)
    536 
    537     def put(self, url, data=None, **kwargs):

/Library/Python/2.7/site-packages/requests/sessions.pyc in request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
    486         }
    487         send_kwargs.update(settings)
--> 488         resp = self.send(prep, **send_kwargs)
    489 
    490         return resp

/Library/Python/2.7/site-packages/requests/sessions.pyc in send(self, request, **kwargs)
    607 
    608         # Send the request
--> 609         r = adapter.send(request, **kwargs)
    610 
    611         # Total elapsed time of the request (approximately)

/Library/Python/2.7/site-packages/requests/adapters.pyc in send(self, request, stream, timeout, verify, cert, proxies)
    421                     decode_content=False,
    422                     retries=self.max_retries,
--> 423                     timeout=timeout
    424                 )
    425 

/Library/Python/2.7/site-packages/requests/packages/urllib3/connectionpool.pyc in urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, **response_kw)
    592                                                   timeout=timeout_obj,
    593                                                   body=body, headers=headers,
--> 594                                                   chunked=chunked)
    595 
    596             # If we're going to release the connection in ``finally:``, then

/Library/Python/2.7/site-packages/requests/packages/urllib3/connectionpool.pyc in _make_request(self, conn, method, url, timeout, chunked, **httplib_request_kw)
    382         try:
    383             try:  # Python 2.7, use buffering of HTTP responses
--> 384                 httplib_response = conn.getresponse(buffering=True)
    385             except TypeError:  # Python 2.6 and older, Python 3
    386                 try:

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/httplib.pyc in getresponse(self, buffering)
   1134 
   1135         try:
-> 1136             response.begin()
   1137             assert response.will_close != _UNKNOWN
   1138             self.__state = _CS_IDLE

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/httplib.pyc in begin(self)
    451         # read until we get a non-100 response
    452         while True:
--> 453             version, status, reason = self._read_status()
    454             if status != CONTINUE:
    455                 break

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/httplib.pyc in _read_status(self)
    407     def _read_status(self):
    408         # Initialize with Simple-Response defaults
--> 409         line = self.fp.readline(_MAXLINE + 1)
    410         if len(line) > _MAXLINE:
    411             raise LineTooLong("header line")

/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/socket.pyc in readline(self, size)
    478             while True:
    479                 try:
--> 480                     data = self._sock.recv(self._rbufsize)
    481                 except error, e:
    482                     if e.args[0] == EINTR:

KeyboardInterrupt: 

In [708]:
plt.plot(tot_reads,np.mean(hfit,1),tot_reads,np.mean(ptrc_e,1))#,tot_reads,[r[0].best_values['C'],r[1].best_values['C'] for r in res])
legend = plt.legend(['H','K'])
#ax = plt.axes()
#ax.set_title('Mean absolute error')
#g.fig.suptitle('THIS IS A TITLE, YOU BET')



In [ ]:
sns_plot.savefig("output.pdf")

In [280]:
args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join('/mnt/vol','comm0','project.conf'),'-e','hedani@chalmers.se']          
create_mount_run('centos/hedani',td," ".join(['/bin/bash -c "']+args0+['make']+[';']+args0+['run"']),envs)


Out[280]:
"Running alignment stage\ncomm00:\nAligning comm00:\n99999 reads; of these:\n  99999 (100.00%) were paired; of these:\n    99690 (99.69%) aligned concordantly 0 times\n    147 (0.15%) aligned concordantly exactly 1 time\n    162 (0.16%) aligned concordantly >1 times\n    ----\n    99690 pairs aligned concordantly 0 times; of these:\n      47312 (47.46%) aligned discordantly 1 time\n    ----\n    52378 pairs aligned 0 times concordantly or discordantly; of these:\n      104756 mates make up the pairs; of these:\n        57 (0.05%) aligned 0 times\n        4319 (4.12%) aligned exactly 1 time\n        100380 (95.82%) aligned >1 times\n99.97% overall alignment rate\nMoving into subjobs\ncomm00: Running pathoscope stage\nRunning without mySQLdb library\ncomm00: Running coverage stage\ncomm00: Sorting\ncomm00: Calculating coverage\ncomm00: Performing piecewise fits\nNC_000913.3\nNC_000913.3.depth\nNC_000913.3.depth: Bin coverage too low, exiting.\nNC_009614.1\nNC_009614.1.depth\nNC_009614.1.depth: Coverage OK (0.813975626304 x).\n\ncntr: 6\nNC_002655.2\nNC_002655.2.depth\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:65: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(iEnd);\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:72: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(np.floor(len(x)/slide));\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:65: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(iEnd);\n/home/Programs/PTR-Pipeline/bin/piecewiseFit.py:72: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n  out=np.zeros(np.floor(len(x)/slide));\nmv: cannot stat 'depth/*.log': No such file or directory\nNC_002655.2.depth: Coverage OK (0.810297289744 x).\n\ncntr: 8\ncomm00: Performing secondary fits\n*\nI/O error *.npy\n"

In [507]:
res=fit_signal(st,1)
err_hfit=((.5-res.best_values['C'])/.5)
err_hfit


Out[507]:
0.8969818062406385

In [464]:
res.best_values['C']


Out[464]:
0.50223664098758436

In [496]:
#depth_file=join(c.conf['output_path'],c.name+'0','npy',c.pop.keys()[0]+'.depth.npy')
depth_file='/Users/Shared/CommTest/comm0/Out/comm00/npy/NC_000913.3.depth.npy'

In [508]:
signal=2**np.load(depth_file)
signal=signal/(signal.sum()/len(signal))

In [509]:
plt.plot(signal)


Out[509]:
[<matplotlib.lines.Line2D at 0x137fd01d0>]

In [429]:
x = np.linspace(-2*np.pi, 2*np.pi, len(signal))

In [339]:
xx=np.linspace(-np.pi, np.pi, 1000)
plt.plot(np.abs(mv_filt(1.2,xx)))


Out[339]:
[<matplotlib.lines.Line2D at 0x12cf47f10>]

In [401]:
y2=[]
x2=np.linspace(0,.5,100)
y2=np.append(pxn2(x2,.5),pxn2(x2,.5)[::-1])
x2=np.linspace(0,1,200)

y=list(triangle(1000, 0.1,10))
x=range(len(y))
x_p = np.linspace(-2*np.pi, 2*np.pi, len(y2))

f=mv_filt(20,x_p)
fft_y=fft(y2)*f
smooth_y=np.roll(ifft(fft_y),0)
smooth_y=smooth_y/np.sum(smooth_y)

x3=np.linspace(0,1,len(y_long))
y_skum=y_long/(np.sum(y_long+np.mean(y_long)))
x_p2 = np.linspace(-2*np.pi, 2*np.pi, len(y_skum))

plt.plot(x2,y2/np.sum(y2),x2,smooth_y,x2,ifft(fft(smooth_y)/f))


Out[401]:
[<matplotlib.lines.Line2D at 0x1325a21d0>,
 <matplotlib.lines.Line2D at 0x1325a2490>,
 <matplotlib.lines.Line2D at 0x1325a2b50>]

In [505]:
xi=np.linspace(0,1,len(signal))
xi2=np.linspace(0,1,10**3)
iss=scipy.interpolate.UnivariateSpline(xi,signal,k=3)
iss.set_smoothing_factor(.00001)
y_long=iss(xi2)
plt.plot(y_long)


Out[505]:
[<matplotlib.lines.Line2D at 0x137ce6410>]

In [208]:
plt.plot(np.abs(ifft(fft(signal,10**4))[0:len(signal)]))


Out[208]:
[<matplotlib.lines.Line2D at 0x121a0eed0>]

In [512]:
signal=y_long
x = np.linspace(-0.00001*np.pi/float(1), 1*np.pi/float(1),10**3)
st=np.roll(ifft(fft(signal,10**3)/mv_filt(1.27,x))[0:len(signal)],np.int64(-125/2))
x2=np.linspace(0,1,len(st))
x1=np.linspace(0,1,len(signal))
plt.plot(x2,st/np.sum(st),x1,signal/np.sum(signal))


Out[512]:
[<matplotlib.lines.Line2D at 0x1383f14d0>,
 <matplotlib.lines.Line2D at 0x1383f1790>]

In [514]:
plt.plot(x2,st/np.sum(st),x1,signal/np.sum(signal))


Out[514]:
[<matplotlib.lines.Line2D at 0x1384a5c10>,
 <matplotlib.lines.Line2D at 0x1384a5ed0>]

In [193]:
plt.plot(np.abs(mv_filt(1*10**1,x)))


Out[193]:
[<matplotlib.lines.Line2D at 0x12058b790>]

In [411]:
np.sqrt(-1)


Out[411]:
nan

In [ ]:
1+np.exp()

In [398]:
np.exp(1)


Out[398]:
2.7182818284590451

In [272]:
create_mount_run('centos/hedani',td," ".join(['/bin/bash -c "']+['which python"']),envs)


Out[272]:
'/usr/bin/python\n'

In [52]:
tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
c.build_index()
for nr in tot_reads:
    c.sample(nr)
    c.write_reads()

c.run_pipeline()
[res,hfit,pfit]=c.compare_fit()


---------------------------------------------------------------------------
NotFound                                  Traceback (most recent call last)
<ipython-input-52-071c0a386248> in <module>()
      1 tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
----> 2 c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
      3 c.build_index()
      4 for nr in tot_reads:
      5     c.sample(nr)

<ipython-input-31-7a930caef3da> in __init__(self, name, acc, growth_param, td, mapper)
      6 
      7         self.create_dirs(True)
----> 8         self.fetch_ref(acc)
      9 
     10         self.pop = self.init_population(acc,growth_param)

<ipython-input-31-7a930caef3da> in fetch_ref(self, acc)
    100         f.close()
    101 
--> 102         create_mount_run('centos/hedani','/Users/Shared/d_py',['-e','hedani@chalmers.se','fetch-references','-s',acc_path])
    103         #args=ptr_pipeline.parse_args(['-e','hedani@chalmers.se','fetch-references','-s',acc_path])
    104         #ptr_pipeline.main(args,self.conf)

<ipython-input-28-257136503bbf> in create_mount_run(image, mount_dir, cmd)
      9         })
     10     )
---> 11     cli.start(container=container.get('Id'))
     12     return cli.logs(container=container.get('Id'))

/Library/Python/2.7/site-packages/docker/utils/decorators.pyc in wrapped(self, resource_id, *args, **kwargs)
     19                 'image or container param is undefined'
     20             )
---> 21         return f(self, resource_id, *args, **kwargs)
     22     return wrapped
     23 

/Library/Python/2.7/site-packages/docker/api/container.pyc in start(self, container, binds, port_bindings, lxc_conf, publish_all_ports, links, privileged, dns, dns_search, volumes_from, network_mode, restart_policy, cap_add, cap_drop, devices, extra_hosts, read_only, pid_mode, ipc_mode, security_opt, ulimits)
    381         url = self._url("/containers/{0}/start", container)
    382         res = self._post_json(url, data=start_config)
--> 383         self._raise_for_status(res)
    384 
    385     @utils.minimum_version('1.17')

/Library/Python/2.7/site-packages/docker/client.pyc in _raise_for_status(self, response, explanation)
    171         except requests.exceptions.HTTPError as e:
    172             if e.response.status_code == 404:
--> 173                 raise errors.NotFound(e, response, explanation=explanation)
    174             raise errors.APIError(e, response, explanation=explanation)
    175 

NotFound: 404 Client Error: Not Found ("{"message":"oci runtime error: exec: \"-e\": executable file not found in $PATH"}")

In [ ]:


In [247]:
community=[]
nrs=np.linspace(10**5,10**5,1)
res=[]
for nr in nrs:
    community.append(run_community_analysis('fs',td,'bowtie2',acc,growth_param,sel,nr,0))
    #res.append(community[i][0].compare_fit())


ERROR:root:Internal Python error in the inspect module.
Below is the traceback from this internal error.

Traceback (most recent call last):
  File "/Library/Python/2.7/site-packages/IPython/core/ultratb.py", line 1132, in get_records
    return _fixed_getinnerframes(etb, number_of_lines_of_context, tb_offset)
  File "/Library/Python/2.7/site-packages/IPython/core/ultratb.py", line 313, in wrapped
    return f(*args, **kwargs)
  File "/Library/Python/2.7/site-packages/IPython/core/ultratb.py", line 358, in _fixed_getinnerframes
    records = fix_frame_records_filenames(inspect.getinnerframes(etb, context))
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 1049, in getinnerframes
    framelist.append((tb.tb_frame,) + getframeinfo(tb, context))
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 1009, in getframeinfo
    filename = getsourcefile(frame) or getfile(frame)
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 454, in getsourcefile
    if hasattr(getmodule(object, filename), '__loader__'):
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 483, in getmodule
    file = getabsfile(object, _filename)
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/inspect.py", line 467, in getabsfile
    return os.path.normcase(os.path.abspath(_filename))
  File "/usr/local/Cellar/python/2.7.11/Frameworks/Python.framework/Versions/2.7/lib/python2.7/posixpath.py", line 364, in abspath
    cwd = os.getcwd()
OSError: [Errno 2] No such file or directory
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/Library/Python/2.7/site-packages/IPython/core/interactiveshell.pyc in run_code(self, code_obj, result)
   2896             if result is not None:
   2897                 result.error_in_exec = sys.exc_info()[1]
-> 2898             self.showtraceback()
   2899         else:
   2900             outflag = 0

/Library/Python/2.7/site-packages/IPython/core/interactiveshell.pyc in showtraceback(self, exc_tuple, filename, tb_offset, exception_only)
   1822                     except Exception:
   1823                         stb = self.InteractiveTB.structured_traceback(etype,
-> 1824                                             value, tb, tb_offset=tb_offset)
   1825 
   1826                     self._showtraceback(etype, value, stb)

/Library/Python/2.7/site-packages/IPython/core/ultratb.pyc in structured_traceback(self, etype, value, tb, tb_offset, number_of_lines_of_context)
   1404         self.tb = tb
   1405         return FormattedTB.structured_traceback(
-> 1406             self, etype, value, tb, tb_offset, number_of_lines_of_context)
   1407 
   1408 

/Library/Python/2.7/site-packages/IPython/core/ultratb.pyc in structured_traceback(self, etype, value, tb, tb_offset, number_of_lines_of_context)
   1312             # Verbose modes need a full traceback
   1313             return VerboseTB.structured_traceback(
-> 1314                 self, etype, value, tb, tb_offset, number_of_lines_of_context
   1315             )
   1316         else:

/Library/Python/2.7/site-packages/IPython/core/ultratb.pyc in structured_traceback(self, etype, evalue, etb, tb_offset, number_of_lines_of_context)
   1196                 structured_traceback_parts += formatted_exception
   1197         else:
-> 1198             structured_traceback_parts += formatted_exception[0]
   1199 
   1200         return structured_traceback_parts

IndexError: string index out of range

In [73]:
c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')


---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-73-fde9aa008fdc> in <module>()
----> 1 c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')

<ipython-input-55-d175fa22a917> in __init__(self, name, acc, growth_param, td, mapper)
     12         self.fetch_ref(acc)
     13 
---> 14         self.pop = self.init_population(acc,growth_param)
     15         self.distribution = self.community_distribution()
     16         self.samples = []

<ipython-input-55-d175fa22a917> in init_population(self, acc, growth_param)
     29         for i,a in enumerate(acc):
     30             pop[a]=population(B = growth_param[i][0],C = growth_param[i][1],
---> 31                              D = growth_param[i][2], l = len(records[a]),
     32                              seq = records[a], cells=growth_param[i][3])
     33         return pop

KeyError: 'NC_000913.3'

In [47]:
tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
c.build_index()


---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-47-5b05c90de0d8> in <module>()
      1 tot_reads=[10**5]#np.linspace(5*10**3,10**4,15)
----> 2 c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
      3 c.build_index()

<ipython-input-29-c0505f730c20> in __init__(self, name, acc, growth_param, td, mapper)
      7         self.args0=['/home/Programs/PTR-Pipeline/ptr_pipeline.py','-c',join(self.d_conf['node_path'],'project.conf'),'-e','hedani@chalmers.se']
      8 
----> 9         save_config(self.d_conf,self.conf)
     10 
     11         self.create_dirs(True)

<ipython-input-43-68021bfe077a> in save_config(lconf, conf)
      3     Config.optionxform = str
      4 
----> 5     cfgfile = open(join(lconf['node_path'],'project.conf'),'w')
      6 
      7     Config.add_section('Project')

IOError: [Errno 2] No such file or directory: '/mnt/vol/comm0/project.conf'

In [2132]:
n=130;
D1_snok=np.random.randint(30,40);C1_snok=np.random.randint(40,70);tau1_snok=np.random.randint(130,250,n)
#tau1_snok[-1]=50;
D2_snok=10.;C2_snok=35.;tau2_snok=np.random.randint(60,200,n)
ab1=[np.random.rand() for x in range(n)];ab2=np.random.randint(100,200,n);
org1=[[np.max([1-(C1_snok+D1_snok)/float(tau1_snok[i]),0]),(C1_snok)/float(tau1_snok[i]), (D1_snok)/float(tau1_snok[i]),ab1[i]] for i in range(tau1_snok.size)]
org2=[[1-(C2_snok+D2_snok)/float(tau2_snok[i]),(C2_snok)/float(tau2_snok[i]), (D1_snok)/float(tau2_snok[i]),ab2[i]/float(ab1[i])] for i in range(tau2_snok.size)]

#same_org[0][0]=0

In [2133]:
print D1_snok,C1_snok


34 44

In [2134]:
same_org=np.concatenate((org1,org2))

In [2135]:
Ri1


Out[2135]:
array([ 0.12517895,  0.28003121,  0.11440445,  0.15001599,  0.07173145,
        0.32649918,  0.39299506,  0.01202917,  0.34304507,  0.33282179,
        0.04329798,  0.26181841,  0.33872681,  0.14846477,  0.24897509,
        0.061663  ,  0.2401119 ,  0.40135043,  0.2886033 ,  0.39829411,
        0.1850574 ,  0.34708688,  0.20622304,  0.16787013,  0.16767945,
        0.37298813,  0.33103408,  0.25134168,  0.35186729,  0.16303172,
        0.30711971,  0.26330621,  0.12609793,  0.20625785,  0.3764792 ,
        0.22850401,  0.0450059 ,  0.0809682 ,  0.18655135,  0.00784618,
        0.29397322,  0.08577948,  0.29319743,  0.10284626,  0.30304625,
        0.28011791,  0.1960176 ,  0.16534839,  0.02916867,  0.0991715 ,
        0.0578688 ,  0.27499659,  0.027967  ,  0.13073409,  0.15571648,
        0.01069903,  0.1112608 ,  0.18262335,  0.20584184,  0.3887354 ,
        0.38976375,  0.38220299,  0.36039617,  0.10669672,  0.22368446,
        0.23371822,  0.03268177,  0.25700886,  0.07782686,  0.37210099,
        0.11653302,  0.06587968,  0.14566878,  0.0068028 ,  0.22986611,
        0.00415846,  0.21185567,  0.23742862,  0.36762303,  0.07634184,
        0.11859662,  0.3574003 ,  0.13858283,  0.08572636,  0.01833882,
        0.19050518,  0.34799022,  0.24642473,  0.10322462,  0.31964564,
        0.07285449,  0.36825552,  0.09071813,  0.18713142,  0.23090128,
        0.19523249,  0.24599339,  0.2281989 ,  0.0326502 ,  0.04275451,
        0.29883877,  0.34733176,  0.25028025,  0.24900968,  0.01261928,
        0.3568365 ,  0.01669861,  0.29840329,  0.10386471,  0.27002372,
        0.04221017,  0.23762403,  0.08465263,  0.40777278,  0.09767022,
        0.12452879,  0.31985767,  0.09076833,  0.15271323,  0.30922279,
        0.2969826 ,  0.29776965,  0.21090022,  0.3868555 ,  0.17569115,
        0.27928653,  0.43143433,  0.0779982 ,  0.01669648,  0.10084281])

In [2136]:
Ri1=np.array([Gekv2(p[1],p[2])*1*p[3] for p in org1])
rs=np.random.randint(2,3,n);
for i in range(n):
    Ri1[i]=Ri1[i]/rs[i];
    org1[i][3]=org1[i][3]/rs[i];


Ki1=np.array([(2**(p[1])-1)/(float(p[1])*np.log(2)) for p in org1])
Ri2=np.array([Gekv2(p[1],p[2])*1*p[3] for p in org2])
Ki2=np.array([(2**(p[1])-1)/(float(p[1])*np.log(2)) for p in org2])

In [2137]:
org1


Out[2137]:
[[0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.46079501691747554],
 [0.59375, 0.22916666666666666, 0.17708333333333334, 0.20116365324132957],
 [0.67088607594936711,
  0.18565400843881857,
  0.14345991561181434,
  0.49654154079404594],
 [0.57837837837837835,
  0.23783783783783785,
  0.1837837837837838,
  0.24764121123035404],
 [0.60999999999999999, 0.22, 0.17, 0.3921142465319019],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.20158666914294182],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.22483740954130388],
 [0.65638766519823788,
  0.19383259911894274,
  0.14977973568281938,
  0.22464293403338875],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.44063274576640726],
 [0.60204081632653061,
  0.22448979591836735,
  0.17346938775510204,
  0.43008967519542796],
 [0.60606060606060608,
  0.2222222222222222,
  0.1717171717171717,
  0.17211265559073496],
 [0.52439024390243905,
  0.2682926829268293,
  0.2073170731707317,
  0.2866343753302563],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.17988289772514116],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.0918123170900288],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.4202347320171348],
 [0.58730158730158732,
  0.2328042328042328,
  0.17989417989417988,
  0.4089650647473905],
 [0.59999999999999998,
  0.22564102564102564,
  0.17435897435897435,
  0.4737860668884134],
 [0.68163265306122445,
  0.17959183673469387,
  0.13877551020408163,
  0.18682768873395322],
 [0.53846153846153844,
  0.2603550295857988,
  0.20118343195266272,
  0.46915131141688116],
 [0.58947368421052637,
  0.23157894736842105,
  0.17894736842105263,
  0.14030413639900796],
 [0.53012048192771077,
  0.26506024096385544,
  0.20481927710843373,
  0.01801941652267519],
 [0.64220183486238525,
  0.2018348623853211,
  0.1559633027522936,
  0.46220745532265606],
 [0.64055299539170507,
  0.20276497695852536,
  0.15668202764976957,
  0.11793696321699298],
 [0.40000000000000002,
  0.3384615384615385,
  0.26153846153846155,
  0.016939443057983616],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.02568227640343973],
 [0.43065693430656937,
  0.32116788321167883,
  0.24817518248175183,
  0.22051022260961034],
 [0.60406091370558368,
  0.2233502538071066,
  0.17258883248730963,
  0.012433187891947306],
 [0.6080402010050252,
  0.22110552763819097,
  0.1708542713567839,
  0.2914623529141575],
 [0.46938775510204078,
  0.29931972789115646,
  0.23129251700680273,
  0.0003985532846665918],
 [0.61194029850746268,
  0.21890547263681592,
  0.1691542288557214,
  0.4267106671672934],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.4075861771070389],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.17460659919835086],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.23454114179864116],
 [0.63888888888888884,
  0.2037037037037037,
  0.1574074074074074,
  0.07579110856008303],
 [0.5357142857142857,
  0.2619047619047619,
  0.20238095238095238,
  0.0652801694669164],
 [0.66808510638297869,
  0.18723404255319148,
  0.14468085106382977,
  0.27763616926454493],
 [0.68674698795180722,
  0.17670682730923695,
  0.13654618473895583,
  0.022855740658188184],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.49074281480303383],
 [0.65638766519823788,
  0.19383259911894274,
  0.14977973568281938,
  0.02739239734518567],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.4917295420246727],
 [0.66233766233766234,
  0.19047619047619047,
  0.1471861471861472,
  0.2497526462054266],
 [0.6080402010050252,
  0.22110552763819097,
  0.1708542713567839,
  0.4770605142372987],
 [0.43478260869565222,
  0.3188405797101449,
  0.2463768115942029,
  0.28254111584207114],
 [0.45833333333333337,
  0.3055555555555556,
  0.2361111111111111,
  0.010230785310108625],
 [0.61576354679802958,
  0.21674876847290642,
  0.16748768472906403,
  0.291575700978587],
 [0.62318840579710144,
  0.21256038647342995,
  0.1642512077294686,
  0.1643615088854778],
 [0.47651006711409394,
  0.2953020134228188,
  0.22818791946308725,
  0.285877263586],
 [0.4135338345864662,
  0.3308270676691729,
  0.2556390977443609,
  0.4079787945674634],
 [0.49677419354838714,
  0.2838709677419355,
  0.21935483870967742,
  0.21611641119296693],
 [0.59999999999999998,
  0.22564102564102564,
  0.17435897435897435,
  0.31412258182920005],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.45168765515216625],
 [0.68421052631578949,
  0.17813765182186234,
  0.13765182186234817,
  0.2223351304214118],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.06362377612636472],
 [0.68032786885245899,
  0.18032786885245902,
  0.13934426229508196,
  0.31736353876696866],
 [0.50632911392405067,
  0.27848101265822783,
  0.21518987341772153,
  0.011060227802662803],
 [0.63207547169811318,
  0.20754716981132076,
  0.16037735849056603,
  0.40062911535554124],
 [0.62318840579710144,
  0.21256038647342995,
  0.1642512077294686,
  0.044284671774929896],
 [0.62318840579710144,
  0.21256038647342995,
  0.1642512077294686,
  0.017877932941876062],
 [0.4388489208633094,
  0.31654676258992803,
  0.2446043165467626,
  0.27840805982487965],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.4303743176937124],
 [0.57608695652173914,
  0.2391304347826087,
  0.18478260869565216,
  0.4789520760689538],
 [0.62857142857142856,
  0.20952380952380953,
  0.1619047619047619,
  0.39675886455683007],
 [0.58288770053475936,
  0.23529411764705882,
  0.18181818181818182,
  0.17842565107350478],
 [0.45454545454545459,
  0.3076923076923077,
  0.23776223776223776,
  0.15424981483099498],
 [0.68292682926829262,
  0.17886178861788618,
  0.13821138211382114,
  0.4159795866799625],
 [0.47297297297297303,
  0.2972972972972973,
  0.22972972972972974,
  0.0400130749241025],
 [0.67088607594936711,
  0.18565400843881857,
  0.14345991561181434,
  0.16761038698281472],
 [0.55932203389830515,
  0.24858757062146894,
  0.192090395480226,
  0.06748110497628335],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.33515107188402804],
 [0.53012048192771077,
  0.26506024096385544,
  0.20481927710843373,
  0.37901848137513133],
 [0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.32363978989470665],
 [0.64055299539170507,
  0.20276497695852536,
  0.15668202764976957,
  0.32810332357826466],
 [0.52727272727272734,
  0.26666666666666666,
  0.20606060606060606,
  0.4715145827317281],
 [0.56666666666666665,
  0.24444444444444444,
  0.18888888888888888,
  0.004128507612861987],
 [0.61764705882352944,
  0.21568627450980393,
  0.16666666666666666,
  0.25827871302554967],
 [0.67226890756302526,
  0.18487394957983194,
  0.14285714285714285,
  0.4987618526306502],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.47396175599898605],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.22916006140325845],
 [0.66808510638297869,
  0.18723404255319148,
  0.14468085106382977,
  0.006042670814164486],
 [0.61764705882352944,
  0.21568627450980393,
  0.16666666666666666,
  0.31867207806094916],
 [0.63033175355450233,
  0.20853080568720378,
  0.16113744075829384,
  0.4533736478201163],
 [0.56906077348066297,
  0.2430939226519337,
  0.1878453038674033,
  0.2480520511184005],
 [0.68163265306122445,
  0.17959183673469387,
  0.13877551020408163,
  0.16669709252320397],
 [0.57608695652173914,
  0.2391304347826087,
  0.18478260869565216,
  0.46919015186765267],
 [0.63551401869158886,
  0.205607476635514,
  0.1588785046728972,
  0.380183644923833],
 [0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.2307091141351154],
 [0.58064516129032251,
  0.23655913978494625,
  0.1827956989247312,
  0.16310104203088738],
 [0.55681818181818188, 0.25, 0.19318181818181818, 0.29090402131285625],
 [0.55932203389830515,
  0.24858757062146894,
  0.192090395480226,
  0.3036843665003522],
 [0.66949152542372881,
  0.1864406779661017,
  0.1440677966101695,
  0.3260796666861614],
 [0.67901234567901236,
  0.18106995884773663,
  0.13991769547325103,
  0.25079441129544927],
 [0.46938775510204078,
  0.29931972789115646,
  0.23129251700680273,
  0.05843929831879591],
 [0.56906077348066297,
  0.2430939226519337,
  0.1878453038674033,
  0.3788527812958234],
 [0.51552795031055898,
  0.2732919254658385,
  0.2111801242236025,
  0.2422078913952499],
 [0.60999999999999999, 0.22, 0.17, 0.4672477014175954],
 [0.56666666666666665,
  0.24444444444444444,
  0.18888888888888888,
  0.42521132206742074],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.10407675396895005],
 [0.58947368421052637,
  0.23157894736842105,
  0.17894736842105263,
  0.47486698927718024],
 [0.68163265306122445,
  0.17959183673469387,
  0.13877551020408163,
  0.07495626881189427],
 [0.42647058823529416, 0.3235294117647059, 0.25, 0.48954133636693964],
 [0.54117647058823537, 0.25882352941176473, 0.2, 0.3871470533255098],
 [0.46575342465753422,
  0.3013698630136986,
  0.2328767123287671,
  0.2675470437722781],
 [0.65333333333333332,
  0.19555555555555557,
  0.1511111111111111,
  0.3255599590818351],
 [0.64864864864864868,
  0.1981981981981982,
  0.15315315315315314,
  0.45573537792144153],
 [0.63720930232558137,
  0.20465116279069767,
  0.15813953488372093,
  0.17392641682893523],
 [0.6637931034482758,
  0.1896551724137931,
  0.14655172413793102,
  0.2968118195696663],
 [0.65333333333333332,
  0.19555555555555557,
  0.1511111111111111,
  0.07989447253978554],
 [0.51552795031055898,
  0.2732919254658385,
  0.2111801242236025,
  0.49371289448667205],
 [0.47999999999999998,
  0.29333333333333333,
  0.22666666666666666,
  0.09733296432041799],
 [0.5185185185185186,
  0.2716049382716049,
  0.20987654320987653,
  0.26059767520839067],
 [0.46938775510204078,
  0.29931972789115646,
  0.23129251700680273,
  0.14626166795275647],
 [0.66233766233766234,
  0.19047619047619047,
  0.1471861471861472,
  0.1628807839145976],
 [0.60999999999999999, 0.22, 0.17, 0.16708772481168682],
 [0.65789473684210531,
  0.19298245614035087,
  0.14912280701754385,
  0.26328747680426245],
 [0.68674698795180722,
  0.17670682730923695,
  0.13654618473895583,
  0.1593009524224957],
 [0.56424581005586594,
  0.24581005586592178,
  0.18994413407821228,
  0.460562307244909],
 [0.61386138613861385,
  0.21782178217821782,
  0.16831683168316833,
  0.28121643054034673],
 [0.48684210526315785,
  0.2894736842105263,
  0.2236842105263158,
  0.1253352606794409],
 [0.45833333333333337,
  0.3055555555555556,
  0.2361111111111111,
  0.010656670378031252],
 [0.56666666666666665,
  0.24444444444444444,
  0.18888888888888888,
  0.44690907408942854],
 [0.53012048192771077,
  0.26506024096385544,
  0.20481927710843373,
  0.31841682872717114],
 [0.65333333333333332,
  0.19555555555555557,
  0.1511111111111111,
  0.4744375311911005],
 [0.68292682926829262,
  0.17886178861788618,
  0.13821138211382114,
  0.31763678516257166],
 [0.49350649350649356,
  0.2857142857142857,
  0.22077922077922077,
  0.3042811363098526],
 [0.61951219512195121,
  0.2146341463414634,
  0.16585365853658537,
  0.10414107380361959],
 [0.40000000000000002,
  0.3384615384615385,
  0.26153846153846155,
  0.4481668632525234],
 [0.55932203389830515,
  0.24858757062146894,
  0.192090395480226,
  0.27911572937343504],
 [0.60406091370558368,
  0.2233502538071066,
  0.17258883248730963,
  0.19864426249618172],
 [0.50943396226415094,
  0.27672955974842767,
  0.2138364779874214,
  0.19862820800678094],
 [0.52727272727272734,
  0.26666666666666666,
  0.20606060606060606,
  0.22704483674014786]]

In [1887]:
np.array([Gekv2(p[1],p[2])*1*p[3] for p in org1])


Out[1887]:
array([ 1.09875756,  0.54429242,  0.14761133])

In [1888]:
[2**(p[2])*((2**(p[1])-1)/(float(p[1])*np.log(2)))*1*p[3]/400. for p in same_org]


Out[1888]:
[0.002746893895910563,
 0.0013607310513527526,
 0.00036902832344279818,
 0.58984513209798939,
 1.0191773861094746,
 5.2424783744028893]

In [939]:
cip


Out[939]:
array([[ 0.23809524,  0.18404908],
       [ 0.23809524,  0.1875    ],
       [ 0.18404908,  0.1875    ]])

In [546]:
[c[1]/float(c[0]) for c in cip]


Out[546]:
[1.3636363636363635, 2.9104477611940296, 2.1343283582089554]

In [545]:
[t[0]/float(t[1]) for t in tip]


Out[545]:
[1.3636363636363635, 2.91044776119403, 2.1343283582089554]

In [2138]:
Rip1=create_pairs(Ri1)
Kip1=create_pairs(Ki1)
cip1=create_pairs([p[1] for p in org1])
tip1=create_pairs([t for t in tau1_snok])
Aip1=create_pairs([p[3] for p in org1])
Rip2=create_pairs(Ri2)
Kip2=create_pairs(Ki2)
cip2=create_pairs([p[1] for p in org2])
tip2=create_pairs([t for t in tau2_snok])
Aip2=create_pairs([p[3] for p in org2])
inde=np.triu_indices(len(Ri1),k=1)
#Aa=np.zeros(len(Rip))

#for i in len(cip):
#    Aa[i]=A_snok(Rip[i],Kip[i],[1,1],cip[i])
    
#M=np.zeros(len(sel[0]))
#for i in len(arrC):
#    M[i][j]=1-

In [1940]:
cip


Out[1940]:
array([[ 0.23809524,  0.18404908],
       [ 0.23809524,  0.1875    ],
       [ 0.18404908,  0.1875    ]])

In [2139]:
M1=np.zeros((len(Rip1),len(Ri1)))
b1=np.zeros(len(Rip1))
for i in range(len(inde[0])):
    cc=cip1[i][1]/cip1[i][0]
    M1[i][inde[0][i]]=cc
    M1[i][inde[1][i]]=-1
    b1[i]=cc*np.log2(Rip1[i][0])-np.log2(Rip1[i][1])+np.log2(Kip1[i][1])-cc*np.log2(Kip1[i][0])

In [2140]:
M2=np.zeros((len(Rip2),len(Ri2)))
b2=np.zeros(len(Rip2))
for i in range(len(inde[0])):
    cc=cip2[i][1]/cip2[i][0]
    M2[i][inde[0][i]]=cc
    M2[i][inde[1][i]]=-1
    b2[i]=cc*np.log2(Rip2[i][0])-np.log2(Rip2[i][1])+np.log2(Kip2[i][1])-cc*np.log2(Kip2[i][0])

In [1807]:
for i in range(3):
    print (Rip1[i][0]/Kip1[i][0])**(cip1[i][1]/cip[i][0])*Kip1[i][1]/Rip1[i][1]


1.38762708215
1.07528402889
1.69467827533

In [1674]:
for i in range(3):
    print (Aip1[i][0]**(cip1[i][1]/cip[i][0]))/Aip1[i][1]


0.449997988656
3.23668895758
0.358153810534

In [1505]:
for i in range(3):
    print cip1[i][1]/cip[i][0]


2.08108108108
1.57142857143
2.1654589372

In [881]:
print Kip1[1]
print Rip1[1]
print cip1[]
dd


[ 1.11527593  1.24003622]
[ 0.3560757   0.92593291]
[ 0.30927835  0.6       ]
Out[881]:
array([ 0.27675277,  0.48214286,  0.56589147])

In [870]:
print b1


[-0.36306836 -2.77405818 -1.85186455]

In [738]:
Z1=np.zeros_like(M1);
Z2=np.zeros_like(M2);
B=np.ones(3*2);

T=np.bmat([[M1,Z1],[Z2,M2]]);
T=np.vstack((T,B));
b=np.vstack((b1,b2));

#b_sum=np.sum([-2/float(3)*x for x in np.array(org1)[:,1]]+Ri1-Ki1)
#b_sum=b_sum+np.sum([-2/float(3)*x for x in np.array(org2)[:,1]]+Ri2-Ki2)

b_sum=np.sum(-np.array(org1)[:,2]+Ri1-Ki1)
b_sum=b_sum+np.sum(-np.array(org2)[:,2]+Ri2-Ki2)


b=np.append(b,b_sum)

In [713]:
aa=np.array(org1)[:,1]

In [726]:
b_sum=np.sum([-2/float(3)*x for x in np.array(org1)[:,1]]+Ri1-Ki1)
b_sum=b_sum+np.sum([-2/float(3)*x for x in np.array(org2)[:,1]]+Ri2-Ki2)

In [727]:
b_sum


Out[727]:
-3.987251671839438

In [761]:
dd


Out[761]:
array([ 0.35869565,  0.28888889,  0.29962547,  0.64130435,  0.71111111,
        0.70037453])

In [511]:
sort_cip=[]*(len(cip*len(cip)))
for i in range(len(cip)):
    for j in range(len(cip)):
        sort_cip[i*len(cip)+j]=cip[i][1]/cip[i][0]*cip[j][1]/cip[j][0]
        
sort_cip=np.sort(sort_cip)


---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-511-17bee1cd483d> in <module>()
      2 for i in range(len(cip)):
      3     for j in range(len(cip)):
----> 4         sort_cip[i*len(cip)+j]=cip[i][1]/cip[i][0]*cip[j][1]/cip[j][0]
      5 
      6 sort_cip=np.sort(sort_cip)

IndexError: list assignment index out of range

In [506]:
F=np.vstack((M[0],M[5],M[9],M[12],M[13],M[14]))

In [1675]:
for i in range(len(cip1)):
    print Aip1[i][1]*Kip1[i][1]/Rip1[i][1]-(Aip1[i][0]*Kip1[i][0]/Rip1[i][0])**(cip1[i][1]/cip1[i][0])


0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
0.0
-1.11022302463e-16
2.22044604925e-16
1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
0.0
0.0
0.0
0.0
1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
0.0
2.22044604925e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
0.0
-1.11022302463e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-1.11022302463e-16
0.0
1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
0.0
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
0.0
0.0
1.11022302463e-16
0.0
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
3.33066907388e-16
2.22044604925e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
2.22044604925e-16
0.0
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
2.22044604925e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
2.22044604925e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
0.0
1.11022302463e-16
0.0
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
1.11022302463e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
3.33066907388e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
2.22044604925e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
1.11022302463e-16
0.0
0.0
-1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
3.33066907388e-16
2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
2.22044604925e-16
0.0
0.0
1.11022302463e-16
0.0
0.0
2.22044604925e-16
1.11022302463e-16
0.0
0.0
0.0
2.22044604925e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
1.11022302463e-16
0.0
2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-2.22044604925e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
1.11022302463e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
3.33066907388e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
2.22044604925e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
0.0
0.0
1.11022302463e-16
1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
0.0
-3.33066907388e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
3.33066907388e-16
2.22044604925e-16
2.22044604925e-16
3.33066907388e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
2.22044604925e-16
0.0
0.0
1.11022302463e-16
0.0
1.11022302463e-16
2.22044604925e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
2.22044604925e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
-1.11022302463e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
0.0
0.0
1.11022302463e-16
1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
1.11022302463e-16
-3.33066907388e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
1.11022302463e-16
0.0
0.0
-1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
0.0
1.11022302463e-16
0.0
0.0
2.22044604925e-16
1.11022302463e-16
0.0
0.0
0.0
2.22044604925e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
0.0
1.11022302463e-16
0.0
0.0
2.22044604925e-16
1.11022302463e-16
0.0
0.0
0.0
2.22044604925e-16
-1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
0.0
0.0
0.0
0.0
1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
3.33066907388e-16
1.11022302463e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
3.33066907388e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
2.22044604925e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
-1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
0.0
-4.4408920985e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-3.33066907388e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
0.0
-3.33066907388e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
0.0
-3.33066907388e-16
0.0
1.11022302463e-16
1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
1.11022302463e-16
-3.33066907388e-16
1.11022302463e-16
0.0
0.0
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
-2.22044604925e-16
0.0
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-2.22044604925e-16
0.0
-4.4408920985e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-3.33066907388e-16
-3.33066907388e-16
0.0
-4.4408920985e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
0.0
-2.22044604925e-16
-2.22044604925e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
0.0
-3.33066907388e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-3.33066907388e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-3.33066907388e-16
-3.33066907388e-16
0.0
-4.4408920985e-16
0.0
-1.11022302463e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
0.0
-3.33066907388e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
0.0
-3.33066907388e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
1.11022302463e-16
-3.33066907388e-16
1.11022302463e-16
0.0
0.0
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
-2.22044604925e-16
0.0
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-1.11022302463e-16
0.0
1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
0.0
-4.4408920985e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
0.0
-3.33066907388e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
-3.33066907388e-16
0.0
1.11022302463e-16
0.0
0.0
2.22044604925e-16
1.11022302463e-16
0.0
0.0
0.0
2.22044604925e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
0.0
0.0
2.22044604925e-16
1.11022302463e-16
0.0
0.0
0.0
2.22044604925e-16
-2.22044604925e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
0.0
-1.11022302463e-16
0.0
1.11022302463e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
0.0
2.22044604925e-16
1.11022302463e-16
0.0
0.0
0.0
2.22044604925e-16
-2.22044604925e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
2.22044604925e-16
0.0
0.0
0.0
0.0
2.22044604925e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
1.11022302463e-16
-3.33066907388e-16
1.11022302463e-16
0.0
0.0
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
-2.22044604925e-16
0.0
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
0.0
0.0
2.22044604925e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
0.0
1.11022302463e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
0.0
2.22044604925e-16
-1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
-1.11022302463e-16
3.33066907388e-16
2.22044604925e-16
1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
0.0
-4.4408920985e-16
1.11022302463e-16
0.0
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-3.33066907388e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
0.0
-3.33066907388e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
0.0
-3.33066907388e-16
4.4408920985e-16
3.33066907388e-16
2.22044604925e-16
-1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
2.22044604925e-16
3.33066907388e-16
1.11022302463e-16
3.33066907388e-16
1.11022302463e-16
2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
3.33066907388e-16
1.11022302463e-16
3.33066907388e-16
1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-4.4408920985e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-3.33066907388e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
0.0
-2.22044604925e-16
-2.22044604925e-16
-2.22044604925e-16
-3.33066907388e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
0.0
-3.33066907388e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
0.0
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-1.11022302463e-16
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
2.22044604925e-16
1.11022302463e-16
3.33066907388e-16
1.11022302463e-16
2.22044604925e-16
3.33066907388e-16
2.22044604925e-16
3.33066907388e-16
2.22044604925e-16
3.33066907388e-16
3.33066907388e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
2.22044604925e-16
3.33066907388e-16
1.11022302463e-16
4.4408920985e-16
1.11022302463e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
1.11022302463e-16
0.0
2.22044604925e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
2.22044604925e-16
0.0
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
-2.22044604925e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
1.11022302463e-16
2.22044604925e-16
0.0
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
1.11022302463e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
-2.22044604925e-16
0.0
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
2.22044604925e-16
1.11022302463e-16
1.11022302463e-16
1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-3.33066907388e-16
-1.11022302463e-16
-2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
-2.22044604925e-16
0.0
-3.33066907388e-16
0.0
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
0.0
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
-1.11022302463e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
0.0
0.0
1.11022302463e-16
2.22044604925e-16
0.0
3.33066907388e-16
0.0
0.0
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
0.0
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
1.11022302463e-16
-1.11022302463e-16
2.22044604925e-16
-1.11022302463e-16
-2.22044604925e-16
1.11022302463e-16
-2.22044604925e-16
3.33066907388e-16
0.0
-3.33066907388e-16

In [1943]:
for i in range(len(cip1)):
    print (cip1[i][1]/cip1[i][0])*np.log2(org1[inde[0][i]][3])-np.log2(org1[inde[1][i]][3])
    print b1[i]


0.0242960306378
0.0242960306378
-0.593614405608
-0.593614405608
1.00682722112
1.00682722112
-1.39951074432
-1.39951074432
-2.7647456932
-2.7647456932
-2.6801539874
-2.6801539874
-0.816833388596
-0.816833388596
-2.71431138142
-2.71431138142
0.949758826131
0.949758826131
-0.431945024915
-0.431945024915
-4.32758429026
-4.32758429026
0.254397200383
0.254397200383
-3.20206889585
-3.20206889585
1.22175010091
1.22175010091
-2.0481735664
-2.0481735664
0.386551000023
0.386551000023
-0.601940146656
-0.601940146656
-2.79234749798
-2.79234749798
-4.62964225429
-4.62964225429
-5.63250540252
-5.63250540252
-0.28050737905
-0.28050737905
-4.4519438523
-4.4519438523
-0.00635582016726
-0.00635582016726
-4.45262640471
-4.45262640471
-1.78119266794
-1.78119266794
-5.51660757361
-5.51660757361
0.709315296884
0.709315296884
0.575700058359
0.575700058359
-3.97414611725
-3.97414611725
-0.615952936147
-0.615952936147
0.98569076397
0.98569076397
-1.42228093728
-1.42228093728
-2.79050214098
-2.79050214098
-2.70522542328
-2.70522542328
-0.838355899755
-0.838355899755
-2.74109223337
-2.74109223337
0.928431088829
0.928431088829
-0.452087887923
-0.452087887923
-4.35722850374
-4.35722850374
0.231627007418
0.231627007418
-3.23228319037
-3.23228319037
1.20080152339
1.20080152339
-2.07495441836
-2.07495441836
0.361479564153
0.361479564153
-0.621497947253
-0.621497947253
-2.81959275777
-2.81959275777
-4.66085702213
-4.66085702213
-5.66291462797
-5.66291462797
-0.301549477013
-0.301549477013
-4.48121981468
-4.48121981468
-0.0285889802792
-0.0285889802792
-4.48384117255
-4.48384117255
-1.80561458475
-1.80561458475
-5.54956862217
-5.54956862217
0.689258148187
0.689258148187
0.554751480831
0.554751480831
-4.0036050544
-4.0036050544
1.56849825064
1.56849825064
-0.794425528942
-0.794425528942
-2.08030503974
-2.08030503974
-2.01391654281
-2.01391654281
-0.244903527486
-0.244903527486
-2.00264865652
-2.00264865652
1.51651285139
1.51651285139
0.103322665611
0.103322665611
-3.53983184005
-3.53983184005
0.859482415761
0.859482415761
-2.39916736006
-2.39916736006
1.77842849906
1.77842849906
-1.3365108415
-1.3365108415
1.05278844461
1.05278844461
-0.0822196504606
-0.0822196504606
-2.06834380096
-2.06834380096
-3.80015457493
-3.80015457493
-4.82442385683
-4.82442385683
0.278656190517
0.278656190517
-3.67397714681
-3.67397714681
0.584458517489
0.584458517489
-3.62313872535
-3.62313872535
-1.13221526078
-1.13221526078
-4.64071498911
-4.64071498911
1.24230525256
1.24230525256
1.13237845651
1.13237845651
-3.19131711986
-3.19131711986
-2.48416035934
-2.48416035934
-3.99164443806
-3.99164443806
-3.87442244649
-3.87442244649
-1.842050148
-1.842050148
-3.99000723545
-3.99000723545
-0.0661799535476
-0.0661799535476
-1.39144276128
-1.39144276128
-5.73967529849
-5.73967529849
-0.830252414636
-0.830252414636
-4.6413155004
-4.6413155004
0.223872455095
0.223872455095
-3.32386942043
-3.32386942043
-0.807717459067
-0.807717459067
-1.53356865416
-1.53356865416
-4.09016524543
-4.09016524543
-6.11654603117
-6.11654603117
-7.08103746903
-7.08103746903
-1.28283983579
-1.28283983579
-5.84649335732
-5.84649335732
-1.06542407634
-1.06542407634
-5.93953018159
-5.93953018159
-2.94452152965
-2.94452152965
-7.08669477857
-7.08669477857
-0.246099470388
-0.246099470388
-0.422177587458
-0.422177587458
-5.37741155669
-5.37741155669
-1.18169255619
-1.18169255619
-1.13920332744
-1.13920332744
0.505991835487
0.505991835487
-1.06829590373
-1.06829590373
2.26061278122
2.26061278122
0.806083710444
0.806083710444
-2.50557973633
-2.50557973633
1.6539079447
1.6539079447
-1.34502579281
-1.34502579281
2.50929998569
2.50929998569
-0.402158088709
-0.402158088709
1.92750165999
1.92750165999
0.600129247842
0.600129247842
-1.11778839928
-1.11778839928
-2.71110765777
-2.71110765777
-3.76348137624
-3.76348137624
1.01279049628
1.01279049628
-2.65257289531
-2.65257289531
1.36014759528
1.36014759528
-2.53409180819
-2.53409180819
-0.280163009528
-0.280163009528
-3.49074237029
-3.49074237029
1.94207582486
1.94207582486
1.86324994313
1.86324994313
-2.16352909179
-2.16352909179
0.0110612352343
0.0110612352343
1.49343356052
1.49343356052
0.160395788216
0.160395788216
3.23911838204
3.23911838204
1.730227889
1.730227889
-1.14551849242
-1.14551849242
2.69859266829
2.69859266829
0.0411904750177
0.0411904750177
3.47040993138
3.47040993138
0.826533603233
0.826533603233
3.07776622266
3.07776622266
1.49743106436
1.49743106436
0.132210085007
0.132210085007
-1.27899018902
-1.27899018902
-2.36832177765
-2.36832177765
1.97819111138
1.97819111138
-1.30940682214
-1.30940682214
2.38019352821
2.38019352821
-1.10197433944
-1.10197433944
0.840301953073
0.840301953073
-1.97850644175
-1.97850644175
2.86228747499
2.86228747499
2.82435988883
2.82435988883
-0.811968230654
-0.811968230654
1.48393807091
1.48393807091
0.148580377852
0.148580377852
3.22970882447
3.22970882447
1.72134108462
1.72134108462
-1.15859718566
-1.15859718566
2.68854671551
2.68854671551
0.0278602684533
0.0278602684533
3.46116765483
3.46116765483
0.814718192869
0.814718192869
3.06670498743
3.06670498743
1.48880238294
1.48880238294
0.120189783134
0.120189783134
-1.29276179315
-1.29276179315
-2.38173798554
-2.38173798554
1.96890757466
1.96890757466
-1.32232304713
-1.32232304713
2.37038450829
2.37038450829
-1.11574594357
-1.11574594357
0.829527278337
0.829527278337
-1.99304848528
-1.99304848528
2.85343848681
2.85343848681
2.81511761228
2.81511761228
-0.824965182054
-0.824965182054
-1.69791074447
-1.69791074447
1.75920005736
1.75920005736
0.332527249027
0.332527249027
-3.2025118871
-3.2025118871
1.11858324919
1.11858324919
-2.05536048494
-2.05536048494
2.01680126581
2.01680126581
-1.03177292946
-1.03177292946
1.33807500057
1.33807500057
0.140327538417
0.140327538417
-1.75832141646
-1.75832141646
-3.44496336619
-3.44496336619
-4.47839887283
-4.47839887283
0.518093121408
0.518093121408
-3.34084750383
-3.34084750383
0.837448482204
0.837448482204
-3.26794751661
-3.26794751661
-0.85431954824
-0.85431954824
-4.26565294352
-4.26565294352
1.47053449732
1.47053449732
1.37075122326
1.37075122326
-2.85610541661
-2.85610541661
3.11138236971
3.11138236971
1.60958832179
1.60958832179
-1.32306351586
-1.32306351586
2.56221750536
2.56221750536
-0.13976887579
-0.13976887579
3.34494478149
3.34494478149
0.666137815017
0.666137815017
2.92760846348
2.92760846348
1.38029555098
1.38029555098
-0.0309671330623
-0.0309671330623
-1.46594157131
-1.46594157131
-2.55044860814
-2.55044860814
1.85216584921
1.85216584921
-1.48474631733
-1.48474631733
2.24703476063
2.24703476063
-1.28892572173
-1.28892572173
0.69403429128
0.69403429128
-2.17591664263
-2.17591664263
2.74216126765
2.74216126765
2.69889473894
2.69889473894
-0.988403597692
-0.988403597692
-1.32893947182
-1.32893947182
-5.64768932532
-5.64768932532
-0.759596522201
-0.759596522201
-4.54756056621
-4.54756056621
0.288875876135
0.288875876135
-3.24076845603
-3.24076845603
-0.729920811545
-0.729920811545
-1.47288081294
-1.47288081294
-4.00562322384
-4.00562322384
-6.01968662896
-6.01968662896
-6.9866776643
-6.9866776643
-1.2175462209
-1.2175462209
-5.75565006705
-5.75565006705
-0.996434596464
-0.996434596464
-5.84267077938
-5.84267077938
-2.86874033932
-2.86874033932
-6.98441666854
-6.98441666854
-0.183862152371
-0.183862152371
-0.357174166418
-0.357174166418
-5.28600049585
-5.28600049585
-3.69189161208
-3.69189161208
0.742682880722
0.742682880722
-2.55415135848
-2.55415135848
1.67097292682
1.67097292682
-1.47388302191
-1.47388302191
0.924184701248
0.924184701248
-0.182541242797
-0.182541242797
-2.20809815792
-2.20809815792
-3.96027049383
-3.96027049383
-4.98040775201
-4.98040775201
0.170720905906
0.170720905906
-3.82414797758
-3.82414797758
0.470413688466
0.470413688466
-3.78325464425
-3.78325464425
-1.25748730095
-1.25748730095
-4.80978844193
-4.80978844193
1.13942225786
1.13942225786
1.02492288427
1.02492288427
-3.34242651832
-3.34242651832
3.57848368421
3.57848368421
1.20873816922
1.20873816922
4.27990966603
4.27990966603
1.86140542309
1.86140542309
4.04658239444
4.04658239444
2.25318807804
2.25318807804
1.18502765897
1.18502765897
-0.0727819751443
-0.0727819751443
-1.19324151768
-1.19324151768
2.79130468412
2.79130468412
-0.178118373097
-0.178118373097
3.23933239753
3.23933239753
0.104233874437
0.104233874437
1.7840192603
1.7840192603
-0.704818048078
-0.704818048078
3.63734041242
3.63734041242
3.63385962348
3.63385962348
0.326390771191
0.326390771191
-3.53963441175
-3.53963441175
0.98770467656
0.98770467656
-2.34737936458
-2.34737936458
0.106443444282
0.106443444282
-0.820447285575
-0.820447285575
-3.09674183601
-3.09674183601
-4.97838543627
-4.97838543627
-5.97224876045
-5.97224876045
-0.515597649047
-0.515597649047
-4.77902596708
-4.77902596708
-0.254753086579
-0.254753086579
-4.80136958669
-4.80136958669
-2.05404355126
-2.05404355126
-5.88486086368
-5.88486086368
0.485229252291
0.485229252291
0.341654634006
0.341654634006
-4.30327249525
-4.30327249525
3.44185120204
3.44185120204
0.790023864013
0.790023864013
3.04358689233
3.04358689233
1.47076835025
1.47076835025
0.095067228922
0.095067228922
-1.32154458705
-1.32154458705
-2.40977799766
-2.40977799766
1.94950488771
1.94950488771
-1.34931808985
-1.34931808985
2.34988355603
2.34988355603
-1.14452873747
-1.14452873747
0.807008097618
0.807008097618
-2.02344150541
-2.02344150541
2.83494401073
2.83494401073
2.79580115948
2.79580115948
-0.852128943796
-0.852128943796
-3.61007000223
-3.61007000223
-1.07564991862
-1.07564991862
-1.74257820767
-1.74257820767
-4.38132884309
-4.38132884309
-6.45013081525
-6.45013081525
-7.40601361353
-7.40601361353
-1.50771172149
-1.50771172149
-6.1593585896
-6.1593585896
-1.30302455934
-1.30302455934
-6.27311496567
-6.27311496567
-3.20551273377
-3.20551273377
-7.43894164847
-7.43894164847
-0.460445438032
-0.460445438032
-0.646050042554
-0.646050042554
-5.69223219666
-5.69223219666
2.30399008346
2.30399008346
0.893821462003
0.893821462003
-0.708656470883
-0.708656470883
-2.24236710405
-2.24236710405
-3.30683735293
-3.30683735293
1.3287718517
1.3287718517
-2.2129466617
-2.2129466617
1.69401468779
1.69401468779
-2.06535125447
-2.06535125447
0.0865718278443
0.0865718278443
-2.99577856881
-2.99577856881
2.24326656364
2.24326656364
2.17782693697
2.17782693697
-1.72115519421
-1.72115519421
-0.903482005596
-0.903482005596
-3.21241448066
-3.21241448066
-5.11091104902
-5.11091104902
-6.10135435739
-6.10135435739
-0.604934111213
-0.604934111213
-4.90332017531
-4.90332017531
-0.349146329622
-0.349146329622
-4.93389519944
-4.93389519944
-2.15772939335
-2.15772939335
-6.02480049672
-6.02480049672
0.400074496865
0.400074496865
0.252715222784
0.252715222784
-4.42834354228
-4.42834354228
-1.95380660003
-1.95380660003
-3.66892983479
-3.66892983479
-4.69658556159
-4.69658556159
0.367115725165
0.367115725165
-3.55090301165
-3.55090301165
0.677925195607
0.677925195607
-3.49191398521
-3.49191398521
-1.02954719984
-1.02954719984
-4.50214900477
-4.50214900477
1.32662412814
1.32662412814
1.22044483767
1.22044483767
-3.06747377135
-3.06747377135
-1.43046267052
-1.43046267052
-2.51588529188
-2.51588529188
1.87608242965
1.87608242965
-1.45147107497
-1.45147107497
2.27230510978
2.27230510978
-1.25344682094
-1.25344682094
0.721792395009
0.721792395009
-2.13845290823
-2.13845290823
2.76495834859
2.76495834859
2.72270502347
2.72270502347
-0.954920385068
-0.954920385068
-1.12233778705
-1.12233778705
2.84036753344
2.84036753344
-0.109857017527
-0.109857017527
3.29117238926
3.29117238926
0.177015849581
0.177015849581
1.84096267091
1.84096267091
-0.627964354044
-0.627964354044
3.68410670283
3.68410670283
3.68270441568
3.68270441568
0.395078760233
0.395078760233
3.61698519859
3.61698519859
0.970654516591
0.970654516591
4.11174954489
4.11174954489
1.32908443894
1.32908443894
2.74232203357
2.74232203357
0.588555624928
0.588555624928
4.42437205174
4.42437205174
4.45587044676
4.45587044676
1.48234349144
1.48234349144
-4.06167271623
-4.06167271623
0.290029335056
0.290029335056
-4.03650949804
-4.03650949804
-1.45562969951
-1.45562969951
-5.07721139943
-5.07721139943
0.976692543298
0.976692543298
0.854960737947
0.854960737947
-3.58143578658
-3.58143578658
3.37460153936
3.37460153936
0.294148166282
0.294148166282
1.9326050534
1.9326050534
-0.504279180465
-0.504279180465
3.75937044676
3.75937044676
3.76131321489
3.76131321489
0.50562238412
0.50562238412
-4.44370300156
-4.44370300156
-1.77421114527
-1.77421114527
-5.5071849591
-5.5071849591
0.715049058056
0.715049058056
0.581688653361
0.581688653361
-3.96572465553
-3.96572465553
1.70246840517
1.70246840517
-0.814883188218
-0.814883188218
3.57036460374
3.57036460374
3.56390711219
3.56390711219
0.228020052191
0.228020052191
-3.11262026653
-3.11262026653
2.17216714758
2.17216714758
2.10356754686
2.10356754686
-1.82558246155
-1.82558246155
4.06622969274
4.06622969274
4.08181064959
4.08181064959
0.95632190166
0.95632190166
-0.165140362831
-0.165140362831
-5.01595295955
-5.01595295955
-4.78372432432
-4.78372432432

In [877]:
np.log2(dd)


Out[877]:
array([-1.85333035, -1.05246742, -0.8214027 ])

In [1508]:
dd=np.array([p[3] for p in org1])
dd


Out[1508]:
array([ 0.685,  0.255,  0.605,  0.515,  0.68 ,  0.435,  0.51 ,  0.365,
        0.575,  0.485])

In [879]:
print M1


[[ 0.76377953 -1.          0.        ]
 [ 1.94        0.         -1.        ]
 [ 0.          2.54       -1.        ]]

In [628]:
tau_snok[0]/float(tau_snok[1])


Out[628]:
0.77300613496932513

In [631]:
tau_snok[0]/float(tau_snok[2])


Out[631]:
0.78749999999999998

In [632]:
tau_snok[1]/float(tau_snok[2])


Out[632]:
1.01875

In [1509]:
viaCellDist=np.dot(M1,np.log2(dd).T)
viaCellDist


Out[1509]:
array([ 1.40102007,  0.29427461,  0.49854699,  0.0838889 ,  0.69438792,
        0.50242643,  0.91287269,  0.27938584,  0.60117694, -0.7636385 ,
       -0.62836045, -1.07665758, -0.5497179 , -0.6495234 , -0.4163002 ,
       -0.99531275, -0.48632815,  0.18508056, -0.23893474,  0.34832098,
        0.18199408,  0.54314305, -0.07519094,  0.29867087, -0.4295401 ,
        0.14399204, -0.00719939,  0.3248429 , -0.28454436,  0.12006166,
        0.60445902,  0.41915893,  0.8167948 ,  0.18724558,  0.52256776,
       -0.14052535,  0.17100525, -0.4320772 , -0.00580551,  0.33314988,
       -0.27657783,  0.12685828, -0.59607403, -0.1457189 ,  0.36281979])

In [1510]:
b1


Out[1510]:
array([ 1.40102007,  0.29427461,  0.49854699,  0.0838889 ,  0.69438792,
        0.50242643,  0.91287269,  0.27938584,  0.60117694, -0.7636385 ,
       -0.62836045, -1.07665758, -0.5497179 , -0.6495234 , -0.4163002 ,
       -0.99531275, -0.48632815,  0.18508056, -0.23893474,  0.34832098,
        0.18199408,  0.54314305, -0.07519094,  0.29867087, -0.4295401 ,
        0.14399204, -0.00719939,  0.3248429 , -0.28454436,  0.12006166,
        0.60445902,  0.41915893,  0.8167948 ,  0.18724558,  0.52256776,
       -0.14052535,  0.17100525, -0.4320772 , -0.00580551,  0.33314988,
       -0.27657783,  0.12685828, -0.59607403, -0.1457189 ,  0.36281979])

In [797]:
print b1


[-0.13128513  0.4296232   0.49356518]

In [613]:
cip[1][1]/cip[1][0]


Out[613]:
0.78750000000000009

In [621]:
tau_snok[2]/float(tau_snok[1])


Out[621]:
0.98159509202453987

In [618]:
tau_snok[1]


Out[618]:
163

In [612]:
-cip[1][1]/cip[1][0]+cip[0][1]/cip[0][0]*cip[2][1]/cip[2][0]


Out[612]:
-1.1102230246251565e-16

In [614]:
cip


Out[614]:
array([[ 0.23809524,  0.18404908],
       [ 0.23809524,  0.1875    ],
       [ 0.18404908,  0.1875    ]])

In [583]:
cip2=create_pairs([p[1] for p in same_org])
cip2


Out[583]:
array([[[ 0.15384615,  0.15384615],
        [ 0.15384615,  0.20979021],
        [ 0.15384615,  0.44776119]],

       [[ 0.20979021,  0.15384615],
        [ 0.20979021,  0.20979021],
        [ 0.20979021,  0.44776119]],

       [[ 0.44776119,  0.15384615],
        [ 0.44776119,  0.20979021],
        [ 0.44776119,  0.44776119]]])

In [578]:
inde


Out[578]:
(array([0, 0, 1]), array([1, 2, 2]))

In [560]:
yy=np.linalg.solve(M[0:6],b.T[0:6])

In [ ]:
Mt

In [1927]:
print np.linalg.pinv(Mt)


[[ 0.44711538  0.375       0.375     ]
 [-0.5         0.41935484  0.41935484]
 [ 0.59615385 -0.5         0.5       ]]

In [1928]:
np.linalg.norm(np.linalg.pinv(Mt))*np.linalg.norm(Mt)


Out[1928]:
3.8000038116125108

In [497]:
inde


Out[497]:
(array([0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]),
 array([1, 2, 3, 4, 5, 2, 3, 4, 5, 3, 4, 5, 4, 5, 5]))

In [489]:
.43*.63-.55*.49


Out[489]:
0.0013999999999999568

In [1892]:
print D1_snok
print C1_snok
print tau1_snok


36
41
[208 186 156]

In [2116]:
M1.shape


Out[2116]:
(8385, 130)

In [2142]:
np.linalg.norm(np.linalg.pinv(Mt))*np.linalg.norm(Mt)


Out[2142]:
1404.5085861661685

In [2249]:
Mt=M1;
bt=b1;
Mt=np.append(Mt,[np.zeros(n)], axis=0);
Mt[-1][-1]=1
bt=np.append(bt,[0], axis=0);
bt[-1]=(np.log2(Ri1[-1])-np.log2(Ki1[-1])-org1[-1][2])
#bt[-1]=np.mean(np.log2(np.array(Rip1)[:,0])-np.log2(Kip1[:,0])-np.array(cip1)[:,0])
bt[-1]=np.log2(Ri1[-1])-np.log2(Ki1[-1])-cip1[-1][0]
bt[-1]=np.log2(Ri1[-1])-np.log2(Ki1[-1])-aa1[-1]
#bt[-1]=np.log2(Ahat[-25])
bt[-1]


Out[2249]:
-2.1995569268971438

In [2238]:
np.abs(np.abs(np.log2(Ri1[-1])-np.log2(Ki1[-1])-org1[-1][2])-abs(bt[-1]))/np.abs(np.log2(Ri1[-1])-np.log2(Ki1[-1])-org1[-1][2])


Out[2238]:
0.028334480029993895

In [2250]:
print np.log2(Ri1[-1])-np.log2(Ki1[-1])-org1[-1][2]
print np.log2(Ri1[-1])-np.log2(Ki1[-1])-org1[-1][1]
print np.log2(Ri1[-1])-np.log2(Ki1[-1])-cip1[-1][0]
print np.mean(np.log2(Ri1)-np.log2(Ki1)-np.array(org1)[:,1])


-2.13895086629
-2.1995569269
-2.20961981998
-2.49371103151

In [2240]:
log2Ahat=np.linalg.lstsq(Mt,bt.T)[0]
log2Ahat


Out[2240]:
array([ -1.16284802,  -2.36564177,  -1.05220777,  -2.06773073,
        -1.40065404,  -2.35363131,  -2.19614945,  -2.19834727,
        -1.24524446,  -1.268311  ,  -2.58907996,  -1.86369207,
        -2.528346  ,  -3.51509855,  -1.32066276,  -1.34286054,
        -1.12897437,  -2.46103614,  -1.15104639,  -2.88600213,
        -5.85454486,  -1.15925912,  -3.12999509,  -5.96039282,
        -5.35157626,  -2.25407526,  -6.38042136,  -1.8288698 ,
       -11.36096697,  -1.27842116,  -1.33597523,  -2.55897227,
        -2.16201715,  -3.76812388,  -3.9967352 ,  -1.89128576,
        -5.49146026,  -1.08043689,  -5.23413352,  -1.09399313,
        -2.04471817,  -1.11800707,  -1.89593104,  -6.68038374,
        -1.82731868,  -2.65336479,  -1.8736463 ,  -1.3686219 ,
        -2.27463559,  -1.72188249,  -1.20007855,  -2.209678  ,
        -4.04278334,  -1.69677531,  -6.56176623,  -1.36683063,
        -4.54535794,  -5.85398543,  -1.91666956,  -1.28482926,
        -1.11639461,  -1.38128469,  -2.54008099,  -2.76658941,
        -1.30606577,  -4.71095226,  -2.61901063,  -3.94586977,
        -1.64000963,  -1.45990086,  -1.67258415,  -1.65386084,
        -1.14523177,  -7.97571947,  -2.00201896,  -1.04559378,
        -1.12026089,  -2.16672473,  -7.41315113,  -1.69887509,
        -1.18862092,  -2.06653383,  -2.62551548,  -1.14610319,
        -1.4419606 ,  -2.16089814,  -2.66992554,  -1.83820304,
        -1.77585263,  -1.65907649,  -2.03657516,  -4.16494454,
        -1.45553937,  -2.10779403,  -1.14774053,  -1.28930364,
        -3.3271733 ,  -1.1270362 ,  -3.77862338,  -1.10402682,
        -1.42786996,  -1.97062866,  -1.66344927,  -1.17877677,
        -2.56996265,  -1.795483  ,  -3.69020494,  -1.08036757,
        -3.42759436,  -2.00183228,  -2.84140359,  -2.66140173,
        -2.63132235,  -1.96914884,  -2.69033385,  -1.17439767,
        -1.87975216,  -3.06192522,  -6.62154389,  -1.21750231,
        -1.71125247,  -1.1201544 ,  -1.6952005 ,  -1.78145826,
        -3.3121694 ,  -1.23481519,  -1.89756184,  -2.38250239,
        -2.39475066,  -2.19955693])

In [2241]:
log2Ahat=np.linalg.solve(Mt,bt.T)
log2Ahat


---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-2241-e4aa116b1efc> in <module>()
----> 1 log2Ahat=np.linalg.solve(Mt,bt.T)
      2 log2Ahat

/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/linalg/linalg.pyc in solve(a, b)
    353     a, _ = _makearray(a)
    354     _assertRankAtLeast2(a)
--> 355     _assertNdSquareness(a)
    356     b, wrap = _makearray(b)
    357     t, result_t = _commonType(a, b)

/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/linalg/linalg.pyc in _assertNdSquareness(*arrays)
    210     for a in arrays:
    211         if max(a.shape[-2:]) != min(a.shape[-2:]):
--> 212             raise LinAlgError('Last 2 dimensions of the array must be square')
    213 
    214 def _assertFinite(*arrays):

LinAlgError: Last 2 dimensions of the array must be square

In [2242]:
Ahat=2**log2Ahat

In [2243]:
np.abs(np.log2(dd))-np.abs(np.log2(Ri1)-np.log2(Ki1)-np.array(org1)[:,2])


Out[2243]:
array([  2.22044605e-16,  -4.44089210e-16,   2.22044605e-16,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,  -2.22044605e-16,
        -2.22044605e-16,   4.44089210e-16,   0.00000000e+00,
         0.00000000e+00,   4.44089210e-16,   2.22044605e-16,
         2.22044605e-16,  -2.22044605e-16,   4.44089210e-16,
        -2.22044605e-16,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         2.22044605e-16,  -1.77635684e-15,  -4.44089210e-16,
        -2.22044605e-16,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   8.88178420e-16,
         2.22044605e-16,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   8.88178420e-16,   0.00000000e+00,
         0.00000000e+00,  -2.22044605e-16,   0.00000000e+00,
         0.00000000e+00,  -2.22044605e-16,   2.22044605e-16,
        -4.44089210e-16,   0.00000000e+00,   2.22044605e-16,
         0.00000000e+00,  -2.22044605e-16,   0.00000000e+00,
         0.00000000e+00,   2.22044605e-16,  -2.22044605e-16,
         2.22044605e-16,   0.00000000e+00,   0.00000000e+00,
         4.44089210e-16,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   4.44089210e-16,  -2.22044605e-16,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
        -2.22044605e-16,   0.00000000e+00,  -4.44089210e-16,
         2.22044605e-16,   4.44089210e-16,   0.00000000e+00,
         8.88178420e-16,  -4.44089210e-16,   0.00000000e+00,
         0.00000000e+00,   4.44089210e-16,   0.00000000e+00,
         0.00000000e+00,   4.44089210e-16,   0.00000000e+00,
         2.22044605e-16,  -2.22044605e-16,   2.22044605e-16,
        -2.22044605e-16,   0.00000000e+00,   0.00000000e+00,
         4.44089210e-16,   0.00000000e+00,   2.22044605e-16,
        -4.44089210e-16,   0.00000000e+00,   4.44089210e-16,
        -2.22044605e-16,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   2.22044605e-16,   0.00000000e+00,
         2.22044605e-16,   0.00000000e+00,   2.22044605e-16,
         0.00000000e+00,  -2.22044605e-16,   0.00000000e+00,
         0.00000000e+00,   4.44089210e-16,   0.00000000e+00,
        -4.44089210e-16,   4.44089210e-16,  -2.22044605e-16,
         0.00000000e+00,   8.88178420e-16,   4.44089210e-16,
        -2.22044605e-16,  -2.22044605e-16,   2.22044605e-16,
         4.44089210e-16,  -4.44089210e-16,   0.00000000e+00,
        -2.22044605e-16,   0.00000000e+00,   0.00000000e+00,
        -4.44089210e-16])

In [2244]:
np.abs(np.abs(np.log2(Ri1)-np.log2(Ki1)-np.array(org1)[:,2])-abs(np.log2(Ahat)))/np.abs(np.log2(Ri1)-np.log2(Ki1)-np.array(org1)[:,2])


Out[2244]:
array([ 0.04029784,  0.02251222,  0.04177576,  0.02684346,  0.0370191 ,
        0.01865524,  0.02001975,  0.02044886,  0.05319322,  0.04191309,
        0.01989504,  0.03382429,  0.02160757,  0.02029801,  0.05591128,
        0.04101712,  0.04758506,  0.01686472,  0.05419266,  0.01857561,
        0.01039658,  0.04119999,  0.01494302,  0.01307444,  0.01296462,
        0.03346627,  0.00801961,  0.02825297,  0.00602387,  0.04049195,
        0.03178215,  0.0163444 ,  0.03342598,  0.01243913,  0.01511827,
        0.02301749,  0.00736717,  0.05207202,  0.0084879 ,  0.06828688,
        0.02162958,  0.0470625 ,  0.03973955,  0.01050447,  0.027705  ,
        0.01854439,  0.03715079,  0.05813051,  0.02919124,  0.03069678,
        0.0466386 ,  0.01866401,  0.01723406,  0.02475167,  0.00973938,
        0.03574389,  0.01074242,  0.00832103,  0.03899896,  0.05631104,
        0.05117272,  0.03570539,  0.0215056 ,  0.02593211,  0.03212416,
        0.01455136,  0.0163745 ,  0.01452604,  0.03987853,  0.04303972,
        0.02767678,  0.02866251,  0.0558774 ,  0.00701445,  0.02509965,
        0.04186705,  0.04001592,  0.01936056,  0.00577337,  0.02971146,
        0.04152841,  0.02746931,  0.01579152,  0.04978022,  0.03349191,
        0.02128931,  0.0205505 ,  0.03189551,  0.03285951,  0.02620943,
        0.02062333,  0.01660449,  0.0394551 ,  0.03036239,  0.0455481 ,
        0.0450299 ,  0.01926706,  0.04898674,  0.01091986,  0.07135332,
        0.04296679,  0.03600855,  0.02745171,  0.03973166,  0.01843175,
        0.0245971 ,  0.01219072,  0.06099823,  0.01983579,  0.03181706,
        0.02452866,  0.01653483,  0.01936992,  0.02278081,  0.01515397,
        0.04994576,  0.02704823,  0.02195811,  0.01059881,  0.04781248,
        0.03648731,  0.04131638,  0.02456886,  0.03782941,  0.0149478 ,
        0.06643372,  0.03068723,  0.02176975,  0.02697124,  0.02833448])

In [2245]:
dd=np.array([p[3] for p in org1])
#aun=np.array(Ahat)/(np.mean(Ahat)/np.mean(dd));
#aun=np.array(Ahat)/(Ahat[0]/dd[0]);
print np.abs((dd-Ahat)/dd)
print Ahat


[ 0.03074045  0.03545753  0.02882317  0.0367742   0.03406367  0.02943513
  0.02943513  0.03007363  0.0426576   0.03474663  0.03440176  0.04138436
  0.03638813  0.04731582  0.04731582  0.03601009  0.03492167  0.02789526
  0.04018491  0.035824    0.04089609  0.03129557  0.0314375   0.05192249
  0.04636648  0.04933609  0.03457334  0.03423188  0.04605844  0.03389711
  0.02812159  0.02812159  0.04731582  0.03158072  0.0404192   0.02906487
  0.02745335  0.03638813  0.03007363  0.04731582  0.02956065  0.03423188
  0.04898748  0.04699509  0.03356881  0.03293094  0.04545448  0.05078156
  0.04373401  0.03492167  0.03638813  0.02767254  0.04636648  0.02800797
  0.0429217   0.0321669   0.03293094  0.03293094  0.04864377  0.04636648
  0.03697032  0.03246822  0.03638813  0.04731582  0.02778345  0.04575447
  0.02882317  0.03840398  0.0426576   0.04089609  0.03074045  0.0314375
  0.04113877  0.03777616  0.03340704  0.02870382  0.02943513  0.02812159
  0.02906487  0.03340704  0.03231686  0.03757143  0.02789526  0.03697032
  0.03187112  0.03074045  0.03658015  0.03861792  0.03840398  0.02894352
  0.02812159  0.04605844  0.03757143  0.04213901  0.03406367  0.03777616
  0.0426576   0.035824    0.02789526  0.04968969  0.03995331  0.04636648
  0.03033686  0.03074045  0.03172525  0.02943513  0.03033686  0.04213901
  0.0451584   0.04188442  0.04605844  0.02956065  0.03406367  0.02994372
  0.02745335  0.03798314  0.03373216  0.04457765  0.04699509  0.03777616
  0.04089609  0.03033686  0.02778345  0.04401165  0.03324682  0.05192249
  0.03840398  0.03457334  0.0426576   0.04113877]
[  4.46629972e-01   1.94030887e-01   4.82229639e-01   2.38534405e-01
   3.78757396e-01   1.95652940e-01   2.18219292e-01   2.17887105e-01
   4.21836412e-01   4.15145509e-01   1.66191677e-01   2.74772197e-01
   1.73337295e-01   8.74681416e-02   4.00350979e-01   3.94238195e-01
   4.57240667e-01   1.81616082e-01   4.50298509e-01   1.35277881e-01
   1.72824929e-02   4.47742408e-01   1.14229320e-01   1.60599051e-02
   2.44914796e-02   2.09631111e-01   1.20033311e-02   2.81485048e-01
   3.80196541e-04   4.12246411e-01   3.96124206e-01   1.69696384e-01
   2.23443634e-01   7.33975706e-02   6.26415970e-02   2.69566709e-01
   2.22282740e-02   4.72885599e-01   2.65686085e-02   4.68462953e-01
   2.42369796e-01   4.60729835e-01   2.68700137e-01   9.74998866e-03
   2.81787851e-01   1.58948930e-01   2.72882862e-01   3.87260994e-01
   2.06664774e-01   3.03152897e-01   4.35251584e-01   2.16182552e-01
   6.06737655e-02   3.08474831e-01   1.05855040e-02   3.87742120e-01
   4.28263358e-02   1.72891958e-02   2.64865242e-01   4.10419375e-01
   4.61245067e-01   3.83876810e-01   1.71933075e-01   1.46951358e-01
   4.04422237e-01   3.81822980e-02   1.62779324e-01   6.48895617e-02
   3.20854333e-01   3.63518109e-01   3.13690958e-01   3.17788575e-01
   4.52117052e-01   3.97254844e-03   2.49650386e-01   4.84445482e-01
   4.60010632e-01   2.22715716e-01   5.86704135e-03   3.08026187e-01
   4.38722036e-01   2.38732382e-01   1.62047034e-01   4.51844044e-01
   3.68066768e-01   2.23617013e-01   1.57134782e-01   2.79669913e-01
   2.92021677e-01   3.16641774e-01   2.43741674e-01   5.57476753e-02
   3.64618742e-01   2.32001490e-01   4.51331529e-01   4.09148470e-01
   9.96370899e-02   4.57855353e-01   7.28653443e-02   4.65216180e-01
   3.71679246e-01   2.55141829e-01   3.15683494e-01   4.41725869e-01
   1.68408557e-01   2.88075126e-01   7.74707255e-02   4.72908320e-01
   9.29375638e-02   2.49682692e-01   1.39525083e-01   1.58065922e-01
   1.61396104e-01   2.55403669e-01   1.54927608e-01   4.43068703e-01
   2.71730393e-01   1.19748109e-01   1.01558592e-02   4.30026564e-01
   3.05394826e-01   4.60044588e-01   3.08811738e-01   2.90889220e-01
   1.00678714e-01   4.24896926e-01   2.68396574e-01   1.91776468e-01
   1.90155206e-01   2.17704491e-01]

In [2246]:
np.log2(dd)


Out[2246]:
array([ -1.11780298,  -2.31355844,  -1.01001368,  -2.01367667,
        -1.35065404,  -2.31052786,  -2.153046  ,  -2.15429441,
        -1.18235138,  -1.2172906 ,  -2.53857491,  -1.80271646,
        -2.47487007,  -3.44516848,  -1.25073269,  -1.28995049,
        -1.07769232,  -2.42021981,  -1.0918748 ,  -2.83337055,
        -5.79430389,  -1.11338756,  -3.08391214,  -5.88346975,
        -5.28308311,  -2.18108256,  -6.32965994,  -1.77861855,
       -11.29293976,  -1.22866992,  -1.29482297,  -2.51782001,
        -2.09208708,  -3.72182758,  -3.93721139,  -1.84873257,
        -5.45129962,  -1.02696095,  -5.19008066,  -1.02406306,
        -2.00142813,  -1.06775581,  -1.82346727,  -6.6109393 ,
        -1.7780576 ,  -2.60505561,  -1.80653221,  -1.29343393,
        -2.21011946,  -1.67060043,  -1.14660261,  -2.16919217,
        -3.97429019,  -1.65579171,  -6.49847509,  -1.31966082,
        -4.49704876,  -5.80567625,  -1.84472712,  -1.21633611,
        -1.06204679,  -1.33366564,  -2.48660506,  -2.69665934,
        -1.26541536,  -4.64338469,  -2.57681654,  -3.88937259,
        -1.57711655,  -1.3996599 ,  -1.6275391 ,  -1.60777789,
        -1.0846257 ,  -7.92016392,  -1.95299935,  -1.00357697,
        -1.07715744,  -2.12557247,  -7.37059793,  -1.64985548,
        -1.14122756,  -2.01128521,  -2.58469915,  -1.09175536,
        -1.39523162,  -2.1158531 ,  -2.6161621 ,  -1.78138486,
        -1.71935545,  -1.61670361,  -1.9954229 ,  -4.09691733,
        -1.40029076,  -2.04568222,  -1.09774053,  -1.23374808,
        -3.26428022,  -1.07440463,  -3.73780705,  -1.03049741,
        -1.36904643,  -1.90213551,  -1.61900482,  -1.13373173,
        -2.52345102,  -1.75237955,  -3.6457605 ,  -1.01825577,
        -3.3609277 ,  -1.94010388,  -2.77337638,  -2.61811168,
        -2.58132235,  -1.92528919,  -2.6501732 ,  -1.11853175,
        -1.83024721,  -2.99613575,  -6.55209944,  -1.16194676,
        -1.65101151,  -1.07570995,  -1.6545501 ,  -1.7165232 ,
        -3.26338891,  -1.15789211,  -1.84106467,  -2.33174097,
        -2.33185757,  -2.13895087])

In [2247]:
aa=np.array(org1)[:,2];
np.abs(aa-np.abs(np.log2(Ri1)-np.log2(Ki1)-np.log2(Ahat)))/aa


Out[2247]:
array([ 0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765,
        0.29411765,  0.29411765,  0.29411765,  0.29411765,  0.29411765])

In [2251]:
aa1=(np.log2(Ri1)-np.log2(Ki1)-np.log2(Ahat))

In [2252]:
aa2=(np.log2(Ri1)-np.log2(Ki1)-np.log2(Ahat))

In [2254]:
k1=aa1/np.array(org1)[:,1]
k2=aa2/np.array(org1)[:,1]
print k1


[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.]

In [ ]:


In [2100]:
aa1/aa2


Out[2100]:
array([ 1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756])

In [2092]:
(np.log2(Ri1)-np.log2(Ki1)-np.log2(Ahat))/np.array(org1)[:,1]


Out[2092]:
array([ 1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756,
        1.05583756,  1.05583756,  1.05583756,  1.05583756,  1.05583756])

In [2078]:
D1_snok/float(C1_snok)


Out[2078]:
0.5660377358490566

In [ ]:


In [2107]:
np.log2(aa1/(np.array(org1)[:,1]))-np.log2(aa2/(np.array(org1)[:,1]))


Out[2107]:
array([ 0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879,
        0.0783879,  0.0783879,  0.0783879,  0.0783879,  0.0783879])

In [2102]:
((aa1-aa2)/(np.array(org1)[:,1]))/()


Out[2102]:
array([ 1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006,
        1.00000006,  1.00000006,  1.00000006,  1.00000006,  1.00000006])

In [2043]:
np.log2(Ri1)-np.log2(Ki1)-np.log2(Ahat)


Out[2043]:
array([ 0.27155172,  0.19364754,  0.21283784,  0.22079439,  0.29166667,
        0.30681818,  0.19207317,  0.24230769,  0.19769874,  0.21774194,
        0.2473822 ,  0.19052419,  0.3375    ,  0.22183099,  0.31711409,
        0.31711409,  0.27155172,  0.32586207,  0.2027897 ,  0.28810976,
        0.20366379,  0.20723684,  0.23275862,  0.21188341,  0.30288462,
        0.328125  ,  0.20633188,  0.28810976,  0.25132979,  0.25      ])

In [2044]:
np.array(org1)[:,2]


Out[2044]:
array([ 0.1954023 ,  0.13934426,  0.15315315,  0.1588785 ,  0.20987654,
        0.22077922,  0.13821138,  0.17435897,  0.14225941,  0.15668203,
        0.17801047,  0.13709677,  0.24285714,  0.15962441,  0.22818792,
        0.22818792,  0.1954023 ,  0.23448276,  0.14592275,  0.20731707,
        0.14655172,  0.14912281,  0.16748768,  0.15246637,  0.21794872,
        0.23611111,  0.14847162,  0.20731707,  0.18085106,  0.17989418])

In [170]:
def null(a, rtol=1e-5):
    u, s, v = np.linalg.svd(a)
    rank = (s > rtol*s[0]).sum()
    return rank, v[rank:].T.copy()

In [594]:
def create_pairs(arr):
    l=len(arr)
    arr=cartesian([arr,arr]).reshape((l,l,2))
    inde=np.triu_indices(l,k=1)
    return arr[inde]

In [74]:
l1=[x.l for x in c.pop.values()]/np.mean([x.l for x in c.pop.values()])

In [75]:
arr1=cartesian([[x.C for x in c.pop.values()],[x.C for x in c.pop.values()]])
arr1=arr1.reshape((len(sel[0]),len(sel[0]),2))
arrC=arr1[inde]

In [160]:
def A_snok(R,K,l,cik):
    return ((K[1]*l[1])/float(R[1]))*(R[0]/(float(K[0])*l[0]))**(cik[1]/float(cik[0]))

In [212]:
1**.2


Out[212]:
1.0

In [68]:
print c.pop.keys()
print c.distribution


['NC_009614.1', 'NC_002655.2', 'NC_000913.3']
[ 0.43919776  0.40942012  0.15138212]

In [21]:
c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')

In [708]:
c.write_reads()

In [731]:
c.collect()


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-731-b901a14963f9> in <module>()
----> 1 c.collect()

<ipython-input-712-2a9bc4d0021b> in collect(self)
    112     def collect(self):
    113         args=ptr_pipeline.parse_args(['collect'])
--> 114         ptr_pipeline.main(args,self.conf)
    115 
    116     def create_dirs(self,from_init=False):

/Users/hedani/Documents/GitRepos/PTRloc/PTR-Pipeline/ptr_pipeline.pyc in main(args, config)
    399 
    400     if(args.subparser_name=='collect'):
--> 401         process = subprocess.Popen(generate_collect_command(config), shell=True)
    402         process.wait()
    403     # prov = provider.get_provider(args['provider_name'], args['project_dir'])

/Users/hedani/Documents/GitRepos/PTRloc/PTR-Pipeline/ptr_pipeline.pyc in generate_collect_command(config)
    339 def generate_collect_command(config):
    340     koremLoc=os.path.join(CODE_DIR,'extra/accLoc.csv')
--> 341     cmd=CODE_DIR+"/bin/PTRMatrix.py {output_path} {ref_path} " + parser_c.min_orics + " {doric_path} " + koremLoc
    342     return cmd.format(**config)
    343 

NameError: global name 'parser_c' is not defined

In [487]:
len(c.samples[0][0])


Out[487]:
100000

In [709]:
c.run_pipeline()

In [512]:
[res,hfit,pfit]=c.compare_fit()


Simulated value: 0.5
Error from this fit: 0.134314081342 value: 0.432842959329
Error from initial PTR fit 0.189275182041 value: 0.405362408979

In [481]:
100*10**5/float(c.pop['NC_000913.3'].l)


Out[481]:
2.154405371191119

In [2255]:
tot_reads=[10**6]#np.linspace(5*10**3,10**4,15)
c=Community('comm0',[acc[s] for s in sel[0]],[growth_param[s] for s in sel[0]],td,'bowtie2')
c.build_index()
for nr in tot_reads:
    c.sample(nr)
    c.write_reads()

c.run_pipeline()
[res,hfit,pfit]=c.compare_fit()


/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py:293: UserWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg)
/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py:293: UserWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  warnings.warn(msg)
Ex
Ex
Ex
---------------------------------------------------------------------------
UnboundLocalError                         Traceback (most recent call last)
<ipython-input-2255-4c45e3923bb7> in <module>()
      7 
      8 c.run_pipeline()
----> 9 [res,hfit,pfit]=c.compare_fit()

<ipython-input-4-2a9bc4d0021b> in compare_fit(self)
     82                     print "Ex"
     83                     pass
---> 84         return [res,err_hfit,err_pfit]
     85 
     86 

UnboundLocalError: local variable 'res' referenced before assignment

In [448]:
np.log2(1.2978)


Out[448]:
0.37606807113368551

In [614]:
c.distribution


Out[614]:
array([ 0.66795389,  0.33204611])

In [616]:
c.pop.keys()


Out[616]:
['NC_007779.1', 'NC_000913.3']

In [570]:
tot_reads[0]*100/float(c.pop['NC_000913.3'].l)


Out[570]:
10.772026855955595

In [343]:
c.write_reads()

In [1]:
pl=plt.plot(tot_reads[0:6],hfit,tot_reads[0:6],pfit)
ax=pl[0]
#pl[0].set_xscale('log')
#plt.legend(['p(x) fit','Linear fit'])
#plt.plot(tot_reads,pfit,label='Linear fit')


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-f7a2c1730d8e> in <module>()
----> 1 pl=plt.plot(tot_reads[0:6],hfit,tot_reads[0:6],pfit)
      2 ax=pl[0]
      3 #pl[0].set_xscale('log')
      4 #plt.legend(['p(x) fit','Linear fit'])
      5 #plt.plot(tot_reads,pfit,label='Linear fit')

NameError: name 'plt' is not defined

In [ ]:
ax.set_yscale('log')

In [270]:
[res,hfit,pfit]=c.compare_fit()


Simulated value: 0.7
Error from this fit: 0.389631075905 value: 0.427258246866
Error from initial PTR fit 0.300151361166 value: 0.489894047184
Simulated value: 0.5
Error from this fit: 0.305636070824 value: 0.347181964588
Error from initial PTR fit -0.437847123705 value: 0.718923561852

In [234]:
res.plot()
pfit


Out[234]:
[0.1343119055981461, 0.13299914267733115]

In [177]:
start = time.time()
print("hello")
end = time.time()
print(end - start)


Out[177]:
u'/Users/hedani/Documents/GitRepos/PTRloc/PTR-Pipeline'

In [106]:
community=[]
nrs=np.linspace(5*10**3,10**6,10)
res=[]
for i in range(9):
    community.append(run_community_analysis('smr',td,'bowtie2',acc,growth_param,comm,nrs[i],0))
    res.append(community[i][0].compare_fit())


---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-106-dcb76e79b9a1> in <module>()
      3 res=[]
      4 for i in range(9):
----> 5     community.append(run_community_analysis('smr',td,'bowtie2',acc,growth_param,comm,nrs[i],0))
      6     res.append(community[i][0].compare_fit())

<ipython-input-83-3eaf4407c66f> in run_community_analysis(exec_string, td, mapper, acc, growth_param, community, nr_samples, topname)
     31         if 's' in exec_string:
     32             # create community and sample
---> 33             comm.append(Community(conf,c,acc,growth_param))
     34             comm[i].sample(nr_samples)
     35             comm[i].write_reads()

<ipython-input-103-3dfe25d04a80> in __init__(self, conf, comm, acc, growth_param)
      5         #self.samples=[]
      6 
----> 7         self.pop = self.init_population(comm,acc,growth_param)
      8         self.distribution = self.community_distribution()
      9         self.samples = []

<ipython-input-103-3dfe25d04a80> in init_population(self, comm, acc, growth_param)
     23         for i in comm:
     24             pop[acc[i]]=population(B = growth_param[i][0],C = growth_param[i][1],
---> 25                              D = growth_param[i][2], l = len(records[acc[i]]),
     26                              seq = records[acc[i]], cells=growth_param[i][3])
     27         return pop

KeyError: 'NC_000913.3'

In [ ]:
cp0=community[0].pop.values()
Gekv(cp0[0].B,cp0[0].C)

In [100]:
#PTRs
for i,c in enumerate(comm):
    print "comm"+str(i)
    for c_sub in c:
        print acc[c_sub]+": "+str(2**growth_param[c_sub][1])


comm0
NC_000913.3: 1.41421356237

In [30]:
cov1=2**(np.load(join(td,'comm1/Data/reads/npy',acc[comm[1][1]]+'.depth.npy')))
acc[comm[1][1]]


Out[30]:
'NC_002655.2'

In [58]:
#print c.pop[0].l/float(c.pop[1].l)*.8
#print c.pop[1].l/float(c.pop[0].l)*.2
c


Out[58]:
[0, 3]

In [32]:
indd=1;
x=np.linspace(0,1,len(cov1))
plt.plot(x*c.pop[indd].l,cov1/(cov1.sum()/len(cov1)*c.pop[indd].l))
plt.plot(x*c.pop[indd].l,piecewise_prob(x*c.pop[indd].l,c.pop[indd].C,c.pop[indd].l))
plt.plot(x*c.pop[indd].l,res.best_fit)


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-32-4d3560460be3> in <module>()
      1 indd=1;
      2 x=np.linspace(0,1,len(cov1))
----> 3 plt.plot(x*c.pop[indd].l,cov1/(cov1.sum()/len(cov1)*c.pop[indd].l))
      4 plt.plot(x*c.pop[indd].l,piecewise_prob(x*c.pop[indd].l,c.pop[indd].C,c.pop[indd].l))
      5 plt.plot(x*c.pop[indd].l,res.best_fit)

TypeError: 'builtin_function_or_method' object has no attribute '__getitem__'

In [24]:
signal=cov1/(cov1.sum()/len(cov1)*c.pop[indd].l)
res=fit_signal(signal,c.pop[indd].l)

In [25]:
err=(c.pop[indd].C-res.best_values['C'])/c.pop[indd].C
from_ptr=np.load(join(td,'comm1/Data/reads/npy',acc[comm[1][1]]+'.depth.best.npy'))
err2=(c.pop[indd].C-from_ptr[2]+from_ptr[3])/c.pop[indd].C
print err, err2


0.0777801701655 0.104943273296

In [26]:
c.pop[indd].C


Out[26]:
0.7

In [45]:
conf=local_conf(join(td,'comm1'),'bowtie2')

In [23]:
import example

In [48]:
example.main(join(conf['ref_path'],'Fasta'),False,conf['data_path'],conf['node_path'],conf['output_path']+"/out.df",None,True)

In [38]:
conf['ref_path']


Out[38]:
'/mnt/comm0/References'

In [111]:
refs=open_records(glob(join(conf['ref_path'],'Fasta','*.fasta')))

In [119]:
glob(join(conf['ref_path'],'Fasta','*.fasta'))


Out[119]:
[]

In [112]:
refs


Out[112]:
[]

In [57]:
import PTRC

In [94]:
import pickle

In [97]:
ind


Out[97]:
({'Bacteroides vulgatus': ['SegalLab|882']}, {})

In [101]:
f=open('/mnt/index.pk','r')
f.close()
f=open('/mnt/index.pk','w')
pickle.dump(ind,f)
f.close()

In [86]:
f=open('/mnt/index.pk','r')
ind=pickle.load(f)
pickle.save

In [93]:
ind=({'Bacteroides vulgatus': ['SegalLab|882']},{})

In [88]:
PTRC.coverage_analysis(args.i1, args.i2, args.m, args.pe, args.gz, \
                          args.outfol, args.qual_format, args.cov_thresh, dbpath = args.db_path_name)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-88-b08053f2169d> in <module>()
----> 1 PTRC.coverage_analysis(args.i1, args.i2, args.m, args.pe, args.gz,                           args.outfol, args.qual_format, args.cov_thresh, dbpath = args.db_path_name)

NameError: name 'args' is not defined

In [23]:
r=glob(conf['data_path']+'/*')

In [32]:
args=namedtuple('args','i1 i2 m pe gz outfol qual_format cov_thresh db_path_name')
args.i1=r[0]
args.i2=r[1]
#args.m=None#conf['data_path']+"/blah.m"
args.pe=True
args.gz=False
args.outfol=conf['output_path']
args.qual_format='offset-33'
args.cov_thresh=5
args.db_path_name=None

In [2]:
cd /mnt


/mnt

In [32]:
ro=pd.read_pickle('comm0/re.df')

In [37]:
ro['reads']


Out[37]:
Series([], Name: reads, dtype: float64)

In [27]:
import pandas
import pickle

In [32]:
cd /Users/hedani/Documents/Projects/MicrobialEcology/CommunitySampling


/Users/hedani/Documents/Projects/MicrobialEcology/CommunitySampling

In [33]:
f=open('comm0/sim_db.pk')
pk=pickle.load(f)

In [34]:
pk


Out[34]:
(defaultdict(list, {'Escherichia coli': ['NC_000913.3']}),
 {'NC_000913.3': 4641652},
 {'NC_000913.3': 'Escherichia coli'},
 {})

In [32]:
mutable_seq


Out[32]:
MutableSeq('AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAG...TTC', IUPACUnambiguousDNA())

In [67]:
mutated


Out[67]:
MutableSeq('AGCTTTTCATTCTGACTGCAACGGGCGATAAGTCTCTGTCTGGATTAAAAAACT...TTC', IUPACUnambiguousDNA())

In [72]:
rec2=rec
rec2.seq=mutated
rec2.id="mut_"+rec2.id
rec2.description="mut_"+rec2.description
rec2.name="mut_"+rec2.name

In [ ]:


In [66]:
mutated=mutateSeq(mutable_seq,.01)


/Library/Python/2.7/site-packages/ipykernel/__main__.py:7: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

In [71]:
with open("/Users/hedani/Documents/Projects/MicrobialEcology/CommunitySampling/comm0/References/Fasta/NC_000913.3_mutated.fasta", "w") as output_handle:
    SeqIO.write(rec2, output_handle, "fasta")

In [65]:
def mutateSeq(mutable_seq,rate):
    nr=['A','C','G','T']
    nt={'A': 0,'C': 1,'G': 2,'T': 3}
    
    prob=np.matrix('0 1 1 1; 1 0 1 1; 1 1 0 1; 1 1 1 0')/float(3)
    prob=prob.tolist()
    rand_ch=np.random.choice(range(len(mutable_seq)),np.round(len(mutable_seq)*rate))
    for i,c in enumerate(rand_ch):
        mutable_seq[c]=np.random.choice(nr,None,False,prob[nt[mutable_seq[c]]])
    
    return mutable_seq

In [46]:
IUPAC.unambiguous_dna.letters


Out[46]:
'GATC'

In [52]:
prob=np.matrix('0 1 1 1; 1 0 1 1; 1 1 0 1; 1 1 1 0')/float(3)
prob.tolist()


Out[52]:
[[0.0, 0.3333333333333333, 0.3333333333333333, 0.3333333333333333],
 [0.3333333333333333, 0.0, 0.3333333333333333, 0.3333333333333333],
 [0.3333333333333333, 0.3333333333333333, 0.0, 0.3333333333333333],
 [0.3333333333333333, 0.3333333333333333, 0.3333333333333333, 0.0]]

In [29]:
# import trimmed multiple alignment for creating MC
handle = open("/Users/Shared/ecol_maf2.fasta-gb", "rU")
seq=[]
for rec in SeqIO.parse(handle, "fasta"):
    seq.append(rec.seq)

In [54]:
for i in range(len(seq[0])):
    seq[i]=


Out[54]:
[Seq('CGGAAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGAGT...GAT', SingleLetterAlphabet()),
 Seq('CGGAAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGAGT...GAT', SingleLetterAlphabet()),
 Seq('CGGAAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGAGT...GAT', SingleLetterAlphabet()),
 Seq('CGGAAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGAGT...GAT', SingleLetterAlphabet()),
 Seq('CGGAAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGAGT...GAT', SingleLetterAlphabet()),
 Seq('CGGAAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGAGT...GAT', SingleLetterAlphabet()),
 Seq('CGGCAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGCGT...GAT', SingleLetterAlphabet()),
 Seq('CGGAAAAGAACTCGCCATCGCAGGTGGCGCAGTAGGCCACACCGCGACCAGAGT...GAT', SingleLetterAlphabet())]

In [9]:
ls comm0/References/Fasta/


NC_000913.3.fasta  NC_002655.2.fasta  NC_009614.1.fasta

In [31]:
from Bio import SeqIO
from Bio.Seq import MutableSeq
from Bio.Alphabet import IUPAC
handle = open("/Users/hedani/Documents/Projects/MicrobialEcology/CommunitySampling/comm0/References/Fasta/NC_000913.3.fasta", "rU")
for rec in SeqIO.parse(handle, "fasta"):
    seq = rec.seq
mutable_seq = MutableSeq(str(rec.seq), IUPAC.unambiguous_dna)
handle.close()

In [ ]:
for rec in SeqIO.parse("sequence.gb", "genbank"):
    if rec.features:
        for feature in rec.features:
            if feature.type == "CDS":
                print feature.location
                print feature.qualifiers["protein_id"]
                print feature.location.extract(rec).seq

In [ ]:
x=None

In [218]:
# Import relevant modules
import pymc
import numpy as np

# Some data

observations=[x for x in samp/float(1000) if x < 1/float(2)]
#observations=samp

#n = 5*np.ones(4,dtype=int)
#x = np.array([-.86,-.3,-.05,.73])

# Priors on unknown parameters
#alpha = pymc.Normal('alpha',mu=0,tau=.01)
#beta = pymc.Normal('beta',mu=0,tau=.01)

# Define prior
probs = pymc.Dirichlet(name="theta",theta=[1/float(2)]*3)
#probs_full = pymc.CompletedDirichlet("theta_full",theta)

# Arbitrary deterministic function of parameters
#@pymc.deterministic
#def theta(a=alpha, b=beta):
#	"""theta = logit^{-1}(a+b)"""
#	return pymc.invlogit(a+b*x)

def PB(x,Bb,Cc):
    return 2**(2*Cc*x)*(-1 + 2**Bb)

def PC(x,Cc):
    return (2**Cc - 2**(1 + 2*Cc*x) + 2**(Cc + 2*Cc*x))/2**Cc

def PD(x,Cc,Dd):
    return 2**(1 - Cc - Dd + 2*Cc*x)*(-1 + 2**Dd)

@pymc.stochastic(dtype=float,observed=True)
def P(probs=probs,value=observations):
    theta = np.hstack((probs,1-probs.sum()))
    #Bb=theta[0]
    #Cc=theta[1]
    #Dd=theta[2]
    return np.prod(PB(value,theta[0],theta[1])+PC(value,theta[1])+PD(value,theta[1],theta[2]))

alpha = pymc.Uniform('alpha',lower=0, upper=1000)
beta = pymc.Uniform('beta',lower=0, upper=1000)

C = pymc.Beta('C',alpha=alpha,beta=beta)

@pymc.stochastic(dtype=float,observed=True)
def pn(C=C,beta=beta,value=observations):
    #theta = np.hstack((probs,1-probs.sum()))
    #Bb=theta[0]
    #Cc=theta[1]
    #Dd=theta[2]
    return np.prod([pxn(x,C,1) for x in value])

#model = pymc.Model([probs,P])

model = pymc.Model(pn)

#graph = pymc.graph.graph(model)
#graph.write_png("graph.png")


#fitter=pymc.MAP(model)
#fitter.fit()

M = pymc.MCMC([C,model])
M.use_step_method(pymc.AdaptiveMetropolis,C)
M.sample(iter=1500, burn=400, thin=2)

print "Maximum posterior probs = ", C.value


 [-----------------80%----------        ] 1201 of 1500 complete in 103.2 sec
/Library/Python/2.7/site-packages/pymc/StepMethods.py:1272: UserWarning: 
Covariance was not positive definite and proposal_sd cannot be computed by 
Cholesky decomposition. The next jumps will be based on the last 
valid covariance matrix. This situation may have arisen because no 
jumps were accepted during the last `interval`. One solution is to 
increase the interval, or specify an initial covariance matrix with 
a smaller variance. For this simulation, each time a similar error 
occurs, proposal_sd will be reduced by a factor .9 to reduce the 
jumps and increase the likelihood of accepted jumps.
  warnings.warn(adjustmentwarning)
 [-----------------100%-----------------] 1500 of 1500 complete in 128.8 secMaximum posterior probs =  0.458280861811

In [219]:
pymc.Matplot.plot(M)


Plotting C

In [92]:
nr_samples=10000
        
samp=inverse_transform_sampling(.6,nr_samples,1000)


/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py:293: UserWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  warnings.warn(msg)
/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py:293: UserWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg)

In [208]:
observations=[x for x in samp/float(1000) if x < 1/float(2)]
[pxn(x,.6,1) for x in observations]


Out[208]:
[1.713822547073222,
 2.201394547941149,
 1.7925554377950603,
 1.946412833079749,
 1.728137172594991,
 2.0717037822746533,
 2.2234773946864945,
 2.0494229299514397,
 1.856284368758959,
 2.0596764670516716,
 1.749833601872147,
 2.3044424474675678,
 2.2834540163954471,
 2.4344868998647438,
 2.166877712490709,
 2.1813446578765125,
 1.951275840002392,
 1.7747523663044473,
 1.744021393459414,
 2.4426001863720925,
 1.870232637303288,
 1.6236241874153752,
 1.847043363703978,
 2.4003050011824101,
 1.8624707100876767,
 2.0613903704684531,
 1.915893991493232,
 2.3430987166439525,
 1.6785525890656403,
 2.1759082690510816,
 2.1904355060573164,
 2.4304403680546796,
 1.7109738821434899,
 1.7747523663044473,
 2.1759082690510816,
 1.6729771460812648,
 2.1977354553145911,
 1.7866013254813549,
 2.4043013594838953,
 2.1489288809409572,
 2.287255825086421,
 1.6925726238540164,
 2.4063020330755549,
 2.366603040263386,
 2.3353159320076466,
 2.2815554824269038,
 2.134676920964667,
 2.2346017495596424,
 1.8332680416739273,
 1.965937882024871,
 2.0803376973656271,
 2.2551411536630392,
 1.7052907493123903,
 1.7109738821434899,
 1.6494873769194809,
 1.9190838356396624,
 2.396315285508984,
 1.9873111042065146,
 1.6222742562804253,
 2.4426001863720925,
 2.2087310198211796,
 2.3903431427788959,
 2.0545433020155759,
 2.0392204372869469,
 2.3178990539648918,
 1.8165709790192808,
 2.0022439180733529,
 2.044315318976373,
 2.2948784423309658,
 1.6141982066970981,
 1.8090317609349087,
 1.946412833079749,
 2.1650761036501951,
 1.9856587925723361,
 2.4243831781419978,
 1.7895759053930809,
 1.6799493506753727,
 2.1258174846791014,
 2.0751530400166631,
 1.8937130477652351,
 1.8624707100876767,
 1.6249752418585832,
 2.2216287271277042,
 2.1686808204907249,
 1.885853671048763,
 1.6467456479985429,
 2.3824034321064569,
 2.2327438328697302,
 1.7166759548526067,
 2.0256960771683747,
 2.3744900938005755,
 1.797034040358231,
 2.0392204372869469,
 1.7310144133341487,
 1.8578290250062042,
 1.7267003461465853,
 2.3607049469584767,
 2.1759082690510816,
 1.6290351543889472,
 1.9127094494109886,
 1.7644490081048305,
 1.843973263090082,
 2.325623794287369,
 1.8090317609349087,
 2.2364612122640923,
 2.2720864638461391,
 1.8332680416739273,
 2.1258174846791014,
 2.1704854288977278,
 1.9989158477081181,
 2.2815554824269038,
 2.1849764608463911,
 1.9351128968157894,
 1.9528995408033487,
 2.2570177076706313,
 1.7238302761217605,
 1.7252647143199236,
 2.201394547941149,
 2.3430987166439525,
 1.8105370972174422,
 2.1614773784629215,
 1.9675737834429252,
 2.3784434718893932,
 1.964303340745619,
 1.8378483624703206,
 1.9545245927225008,
 2.4183410840737238,
 1.744021393459414,
 2.2853541301770872,
 2.0890075946710644,
 2.2160919415753586,
 2.3353159320076466,
 2.0256960771683747,
 2.1722915389602409,
 2.4043013594838953,
 2.1596802596263758,
 1.8317438061843694,
 1.8717889002205601,
 1.9610383340810176,
 2.3923322014224384,
 1.8531989080209592,
 1.6426415970376906,
 2.2664239265584767,
 2.4163304008444264,
 2.0477189766331381,
 2.1364532339570679,
 1.7123976222451709,
 1.6440084762979756,
 1.6494873769194809,
 2.0324470079915362,
 1.7425713603807635,
 2.3372592012633575,
 2.1382310250597234,
 1.6399112469687223,
 2.044315318976373,
 1.713822547073222,
 1.6771569887665114,
 2.3082792005888151,
 2.4123140482277061,
 2.2087310198211796,
 1.6939810518520932,
 2.1795310207722207,
 2.0596764670516716,
 2.3863699854105977,
 2.1813446578765125,
 1.8795896613964591,
 1.9095302005800006,
 1.6785525890656403,
 1.9415619458358537,
 2.0528350914623026,
 1.9206807487450022,
 1.7024562656834648,
 2.2068946128524392,
 2.013935860541598,
 2.4143213893590665,
 2.2739771184299546,
 2.3314342392345435,
 1.7324548296155262,
 1.7310144133341487,
 2.3294958130315937,
 2.1777188915816517,
 1.7024562656834648,
 2.4123140482277061,
 1.6236241874153752,
 2.3159718813316736,
 1.9254794654502114,
 2.134676920964667,
 2.2032263781813621,
 1.6701963731637748,
 2.2346017495596424,
 2.0172889382925381,
 2.287255825086421,
 2.4264005622695906,
 2.0786080405579432,
 1.7252647143199236,
 1.6701963731637748,
 2.3903431427788959,
 1.9079425583836909,
 1.9757737322612936,
 1.6290351543889472,
 2.1560905032507582,
 2.3159718813316736,
 2.0426156122830812,
 2.0122614122653224,
 1.6688077204798046,
 2.0665406425338078,
 2.0648224570222724,
 1.8968659620221067,
 1.7762291779266652,
 2.3705432868993839,
 1.8060248420780252,
 1.6276807244094718,
 1.9873111042065146,
 2.4344868998647438,
 2.4123140482277061,
 2.3101999719682347,
 1.7718024256390705,
 1.6577399829425863,
 1.7940470637561052,
 2.2986992719794448,
 2.1813446578765125,
 2.0803376973656271,
 2.3392040875574089,
 2.2758693462702522,
 2.0528350914623026,
 2.0528350914623026,
 1.7732767825488778,
 1.6426415970376906,
 2.044315318976373,
 2.4123140482277061,
 1.8241416171210074,
 2.3450484621309515,
 1.6385477742710401,
 2.3567810533188873,
 1.7910650520177611,
 1.6399112469687223,
 2.2626567440304379,
 2.2588958232012959,
 2.3489528217370874,
 2.2739771184299546,
 2.3198278302420858,
 2.0768798218401789,
 1.8150606259816562,
 1.9431775634855639,
 2.0528350914623026,
 2.1275864255340031,
 1.9561509968841435,
 1.6263274205447789,
 1.7367832740680103,
 1.9955933091625873,
 2.4284196251081864,
 2.3430987166439525,
 2.3626693429883048,
 1.7267003461465853,
 2.0290687349474035,
 2.1169948172415651,
 2.30061207168339,
 2.2513927250249601,
 2.2032263781813621,
 2.325623794287369,
 2.1258174846791014,
 2.0579639886269554,
 2.3372592012633575,
 2.396315285508984,
 2.2476505269136178,
 2.0768798218401789,
 2.2253276005634137,
 1.6426415970376906,
 2.1152346825758759,
 2.134676920964667,
 2.2216287271277042,
 2.1886143105419724,
 1.6660338778576991,
 2.2383222222695554,
 1.6331052103906272,
 1.8135515286980448,
 2.2588958232012959,
 1.8749053121375208,
 2.1222840140021297,
 2.0055775294684253,
 2.1417910465166972,
 1.7123976222451709,
 1.9989158477081181,
 1.7252647143199236,
 2.321758211497686,
 1.7166759548526067,
 2.2087310198211796,
 1.9724896598526351,
 2.287255825086421,
 2.1258174846791014,
 1.7703292945549882,
 2.1578846349697569,
 2.013935860541598,
 2.1507170533119009,
 2.0005791908379988,
 2.2364612122640923,
 1.839377677699364,
 1.7109738821434899,
 2.0751530400166631,
 2.1525067136615572,
 2.2664239265584767,
 1.8045232574233308,
 1.9528995408033487,
 2.1258174846791014,
 2.4183410840737238,
 2.0665406425338078,
 2.3140463110091143,
 2.4304403680546796,
 1.7880879968907495,
 2.1134760113403361,
 2.2197815966080348,
 2.0959695204602458,
 1.8180825888558607,
 2.3101999719682347,
 1.8686776683109168,
 2.1099630542937522,
 1.6426415970376906,
 2.192258217030679,
 2.1507170533119009,
 2.30061207168339,
 1.9318964096400044,
 2.2383222222695554,
 2.4143213893590665,
 1.959407866436766,
 1.9545245927225008,
 1.8516581013955205,
 2.4426001863720925,
 1.7985293930658586,
 2.3353159320076466,
 1.7310144133341487,
 1.7223970305596761,
 1.6577399829425863,
 1.9545245927225008,
 1.8686776683109168,
 2.2986992719794448,
 2.1596802596263758,
 1.8780269126861882,
 2.1382310250597234,
 2.0562529340095228,
 1.7762291779266652,
 2.30061207168339,
 2.4183410840737238,
 2.210568954905733,
 1.7310144133341487,
 2.2142494145096552,
 1.7152486576134809,
 1.9016051787740045,
 2.3705432868993839,
 2.4243831781419978,
 2.3372592012633575,
 2.3275589984958649,
 1.8045232574233308,
 1.6799493506753727,
 1.8516581013955205,
 2.183159804149299,
 1.7940470637561052,
 1.9111191638893235,
 1.9807100957984993,
 1.8889934999329971,
 2.0768798218401789,
 1.6996264934490304,
 1.6660338778576991,
 2.2739771184299546,
 1.7688573882774419,
 1.9561509968841435,
 2.2701973812107501,
 2.4365126915279181,
 1.7985293930658586,
 1.9367231480658114,
 2.2834540163954471,
 2.2050597327289885,
 2.1012061866568619,
 1.9561509968841435,
 1.9610383340810176,
 1.939947671460396,
 1.8165709790192808,
 1.7910650520177611,
 1.9675737834429252,
 1.9972538875329251,
 2.4446327292846441,
 2.1169948172415651,
 2.044315318976373,
 1.9922762932419316,
 2.2234773946864945,
 2.3159718813316736,
 2.1222840140021297,
 1.6168857461557915,
 2.4183410840737238,
 2.2796585269579497,
 2.0648224570222724,
 2.4446327292846441,
 2.192258217030679,
 1.6897592799154439,
 2.087270729443846,
 2.0751530400166631,
 2.3607049469584767,
 1.7527469658144854,
 2.0734276938927159,
 1.870232637303288,
 1.7425713603807635,
 1.9143010582452387,
 1.965937882024871,
 2.0005791908379988,
 2.2588958232012959,
 1.7469250804368135,
 1.9286852688136598,
 2.3392040875574089,
 1.6331052103906272,
 1.6855480295594616,
 1.6141982066970981,
 2.0977136241483993,
 2.2815554824269038,
 2.1117188023182041,
 2.2495208478011124,
 2.2308874609089528,
 2.1525067136615572,
 1.6605000167463213,
 1.6646486860004741,
 2.3943229152067076,
 2.0977136241483993,
 2.3843858839370027,
 2.3567810533188873,
 2.106455936380021,
 2.325623794287369,
 2.3314342392345435,
 1.783631689828,
 1.8516581013955205,
 2.1777188915816517,
 1.749833601872147,
 1.6605000167463213,
 2.0786080405579432,
 2.0477189766331381,
 2.1417910465166972,
 1.8485803303017119,
 2.3509074385574649,
 2.2570177076706313,
 1.7880879968907495,
 2.3314342392345435,
 1.8075276762353043,
 2.2346017495596424,
 2.2701973812107501,
 2.1704854288977278,
 2.132902084853574,
 2.0717037822746533,
 2.3333742784458216,
 1.9367231480658114,
 1.8795896613964591,
 1.9190838356396624,
 1.9873111042065146,
 1.8609221946005545,
 1.9254794654502114,
 1.7615161933753989,
 2.3101999719682347,
 1.8075276762353043,
 2.396315285508984,
 2.3044424474675678,
 2.1813446578765125,
 1.9174882502570421,
 2.1099630542937522,
 1.713822547073222,
 2.0055775294684253,
 2.0613903704684531,
 1.9708496712235921,
 1.959407866436766,
 2.0189675700870375,
 1.9626701584743083,
 2.4304403680546796,
 2.1400102955026021,
 1.9692110461315844,
 2.2796585269579497,
 2.0648224570222724,
 2.2457817610684891,
 1.7866013254813549,
 1.7940470637561052,
 2.2142494145096552,
 2.3275589984958649,
 1.6729771460812648,
 1.6605000167463213,
 1.9206807487450022,
 2.3923322014224384,
 2.1650761036501951,
 2.2364612122640923,
 2.3178990539648918,
 2.0959695204602458,
 1.9545245927225008,
 1.6660338778576991,
 2.3159718813316736,
 2.3198278302420858,
 2.3430987166439525,
 2.0768798218401789,
 2.0924856621827836,
 2.1099630542937522,
 1.9873111042065146,
 1.6953906518349915,
 2.1867946292246474,
 2.0734276938927159,
 2.1258174846791014,
 2.044315318976373,
 2.3353159320076466,
 2.0768798218401789,
 1.7585882534834925,
 1.8332680416739273,
 2.2439145489728172,
 2.1382310250597234,
 1.9626701584743083,
 2.4223674713296774,
 2.0290687349474035,
 1.6897592799154439,
 2.0392204372869469,
 2.1614773784629215,
 2.3101999719682347,
 1.9724896598526351,
 2.2570177076706313,
 2.2739771184299546,
 2.1904355060573164,
 1.6632646458344045,
 1.8180825888558607,
 1.8135515286980448,
 2.1012061866568619,
 2.063105700063069,
 1.9577787544135095,
 2.2532661598802224,
 1.7324548296155262,
 1.9561509968841435,
 2.3121223416650016,
 1.7940470637561052,
 2.0768798218401789,
 2.4284196251081864,
 2.3314342392345435,
 1.6813472745620617,
 2.0924856621827836,
 2.4203534404380567,
 1.691165366866338,
 2.0240118510859415,
 2.2758693462702522,
 1.8686776683109168,
 2.2253276005634137,
 1.9889647907666099,
 1.9351128968157894,
 1.890565374188786,
 2.2050597327289885,
 1.8060248420780252,
 1.90318754749327,
 1.9577787544135095,
 2.3943229152067076,
 1.6799493506753727,
 1.7673867057880892,
 1.9335039843782302,
 1.7166759548526067,
 2.2197815966080348,
 1.7181044397780765,
 1.6168857461557915,
 2.3824034321064569,
 2.2364612122640923,
 1.6688077204798046,
 2.1886143105419724,
 1.9972538875329251,
 1.6771569887665114,
 2.0613903704684531,
 2.2142494145096552,
 1.9335039843782302,
 1.948032487262553,
 2.1614773784629215,
 2.4324627925071254,
 2.3121223416650016,
 2.3705432868993839,
 1.7223970305596761,
 2.3236901990672303,
 1.8135515286980448,
 2.2087310198211796,
 1.9906198533967292,
 1.959407866436766,
 2.2364612122640923,
 1.951275840002392,
 1.6869506122335969,
 2.3025264630712718,
 1.9302901714889078,
 2.3784434718893932,
 1.9577787544135095,
 2.0055775294684253,
 1.6968014247779475,
 2.4243831781419978,
 2.3804226285473944,
 1.6771569887665114,
 2.0648224570222724,
 1.7152486576134809,
 1.9016051787740045,
 1.9127094494109886,
 1.8120436861243752,
 1.7425713603807635,
 2.0189675700870375,
 2.0890075946710644,
 2.0838013300384186,
 2.0172889382925381,
 2.3548215529943701,
 2.2197815966080348,
 2.4324627925071254,
 2.0699813039697834,
 2.4284196251081864,
 1.8764654632933062,
 2.2197815966080348,
 2.3824034321064569,
 1.7732767825488778,
 2.0358309056498398,
 1.6222742562804253,
 1.8686776683109168,
 2.2179360018495458,
 2.0734276938927159,
 2.2834540163954471,
 2.0392204372869469,
 1.8485803303017119,
 1.959407866436766,
 1.9577787544135095,
 2.0528350914623026,
 2.4304403680546796,
 1.6925726238540164,
 2.2216287271277042,
 2.087270729443846,
 1.6785525890656403,
 1.6536085316901676,
 2.1704854288977278,
 1.7806669902148173,
 1.948032487262553,
 2.1904355060573164,
 1.8135515286980448,
 2.0734276938927159,
 1.6317473964412894,
 2.4103083760615611,
 2.3764659607631691,
 2.0511283011678727,
 2.1777188915816517,
 1.8485803303017119,
 1.7324548296155262,
 2.30061207168339,
 1.8984443871318231,
 2.192258217030679,
 2.1082087660522606,
 2.0426156122830812,
 2.3725158696346043,
 1.9989158477081181,
 2.1995642407409925,
 1.6453764929694381,
 1.8952888492636204,
 1.8105370972174422,
 2.3411505922353766,
 1.8686776683109168,
 2.0734276938927159,
 2.3784434718893932,
 2.2457817610684891,
 1.6968014247779475,
 1.9063562362019808,
 2.4203534404380567,
 1.7483787363444852,
 2.192258217030679,
 1.9741310131533445,
 1.7396749100275173,
 2.3314342392345435,
 2.063105700063069,
 2.2439145489728172,
 1.8937130477652351,
 1.9286852688136598,
 1.9351128968157894,
 2.4143213893590665,
 2.0977136241483993,
 1.8424401269498638,
 1.9367231480658114,
 2.1777188915816517,
 1.959407866436766,
 1.9840078547209199,
 2.1650761036501951,
 2.0324470079915362,
 2.3528636818587483,
 1.9016051787740045,
 2.3236901990672303,
 2.2701973812107501,
 2.1977354553145911,
 2.4103083760615611,
 2.0172889382925381,
 1.6577399829425863,
 2.1258174846791014,
 2.3392040875574089,
 2.287255825086421,
 1.7880879968907495,
 1.9367231480658114,
 1.8811537105053111,
 1.8378483624703206,
 2.0105883561768949,
 1.8165709790192808,
 2.106455936380021,
 1.9528995408033487,
 2.4103083760615611,
 2.4183410840737238,
 1.8764654632933062,
 2.2720864638461391,
 2.2439145489728172,
 1.7673867057880892,
 2.1777188915816517,
 1.843973263090082,
 2.2420488893347654,
 2.4063020330755549,
 2.1777188915816517,
 1.9955933091625873,
 2.4264005622695906,
 1.7324548296155262,
 2.3430987166439525,
 2.2607755015544111,
 1.797034040358231,
 2.4243831781419978,
 1.8165709790192808,
 1.9143010582452387,
 2.192258217030679,
 2.402302349313854,
 1.7747523663044473,
 2.2929704097429822,
 2.2142494145096552,
 1.8952888492636204,
 2.003910030565927,
 1.6453764929694381,
 2.4123140482277061,
 1.797034040358231,
 2.1940824447231067,
 1.7629819908784712,
 2.2327438328697302,
 2.3923322014224384,
 2.2626567440304379,
 1.8060248420780252,
 2.106455936380021,
 1.7644490081048305,
 2.2383222222695554,
 2.227179346038529,
 2.2777631486761734,
 1.8889934999329971,
 1.6467456479985429,
 2.3528636818587483,
 1.6841466130357414,
 2.0959695204602458,
 2.0223290253209236,
 1.9955933091625873,
 1.8655716077985363,
 2.2068946128524392,
 2.0409173187782623,
 1.964303340745619,
 2.3646353736365393,
 1.7469250804368135,
 2.1117188023182041,
 1.6344641542076912,
 2.1169948172415651,
 1.7556651803263019,
 2.3314342392345435,
 2.1117188023182041,
 2.0358309056498398,
 2.3025264630712718,
 1.6605000167463213,
 1.9955933091625873,
 2.1047045640643343,
 2.1258174846791014,
 2.1489288809409572,
 1.7542054662447801,
 2.2383222222695554,
 1.6331052103906272,
 2.2910639635520336,
 2.0579639886269554,
 2.4426001863720925,
 2.2327438328697302,
 1.659119425911592,
 2.3705432868993839,
 1.8655716077985363,
 2.4426001863720925,
 1.9496534891962694,
 2.4344868998647438,
 1.8842857142480207,
 2.3353159320076466,
 1.7238302761217605,
 1.8686776683109168,
 2.0942268668735942,
 2.087270729443846,
 2.2815554824269038,
 2.3784434718893932,
 2.4183410840737238,
 2.3626693429883048,
 2.0768798218401789,
 2.4284196251081864,
 2.3392040875574089,
 2.4304403680546796,
 1.8241416171210074,
 2.3587421841879848,
 1.6536085316901676,
 2.4043013594838953,
 2.3607049469584767,
 2.0890075946710644,
 1.7267003461465853,
 2.120519481735418,
 2.2476505269136178,
 2.3824034321064569,
 2.4405693333782938,
 2.2720864638461391,
 1.7483787363444852,
 2.1099630542937522,
 2.0665406425338078,
 1.946412833079749,
 2.3275589984958649,
 1.7600516145814959,
 2.1813446578765125,
 2.1222840140021297,
 1.9254794654502114,
 1.6249752418585832,
 2.0977136241483993,
 2.2532661598802224,
 2.0838013300384186,
 2.2142494145096552,
 1.7038729180836609,
 1.7762291779266652,
 1.6358242288326703,
 1.9206807487450022,
 1.7777072184372671,
 2.227179346038529,
 1.8717889002205601,
 1.7267003461465853,
 1.6481159423325433,
 2.1134760113403361,
 2.3903431427788959,
 1.6358242288326703,
 2.2986992719794448,
 2.3101999719682347,
 2.4243831781419978,
 2.3923322014224384,
 2.1886143105419724,
 1.7615161933753989,
 2.4344868998647438,
 1.8060248420780252,
 2.2929704097429822,
 2.1240500145760866,
 1.8780269126861882,
 2.368572344230178,
 1.6168857461557915,
 1.9095302005800006,
 2.3705432868993839,
 1.9143010582452387,
 2.0172889382925381,
 2.2758693462702522,
 1.959407866436766,
 2.0717037822746533,
 1.8271786992797647,
 2.2588958232012959,
 2.2346017495596424,
 2.2439145489728172,
 1.7659172460694348,
 2.1596802596263758,
 2.0838013300384186,
 1.6399112469687223,
 1.8889934999329971,
 2.2645395519309179,
 2.174099151927825,
 2.0392204372869469,
 1.7718024256390705,
 1.6331052103906272,
 2.3469998300453105,
 2.2626567440304379,
 1.8937130477652351,
 2.2401847808635744,
 2.0358309056498398,
 1.783631689828,
 1.964303340745619,
 1.8952888492636204,
 2.0156117021641919,
 2.3198278302420858,
 2.1813446578765125,
 1.9222789906778879,
 1.8075276762353043,
 2.0256960771683747,
 1.8795896613964591,
 1.915893991493232,
 1.843973263090082,
 1.9174882502570421,
 2.3063600261978547,
 1.7052907493123903,
 1.6785525890656403,
 1.7688573882774419,
 1.843973263090082,
 1.6344641542076912,
 1.8749053121375208,
 2.3275589984958649,
 1.7038729180836609,
 2.1867946292246474,
 1.733896444499901,
 2.3923322014224384,
 2.2834540163954471,
 1.7310144133341487,
 1.6953906518349915,
 2.2160919415753586,
 1.8968659620221067,
 2.4043013594838953,
 1.8075276762353043,
 1.9939341114482292,
 1.9174882502570421,
 1.6522336706482783,
 1.9111191638893235,
 2.0803376973656271,
 1.6331052103906272,
 1.783631689828,
 2.0460164400340841,
 1.7910650520177611,
 1.7851158901363389,
 2.2420488893347654,
 1.7659172460694348,
 2.1759082690510816,
 1.959407866436766,
 2.4063020330755549,
 2.1222840140021297,
 2.2457817610684891,
 2.0256960771683747,
 2.1759082690510816,
 2.1152346825758759,
 2.0358309056498398,
 2.0699813039697834,
 1.7851158901363389,
 2.4183410840737238,
 1.7367832740680103,
 2.2050597327289885,
 1.9302901714889078,
 1.7585882534834925,
 1.733896444499901,
 1.6182311933894769,
 2.1117188023182041,
 2.0613903704684531,
 1.7152486576134809,
 2.2050597327289885,
 2.1471421953115737,
 2.2796585269579497,
 2.3333742784458216,
 2.0648224570222724,
 2.30061207168339,
 2.2327438328697302,
 2.0820687934598978,
 2.2124084193776792,
 2.2683098692168224,
 2.0324470079915362,
 2.3744900938005755,
 1.7955399309328812,
 1.6563616868849336,
 1.6827463616928642,
 2.0665406425338078,
 2.3121223416650016,
 2.120519481735418,
 2.3294958130315937,
 1.7542054662447801,
 1.8749053121375208,
 2.2457817610684891,
 2.4223674713296774,
 1.7038729180836609,
 1.839377677699364,
 1.9610383340810176,
 1.964303340745619,
 1.8241416171210074,
 1.9610383340810176,
 1.9270817005039056,
 2.1686808204907249,
 2.3063600261978547,
 2.0613903704684531,
 2.0545433020155759,
 1.6618817564019388,
 2.1904355060573164,
 1.9626701584743083,
 1.9335039843782302,
 1.9238785625440689,
 1.6757625488124366,
 1.7052907493123903,
 2.1152346825758759,
 1.7010407911316876,
 2.3784434718893932,
 2.3705432868993839,
 1.8733464581394381,
 1.9047712329373705,
 1.7010407911316876,
 2.2701973812107501,
 2.2645395519309179,
 2.2607755015544111,
 1.7615161933753989,
 1.8686776683109168,
 1.9190838356396624,
 1.7542054662447801,
 2.1904355060573164,
 2.120519481735418,
 2.1777188915816517,
 1.6618817564019388,
 1.8609221946005545,
 2.1152346825758759,
 1.7542054662447801,
 2.2476505269136178,
 1.8578290250062042,
 1.8531989080209592,
 2.4223674713296774,
 2.003910030565927,
 1.667420222364427,
 2.106455936380021,
 2.0290687349474035,
 2.2179360018495458,
 2.3824034321064569,
 1.959407866436766,
 2.3725158696346043,
 ...]

In [205]:
pxn2(.5,.6)


Out[205]:
1.6128561126127094

In [228]:
import pymc
import numpy as np


n = 5*np.ones(4,dtype=int)
x = np.array([-0.86,-0.3,-0.05,.73])


alpha = pymc.Normal('alpha',mu=0,tau=.01)
beta = pymc.Normal('beta',mu=0,tau=.01)


@pymc.deterministic
def theta(a=alpha, b=beta):
	    """theta = logit^{−1}(a+b)"""
            return pymc.invlogit(a+b*x)


d = pymc.Binomial('d', n=n, p=theta, value=np.array([0.,1.,3.,5.]),
                         observed=True)

In [230]:
S = pymc.MCMC([theta],db='pickle')
S.sample(iter = 10000, burn = 5000, thin = 2)
pymc.Matplot.plot(S)


 [-----------------100%-----------------] 10000 of 10000 complete in 0.1 secPlotting theta_0
Plotting theta_1
Plotting theta_2
Plotting theta_3

In [226]:
pymc.Matplot.plot(S)

In [ ]: