def dynamically_bin(hic1, hic2, chrom, binbounds):
    unbinned1, map1 = hic1.cis_heatmap(chrom, start=binbounds[0, 0], stop=binbounds[-1, 1], datatype='fend',
                                 arraytype='full', returnmapping=True)
    unbinned2, map2 = hic2.cis_heatmap(chrom, start=binbounds[0, 0], stop=binbounds[-1, 1], datatype='fend',
                                 arraytype='full', returnmapping=True)
    map1[:, 2] = (map1[:, 0] + map1[:, 1])
    map2[:, 2] = (map2[:, 0] + map2[:, 1])
    allmap = numpy.vstack((map1, map2))
    allmap = allmap[numpy.argsort(allmap[:, 2]), :]
    indices1 = numpy.searchsorted(allmap[:, 2], map1[:, 2])
    indices1_1 = (indices1.reshape(-1, 1) * allmap.shape[0] + indices1.reshape(1, -1)).ravel()
    indices2 = numpy.searchsorted(allmap[:, 2], map2[:, 2])
    indices2_1 = (indices2.reshape(-1, 1) * allmap.shape[0] + indices2.reshape(1, -1)).ravel()
    unbinned = numpy.zeros((allmap.shape[0], allmap.shape[0], 2), dtype=numpy.float32)
    unbinned[:, :, 0] += numpy.bincount(indices1_1, minlength=allmap.shape[0] ** 2,
                                        weights=unbinned1[:, :, 0].ravel()).reshape(allmap.shape[0], -1)
    unbinned[:, :, 1] += numpy.bincount(indices1_1, minlength=allmap.shape[0] ** 2,
                                        weights=unbinned1[:, :, 1].ravel()).reshape(allmap.shape[0], -1)
    unbinned[:, :, 0] += numpy.bincount(indices2_1, minlength=allmap.shape[0] ** 2,
                                        weights=unbinned2[:, :, 0].ravel()).reshape(allmap.shape[0], -1)
    unbinned[:, :, 1] += numpy.bincount(indices2_1, minlength=allmap.shape[0] ** 2,
                                        weights=unbinned2[:, :, 1].ravel()).reshape(allmap.shape[0], -1)
    indices = numpy.triu_indices(allmap.shape[0], 1)
    unbinned = unbinned[indices[0], indices[1], :]
    binned, binmap = hic1.cis_heatmap(chrom, binbounds=binbounds, datatype='fend', arraytype='full',
                                      returnmapping=True)
    binned += hic2.cis_heatmap(chrom, binbounds=binbounds, datatype='fend', arraytype='full')
    indices = numpy.triu_indices(binbounds.shape[0], 1)
    upper = binned[indices[0], indices[1], :]
    hifive.hic_binning.dynamically_bin_cis_array(unbinned, allmap, upper, binmap,
                                                 expansion_binsize=0, minobservations=25)
    binned[indices[0], indices[1], :] = upper
    binned[indices[1], indices[0], :] = upper
    return binned

    def mvn_msr(self, corrDS=None, abstol=1e-12, reltol=1e-12, intLb=-10, intUb=10):
        systype = self.systype
        beta = self.beta
        nls = len(self.comps)
        if corrDS is None:
            correl = self.syscorrDS[np.triu_indices(nls, 1)]
        else:
            correl = corrDS[np.triu_indices(nls, 1)]
        if corrDS is None:
            corrDS = self.syscorrDS
        i = 1
        n = 10000
        syspf0 = 0.0
        dpf = 1.0
        # while i!=0:
        # n +=10000
        # v,res,i = stats.mvn.mvndst(intLb*np.ones(nls), beta, np.zeros(nls, dtype=int), correl, [nls*n,1e-12, 1e-12])
        while i != 0:
            n += 10000
            res, i = stats.mvn.mvnun(-10 * np.ones(nls), beta, np.zeros(nls), corrDS, [nls * n, abstol, reltol])
        # if abs(res-res1)/(0.5*(res+res1))>1e-3:
        # print 'warning: abnormal difference between mvnun and mvndst results'
        if systype.lower() == "series":
            syspf = 1.0 - res
            sysbeta = -stats.norm.ppf(syspf)
            results = ReliabilityResults(sysbeta, syspf)
        elif systype.lower() == "parallel":
            syspf = res
            sysbeta = -stats.norm.ppf(syspf)
            results = ReliabilityResults(sysbeta, syspf)
        else:
            print ("mvn_msr only supports series or parallel system")
            sys.exit(0)

        return results

 def normalization(self):
     """
     Normalize the equilibrium steady state correlations according to
     Eq 76 in Lorenzo's writeup
     """
     N = self.latsize
     #First disconnect
     self.disconnect(self.steady_state)
     norm_1 = N+np.sum(self.steady_state[2*N:3*N])
     sxxpsyy = self.steady_state[3*N:].reshape(3,3,N,N)[0,0,:,:] +\
         self.steady_state[3*N:].reshape(3,3,N,N)[1,1,:,:]
     sxymsyx = self.steady_state[3*N:].reshape(3,3,N,N)[0,1,:,:] -\
         self.steady_state[3*N:].reshape(3,3,N,N)[1,0,:,:]
     norms = []
     for kvec in self.kvecs:
         argmat = np.zeros((N,N))
         for (m,n) in combinations(np.arange(N),2):
             argmat[m,n] = kvec.dot(self.atoms[m].coords-self.atoms[n].coords)
         norm_2 = np.sum(\
                 np.cos(argmat[np.triu_indices(N, k=1)]) *\
                     sxxpsyy[np.triu_indices(N, k=1)] +\
                 np.sin(argmat[np.triu_indices(N, k=1)]) *\
                     sxymsyx[np.triu_indices(N, k=1)])
         norms.append(0.5*(norm_1+norm_2))
     #Reconnect before exit
     self.reconnect(self.steady_state)
     return np.array(norms).flatten()

def get_query_clusters(points, k):
    '''
    points [n,m] - array for n points with dimention m - encoded query
    '''
    # normalize input
    points = normalize(points.astype(np.float))
    # get similarity matrix (cosine distance)
    dist = points.dot(points.T)
    # initialize variables
    n_pt = len(points)
    cluster_old, cluster_new = np.ones(n_pt), np.zeros(n_pt)
    # special case, no clustering
    if k==1 or n_pt==1:
        return np.zeros(n_pt), 1 if n_pt==1 else np.mean(dist[np.triu_indices(n_pt,k=1)])
    # randomly choose k starting centroids
    centroids = points[np.random.permutation(n_pt)[:k]]
    while not np.array_equal(cluster_old, cluster_new):
        cluster_old = cluster_new
        # get cluster index for each point
        cluster_new = np.argmax(points.dot(centroids.T), axis=1)
        # get new centroids, and within class mean distance/similarity
        centroids, in_dist = [], []
        for c in np.unique(cluster_new):
            pid = cluster_new==c
            # set new centroid as the one who has minimum total distance to rest of the points in the cluster
            cid = np.argmax(np.sum(dist[np.ix_(pid, pid)], axis=1))
            centroids.append(points[pid][cid])
            in_dist.append(1 if sum(pid)==1 else np.mean(dist[np.ix_(pid,pid)][np.triu_indices(sum(pid),k=1)]))
        centroids = np.array(centroids)
        # traditional way to get new centroid, not working well for cosine distance
#         centroids = normalize([np.mean(points[cluster_new==c], axis=0) for c in np.unique(cluster_new)])

    return cluster_new, np.mean(in_dist), centroids

def example_one():
    """
    Generates a set of sample data for the
    examples page of the hetaira web tool.
    """
    
    np.random.seed(5)
    ids = ['Pr'] + list(ascii_lowercase) + ['Sp']
    
    # make some data where all activities are the same
    data = np.ones((26,26))
    
    # make some random activites to pull from
    y = np.random.uniform(1000, 2500, (26,26))
    
    # this will replace the ones with numbers from the uniform
    # distribution, increasing by one at each column
    # using the upper triangular matrix
    data[np.triu_indices(26)] = y[np.triu_indices(26)]

    # stack a perfectly promiscuous and a perfectly (almost)
    # specific column on either side of the data
    data = np.hstack((np.full((26,1), 1e-10), data, np.ones((26,1))))
    data[0,0] = 100
    descriptors = None
    example = Promiscuity(ids, np.fliplr(data), descriptors)
    return example.hetaira_results()

def test_syrk(k, n, dtype, rng):
    tols = tolerances[dtype]

    A = np.zeros((n, k), dtype=dtype)
    C = np.zeros((n, n), dtype=dtype)
    D = np.zeros((k, k), dtype=dtype)
    A[...] = rng.uniform(-1, 1, size=A.shape)
    C[...] = rng.uniform(-1, 1, size=C.shape)
    D[...] = rng.uniform(-1, 1, size=D.shape)

    clA, clC, clD = map(to_ocl, [A, C, D])
    a = 0.9
    b = 0.5

    try:
        blas.setup()

        # normal syrk
        up = np.triu_indices(n)
        event = blas.syrk(queue, clA, clC, alpha=a, beta=b)
        assert np.allclose(clC.get()[up], (a*np.dot(A, A.T) + b*C)[up], **tols)
        assert isinstance(event, cl.Event)

        # transposed syrk
        up = np.triu_indices(k)
        blas.syrk(queue, clA, clD, transA=True, alpha=a, beta=b)
        assert np.allclose(clD.get()[up], (a*np.dot(A.T, A) + b*D)[up], **tols)
    finally:
        blas.teardown()

    def loss_function(mapping12):
        """Computes the loss function of a given mapping.

        Using the graph kernel of two sets of distance.
        """
        
        global tractography1, tractography2
        global dm1_all, dm1_all_small, dm2_all, dm2_all_small
        global kdt1, kdt2
        k = 10
        
        radius = 150  
        loss = 0.0
        for sid in np.arange(len(tractography1)):               
            #idx1 = kdt1.query_radius(dm1_all_small[sid], radius)[0]
            idx1 = kdt1.query(dm1_all_small[sid], k)[1][0]            
            dm_small1 = dm1_all[idx1][:,idx1]
            e1 = dm_small1[np.triu_indices(dm_small1.shape[0],1)]
                    
            #idx2 = kdt2.query_radius(dm2_all_small[mapping12[sid]], radius)[0]
            idx2 = kdt2.query(dm2_all_small[mapping12[sid]], k)[1][0]
            dm_small2 = dm2_all[idx2][:,idx2]
            e2 = dm_small2[np.triu_indices(dm_small2.shape[0],1)]
            
            #loss = loss + Graph_KN(e1, e2, weight=1., num_bins = 128)
            #similarity = similarity + Pyramid_KN(e1, e2, weight=1., num_bins = 128)
            loss = loss +  Pyramid_KN(e1, e2, weight=1., num_bins = 128)
            
        return loss

def sim_matrix_within_group_means(matrix, n1):
    """
    Computes the mean of the upper triangle (k=1) for the blocks
    (0,n-1)*(0,n-1) and (n,2n-1)*(n,2n-1), and their difference
    (for convenience).

    Parameters
    ----------
    matrix : 2D symmetric numpy array
        1 or 2 dimensional numpy array, the n1 first indices in
        the zeroth axis of the array, should correspond to the
        values of the first group.
        The value of ``matrix[i][j]`` should correspond to
    n1 : int
        the number of elements in the first group

    Returns
    -------
    mean1 : float
        the average similarity between members in the first group
    mean2: float
        the average similarity between members in the second group
    mean1-mean2: float
        just mean1-mean2 (as a convenience for stat. testing)
    """
    n2 = matrix.shape[0] - n1
    indices1 = np.triu_indices(n1, k=1)
    indices2base = np.triu_indices(n2, k=1)
    indices2I = indices2base[0].copy() + n1
    indices2J = indices2base[1].copy() + n1
    indices2 = (indices2I, indices2J)
    mean1 = np.average(matrix[indices1])
    mean2 = np.average(matrix[indices2])
    return mean1, mean2, mean1 - mean2

    def test_triu_indices(self):
        iu1 = triu_indices(4)
        iu2 = triu_indices(4, 2)

        a = np.array([[1, 2, 3, 4],
                      [5, 6, 7, 8],
                      [9, 10, 11, 12],
                      [13, 14, 15, 16]])

        # Both for indexing:
        yield (assert_array_equal, a[iu1],
               array([1, 2,  3,  4,  6, 7, 8, 11, 12, 16]))

        # And for assigning values:
        a[iu1] = -1
        yield (assert_array_equal, a,
               array([[-1, -1, -1, -1],
                      [ 5, -1, -1, -1],
                      [ 9, 10, -1, -1],
                      [13, 14, 15, -1]])  )

        # These cover almost the whole array (two diagonals right of the main one):
        a[iu2] = -10
        yield ( assert_array_equal, a,
                array([[ -1,  -1, -10, -10],
                       [  5,  -1,  -1, -10],
                       [  9,  10,  -1,  -1],
                       [ 13,  14,  15,  -1]]) )

def calcH2Continuous(XXT, phe, keepArr, prev, h2coeff):
	t = stats.norm(0,1).isf(prev)
	phit = stats.norm(0,1).pdf(t)
	
	K1 = 1 - prev
	K2 = 1 - K1
	P = np.sum(phe<t) / float(phe.shape[0])	
	P2 = 1.0
	P1 = K2*P2*P / (K1*(1-P))
	R = P2 / P1
	
	XXT = XXT[np.ix_(keepArr, keepArr)]
	phe = phe[keepArr]
	
	xCoeff = (((R-1)*phit*t + K1 + R*K2)**2 * (K1+R*K2)**2 - ((R-1)*phit)**4) / (K1 + R*K2)**4
	x = (xCoeff * h2coeff) * XXT 
	pheMean = 0
	pheVar = 1	
	y = np.outer((phe-pheMean) / np.sqrt(pheVar), (phe-pheMean)/np.sqrt(pheVar))
	y -= ((R-1)*phit / (K1+R*K2))**2
	
	y = y[np.triu_indices(y.shape[0], 1)]
	x = x[np.triu_indices(x.shape[0], 1)]
	
	slope, intercept, rValue, pValue, stdErr = stats.linregress(x,y)
	return slope

    def __init__(self, endmembers, alphas, energy_interaction, volume_interaction=None, entropy_interaction=None):

        self.n_endmembers = len(endmembers)

        # Create array of van Laar parameters
        self.alphas = np.array(alphas)

        # Create 2D arrays of interaction parameters
        self.We = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
        self.We[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in energy_interaction
                                                                for i in row])

        if entropy_interaction is not None:
            self.Ws = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
            self.Ws[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in entropy_interaction
                                                                        for i in row])
        else:
            self.Ws = np.zeros((self.n_endmembers, self.n_endmembers))

        if volume_interaction is not None:
            self.Wv = np.triu(2. / (self.alphas[:, np.newaxis] + self.alphas), 1)
            self.Wv[np.triu_indices(self.n_endmembers, 1)] *= np.array([i for row in volume_interaction
                                                                        for i in row])
        else:
            self.Wv = np.zeros((self.n_endmembers, self.n_endmembers))


        # initialize ideal solution model
        IdealSolution.__init__(self, endmembers)

def scoring2B_behavior():
	t_clusters = np.zeros((600,3))
	t_clusters[0:200,0] = 1
	t_clusters[200:400,1] = 1
	t_clusters[400:,2] = 1
	t_ccm = np.dot(t_clusters,t_clusters.T)

	n_uniq = len(np.triu_indices(t_ccm.shape[0],k=1)[0])
	res = []
	concentrations = [1000,100,50,25,10,5,3,1]
	for c in concentrations:
		for i in range(50):
			ccm = np.copy(t_ccm)
			ccm[np.triu_indices(t_ccm.shape[0],k=1)] -= np.random.beta(1,c,n_uniq)
			#ccm[np.tril_indices(t_ccm.shape[0],k=-1)] = ccm[np.triu_indices(t_ccm.shape[0],k=1)]
			ccm[np.tril_indices(t_ccm.shape[0],k=-1)] = 0
			ccm = ccm + ccm.T
			np.fill_diagonal(ccm,1)
			ccm = np.abs(ccm)
			res.append([c,calculate2(ccm,t_ccm)])
	res = [map(str,x) for x in res]
	res = ['\t'.join(x) for x in res]
	f = open('scoring2B_beta.tsv', 'w')
	f.write('\n'.join(res))
	f.close()	

def calcH2Binary(XXT, phe, probs, thresholds, keepArr, prev, h2coeff):
	K = prev
	P = np.sum(phe>0) / float(phe.shape[0])
	
	XXT = XXT[np.ix_(keepArr, keepArr)]
	phe = phe[keepArr]
	
	if (thresholds is None):
		t = stats.norm(0,1).isf(K)
		phit = stats.norm(0,1).pdf(t)
		xCoeff = P*(1-P) / (K**2 * (1-K)**2) * phit**2 * h2coeff
		y = np.outer((phe-P) / np.sqrt(P*(1-P)), (phe-P) / np.sqrt(P*(1-P)))
		x = xCoeff * XXT
		
	else:
		probs = probs[keepArr]
		thresholds = thresholds[keepArr]
		Ki = K*(1-P) / (P*(1-K)) * probs / (1 + K*(1-P) / (P*(1-K))*probs - probs)
		phit = stats.norm(0,1).pdf(thresholds)	
		probsInvOuter = np.outer(probs*(1-probs), probs*(1-probs))
		y = np.outer(phe-probs, phe-probs) / np.sqrt(probsInvOuter)	
		sumProbs = np.tile(np.column_stack(probs).T, (1,probs.shape[0])) + np.tile(probs, (probs.shape[0], 1))
		Atag0 = np.outer(phit, phit) * (1 - (sumProbs)*(P-K)/(P*(1-K)) + np.outer(probs, probs)*(((P-K)/(P*(1-K)))**2)) / np.sqrt(probsInvOuter)
		B0 = np.outer(Ki + (1-Ki)*(K*(1-P))/(P*(1-K)), Ki + (1-Ki)*(K*(1-P))/(P*(1-K)))
		x = (Atag0 / B0 * h2coeff) * XXT	
	
	y = y[np.triu_indices(y.shape[0], 1)]
	x = x[np.triu_indices(x.shape[0], 1)]
	
	slope, intercept, rValue, pValue, stdErr = stats.linregress(x,y)
	return slope

def calcH2Continuous_twotails(XXT, phe, keepArr, prev, h2coeff):

	print 'computing h2 for a two-tails ascertained study...'
	
	XXT = XXT[np.ix_(keepArr, keepArr)]
	phe = phe[keepArr]	

	t1 = stats.norm(0,1).ppf(prev)
	t2 = stats.norm(0,1).isf(prev)
	phit1 = stats.norm(0,1).pdf(t1)
	phit2 = stats.norm(0,1).pdf(t2)
	
	K1 = prev
	K2 = prev
	
	xCoeff = ((phit2*t2 - phit1*t1 + K1 + K2)**2 * (K1+K2)**2 - (phit2-phit1)**4) / (K1 + K2)**4
	intersect = ((phit2-phit1) / (K1+K2))**2
		
	pheMean = 0
	pheVar = 1
	
	x = (xCoeff * h2coeff) * XXT 
	y = np.outer((phe-pheMean)/np.sqrt(pheVar), (phe-pheMean)/np.sqrt(pheVar))
	y -= intersect
	
	y = y[np.triu_indices(y.shape[0], 1)]
	x = x[np.triu_indices(x.shape[0], 1)]
	
	slope, intercept, rValue, pValue, stdErr = stats.linregress(x,y)	
	return slope

    def loss_function2(mapping12):
        """Computes the loss function of a given mapping.

        Using the graph kernel of two sets of distance.
        """
        global dis_1, dis_2
        global kdt_1, kdt_2
        global dm1, dm2, dm1_all, dm2_all
        
        k = 15
        
        radius = 100  
        similarity = 0.0
        for sid in np.arange(len(pro_1)):               
            idx1 = kdt_1.query_radius(dm1[sid], radius)[0]
            #idx1 = kdt_1.query(dm1[sid], k)[1][0]            
            dm_small1 = dm1_all[idx1][:,idx1]
            e1 = dm_small1[np.triu_indices(dm_small1.shape[0],1)]
                    
            idx2 = kdt_2.query_radius(dis_2[mapping12[sid]], radius)[0]
            #idx2 = kdt_2.query(dis_2[mapping12[sid]], k)[1][0]
            dm_small2 = dm2_all[idx2][:,idx2]
            e2 = dm_small2[np.triu_indices(dm_small2.shape[0],1)]
            
            #loss = loss + Graph_KN(e1, e2, weight=1., num_bins = 128)
            similarity = similarity + Pyramid_KN(e1, e2, weight=1., num_bins = 128)
            
        return 1./similarity

def insertCartesainData(gt_row, xyzuvw_mean, xyzuvw_cov):
    dim = 6
    # CART_COL_NAMES = ['X', 'Y', 'Z', 'U', 'V', 'W',
    #                   'dX', 'dY', 'dZ', 'dU', 'dV', 'dW',
    #                   'c_XY', 'c_XZ', 'c_XU', 'c_XV', 'c_XW',
    #                           'c_YZ', 'c_YU', 'c_YV', 'c_YW',
    #                                   'c_ZU', 'c_ZV', 'c_ZW',
    #                                           'c_UV', 'c_UW',
    #                                                   'c_VW']

    # fill in cartesian mean
    try:
        for col_ix, col_name in enumerate(CART_COL_NAMES[:6]):
            gt_row[col_name] = xyzuvw_mean[col_ix]
    except IndexError:
        import pdb; pdb.set_trace()

    # fill in standard deviations
    xyzuvw_stds = np.sqrt(xyzuvw_cov[np.diag_indices(dim)])
    for col_ix, col_name in enumerate(CART_COL_NAMES[6:12]):
        gt_row[col_name] = xyzuvw_stds[col_ix]

    correl_matrix = xyzuvw_cov / xyzuvw_stds / xyzuvw_stds.reshape(6, 1)
    # fill in correlations
    for col_ix, col_name in enumerate(CART_COL_NAMES[12:]):
        gt_row[col_name] = correl_matrix[
            np.triu_indices(dim, k=1)[0][col_ix],
            np.triu_indices(dim, k=1)[1][col_ix]
        ]

    def __init__(self, eta=None, n=None, p=None, transform='interval', *args, **kwargs):
        if (p is not None) and (n is not None) and (eta is None):
            warnings.warn('Parameters to LKJCorr have changed: shape parameter n -> eta '
                          'dimension parameter p -> n. Please update your code. '
                          'Automatically re-assigning parameters for backwards compatibility.',
                          DeprecationWarning)
            self.n = p
            self.eta = n
            eta = self.eta
            n = self.n
        elif (n is not None) and (eta is not None) and (p is None):
            self.n = n
            self.eta = eta
        else:
            raise ValueError('Invalid parameter: please use eta as the shape parameter and '
                             'n as the dimension parameter.')

        n_elem = int(n * (n - 1) / 2)
        self.mean = np.zeros(n_elem, dtype=theano.config.floatX)

        if transform == 'interval':
            transform = transforms.interval(-1, 1)

        super(LKJCorr, self).__init__(shape=n_elem, transform=transform,
                                      *args, **kwargs)
        warnings.warn('Parameters in LKJCorr have been rename: shape parameter n -> eta '
                      'dimension parameter p -> n. Please double check your initialization.',
                      DeprecationWarning)
        self.tri_index = np.zeros([n, n], dtype='int32')
        self.tri_index[np.triu_indices(n, k=1)] = np.arange(n_elem)
        self.tri_index[np.triu_indices(n, k=1)[::-1]] = np.arange(n_elem)

    def find_stationary_var(amat=None, bmat=None, cmat=None):
        """Find fixed point of H = CC' + AHA' + BHB' given A, B, C.

        Parameters
        ----------
        amat, bmat, cmat : (nstocks, nstocks) arrays
            Parameter matrices

        Returns
        -------
        (nstocks, nstocks) array
            Unconditional variance matrix

        """
        nstocks = amat.shape[0]
        kwargs = {'amat': amat, 'bmat': bmat, 'ccmat': cmat.dot(cmat.T)}
        fun = partial(ParamGeneric.fixed_point, **kwargs)
        try:
            with np.errstate(divide='ignore', invalid='ignore'):
                hvar = np.eye(nstocks)
                sol = sco.fixed_point(fun, hvar[np.tril_indices(nstocks)])
                hvar[np.tril_indices(nstocks)] = sol
                hvar[np.triu_indices(nstocks, 1)] \
                    = hvar.T[np.triu_indices(nstocks, 1)]
                return hvar
        except RuntimeError:
            # warnings.warn('Could not find stationary varaince!')
            return None

def angles_and_contour_lengths(spline,deriv,
                               min_change_px=0,max_change_px=np.inf):
    """
    gets Cos(Theta(i)) and L(i), where i runs along the spline order given,
    and L is the contour length between segments chosen at index i

    Args:
        spline: tuple of x_spline,y_spline -- x and y values of the line, size N
        deriv: the continuous derivative of spline, size N
        <min/max>_change_px: the minimum and maximum pixel changes
    Returns:
        tuple of angle_info object, L0_px
    """
    # get the x and y coordinates of the spline
    x_spline,y_spline = spline
    x_deriv,y_deriv = deriv
    deriv_unit_vector = np.array((x_deriv,y_deriv))
    deriv_unit_vector /= np.sqrt(np.sum(np.abs(deriv_unit_vector**2),axis=0))
    assert ((np.sum(deriv_unit_vector**2,axis=0) -1) < 1e-6).all() , \
        "Unit vectors not correct"
    # POST: unit vector are normalized, |v| = 1
    dx_spline = np.array([0] + list(np.diff(x_spline)))
    dy_spline = np.array([0] + list(np.diff(y_spline)))
    # d_spline(i) is the change from i-i to i (zero if i=0)
    d_spline = np.sqrt(dx_spline**2 + dy_spline**2)
    assert (dx_spline <= d_spline).all()
    contour_lengths = np.cumsum(d_spline)
    L0 = contour_lengths[-1]
    n = x_spline.shape[0]
    contour_length_matrix = _difference_matrix(contour_lengths,contour_lengths)
    dx_deriv = deriv_unit_vector[0, :]
    dy_deriv = deriv_unit_vector[1, :]
    angle2 = np.arctan2(dy_deriv, dx_deriv)
    angle_diff_matrix = _difference_matrix(angle2.T, angle2.T)
    # normalize to 0 to 2*pi
    where_le_0 = np.where(angle_diff_matrix < 0)
    angle_diff_matrix[where_le_0] += 2 * np.pi
    assert ((angle_diff_matrix >= 0) & (angle_diff_matrix <= 2*np.pi)).all()
    # POST: angles calculated correctly...
    # only look at the upper triangular part
    idx_upper_tri = np.triu_indices(n)
    idx_upper_tri_no_diag =np.triu_indices(n,k=1)
    # upper diagonal should have >0 contour length
    assert (contour_length_matrix[idx_upper_tri_no_diag] > 0).all() , \
        "Contour lengths should be positive"
    # POST: contour lengths and angles make sense; we only want upper triangular
    # (*including* the trivial 0,0 point along the diagonal)
    contour_length_matrix_check_valid = contour_length_matrix[idx_upper_tri]
    # POST: matrix is filled in, determine where the value are valid
    ok_idx = np.where( (contour_length_matrix_check_valid > min_change_px) &
                       (contour_length_matrix_check_valid < max_change_px))
    sanit = lambda x: x[idx_upper_tri][ok_idx].flatten()
    sort_idx = np.argsort(sanit(contour_length_matrix))
    sanit_and_sort = lambda x: sanit(x)[sort_idx]
    # return everything sorted as per sort_idx
    flat_L = sanit_and_sort(contour_length_matrix)
    flat_angle = np.arccos(np.cos(sanit_and_sort(angle_diff_matrix)))
    to_ret = angle_info(theta=flat_angle, L_px=flat_L)
    return to_ret,L0

def coefs2mats(coefs, n=8):
    const = coefs[0]
    jac = coefs[1:n+1]
    hes = np.zeros((n,n))
    hes[np.triu_indices(n)] = hes.T[np.triu_indices(n)] = coefs[n+1:]
    
    hes[np.diag_indices(n)] *= 2
    return const, jac, hes