def calc_eff(ns0, nb0, ns1, nb1, alpha, sum_axes=None):

    if sum_axes:
        ns0 = np.apply_over_axes(np.sum, ns0, axes=sum_axes)
        nb0 = np.apply_over_axes(np.sum, nb0, axes=sum_axes)
        ns1 = np.apply_over_axes(np.sum, ns1, axes=sum_axes)
        nb1 = np.apply_over_axes(np.sum, nb1, axes=sum_axes)

    shape = np.broadcast(ns0, nb0, ns1, nb1).shape
    eff = np.zeros(shape)
    eff_var = np.zeros(shape)

    s0 = ns0 - alpha * nb0
    s1 = ns1 - alpha * nb1
    mask = (s0 * np.ones(shape) > 0)
    mask &= (s1 * np.ones(shape) > 0)

    s0[s0 <= 0] = 1.0
    eff = s1 / s0
    eff_var = (((ns0 - ns1 + alpha**2 * (nb0 - nb1)) * eff**2 +
                (ns1 + alpha**2 * nb1) * (1 - eff)**2) / s0**2)

    eff[~mask] = 0.0
    eff_var[~mask] = 0.0
    return eff, eff_var

def sort_and_avg(fr_array,sort_ind):
	#sort_ind 2D array, with col 0 = trial #, col 1 = sorting param

	unit_fr_mean = np.squeeze(np.apply_over_axes(np.mean,fr_array,(0,2)))
	unit_fr_std = np.squeeze(np.apply_over_axes(np.std,fr_array,(0,2)))

	zscore_array_all = []
	for i in range(int(np.max(sort_ind[:,1]) + 1)):
		temp = sort_ind[sort_ind[:,1] == i]
		for j in range(np.shape(temp)[0]):
			cond_trial_nums = temp[:,0]
			cond_trial_nums = cond_trial_nums.astype(int)

			fr_array_cond = fr_array[cond_trial_nums,:,:]
			mean_fr_cond = np.mean(fr_array_cond,axis=0)
			
			#gaussian_smoothed = gaussian_filter(mean_fr_cond,sigma=sigma_val)
			zscore_array = np.zeros((np.shape(mean_fr_cond)))
			for k in range(np.shape(mean_fr_cond)[0]):
				zscore_array[k,:] = (mean_fr_cond[k,:] - unit_fr_mean[k]) / unit_fr_std[k]

		zscore_array_all.append(zscore_array)

	zscore_array_all = np.asarray(zscore_array_all)
	dims = np.shape(zscore_array_all)
	zscore_array_all = zscore_array_all.reshape((dims[1],dims[0],dims[2]))



	return(zscore_array_all)

    def reject_parts(self):
        ok = self

        """
        Hall = self._parts * np.log(self._parts) + \
            (1 - self._parts) * np.log(1 - self._parts)
        H = -np.apply_over_axes(np.mean, Hall, [1, 2, 3]).ravel()
        """

        th1 = self._parts
        th0 = np.apply_over_axes(np.mean, self._parts, [1, 2])

        M1 = th1 * np.log(th1 / th0) +\
            (1 - th1) * np.log((1 - th1) / (1 - th0))
        S1 = np.log((th1 / (1 - th1) * ((1 - th0) / th0)))**2 * th1 * (1 - th1)
        mu1 = np.apply_over_axes(np.sum, M1, [1, 2, 3]).ravel()
        sigma1 = np.sqrt(np.apply_over_axes(np.sum, S1, [1, 2, 3])).ravel()

        M1 = th0 * np.log(th1 / th0) +\
            (1 - th1) * np.log((1 - th0) / (1 - th0))
        S1 = np.log((th1 / (1 - th1) * ((1 - th0) / th0)))**2 * th0 * (1 - th0)
        mu0 = np.apply_over_axes(np.sum, M1, [1, 2, 3]).ravel()
        # sigma0 = np.sqrt(np.apply_over_axes(np.sum, S1, [1, 2, 3])).ravel()

        ok = ((mu1 - mu0) / sigma1) > self._settings.get('reject_entropy', 1.0)

        print(ok.shape)
        print(ok.sum())

        self._parts = self._parts[ok]
        self._num_parts = self._parts.shape[0]

    def lnl_null(ns,nc,mub,alpha=None,data_axes=0,sum_lnl=True):
        """
        Log-likelihood for null hypothesis.

        Parameters
        ----------
        ns: Vector of observed counts in signal region.

        nc: Vector of observed counts in control region(s).
        """       
        lnls = poisson_lnl(ns,mub)
        lnlc = np.zeros(nc.shape)

        if alpha: 
            # model amplitude for counts in control region
            muc = np.apply_over_axes(np.sum,mub,data_axes)/alpha
            lnlc = poisson_lnl(nc,muc)

        if sum_lnl: 
            lnls = np.apply_over_axes(np.sum,lnls,data_axes)
            lnls = np.squeeze(lnls,data_axes)

            lnlc = np.apply_over_axes(np.sum,lnlc,data_axes)
            lnlc = np.squeeze(lnlc,data_axes)
            return lnls+lnlc
        else:
            return lnls

def pool(theta, size):
    w, h = theta.shape[0]//size, theta.shape[1]//size
    feat = np.zeros((w, h, theta.shape[-1]))
    for x, y in gv.multirange(w, h):
        feat[x,y] = np.apply_over_axes(np.sum, theta[x*size:(x+1)*size, y*size:(y+1)*size], [0, 1])[0,0]

    feat /= np.apply_over_axes(np.sum, feat, [-1]) + 1e-8
    
    return feat

 def verify_cumsum(h):
     for op in '<', '>':
         for kind in 'bincontent', 'binerror':
             func = lambda arr, axis: d.histfuncs.cumsum(arr, operator=op, axis=axis)
             if kind == 'bincontent':
                 cum = d.histfuncs.cumulative_bincontent(h, op)
                 cum_full = n.apply_over_axes(func, h._h_bincontent, range(h.ndim-1, -1, -1))[h._h_visiblerange]
             else:
                 cum = d.histfuncs.cumulative_binerror(h, op)
                 cum_full = n.sqrt(n.apply_over_axes(func, h._h_squaredweights, range(h.ndim-1, -1, -1))[h._h_visiblerange])
             assert((cum == cum_full).all())

def two_way(cells):
	dt = cells.dtype
	cells = cells.astype('float64')  # Make sure we don't overflow
	total = np.apply_over_axes(np.sum, cells, [1,2]).ravel()
	chi_sq = np.zeros(cells.shape, dtype='float64')
	for i in range(2):
		for j in range(2):
			exp = np.sum(cells[:,i,:], 1).ravel() * np.sum(cells[:,:,j], 1).ravel() / total
			chi_sq[:,i,j] = (cells[:,i,j] - exp)**2 / exp
	chi_sq = np.apply_over_axes(np.sum, chi_sq, [1,2]).ravel()
	return special.chdtrc(1, chi_sq).astype(dt)

def two_way(cells):
    """ Two-way chi-square test of independence. 
	Takes a 3D array as input: N(voxels) x 2 x 2, where the last two dimensions
	are the contingency table for each of N voxels. Returns an array of p-values.
	"""
    # dt = cells.dtype
    cells = cells.astype("float64")  # Make sure we don't overflow
    total = np.apply_over_axes(np.sum, cells, [1, 2]).ravel()
    chi_sq = np.zeros(cells.shape, dtype="float64")
    for i in range(2):
        for j in range(2):
            exp = np.sum(cells[:, i, :], 1).ravel() * np.sum(cells[:, :, j], 1).ravel() / total
            chi_sq[:, i, j] = (cells[:, i, j] - exp) ** 2 / exp
    chi_sq = np.apply_over_axes(np.sum, chi_sq, [1, 2]).ravel()
    return special.chdtrc(1, chi_sq)  # .astype(dt)

    def int_flux_threshold(self, skydir, fn, ts_thresh, min_counts):
        """Compute the integral flux threshold for a point source at
        position ``skydir`` with spectral parameterization ``fn``.

        """

        ebins = 10**np.linspace(np.log10(self.ebins[0]),
                                np.log10(self.ebins[-1]), 33)
        ectr = np.sqrt(ebins[0] * ebins[-1])

        sig, bkg, bkg_fit = self.compute_counts(skydir, fn, ebins)

        norms = irfs.compute_norm(sig, bkg, ts_thresh,
                                  min_counts, sum_axes=[1, 2, 3], bkg_fit=bkg_fit,
                                  rebin_axes=[4, 10, 1])

        npred = np.squeeze(np.apply_over_axes(np.sum, norms * sig, [1, 2, 3]))
        npred = np.array(npred, ndmin=1)
        flux = np.squeeze(norms) * fn.flux(ebins[0], ebins[-1])
        eflux = np.squeeze(norms) * fn.eflux(ebins[0], ebins[-1])
        dnde = np.squeeze(norms) * fn.dnde(ectr)
        e2dnde = ectr**2 * dnde

        o = dict(e_min=self.ebins[0], e_max=self.ebins[-1], e_ref=ectr,
                 npred=npred, flux=flux, eflux=eflux,
                 dnde=dnde, e2dnde=e2dnde)

        sig, bkg, bkg_fit = self.compute_counts(skydir, fn)

        npred = np.squeeze(np.apply_over_axes(np.sum, norms * sig,
                                              [2, 3]))
        flux = np.squeeze(np.squeeze(norms, axis=(1, 2, 3))[:, None] *
                          fn.flux(self.ebins[:-1], self.ebins[1:]))
        eflux = np.squeeze(np.squeeze(norms, axis=(1, 2, 3))[:, None] *
                           fn.eflux(self.ebins[:-1], self.ebins[1:]))
        dnde = np.squeeze(np.squeeze(norms, axis=(1, 2, 3))
                          [:, None] * fn.dnde(self.ectr))
        e2dnde = ectr**2 * dnde

        o['bins'] = dict(npred=npred,
                         flux=flux,
                         eflux=eflux,
                         dnde=dnde,
                         e2dnde=e2dnde,
                         e_min=self.ebins[:-1], e_max=self.ebins[1:],
                         e_ref=self.ectr)

        return o

    def _extract(self, phi, data):
        X = phi(data)
        XX = X[:, np.newaxis, np.newaxis]
        theta = self._models[np.newaxis]

        S = self._settings.get('standardize')
        if S:
            llh = XX * logit(theta)
            bb = np.apply_over_axes(np.sum, llh, [-3, -2, -1])[..., 0, 0, 0]
            bb = (bb - self._means) / self._sigmas
            yhat = np.argmax(bb.max(-1), axis=1)
        else:
            llh = XX * np.log(theta) + (1 - XX) * np.log(1 - theta)
            bb = np.apply_over_axes(np.sum, llh, [-3, -2, -1])[..., 0, 0, 0]
            yhat = np.argmax(bb.max(-1), axis=1)
        return yhat

    def to_wcs(self, sum_bands=False, normalize=True, proj='AIT', oversample=2,
               width_pix=None, hpx2wcs=None):

        from .wcsnd import WcsNDMap

        if sum_bands and self.geom.nside.size > 1:
            map_sum = self.sum_over_axes()
            return map_sum.to_wcs(sum_bands=False, normalize=normalize, proj=proj,
                                  oversample=oversample, width_pix=width_pix)

        # FIXME: Check whether the old mapping is still valid and reuse it
        if hpx2wcs is None:
            hpx2wcs = self.make_wcs_mapping(oversample=oversample, proj=proj,
                                            width_pix=width_pix)

        # FIXME: Need a function to extract a valid shape from npix property

        if sum_bands:
            hpx_data = np.apply_over_axes(np.sum, self.data,
                                          axes=np.arange(self.data.ndim - 1))
            hpx_data = np.squeeze(hpx_data)
            wcs_shape = tuple([t.flat[0] for t in hpx2wcs.npix])
            wcs_data = np.zeros(wcs_shape).T
            wcs = hpx2wcs.wcs.to_image()
        else:
            hpx_data = self.data
            wcs_shape = tuple([t.flat[0] for t in
                               hpx2wcs.npix]) + self.geom._shape
            wcs_data = np.zeros(wcs_shape).T
            wcs = hpx2wcs.wcs.to_cube(self.geom.axes)

        # FIXME: Should reimplement instantiating map first and fill data array
        hpx2wcs.fill_wcs_map_from_hpx_data(hpx_data, wcs_data, normalize)
        return WcsNDMap(wcs, wcs_data)

def downsample_rect(img, start_row, start_col, end_row, end_col, width, output, start_idx):
    """
    .. todo::

        WRITEME

    Parameters
    ----------
    img : WRITEME
        numpy matrix in topological order
        (batch size, rows, cols, channels)
    start_row : WRITEME
        row index of top-left corner of rectangle to average pool
    start_col : WRITEME
        col index of top-left corner of rectangle to average pool
    end_row : WRITEME
        row index of bottom-right corner of rectangle to average pool
    end_col : WRITEME
        col index of bottom-right corner of rectangle to average pool
    width : WRITEME
        take the mean over rectangular block of this width
    output : WRITEME
        dense design matrix, of shape (batch size, rows*cols*channels)
    start_idx : WRITEME
        column index where to start writing the output
    """
    idx = start_idx

    for i in xrange(start_row, end_row - width + 1, width):
        for j in xrange(start_col, end_col - width + 1, width):
            block = img[:, i:i+width, j:j+width]
            output[:,idx] = numpy.apply_over_axes(numpy.mean, block, axes=[1,2])[:,0,0]
            idx += 1

    return idx

def poisson_ul(nc,mus,mub,data_axes=1):
    """Test statistic for discovery with known background."""


    # MLE for signal norm under signal hypothesis
    snorm = np.apply_over_axes(np.sum,nc-mub,data_axes)
    snorm[snorm<0] = 0

    x = np.linspace(-3,3,50)

    mutot = snorm*mus+mub

    deltas = 10**x#*np.sum(mub)

    smutot = deltas[np.newaxis,np.newaxis,:]*mus[...,np.newaxis] + mutot[...,np.newaxis]

    lnl = nc[...,np.newaxis]*np.log(smutot) - smutot

    lnl = np.sum(lnl,axis=data_axes)

    ul = np.zeros(lnl.shape[0])

    for i in range(lnl.shape[0]):
        
        dlnl = -2*(lnl[i]-lnl[i][0])
        deltas_root = find_root(deltas,dlnl,2.72)
        ul[i] = snorm[i][0] + deltas_root

    return ul

def median_absolute_deviation(array, c=scipy.stats.norm.ppf(3/4.), axis=0,
                              center=np.median):
    """ The Median Absolute Deviation along given axis of an array.

    Parameters
    ----------
    array: array-like
        input array.
    c: float (optional, default scipy.stats.norm.ppf(3/4.) ~ .6745
        the normalization constant.
    axis: int (optional default 0)
        axes over which the callable fucntion `center` is applied.
    center: callable or float (default `np.median`)
        If a callable is provided then the array is centerd.
        Otherwise, a float represented the center is provided.

    Returns
    -------
    mad: float
        `mad` = median(abs(`array` - center)) / `c`
    """
    # Convert array-like object
    array = np.asarray(array)

    # Compute the center if a callable is passed in parameters
    if callable(center):
        center = np.apply_over_axes(center, array, axis)

    # Compute the median absolute deviation
    return np.median((np.fabs(array - center)) / c, axis=axis)

    def psf(self,emin,emax,cthmin,cthmax):
        """Return energy- and livetime-weighted PSF density vector as
        a function of angular offset for a bin in energy and
        inclination angle."""
        
        logemin = np.log10(emin)
        logemax = np.log10(emax)

        ilo = np.argwhere(self._energy > emin)[0,0]
        ihi = np.argwhere(self._energy < emax)[-1,0]+1
        
        jlo = np.argwhere(self._ctheta_axis.center > cthmin)[0,0]
        jhi = np.argwhere(self._ctheta_axis.center < cthmax)[-1,0] +1
        
        weights = (self._energy[ilo:ihi,np.newaxis]*
                   self._exp[ilo:ihi,jlo:jhi]*
                   self._wfn(self._energy[ilo:ihi,np.newaxis]))
        
        wsum = np.sum(weights)
        psf = np.apply_over_axes(np.sum,
                                 self._psf[:,ilo:ihi,jlo:jhi]*
                                 weights[np.newaxis,...],
                                 [1,2])
        psf = np.squeeze(psf)                   
        psf *= (1./wsum)
        return self._dtheta, psf

    def deviance(self, X, y, weights=None):
        '''Calculate the normal deviance (i.e. sum of squared errors) for
        every lambda.  The model must already be fit to call this method.
        '''
        self._check_if_fit()
        if weights is not None and weights.shape[0] != X.shape[0]:
            raise ValueError("The weights vector must have the same length as X.")

        # We normalise responses by default
        resp_weights = 1.0 / np.apply_along_axis(np.nanstd, 0, np.array(y))

        y_hat = self.predict(X)
        # Take the response y, and repeat it to produce a matrix
        # of the same dimensions as y_hat
        a = np.array(y)
        y_stacked = np.tile(a.reshape(a.shape + (1,)), (1, 1, y_hat.shape[-1]))
        rw_stacked = np.tile(resp_weights.reshape(1, len(resp_weights), 1), (y_hat.shape[0], 1, y_hat.shape[2]))
        if weights is None:
            sq_residuals = ((y_stacked - y_hat) * rw_stacked)**2
            normfac = X.shape[0] * y.shape[1]
        else:
            w = np.array(weights)
            w_stacked = np.tile(w.reshape((y_hat.shape[0], 1, 1)), (1,) + y_hat.shape[1:])
            sq_residuals = w_stacked * ((y_stacked - y_hat) * rw_stacked)**2
            normfac = np.sum(weights) * y.shape[1]
        return np.apply_over_axes(np.sum, sq_residuals, [0, 1]).ravel() / normfac

def cut_positions(filename, blurred, *positions):
	blurred = int(blurred)
	pos = eval("".join(positions))
	root = nc.open(filename)[0]
	lat = nc.getvar(root, 'lat')
	lon = nc.getvar(root, 'lon')
	data = nc.getvar(root, 'data')
	root_cut = nc.clonefile(root, 'cut_positions.' + filename, ['lat', 'lon', 'data'])[0]
	nc.getdim(root_cut, 'northing_cut', len(pos))
	nc.getdim(root_cut, 'easting_cut', 2)
	lat_cut = nc.getvar(root_cut, 'lat', 'f4', ('northing_cut','easting_cut',),4)
	lon_cut = nc.getvar(root_cut, 'lon', 'f4', ('northing_cut','easting_cut',),4)
	data_cut = nc.getvar(root_cut, 'data', 'f4', ('timing','northing_cut','easting_cut',),4)
	ix = 0
	for i in range(len(pos)):
		show("\rCutting data: processing position %d / %d " % (i+1, len(pos)))
		x, y = statistical_search_position(pos[i], lat, lon)
		if x and y:
			lat_cut[ix,0] = lat[x,y]
			lon_cut[ix,0] = lon[x,y]
			data_cut[:,ix,0] = np.apply_over_axes(np.mean, data[:,x-blurred:x+blurred,y-blurred:y+blurred], axes=[1,2]) if blurred > 0 else data[:,x,y]
			lat_cut[ix,1], lon_cut[ix,1], data_cut[:,ix,1] = lat_cut[ix,0], lon_cut[ix,0], data_cut[:,ix,0]
			ix += 1
	nc.close(root)
	nc.close(root_cut)

def genLocs(locs,predlocs,conf):
    dlocs = np.apply_over_axes(np.sum,(locs-predlocs)**2,axes=[1,2])
    dlocs = old_div(np.sqrt(dlocs),conf.n_classes)
    close = np.reshape(dlocs < (old_div(conf.gen_minlen,2)),[-1])
    newlocs = copy.deepcopy(predlocs)
    newlocs[close,...] = genFewMovedNegSamples(newlocs[close,...],conf,nmove=3)
    return newlocs

def _cumsum_with_overflow(bincontent, overflow, func):
    """
    Emulate the result of a cumulative sum over the entire histogram backing
    array given only disjoint views into the visible *bincontent* and the
    *overflow* (a slice in each dimension).
    """
    cum = bincontent
    ndim = bincontent.ndim
    # TODO: take axes to sum over as a parameter
    axes = range(ndim-1, -1, -1)
    for i, axis in enumerate(axes):
        # overflow should be a slab with one trivial dimension
        oflow = overflow[axis]
        assert(oflow.ndim == cum.ndim)
        assert(oflow.shape[axis] == 1)
    
        # sum over the visible part of the array
        cum = func(cum, axis=axis)
        # apply all the sums taken so far to the overflow slab, then add only
        # the part that is either in the overflow bin for this dimension, or
        # in a visible bin of another dimension
        idx = [slice(1,-1)]*ndim
        idx[axis] = slice(0,1)
        cum += n.apply_over_axes(func, oflow, axes[:i])[idx]
        
    return cum

def main():
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument('parts', metavar='<parts file>',
                        type=argparse.FileType('rb'),
                        help='Filename of parts file')
    args = parser.parse_args()

    feat_net = pnet.Layer.load(args.parts)

    layers = feat_net.layers

    S = 11

    layers += [
        pnet.PoolingLayer(shape=(S, S), strides=(S, S), operation='sum'),
        # pnet.MixtureClassificationLayer(n_components=10, min_prob=1e-5,
        # settings=dict( standardize=False,),)
        pnet.SVMClassificationLayer(C=100.0, settings=dict(standardize=True)),
    ]

    net = pnet.PartsNet(layers)

    limit = None
    error_rate, conf_mat = pnet.norb.train_and_test(net,
                                                    samples_per_class=None,
                                                    seed=0, limit=limit)

    print('Error rate: {:.02f}'.format(error_rate * 100))
    np.set_printoptions(precision=2, suppress=True)
    print('Confusion matrix:')

    norm_conf = conf_mat / np.apply_over_axes(np.sum, conf_mat, [1])
    print(norm_conf)