def points_to_basis_dists(self, points):
        assert is_mat(points)
        assert is_float(points)
        (N, D) = points.shape
        assert D == self.grid.get_dim()

        G = self.grid

        # Get indices
        cell_coords = G.points_to_cell_coords(points)

        # Get rel distances
        rel_dist = G.points_to_low_vertex_rel_distance(points, cell_coords)
        assert (N, D) == rel_dist.shape

        # Get the vertices
        vertices = self.grid.cell_coords_to_vertex_indices(cell_coords)
        assert (N, 2 ** D) == vertices.shape

        # Calculate multilinear interp weights from distances
        weights = np.empty((N, 2 ** D))
        for (i, bin_vertex) in enumerate(itertools.product([0, 1], repeat=D)):
            vert_mask = np.array(bin_vertex, dtype=bool)
            weights[:, i] = np.product(rel_dist[:, vert_mask], axis=1) * np.product(
                1.0 - rel_dist[:, ~vert_mask], axis=1
            )

        point_dist = self.convert_to_sparse_matrix(cell_coords, vertices, weights)
        return point_dist

    def get_dataset_slice(self, in_dataset_obj, dataset_info, in_data_idx, out_shape, inp_filename=""):
        """Copys dataset values from one dataset object to another, but only certain indexes along a
        specific dimension of the data"""

        # Determine how to extact data other than the splice dimension
        in_dataset_indexes = dataset_info.input_data_indexes(in_dataset_obj, in_data_idx)

        # Obtain selected data for copying into output dataset
        try:
            if len(in_dataset_indexes) == 1 and not isinstance(in_dataset_indexes[0], slice):
                in_data = in_dataset_obj[:][numpy.array(in_dataset_indexes[0])]
            else:
                in_data = in_dataset_obj[:][tuple(in_dataset_indexes)]
        except IOError as exc:
            raise IOError("Can not read dataset %s from file %s: %s" % (dataset_info.inp_name, inp_filename, exc))

        # Set sliced data into output dataset
        if numpy.product(in_data.shape) > numpy.product(out_shape):
            self.logger.warning("Dataset %s requires destructive resizing" % (dataset_info.out_name))
            self.logger.debug("At indexes %s resizing source data of shape %s to %s." % (in_data_idx, in_data.shape, out_shape)) 
            stored_data = numpy.resize(in_data, out_shape)
        else:
            stored_data = in_data.reshape(out_shape)

        return stored_data

def fftconvolve(in1, in2, in3=None, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    copied from scipy, but here used to try out inverse filter
    doesn't work or I can't get it to work
    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1+s2-1

    # Always use 2**n-sized FFT
    fsize = 2**np.ceil(np.log2(size))
    IN1 = fftn(in1,fsize)
    #IN1 *= fftn(in2,fsize)
    IN1 /= fftn(in2,fsize)  # use inverse filter
    # note the inverse is elementwise not matrix inverse
    # is this correct, NO  doesn't seem to work
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if product(s1,axis=0) > product(s2,axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret,osize)
    elif mode == "valid":
        return _centered(ret,abs(s2-s1)+1)

def maxprod(data, num):
    
    def diagdirection(largest):
        for i in xrange(-Len+num, Len-num+1):
            diag = data.diagonal(i)
            size = diag.size
            for j in xrange(size-num+1):
                dp = np.product(diag[j:j+num])
                if dp > largest:
                    largest = dp
        return largest
    
    largest = 0
    Len = data.shape[0]
    for i in xrange(Len):
        row = data[i,:]
        col = data[:,i]
        for j in xrange(Len-num+1):
            rp = np.product(row[j:j+num])
            cp = np.product(col[j:j+num])
            if rp > largest:
                largest = rp
            if cp > largest:
                largest = cp
    
    largest = diagdirection(largest)
    data = data[:,::-1]
    largest = diagdirection(largest)
    
    return largest

def fftconvolve(in1, in2, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1+s2-1
    IN1 = fftn(in1,size)
    IN1 *= fftn(in2,size)
    ret = ifftn(IN1)
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if product(s1,axis=0) > product(s2,axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret,osize)
    elif mode == "valid":
        return _centered(ret,abs(s2-s1)+1)

def weighted_variance(image, mask, binary_image):
    """Compute the log-transformed variance of foreground and background
    
    image - intensity image used for thresholding
    
    mask - mask of ignored pixels
    
    binary_image - binary image marking foreground and background
    """
    if not np.any(mask):
        return 0
    #
    # Clamp the dynamic range of the foreground
    #
    minval = np.max(image[mask])/256
    if minval == 0:
        return 0
    
    fg = np.log2(np.maximum(image[binary_image & mask], minval))
    bg = np.log2(np.maximum(image[(~ binary_image) & mask], minval))
    nfg = np.product(fg.shape)
    nbg = np.product(bg.shape)
    if nfg == 0:
        return np.var(bg)
    elif nbg == 0:
        return np.var(fg)
    else:
        return (np.var(fg) * nfg + np.var(bg)*nbg) / (nfg+nbg)

    def predict(self, context):
        if self.dictionary is None or self.parameters is None:
            print('Train before predict!')
            return
        context = context[-self.context_size:]
        input = []
        for word in context:
            if word in self.dictionary:
                input.append(self.dictionary[word])
            else:
                input.append(0)
        W_size = np.product(self.W_shape)
        U_size = np.product(self.U_shape)
        H_size = np.product(self.H_shape)
        split_indices = [W_size, W_size + U_size, W_size + U_size + H_size]
        W, U, H, C = np.split(self.parameters, split_indices)
        W = W.reshape(self.W_shape)
        U = U.reshape(self.U_shape)
        H = H.reshape(self.H_shape)
        C = C.reshape(self.C_shape)

        x = np.concatenate([C[input[i]] for i in range(self.context_size)])
        x = np.append(x, 1.)    # Append bias term
        x = x.reshape(-1, 1)
        y = W.dot(x) + U.dot(np.tanh(H.dot(x)))

        # You don't want to predict unknown words (index 0)
        prediction = np.argmax(y[1:]) + 1
        return self.reverse_dictionary[prediction]

def SceneTrain(NumSub, ImagPerClassToRead, DataPath, TrainImgs, ResizeAmt):
	# Read Images and Generate Image Features XNow
	for i in range(0,NumSub):
		for j in range(0, ImagPerClassToRead):
			InitAll.tic()
			#print DataPath+str(i+1)+'frame'+str(TrainImgs[i][j])+'.png'
			XNow = cv2.imread(DataPath+str(i+1)+'frame'+str(TrainImgs[i][j])+'.png',0)
			XNow = cv2.resize(XNow, ResizeAmt, interpolation = cv2.INTER_CUBIC)
			XNow = InitAll.ComputeGIST(XNow)
			#print("Sub " + str(i+1) + " Image " + str(j+1))
			if(i==0 and j==0):
				X = np.reshape(XNow, (1,np.product(XNow.shape)))
			else:
				X = np.vstack((X,np.reshape(XNow, (1,np.product(XNow.shape)))))
			InitAll.toc()
		print "Subject " + str(i+1) + " done...."

	# Now Generate Class Labels Y
	# Class labels start from 1 and not 0
	Y = [i for i in range(1,NumSub+1)]*ImagPerClassToRead
	Y = list(np.sort(Y))

	SVMModel = svm.SVC()
	SVMModel.fit(X, Y)
	# Saving the objects:
	with open('SceneTrainedSVMModel'+strftime("%Y-%m-%d %H:%M:%S", gmtime())+'.pickle', 'w') as f:
		pickle.dump([X, Y, SVMModel], f)
	return SVMModel

   def _MakeUnitCellPhasesForT_Restricted(self, ijkRow, iOpTs):
      # <r + T|op| s> = <r + T + S|op|s + S> for any unit-cell translation S,T
      # we know the op is only non-zero for a limited range of T, which we
      # have given via iOpTs.
      #
      # we have given the ijkRow (lhs T), and are now looking for the ijkCol
      # (right T) for which <r+Tr| op |r+Tc> is non-zero. That means that
      # Tr-Tc must lie within the Tc supplied in iOpTs.

      # FIXME: do this properly.
      if 0:
         I = []
         PF = []
         for ijkOp in (self.iTs[o] for o in iOpTs):
            for iTCol,ijkCol in enumerate(self.iTs):
               if np.all((ijkCol - ijkRow) % self.Size == ijkOp):
                  I.append(iTCol)
                  PF.append(np.product(self.PhaseShift**((ijkCol - ijkRow) / self.Size)))
         return np.array(PF), np.array(I)
      else:
         # well... that's not properly, but atm it's not the main problem in this form.
         I = []
         PF = []
         for ijkOp in (self.iTs[o] for o in iOpTs):
            ijkCol = ijkOp + ijkRow
            PF.append(np.product(self.PhaseShift**((ijkCol) / self.Size)))
            I.append(self.Fixme_iTsLookup[tuple(ijkCol % self.Size)])
         return np.array(PF), np.array(I)

    def get_srcp(self, img):
        import sourcecounts as sc
        fwsig = const.fwsig
        cutoff = 5.0
        spin = -0.80
        freq = img.frequency
        bm = (img.beam[0], img.beam[1])
        cdelt = img.wcs_obj.acdelt[:2]
        x = 2.0*pi*N.product(bm)/abs(N.product(cdelt))/(fwsig*fwsig)*img.omega

        smin_L = img.clipped_rms*cutoff*((1.4e9/freq)**spin)
        scflux = sc.s
        scnum = sc.n
        index = 0
        for i,s in enumerate(scflux):
            if s < smin_L:
                index = i
                break
        n1 = scnum[index]; n2 = scnum[-1]
        s1 = scflux[index]; s2 = scflux[-1]
        alpha = 1.0-log(n1/n2)/log(s1/s2)
        A = (alpha-1.0)*n1/(s1**(1.0-alpha))
        source_p = x*A*((cutoff*img.clipped_rms)**(1.0-alpha)) \
                     /((1.0-alpha)*(1.0-alpha))

        return source_p

    def check_meanmap(self, img, mean):
        """Calculates the statistics of the mean map and decides, when
        mean_map=None, whether to take the map (if variance
        is significant) or a constant value
        """
    	from math import sqrt

        mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Rmsimage.Checkmean ")
        cdelt = img.wcs_obj.acdelt[:2]
        bm = (img.beam[0], img.beam[1])
        fw_pix = sqrt(N.product(bm)/abs(N.product(cdelt)))
    	if img.masked:
            unmasked = N.where(~img.mask_arr)
            stdsub = N.std(mean[unmasked])
            maxmean = N.max(mean[unmasked])
        else:
            stdsub = N.std(mean)
            maxmean = N.max(mean)
        rms_expect = img.clipped_rms/img.rms_box[0]*fw_pix
        mylog.debug('%s %10.6f %s' % ('Standard deviation of mean image = ', stdsub*1000.0, 'mJy'))
        mylog.debug('%s %10.6f %s' % ('Expected standard deviation = ', rms_expect*1000.0, 'mJy'))

        # For mean map, use a higher threshold than for the rms map, as radio images
        # should rarely, if ever, have significant variations in the mean
        if stdsub > 5.0*rms_expect:
          img.mean_map_type = 'map'
          mylogger.userinfo(mylog, 'Variation in mean image significant')
        else:
          if img.confused:
            img.mean_map_type = 'zero'
          else:
            img.mean_map_type = 'const'
          mylogger.userinfo(mylog, 'Variation in mean image not significant')

        return img

    def check_rmsmap(self, img, rms):
        """Calculates the statistics of the rms map and decides, when
        rms_map=None, whether to take the map (if variance
        is significant) or a constant value
        """
    	from math import sqrt

        mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Rmsimage.Checkrms  ")
        cdelt = img.wcs_obj.acdelt[:2]
    	bm = (img.beam[0], img.beam[1])
    	fw_pix = sqrt(N.product(bm)/abs(N.product(cdelt)))
    	if img.masked:
    	    unmasked = N.where(~img.mask_arr)
            stdsub = N.std(rms[unmasked])
            maxrms = N.max(rms[unmasked])
        else:
            stdsub = N.std(rms)
            maxrms = N.max(rms)

    	rms_expect = img.clipped_rms/sqrt(2)/img.rms_box[0]*fw_pix
        mylog.debug('%s %10.6f %s' % ('Standard deviation of rms image = ', stdsub*1000.0, 'mJy'))
        mylog.debug('%s %10.6f %s' % ('Expected standard deviation = ', rms_expect*1000.0, 'mJy'))
    	if stdsub > 1.1*rms_expect:
            img.use_rms_map = True
            mylogger.userinfo(mylog, 'Variation in rms image significant')
        else:
            img.use_rms_map = False
            mylogger.userinfo(mylog, 'Variation in rms image not significant')

        return img

def df_fromdict(data, repeat=1):
    """
    Produces a factorial DataFrame from a dict or list of tuples.

    For example, suppose you want to generate a DataFrame like this::

           a    b
        0  one  0
        1  one  1
        2  two  0
        3  two  1

    This function generates such output simply by providing the following:
    df_fromdict([('a', ['one', 'two']), ('b', [0, 1])])

    :Args:
        data: dict or a list of tuples
            Data used to produce a DataFrame. Keys specify column names, and
            values specify possible (unique) values.
    :Kwargs:
        repeat: int (default: 1)
            How many times everything should be repeated. Useful if you want to
            simulate multiple samples of each condition, for example.
    :Returns:
        pandas.DataFrame with data.items() column names
    """
    data = OrderedDict(data)
    count = map(len, data.values())
    df = {}
    for i, (key, vals) in enumerate(data.items()):
        rep = np.repeat(vals, np.product(count[i+1:]))
        tile = np.tile(rep, np.product(count[:i]))
        df[key] = np.repeat(tile, repeat)
    df = pandas.DataFrame(df, columns=data.keys())
    return df

def fftconvolve(in1, in2, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1 + s2 - 1

    # Always use 2**n-sized FFT
    fsize = 2 ** np.ceil(np.log2(size))
    IN1 = fftn(in1, fsize)
    IN1 *= fftn(in2, fsize)
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if product(s1, axis=0) > product(s2, axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret, osize)
    elif mode == "valid":
        return _centered(ret, abs(s2 - s1) + 1)

def plot_hist_rmsRho_Levs012(rho_L0, rho_L1, rho_L2, tSim, fname):
    ''''''
    import pylab as py
    #make 6 hists
    rms = rho_L0
    lab = r"$\rho$, Lev 0"
    rms1d = np.reshape(rms, np.product(rms.shape))
    logrms= np.log10(rms1d)
    cnt,bins,patches = py.hist(logrms, bins=100, color="blue", alpha=0.5, label=lab)
    rms = rho_L1
    lab = r"$\rho$, Lev 1"
    rms1d = np.reshape(rms, np.product(rms.shape))
    logrms= np.log10(rms1d)
    cnt,bins,patches = py.hist(logrms, bins=100, color="red", alpha=0.5, label=lab)
    rms = rho_L2
    lab = r"$\rho$, Lev 2"
    rms1d = np.reshape(rms, np.product(rms.shape))
    logrms= np.log10(rms1d)
    #plot quantities
    Tratio = tSim / ic.tCr
    #plot
    cnt,bins,patches = py.hist(logrms, bins=100, color="green", alpha=0.5,label=lab)
    py.vlines(np.log10( ic.rho0 ), 0, cnt.max(), colors="black", linestyles='dashed',label=r"$\rho_{0}$ = %2.2g [g/cm^3]" % ic.rho0)
    #py.xlim([-13.,-9.])
    py.xlabel("Log10 Density [g/cm^3]")
    py.ylabel("count")
    py.title(r"$T/T_{\rmCross}$ = %g" % Tratio)
    py.legend(loc=0, fontsize="small")
    py.savefig(fname,format="pdf")
    py.close()

 def fit(self, n_components, data):
     if len(data.shape) == 2:
         self.reshape = None
     else:
         assert n_components == np.product(data.shape[1:]), \
             'ZCA whitening components should be %d for convolutional data'\
             % np.product(data.shape[1:])
         self.reshape = data.shape[1:]
     data = self._flatten_data(data)
     assert len(data.shape) == 2
     n, m = data.shape
     self.mean = np.mean(data, axis=0)
     bias_filter = self.filter_bias * np.identity(m, 'float64')
     cov = np.cov(data, rowvar=0, bias=1) + bias_filter
     eigs, eigv = np.linalg.eig(cov.astype(np.float64))
     assert not np.isnan(eigs).any()
     assert not np.isnan(eigv).any()
     print 'eigenvals larger than bias', np.sum(eigs > 0.1)/3072.
     print 'min eigenval: ', eigs.min(), 'max eigenval: ', eigs.max()
     assert eigs.min() > 0
     if self.n_components:
         eigs = eigs[-self.n_components:]
         eigv = eigv[:, -self.n_components:]
     sqrt_eigs = np.sqrt(eigs)
     self.P = np.dot(eigv * (1.0 / sqrt_eigs), eigv.T)
     assert not np.isnan(self.P).any()
     self.P_inv = np.dot(eigv * sqrt_eigs, eigv.T)
     self.P = np.float32(self.P)
     self.P_inv = np.float32(self.P_inv)
     self.is_fit = True

def weighted_variance(image,mask,threshold):
    """Compute the log-transformed variance of foreground and background"""
    if not np.any(mask):
        return 0
    #
    # Clamp the dynamic range of the foreground
    #
    minval = np.max(image[mask])/256
    if minval == 0:
        return 0
    clamped_image = image[mask]
    clamped_image[clamped_image < minval] = minval
    
    if isinstance(threshold,np.ndarray):
        threshold = threshold[mask]
    fg = np.log2(clamped_image[clamped_image >=threshold])
    bg = np.log2(clamped_image[clamped_image < threshold])
    nfg = np.product(fg.shape)
    nbg = np.product(bg.shape)
    if nfg == 0:
        return np.var(bg)
    elif nbg == 0:
        return np.var(fg)
    else:
        return (np.var(fg) * nfg + np.var(bg)*nbg) / (nfg+nbg)

def main():

    numberString = "08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08\
 49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00\
 81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65\
 52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91\
 22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80\
 24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50\
 32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70\
 67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21\
 24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72\
 21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95\
 78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92\
 16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57\
 86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58\
 19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40\
 04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66\
 88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69\
 04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36\
 20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16\
 20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54\
 01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48"

    array = np.array(map(int, numberString.split(' '))).reshape(20, 20)
    maxProduct = 0

    for i in xrange(array.shape[0]):
        for j in xrange(array.shape[1]):

            try:
                if np.product(array[i:i+4, j]) > maxProduct:
                    maxProduct = np.product(array[i:i+4, j])
            except IndexError:
                pass

            try:
                if np.product(array[i, j:j+4]) > maxProduct:
                    maxProduct = np.product(array[i, j:j+4])
            except IndexError:
                pass

            try:
                product = 1
                for k in range(4):
                    product *= array[i+k, j+k]
                if product > maxProduct:
                    maxProduct = product
            except IndexError:
                pass

            try:
                product = 1
                for k in range(4):
                    product *= array[i+k, j-k]
                if product > maxProduct:
                    maxProduct = product
            except IndexError:
                pass

    print maxProduct

    def modify(self, name, value):
        """ Change the value of an attribute while preserving its type.

        Differs from __setitem__ in that if the attribute already exists, its
        type is preserved.  This can be very useful for interacting with
        externally generated files.

        If the attribute doesn't exist, it will be automatically created.
        """
        with phil:
            if not name in self:
                self[name] = value
            else:
                value = numpy.asarray(value, order='C')

                attr = h5a.open(self._id, self._e(name))

                if attr.get_space().get_simple_extent_type() == h5s.NULL:
                    raise IOError("Empty attributes can't be modified")

                # Allow the case of () <-> (1,)
                if (value.shape != attr.shape) and not \
                   (numpy.product(value.shape) == 1 and numpy.product(attr.shape) == 1):
                    raise TypeError("Shape of data is incompatible with existing attribute")
                attr.write(value)

def apply_along_axis_with_idx(func1d,axis,arr,*args):
    """ Execute func1d(arr[i], i, *args) where func1d takes 1-D arrays
        and arr is an N-d array.  i varies so as to apply the function
        along the given axis for each 1-d subarray in arr.
    """
    arr = np.asarray(arr)
    nd = arr.ndim
    if axis < 0:
        axis += nd
    if (axis >= nd):
        raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
            % (axis,nd))
    ind = [0]*(nd-1)
    i = np.zeros(nd,'O')
    indlist = range(nd)
    indlist.remove(axis)
    i[axis] = slice(None,None)
    outshape = np.asarray(arr.shape).take(indlist)
    i.put(indlist, ind)
    res = func1d(arr[tuple(i.tolist())], tuple(i.tolist()), *args)
    #  if res is a number, then we have a smaller output array
    if isscalar(res):
        outarr = np.zeros(outshape,np.asarray(res).dtype)
        outarr[tuple(ind)] = res
        Ntot = np.product(outshape)
        k = 1
        while k < Ntot:
            # increment the index
            ind[-1] += 1
            n = -1
            while (ind[n] >= outshape[n]) and (n > (1-nd)):
                ind[n-1] += 1
                ind[n] = 0
                n -= 1
            i.put(indlist,ind)
            res = func1d(arr[tuple(i.tolist())], tuple(i.tolist()), *args)
            outarr[tuple(ind)] = res
            k += 1
        return outarr
    else:
        Ntot = np.product(outshape)
        holdshape = outshape
        outshape = list(arr.shape)
        outshape[axis] = len(res)
        outarr = np.zeros(outshape,np.asarray(res).dtype)
        outarr[tuple(i.tolist())] = res
        k = 1
        while k < Ntot:
            # increment the index
            ind[-1] += 1
            n = -1
            while (ind[n] >= holdshape[n]) and (n > (1-nd)):
                ind[n-1] += 1
                ind[n] = 0
                n -= 1
            i.put(indlist, ind)
            res = func1d(arr[tuple(i.tolist())], tuple(i.tolist()), *args)
            outarr[tuple(i.tolist())] = res
            k += 1
        return outarr