def LambdaEstimate(Array=scipy.array,
                   Filter=True):
    #===========================================================================
    # copied from R-GenABEL.lamdaest
    #===========================================================================
    Array      = Array.astype(float)
    Estimate   = None
    Ntp        = len(Array)
    QChi2Array = None

    if(Array.max()<=1.0):
#       Convert to quantile function of PValObsArray (df=1) if input are p-values.
        QChi2Array = scipy.stats.chi2.isf(Array,
                                          1)
    else:
        QChi2Array = Array

    if(Filter):
        FilterArray        = (QChi2Array>=1.0e-8)
        QChi2FilteredArray = scipy.compress(FilterArray,QChi2Array)
    else:
        QChi2FilteredArray = QChi2Array
    QChi2FilteredArray.sort()

    PPointsArray = GetPPointsArray(len(QChi2FilteredArray))
    PPointsArray = scipy.sort(scipy.stats.chi2.ppf(PPointsArray,1)) # df=1

    FilterArray          = (PPointsArray!=0.0)
    FilterArray         *= (QChi2FilteredArray!=0.0)
    PPointsArray         = scipy.compress(FilterArray,PPointsArray)
    QChi2FilteredArray   = scipy.compress(FilterArray,QChi2FilteredArray)

#   Fit PPointsArray,QChi2FilteredArray to the linear model.
    P0           = [1.0]
    PBest        = scipy.optimize.leastsq(Residuals,
                                          P0,
                                          args=(PPointsArray,QChi2FilteredArray),
                                          full_output=1,
                                          maxfev=100)
    Estimate = None
    if(type(PBest[0])==scipy.float64):
        Estimate = PBest[0]
    else:
        Estimate     = PBest[0][0]
#   Error estimation of parameter.
    Chi2 = scipy.power(PBest[2]['fvec'],2.0).sum()
    Dof  = len(QChi2FilteredArray)-len(P0)-1
    SE   = scipy.real(scipy.sqrt(PBest[1][0,0])*scipy.sqrt(Chi2/float(Dof)))

    return Estimate,SE

    def downsample(self,freq,factor,ranges=[]):
        if shape(factor) and ranges:
#            pdb.set_trace()
            mask=CreateVariableMask(freq, ranges, factor)
        else:
            mask=CreateDownSampleMask(freq, factor)
#            mask=CreateLogMask(freq)
        dsfreq=compress(mask,freq)
        mylist=['mag','phase','coh']
        for item in mylist:
            tempvect=getattr(self,item)
            if tempvect is not None:
                tempvect=colwise(tempvect)
                setattr(self,item,compress(mask,tempvect,0))
        return dsfreq

def null(A, eps=1e-12):
    '''Compute a base of the null space of A.'''
    u, s, vh = np.linalg.svd(A)
    padding = max(0,np.shape(A)[1]-np.shape(s)[0])
    null_mask = np.concatenate(((s <= eps), np.ones((padding,),dtype=bool)),axis=0)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

    def __init__(self, idx, theta, diff_matrix, log_file):
        self.log_file = log_file
        self.log("create beta form")
        self.idx = idx
        self.theta = theta
        self.theta_norm = sp.linalg.norm(theta, ord=None)
        self.diff_matrix = diff_matrix

        # Find the null space for the subsetted diff matrix
        start = time.time()
        zero_theta_idx = self._get_zero_theta_indices(diff_matrix * theta)
        u, s, v = sp.linalg.svd(diff_matrix[zero_theta_idx,:])
        self.log("SVD done %f" % (time.time() - start))
        null_mask = np.ones(v.shape[1])
        null_mask[:s.size] = s <= self.eps
        null_space = sp.compress(null_mask, v, axis=0)
        null_matrix = np.matrix(sp.transpose(null_space))
        start = time.time()
        beta, istop, itn, normr, normar, norma, conda, normx = sp.sparse.linalg.lsmr(null_matrix, theta.A1, atol=self.eps, btol=self.eps)
        self.log("sp.sparse.linalg.lsmr done %f, istop %d, itn %d" % ((time.time() - start), istop, itn))
        self.beta = np.matrix(beta).T
        self.u = null_matrix

        # Check that we reformulated theta but it is still very close to the original theta
        # assert(res.size == 0 or res < self.eps)
        if sp.linalg.norm(self.u * self.beta - self.theta, ord=2) > self.eps:
            self.log("Warning: Reformulation is off: diff %f" % sp.linalg.norm(self.u * self.beta - self.theta, ord=2))
        self.log("create beta form success")

def nullspace(A, atol=1e-9):
    '''Compute an approximate basis for the nullspace of A using the singular
    value decomposition of `A`.

    Parameters
    ----------
    'A' = ndarray;  A should be at most 2-D.  A 1-D array with length k will be
     treated  as a 2-D with shape (1, k)
    'atol' = float; The absolute tolerance for a zero singular value.  Singular
     values smaller than `atol` are considered to be zero.
    'rtol' = float; The relative tolerance.  Singular values less than
     rtol*smax are considered to be zero, where smax is the largest singular
     value.

    If both `atol` and `rtol` are positive, the combined tolerance is the
    maximum of the two; that is::
    tol = max(atol, rtol * smax)
    Singular values smaller than `tol` are considered to be zero.

    Returns
    -------
    'ns' = ndarray; If `A` is an array with shape (m, k), then `ns` will be an
    array with shape (k, n), where n is the estimated dimension of the
    nullspace of `A`.  The columns of `ns` are a basis for the
    nullspace; each element in numpy.dot(A, ns) will be approximately
    zero.
    '''

    # singular value decomposition
    u, s, vh = sp_la.svd(A)
    null_mask = (s <= atol)
    null_space = sp.compress(null_mask, vh, axis=0)
    return sp.transpose(null_space)

def null(mat, eps=1e-12):
    '''returns null space of matrix mat'''
    u, s, vh = scipy.linalg.svd(mat)   # , full_matrices=False)
    padding = max(0, np.shape(mat)[1]-np.shape(s)[0])
    null_mask = np.concatenate(((s <= eps), np.ones((padding, ), dtype=bool)), axis=0)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

def _sampling_matrix(hessian, cutoff=0, temperature=1, step_scale=1):    
    # basically need SVD of hessian - singular values and eigenvectors
    # hessian = u * diag(singVals) * vh
    u, sing_vals, vh = scipy.linalg.svd(0.5 * hessian)

    # scroll through the singular values and find the ones whose inverses will
    # be huge and set them to zero also, load up the array of singular values 
    # that we store
    # cutoff = (1.0/_.singVals[0])*1.0e03
    # double cutoff = _.singVals[0]*1.0e-02
    cutoff_sing_val = cutoff * max(sing_vals)

    D = 1.0/scipy.maximum(sing_vals, cutoff_sing_val)

    ## now fill in the sampling matrix ("square root" of the Hessian)
    ## note that sqrt(D[i]) is taken here whereas Kevin took sqrt(D[j])
    ## this is because vh is the transpose of his PT -JJW
    samp_mat = scipy.transpose(vh) * scipy.sqrt(D)

    # Divide the sampling matrix by an additional factor such
    # that the expected quadratic increase in cost will be about 1.
    cutoff_vals = scipy.compress(sing_vals < cutoff_sing_val, sing_vals)
    if len(cutoff_vals):
        scale = scipy.sqrt(len(sing_vals) - len(cutoff_vals)
                           + sum(cutoff_vals)/cutoff_sing_val)
    else:
        scale = scipy.sqrt(len(sing_vals))

    samp_mat /= scale
    samp_mat *= step_scale
    samp_mat *= scipy.sqrt(temperature)

    return samp_mat

 def SetMeanMetaboliteConcentrationAndStdExcludingMissingValues(self):
     boFloats  = scipy.array(list(type(Entry)==float for Entry in self.DataArray))
     TmpDataArray = scipy.compress(boFloats,scipy.array(self.DataArray))
     TmpDataArray = scipy.array(TmpDataArray,dtype=float)
     self.MeanMetaboliteConcentrationExcludingMissingValues = TmpDataArray.mean()
     self.StdMetaboliteConcentrationExlcudingMissingValues  = TmpDataArray.std()
     return

def null(A, eps=1e-15):
  import scipy
  from scipy import matrix
  A = matrix(A)
  u, s, vh = scipy.linalg.svd(A)
  null_mask = (s <= eps)
  null_space = scipy.compress(null_mask, vh, axis=0)
  return scipy.transpose(null_space)

def null(A, eps=1e-15):
    """
    Compute nullspace of A.  Thanks Robert Kern and Ryan Krauss:
    http://stackoverflow.com/questions/5889142/python-numpy-scipy-finding-the-null-space-of-a-matrix
    """
    u, s, vh = la.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

def null(A, eps=1e-15):
    """Returns the null-space of the matrix A
    Implementation from
    http://stackoverflow.com/questions/5889142/python-numpy-scipy-finding-the-null-space-of-a-matrix
    """
    u, s, vh = linalg.svd(A)
    null_mask = (s < eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

def main():
    parser = OptionParser(usage=usage)
    parser.add_option(
        "-e", "--epsilon", type=float, dest="eps", default=None, help="set drop-off tolerance [default: %default]"
    )
    parser.add_option("-o", metavar="figname", dest="figname", default=None, help="save the figure to figname")
    parser.add_option(
        "-t",
        "--transparent",
        action="store_true",
        dest="transparent",
        default=False,
        help="save the figure as transparent",
    )
    parser.add_option(
        "-n", "--no-show", action="store_true", dest="no_show", default=False, help="do not show the figure"
    )
    (options, args) = parser.parse_args()
    if len(args) < 1:
        print usage
        return
    filename = args[0]
    print filename + ":"

    fd = open(filename, "r")
    n_row, n_col = map(int, fd.readline().split())
    n_item = int(fd.readline())
    print n_row, n_col, n_item

    ij = nm.zeros((n_item, 2), nm.int32)
    val = nm.zeros((n_item,), nm.float64)
    for ii, row in enumerate(fd.readlines()):
        aux = row.split()
        ij[ii] = int(aux[0]), int(aux[1])
        val[ii] = float(aux[2])

    if options.eps is not None:
        print "using", options.eps
        ij = nm.compress(nm.absolute(val) > options.eps, ij, 0)
        n_item = ij.shape[0]
    else:
        print "showing all"

    print n_item
    if n_item:
        plot(ij[:, 1] + 0.5, ij[:, 0] + 0.5, linestyle="None", marker=",", markersize=0.5, markeredgewidth=0.1)
    axis([-0.5, n_row + 0.5, -0.5, n_col + 0.5])
    axis("image")
    xlabel("%d x %d: %d nnz, %.2f\%% fill" % (n_row, n_col, n_item, 100.0 * n_item / float(n_row * n_col)))
    gca().set_ylim(gca().get_ylim()[::-1])

    if options.figname is not None:
        savefig(options.figname, transparent=options.transparent)

    if not options.no_show:
        show()

def get_intervals(traj):
    # We want to break up our integrals when events fire, so first we figure out
    #  when they fired by looking for duplicated times in the trajectory
    times = traj.get_times()
    eventIndices = scipy.compress(scipy.diff(times) == 0, 
                                  scipy.arange(len(times)))
    intervals = zip([0] + list(eventIndices + 1), 
                    list(eventIndices + 1) + [len(times)])

    return intervals

def null(A, eps=1e-12,sparse=True):
    if sparse == True:
        X = scipy.sparse.csc_matrix(A)
        n=X.shape[1]
        u, s, vh = svds(X, n-1, which='SM')
    else:
        u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

    def FilterSEs(self,
                  XmlObj=lxml.etree._ElementTree,
                  DCs=DataContainer.DataContainers,
                  ColumnTag=str,
                  boDryRun=False):

        FiltersTag = None
        if(boDryRun):
            FiltersTag = 'DryRunFilters'
        else:
            FiltersTag = 'Filters'

        FilterTags  = XmlObj.getroot().find('MtbGWAColumns').find(ColumnTag).find(FiltersTag).text.split(',')

        self.InitFilterReportDictDict(ColumnTag)

        FilterArray = None
        for Tag in FilterTags:
            Operator     = XmlObj.getroot().find('QCFilters').find(Tag).find('Operator').text
            CompareValue = XmlObj.getroot().find('QCFilters').find(Tag).find('Compare').text
            ValueType    = XmlObj.getroot().find('QCFilters').find(Tag).find('CompareType').text

            FFunction = FilterFunction.FilterFunction(OperatorString=Operator,
                                                      CompareString=CompareValue,
                                                      CompareType=ValueType)

            InitLength  = len(DCs.DataContainers[ColumnTag].GetDataArray())
            FinalLength = None
            FilterArray = FFunction.Run(DataArray=DCs.DataContainers[ColumnTag].GetDataArray())

            if(not boDryRun):
                for Key in DCs.DataContainers.iterkeys():
                    DataArray = scipy.compress(FilterArray,
                                               DCs.DataContainers[Key].GetDataArray())
                    DCs.DataContainers[Key].ReplaceDataArray(DataArray)
                FinalLength = len(DCs.DataContainers[ColumnTag].GetDataArray())
            else:
                CounterDict = collections.defaultdict(int)
                for Entry in FilterArray:
                    CounterDict[Entry] += 1
                Difference  = CounterDict[False]
                FinalLength = InitLength-Difference


            self.SetFilterReportDictDict(ParentTag=ColumnTag,
                                         ChildTag=Tag,
                                         Value='Column Tag   = '+ColumnTag+'\n'+\
                                               'Filter Tag   = '+Tag+'\n'+\
                                               'Start Length = '+str(InitLength)+'\n'+\
                                               'Final Length = '+str(FinalLength)+'\n'+\
                                               'Difference   = '+str(InitLength-FinalLength))

        return DCs,\
               FilterTags

def CompleteBase(V, B, eps=1e-4):
    tbase = append(V, B, axis=1)
    p, l, u = lu(tbase)
    echelon = zeros(u.shape[1], int)
    
    for row in u:
        tmp = nonzero(abs(row) > eps)[0]
        if tmp.size:
            echelon[tmp[0]] = 1
  
    return compress(echelon, tbase, axis=1)

def null(A, eps=1e-6):
    u, s, vh = numpy.linalg.svd(A,full_matrices=1,compute_uv=1)
    null_rows = [];
    rank = numpy.linalg.matrix_rank(A)
    for i in range(A.shape[1]):
        if i<rank:
            null_rows.append(False);
        else:
            null_rows.append(True);
    null_space = scipy.compress(null_rows, vh, axis=0)
    return null_space.T

def null(A, eps=1e-10):
    """
    null-space of a Matrix or 2d-array
    """
    n, m = A.shape
    if n > m :
        return null(A.T, eps).T
        #return null(scipy.transpose(A), eps)
    u, s, vh = sc.linalg.svd(A)
    s=sc.append(s,sc.zeros(m))[0:m]
    null_mask = (s <= eps)
    null_space = sc.compress(null_mask, vh, axis=0)
    return null_space.T

def compute_stationary_dist(Pt):
    eps = 1e-15
    u,s,vh = linalg.svd((np.eye(Pt.shape[0]) - Pt).T)
    null_mask = (s<eps)
    null_space = scipy.compress(null_mask, vh, axis = 0)
    #print 'nullspace:',null_space
    Pi = null_space[0]/null_space[0].sum()
    if (Pi<0).sum()>0 or null_space.shape[0]>1:
        print Pt
        print null_space
        print Pi
    assert (Pi<0).sum()==0
    return Pi

def _sampling_matrix(hessian, cutoff=None, diag = None, temperature=1, step_scale=1):
    # basically need SVD of hessian - singular values and eigenvectors
    # hessian = u * diag(singVals) * vh
    #u, sing_vals, vh = scipy.linalg.svd(hessian)

    # scroll through the singular values and find the ones whose inverses will
    # be huge and set them to zero also, load up the array of singular values 
    # that we store
    # cutoff = (1.0/_.singVals[0])*1.0e03
    # double cutoff = _.singVals[0]*1.0e-02
    # when cutoff is set to zero it means that all values are included
    # cutoff*(sloppiest eigenvalue)
    
    if cutoff:
        u, sing_vals, vh = scipy.linalg.svd(hessian)
        cutoff_sing_val = cutoff * max(sing_vals)
        #when cutoff is set to zero it means that all values are included
        D = 1.0/scipy.maximum(sing_vals, cutoff_sing_val)
        samp_mat = scipy.transpose(vh)*scipy.sqrt(D)
    # instead of cutoff use another method, adding diagonal term to hessian
    elif diag is not None:
        u, sing_vals, vh = scipy.linalg.svd(hessian+diag)
        D = 1.0/sing_vals
        samp_mat = scipy.transpose(vh)*scipy.sqrt(D)
        cutoff_sing_val = diag[0,0]
    else: 
        u, sing_vals, vh = scipy.linalg.svd(hessian)
        D = 1.0/sing_vals
        samp_mat = scipy.transpose(vh)*scipy.sqrt(D)
        cutoff_sing_val = 0

    ## now fill in the sampling matrix ("square root" of the Hessian)
    ## note that sqrt(D[i]) is taken here whereas Kevin took sqrt(D[j])
    ## this is because vh is the transpose of his PT -JJW
    #samp_mat = scipy.transpose(vh) * scipy.sqrt(D)

    # Divide the sampling matrix by an additional factor such
    # that the expected quadratic increase in cost will be about 1.
    cutoff_vals = scipy.compress(sing_vals < cutoff_sing_val, sing_vals)
    if len(cutoff_vals):
        scale = scipy.sqrt(len(sing_vals) - len(cutoff_vals)
                           + sum(cutoff_vals)/cutoff_sing_val)
    else:
        scale = scipy.sqrt(len(sing_vals))

    samp_mat /= scale
    samp_mat *= step_scale
    samp_mat *= scipy.sqrt(temperature)

    return samp_mat