def cca(x_tn,y_tm, reg=0.00000001):
    x_tn = x_tn-x_tn.mean(axis=0)
    y_tm = y_tm-y_tm.mean(axis=0)
    N = x_tn.shape[1]
    M = y_tm.shape[1]
    xy_tq = c_[x_tn,y_tm]
    cqq = cov(xy_tq,rowvar=0)
    cxx = cqq[:N,:N]+reg*np.eye(N)+0.000000001*np.ones((N,N))
    cxy = cqq[:N,N:(N+M)]+0.000000001*np.ones((N,N))
    cyx = cqq[N:(N+M),:N]+0.000000001*np.ones((N,N))
    cyy = cqq[N:(N+M),N:(N+M)]+reg*np.eye(N)+0.000000001*np.ones((N,N))
    
    K = min(N,M)
    
    xldivy = lstsq(cxx,cxy)[0]
    yldivx = lstsq(cyy,cyx)[0]
    #print xldivy
    #print dot(np.linalg.inv(cxx),cxy)
    _,vecs = eig(dot(xldivy,yldivx))
    a_nk = vecs[:,:K]
    #print normr(vecs.T)
    b_mk = dot(yldivx,a_nk)
    
    u_tk = dot(x_tn,a_nk)
    v_tk = dot(y_tm,b_mk)
    
    return a_nk,b_mk,u_tk,v_tk

    def initialize(self):
        S = sum([np.dot(unit.X.T, unit.X) for unit in self.units])
        Y = sum([np.dot(unit.X.T, unit.Y) for unit in self.units])
        self.a = L.lstsq(S, Y)[0]

        D = 0
        t = 0
        sigmasq = 0
        for unit in self.units:
            unit.r = unit.Y - np.dot(unit.X, self.a)
            if self.q > 1:
                unit.b = L.lstsq(unit.Z, unit.r)[0]
            else:
                Z = unit.Z.reshape((unit.Z.shape[0], 1))
                unit.b = L.lstsq(Z, unit.r)[0]

            sigmasq += (np.power(unit.Y, 2).sum() -
                        (self.a * np.dot(unit.X.T, unit.Y)).sum() -
                        (unit.b * np.dot(unit.Z.T, unit.r)).sum())
            D += np.multiply.outer(unit.b, unit.b)
            t += L.pinv(np.dot(unit.Z.T, unit.Z))

        #TODO: JP added df_resid check
        self.df_resid = (self.N - (self.m - 1) * self.q - self.p)
        sigmasq /= (self.N - (self.m - 1) * self.q - self.p)
        self.sigma = np.sqrt(sigmasq)
        self.D = (D - sigmasq * t) / self.m

 def _calculate_log_likelihood(self):
     #if self.m == None:
     #    Give error message
     R = zeros((self.n, self.n))
     X,Y = array(self.X), array(self.Y)
     thetas = 10.**self.thetas
     for i in range(self.n):
         for j in arange(i+1,self.n):
             R[i,j] = (1-self.nugget)*e**(-sum(thetas*(X[i]-X[j])**2.)) #weighted distance formula
     R = R + R.T + eye(self.n)
     self.R = R
     one = ones(self.n)
     try:
         self.R_fact = cho_factor(R)
         rhs = vstack([Y, one]).T
         R_fact = (self.R_fact[0].T,not self.R_fact[1])
         cho = cho_solve(R_fact, rhs).T
         
         self.mu = dot(one,cho[0])/dot(one,cho[1])
         self.sig2 = dot(Y-dot(one,self.mu),cho_solve(self.R_fact,(Y-dot(one,self.mu))))/self.n
         #self.log_likelihood = -self.n/2.*log(self.sig2)-1./2.*log(abs(det(self.R)+1.e-16))-sum(thetas)
         self.log_likelihood = -self.n/2.*log(self.sig2)-1./2.*log(abs(det(self.R)+1.e-16))
     except (linalg.LinAlgError,ValueError):
         #------LSTSQ---------
         self.R_fact = None #reset this to none, so we know not to use cholesky
         #self.R = self.R+diag([10e-6]*self.n) #improve conditioning[Booker et al., 1999]
         rhs = vstack([Y, one]).T
         lsq = lstsq(self.R.T,rhs)[0].T
         self.mu = dot(one,lsq[0])/dot(one,lsq[1])
         self.sig2 = dot(Y-dot(one,self.mu),lstsq(self.R,Y-dot(one,self.mu))[0])/self.n
         self.log_likelihood = -self.n/2.*log(self.sig2)-1./2.*log(abs(det(self.R)+1.e-16))

def dfa(x, ave=None, l=None):
    x = np.array(x)

    if ave is None:
        ave = np.mean(x)

    y = np.cumsum(x)
    y -= ave

    if l is None:
        l = np.floor(len(x) * 1 / (2 ** np.array(range(4, int(np.log2(len(x))) - 4))))

    f = np.zeros(len(l))  # f(n) of different given box length n

    for i in xrange(0, len(l)):
        n = int(l[i])  # for each box length L[i]
        if n == 0:
            print "time series is too short while the box length is too big"
            print "abort"
            exit()
        for j in xrange(0, len(x), n):  # for each box
            if j + n < len(x):
                c = range(j, j + n)
                c = np.vstack([c, np.ones(n)]).T  # coordinates of time in the box
                y = y[j:j + n]  # the value of data in the box
                f[i] += lstsq(c, y)[1]  # add residue in this box
        f[i] /= ((len(x) / n) * n)
    f = np.sqrt(f)

    alpha = lstsq(np.vstack([np.log(l), np.ones(len(l))]).T, np.log(f))[0][0]

    return alpha

def simplex(A,b,c,basis):
    B=A[:,basis]
    xx,resid,rank,s = lin.lstsq(B,b)
    x = np.zeros(c.shape[0])    
    v = np.zeros(c.shape[0])
    x[basis[0]] = xx[0,0]
    cost = c[basis].T*x[basis]     # cost at starting corner
    for iteration in np.arange(100):
        y = lin.lstsq(B.T,c[basis])           # this y may not be feasible
        y = np.array([y[0]]).T
        idx = (c - np.dot(A.T,y).T).argmin()
        rmin = (c - np.dot(A.T,y).T).min()
        if rmin >= -0.00000001:      # optimality is reached, r>=0
            break                 # current x and y are optimal
         
        print B
        print A 
        print A[:,idx]
        print idx
        v[basis] = lin.lstsq(B,A[:,idx])[0]
        tmp = x[basis] / np.max(v[basis],.000001)
        out = tmp.argmin()
        minratio = tmp.min()
        if v[out] == .000001:  # out = index of first x to reach 0
            break      # break when that edge is extremely short
            
        cost = cost + minratio*rmin  # lower cost at end of step
        x[basis] = x[basis] - minratio*v[basis]   # update old x
        x[idx] = minratio      # find new positive component of x
        basis[out] = idx      # replace old index by new in basis
        print basis
    
    return x,y,cost

    def initialize(self):
        S = sum([N.dot(unit.X.T, unit.X) for unit in self.units])
        Y = sum([N.dot(unit.X.T, unit.Y) for unit in self.units])
        self.a = L.lstsq(S, Y)[0]

        D = 0
        t = 0
        sigmasq = 0
        for unit in self.units:
            unit.r = unit.Y - N.dot(unit.X, self.a)
            if self.q > 1:
                unit.b = L.lstsq(unit.Z, unit.r)[0]
            else:
                Z = unit.Z.reshape((unit.Z.shape[0], 1))
                unit.b = L.lstsq(Z, unit.r)[0]

            sigmasq += (N.power(unit.Y, 2).sum() -
                        (self.a * N.dot(unit.X.T, unit.Y)).sum() -
                        (unit.b * N.dot(unit.Z.T, unit.r)).sum())
            D += N.multiply.outer(unit.b, unit.b)
            t += L.pinv(N.dot(unit.Z.T, unit.Z))

        sigmasq /= (self.N - (self.m - 1) * self.q - self.p)
        self.sigma = N.sqrt(sigmasq)
        self.D = (D - sigmasq * t) / self.m

def con2vert(A, b):
    """
    Convert sets of constraints to a list of vertices (of the feasible region).
    If the shape is open, con2vert returns False for the closed property.
    """
    # Python implementation of con2vert.m by Michael Kleder (July 2005),
    #  available: http://www.mathworks.com/matlabcentral/fileexchange/7894
    #  -con2vert-constraints-to-vertices
    # Author: Michael Kelder (Original)
    #         Andre Campher (Python implementation)
    c = linalg.lstsq(mat(A), mat(b))[0]
    btmp = mat(b)-mat(A)*c
    D = mat(A)/matlib.repmat(btmp, 1, A.shape[1])

    fmatv = qhull(D, "Ft") #vertices on facets

    G  = zeros((fmatv.shape[0], D.shape[1]))
    for ix in range(0, fmatv.shape[0]):
        F = D[fmatv[ix, :], :].squeeze()
        G[ix, :] = linalg.lstsq(F, ones((F.shape[0], 1)))[0].transpose()

    V = G + matlib.repmat(c.transpose(), G.shape[0], 1)
    ux = uniqm(V)

    eps = 1e-13
    Av = dot(A, ux.T)
    bv = tile(b, (1, ux.shape[0]))
    closed = sciall(Av - bv <= eps)

    return ux, closed

def solve_and_check_impl(A,Ah,B,Bh):
    Am=numpy.matrix(A)
    Amh=numpy.matrix(Ah)
    Bm=numpy.matrix(B)
    Bmh=numpy.matrix(Bh)
    Ar,Aresidue=linalg.lstsq(Am,Amh)[0:2]
    Br,Bresidue=linalg.lstsq(Bm,Bmh)[0:2]
    A11,A12,Tx=Ar.transpose()[0].tolist()[0]
    A21,A22,Ty=Br.transpose()[0].tolist()[0]
    mA=numpy.matrix([[A11,A12],[A21,A22]])
    T=numpy.matrix([[Tx],[Ty]])
    print "Determinant:",linalg.det(mA)
    error=0
    for m,mh in zip(A,Ah):
        #print "m:",m,"mh:",mh        
        x,y,one=m
        deg=mh[0]
        X=numpy.matrix([[x],[y]])        
        Y=mA*X+T
        lat,lon=Y[0,0],Y[1,0]
        #print "Mapped lat",lat,"correct",deg
        error=max(error,(deg-lat))
    for m,mh in zip(B,Bh):
        #print "m:",m,"mh:",mh        
        x,y,one=m
        deg=mh[0]
        X=numpy.matrix([[x],[y]])        
        Y=mA*X+T
        lat,lon=Y[0,0],Y[1,0]
        #print "Mapped lon",lon,"correct",deg
        error=max(error,(deg-lon))
    Ai=linalg.inv(mA)
        
    return error,mA,T

def multiLinearExampleWithLeastSQR():
  from numpy.linalg import lstsq
  X = [[1,6,2],[1,8,1],[1,10,0],[1,14,2],[1,18,0]]
  y =[[7],[9],[13],[17.5],[18]]  
  print "\n Least sqr:"
  print lstsq(X,y)[0]
  print lstsq(X,y)[0]  

def calibrate_peak_ensemble(models_forecasts, measurements, forecast_len = 48, peak_level = 80):
    T_predictors = []
    T_target = []
    H_predictors = []
    H_target = []

    rng = min(len(model) for model in models_forecasts)
    for tm in range(rng):
        msm = measurements[tm * 6: tm * 6 + forecast_len]
        measured_peaks = detect_peaks(msm, peak_level)
        if not measured_peaks:
            continue

        forecasts = [prd[tm] for prd in models_forecasts]
        forecasts_peaks = [detect_peaks(fcst, peak_level) for fcst in forecasts]
        
        forecasts_peaks_cor = [list(map(lambda x: find_corresponding_peak(x, forecast_peaks),
                                      measured_peaks)) for forecast_peaks in forecasts_peaks]
        for measured, *corresponding in zip(measured_peaks, *forecasts_peaks_cor):
            if all(corresponding):
                H_predictors.append([peak[3] for peak in corresponding] + [1])
                T_predictors.append([peak[2] for peak in corresponding] + [1])
                H_target.append(measured[3])
                T_target.append(measured[2])

    print(H_predictors, H_target)

    H_coefs = lstsq(H_predictors, H_target)[0]
    T_coefs = lstsq(T_predictors, T_target)[0]
    return list(T_coefs), list(H_coefs)

def getOrthColumns(m):
    '''
    Constructs the orthogonally complementing columns of the input.

    Input of the form pxr is assumed to have r<=p,
    and have either full column rank r or rank 0 (scalar or matrix)
    Output is of the form px(p-r), except:
    a) if M square and full rank p, returns scalar 0
    b) if rank(M)=0 (zero matrix), returns I_p
    (Note you cannot pass scalar zero, because dimension info would be
    missing.)
    Return type is as input type.
    '''
    if type(m) == type(asarray(m)):
        m = mat(m)
        output = 'array'
    else: output = 'matrix'
    p, r = m.shape
    # first catch the stupid input case
    if p < r: raise ValueError, 'need at least as many rows as columns'
    # we use lstsq(M, ones) just to exploit its rank-finding algorithm,
    rk = lstsq(m, ones(p).T)[2]
    # first the square and full rank case:
    if rk == p: result = zeros((p,0))   # note the shape! hopefully octave-like
    # then the zero-matrix case (within machine precision):
    elif rk == 0: result = eye(p)
    # now the rank-deficient case:
    elif rk < r:
        raise ValueError, 'sorry, matrix does not have full column rank'
    # (what's left should be ok)
    else:
        # we have to watch out for zero rows in M,
        # if they are in the first p-r positions!
        # so the (probably inefficient) algorithm:
            # 1. check the rank of each row
            # 2. if zero, then also put a zero row in c
            # 3. if not, put the next unit vector in c-row
        idr = eye(r)
        idpr = eye(p-r)
        c = empty([0,r])    # starting point  
        co = empty([0, p-r]) # will hold orth-compl.
        idrcount = 0
        for row in range(p):
            # (must be ones() instead of 1 because of 2d-requirement
            if lstsq( m[row,:], ones(1) )[2] == 0 or idrcount >= r:
                c = r_[ c, zeros(r) ]
                co = r_[ co, idpr[row-idrcount, :] ]
            else:     # row is non-zero, and we haven't used all unit vecs 
                c = r_[ c, idr[idrcount, :] ] 
                co = r_[ co, zeros(p-r) ]
                idrcount += 1
        # earlier non-general (=bug) line: c = mat(r_[eye(r), zeros((p-r, r))])
        # and:  co = mat( r_[zeros((r, p-r)), eye(p-r)] )
        # old:
        # result = ( eye(p) - c * (M.T * c).I * M.T ) * co
        result = co - c * solve(m.T * c, m.T * co)
    if output == 'array': return result.A
    else: return result

def my_nmf(document_term_mat, n_components=15, n_iterations=50, eps=1e-6):
    n_rows, n_cols = document_term_mat.shape
    W = rand(n_rows*n_components).reshape([n_rows, n_components])
    H = rand(n_components*n_cols).reshape([n_components, n_cols])
    # linalg.lstsq doesn't work on sparse mats
    dense_document_term_mat = document_term_mat.todense()
    for i in range(n_iterations):
        H = linalg.lstsq(W, dense_document_term_mat)[0].clip(eps)
        W = linalg.lstsq(H.T, dense_document_term_mat.T)[0].clip(eps).T
    return array(W), array(H)

def glrd_diverse(V, G, F, r, err_V, err_F):
    # diversity threshold is 0.5
    for k in xrange(r):
        G_copy = np.copy(G)  # create local copies for excluding the k^th col and row of G and F resp.
        F_copy = np.copy(F)
        G_copy[:, k] = 0.0
        F_copy[k, :] = 0.0

        R = V - np.dot(G_copy, F_copy)  # compute residual

        # Solve for optimal G(.)(k) with diversity constraints
        F_k = F[k, :]
        x_star_G = linalg.lstsq(R.T, F_k.T)[0].T
        x_G = cvx.Variable(x_star_G.shape[0])

        objective_G = cvx.Minimize(cvx.norm2(x_star_G - x_G))

        constraints_G = [x_G >= 0]
        for j in xrange(r):
            if j != k:
                constraints_G += [x_G.T * G[:, j] <= err_V]

        prob_G = cvx.Problem(objective_G, constraints_G)
        result = prob_G.solve(solver='SCS')
        if not np.isinf(result):
            G_k_min = np.asarray(x_G.value)
            G[:, k] = G_k_min[:, 0]
        else:
            print result

        # Solve for optimal F(k)(.) with diversity constraints
        G_k = G[:, k]
        x_star_F = linalg.lstsq(R, G_k)[0]
        x_F = cvx.Variable(x_star_F.shape[0])
        objective_F = cvx.Minimize(cvx.norm2(x_star_F - x_F))

        constraints_F = [x_F >= 0]
        for j in xrange(r):
            if j != k:
                constraints_F += [x_F.T * F[j, :] <= err_F]

        prob_F = cvx.Problem(objective_F, constraints_F)
        result = prob_F.solve(solver='SCS')
        if not np.isinf(result):
            F_k_min = np.asarray(x_F.value)
            F[k, :] = F_k_min[0, :]
        else:
            print result

    return G, F

    def add_data(self, v):
        mask_c = np.isnan(v)
        mask = 1 - mask_c
        U = self.U

        n = self.n
        d = self.d

        Ov = v[mask==1]
        OU = U[mask==1,:]

        w, _, __, ___ = la.lstsq(OU, Ov)
        p = U.dot(w)
        r = np.zeros((n,))
        r[mask==1] = Ov - p[mask==1]

        sigma = la.norm(r) * la.norm(p)
        eta = self.eta0 / self.it

        pw = la.norm(p) * la.norm(w)
        rw = la.norm(r) * la.norm(w)

        if pw == 0 or rw == 0: return

        U = U + (np.cos(sigma * eta) - 1.0) * np.outer(p, w) / pw \
              + np.sin(sigma * eta) * np.outer(r, w) / rw

        self.U = U
        self.it += 1.0

def calibrate_ensemble(models_forecasts, measurements, forecast_len = 48):
    """Calculates coefficient for models in ensembles usulg OLS.
       Returns coefficients for all possible ensembles obtained by models combinations.
    """
    models_count = len(models_forecasts)

    predictors = [list() for mdl in range(models_count)]
    target = list()

    rng = min(len(model) for model in models_forecasts)
    for tm in range(rng):
        msm = measurements[tm * 6: tm * 6 + forecast_len]
        msm_len = len(msm)
        for current_prediction, predictor \
                in zip([prd[tm] for prd in models_forecasts], predictors):
            predictor.extend(current_prediction[:msm_len])
        target.extend(msm)

    ensembles = list()
    for ens_map in reversed(list(product([1,0], repeat = models_count))):
        ensemble_predictors = \
                    [[a*b for a,b in zip(point, ens_map)] for point in zip(*predictors)]

        ensemble_predictors = [pred + [1] for pred in ensemble_predictors]
        coefs = list(lstsq(ensemble_predictors, target)[0])
        ensembles.append(coefs)
    return ensembles

def calibrate(x, y, z, sensor_type):
  H = numpy.array([x, y, z, -y**2, -z**2, numpy.ones([len(x), 1])])
  H = numpy.transpose(H)
  w = x**2

  (X, residues, rank, shape) = linalg.lstsq(H, w)

  OSx = X[0] / 2
  OSy = X[1] / (2 * X[3])
  OSz = X[2] / (2 * X[4])

  A = X[5] + OSx**2 + X[3] * OSy**2 + X[4] * OSz**2
  B = A / X[3]
  C = A / X[4]

  SCx = numpy.sqrt(A)
  SCy = numpy.sqrt(B)
  SCz = numpy.sqrt(C)

  # type conversion from numpy.float64 to standard python floats
  offsets = [OSx, OSy, OSz]
  scale = [SCx, SCy, SCz]

  offsets = map(numpy.asscalar, offsets)
  scale = map(numpy.asscalar, scale)

  #misalignment matrix
  if(sensor_type == "mag"):
    mis_matrix = calibrate_misalignment(x_file, y_file, z_file)
    return (offsets, scale, mis_matrix)
  else:
    return (offsets, scale)

  def run(self, genes):
    """ Computes cross validations for given set of clusters. """
    from numpy import dot, array
    from numpy.linalg import lstsq

    # creates boolean array from genes.
    genes = array([i in genes for i in xrange(self.pis.shape[1])])
    # selects pis for given set of clusters.
    pis = self.pis[:,genes]

    # loops over cross-validation sets.
    scores, fits = [], []
    for croset, fitset in zip(self.crosets, self.fitsets): 
      # pis, energies in cross-validation set.
      cropis = pis[croset, :]
      croene = self.energies[croset]
      # pis, energies in fitting set.
      fitpis = pis[fitset,:]
      fitene = self.energies[fitset]
      # performs fitting.
      try: interactions = lstsq(fitpis, fitene)
      except:
        print "Encountered error in least-square fit."
        print fitpis.shape, fitene.shape, self.pis.shape, genes
        raise
      else: interactions = interactions[0]
      scores.append(dot(cropis, interactions) - croene)
      fits.append(dot(fitpis, interactions) - fitene)
      
    return array(scores), array(fits)

 def trainClassifer(self,labels,vectors,ilog=None):
     '''
     Train the polynomial.  Do not call this function
     manually, instead call the train function on the super
     class.
     '''
     #build matrix
     matrix = []
     for each in vectors:
         if len(each) != 2:
             raise ValueError("ERROR: Vector length=%d.  Polynomial2D only predicts for vectors of length 2."%len(each))
         x,y = each
         matrix.append(self.buildRow(x,y))
     
     matrix = array(matrix)
     labels = array(labels)
     
     x,resids,rank,s = lstsq(matrix,labels)
     
     self.x = x
     self.resids = resids
     self.rank = rank
     self.s = s
     
     if rank != matrix.shape[1]:
         print "WARNING: Polynomial is not fully constrained."

def __loglinregression(rs, zs):
    coef = linalg.lstsq(c_[log(rs), (1,)*len(rs)], log(zs))[0]
    a, b = coef
    #print 'Regression: log(z) = %f*log(r) + %f' % (a,b)
    if a > -1.0:
       print 'Warning: slope is > -1.0'
    return a, b

def least_squares_fit(M, S):
    """
    Least squares fit to fit two datastreams to create calibration matrix C.
    Args:

    Returns:
        Calibration matrix C.
    """
    S_squared = S * S
    S_cubed = S * S * S
    S_24 = np.zeros((S.shape[0], 8, 3))

    for i in xrange(S.shape[0]):
        for j in xrange(S.shape[1]):
            thirds = np.zeros((3, ))
            thirds[0], thirds[1], thirds[2] = S[i][j], S_squared[i][
                j], S_cubed[i][j]
            S_24[i][j] = thirds

    S_24 = S_24.reshape((S_24.shape[0], 24))
    print "=== S_24 shape:{0} ===".format(S_24.shape)

    #C is the vector we are solving for in the S_24 * C = M equation.
    C = lstsq(S_24, M)[0]
    print "=== least squares fit DONE ==="
    print "=== C shape: {0} ===".format(C.shape)
    return C