def findmaxima(c, var = 'x'):
    """
    Creates a RK-4 approximation of Rossler curve with the given c value
    and returns local maxima along the x(t) curve from T0 to T
    """
    T0 = 250
    e = 3E-4

    Rc = Rossler(c, dt=0.01, T0 = T0)
    Rc.run()

    if var == 'x':
        var = Rc.x
    elif var == 'y':
        var = Rc.y
    else:
        var = Rc.z

    #using only valuse where t > T0
    initial_index = np.where(Rc.t == T0)[0][0]

    #moving back 1 to get a more exact diff
    usable_x = var[initial_index - 1:]
    x_diff = np.diff(usable_x)
    usable_x = usable_x[1:]

    np.set_printoptions(threshold='nan')


    critical_points = usable_x[np.abs(x_diff) < e]
    return critical_points[critical_points > 0]

def main():
    
    # read data from Bingo-F protocols
    in_dt = InputData(BingoParser(options["bingo_data_dir"]))

    phs = in_dt.GetPhotos()
    pts = in_dt.GetPoints()
    cam_m, distor = in_dt.GetCameras().values()[0].GetParams()
    

    # compute relative orientations and merged them into common system
    ros = RelativeOrientation(in_dt)

    np.set_printoptions(precision=3, suppress=True)

    #PlotScene(pts, phs)
    #PlotRelOrs(pts, phs)

    # performs helmert transformation into world coordinate system 
    HelmertTransform(pts, phs)

    # run bundle block adjustment
    ph_errs, pt_errs = RunBundleBlockAdjutment(in_dt, options['protocol_dir'], flags['f'])   

    PlotScene(pts, phs)
    return 

    def info(self):
        """

        """
        print '{0:4} , {1:15}, {2:5}, {3:5}, {4:7}, {5:6}, {6:8}, {7:9}'.format('type', 'p', 'value', 'std', 'runable' , 'usable' , 'obsolete' , 'evaluated')
        np.set_printoptions(precision=3)
        print '{0:4} , {1:15}, {2:5}, {3:5}, {4:7}, {5:6}, {6:8}, {7:9}'.format(self.type, self.p, self.value, self.std, self.runable, self.usable , self.obsolete , self.evaluated)

def _test():
    import doctest

    start_suppress = np.get_printoptions()["suppress"]
    np.set_printoptions(suppress=True)
    doctest.testmod()
    np.set_printoptions(suppress=start_suppress)

def generate_eventdir_matrix(fileName, header=True, direction=None):
    '''write out directed event matrix from fast5, False if not present'''
    try: # check to make sure the file actually exists
        h5File = h5py.File(fileName, 'r')
        h5File.close()
    except:
        return False
    with h5py.File(fileName, 'r') as h5File:
      rowData = get_telemetry(h5File, "000", fileName)
      channel = str(rowData['channel'])
      mux = str(rowData['mux'])
      runID = rowData['runID']
      dir = "complement" if (direction=="r") else "template"
      eventLocation = "/Analyses/Basecall_1D_000/BaseCalled_%s/Events/" % (dir)
      if(not eventLocation in h5File):
          return False
      readName = rowData['read']
      sampleRate = str(int(rowData['sampleRate']))
      rawStart = str(rowData['rawStart'])
      outData = h5File[eventLocation]
      dummy = h5File[eventLocation]['move'][0] ## hack to get order correct
      numpy.set_printoptions(precision=15)
      headers = h5File[eventLocation].dtype
      if(header):
          sys.stdout.write("runID,channel,mux,read,sampleRate,rawStart,"+",".join(headers.names)+"\n")
      for line in outData:
        res=[repr(x) for x in line]
        # data seems to be normalised, but just in case it isn't in the future,
        # here's the formula for calculation:
        # pA = (raw + offset)*range/digitisation
        # (using channelMeta[("offset", "range", "digitisation")])
        # - might also be useful to know start_time from outMeta["start_time"]
        #   which should be subtracted from event/start
        sys.stdout.write(",".join((runID,channel,mux,readName,sampleRate,rawStart)) + "," + ",".join(res) + "\n")

def align_se3(model, data, precision=False):
    """Align two trajectories using the method of Horn (closed-form).

    Input:
    model -- first trajectory (3xn)
    data -- second trajectory (3xn)

    Output:
    R -- rotation matrix (3x3)
    t -- translation vector (3x1)
    t_error -- translational error per point (1xn)

    """
    if not precision:
        np.set_printoptions(precision=3, suppress=True)
    model_zerocentered = model - model.mean(1).reshape(model.shape[0], 1)
    data_zerocentered = data - data.mean(1).reshape(data.shape[0], 1)

    W = np.zeros((3, 3))
    for column in range(model.shape[1]):
        W += np.outer(model_zerocentered[:, column], data_zerocentered[:, column])
    U, d, Vh = np.linalg.linalg.svd(W.transpose())
    S = np.matrix(np.identity(3))
    if np.linalg.det(U) * np.linalg.det(Vh) < 0:
        S[2, 2] = -1
    R = U * S * Vh
    t = data.mean(1).reshape(data.shape[0], 1) - R * model.mean(1).reshape(model.shape[0], 1)

    model_aligned = R * model + t
    alignment_error = model_aligned - data
    t_error = np.sqrt(np.sum(np.multiply(alignment_error, alignment_error), 0)).A[0]

    return R, t, t_error

	def general_mix(self, mix, **kwargs):
		"""
			simple mix data test
		"""
		npr.seed(122351)
		self.mix = mix(K = self.nClass,**kwargs)
		self.mix.set_data(self.Y)
		self.mu_sample = list()
		self.sigma_sample = list()
		self.p_sample = list()
		for k in range(self.nClass):
			self.mu_sample.append(np.zeros_like(self.Thetas[k])) 
			self.sigma_sample.append(np.zeros_like(self.Sigmas[k])) 
			self.p_sample.append(np.zeros_like(self.P[k])) 
			
			
		for i in range(self.sim):  # @UnusedVariable
			self.mix.sample()
			for k in range(self.nClass):
				self.mu_sample[k] += self.mix.mu[k]
				self.sigma_sample[k] += self.mix.sigma[k]
				self.p_sample[k] += self.mix.p[k]
		np.set_printoptions(precision=2)
			
		self.compare_class("MCMC:")

def branch_PSSM(peak_branch_df, fa_dict):
    base_dict = {"A":0, "C":1, "T":2, "G":3}
    nuc_prob = gc_content(fa_dict)
    
    pos_matrix_branch = np.zeros([4,5])
    counter = 0
    if type(peak_branch_df) == str:
        with open(peak_branch_df) as f:
            for line in f:
                counter += 1
                seq = line.strip()
                for a, base in enumerate(seq):
                    pos_matrix_branch[base_dict[base],a] += 1
    else:
        for seq in peak_branch_df['Branch seq']:
            counter += 1
            seq = seq[2:7]
            for a, base in enumerate(seq):
                pos_matrix_branch[base_dict[base],a] += 1

    float_formatter = lambda x: "%.1f" % x
    np.set_printoptions(formatter={'float_kind':float_formatter})
    
    a = 0
    while a < 4:
        b = 0
        while b < 5:
            if pos_matrix_branch[a,b] == 0: pos_matrix_branch[a,b] += 1
            pos_matrix_branch[a,b] = np.log2((pos_matrix_branch[a,b]/float(counter))/nuc_prob[a])
            b += 1
        a += 1
    
    return pos_matrix_branch

def main():
    numpy.set_printoptions(precision=1, linewidth=284, threshold=40, edgeitems=13)

    X = []
    Y = []
    order = 2

    coeffs = raw_readings_norm_coeffs = {"temp": 100,
                                        "light": 100,
                                        "humidity" : 100, "pressure": 1e5, 
                                        "audio_p2p": 100, "motion" : 1}

    data_provider = dataprovider.DataProvider(order=order, debug=True,
                                                      start_time = 1379984887,
                                                      stop_time = 1379984887+3600*24*2,
                                                      device_list = ["17030002"],
                                                      eliminate_const_one=True, device_groupping="dict",
                                                  raw_readings_norm_coeffs = coeffs)
    f = open("2_order_knode.p", "rb")
    knode = pickle.load(f)
    f.close()
    

    for data in data_provider:
        print data
        t, d = data["17030002"]
        X.append(t)
        l = knode.label(d)[0]
        Y.append(l)


    plt.plot(X, Y, 'ro')
    plt.show()

def generate_PSSM(seq_list, fasta_dict):
    #Populate gene dictionary and build genome
    genome = fasta_dict
    nuc_prob = gc_content(fasta_dict)

    base_dict = {"A":0, "C":1, "T":2, "G":3}
    
    #First generate a consensus matrix for the sequence, where 1st row is A counts, second row is C, third row is T, fourth row is G.
    PSSM = np.zeros([4,len(seq_list[0])])

    counter = 0
    for seq in seq_list:
        counter += 1
        for a, base in enumerate(seq):
            PSSM[base_dict[base],a] += 1

    float_formatter = lambda x: "%.1f" % x
    np.set_printoptions(formatter={'float_kind':float_formatter})
    
    a = 0
    while a < 4:
        b = 0
        while b < len(seq_list[0]):
            if PSSM[a,b] == 0: PSSM[a,b] += 1
            PSSM[a,b] = np.log2((PSSM[a,b]/float(counter))/nuc_prob[a])
            b += 1
        a += 1
    
    return PSSM

def main(file,N,db_file):
    
    original_db = list(SeqIO.parse(db_file, 'fasta'))
    original_db_dict = defaultdict(Bio.SeqRecord.SeqRecord)
    for i in original_db:
        original_db_dict[i.id] = i
        
    
    np.set_printoptions(suppress=True)
    with open(file) as f:
        lines = f.readlines()
    lines = [l.strip().split() for l in lines if l[0] != "#"]
    mat = np.array([ [v[2],v[3],v[4],v[5],v[6],v[7],v[8],v[9],v[10],v[11]] for v in lines]  ,dtype=float)
    id = np.array([ [l[0],l[1]] for l in lines],dtype=str)
    
    gaps = xrange(100,N-1,-1)
    for i in gaps:
        hit = (mat[:,0] >= i) & (mat[:,0] < i+1) &(mat[:,1] > 370)
        id_result = id[hit,:]
        mat_result = mat[hit,:]
        np.savetxt('.'.join([file,str(i)]),np.hstack((id_result,np.asarray(mat_result,dtype='str'))),delimiter='\t',fmt='%s' )
        for result in set(id_result[:,1]):
            try:
                del original_db_dict[result]
            except:
                1
        SeqIO.write(original_db_dict.values(), "%s.%s" %(db_file,i), "fasta")

    def fo_reverse(self, xs_bar):

        numpy.set_printoptions(precision=None, threshold=None, edgeitems=None, linewidth=200, suppress=None, nanstr=None, infstr=None, formatter=None)
        self.xs_bar = xs_bar.copy()


        self.x0_bar = numpy.zeros(self.x0.shape)
        self.f_bar  = numpy.zeros(self.f.shape)
        self.p_bar  = numpy.zeros(self.p.shape)
        self.q_bar  = numpy.zeros(self.q.shape)
        self.u_bar  = numpy.zeros(self.u.shape)

        for i in range(self.M-1)[::-1]:
            h = self.ts[i+1] - self.ts[i]
            self.update_u(i)

            self.xs_bar[i,:] += self.xs_bar[i + 1, :]

            self.f_bar[:] = h*self.xs_bar[i+1, :]
            self.backend.ffcn_bar(self.ts[i:i+1],
                              self.xs[i, :], self.xs_bar[i,:],
                              self.f, self.f_bar,
                              self.p, self.p_bar,
                              self.u, self.u_bar)

            self.xs_bar[i + 1, :] = 0
            self.update_u_bar(i)

        self.x0_bar[:] += self.xs_bar[0, :]
        self.xs_bar[0, :] = 0.

    def get_occ(self, mo_energy=None, mo_coeff=None):
        '''Label the occupancies for each orbital

        Kwargs:
            mo_energy : 1D ndarray
                Obital energies

            mo_coeff : 2D ndarray
                Obital coefficients

        Examples:

        >>> from pyscf import gto, scf
        >>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
        >>> mf = scf.hf.SCF(mol)
        >>> mf.get_occ(numpy.arange(mol.nao_nr()))
        array([2, 2, 2, 2, 2, 0])
        '''
        if mo_energy is None: mo_energy = self.mo_energy
        mo_occ = numpy.zeros_like(mo_energy)
        nocc = self.mol.nelectron // 2
        mo_occ[:nocc] = 2
        if nocc < mo_occ.size:
            logger.info(self, 'HOMO = %.12g, LUMO = %.12g,',
                        mo_energy[nocc-1], mo_energy[nocc])
            if mo_energy[nocc-1]+1e-3 > mo_energy[nocc]:
                logger.warn(self, '!! HOMO %.12g == LUMO %.12g',
                            mo_energy[nocc-1], mo_energy[nocc])
        else:
            logger.info(self, 'HOMO = %.12g,', mo_energy[nocc-1])
        if self.verbose >= logger.DEBUG:
            numpy.set_printoptions(threshold=len(mo_energy))
            logger.debug(self, '  mo_energy = %s', mo_energy)
            numpy.set_printoptions()
        return mo_occ

	def general_mix_AMCMC(self, mix, **kwargs):
		"""
			simple mix data test using AMCMC
		"""
		self.mix = mix(K = self.nClass,**kwargs)
		self.mix.set_data(self.Y)
		for i in range(100):#np.int(np.ceil(0.1*self.sim))):  # @UnusedVariable
			self.mix.sample()
		self.mu_sample = list()
		self.sigma_sample = list()
		self.p_sample = list()
		
		for k in range(self.nClass):
			self.mu_sample.append(np.zeros_like(self.Thetas[k])) 
			self.sigma_sample.append(np.zeros_like(self.Sigmas[k])) 
			self.p_sample.append(np.zeros_like(self.P[k])) 
			
		self.mix.set_AMCMC(1200)
		sim_m = 2.
		for i in range(np.int(np.ceil(sim_m*self.sim))):  # @UnusedVariable
			self.mix.sample()
			for k in range(self.nClass):
				
				self.mu_sample[k] += self.mix.mu[k]/sim_m
				self.sigma_sample[k] += self.mix.sigma[k]/sim_m
				self.p_sample[k] += self.mix.p[k]	/sim_m	
		np.set_printoptions(precision=2)
			
		self.compare_class("AMCMC:")

def work_with_simple_bag_of_words():
    count = CountVectorizer()
    docs = np.array([
        'The sun is shining',
        'The weather is sweet',
        'The sun is shining and the weather is sweet',
    ])
    bag = count.fit_transform(docs)
    print(count.vocabulary_)
    print(bag.toarray())

    np.set_printoptions(precision=2)
    tfidf = TfidfTransformer(use_idf=True, norm='l2', smooth_idf=True)
    print(tfidf.fit_transform(bag).toarray())

    tf_is = 2
    n_docs = 3
    idf_is = np.log((n_docs+1) / (3+1))
    tfidf_is = tf_is * (idf_is + 1)
    print("tf-idf of term 'is' = %.2f" % tfidf_is)

    tfidf = TfidfTransformer(use_idf=True, norm=None, smooth_idf=True)
    raw_tfidf = tfidf.fit_transform(bag).toarray()[-1]
    print(raw_tfidf)

    l2_tfidf = raw_tfidf / np.sqrt(np.sum(raw_tfidf**2))
    print(l2_tfidf)

def get_response_content(fs):
    # make the laplacian matrix for the graph
    weight = 1 / fs.edge_length
    n = fs.nvertices
    L = np.zeros((n,n))
    # set the diagonal
    for i in range(n):
        L[i,i] = 2 * weight
    L[0,0] = weight
    L[-1,-1] = weight
    # complete the tridiagonal
    for i in range(n-1):
        L[i+1,i] = -weight
        L[i,i+1] = -weight
    # define other matrices
    L_pinv = np.linalg.pinv(L)
    HDH = -2*L_pinv
    v = np.diag(HDH)
    e = np.ones(n)
    D = HDH - (np.outer(v, e) + np.outer(e, v))/2
    # show some matrices
    out = StringIO()
    np.set_printoptions(linewidth=300)
    print >> out, 'Laplacian matrix:'
    print >> out, L
    print >> out
    print >> out, 'HDH:'
    print >> out, HDH
    print >> out
    print >> out, 'EDM:'
    print >> out, D
    print >> out
    return out.getvalue()

def main():
    np.set_printoptions(threshold=np.nan)
    
    feature_names = get_feature_names()
    x_train, y_train,  = get_data("training.dat")

    clf = linear_model.LinearRegression()
    clf.fit (x_train, y_train)
    
    w = clf.coef_.reshape(clf.coef_.shape[1],1) 
    y_hat_train = x_train.dot(w)

    rmse_our_train, rmse_oba_train = get_rmse(y_train, y_hat_train)
    
    x_test, y_test = get_data("test.dat")
    y_hat_test = x_test.dot(w)

    rmse_our_test, rmse_oba_test = get_rmse(y_test, y_hat_test)
 
    print "RMSE OUR Train ", rmse_our_train
    print "RMSE OBA Train ", rmse_oba_train
    print "RMSE OUR Test ", rmse_our_test
    print "RMSE OBA Test ", rmse_oba_test

    save_scatter_plot(y_train, y_hat_train, "train")
    save_scatter_plot(y_test, y_hat_test, "test")
    
    build_output_files(y_hat_train, y_hat_test, y_train, y_test)
    print_weights(w, feature_names);
    report_range(y_train)
    report_range(y_test)

    def test_0d_arrays(self):
        assert_equal(repr(np.datetime64('2005-02-25')[...]),
                     "array('2005-02-25', dtype='datetime64[D]')")

        x = np.array(1)
        np.set_printoptions(formatter={'all':lambda x: "test"})
        assert_equal(repr(x), "array(test)")

def writeOVERFLOW(filename, X, Y):
    f = open(filename, 'w')
    print 'Writing ', filename
   
    # Overflow requires 3 spanwise ndoes

    ni, nj = X.shape; nk = 3
    
    npy.set_printoptions( precision=16, threshold = ni )
    
    f.write('1'+'\n')
    f.write(str(ni) + ' ' + str(nj) + ' ' + str(nk) + '\n')
    Write3DArray(f,X)
    Write3DArray(f,X)
    Write3DArray(f,X)

    Z = npy.ones(X.shape)*0
    Write3DArray(f,Z)
    Z = npy.ones(X.shape)*-0.5
    Write3DArray(f,Z)
    Z = npy.ones(X.shape)*-1
    Write3DArray(f,Z)

    Write3DArray(f,Y)
    Write3DArray(f,Y)
    Write3DArray(f,Y)
               
    f.close()

def train_network(model, num_epochs = 100, minibatch_size = 256, lr = 0.01, mom = 0.9, wd = 0.0000):
    np.set_printoptions(linewidth=200)
    owl.set_device(owl.create_gpu_device(0))
    count = 0
    # load data
    (train_data, test_data) = imageio.load_mb_from_mat("mnist_all.mat", minibatch_size)
    num_test_samples = test_data[0].shape[0]
    (test_samples, test_labels) = map(lambda npdata : owl.from_nparray(npdata), test_data)
    for i in xrange(num_epochs):
        print "---Epoch #", i
        for (mb_samples, mb_labels) in train_data:
            num_samples = mb_samples.shape[0]
            data = owl.from_nparray(mb_samples).reshape([28, 28, 1, num_samples])
            label = owl.from_nparray(mb_labels)
            out, weightgrad, biasgrad = train(model, data, label)
            for k in range(len(model.weights)):
                model.weightdelta[k] = mom * model.weightdelta[k] - lr / num_samples * weightgrad[k] - wd * model.weights[k]
                model.biasdelta[k] = mom * model.biasdelta[k] - lr / num_samples * biasgrad[k]
                model.weights[k] += model.weightdelta[k]
                model.bias[k] += model.biasdelta[k]

            count = count + 1
            if (count % 1) == 0:
                print_training_accuracy(out, label, num_samples)
            if count == 100:
                sys.exit()