FWHM calculation on a 1D array by using least square fitting of
        a Moffat function on the data.  arr1d is a 1D array cut in either
        X or Y direction on the object.
350139
def calc_fwhm_moffat(self, arr1d, medv=None, moffat_fn=None):
        
        if moffat_fn is None:
            moffat_fn = self.moffat

        N = len(arr1d)
        X = np.array(list(range(N)))
        Y = arr1d
        # Fitting works more reliably if we do the following
        # a. subtract sky background
        if medv is None:
            medv = get_median(Y)
        Y = Y - medv
        maxv = Y.max()
        # b. clamp to 0..max (of the sky subtracted field)
        Y = Y.clip(0, maxv)

        # Fit a moffat
        p0 = [0, N - 1, 2, maxv]              # Inital guess
        # Distance to the target function
        errfunc = lambda p, x, y: moffat_fn(x, p) - y  # noqa
        # Least square fit to the gaussian
        with self.lock:
            # NOTE: without this mutex, optimize.leastsq causes a fatal error
            # sometimes--it appears not to be thread safe.
            # The error is:
            # "SystemError: null argument to internal routine"
            # "Fatal Python error: GC object already tracked"
            p1, success = optimize.leastsq(errfunc, p0[:], args=(X, Y))

        if not success:
            raise IQCalcError("FWHM moffat fitting failed")

        mu, width, power, maxv = p1
        width = np.abs(width)
        self.logger.debug("mu=%f width=%f power=%f maxv=%f" % (
            mu, width, power, maxv))

        fwhm = 2.0 * width * np.sqrt(2.0 ** (1.0 / power) - 1.0)

        # some routines choke on numpy values and need "pure" Python floats
        # e.g. when marshalling through a remote procedure interface
        fwhm = float(fwhm)
        mu = float(mu)
        width = float(width)
        power = float(power)
        maxv = float(maxv)

        res = Bunch.Bunch(fwhm=fwhm, mu=mu, width=width, power=power,
                          maxv=maxv, fit_fn=moffat_fn,
                          fit_args=[mu, width, power, maxv])
        return res
443128
def _get_2e_asym(self):
        
        
        # get needed PPG parameters for splitting 1D histos into 2D histos
        try:        
            # get frequency vector
            freq = np.arange(self.ppg['freq_start'].mean,\
                             self.ppg['freq_stop'].mean+\
                                    self.ppg['freq_incr'].mean,\
                             self.ppg['freq_incr'].mean)
                             
            # number of dwelltimes per frequency bin 
            ndwell = 2*int(self.ppg['ndwell_per_f'].mean)-1
            
            # number of RF on delays for the start bin. 
            start_bin = int(self.ppg['rf_on_delay'].mean)
            
            # get bin centers in ms
            time = self.ppg['rf_on_ms'].mean*(np.arange(ndwell)+0.5-ndwell/2.)
            
            # get the time and index of the middle time 
            mid_time_i = int(np.floor(ndwell/2.))
            mid_time = time[mid_time_i]
        
            # beam off time after pulse in ms
            beam_off = int(self.ppg['beam_off_ms'].mean)
        
        except KeyError:
            raise RuntimeError("Not all dictionary variables read out to "+\
                               "proper locations")
            
        # setup output
        out = bdict()
        out['freq'] = freq
        out['time'] = time
            
        # get data
        data = np.array(self._get_area_data()) # [[fp], [bfm], [bp], [bm]]
        
        # discared initial bad bins, and beam-off trailing bins
        data = data[:,start_bin:len(freq)*ndwell+start_bin]
        
        # split data by frequency
        nsplit = len(data[0])/ndwell
        
        fp = np.array(np.split(data[0],nsplit))
        fm = np.array(np.split(data[1],nsplit))
        bp = np.array(np.split(data[2],nsplit))
        bm = np.array(np.split(data[3],nsplit))
        
        # get raw asymmetries 
        asym_p_2cntr = (bp-fp)/(bp+fp)      # two counter
        asym_m_2cntr = (bm-fm)/(bm+fm)      # two counter
        r = np.sqrt(bp*fm/(bm*fp))
        asym_4cntr = (r-1)/(r+1)            # four counter
        
        # get raw asymmetry errors
        asym_p_2cntr_err = 2*np.sqrt(bp*fp)/((bp+fp)**1.5)
        asym_m_2cntr_err = 2*np.sqrt(bm*fm)/((bm+fm)**1.5)
        asym_4cntr_err = r*np.sqrt(1./bp+1./bm+1./fp+1./fm)/((r+1)**2)

        # save to output
        out['raw_p'] = np.array([asym_p_2cntr,asym_p_2cntr_err])
        out['raw_n'] = np.array([asym_m_2cntr,asym_m_2cntr_err])
        out['raw_c'] = np.array([asym_4cntr,asym_4cntr_err])

        # wrap asymmetry arrays into one for calculations [p,m,4]
        # indexing is now [pm4][freq][time bin]
        asym     = np.array([asym_p_2cntr,    asym_m_2cntr,    asym_4cntr])
        asym_err = np.array([asym_p_2cntr_err,asym_m_2cntr_err,asym_4cntr_err])

        # compute differenced asymmetries via slopes from weighted least squares 
        # minimization.
        if ndwell >= 5:
            
            # calculate needed components element-wise
            w = asym_err**-2
            x = time
            y = asym
            
            wx = w*x
            wy = w*y
            wxy = w*x*y
            wxx = w*x*x
        
            # sum over values i < mid_time_i within each asymmetry and frequency 
            # Indexing: [pm4][freq]
            w_pre   = np.sum(w  [:,:,:mid_time_i],2)    
            wx_pre  = np.sum(wx [:,:,:mid_time_i],2)
            wy_pre  = np.sum(wy [:,:,:mid_time_i],2)
            wxy_pre = np.sum(wxy[:,:,:mid_time_i],2)
            wxx_pre = np.sum(wxx[:,:,:mid_time_i],2)
            
            # sum over values i > mid_time_i
            w_pst   = np.sum(w  [:,:,-mid_time_i:],2)
            wx_pst  = np.sum(wx [:,:,-mid_time_i:],2)
            wy_pst  = np.sum(wy [:,:,-mid_time_i:],2)
            wxy_pst = np.sum(wxy[:,:,-mid_time_i:],2)
            wxx_pst = np.sum(wxx[:,:,-mid_time_i:],2)
            
            # calculate slopes and intercepts
            delta_pre = w_pre*wxx_pre - wx_pre**2
            delta_pst = w_pst*wxx_pst - wx_pst**2
            
            sl_pre = (w_pre*wxy_pre - wx_pre*wy_pre)/delta_pre
            sl_pst = (w_pst*wxy_pst - wx_pst*wy_pst)/delta_pst
            dsl_pre = np.sqrt(w_pre/delta_pre)
            dsl_pst = np.sqrt(w_pst/delta_pst)
            
            intr_pre = (wy_pre*wxx_pre - wx_pre*wxy_pre)/delta_pre
            intr_pst = (wy_pst*wxx_pst - wx_pst*wxy_pst)/delta_pst
            dintr_pre = np.sqrt(wxx_pre/delta_pre)
            dintr_pst = np.sqrt(wxx_pst/delta_pst)
            
            # extrapolate to middle time bin
            asym_slopes = intr_pst-intr_pre+(sl_pst-sl_pre)*mid_time
            asym_slopes_err = np.sqrt(dintr_pre**2 + dintr_pst**2 + \
                                      (dsl_pre**2 + dsl_pst**2) * mid_time**2)
            
            # save to output        
            out['sl_p'] = np.array([asym_slopes[0],asym_slopes_err[0]])
            out['sl_n'] = np.array([asym_slopes[1],asym_slopes_err[1]])
            out['sl_c'] = np.array([asym_slopes[2],asym_slopes_err[2]])
        
        # calculate asymmetry using differences
        asym_diff = asym[:,:,mid_time_i+1] - asym[:,:,mid_time_i-1]
        asym_diff_err = np.sqrt(asym_err[:,:,mid_time_i+1]**2+\
                                asym_err[:,:,mid_time_i-1]**2)
        
        # save to output
        out['dif_p'] = np.array([asym_diff[0],asym_diff_err[0]])
        out['dif_n'] = np.array([asym_diff[1],asym_diff_err[1]])
        out['dif_c'] = np.array([asym_diff[2],asym_diff_err[2]])   
            
        return out
676362
private void tred2() {
    // This is derived from the Algol procedures tred2 by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    System.arraycopy(V[n - 1], 0, d, 0, n);

    // Householder reduction to tridiagonal form.
    for(int i = n - 1; i > 0; i--) {
      double[] V_i = V[i], V_im1 = V[i - 1];
      // Scale to avoid under/overflow.
      double scale = 0.;
      for(int k = 0; k < i; k++) {
        scale += abs(d[k]);
      }
      if(scale < Double.MIN_NORMAL) {
        e[i] = d[i - 1];
        for(int j = 0; j < i; j++) {
          d[j] = V_im1[j];
          V_i[j] = V[j][i] = 0.;
        }
        d[i] = 0;
        continue;
      }
      // Generate Householder vector.
      double h = 0.;
      for(int k = 0; k < i; k++) {
        double dk = d[k] /= scale;
        h += dk * dk;
      }
      {
        double f = d[i - 1], g = FastMath.sqrt(h);
        g = (f > 0) ? -g : g;
        e[i] = scale * g;
        h -= f * g;
        d[i - 1] = f - g;
        Arrays.fill(e, 0, i, 0.);
      }

      // Apply similarity transformation to remaining columns.
      for(int j = 0; j < i; j++) {
        final double[] Vj = V[j];
        double dj = Vj[i] = d[j], ej = e[j] + Vj[j] * dj;
        for(int k = j + 1; k <= i - 1; k++) {
          final double Vkj = V[k][j];
          ej += Vkj * d[k];
          e[k] += Vkj * dj;
        }
        e[j] = ej;
      }
      double sum = 0.;
      for(int j = 0; j < i; j++) {
        sum += (e[j] /= h) * d[j];
      }
      double hh = sum / (h + h);
      for(int j = 0; j < i; j++) {
        e[j] -= hh * d[j];
      }
      for(int j = 0; j < i; j++) {
        double dj = d[j], ej = e[j];
        for(int k = j; k <= i - 1; k++) {
          V[k][j] -= (dj * e[k] + ej * d[k]);
        }
        d[j] = V_im1[j];
        V_i[j] = 0.;
      }
      d[i] = h;
    }

    // Accumulate transformations.
    tred2AccumulateTransformations();
    System.arraycopy(V[n - 1], 0, d, 0, n);
    Arrays.fill(V[n - 1], 0.);
    V[n - 1][n - 1] = 1.;
    e[0] = 0.;
  }
531223
function estimateBandwidth(array) {
  var n = array.length,
      q = quartiles(array),
      h = (q[2] - q[0]) / 1.34;
  return 1.06 * Math.min(Math.sqrt(variance(array)), h) * Math.pow(n, -0.2);
}
1096009
def F(self, x):
        
        if isinstance(x, np.ndarray):
            nfwvals = np.ones_like(x)
            inds1 = np.where(x < 1)
            inds2 = np.where(x > 1)
            nfwvals[inds1] = (1 - x[inds1] ** 2) ** -.5 * np.arctanh((1 - x[inds1] ** 2) ** .5)
            nfwvals[inds2] = (x[inds2] ** 2 - 1) ** -.5 * np.arctan((x[inds2] ** 2 - 1) ** .5)
            return nfwvals

        elif isinstance(x, float) or isinstance(x, int):
            if x == 1:
                return 1
            if x < 1:
                return (1 - x ** 2) ** -.5 * np.arctanh((1 - x ** 2) ** .5)
            else:
                return (x ** 2 - 1) ** -.5 * np.arctan((x ** 2 - 1) ** .5)
361911
public static void FFTShift1D(double[] x, Direction direction) {

        if (x.length == 1)
            return;

        double[] temp = x.clone();
        int move = x.length / 2;

        if (direction == Direction.Forward) {
            int c = 0;
            for (int i = x.length - move; i < x.length; i++)
                x[c++] = temp[i];

            for (int i = 0; i < x.length - move; i++)
                x[c++] = temp[i];

        } else {
            int c = 0;
            for (int i = move; i < x.length; i++)
                x[c++] = temp[i];

            for (int i = 0; i < move; i++)
                x[c++] = temp[i];

        }
    }
329155
def array_map(f, ar):
    ""
    flat_ar = ravel(ar)
    out = zeros(len(flat_ar), flat_ar.typecode())
    for i in range(len(flat_ar)):
        out[i] = f(flat_ar[i])
    out.shape = ar.shape
    return out
1330118
protected DoubleMatrix2D getHessF0(DoubleMatrix1D X) {
		double[][] hess = request.getF0().hessian(X.toArray());
		if(hess == FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER){
			return F2.make(X.size(), X.size());
		}else{
			return F2.make(hess);
		}
	}
170401
def asarray2d(a):
    """"""
    arr = np.asarray(a)
    if arr.ndim == 1:
        arr = arr.reshape(-1, 1)
    return arr
910897
def obfn_reg(self):
        

        rl1 = np.linalg.norm((self.Wl1 * self.obfn_g1var()).ravel(), 1)
        rtv = np.sum(np.sqrt(np.sum(self.obfn_g0var()**2, axis=-1)))
        return (self.lmbda*rl1 + self.mu*rtv, rl1, rtv)
510473
private double calcSD(double[] arr, double mean) {
        double sd = 0;
        for (double resid : arr) {
            sd += Math.pow(resid - mean, 2);
        }
        sd = Math.sqrt(sd/arr.length);
        return sd;
    }
38496
def obfn_dfd(self):
        
        """

        DXBf = sl.dot(self.B, sl.inner(self.Df, self.obfn_fvarf(),
                                       axis=self.cri.axisM),
                       axis=self.cri.axisC)
        Ef = DXBf - self.Sf
        return sl.rfl2norm2(Ef, self.S.shape, axis=self.cri.axisN) / 2.0
669111
def MeanVarEstEQ(y, x, covariates, tol=1e-8):
    

    pcount = covariates.shape[0] + 2
    N = y.shape[0]
    beta_count = pcount * 2
    X = numpy.ones((N, pcount))
    X[:, 1] = x

    aprime = [1, 0]
    idx = 2
    for cov in covariates:
        aprime.append(-numpy.mean(cov)/numpy.std(cov))
        X[:, idx] = cov
        idx += 1
    aprime = numpy.matrix(aprime)

    # http://stackoverflow.com/questions/7267226/range-for-floats
    def frange(x, y, jump):
        while x < y:
            yield x
            x += jump
        yield 0

    class PhiReturn(object):
        def __init__(self, phi, dtheta):
            self.phi = phi
            self.dtheta = dtheta

    def dot_diag(x, y):
        """This is a conveniance function to perform the following \
        in a more efficient manner:

        :param x: arr1
        :param y: arr2

        x.dot(numpy.diag(y))

        y must be a single dimensioned array

        """
        if len(y.shape) != 1:
            print >> sys.stderr, "You can't pass an array of shape %s to dot_diag" % (y.shape)
            sys.exit(1)

        result = numpy.empty(x.shape)
        for i in range(0, x.shape[0]):
            result[i] = x[i] * y

        return result

    def diag_dot(x, y):
        """diagonal multiply
        """
        result = numpy.empty(y.shape)
        for i in range(0, y.shape[0]):
            result[i] = y[i] * x[i]
        return result

    def Phi(theta):
        MM = y - numpy.dot(X, theta[0:pcount])
        SS = numpy.exp(numpy.dot(-X, theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        Phi = numpy.hstack((DD1.dot(X), (DD2-0.5).dot(X)))

        tX = X.transpose()

        t1 = numpy.empty(tX.shape)
        t2 = numpy.empty(tX.shape)
        t3 = numpy.empty(tX.shape)
        for i in range(0,tX.shape[0]):
            t1[i] = tX[i]*-SS
            t2[i] = tX[i]*-DD1
            t3[i] = tX[i]*-DD2
        t1 = t1.dot(X)
        t2 = t2.dot(X)
        t3 = t3.dot(X)

        dtheta = numpy.hstack((numpy.vstack((t1, t2)), numpy.vstack((t2, t3))))
        return PhiReturn(Phi, dtheta)



    class MvSolveReturn(object):
        def __init__(self, theta, dtheta):
            self.theta = theta
            self.dtheta = dtheta

    def MVsolve(theta_new):
        solution_found = False

        mvsolve_iterations = 0

        while not solution_found:
            theta_old = theta_new.copy()
            tmp = Phi(theta_old)
            #print "ITR", mvsolve_iterations, numpy.mean(theta_new),numpy.sum(numpy.absolute(theta_new - theta_old))
            theta_new = theta_old - scipy.linalg.solve(tmp.dtheta, tmp.phi)
            mvsolve_iterations += 1

            if (numpy.sum(numpy.absolute(theta_new - theta_old)) < tol):
                tmp = Phi(theta_new)
                solution_found = True
            else:
                if mvsolve_iterations > 25000:
                    #print >> sys.stderr, mvsolve_iterations, "failures"
                    raise UnsolvedLocus("")

        return MvSolveReturn(theta_new, tmp.dtheta), mvsolve_iterations

    def MVcalcB(theta):
        MM = y - X.dot(theta[0:pcount])
        SS = numpy.exp(-X.dot(theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        AA = numpy.hstack((diag_dot(DD1, X), diag_dot(DD2-0.5,X)))

        return numpy.transpose(AA).dot(AA)/ N


    mod = None
    itr = 0
    total_iterations = 0
    for i in frange(0.00, 1.0, 0.05):
        theta=numpy.empty((beta_count))
        theta[:] = i

        try:
            mod, iterations = MVsolve(theta)
            total_iterations += iterations
            if i > 0.05:
                print >> sys.stderr, "Completed: ", total_iterations, itr
            break
        except exceptions.ValueError as e:
            pass
        except numpy.linalg.linalg.LinAlgError as e:
            pass
        except Exception as inst:
            #print type(inst)
            pass
        itr += 1

    if not mod:
        raise UnsolvedLocus("")

    try:
        ainv = scipy.linalg.inv(mod.dtheta) * N
    except:
        raise ValueError("Singular Matrix Encountered")
    B = MVcalcB(mod.theta)
    V = ainv.dot(B).dot(ainv.transpose())



    # Focus on the two parameters of interest
    theta2 = numpy.array([mod.theta[1], mod.theta[pcount+1]])
    V2 = V[1:beta_count:pcount,1:beta_count:pcount]
    pvalt = 1 - scipy.stats.chi2.cdf(theta2.dot(scipy.linalg.inv(V2)).dot(theta2) * N, 2)


    ## From Chun's updated code:
    theta= mod.theta
    se = numpy.sqrt(numpy.diag(V)/N)
    pval = 2*scipy.stats.norm.cdf(-numpy.absolute(mod.theta/se))
    return pvalt, theta, pval, se, V/N
141532
def _fD(self, neg_pairs, A):
    
    diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
    return np.log(np.sum(np.sqrt(np.sum(np.dot(diff, A) * diff, axis=1))) + 1e-6)
1351328
def bb2hw(a:Collection[int])->np.ndarray:
    ""
    return np.array([a[1],a[0],a[3]-a[1],a[2]-a[0]])
1145532
def multidimensional_discretization(rho, sigma, N=3, method='rouwenhorst', m=2):
    

    # rho is assumed to be a scalar
    # sigma is a positive symmetric matrix
    # N number of points in each non-degenerate dimension
    # m : standard deviations to approximate

    import scipy.linalg
    from itertools import product

    d = sigma.shape[1]

    sigma = sigma.copy()

    zero_columns = np.where(sigma.sum(axis=0)==0)[0]
    for i in zero_columns:
        sigma[i,i] = 1

    L = scipy.linalg.cholesky(sigma)

    N = int(N)

    if method=='tauchen':
        [nodes_1d, probas_1d] = tauchen(N, 0, rho, 1, m=m)
    elif method=='rouwenhorst':
        [nodes_1d, probas_1d] = rouwenhorst(rho, 1, N)

    markov_nodes = np.array( list( product( *([nodes_1d]*d))) ).T
    markov_indices = np.array( list( product( *([range(N)]*d)  )  ) ).T

    markov_nodes = np.dot(L, markov_nodes)


    transition_matrix = 1
    for i in range(d):
        transition_matrix = np.kron(transition_matrix, probas_1d)

    markov_nodes = np.ascontiguousarray(markov_nodes.T)

    for i in zero_columns:
        markov_nodes[:,i] = 0

    return [markov_nodes, transition_matrix]
166895
def cwt_haar(x: np.ndarray, scale=10):
    
    next_power_two = 2 ** int(np.log2(len(x)))

    x = x[0:next_power_two]
    num_data = len(x)

    # get FFT of x (eq. (3) in paper)
    x_hat = np.fft.fft(x)

    # Get omega (eq. (5) in paper)
    f = (2.0 * np.pi / num_data)
    omega = f * np.concatenate((np.arange(0, num_data // 2), np.arange(num_data // 2, num_data) * -1))

    # get psi hat (eq. (6) in paper)
    psi_hat = np.sqrt(2.0 * np.pi * scale) * normalized_haar_wavelet(scale * omega, scale)

    # get W (eq. (4) in paper)
    W = np.fft.ifft(x_hat * psi_hat)

    return W[2 * scale:-2 * scale]