def perform(self, node, inp, out):
        x, = inp
        z, = out
        if len(x.shape) != 4:
            raise NotImplementedError(
                'Pool requires 4D input for now')
        z_shape = self.out_shape(x.shape, self.ds, self.ignore_border, self.st,
                                 self.padding)
        if not self.ignore_border:
            assert z_shape[2] > 0
            assert z_shape[3] > 0
        if (z[0] is None) or (z[0].shape != z_shape):
            z[0] = numpy.empty(z_shape, dtype=x.dtype)
        zz = z[0]
        # number of pooling output rows
        pr = zz.shape[-2]
        # number of pooling output cols
        pc = zz.shape[-1]
        ds0, ds1 = self.ds
        st0, st1 = self.st
        pad_h = self.padding[0]
        pad_w = self.padding[1]
        img_rows = x.shape[-2] + 2 * pad_h
        img_cols = x.shape[-1] + 2 * pad_w
        inc_pad = self.mode == 'average_inc_pad'

        # pad the image
        if self.padding != (0, 0):
            y = numpy.zeros(
                (x.shape[0], x.shape[1], img_rows, img_cols),
                dtype=x.dtype)
            y[:, :, pad_h:(img_rows - pad_h), pad_w:(img_cols - pad_w)] = x
        else:
            y = x
        func = numpy.max
        if self.mode == 'sum':
            func = numpy.sum
        elif self.mode != 'max':
            func = numpy.average

        for n in xrange(x.shape[0]):
            for k in xrange(x.shape[1]):
                for r in xrange(pr):
                    row_st = r * st0
                    row_end = builtins.min(row_st + ds0, img_rows)
                    if not inc_pad:
                        row_st = builtins.max(row_st, self.padding[0])
                        row_end = builtins.min(row_end, x.shape[-2] + pad_h)
                    for c in xrange(pc):
                        col_st = c * st1
                        col_end = builtins.min(col_st + ds1, img_cols)
                        if not inc_pad:
                            col_st = builtins.max(col_st, self.padding[1])
                            col_end = builtins.min(col_end,
                                                   x.shape[-1] + pad_w)
                        zz[n, k, r, c] = func(y[
                            n, k, row_st:row_end, col_st:col_end])

    def perform(self, node, inp, out):
        if self.mode == 'average_exc_pad' and self.padding != (0, 0):
            raise NotImplementedError()
        x, gz = inp
        gx_stg, = out
        z_shape = self.out_shape(x.shape, self.ds, self.ignore_border, self.st,
                                 self.padding)
        if (gx_stg[0] is None) or (gx_stg[0].shape != z_shape):
            gx_stg[0] = numpy.empty(z_shape, dtype=x.dtype)
        zz = gx_stg[0]
        # number of pooling output rows
        pr = zz.shape[-2]
        # number of pooling output cols
        pc = zz.shape[-1]
        ds0, ds1 = self.ds
        st0, st1 = self.st
        pad_h = self.padding[0]
        pad_w = self.padding[1]
        img_rows = x.shape[-2] + 2 * pad_h
        img_cols = x.shape[-1] + 2 * pad_w
        inc_pad = self.mode == 'average_inc_pad'
        sum_mode = self.mode == 'sum'

        # pad the image
        if self.padding != (0, 0):
            y = numpy.zeros(
                (x.shape[0], x.shape[1], img_rows, img_cols),
                dtype=x.dtype)
            y[:, :, pad_h:(img_rows - pad_h), pad_w:(img_cols - pad_w)] = x
        else:
            y = x
        gx = numpy.zeros_like(y)
        for n in xrange(x.shape[0]):
            for k in xrange(x.shape[1]):
                for r in xrange(pr):
                    if sum_mode or inc_pad:
                        row_st = r * st0
                    else:
                        row_st = builtins.max(r * st0, self.padding[0])
                    row_end = builtins.min(row_st + ds0, img_rows)
                    for c in xrange(pc):
                        if sum_mode or inc_pad:
                            col_st = c * st1
                        else:
                            col_st = builtins.max(c * st1,
                                                  self.padding[1])
                        col_end = builtins.min(col_st + ds1, img_cols)
                        if sum_mode:
                            val = gz[n, k, r, c]
                        else:
                            val = gz[n, k, r, c] / ((row_end - row_st) *
                                                    (col_end - col_st))
                        gx[n, k, row_st:row_end, col_st:col_end] += val
        # unpad the image
        gx = gx[:, :, pad_h:(img_rows - pad_h), pad_w:(img_cols - pad_w)]
        gx_stg[0] = gx

    def perform(self, node, inp, out):
        x, = inp
        z, ind = out
        ind = numpy.zeros_like(x)
        if len(x.shape) != 4:
            raise NotImplementedError('Pool requires 4D input for now')
        z_shape = self.out_shape(x.shape, self.ds, self.ignore_border, self.st,
                                 self.padding)
        if (z[0] is None) or (z[0].shape != z_shape):
            z[0] = numpy.empty(z_shape, dtype=x.dtype)
        zz = z[0]
        # number of pooling output rows
        pr = zz.shape[-2]
        # number of pooling output cols
        pc = zz.shape[-1]
        ds0, ds1 = self.ds
        st0, st1 = self.st
        pad_h = self.padding[0]
        pad_w = self.padding[1]
        img_rows = x.shape[-2] + 2 * pad_h
        img_cols = x.shape[-1] + 2 * pad_w
        inc_pad = 0

        # pad the image
        if self.padding != (0, 0):
            y = numpy.zeros(
                (x.shape[0], x.shape[1], img_rows, img_cols),
                dtype=x.dtype)
            y[:, :, pad_h:(img_rows - pad_h), pad_w:(img_cols - pad_w)] = x
        else:
            y = x

        for n in xrange(x.shape[0]):
            for k in xrange(x.shape[1]):
                for r in xrange(pr):
                    row_st = r * st0
                    row_end = builtins.min(row_st + ds0, img_rows)
                    if not inc_pad:
                        row_st = builtins.max(row_st, self.padding[0])
                        row_end = builtins.min(row_end, x.shape[-2] + pad_h)
                    for c in xrange(pc):
                        col_st = c * st1
                        col_end = builtins.min(col_st + ds1, img_cols)
                        if not inc_pad:
                            col_st = builtins.max(col_st, self.padding[1])
                            col_end = builtins.min(col_end,
                                                   x.shape[-1] + pad_w)
                        cur_max = y[n, k, row_st, col_st]
                        max_r, max_c = row_st, col_st
                        for rr in xrange(row_st, row_end):
                            for cc in xrange(col_st, col_end):
                                if y[n, k, rr, cc] > cur_max:
                                    cur_max = y[n, k, rr, cc]
                                    max_r, max_c = rr, cc
                        zz[n, k, r, c] = cur_max
                        ind[n, k, max_r, max_c] = 1

    def numpy_max_pool_2d_stride_padding(
            x, ds, ignore_border=True, st=None, padding=(0, 0), mode='max'):
        assert ignore_border
        pad_h = padding[0]
        pad_w = padding[1]
        h = x.shape[-2]
        w = x.shape[-1]
        assert ds[0] > pad_h
        assert ds[1] > pad_w

        def pad_img(x):
            y = numpy.zeros(
                (x.shape[0], x.shape[1],
                 x.shape[2] + pad_h * 2, x.shape[3] + pad_w * 2),
                dtype=x.dtype)
            y[:, :, pad_h:(x.shape[2] + pad_h), pad_w:(x.shape[3] + pad_w)] = x

            return y
        img_rows = h + 2 * pad_h
        img_cols = w + 2 * pad_w
        out_r = (img_rows - ds[0]) // st[0] + 1
        out_c = (img_cols - ds[1]) // st[1] + 1
        out_shp = list(x.shape[:-2])
        out_shp.append(out_r)
        out_shp.append(out_c)
        ds0, ds1 = ds
        st0, st1 = st
        output_val = numpy.zeros(out_shp)
        y = pad_img(x)
        func = numpy.max
        if mode == 'sum':
            func = numpy.sum
        elif mode != 'max':
            func = numpy.average
        inc_pad = mode == 'average_inc_pad'

        for k in numpy.ndindex(*x.shape[:-2]):
            for i in range(output_val.shape[-2]):
                ii_st = i * st[0]
                ii_end = builtins.min(ii_st + ds[0], img_rows)
                if not inc_pad:
                    ii_st = builtins.max(ii_st, pad_h)
                    ii_end = builtins.min(ii_end, h + pad_h)
                for j in range(output_val.shape[-1]):
                    jj_st = j * st[1]
                    jj_end = builtins.min(jj_st + ds[1], img_cols)
                    if not inc_pad:
                        jj_st = builtins.max(jj_st, pad_w)
                        jj_end = builtins.min(jj_end, w + pad_w)
                    patch = y[k][ii_st:ii_end, jj_st:jj_end]
                    output_val[k][i, j] = func(patch)
        return output_val

    def restricted(self, min=None, max=None):
        """
        Returns the power series restricted to the coefficients starting at
        min and going up to, but not including max. If min is not
        specified, then it is assumed to be zero. If max is not specified,
        then it is assumed to be infinity.

        EXAMPLES::

            sage: L = LazyPowerSeriesRing(QQ)
            sage: a = L([1])
            sage: a.restricted().coefficients(10)
            [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
            sage: a.restricted(min=2).coefficients(10)
            [0, 0, 1, 1, 1, 1, 1, 1, 1, 1]
            sage: a.restricted(max=5).coefficients(10)
            [1, 1, 1, 1, 1, 0, 0, 0, 0, 0]
            sage: a.restricted(min=2, max=6).coefficients(10)
            [0, 0, 1, 1, 1, 1, 0, 0, 0, 0]
        """
        from six.moves import builtins
        if ((min is None and max is None) or
            (max is None and self.get_aorder() >= min)):
            return self

        return self._new(partial(self._restricted_gen, min, max),
                         lambda ao: builtins.max(ao, min), self)

def max(xs, y_from_x=lambda x: x):
    '''
    >>> range(10) | max()
    9
    >>> range(10) | max(lambda x: -x)
    0
    '''
    return __builtins__.max(y_from_x(x) for x in xs)

    def perform(self, node, inp, out):
        assert self.mode == 'max'
        x, maxout, gz, ws, stride, pad = inp
        gx_stg, = out
        assert ws.shape == stride.shape == pad.shape == (2,)
        # number of pooling output rows
        pr = maxout.shape[-2]
        # number of pooling output cols
        pc = maxout.shape[-1]
        ws0, ws1 = ws
        st0, st1 = stride
        pad_h = pad[0]
        pad_w = pad[1]
        img_rows = x.shape[-2] + 2 * pad_h
        img_cols = x.shape[-1] + 2 * pad_w

        # pad the image
        if (pad_h, pad_w) != (0, 0):
            y = numpy.zeros(
                (x.shape[0], x.shape[1], img_rows, img_cols),
                dtype=x.dtype)
            y[:, :, pad_h:(img_rows - pad_h), pad_w:(img_cols - pad_w)] = x
        else:
            y = x
        gx = numpy.zeros_like(y)
        for n in xrange(x.shape[0]):
            for k in xrange(x.shape[1]):
                for r in xrange(pr):
                    row_st = builtins.max(r * st0, pad_h)
                    row_end = builtins.min(row_st + ws0, img_rows)
                    for c in xrange(pc):
                        col_st = builtins.max(c * st1, pad_w)
                        col_end = builtins.min(col_st + ws1, img_cols)
                        for row_ind in xrange(row_st, row_end):
                            for col_ind in xrange(col_st, col_end):
                                if (maxout[n, k, r, c] == y[n, k, row_ind, col_ind]):
                                    gx[n, k, row_ind, col_ind] += gz[n, k, r, c]
        # unpad the image
        gx = gx[:, :, pad_h:(img_rows - pad_h), pad_w:(img_cols - pad_w)]
        gx_stg[0] = gx

    def numpy_max_pool_nd_stride_padding(
            input, ds, ignore_border=True, st=None, padding=None, mode='max'):
        assert ignore_border
        nd = len(ds)
        if padding is None:
            padding = (0,) * nd
        if st is None:
            st = (0,) * nd
        assert len(padding) == len(ds) == len(st)
        assert all(ds[i] > padding[i] for i in range(nd))

        def pad_img(x):
            # initialize padded input
            y = numpy.zeros(
                x.shape[0:-nd] +
                tuple(x.shape[-nd + i] + padding[i] * 2 for i in range(nd)),
                dtype=x.dtype)
            # place the unpadded input in the center
            block = ((slice(None),) * (len(x.shape) - nd) +
                     tuple(slice(padding[i], x.shape[-nd + i] + padding[i])
                           for i in range(nd)))
            y[block] = x
            return y

        pad_img_shp = list(input.shape[:-nd])
        out_shp = list(input.shape[:-nd])
        for i in range(nd):
            padded_size = input.shape[-nd + i] + 2 * padding[i]
            pad_img_shp.append(padded_size)
            out_shp.append((padded_size - ds[i]) // st[i] + 1)
        output_val = numpy.zeros(out_shp)
        padded_input = pad_img(input)
        func = numpy.max
        if mode == 'sum':
            func = numpy.sum
        elif mode != 'max':
            func = numpy.average
        inc_pad = mode == 'average_inc_pad'

        for l in numpy.ndindex(*input.shape[:-nd]):
            for r in numpy.ndindex(*output_val.shape[-nd:]):
                region = []
                for i in range(nd):
                    r_st = r[i] * st[i]
                    r_end = builtins.min(r_st + ds[i], pad_img_shp[-nd + i])
                    if not inc_pad:
                        r_st = builtins.max(r_st, padding[i])
                        r_end = builtins.min(r_end, input.shape[-nd + i] + padding[i])
                    region.append(slice(r_st, r_end))
                patch = padded_input[l][region]
                output_val[l][r] = func(patch)
        return output_val

    def perform(self, node, inp, out):
        assert self.mode == "max"
        x, maxout, gz = inp
        gx_stg, = out
        # number of pooling output rows
        pr = maxout.shape[-2]
        # number of pooling output cols
        pc = maxout.shape[-1]
        ds0, ds1 = self.ds
        st0, st1 = self.st
        pad_h = self.padding[0]
        pad_w = self.padding[1]
        img_rows = x.shape[-2] + 2 * pad_h
        img_cols = x.shape[-1] + 2 * pad_w

        # pad the image
        if self.padding != (0, 0):
            y = numpy.zeros((x.shape[0], x.shape[1], img_rows, img_cols), dtype=x.dtype)
            y[:, :, pad_h : (img_rows - pad_h), pad_w : (img_cols - pad_w)] = x
        else:
            y = x
        gx = numpy.zeros_like(y)
        for n in xrange(x.shape[0]):
            for k in xrange(x.shape[1]):
                for r in xrange(pr):
                    row_st = builtins.max(r * st0, self.padding[0])
                    row_end = builtins.min(row_st + ds0, img_rows)
                    for c in xrange(pc):
                        col_st = builtins.max(c * st1, self.padding[1])
                        col_end = builtins.min(col_st + ds1, img_cols)
                        for row_ind in xrange(row_st, row_end):
                            for col_ind in xrange(col_st, col_end):
                                if maxout[n, k, r, c] == y[n, k, row_ind, col_ind]:
                                    gx[n, k, row_ind, col_ind] += gz[n, k, r, c]
        # unpad the image
        gx = gx[:, :, pad_h : (img_rows - pad_h), pad_w : (img_cols - pad_w)]
        gx_stg[0] = gx

    def numpy_pool_2d_stride_padding(
            x, ds, ignore_border=True, st=None, padding=(0, 0), mode='max'):
        assert (ignore_border is False)

        in_h = x.shape[-2]
        in_w = x.shape[-1]
        kernel_h = ds[0]
        kernel_w = ds[1]
        stride_h = st[0]
        stride_w = st[1]
        pad_h = padding[0]
        pad_w = padding[1]

        assert ds[0] > pad_h
        assert ds[1] > pad_w

        def pad_img(x):
            y = numpy.zeros(
                (x.shape[0], x.shape[1],
                 x.shape[2] + pad_h * 2, x.shape[3] + pad_w * 2),
                dtype=x.dtype)
            y[:, :, pad_h:(x.shape[2] + pad_h), pad_w:(x.shape[3] + pad_w)] = x

            return y

        h = in_h + 2 * pad_h
        w = in_w + 2 * pad_w

        out_h = int(math.ceil((float)(h - kernel_h) / stride_h)) + 1
        out_w = int(math.ceil((float)(w - kernel_w) / stride_w)) + 1

        out_shp = list(x.shape[:-2])
        out_shp.extend([out_h, out_w])

        output_val = numpy.zeros(out_shp)

        y = pad_img(x)
        func = numpy.max
        if mode == 'sum':
            func = numpy.sum
        elif mode != 'max':
            func = numpy.average
        inc_pad = mode == 'average_inc_pad'

        for k in numpy.ndindex(*x.shape[:-2]):
            for i in range(output_val.shape[-2]):
                ii_st = i * st[0]
                if ii_st > h:
                    print ('ii_st > h!!!')
                    continue
                ii_end = builtins.min(ii_st + ds[0], h)
                if not inc_pad:
                    ii_st = builtins.max(ii_st, pad_h)
                    ii_end = builtins.min(ii_end, in_h + pad_h)
                for j in range(output_val.shape[-1]):
                    jj_st = j * st[1]
                    if jj_st > w:
                        print ('jj_st > w!!!')
                        continue
                    jj_end = builtins.min(jj_st + ds[1], w)
                    if not inc_pad:
                        jj_st = builtins.max(jj_st, pad_w)
                        jj_end = builtins.min(jj_end, in_w + pad_w)
                    patch = y[k][ii_st:ii_end, jj_st:jj_end]
                    output_val[k][i, j] = func(patch)
        return output_val

    def perform(self, node, inp, out):
        if self.mode not in ('max', 'sum') and self.padding != (0, 0):
            raise NotImplementedError()
        x, maxout, gz = inp
        gx_stg, = out
        # number of pooling output rows
        pr = maxout.shape[-2]
        # number of pooling output cols
        pc = maxout.shape[-1]
        ds0, ds1 = self.ds
        st0, st1 = self.st
        pad_h = self.padding[0]
        pad_w = self.padding[1]
        img_rows = x.shape[-2] + 2 * pad_h
        img_cols = x.shape[-1] + 2 * pad_w
        inc_pad = self.mode == 'average_inc_pad'
        sum_mode = self.mode == 'sum'

        # pad the image
        if self.padding != (0, 0):
            y = numpy.zeros(
                (x.shape[0], x.shape[1], img_rows, img_cols),
                dtype=x.dtype)
            y[:, :, pad_h:(img_rows-pad_h), pad_w:(img_cols-pad_w)] = x
        else:
            y = x
        gx = numpy.zeros_like(y)
        if self.mode == 'max':
            for n in xrange(x.shape[0]):
                for k in xrange(x.shape[1]):
                    for r in xrange(pr):
                        row_st = builtins.max(r * st0, self.padding[0])
                        row_end = builtins.min(row_st + ds0, img_rows)
                        for c in xrange(pc):
                            col_st = builtins.max(c * st1, self.padding[1])
                            col_end = builtins.min(col_st + ds1, img_cols)
                            for row_ind in xrange(row_st, row_end):
                                for col_ind in xrange(col_st, col_end):
                                    if (maxout[n, k, r, c] == y[n, k, row_ind, col_ind]):
                                        gx[n, k, row_ind, col_ind] += gz[n, k, r, c]
        else:
            for n in xrange(x.shape[0]):
                for k in xrange(x.shape[1]):
                    for r in xrange(pr):
                        if sum_mode or inc_pad:
                            row_st = r * st0
                        else:
                            row_st = builtins.max(r * st0, self.padding[0])
                        row_end = builtins.min(row_st + ds0, img_rows)
                        for c in xrange(pc):
                            if sum_mode or inc_pad:
                                col_st = c * st1
                            else:
                                col_st = builtins.max(c * st1,
                                                         self.padding[1])
                            col_end = builtins.min(col_st + ds1, img_cols)
                            if sum_mode:
                              val = gz[n, k, r, c]
                            else:
                              val = gz[n, k, r, c] / ((row_end - row_st) *
                                                      (col_end - col_st))
                            gx[n, k, row_st:row_end, col_st:col_end] += val
        # unpad the image
        gx = gx[:, :, pad_h:(img_rows-pad_h), pad_w:(img_cols-pad_w)]
        gx_stg[0] = gx