def createSimilarAR(data):
    """
    creates an AR-process that is similar to a given data set.
    data must be given in n x d-format
    """
    # step 1: get "average" fit matrix
    l_A = []
    for rep in arange(100):
        idx = randint(0,data.shape[0]-1,data.shape[0]-1)
        idat = data[idx,:]
        odat = data[idx+1,:]
        l_A.append(lstsq(idat,odat)[0])

    sysmat = meanMat(l_A).T
    # idea: get "dynamic noise" from input data as difference of
    # expected vs. predicted data:
    # eta_i = (sysmat*(data[:,i-1]).T - data[:,i])
    # however, in order to destroy any possible correlations in the
    # input noise (they would also occur in the output), the
    # noise per section has to be permuted.
    prediction = dot(sysmat,data[:-1,:].T)
    dynNoise = data[1:,:].T - prediction
    res = [zeros((dynNoise.shape[0],1)), ]
    for nidx in permutation(dynNoise.shape[1]):
        res.append( dot(sysmat,res[-1]) + dynNoise[:,nidx][:,newaxis] )
    
    return hstack(res).T

def reduceDim(fullmat,n=1):
    """
    reduces the dimension of a d x d - matrix to a (d-n)x(d-n) matrix, 
    keeping the largest eigenvalues unchanged.
    """
    u,s,v = svd(fullmat)
    return dot(u[:-n,:-n],dot(diag(s[:-n]),v[:-n,:-n]))

    def Strain_stress(self):
        '''
        Strain is obtained by Strain_Matrix (B)X Displacement (u). Stress is Stiffess_Matrix(E) x Strain.
        Von mises stress is also computed to study convergence.
        '''
        
        # copy displacement vector after reshaping with two columns for x and y         
        d = self.displacement.reshape(self.n_nodes,2) 
        
        #stress-strain calculation for each element
        for i in range(self.n_el): 
            
            el = self.element[i] # present element
            
            #Displacement formatted for an element
            disp=py.array([d[el[0]][0],d[el[0]][1],d[el[1]][0],d[el[1]][1],d[el[2]][0],d[el[2]][1]])
        
            #Element Strain vector = Product of Strain Matrix and Displacement
            [J,B] = self.B(el)
            strain = py.dot(B,disp.T)
            self.strain_vector[i] = strain
        
            #Element Stress vector = Product of Element K-Matrix and Strain Vector
            stress = py.dot(self.E_matrix,strain)
            self.stress_vector[i] = stress

            #von-mises stress for plotting
            self.von[i] = py.math.sqrt(0.5*((stress[0]-stress[1])**2 + stress[0]**2 + stress[1]**2 + 6*(stress[2])**2))

def getApices(y):
    """ 
    returns the time (in frames) and position of initial and final apex
    height from a given trajectory y, which are obtained by fitting a cubic
    spline 
    
    ==========
    parameter:
    ==========
    y : *array* (1D)
        the trajectory. Should ideally start ~1 frame before an apex and end ~1
        frame behind an apex

    ========
    returns:
    ========

    [x0, xF], [y0, yF] : location (in frames) and value of the first and final
    apices

    """
    # the math behind here is: fitting a 2nd order polynomial and finding 
    # the root of its derivative. Here, only the results are applied, that's
    # why it appears like "magic" numbers
    c = dot(array([[.5, -1, .5], [-1.5, 2., -.5], [1., 0., 0.]]), y[:3])
    x0 = -1. * c[1] / (2. * c[0])
    y0 = polyval(c, x0)

    c = dot(array([[.5, -1, .5], [-1.5, 2., -.5], [1., 0., 0.]]), y[-3:])
    xF = -1. * c[1] / (2. * c[0])
    yF = polyval(c, xF)
    xF += len(y) - 3

    return [x0, xF], [y0, yF]

def rotate_molecule(coords, rotp = m.array((0.,0.,0.)), phi = 0., \
        theta = 0., psi = 0.):
    """Rotate a molecule via Euler angles.

    See http://mathworld.wolfram.com/EulerAngles.html for definition.
    Input arguments:
    coords: Atom coordinates, as Nx3 2d pylab array.
    rotp: The point to rotate about, as a 1d 3-element pylab array
    phi: The 1st rotation angle around z axis.
    theta: Rotation around x axis.
    psi: 2nd rotation around z axis.

    """
# First move the molecule to the origin
# In contrast to MATLAB, numpy broadcasts the smaller array to the larger
# row-wise, so there is no need to play with the Kronecker product.
    rcoords = coords - rotp
# First Euler rotation about z in matrix form
    D = m.array(((m.cos(phi), m.sin(phi), 0.), (-m.sin(phi), m.cos(phi), 0.), \
            (0., 0., 1.)))
# Second Euler rotation about x:
    C = m.array(((1., 0., 0.), (0., m.cos(theta), m.sin(theta)), \
            (0., -m.sin(theta), m.cos(theta))))
# Third Euler rotation, 2nd rotation about z:
    B = m.array(((m.cos(psi), m.sin(psi), 0.), (-m.sin(psi), m.cos(psi), 0.), \
            (0., 0., 1.)))
# Total Euler rotation
    A = m.dot(B, m.dot(C, D))
# Do the rotation
    rcoords = m.dot(A, m.transpose(rcoords))
# Move back to the rotation point
    return m.transpose(rcoords) + rotp

def ssc(signal,samplerate=16000,winlen=0.025,winstep=0.01,
          nfilt=26,nfft=512,lowfreq=0,highfreq=None,preemph=0.97):
    """Compute Spectral Subband Centroid features from an audio signal.

    :param signal: the audio signal from which to compute features. Should be an N*1 array
    :param samplerate: the samplerate of the signal we are working with.
    :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds)    
    :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds)    
    :param nfilt: the number of filters in the filterbank, default 26.
    :param nfft: the FFT size. Default is 512.
    :param lowfreq: lowest band edge of mel filters. In Hz, default is 0.
    :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2
    :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97. 
    :returns: A numpy array of size (NUMFRAMES by nfilt) containing features. Each row holds 1 feature vector. 
    """          
    highfreq= highfreq or samplerate/2
    signal = sigproc.preemphasis(signal,preemph)
    frames = sigproc.framesig(signal, winlen*samplerate, winstep*samplerate)
    pspec = sigproc.powspec(frames,nfft)
    pspec = pylab.where(pspec == 0,pylab.finfo(float).eps,pspec) # if things are all zeros we get problems
    
    fb = get_filterbanks(nfilt,nfft,samplerate,lowfreq,highfreq)
    feat = pylab.dot(pspec,fb.T) # compute the filterbank energies
    R = pylab.tile(pylab.linspace(1,samplerate/2,pylab.size(pspec,1)),(pylab.size(pspec,0),1))
    
    return pylab.dot(pspec*R,fb.T) / feat

def varRed(idat,odat,A,bootstrap = None):
    """
    computed the variance reduction when using A*idat[:,x].T as predictor for odat[:,x].T
    if bootstrap is an integer > 1, a bootstrap with the given number of iterations
    will be performed.
    returns
    tVred, sVred: the total relative variance after prediction (all coordinates)
       and the variance reduction for each coordinate separately. These data are
       scalar and array or lists of scalars and arrays when a bootstrap is performed.
       
    Note: in the bootstrapped results, the first element refers to the "full" 
    data variance reduction.
    """
    
    nBoot = bootstrap if type(bootstrap) is int else 0
    if nBoot < 2:
        nBoot = 0        
    
    odat_pred = dot(A,idat.T)
    rdiff = odat_pred - odat.T # remaining difference
    rvar = var(rdiff,axis=1)/var(odat.T,axis=1) # relative variance
    trvar = var(rdiff.flat)/var(odat.T.flat)    # total relative variance
    
    if nBoot > 0:
        rvar = [rvar,]
        trvar = [trvar,]
    for rep in range(nBoot-1):
        indices = randint(0,odat.T.shape[1],odat.T.shape[1])
        odat_pred = dot(A,idat[indices,:].T)
        rdiff = odat_pred - odat[indices,:].T # remaining difference
        rvar.append( var(rdiff,axis=1)/var(odat.T,axis=1) ) # relative variance
        trvar.append (var(rdiff.flat)/var(odat.T.flat) )    # total relative variance
        
    return trvar, rvar

def fBM_nd(dims, H, return_mat = False, use_eig_ev = True):
    """
    creates fractional Brownian motion
    parameters: dims is a tuple of the shape of the sample path (nxd); 
                H: Hurst exponent
    this is the slow version of fBM. It might, however, be more precise than
    fBM, however - sometimes, the matrix square root has a problem, which might
    induce inaccuracy    
    use_eig_ev: use eigenvalue decomposition for matrix square root computation
    (faster)
    """
    n = dims[0]
    d = dims[1]
    Gamma = zeros((n,n))
    print ('building ...\n')
    for t in arange(n):
        for s in arange(n):
            Gamma[t,s] = .5*((s+1)**(2.*H) + (t+1)**(2.*H) - abs(t-s)**(2.*H))
    print('rooting ...\n')    
    if use_eig_ev:
        ev,ew = eig(Gamma.real)
        Sigma = dot(ew, dot(diag(sqrt(ev)),ew.T) )
    else:
        Sigma = sqrtm(Gamma)
    if return_mat:
        return Sigma
    v = randn(n,d)
    return dot(Sigma,v)

	def Q_calc(self,X):

		"""
			calculates Q (n_x by n_theta) matrix of the IDE model at  each time step
	
			Arguments
			----------
			X: list of ndarray
				state vectors

			Returns
			---------
			Q : list of ndarray (n_x by n_theta)
		"""

		Q=[]	
		T=len(X)
		Psi=self.model.Gamma_inv_psi_conv_Phi
		Psi_T=pb.transpose(self.model.Gamma_inv_psi_conv_Phi,(0,2,1))

		for t in range(T):

			firing_rate_temp=pb.dot(X[t].T,self.model.Phi_values)
			firing_rate=self.model.act_fun.fmax/(1.+pb.exp(self.model.act_fun.varsigma*(self.model.act_fun.v0-firing_rate_temp)))	

			#calculate q
			g=pb.dot(firing_rate,Psi_T)

			g *=(self.model.spacestep**2)	
			q=self.model.Ts*g
			q=q.reshape(self.model.nx,self.model.n_theta)
			Q.append(q)
		return Q

    def estimate_kernel(self, X, P, M):

        """estimate the ide model's kernel weights from data stored in the ide object"""

        # form Xi variables
        Xi_0 = pb.zeros([self.model.nx, self.model.nx])
        Xi_1 = pb.zeros([self.model.nx, self.model.nx])
        for t in range(1, len(X)):
            Xi_0 += pb.dot(X[t - 1, :].reshape(self.model.nx, 1), X[t, :].reshape(self.model.nx, 1).T) + M[
                t, :
            ].reshape(self.model.nx, self.model.nx)
            Xi_1 += pb.dot(X[t - 1, :].reshape(self.model.nx, 1), X[t - 1, :].reshape(self.model.nx, 1).T) + P[
                t - 1, :
            ].reshape(self.model.nx, self.model.nx)

            # form Upsilon and upsilons
        Upsilon = pb.zeros([self.model.ntheta, self.model.ntheta])
        upsilon0 = pb.zeros([1, self.model.ntheta])
        upsilon1 = pb.zeros([1, self.model.ntheta])
        for i in range(self.model.nx):
            for j in range(self.model.nx):
                Upsilon += Xi_1[i, j] * self.model.Delta_Upsilon[j, i]
                upsilon0 += Xi_0[i, j] * self.model.Delta_upsilon[j, i]
                upsilon1 += Xi_1[i, j] * self.model.Delta_upsilon[j, i]
        upsilon1 = upsilon1 * self.model.xi
        Upsilon = Upsilon * self.model.Ts * self.model.varsigma

        weights = pb.dot(pb.inv(Upsilon.T), upsilon0.T - upsilon1.T)

        return weights

def main():
    mu = pl.array([[0], [12], [24], [36]])
    Sigma = pl.array([[3.01602775,  1.02746769, -3.60224613, -2.08792829],
                      [1.02746769,  5.65146472, -3.98616664,  0.48723704],
                      [-3.60224613, -3.98616664, 13.04508284, -1.59255406],
                      [-2.08792829,  0.48723704, -1.59255406,  8.28742469]])

    # The data matrix is created for above mu and Sigma.
    d, U = pl.eig(Sigma)
    L = pl.diagflat(d)
    A = pl.dot(U, pl.sqrt(L))
    X = pl.randn(4, 1000)

    # Y is the data matrix of random samples.
    Y = pl.dot(A, X) + pl.tile(mu, 1000)

    pl.figure(1)
    pl.clf()
    pl.plot(X[0], Y[1], '+', color='#0000FF', label='i=0,j=1')
    pl.plot(X[0], Y[2], '+', color='#FF0000', label='i=0,j=2')
    pl.plot(X[0], Y[3], '+', color='#00FF00', label='i=0,j=3')
    pl.plot(X[1], Y[0], 'x', color='#FFFF00', label='i=1,j=0')
    pl.plot(X[1], Y[2], 'x', color='#00FFFF', label='i=1,j=2')
    pl.plot(X[1], Y[3], 'x', color='#444444', label='i=1,j=3')
    pl.plot(X[2], Y[0], '.', color='#774411', label='i=2,j=0')
    pl.plot(X[2], Y[1], '.', color='#222222', label='i=2,j=1')
    pl.plot(X[2], Y[3], '.', color='#AAAAAA', label='i=2,j=3')
    pl.plot(X[3], Y[0], '+', color='#FFAA22', label='i=3,j=0')
    pl.plot(X[3], Y[1], '+', color='#22AAFF', label='i=3,j=1')
    pl.plot(X[3], Y[2], '+', color='#FFDD00', label='i=3,j=2')
    pl.legend()
    pl.savefig('fig21.png')

    def Global_Stiffness(self):
        '''
        Generates Global Stiffness Matrix for the plane structure
        '''
        elem = self.element;
        B = py.zeros((6,6))
        for i in range (0,py.size(elem,0)): 
            #for each element find the stifness matrix
            K = py.zeros((self.n_nodes*2,self.n_nodes*2))            
            el = elem[i]
            
            #nodes formatted for input            
            [node1, node2, node3] = el;
            node1x = 2*(node1-1);node2x = 2*(node2-1);node3x = 2*(node3-1);
            node1y = 2*(node1-1)+1;node2y = 2*(node2-1)+1;node3y = 2*(node3-1)+1;
            
            #Area, Strain Matrix and E Matrix multiplied to get element stiffness            
            [J,B] = self.B(el)
            local_k =0.5*abs(J)*py.dot(py.transpose(B),py.dot(self.E_matrix,B))
            if self.debug:            
                print 'K for elem', el, '\n', local_k
                
            #Element K-Matrix converted into Global K-Matrix format 
            K[py.ix_([node1x,node1y,node2x,node2y,node3x,node3y],[node1x,node1y,node2x,node2y,node3x,node3y])] = K[py.ix_([node1x,node1y,node2x,node2y,node3x,node3y],[node1x,node1y,node2x,node2y,node3x,node3y])]+local_k

            #Adding contibution into Global Stiffness           
            self.k_global = self.k_global + K
            
        if self.debug:            
                print 'Global Stiffness','\n', self.k_global        

def renormalize(x_unpurt,x_before,x_purt,epsilon,N):
    # BEFORE ANYTHING: make sure particles near boundaries are shuffeled into places where where the
    # seam is not between any purturbed and fudicial trajectories.
    x_unpurt,x_purt = shuff(x_unpurt,x_purt,N)

    # The trajectory we are going to be returning is going to be the new one for the next run. lets
    # call it
    x_new = pl.copy(x_unpurt)
    # copied it because we are going to add the small amounts to it to purturb it.

    # lets find a vector pointing in the direction of the trajectories path. For this we need the
    # fiducual point at t-dt, which is given to us in the function as x_before. find the vector
    # between x_before and x_unpurt
    traj_vec = x_unpurt-x_before 
    # normalize it
    traj_vec = traj_vec/pl.sqrt(pl.dot(traj_vec,traj_vec))
    print('traj_vec magnitude (should be 1): ' + str(pl.sqrt(pl.dot(traj_vec,traj_vec))))

    # Now lets see how close the vector pointing from the fidicial to the perturbed trajectorie is
    # to orthogonal with the trajectory... should get closer to 1 as we check more because it should
    # be aligning itself with the axis of greatest expansion and that should be orthogonal.
    # First normalize the difference vector
    diff_vec = x_unpurt - x_purt
    # normalize it
    diff_vec = diff_vec/pl.sqrt(pl.dot(diff_vec,diff_vec))
    print('diff_vec magnitude (should be 1): ' + str(pl.sqrt(pl.dot(diff_vec,diff_vec))))
    print('normalized(x_unpurt-x_purt)dot(traj_vec)  (should get close to 0): '+ str(pl.dot(diff_vec,traj_vec)))

    # for now lets just return a point moved back along the difference vector. no gram shmidt or
    # anything.
    return x_new + epsilon*diff_vec

def plotEnsemble2D(ens,v1,v2,colordata=None,hess=None,\
		   size=50,labelBest=True,ensembleAlpha=0.75,contourAlpha=1.0):
	"""
	Plots a 2-dimensional projection of a given parameter
	ensemble, along given directions:
	     -- If v1 and v2 are scalars, project onto plane given by
		those two bare parameter directions.
	     -- If v1 and v2 are vectors, project onto those two vectors.
	
	When given colordata (either a single color, or an array
	of different colors the length of ensemble size), each point
	will be assigned a color based on the colordata.
	
	With labelBest set, the first point in the ensemble is
	plotted larger (to show the 'best fit' point for a usual 
	parameter ensemble).
	
	If a Hessian is given, cost contours will be plotted
	using plotContours2D.
	"""
	if pylab.shape(v1) is ():
		xdata = pylab.transpose(ens)[v1]
		ydata = pylab.transpose(ens)[v2]
		
		# label axes
		param1name, param2name = '',''
		try:
		    paramLabels = ens[0].keys()
		except:
		    paramLabels = None
		if paramLabels is not None:
		    param1name = ' ('+paramLabels[param1]+')'
		    param2name = ' ('+paramLabels[param2]+')'
		pylab.xlabel('Parameter '+str(v1)+param1name)
		pylab.ylabel('Parameter '+str(v2)+param2name)
	else:
		xdata = pylab.dot(ens,v1)
		ydata = pylab.dot(ens,v2)

	if colordata==None:
		colordata = pylab.ones(len(xdata))
		
	if labelBest: # plot first as larger circle
		if pylab.shape(colordata) is (): # single color
		    colordata0 = colordata
		    colordataRest = colordata
		else: # specified colors
		    colordata0 = [colordata[0]]
		    colordataRest = colordata[1:]
		scatterColors(xdata[1:],ydata[1:],colordataRest,		\
				size,alpha=ensembleAlpha)
		scatterColors([xdata[0]],[ydata[0]],colordata0,			\
				size*4,alpha=ensembleAlpha)
	else:
		scatterColors(xdata,ydata,colordata,size,alpha=ensembleAlpha)
		
	if hess is not None:
		plotApproxContours2D(hess,param1,param2,pylab.array(ens[0]),	\
			alpha=contourAlpha)

def brute_force_Z_vis_gauss(W):
    v,h = W.v, W.h
    Z = zeros(2**h)
    for i in xrange(2**h):
        H = int_to_bin(i, h)
        b = W.T() * H
        Z[i] = .5 * dot(b,b.T) + dot(H,W[2])
    return log_sum_exp(Z)

def simCSLIP_xp(x0, x0R, x0L, p0R, p0L, AR, AL, SLIP_param0, n=50):
    """
    simulates the controlled 2step-SLIP, using [x,p]-referenced control
    input:
        x0 - initial (augmented) state, e.g. [x0L, p0R].T
        x0R - reference right apex (y, vx, vz)
        x0L - reference left apex     -"-
        p0R - reference right parameters
        p0L - reference left parameters
        AR - parameter control right leg
        AL - parameter control left leg
        SLIP_param0: dict, containing {'m': ..., 'g': ... }
        n - number of strides to simulate at most
    """
    res = []
    refStateL = hstack([x0L, squeeze(sp_d2a(p0R))])[:,newaxis]
    refStateR = hstack([x0R, squeeze(sp_d2a(p0L))])[:,newaxis]
    currState = array(x0)
    slip_params = copy.deepcopy(SLIP_param0)
    if currState.ndim == 1:
        currState = currState[:,newaxis]
    elif currState.shape[0] == 1:
        currState = currState.T
    for step in range(n):
        #print 'AL: ', AL.shape, 'p0L: ', sp_d2a(p0L).shape
        pL = sp_d2a(p0L) + dot(AL, currState - refStateL)
        #print 'pL changed:', not allclose(pL,sp_d2a(p0L))
        slip_params.update(sp_a2d(pL))
        try:
            resL = sl.SLIP_step3D(currState[:3,0], slip_params)
        except ValueError:
            print 'simulation aborted (l1)\n'
            break
        if resL['sim_fail']:
            print 'simulation aborted (l2)\n'
            break
        res.append(resL)
        currState = hstack([resL['y'][-1],
                            resL['vx'][-1],
                            resL['vz'][-1],
                            squeeze(pL)])[:,newaxis]
        pR = sp_d2a(p0R) + dot(AR, currState - refStateR)
        #print 'pR changed:', not allclose(pR,sp_d2a(p0R))
        slip_params.update(sp_a2d(pR))
        try:
            resR = sl.SLIP_step3D(currState[:3,0], slip_params)
        except ValueError:
            print 'simulation aborted (r1)\n'
            break
        if resR['sim_fail']:
            print 'simulation aborted (r2)\n'
            break
        res.append(resR)
        currState = hstack([resR['y'][-1],
                            resR['vx'][-1],
                            resR['vz'][-1],
                            squeeze(pR)])[:,newaxis]
    return res

	def  state_equation(self,x):

		'''state equation for sigma points propogation '''
		firing_rate_temp=pb.dot(x.T,self.Phi_values)
 		firing_rate=self.act_fun.fmax/(1.+pb.exp(self.act_fun.varsigma*(self.act_fun.v0-firing_rate_temp))).T	
		g=pb.dot(self.Gamma_inv_Psi_conv_Phi,firing_rate)
		g *=(self.spacestep**2)
		x=self.Ts*g+self.xi*x
		return x

    def _get_eff(window_positions,window_steps,window_length):
        s2 = 0
        dispVec = window_steps[window_length-1,:] - window_steps[0,:]
        disp = pl.dot(dispVec,dispVec)
        for i in range(window_length-1):
            s2 += pl.dot(window_steps[i,:],window_steps[i,:])

        efficiency = disp/(window_length*s2+epsilon)
    
        return efficiency

def JointValuesFromTransform(robot, Tdesired):
    """
    Find the dummy joint values such that the Transform of the baselink is
    Tdesired. T0 is the initial Transform of the baselink.
    """
    R = dot(inv(robot.baselinkinittransform[0:3, 0:3]), Tdesired[0:3, 0:3])
    [h1, h2, h3] = AnglesFromRot(R)
    offset = dot(R, robot.baselinkinittransform[0:3, 3])
    [s1, s2, s3] = Tdesired[0:3, 3] - offset
    return [s1, s2, s3, h1, h2, h3]

def find_factors(idat, odat, k = None):
    """
    A routine to compute the main predictors (linear combinations of
    coordinates) in idat to predict odat.

    *Parameters*
        idat: d x n data matrix,
            with n measurements, each of dimension d
        odat: q x n data matrix
             with n measurements, each of dimension q

    *Returns*
        **Depending on whether or not** *k* **is provided, the returned
        value is different**

      * if k is given, compute the first k regressors and return an orthogonal
        matrix that contains the regressors in its colums,
        i.e. reg[0,:] is the first regressor

      * if k is not given or None, return a d-dimensional vector v(k)
        explaining which fraction of the total predictable variance can be
        explained using only k regressors.

    **NOTE**


    #. idat and odat must have zero mean
    #. To interpret the regressors, it is advisable to have the
       for columns of idat having the same variance
    """
    # transform into z-scores
    u, s, v = svd(idat, full_matrices = False)
    su = dot(diag(1./s), u.T)

    z = dot(su,idat)
    # ! Note that the covariance of z is *NOT* 1, but 1/n; z*z.T = 1 !

    # least-squares regression:
    A = dot(odat, pinv(z))

    uA, sigma_A, vA = svd(A, full_matrices = False)
    if k is None:
        vk = cumsum(sigma_A**2) / sum(sigma_A**2)
        return vk

    else:
    # choose k predictors
        sigma_A1 = sigma_A.copy()
        sigma_A1[k:] = 0
        A1 = reduce(dot, [uA, diag(sigma_A1), vA])
        B = dot(A1, su)
        uB, sigma_B, vB = svd(B, full_matrices = False)
        regs = vB[:k,:].T
        return regs