Python scipy.zeros() Examples

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm
    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp = airtemp + 237.3
        b = 0.6108*(math.exp((17.27*airtemp)/temp))
        delta = (4098*b)/(temp**2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp = airtemp[i] + 237.3
            b = 0.6108*(math.exp(17.27*airtemp[i])/temp)
            delta[i] = (4098*b)/(temp**2)
    return delta 

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm

    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp = airtemp + 237.3
        b = 0.6108 * (math.exp((17.27 * airtemp) / temp))
        delta = (4098 * b) / (temp ** 2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp = airtemp[i] + 237.3
            b = 0.6108 * (math.exp(17.27 * airtemp[i]) / temp)
            delta[i] = (4098 * b) / (temp ** 2)
    return delta 

def gap(data, refs=None, nrefs=20, ks=range(1,11), method=None):
    shape = data.shape
    if refs is None:
        tops = data.max(axis=0)
        bots = data.min(axis=0)
        dists = scipy.matrix(scipy.diag(tops-bots))

        rands = scipy.random.random_sample(size=(shape[0], shape[1], nrefs))
        for i in range(nrefs):
            rands[:, :, i] = rands[:, :, i]*dists+bots
    else:
        rands = refs
    gaps = scipy.zeros((len(ks),))
    for (i, k) in enumerate(ks):
        g1 = method(n_clusters=k).fit(data)
        (kmc, kml) = (g1.cluster_centers_, g1.labels_)
        disp = sum([euclidean(data[m, :], kmc[kml[m], :]) for m in range(shape[0])])

        refdisps = scipy.zeros((rands.shape[2],))
        for j in range(rands.shape[2]):
            g2 = method(n_clusters=k).fit(rands[:, :, j])
            (kmc, kml) = (g2.cluster_centers_, g2.labels_)
            refdisps[j] = sum([euclidean(rands[m, :, j], kmc[kml[m],:]) for m in range(shape[0])])
        gaps[i] = scipy.log(scipy.mean(refdisps))-scipy.log(disp)
    return gaps 

def coupling_optim_garrick(y,t):
	creation=s.zeros(n_bin)
	destruction=s.zeros(n_bin)
	#now I try to rewrite this in a more optimized way
	destruction = -s.dot(s.transpose(kernel),y)*y #much more concise way to express\
	#the destruction of k-mers 
	
	for k in xrange(n_bin):
		kyn = (kernel*f_garrick[:,:,k])*y[:,s.newaxis]*y[s.newaxis,:]
		creation[k] = s.sum(kyn)
	creation=0.5*creation
	out=creation+destruction
	return out



#Now I work with the function for espressing smoluchowski equation when a uniform grid is used 

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm
    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp = airtemp + 237.3
        b = 0.6108*(math.exp((17.27*airtemp)/temp))
        delta = (4098*b)/(temp**2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp = airtemp[i] + 237.3
            b = 0.6108*(math.exp(17.27*airtemp[i])/temp)
            delta[i] = (4098*b)/(temp**2)
    return delta 

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm
    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp = airtemp + 237.3
        b = 0.6108*(math.exp((17.27*airtemp)/temp))
        delta = (4098*b)/(temp**2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp = airtemp[i] + 237.3
            b = 0.6108*(math.exp(17.27*airtemp[i])/temp)
            delta[i] = (4098*b)/(temp**2)
    return delta 

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm
    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp = airtemp + 237.3
        b = 0.6108*(math.exp((17.27*airtemp)/temp))
        delta = (4098*b)/(temp**2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp = airtemp[i] + 237.3
            b = 0.6108*(math.exp(17.27*airtemp[i])/temp)
            delta[i] = (4098*b)/(temp**2)
    return delta 

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm
    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp = airtemp + 237.3
        b = 0.6108*(math.exp((17.27*airtemp)/temp))
        delta = (4098*b)/(temp**2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp = airtemp[i] + 237.3
            b = 0.6108*(math.exp(17.27*airtemp[i])/temp)
            delta[i] = (4098*b)/(temp**2)
    return delta 

def coupling_optim(y,t):
	creation=s.zeros(n_bin)
	destruction=s.zeros(n_bin)
	#now I try to rewrite this in a more optimized way
	destruction = -s.dot(s.transpose(kernel),y)*y #much more concise way to express\
	#the destruction of k-mers 
	kyn = kernel*y[:,s.newaxis]*y[s.newaxis,:]
	for k in xrange(n_bin):
		creation[k] = s.sum(kyn[s.arange(k),k-s.arange(k)-1])
	creation=0.5*creation
	out=creation+destruction
	return out


#Now I go for the optimal optimization of the chi_{i,j,k} coefficients used by Garrick for
# dealing with a non-uniform grid. 

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm
    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp_kelvin = airtemp + 237.3
        b = 0.6108*(math.exp((17.27*airtemp)/temp_kelvin))
        delta = (4098*b)/(temp_kelvin**2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp_kelvin = airtemp[i] + 237.3
            b = 0.6108*(math.exp(17.27*airtemp[i])/temp_kelvin)
            delta[i] = (4098*b)/(temp_kelvin**2)
    return delta 

def delta_calc(airtemp):
    """
    Calculates slope of saturation vapour pressure curve at air temperature [kPa/Celsius]
    http://www.fao.org/docrep/x0490e/x0490e07.htm
    :param airtemp: Temperature in Celsius
    :return: slope of saturation vapour pressure curve [kPa/Celsius]
    """
    l = sp.size(airtemp)
    if l < 2:
        temp = airtemp + 237.3
        b = 0.6108*(math.exp((17.27*airtemp)/temp))
        delta = (4098*b)/(temp**2)
    else:
        delta = sp.zeros(l)
        for i in range(0, l):
            temp = airtemp[i] + 237.3
            b = 0.6108*(math.exp(17.27*airtemp[i])/temp)
            delta[i] = (4098*b)/(temp**2)
    return delta 

def calculate_daily_extraterrestrial_irradiation(doy, latitude):
    lat = latitude
    l = np.size(doy)
    s = 0.0820  # MJ m-2 min-1
    lat_rad = lat * (math.pi / 180)
    if l < 2:
        day = doy
        dr = 1 + (0.033 * math.cos((2 * math.pi * day) / 365))  # inverse relative distance Earth-Sun
        dt = 0.409 * math.sin(((2 * math.pi * day) / 365) - 1.39)  # solar declination in radian
        ws = math.acos(-math.tan(lat_rad) * math.tan(dt))   # sunset hour angle in radian
        rext = ((24* 60) / math.pi) * s * dr * ((ws * math.sin(lat_rad) * math.sin(dt)) + (math.cos(lat_rad) * math.cos(dt) * math.sin(ws)))  # MJm-2day-1
    else:
        rext = np.zeros(l)
        for i in range(0, l):
            day = doy[i]
            dr = 1 + (0.033 * math.cos((2 * math.pi * day) / 365))  # inverse relative distance Earth-Sun
            dt = 0.409 * math.sin(((2 * math.pi * day) / 365) - 1.39)  # solar declination in radian
            ws = math.acos(-math.tan(lat_rad) * math.tan(dt))  # sunset hour angle in radian
            rext[i] = ((24 * 60) / math.pi) * s * dr * ((ws * math.sin(lat_rad) * math.sin(dt)) + (math.cos(lat_rad) * math.cos(dt) * math.sin(ws)))  # MJm-2day-1
    return rext 

def impz(b,a):
    """Pseudo implementation of the impz method of MATLAB"""
#% Compute time vector
# M = 0;  NN = [];
# if isempty(N)
#   % if not specified, determine the length
#   if isTF
#     N = impzlength(b,a,.00005);
#   else
#     N  = impzlength(b,.00005);
#   end
    p = np.roots(a)
    N = stableNmarginal_length(p, 0.00005, 0)
    N = len(b) * len(b) * len(b) # MATLAB AUTOFINDS THE SIZE HERE... 
    #TODO: Implement some way of finding the autosieze of this... I used a couple of examples... matlab gave 43 as length we give 64
    x = zeros(N)
    x[0] = 1
    h = lfilter(b,a, x)
    return h 

def __MR_W_D_matrix(self,img,labels):
        s = sp.amax(labels)+1
        vect = self.__MR_superpixel_mean_vector(img,labels)
        
        adj = self.__MR_get_adj_loop(labels)
        
        W = sp.spatial.distance.squareform(sp.spatial.distance.pdist(vect))
        
        W = sp.exp(-1*W / self.weight_parameters['delta'])
        W[adj.astype(np.bool)] = 0
        

        D = sp.zeros((s,s)).astype(float)
        for i in range(s):
            D[i, i] = sp.sum(W[i])

        return W,D 

def star_graph_lumped(N, tau, gamma, I0, tmin, tmax, tcount):
    times = scipy.linspace(tmin, tmax, tcount)
    #    [[central node infected] + [central node susceptible]]
    #X = [Y_1^1, Y_1^2, ..., Y_1^{N}, Y_2^0, Y_2^1, ..., Y_2^{N-1}]
    X0 = scipy.zeros(2*N)  #length 2*N of just 0 entries
    X0[I0]=I0*1./N #central infected, + I0-1 periph infected prob
    X0[N+I0] = 1-I0*1./N #central suscept + I0 periph infected
    X = EoN.my_odeint(star_graph_dX, X0, times, args = (tau, gamma, N))
    #X looks like [[central susceptible,k periph] [ central inf, k-1 periph]] x T

    central_inf = X[:,:N]
    central_susc = X[:,N:]

    I = scipy.array([ sum(k*central_susc[t][k] for k in range(N))
            + sum((k+1)*central_inf[t][k] for k in range(N))
            for t in range(len(X))])
    S = N-I
    return times, S, I 

def _flat_field(X, uniformity_thresh):
    """."""

    Xhoriz = _low_frequency_horiz(X, sigma=4.0)
    Xhorizp = _low_frequency_horiz(X, sigma=3.0)
    nl, nb, nc = X.shape
    FF = s.zeros((nb, nc))
    use_ff = s.ones((X.shape[0], X.shape[2])) > 0
    for b in range(nb):
        xsub = Xhoriz[:, b, :]
        xsubp = Xhorizp[:, b, :]
        mu = xsub.mean(axis=0)
        dists = abs(xsub - mu)
        distsp = abs(xsubp - mu)
        thresh = _percentile(dists.flatten(), 90.0)
        uthresh = dists * uniformity_thresh
        #use       = s.logical_and(dists<thresh, abs(dists-distsp) < uthresh)
        use = dists < thresh
        FF[b, :] = ((xsub*use).sum(axis=0)/use.sum(axis=0)) / \
            ((X[:, b, :]*use).sum(axis=0)/use.sum(axis=0))
        use_ff = s.logical_and(use_ff, use)
    return FF, Xhoriz, Xhorizp, s.array(use_ff) 

def star_graph_lumped(N, tau, gamma, I0, tmin, tmax, tcount):
    times = scipy.linspace(tmin, tmax, tcount)
    #    [[central node infected] + [central node susceptible]]
    #X = [Y_1^1, Y_1^2, ..., Y_1^{N}, Y_2^0, Y_2^1, ..., Y_2^{N-1}]
    X0 = scipy.zeros(2*N)  #length 2*N of just 0 entries
    #X0[I0]=I0*1./N #central infected, + I0-1 periph infected prob
    X0[N+I0] = 1#-I0*1./N #central suscept + I0 periph infected
    X = EoN.my_odeint(star_graph_dX, X0, times, args = (tau, gamma, N))
    #X looks like [[central susceptible,k periph] [ central inf, k-1 periph]] x T

    central_susc = X[:,N:]
    central_inf = X[:,:N]
    print(central_susc[-1][:])
    print(central_inf[-1][:])
    I = scipy.array([ sum(k*central_susc[t][k] for k in range(N))
              + sum((k+1)*central_inf[t][k] for k in range(N))
              for t in range(len(X))])
    S = N-I
    return times, S, I 

def compute_score(self, img_idx, y, x, mask):
    " Compute the score for deck or met with idx "
    qtrwin = self.winsize/2
    if mask==0:
      mask_file = self.datafiles[img_idx].split('.')[0] + '.jpg'
    elif mask==1:
      mask_file = self.datafiles[img_idx].split('.')[0] + '.msk.jpg'
    else:
      mask_file = self.datafiles[img_idx].split('.')[0] + '.shadow.jpg'
      
    selections_pad = N.zeros((self.height[img_idx] + self.winsize, 
                              self.width[img_idx] + self.winsize))
    mask_img  = cv2.imread(mask_file, 0)
    selections_pad[qtrwin:self.height[img_idx]+qtrwin,
                   qtrwin:self.width[img_idx]+qtrwin] = mask_img
    csel_mask = selections_pad[y:y+self.winsize, x:x+self.winsize]
    
    # Matches are pixels with intensity 255, so divide by this
    # to get number of matching pixels.
    return (csel_mask.sum()/255) 

def prepare_constants(ldpred_n,ns,m,p,h2,sampl_var_shrink_factor):
    Mp = m * p
    hdmp = (h2 / Mp)
    const_dict = {'Mp':Mp, 'hdmp':hdmp}
    rv_scalars = sp.zeros(m)
    if ldpred_n is not None:
        hdmpn = hdmp + 1.0 / ldpred_n
        hdmp_hdmpn = (hdmp / hdmpn)
        c_const = (p / sp.sqrt(hdmpn))
        d_const = (1.0 - p) / (sp.sqrt(1.0 / ldpred_n))
        const_dict['n']=ldpred_n
        const_dict['hdmpn']=hdmpn
        const_dict['hdmp_hdmpn']=hdmp_hdmpn
        const_dict['c_const']=c_const
        const_dict['d_const']=d_const
        rv_scalars[:]=sampl_var_shrink_factor* sp.sqrt((hdmp_hdmpn) * (1.0 / ldpred_n))
    else:
        snp_dict = {}
        for i in range(m):
            ni = ns[i]
            hdmpn_i = hdmp + 1.0 / ni
            hdmp_hdmpn_i = (hdmp / hdmpn_i)
            c_const_i = (p / sp.sqrt(hdmpn_i))
            d_const_i = (1.0 - p) / (sp.sqrt(1.0 / ni))
            snp_dict[i]={'n':ni, 'hdmpn':hdmpn_i, 
                           'hdmp_hdmpn':hdmp_hdmpn_i, 
                           'c_const':c_const_i, 
                           'd_const':d_const_i}
            rv_scalars[i]=sampl_var_shrink_factor* sp.sqrt((hdmp_hdmpn_i) * (1.0 / ni))
        const_dict['snp_dict']=snp_dict
    const_dict['rv_scalars']=rv_scalars
    return const_dict 

def FlipY(A, order):
    """Change sign of all coefficients with odd y power"""
    terms = (order+1)*(order+2)//2
    AF = scipy.zeros(terms)
    k = 0
    for i in range(order+1):
        for j in range (i+1):
            AF[k] = (-1)**(j)*A[k]
            k += 1
    return AF 

def smart_ld_pruning(scores, ld_table, max_ld=0.5, verbose=False, reverse=False):
    """
    Prunes SNPs in LD, but with smaller scores (p-values or betas)
    
    If using betas, set reversed to True.
    """
    if verbose:
        print('Performing smart LD pruning')
    t0 = time.time()
    if type(scores) == type([]):
        l = list(zip(scores, list(range(len(scores)))))
    else:
        l = list(zip(scores.tolist(), list(range(len(scores)))))
    l.sort(reverse=reverse)
    l = list(map(list, list(zip(*l))))
    rank_order = l[1]
    indices_to_keep = []
    remaining_indices = set(rank_order)
    for i in rank_order:
        if len(remaining_indices) == 0:
            break
        elif not (i in remaining_indices):
            continue
        else:
            indices_to_keep.append(i)
            for j in ld_table[i]:
                if ld_table[i][j] > max_ld and j in remaining_indices:
                    remaining_indices.remove(j)
    indices_to_keep.sort()
    pruning_vector = sp.zeros(len(scores), dtype='bool8')
    pruning_vector[indices_to_keep] = 1
    t1 = time.time()
    t = (t1 - t0)
    if verbose:
        print('\nIt took %d minutes and %0.2f seconds to LD-prune' % (t / 60, t % 60))
    return pruning_vector 

def get_snp_lrld_status(chromosome, positions, lrld_dict):
    snp_lrld = sp.zeros(len(positions))
    for snp_i in range(len(positions)):
        snp_lrld[snp_i] = is_in_lrld(chromosome, positions[snp_i], lrld_dict)
    return snp_lrld 

def _qqplot_bar(M=1000000, alphalevel = 0.05,distr = 'log10'):
    '''
    calculate error bars for a QQ-plot
    --------------------------------------------------------------------
    Input:
    -------------   ----------------------------------------------------
    M               number of points to compute error bars
    alphalevel      significance level for the error bars (default 0.05)
    distr           space in which the error bars are implemented
                    Note only log10 is implemented (default 'log10')
    --------------------------------------------------------------------
    Returns:
    -------------   ----------------------------------------------------
    betaUp          upper error bars
    betaDown        lower error bars
    theoreticalPvals    theoretical P-values under uniform
    --------------------------------------------------------------------
    '''


    #assumes 'log10'

    mRange=10**(sp.arange(sp.log10(0.5),sp.log10(M-0.5)+0.1,0.1));#should be exp or 10**?
    numPts=len(mRange);
    betaalphaLevel=sp.zeros(numPts);#down in the plot
    betaOneMinusalphaLevel=sp.zeros(numPts);#up in the plot
    betaInvHalf=sp.zeros(numPts);
    for n in xrange(numPts):
        m=mRange[n]; #numplessThanThresh=m;
        betaInvHalf[n]=st.beta.ppf(0.5,m,M-m);
        betaalphaLevel[n]=st.beta.ppf(alphalevel,m,M-m);
        betaOneMinusalphaLevel[n]=st.beta.ppf(1-alphalevel,m,M-m);
        pass
    betaDown=betaInvHalf-betaalphaLevel;
    betaUp=betaOneMinusalphaLevel-betaInvHalf;

    theoreticalPvals=mRange/M;
    return betaUp, betaDown, theoreticalPvals 

def FlipX(A, order):
    """Change sign of all coefficients with odd x power"""
    terms = (order+1)*(order+2)//2
    AF = scipy.zeros(terms)
    k = 0
    for i in range(order+1):
        for j in range (i+1):
            AF[k] = (-1)**(i-j)*A[k]
            k += 1
    return  AF 

def polyfit(u, x, y, order):
    '''Fit polynomial to a set of u values on an x, y grid
    u is a function u(x, y) being a polynomial of the form
    u = a[i, j] x**(i-j) y**j. x and y can be on a grid or be arbitrary values'''

    # First set up x and y powers for each coefficient
    px = []
    py = []
    for i in range(order+1):
        for j in range(i+1):
            px.append(i-j)
            py.append(j)
    terms = len(px)
    #print terms, ' terms for order ', order
    #print px
    #print py

    # Make up matrix and vector
    vector = scipy.zeros((terms))
    mat = scipy.zeros((terms, terms))
    for i in range(terms):
        vector[i] = (u*x**px[i]*y**py[i]).sum()
        for j in range(terms):
            mat[i, j] =  (x**px[i]*y**py[i]*x**px[j]*y**py[j]).sum()

    #print 'Vector', vector
    #print 'Matrix'
    #print mat
    imat = linalg.inv(mat)
    #print 'Inverse'
    #print imat
    # Check that inversion worked
    #print scipy.dot(mat, imat)
    coeffs = scipy.dot(imat, vector)
    return coeffs 

def ShiftCoeffs(a, xshift, yshift, order, verbose=False):
    '''Calculate coefficients of polynomial when shifted to new origin'''

    # First place in triangular layout
    at = scipy.zeros([order+1, order+1])
    atshift = scipy.zeros([order+1, order+1])
    ashift = scipy.zeros([len(a)]) # Copy shape of a
    k = 0
    for p in range(order+1):
        for q in range(p+1):
            at[p, q] = a[k]
            k+=1

    # Apply shift
    for p in range(order+1):
        for q in range(p+1):
            if verbose: print("A'%1d%1d" %(p, q))
            for i in range(p, order+1):
                for j in range(q, i+1-(p-q)):
                    f = choose(j, q)*choose(i-j, p-q)
                    atshift[p, q] = atshift[p, q] + f*xshift**((i-j)-(p-q))*yshift**(j-q)*at[i, j]
                    if verbose: print('%2d A(%1d, %1d) x^%1d y^%1d' %(f, i , j, i-j-(p-q), (j-q)))
            if verbose: print()

    # Put back in linear layout
    k = 0
    for p in range(order+1):
        for q in range(p+1):
            ashift[k] = atshift[p, q]
            k+=1

    return ashift 

def reorder(A, B, verbose=False) :
    '''Reorder Sabatke coefficients to my convention'''
    order = 5
    terms = (order+1)*(order+2)//2
    Aarray = scipy.zeros((order+1, order+1))
    Barray = scipy.zeros((order+1, order+1))
    #Carray = scipy.zeros((order+1, order+1))
    #Darray = scipy.zeros((order+1, order+1))
    k1 = 0
    for i in range(order+1):
        for j in range(order+1-i):
            Aarray[j, i] = A[k1]
            Barray[j, i] = B[k1]
            #Carray[j, i] = C[k1]
            #Darray[j, i] = D[k1]
            k1 += 1

    A2 = scipy.zeros((terms))
    B2 = scipy.zeros((terms))
    #C2 = scipy.zeros((terms))
    #D2 = scipy.zeros((terms))
    k2 = 0
    for i in range(order+1):
        for j in range(i+1):
            A2[k2] = Aarray[j, i-j]
            B2[k2] = Barray[j, i-j]
            #C2[k2] = Carray[j, i-j]
            #D2[k2] = Darray[j, i-j]
            k2 += 1

    if verbose:
        print('A')
        triangle(A2, order)
        print('\nB')
        triangle(B2, order)
        #print '\nC'
        #polyfit(C2, order)
        #print '\nD'
        #triangle(D2, order)

    return (A2, B2) #, C2, D2) 

def flatten(A, order):
    '''Convert triangular layout to linear array'''
    terms = (order+1)*(order+2)//2
    AF = scipy.zeros(terms)
    k = 0
    for i in range(order+1):
        for j in range(i+1):
            AF[k] = A[i, j]
            k += 1
    return AF 

def triangulate(A, order):
    '''Convert linear array to 2-D array with triangular coefficient layout'''
    AT = scipy.zeros((order+1, order+1))
    k = 0
    for i in range(order+1):
        for j in range(i+1):
            AT[i, j] = A[k]
            k += 1
    return AT 

def dLs_dsurface(self, x_surface, geom):
        """Partial derivative of surface emission with respect to state vector, 
        calculated at x_surface.  We append a column of zeros to handle 
        the extra glint parameter"""

        dLs_dsurface = super().dLs_dsurface(x_surface, geom)
        dLs_dglint = s.zeros((dLs_dsurface.shape[0],1))
        dLs_dsurface = s.hstack([dLs_dsurface, dLs_dglint]) 
        return dLs_dsurface