Python numpy.asarray_chkfinite函数代码示例

    def __init__(self, correct, student_ids=None, item_ids=None, student_idx=None,
                 item_idx=None, is_held_out=None, num_students=None, num_items=None,
                 **bn_learner_kwargs):
        """
        :param np.ndarray[bool] correct: a 1D array of correctness values
        :param np.ndarray|None student_ids: student identifiers for each interaction; if no student
            indices provided, sort order of these ids determines theta indices.
        :param np.ndarray|None item_ids: item identifiers for each interaction; if no item indices
            are provided, sort order of these ids determines item indices.
        :param np.ndarray[int]|None student_idx: a 1D array mapping `correct` to student index
        :param np.ndarray[int]|None item_idx: a 1D array mapping `correct` to item index
        :param np.ndarray[bool] is_held_out: a 1D array indicating whether the interaction should be
            held out from training (if not all zeros, a held_out test node will be added to learner)
        :param int|None num_students: optional number of students. Default is one plus
            the maximum index.
        :param int|None num_items: optional number of items. Default is one plus
            the maximum index.
        :param bn_learner_kwargs: arguments to be passed on to the BayesNetLearner init
        """
        # convert pandas Series to np.ndarray and check argument dimensions
        correct = np.asarray_chkfinite(correct, dtype=bool)
        student_ids, student_idx = check_and_set_idx(student_ids, student_idx, 'student')
        item_ids, item_idx = check_and_set_idx(item_ids, item_idx, 'item')

        if len(correct) != len(student_idx) or len(correct) != len(item_idx):
            raise ValueError("number of elements in correct ({}), student_idx ({}), and item_idx"
                             "({}) must be the same".format(len(correct), len(student_idx),
                                                            len(item_idx)))
        if is_held_out is not None and (
                len(is_held_out) != len(correct) or is_held_out.dtype != bool):
            raise ValueError("held_out ({}) must be None or an array of bools the same length as "
                             "correct ({})".format(len(is_held_out), len(correct)))

        self.num_students = set_or_check_min(num_students, np.max(student_idx) + 1, 'num_students')
        self.num_items = set_or_check_min(num_items, np.max(item_idx) + 1, 'num_items')

        theta_node = DefaultGaussianNode(THETAS_KEY, self.num_students, ids=student_ids)
        offset_node = DefaultGaussianNode(OFFSET_COEFFS_KEY, self.num_items, ids=item_ids)
        nodes = [theta_node, offset_node]

        # add response nodes (train/test if there is held-out data; else just the train set)
        if is_held_out is not None and np.sum(is_held_out):
            if np.sum(is_held_out) == len(is_held_out):
                raise ValueError("some interactions must be not held out")
            is_held_out = np.asarray_chkfinite(is_held_out, dtype=bool)
            node_names = (TRAIN_RESPONSES_KEY, TEST_RESPONSES_KEY)
            response_idxs = (np.logical_not(is_held_out), is_held_out)
        else:
            node_names = (TRAIN_RESPONSES_KEY,)
            response_idxs = (np.ones_like(correct, dtype=bool),)
        for node_name, response_idx in zip(node_names, response_idxs):
            cpd = OnePOCPD(item_idx=item_idx[response_idx], theta_idx=student_idx[response_idx],
                           num_thetas=self.num_students, num_items=self.num_items)
            param_nodes = {THETAS_KEY: theta_node, OFFSET_COEFFS_KEY: offset_node}
            nodes.append(Node(name=node_name, data=correct[response_idx], cpd=cpd,
                              solver_pars=SolverPars(learn=False), param_nodes=param_nodes,
                              held_out=(node_name == TEST_RESPONSES_KEY)))

        # store leaf nodes for learning
        super(OnePOLearner, self).__init__(nodes=nodes, **bn_learner_kwargs)

def check_and_set_idx(ids, idx, prefix):
    """ Reconciles passed-in IDs and indices and returns indices, as well as unique IDs
    in the order specified by the indices.  If only IDs supplied, returns the sort-arg
    as the index.  If only indices supplied, returns None for IDs.  If both supplied,
    checks that the correspondence is unique and returns unique IDs in the sort order of
    the associated index.
    :param np.ndarray ids: array of IDs
    :param np.ndarray[int] idx: array of indices
    :param str prefix: variable name (for error logging)
    :return: unique IDs and indices (passed in or derived from the IDs)
    :rtype: np.ndarray, np.ndarray
    """
    if ids is None and idx is None:
        raise ValueError('Both {}_ids and {}_idx cannot be None'.format(prefix, prefix))
    if ids is None:
        return None, np.asarray_chkfinite(idx)
    if idx is None:
        return np.unique(ids, return_inverse=True)
    else:
        ids = np.asarray(ids)
        idx = np.asarray_chkfinite(idx)
        if len(idx) != len(ids):
            raise ValueError('{}_ids ({}) and {}_idx ({}) must have the same length'.format(
                prefix, len(ids), prefix, len(idx)))
        uniq_idx, idx_sort_index = np.unique(idx, return_index=True)
        # make sure each unique index corresponds to a unique id
        if not all(len(set(ids[idx == i])) == 1 for i in uniq_idx):
            raise ValueError("Each index must correspond to a unique {}_id".format(prefix))
        return ids[idx_sort_index], idx

def orthogonal_procrustes(A, B, check_finite=True):
    """
    Compute the matrix solution of the orthogonal Procrustes problem.

    Given matrices A and B of equal shape, find an orthogonal matrix R
    that most closely maps A to B [1]_.
    Note that unlike higher level Procrustes analyses of spatial data,
    this function only uses orthogonal transformations like rotations
    and reflections, and it does not use scaling or translation.

    Parameters
    ----------
    A : (M, N) array_like
        Matrix to be mapped.
    B : (M, N) array_like
        Target matrix.
    check_finite : bool, optional
        Whether to check that the input matrices contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.

    Returns
    -------
    R : (N, N) ndarray
        The matrix solution of the orthogonal Procrustes problem.
        Minimizes the Frobenius norm of dot(A, R) - B, subject to
        dot(R.T, R) == I.
    scale : float
        Sum of the singular values of ``dot(A.T, B)``.

    Raises
    ------
    ValueError
        If the input arrays are incompatibly shaped.
        This may also be raised if matrix A or B contains an inf or nan
        and check_finite is True, or if the matrix product AB contains
        an inf or nan.

    References
    ----------
    .. [1] Peter H. Schonemann, "A generalized solution of the orthogonal
           Procrustes problem", Psychometrica -- Vol. 31, No. 1, March, 1996.

    """
    if check_finite:
        A = np.asarray_chkfinite(A)
        B = np.asarray_chkfinite(B)
    else:
        A = np.asanyarray(A)
        B = np.asanyarray(B)
    if A.ndim != 2:
        raise ValueError('expected ndim to be 2, but observed %s' % A.ndim)
    if A.shape != B.shape:
        raise ValueError('the shapes of A and B differ (%s vs %s)' % (
            A.shape, B.shape))
    # Be clever with transposes, with the intention to save memory.
    u, w, vt = svd(B.T.dot(A).T)
    R = u.dot(vt)
    scale = w.sum()
    return R, scale

    def __init__( self, griddata, lo, hi, maps=None, copy=True, verbose=1,
            order=1, prefilter=False ):
        griddata = np.asanyarray( griddata )
        dim = griddata.ndim # - (griddata.shape[-1] == 1) # ??
        assert dim >= 2, griddata.shape
        self.dim = dim
        if np.isscalar(lo):
            lo *= np.ones(dim)
        if np.isscalar(hi):
            hi *= np.ones(dim)
        self.loclip = lo = np.asarray_chkfinite( lo ).copy()
        self.hiclip = hi = np.asarray_chkfinite( hi ).copy()
        assert lo.shape == (dim,), lo.shape
        assert hi.shape == (dim,), hi.shape
        self.copy = copy
        self.verbose = verbose
        self.order = order
        if order > 1 and 0 < prefilter < 1: # 1/3: Mitchell-Netravali = 1/3 B + 2/3 fit
            exactfit = spline_filter( griddata ) # see Unser
            griddata += prefilter * (exactfit - griddata)
            prefilter = False
        self.griddata = griddata
        self.prefilter = (prefilter == True)

        if maps is None:
            maps = [None,] * len(lo)

        self.maps = maps
        self.nmap = 0
        if len(maps) > 0:
            assert len(maps) == dim, "maps must have len %d, not %d" % (
                    dim, len(maps))
            # linear maps (map None): Xcol -= lo *= scale -> [0, n-1]
            # nonlinear: np.interp e.g. [50 52 62 63] -> [0 1 2 3]
            self._lo = np.zeros(dim)
            self._scale = np.ones(dim)

            for j, (map, n, l, h) in enumerate( zip( maps, griddata.shape, lo, hi )):
                ## print "test: j map n l h:", j, map, n, l, h
                if map is None or callable(map):
                    self._lo[j] = l
                    if h > l:
                        self._scale[j] = (n - 1) / (h - l) # _map lo -> 0, hi -> n - 1
                    else:
                        self._scale[j] = 0 # h <= l: X[:,j] -> 0
                    continue
                self.maps[j] = map = np.asanyarray(map)
                self.nmap += 1
                assert len(map) == n, "maps[%d] must have len %d, not %d" % (
                    j, n, len(map) )
                mlo, mhi = map.min(), map.max()
                if not (l <= mlo <= mhi <= h):
                    print "Warning: Intergrid maps[%d] min %.3g max %.3g " \
                        "are outside lo %.3g hi %.3g" % (
                        j, mlo, mhi, l, h )

def cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True):
    """Solve the linear equations A x = b, given the Cholesky factorization of A.

    Parameters
    ----------
    (c, lower) : tuple, (array, bool)
        Cholesky factorization of a, as given by cho_factor
    b : array
        Right-hand side
    overwrite_b : bool, optional
        Whether to overwrite data in b (may improve performance)
    check_finite : bool, optional
        Whether to check that the input matrices contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.

    Returns
    -------
    x : array
        The solution to the system A x = b

    See also
    --------
    cho_factor : Cholesky factorization of a matrix

    Examples
    --------
    >>> from scipy.linalg import cho_factor, cho_solve
    >>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])
    >>> c, low = cho_factor(A)
    >>> x = cho_solve((c, low), [1, 1, 1, 1])
    >>> np.allclose(A @ x - [1, 1, 1, 1], np.zeros(4))
    True

    """
    (c, lower) = c_and_lower
    if check_finite:
        b1 = asarray_chkfinite(b)
        c = asarray_chkfinite(c)
    else:
        b1 = asarray(b)
        c = asarray(c)
    if c.ndim != 2 or c.shape[0] != c.shape[1]:
        raise ValueError("The factored matrix c is not square.")
    if c.shape[1] != b1.shape[0]:
        raise ValueError("incompatible dimensions.")

    overwrite_b = overwrite_b or _datacopied(b1, b)

    potrs, = get_lapack_funcs(('potrs',), (c, b1))
    x, info = potrs(c, b1, lower=lower, overwrite_b=overwrite_b)
    if info != 0:
        raise ValueError('illegal value in %d-th argument of internal potrs'
                         % -info)
    return x

def cho_solve_banded(cb_and_lower, b, overwrite_b=False, check_finite=True):
    """Solve the linear equations A x = b, given the Cholesky factorization of A.

    Parameters
    ----------
    (cb, lower) : tuple, (array, bool)
        `cb` is the Cholesky factorization of A, as given by cholesky_banded.
        `lower` must be the same value that was given to cholesky_banded.
    b : array
        Right-hand side
    overwrite_b : bool
        If True, the function will overwrite the values in `b`.
    check_finite : boolean, optional
        Whether to check that the input matrices contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.

    Returns
    -------
    x : array
        The solution to the system A x = b

    See also
    --------
    cholesky_banded : Cholesky factorization of a banded matrix

    Notes
    -----

    .. versionadded:: 0.8.0

    """
    (cb, lower) = cb_and_lower
    if check_finite:
        cb = asarray_chkfinite(cb)
        b = asarray_chkfinite(b)
    else:
        cb = asarray(cb)
        b = asarray(b)

    # Validate shapes.
    if cb.shape[-1] != b.shape[0]:
        raise ValueError("shapes of cb and b are not compatible.")

    pbtrs, = get_lapack_funcs(('pbtrs',), (cb, b))
    x, info = pbtrs(cb, b, lower=lower, overwrite_b=overwrite_b)
    if info > 0:
        raise LinAlgError("%d-th leading minor not positive definite" % info)
    if info < 0:
        raise ValueError('illegal value in %d-th argument of internal pbtrs'
                                                                    % -info)
    return x

def _cholesky(a, lower=False, overwrite_a=False, clean=True,
              check_finite=True):
    """Common code for cholesky() and cho_factor()."""

    a1 = asarray_chkfinite(a) if check_finite else asarray(a)
    a1 = atleast_2d(a1)

    # Dimension check
    if a1.ndim != 2:
        raise ValueError('Input array needs to be 2 dimensional but received '
                         'a {}d-array.'.format(a1.ndim))
    # Squareness check
    if a1.shape[0] != a1.shape[1]:
        raise ValueError('Input array is expected to be square but has '
                         'the shape: {}.'.format(a1.shape))

    # Quick return for square empty array
    if a1.size == 0:
        return a1.copy(), lower

    overwrite_a = overwrite_a or _datacopied(a1, a)
    potrf, = get_lapack_funcs(('potrf',), (a1,))
    c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean)
    if info > 0:
        raise LinAlgError("%d-th leading minor of the array is not positive "
                          "definite" % info)
    if info < 0:
        raise ValueError('LAPACK reported an illegal value in {}-th argument'
                         'on entry to "POTRF".'.format(-info))
    return c, lower

def schur(a,output='real',lwork=None,overwrite_a=0):
    """Compute Schur decomposition of matrix a.

    Description:

      Return T, Z such that a = Z * T * (Z**H) where Z is a
      unitary matrix and T is either upper-triangular or quasi-upper
      triangular for output='real'
    """
    if not output in ['real','complex','r','c']:
        raise ValueError, "argument must be 'real', or 'complex'"
    a1 = asarray_chkfinite(a)
    if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
        raise ValueError, 'expected square matrix'
    typ = a1.dtype.char
    if output in ['complex','c'] and typ not in ['F','D']:
        if typ in _double_precision:
            a1 = a1.astype('D')
            typ = 'D'
        else:
            a1 = a1.astype('F')
            typ = 'F'
    overwrite_a = overwrite_a or (_datanotshared(a1,a))
    gees, = get_lapack_funcs(('gees',),(a1,))
    if lwork is None or lwork == -1:
        # get optimal work array
        result = gees(lambda x: None,a,lwork=-1)
        lwork = result[-2][0]
    result = gees(lambda x: None,a,lwork=result[-2][0],overwrite_a=overwrite_a)
    info = result[-1]
    if info<0: raise ValueError,\
       'illegal value in %-th argument of internal gees'%(-info)
    elif info>0: raise LinAlgError, "Schur form not found.  Possibly ill-conditioned."
    return result[0], result[-3]

def lu(a,permute_l=0,overwrite_a=0):
    """Return LU decompostion of a matrix.

    Inputs:

      a     -- An M x N matrix.
      permute_l  -- Perform matrix multiplication p * l [disabled].

    Outputs:

      p,l,u  -- LU decomposition matrices of a [permute_l=0]
      pl,u   -- LU decomposition matrices of a [permute_l=1]

    Definitions:

      a = p * l * u    [permute_l=0]
      a = pl * u       [permute_l=1]

      p   -  An M x M permutation matrix
      l   -  An M x K lower triangular or trapezoidal matrix
             with unit-diagonal
      u   -  An K x N upper triangular or trapezoidal matrix
             K = min(M,N)
    """
    a1 = asarray_chkfinite(a)
    if len(a1.shape) != 2:
        raise ValueError, 'expected matrix'
    overwrite_a = overwrite_a or (_datanotshared(a1,a))
    flu, = get_flinalg_funcs(('lu',),(a1,))
    p,l,u,info = flu(a1,permute_l=permute_l,overwrite_a = overwrite_a)
    if info<0: raise ValueError,\
       'illegal value in %-th argument of internal lu.getrf'%(-info)
    if permute_l:
        return l,u
    return p,l,u

def rq(a, overwrite_a=False, lwork=None):
    """Compute RQ decomposition of a square real matrix.

    Calculate the decomposition :lm:`A = R Q` where Q is unitary/orthogonal
    and R upper triangular.

    Parameters
    ----------
    a : array, shape (M, M)
        Square real matrix to be decomposed
    overwrite_a : boolean
        Whether data in a is overwritten (may improve performance)
    lwork : integer
        Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
        is computed.
    econ : boolean

    Returns
    -------
    R : double array, shape (M, N) or (K, N) for econ==True
        Size K = min(M, N)
    Q : double or complex array, shape (M, M) or (M, K) for econ==True

    Raises LinAlgError if decomposition fails

    """
    # TODO: implement support for non-square and complex arrays
    a1 = asarray_chkfinite(a)
    if len(a1.shape) != 2:
        raise ValueError('expected matrix')
    M,N = a1.shape
    if M != N:
        raise ValueError('expected square matrix')
    if issubclass(a1.dtype.type, complexfloating):
        raise ValueError('expected real (non-complex) matrix')
    overwrite_a = overwrite_a or (_datanotshared(a1, a))
    gerqf, = get_lapack_funcs(('gerqf',), (a1,))
    if lwork is None or lwork == -1:
        # get optimal work array
        rq, tau, work, info = gerqf(a1, lwork=-1, overwrite_a=1)
        lwork = work[0]
    rq, tau, work, info = gerqf(a1, lwork=lwork, overwrite_a=overwrite_a)
    if info < 0:
        raise ValueError('illegal value in %d-th argument of internal geqrf'
                                                                    % -info)
    gemm, = get_blas_funcs(('gemm',), (rq,))
    t = rq.dtype.char
    R = special_matrices.triu(rq)
    Q = numpy.identity(M, dtype=t)
    ident = numpy.identity(M, dtype=t)
    zeros = numpy.zeros

    k = min(M, N)
    for i in range(k):
        v = zeros((M,), t)
        v[N-k+i] = 1
        v[0:N-k+i] = rq[M-k+i, 0:N-k+i]
        H = gemm(-tau[i], v, v, 1+0j, ident, trans_b=2)
        Q = gemm(1, Q, H)
    return R, Q

    def __init__(self, data, segment_iterator, uri=None):
        # make sure data does not contain NaN nor inf
        data = np.asarray_chkfinite(data)

        # make sure segment_iterator is actually one of those
        if not isinstance(segment_iterator, (SlidingWindow, Timeline)):
            raise TypeError("segment_iterator must 'Timeline' or "
                            "'SlidingWindow'.")

        # make sure it iterates the correct number of segments
        try:
            N = len(segment_iterator)
        except Exception:
            # an exception is raised by `len(sliding_window)`
            # in case sliding window has infinite end.
            # this is acceptable, no worry...
            pass
        else:
            n = data.shape[0]
            if n != N:
                raise ValueError("mismatch between number of segments (%d) "
                                 "and number of feature vectors (%d)." % (N, n))

        super(BasePrecomputedSegmentFeature, self).__init__(uri=uri)
        self.__data = data
        self._segment_iterator = segment_iterator

    def find_disjoint_biclusters(self, biclusters_number=50):
        data = np.asarray_chkfinite(self.matrix)
        data[data == 0] = 0.000001
        coclustering = SpectralCoclustering(n_clusters=biclusters_number, random_state=0)
        coclustering.fit(data)

        biclusters = set()
        for i in range(biclusters_number):
            rows, columns = coclustering.get_indices(i)
            row_set = set(rows)
            columns_set = set(columns)
            if len(row_set) > 0 and len(columns_set) > 0:
                density = self._calculate_box_cluster_density(row_set, columns_set)
                odd_columns = set()
                for column in columns_set:
                    col_density = self._calculate_column_density(column, row_set)
                    if col_density < density / 4:
                        odd_columns.add(column)
                columns_set.difference_update(odd_columns)
                if len(columns_set) == 0:
                    continue

                odd_rows = set()
                for row in row_set:
                    row_density = self._calculate_row_density(row, columns_set)
                    if row_density < density / 4:
                        odd_rows.add(row)
                row_set.difference_update(odd_rows)

                if len(row_set) > 0 and len(columns_set) > 0:
                    density = self._calculate_box_cluster_density(row_set, columns_set)
                    biclusters.add(Bicluster(row_set, columns_set, density))

        return biclusters

def det(a, overwrite_a=False):
    """Compute the determinant of a matrix

    Parameters
    ----------
    a : array, shape (M, M)

    Returns
    -------
    det : float or complex
        Determinant of a

    Notes
    -----
    The determinant is computed via LU factorization, LAPACK routine z/dgetrf.
    """
    a1 = asarray_chkfinite(a)
    if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
        raise ValueError('expected square matrix')
    overwrite_a = overwrite_a or _datacopied(a1, a)
    fdet, = get_flinalg_funcs(('det',), (a1,))
    a_det, info = fdet(a1, overwrite_a=overwrite_a)
    if info < 0:
        raise ValueError('illegal value in %d-th argument of internal '
                                                        'det.getrf' % -info)
    return a_det

def orthogonal_procrustes(A, ref_matrix, reflection=False):
	# Adaptation of scipy.linalg.orthogonal_procrustes -> https://github.com/scipy/scipy/blob/v0.16.0/scipy/linalg/_procrustes.py#L14
	# Info here: http://compgroups.net/comp.soft-sys.matlab/procrustes-analysis-without-reflection/896635
	# goal is to find unitary matrix R with det(R) > 0 such that ||A*R - ref_matrix||^2 is minimized
	from scipy.linalg.decomp_svd import svd # Singular Value Decomposition, factors matrices
	from scipy.linalg import det
	import numpy as np

	A = np.asarray_chkfinite(A)
	ref_matrix = np.asarray_chkfinite(ref_matrix)

	if A.ndim != 2:
		raise ValueError('expected ndim to be 2, but observed %s' % A.ndim)
	if A.shape != ref_matrix.shape:
		raise ValueError('the shapes of A and ref_matrix differ (%s vs %s)' % (A.shape, ref_matrix.shape))


	u, w, vt = svd(ref_matrix.T.dot(A).T)

	# Goal: minimize ||A*R - ref||^2, switch to trace
	# trace((A*R-ref).T*(A*R-ref)), now we distribute
	# trace(R'*A'*A*R) + trace(ref.T*ref) - trace((A*R).T*ref) - trace(ref.T*(A*R)), trace doesn't care about order, so re-order
	# trace(R*R.T*A.T*A) + trace(ref.T*ref) - trace(R.T*A.T*ref) - trace(ref.T*A*R), simplify
	# trace(A.T*A) + trace(ref.T*ref) - 2*trace(ref.T*A*R)
	# Thus, to minimize we want to maximize trace(ref.T * A * R) 

	# u*w*v.T = (ref.T*A).T
	# ref.T * A = w * u.T * v
	# trace(ref.T * A * R) = trace (w * u.T * v * R)
	# differences minimized when trace(ref.T * A * R) is maximized, thus when trace(u.T * v * R) is maximized
	# This occurs when u.T * v * R = I (as u, v and R are all unitary matrices so max is 1)
	# R is a rotation matrix so R.T = R^-1
	# u.T * v * I = R^-1 = R.T
	# R = u * v.T
	# Thus, R = u.dot(vt)

	R = u.dot(vt) # Get the rotation matrix, including reflections
	if not reflection and det(R) < 0: # If we don't want reflection
		# To remove reflection, we change the sign of the rightmost column of u (or v) and the scalar associated
		# with that column
		u[:,-1] *= -1
		w[-1] *= -1
		R = u.dot(vt)
	
	scale = w.sum() # Get the scaled difference

	return R,scale

def qr_old(a, overwrite_a=False, lwork=None, check_finite=True):
    """Compute QR decomposition of a matrix.

    Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
    and R upper triangular.

    Parameters
    ----------
    a : array, shape (M, N)
        Matrix to be decomposed
    overwrite_a : boolean
        Whether data in a is overwritten (may improve performance)
    lwork : integer
        Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
        is computed.
    check_finite : boolean, optional
        Whether to check that the input matrix contains only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.

    Returns
    -------
    Q : float or complex array, shape (M, M)
    R : float or complex array, shape (M, N)
        Size K = min(M, N)

    Raises LinAlgError if decomposition fails

    """
    if check_finite:
        a1 = numpy.asarray_chkfinite(a)
    else:
        a1 = numpy.asarray(a)
    if len(a1.shape) != 2:
        raise ValueError('expected matrix')
    M,N = a1.shape
    overwrite_a = overwrite_a or (_datacopied(a1, a))
    geqrf, = get_lapack_funcs(('geqrf',), (a1,))
    if lwork is None or lwork == -1:
        # get optimal work array
        qr, tau, work, info = geqrf(a1, lwork=-1, overwrite_a=1)
        lwork = work[0]
    qr, tau, work, info = geqrf(a1, lwork=lwork, overwrite_a=overwrite_a)
    if info < 0:
        raise ValueError('illegal value in %d-th argument of internal geqrf'
                                                                    % -info)
    gemm, = get_blas_funcs(('gemm',), (qr,))
    t = qr.dtype.char
    R = numpy.triu(qr)
    Q = numpy.identity(M, dtype=t)
    ident = numpy.identity(M, dtype=t)
    zeros = numpy.zeros
    for i in range(min(M, N)):
        v = zeros((M,), t)
        v[i] = 1
        v[i+1:M] = qr[i+1:M, i]
        H = gemm(-tau[i], v, v, 1+0j, ident, trans_b=2)
        Q = gemm(1, Q, H)
    return Q, R