1"""
2TODO: implement Images2Neibs.infer_shape() methods
3
4"""
5from __future__ import absolute_import, print_function, division
6
7import numpy as np
8
9import theano
10from theano import Op, Apply
11from theano.gof import EnumList
12import theano.tensor as T
13from theano.gradient import grad_not_implemented
14from theano.gradient import grad_undefined
15
16
17class Images2Neibs(Op):
18    """
19    Reshapes the input as a 2D tensor where each row is an pooling
20    example.
21
22    Parameters
23    ----------
24    mode : {'valid', 'ignore_borders', 'wrap_centered'}
25        - 'valid' :
26            Requires an input that is a multiple of the pooling factor
27            (in each direction).
28        - 'half' :
29            Equivalent to 'valid' if we pre-pad with zeros the input on
30            each side by (neib_shape[0]//2, neib_shape[1]//2)
31        - 'full' :
32            Equivalent to 'valid' if we pre-pad with zeros the input on
33            each side by (neib_shape[0] - 1, neib_shape[1] - 1)
34        - 'ignore_borders' :
35            Same as valid, but will ignore the borders if the shape(s)
36            of the input is not a multiple of the pooling factor(s).
37        - 'wrap_centered' :
38            ?? TODO comment
39
40    """
41
42    __props__ = ("mode",)
43    BORDER_MODE = EnumList(('MODE_VALID', 'valid'),
44                           ('MODE_HALF', 'half'),
45                           ('MODE_FULL', 'full'),
46                           ('MODE_WRAP_CENTERED', 'wrap_centered'),
47                           ('MODE_IGNORE_BORDERS', 'ignore_borders'))
48    params_type = BORDER_MODE
49
50    def get_params(self, node):
51        return self.mode
52
53    def __init__(self, mode='valid'):
54        implemented_modes = self.BORDER_MODE.get_aliases()
55        if mode not in implemented_modes:
56            raise NotImplementedError("Only modes %s have been implemented for %s"
57                                      % (', '.join(implemented_modes), type(self).__name__))
58        self.mode = mode
59
60    def __str__(self):
61        return self.__class__.__name__ + "{%s}" % self.mode
62
63    def __setstate__(self, d):
64        self.__dict__.update(d)
65        if not hasattr(self, "mode"):
66            self.mode = 'valid'
67
68    def make_node(self, ten4, neib_shape, neib_step=None):
69        """
70        Parameters
71        ----------
72        ten4 : a list of lists of images
73            ten4 is of shape (list 1 dim, list 2 dim, row, col).
74        neib_shape
75            (r,c) where r is the height of the neighborhood in rows and c is
76            the width of the neighborhood in columns.
77        neib_step
78            (dr,dc) where dr is the number of rows to skip between patch and dc
79            is the number of columns. When None, this is the same as neib_shape
80            (patch are disjoint).
81
82        Returns
83        -------
84        matrix
85            A 2D matrix, written using the following pattern::
86
87                idx = 0
88                for i in xrange(list 1 dim)
89                    for j in xrange(list 2 dim)
90                        for k in <image column coordinates>
91                            for l in <image row coordinates>
92                                output[idx,:]
93                                     = flattened version of ten4[i,j,l:l+r,k:k+c]
94                                idx += 1
95
96            .. note:: The op isn't necessarily implemented internally with these
97                for loops, they're just the easiest way to describe the output
98                pattern.
99
100        """
101        ten4 = T.as_tensor_variable(ten4)
102        neib_shape = T.as_tensor_variable(neib_shape)
103        if neib_step is None:
104            neib_step = neib_shape
105        else:
106            neib_step = T.as_tensor_variable(neib_step)
107
108        assert ten4.ndim == 4
109        assert neib_shape.ndim == 1
110        assert neib_step.ndim == 1
111
112        return Apply(self, [ten4, neib_shape, neib_step],
113                     [T.matrix(dtype=ten4.type.dtype)])
114
115    def grad(self, inp, grads):
116        x, neib_shape, neib_step = inp
117        gz, = grads
118
119        if self.mode in ['valid', 'ignore_borders']:
120            if (neib_shape is neib_step or
121                neib_shape == neib_step or
122                # Theano Constant == do not compare the data
123                # the equals function do that.
124                (hasattr(neib_shape, "equals") and
125                 neib_shape.equals(neib_step))):
126                return [neibs2images(gz, neib_shape, x.shape, mode=self.mode),
127                        grad_undefined(self, 1, neib_shape),
128                        grad_undefined(self, 2, neib_step)]
129
130        if self.mode in ['valid']:
131            # Iterate over neighborhood positions, summing contributions.
132            def pos2map(pidx, pgz, prior_result, neib_shape, neib_step):
133                '''
134                Helper function that adds gradient contribution from a single
135                neighborhood position i,j.
136                pidx = Index of position within neighborhood.
137                pgz  = Gradient of shape (batch_size*num_channels*neibs)
138                prior_result  = Shape (batch_size, num_channnels, rows, cols)
139                neib_shape = Number of rows, cols in a neighborhood.
140                neib_step  = Step sizes from image2neibs.
141                '''
142                nrows, ncols = neib_shape
143                rstep, cstep = neib_step
144                batch_size, num_channels, rows, cols = prior_result.shape
145                i = pidx // ncols
146                j = pidx - (i * ncols)
147                # This position does not touch some img pixels in valid mode.
148                result_indices = prior_result[:, :,
149                                              i:(rows - nrows + i + 1):rstep,
150                                              j:(cols - ncols + j + 1):cstep]
151                newshape = (batch_size, num_channels) + \
152                           ((rows - nrows) // rstep + 1,) + \
153                           ((cols - ncols) // cstep + 1,)
154                return T.inc_subtensor(result_indices, pgz.reshape(newshape))
155            indices = T.arange(neib_shape[0] * neib_shape[1])
156            pgzs = gz.dimshuffle((1, 0))
157            result, _ = theano.scan(fn=pos2map,
158                                    sequences=[indices, pgzs],
159                                    outputs_info=T.zeros(x.shape),
160                                    non_sequences=[neib_shape, neib_step])
161            grad_input = result[-1]
162            return [grad_input,
163                    grad_undefined(self, 1, neib_shape),
164                    grad_undefined(self, 2, neib_step)]
165
166        return [grad_not_implemented(self, 0, x),
167                grad_undefined(self, 1, neib_shape),
168                grad_undefined(self, 2, neib_step)]
169
170    def c_code_cache_version(self):
171        return (10,)
172
173    def perform(self, node, inp, out_, params):
174        ten4, neib_shape, neib_step = inp
175        z, = out_
176        # GpuImages2Neibs should not run this perform in DebugMode
177        if type(self) != Images2Neibs:
178            raise theano.gof.utils.MethodNotDefined()
179
180        def CEIL_INTDIV(a, b):
181            if a % b:
182                return (a // b) + 1
183            else:
184                return a // b
185
186        grid_c = -1  # number of patch in height
187        grid_d = -1  # number of patch in width
188        assert ten4.ndim == 4
189        assert neib_shape.ndim == 1
190        assert neib_shape.shape[0] == 2
191        assert neib_step.ndim == 1
192        assert neib_step.shape[0] == 2
193        c, d = neib_shape
194        step_x, step_y = neib_step
195        mode = self.mode
196        if step_x <= 0 or step_y <= 0:
197            raise ValueError(
198                "neib_step wrong step ; values <= 0. Got " + str(neib_step))
199        if c <= 0 or d <= 0:
200            raise ValueError(
201                "neib_shape values <=0. Got " + str(neib_shape))
202
203        if mode == "wrap_centered":
204            if (c % 2 != 1) or (d % 2 != 1):
205                raise TypeError(
206                    "Images2Neibs:"
207                    " in mode wrap_centered need patch with odd shapes")
208
209            if (ten4.shape[2] < c) or (ten4.shape[3] < d):
210                raise TypeError(
211                    "Images2Neibs: in wrap_centered mode, don't support"
212                    " image shapes smaller then the patch shapes:"
213                    " neib_shape=(%d,%d), ten4[2:]=[%d,%d]" %
214                    (c, d, ten4.shape[2], ten4.shape[3]))
215            grid_c = CEIL_INTDIV(ten4.shape[2], step_x)
216            grid_d = CEIL_INTDIV(ten4.shape[3], step_y)
217        elif mode == "valid":
218            if (ten4.shape[2] < c) or (((ten4.shape[2] - c) % step_x) != 0):
219                raise TypeError(
220                    "neib_shape[0]=%d, neib_step[0]=%d and"
221                    " ten4.shape[2]=%d not consistent" %
222                    (c, step_x, ten4.shape[2]))
223            if (ten4.shape[3] < d) or (((ten4.shape[3] - d) % step_y) != 0):
224                raise TypeError(
225                    "neib_shape[1]=%d, neib_step[1]=%d and"
226                    " ten4.shape[3]=%d not consistent" %
227                    (d, step_y, ten4.shape[3]))
228            # number of patch in height
229            grid_c = 1 + ((ten4.shape[2] - c) // step_x)
230            # number of patch in width
231            grid_d = 1 + ((ten4.shape[3] - d) // step_y)
232        elif mode == "ignore_borders":
233            # number of patch in height
234            grid_c = 1 + ((ten4.shape[2] - c) // step_x)
235            # number of patch in width
236            grid_d = 1 + ((ten4.shape[3] - d) // step_y)
237        elif mode == "half":
238            # This is equivalent to 'valid' with padding (c // 2, d // 2) on both sides
239            # Thus the expanded image will have size (h + 2 * (c // 2), w + 2 * (d // 2))
240            # Plugging these in the equation for 'valid' we get
241            # h + 2 * (c // 2) - c  = h - (c % 2)
242            # w + 2 * (d // 2) - c  = w - (d % 2)
243            if (ten4.shape[2] < c) or (((ten4.shape[2] - (c % 2)) % step_x) != 0):
244                raise TypeError(
245                    "neib_shape[0]=%d, neib_step[0]=%d and"
246                    " ten4.shape[2]=%d not consistent" %
247                    (c, step_x, ten4.shape[2]))
248            if (ten4.shape[3] < d) or (((ten4.shape[3] - (d % 2)) % step_y) != 0):
249                raise TypeError(
250                    "neib_shape[1]=%d, neib_step[1]=%d and"
251                    " ten4.shape[3]=%d not consistent" %
252                    (d, step_y, ten4.shape[3]))
253            # number of patch in height
254            grid_c = 1 + ((ten4.shape[2] - (c % 2)) // step_x)
255            # number of patch in width
256            grid_d = 1 + ((ten4.shape[3] - (d % 2)) // step_y)
257        elif mode == "full":
258            # This is equivalent to 'valid' with padding (c - 1, d - 1) on both sides
259            # Thus the expanded image will have size (h + 2 * (c - 1), w + 2 * (d - 1))
260            # Plugging these in the equation for 'valid' we get
261            # h + 2 * (c - 1) - c  = h + c - 2
262            # w + 2 * (d - 1) - c  = w + d - 2
263            if (ten4.shape[2] < c) or (((ten4.shape[2] + c - 2) % step_x) != 0):
264                raise TypeError(
265                    "neib_shape[0]=%d, neib_step[0]=%d and"
266                    " ten4.shape[2]=%d not consistent" %
267                    (c, step_x, ten4.shape[2]))
268            if (ten4.shape[3] < d) or (((ten4.shape[3] + d - 2) % step_y) != 0):
269                raise TypeError(
270                    "neib_shape[1]=%d, neib_step[1]=%d and"
271                    " ten4.shape[3]=%d not consistent" %
272                    (d, step_y, ten4.shape[3]))
273            # number of patch in height
274            grid_c = 1 + ((ten4.shape[2] + c - 2) // step_x)
275            # number of patch in width
276            grid_d = 1 + ((ten4.shape[3] + d - 2) // step_y)
277        else:
278            raise TypeError("Images2Neibs: unknow mode '%s'" % mode)
279        z_dim0 = grid_c * grid_d * ten4.shape[1] * ten4.shape[0]
280        z_dim1 = c * d
281        z[0] = np.empty((z_dim0, z_dim1), dtype=node.outputs[0].dtype)
282
283        nb_batch = ten4.shape[0]
284        nb_stack = ten4.shape[1]
285        height = ten4.shape[2]
286        width = ten4.shape[3]
287
288        wrap_centered_half_idx_shift_x = c // 2
289        wrap_centered_half_idx_shift_y = d // 2
290        for n in range(nb_batch):
291            for s in range(nb_stack):
292                # loop over the number of patch in height
293                for a in range(grid_c):
294                    # loop over the number of patch in width
295                    for b in range(grid_d):
296                        z_row = b + grid_d * (a + grid_c * (s + nb_stack * n))
297                        for i in range(c):
298                            ten4_2 = i + a * step_x
299                            if mode == "wrap_centered":
300                                ten4_2 -= wrap_centered_half_idx_shift_x
301                                if ten4_2 < 0:
302                                    ten4_2 += height
303                                elif ten4_2 >= height:
304                                    ten4_2 -= height
305                            elif mode == "half":
306                                ten4_2 -= wrap_centered_half_idx_shift_x
307                            elif mode == "full":
308                                ten4_2 -= c - 1
309                            if ten4_2 < 0 or ten4_2 >= height:
310                                z[0][z_row, d * i: d * i + d] = 0
311                            else:
312                                for j in range(d):
313                                    ten4_3 = j + b * step_y
314                                    if mode == "wrap_centered":
315                                        ten4_3 -= wrap_centered_half_idx_shift_y
316                                        if ten4_3 < 0:
317                                            ten4_3 += width
318                                        elif ten4_3 >= width:
319                                            ten4_3 -= width
320                                    elif mode == "half":
321                                        ten4_3 -= wrap_centered_half_idx_shift_y
322                                    elif mode == "full":
323                                        ten4_3 -= d - 1
324                                    z_col = j + d * i
325                                    if ten4_3 < 0 or ten4_3 >= width:
326                                        z[0][z_row, z_col] = 0
327                                    else:
328                                        z[0][z_row, z_col] = ten4[n, s, ten4_2, ten4_3]
329
330    def infer_shape(self, node, input_shape):
331        in_shape = input_shape[0]
332        c, d = node.inputs[1]
333        step_x, step_y = node.inputs[2]
334        if self.mode == 'wrap_centered':
335            grid_c = T.ceil_intdiv(in_shape[2], step_x)
336            grid_d = T.ceil_intdiv(in_shape[3], step_y)
337        elif self.mode == 'valid':
338            grid_c = 1 + ((in_shape[2] - c) // step_x)
339            grid_d = 1 + ((in_shape[3] - d) // step_y)
340        elif self.mode == 'ignore_borders':
341            grid_c = 1 + ((in_shape[2] - c) // step_x)
342            grid_d = 1 + ((in_shape[3] - d) // step_y)
343        elif self.mode == 'half':
344            grid_c = 1 + ((in_shape[2] - (c % 2)) // step_x)
345            grid_d = 1 + ((in_shape[3] - (d % 2)) // step_y)
346        elif self.mode == 'full':
347            grid_c = 1 + ((in_shape[2] + c - 2) // step_x)
348            grid_d = 1 + ((in_shape[3] + d - 2) // step_y)
349        else:
350            raise TypeError("Images2Neibs: unknow mode '%s'" % self.mode)
351        z_dim0 = grid_c * grid_d * in_shape[1] * in_shape[0]
352        z_dim1 = c * d
353        return [(z_dim0, z_dim1)]
354
355    def c_code(self, node, name, inp, out, sub):
356        return """
357#ifndef CEIL_INTDIV
358#define CEIL_INTDIV(a, b) ((a/b) + ((a %% b) ? 1: 0))
359#endif
360
361        int grid_c = -1; //number of patch in height
362        int grid_d = -1; //number of patch in width
363        {
364        if (PyArray_NDIM(%(ten4)s) != 4)
365        {
366            PyErr_Format(PyExc_TypeError, "ten4 wrong rank");
367            %(fail)s;
368        }
369        if (PyArray_NDIM(%(neib_shape)s) != 1)
370        {
371            PyErr_Format(PyExc_TypeError, "neib_shape wrong rank");
372            %(fail)s;
373        }
374        if ( (PyArray_DIMS(%(neib_shape)s))[0] != 2)
375        {
376            PyErr_Format(PyExc_TypeError, "neib_shape wrong shape ; has to"
377                                          " contain 2 elements");
378            %(fail)s;
379        }
380        if (PyArray_NDIM(%(neib_step)s) != 1)
381        {
382            PyErr_Format(PyExc_TypeError, "neib_step wrong rank");
383            %(fail)s;
384        }
385        if ( (PyArray_DIMS(%(neib_step)s))[0] != 2)
386        {
387            PyErr_Format(PyExc_TypeError,
388                         "neib_step wrong step ; has to contain 2 elements");
389            %(fail)s;
390        }
391
392        // (c,d) = neib_shape
393        const npy_intp c = (npy_intp) *(dtype_%(neib_shape)s*) PyArray_GETPTR1(%(neib_shape)s, 0);
394        const npy_intp d = (npy_intp) *(dtype_%(neib_shape)s*) PyArray_GETPTR1(%(neib_shape)s, 1);
395        // (step_x,step_y) = neib_step
396        const dtype_%(neib_step)s step_x = *(dtype_%(neib_step)s*) PyArray_GETPTR1(%(neib_step)s, 0);
397        const dtype_%(neib_step)s step_y = *(dtype_%(neib_step)s*) PyArray_GETPTR1(%(neib_step)s, 1);
398
399        if (step_x <=0 || step_y <=0)
400        {
401            PyErr_Format(PyExc_ValueError,
402                         "neib_step wrong step ; values <= 0. Got %%lld %%lld.",
403                         (long long) step_x, (long long) step_y);
404            %(fail)s;
405        }
406
407        if (c <=0 || d <=0)
408        {
409            PyErr_Format(PyExc_ValueError,
410                         "neib_shape values <= 0. Got %%lld %%lld.",
411                         (long long)c, (long long)d);
412            %(fail)s;
413        }
414
415        if (%(mode)s == MODE_WRAP_CENTERED) {
416            if (c%%2!=1 || d%%2!=1){
417                PyErr_Format(PyExc_TypeError,
418                             "Images2Neibs: in mode wrap_centered"
419                             " need patch with odd shapes");
420                %(fail)s;
421            }
422            if ( (PyArray_DIMS(%(ten4)s))[2] < c ||
423                 (PyArray_DIMS(%(ten4)s))[3] < d)
424            {
425                PyErr_Format(PyExc_TypeError,
426                    "Images2Neibs: in wrap_centered mode, don't support image"
427                    " shapes smaller then the patch shapes:"
428                    " neib_shape=(%%ld,%%ld), ten4[2:]=[%%ld,%%ld]",
429                    (long int)c, (long int)d,
430                    (long int)(PyArray_DIMS(%(ten4)s)[2]),
431                    (long int)(PyArray_DIMS(%(ten4)s)[3]));
432                %(fail)s;
433            }
434            grid_c = CEIL_INTDIV(((PyArray_DIMS(%(ten4)s))[2]),step_x);
435            grid_d = CEIL_INTDIV(((PyArray_DIMS(%(ten4)s))[3]),step_y);
436
437        } else if (%(mode)s == MODE_VALID) {
438            if ( ((PyArray_DIMS(%(ten4)s))[2] < c) ||
439                 ( (((PyArray_DIMS(%(ten4)s))[2]-c) %% step_x)!=0))
440            {
441                PyErr_Format(PyExc_TypeError,
442                             "neib_shape[0]=%%ld, neib_step[0]=%%ld and"
443                             " ten4.shape[2]=%%ld not consistent",
444                             (long int)c, (long int)step_x,
445                             (long int)(PyArray_DIMS(%(ten4)s)[2]));
446                %(fail)s;
447            }
448            if ( ((PyArray_DIMS(%(ten4)s))[3] < d) ||
449                 ( (((PyArray_DIMS(%(ten4)s))[3]-d) %% step_y)!=0))
450            {
451                PyErr_Format(PyExc_TypeError,
452                             "neib_shape[1]=%%ld, neib_step[1]=%%ld and"
453                             " ten4.shape[3]=%%ld not consistent",
454                             (long int)d, (long int)step_y,
455                             (long int)(PyArray_DIMS(%(ten4)s)[3]));
456                %(fail)s;
457            }
458            //number of patch in height
459            grid_c = 1+(((PyArray_DIMS(%(ten4)s))[2]-c)/step_x);
460            //number of patch in width
461            grid_d = 1+(((PyArray_DIMS(%(ten4)s))[3]-d)/step_y);
462        } else if (%(mode)s == MODE_IGNORE_BORDERS) {
463            //number of patch in height
464            grid_c = 1+(((PyArray_DIMS(%(ten4)s))[2]-c)/step_x);
465            //number of patch in width
466            grid_d = 1+(((PyArray_DIMS(%(ten4)s))[3]-d)/step_y);
467        } else if (%(mode)s == MODE_HALF) {
468            if ( ((PyArray_DIMS(%(ten4)s))[2] < c) ||
469                 ( (((PyArray_DIMS(%(ten4)s))[2]-(c%%2)) %% step_x)!=0))
470            {
471                PyErr_Format(PyExc_TypeError,
472                             "neib_shape[0]=%%ld, neib_step[0]=%%ld and"
473                             " ten4.shape[2]=%%ld not consistent",
474                             (long int)c, (long int)step_x,
475                             (long int)(PyArray_DIMS(%(ten4)s)[2]));
476                %(fail)s;
477            }
478            if ( ((PyArray_DIMS(%(ten4)s))[3] < d) ||
479                 ( (((PyArray_DIMS(%(ten4)s))[3]-(d%%2)) %% step_y)!=0))
480            {
481                PyErr_Format(PyExc_TypeError,
482                             "neib_shape[1]=%%ld, neib_step[1]=%%ld and"
483                             " ten4.shape[3]=%%ld not consistent",
484                             (long int)d, (long int)step_y,
485                             (long int)(PyArray_DIMS(%(ten4)s)[3]));
486                %(fail)s;
487            }
488            //number of patch in height
489            grid_c = 1+(((PyArray_DIMS(%(ten4)s))[2]-(c%%2))/step_x);
490            //number of patch in width
491            grid_d = 1+(((PyArray_DIMS(%(ten4)s))[3]-(d%%2))/step_y);
492        } else if (%(mode)s == MODE_FULL) {
493            if ( ((PyArray_DIMS(%(ten4)s))[2] < c) ||
494                 ( (((PyArray_DIMS(%(ten4)s))[2]+c-2) %% step_x)!=0))
495            {
496                PyErr_Format(PyExc_TypeError,
497                             "neib_shape[0]=%%ld, neib_step[0]=%%ld and"
498                             " ten4.shape[2]=%%ld not consistent",
499                             (long int)c, (long int)step_x,
500                             (long int)(PyArray_DIMS(%(ten4)s)[2]));
501                %(fail)s;
502            }
503            if ( ((PyArray_DIMS(%(ten4)s))[3] < d) ||
504                 ( (((PyArray_DIMS(%(ten4)s))[3]+d-2) %% step_y)!=0))
505            {
506                PyErr_Format(PyExc_TypeError,
507                             "neib_shape[1]=%%ld, neib_step[1]=%%ld and"
508                             " ten4.shape[3]=%%ld not consistent",
509                             (long int)d, (long int)step_y,
510                             (long int)(PyArray_DIMS(%(ten4)s)[3]));
511                %(fail)s;
512            }
513            //number of patch in height
514            grid_c = 1+(((PyArray_DIMS(%(ten4)s))[2]+c-2)/step_x);
515            //number of patch in width
516            grid_d = 1+(((PyArray_DIMS(%(ten4)s))[3]+d-2)/step_y);
517        } else {
518            PyErr_Format(PyExc_TypeError,
519                         "Images2Neibs: unknow mode %%d", %(mode)s);
520            %(fail)s;
521        }
522
523        // new dimensions for z
524        const npy_intp z_dim1 = c * d;
525        const npy_intp z_dim0 =  grid_c
526                            * grid_d
527                            * (PyArray_DIMS(%(ten4)s))[1]
528                            * (PyArray_DIMS(%(ten4)s))[0];
529
530        if ((NULL == %(z)s)
531            || ((PyArray_DIMS(%(z)s))[0] != z_dim0 )
532            || ((PyArray_DIMS(%(z)s))[1] != z_dim1 )
533        )
534        {
535            Py_XDECREF(%(z)s);
536            npy_intp dims[2];
537            dims[0] = z_dim0;
538            dims[1] = z_dim1;
539
540            %(z)s = (PyArrayObject*) PyArray_EMPTY(2,
541                dims,
542                PyArray_TYPE((PyArrayObject*) py_%(ten4)s),
543                0);
544
545            if (!%(z)s)
546            {
547                PyErr_SetString(PyExc_MemoryError, "failed to alloc z output");
548                %(fail)s;
549            }
550        }
551        }
552
553        { // NESTED SCOPE
554
555        const int nb_batch = (PyArray_DIMS(%(ten4)s))[0];
556        const int nb_stack = (PyArray_DIMS(%(ten4)s))[1];
557        const int height = (PyArray_DIMS(%(ten4)s))[2];
558        const int width = (PyArray_DIMS(%(ten4)s))[3];
559
560        // (c,d) = neib_shape
561        const npy_intp c = (npy_intp) *(dtype_%(neib_shape)s*) PyArray_GETPTR1(%(neib_shape)s, 0);
562        const npy_intp d = (npy_intp) *(dtype_%(neib_shape)s*) PyArray_GETPTR1(%(neib_shape)s, 1);
563        // (step_x,step_y) = neib_step
564        const npy_intp step_x = (npy_intp) *(dtype_%(neib_step)s*) PyArray_GETPTR1(%(neib_step)s, 0);
565        const npy_intp step_y = (npy_intp) *(dtype_%(neib_step)s*) PyArray_GETPTR1(%(neib_step)s, 1);
566
567        const int wrap_centered_half_idx_shift_x = c/2;
568        const int wrap_centered_half_idx_shift_y = d/2;
569        // Oh this is messed up...
570        for (int n = 0; n < nb_batch; n++)              // loop over batches
571            for (int s = 0; s < nb_stack; s++)          // loop over stacks
572                for (int a = 0; a < grid_c; a++)        // loop over the number of patch in height
573                    for (int b = 0; b < grid_d; b++)    // loop over the number of patch in width
574                    {
575                        int z_row = b + grid_d*(a + grid_c*(s + nb_stack*n));
576                        for (int i = 0; i < c; i++)     // loop over c
577                        {
578                            int ten4_2 = i + a * step_x;
579                            if (%(mode)s == MODE_WRAP_CENTERED) {
580                                ten4_2 -= wrap_centered_half_idx_shift_x;
581                                if ( ten4_2 < 0 ) ten4_2 += height;
582                                else if (ten4_2 >= height) ten4_2 -= height;
583                            } else if (%(mode)s == MODE_HALF) {
584                                ten4_2 -= wrap_centered_half_idx_shift_x;
585                            } else if (%(mode)s == MODE_FULL) {
586                                ten4_2 -= c - 1;
587                            }
588                            if (ten4_2 < 0 | ten4_2 >= height) {
589                                dtype_%(z)s* curr_z = (dtype_%(z)s*) PyArray_GETPTR2(%(z)s, z_row, d * i);
590                                memset(curr_z, 0, d*sizeof(*curr_z));
591                            } else {
592                                for (int j = 0; j < d; j++)  // loop over d
593                                {
594                                    int ten4_3 = j + b * step_y;
595                                    if (%(mode)s == MODE_WRAP_CENTERED) {
596                                        ten4_3 -= wrap_centered_half_idx_shift_y;
597                                        if ( ten4_3 < 0 ) ten4_3 += width;
598                                        else if (ten4_3 >= width) ten4_3 -= width;
599                                    } else if (%(mode)s == MODE_HALF) {
600                                        ten4_3 -= wrap_centered_half_idx_shift_y;
601                                    } else if (%(mode)s == MODE_FULL) {
602                                        ten4_3 -= d - 1;
603                                    }
604                                    int z_col = j + d * i;
605                                    dtype_%(z)s* curr_z = (dtype_%(z)s*) PyArray_GETPTR2(%(z)s, z_row, z_col);
606                                    if (ten4_3 < 0 | ten4_3 >= width) {
607                                        *curr_z = 0;
608                                    } else {
609                                        *curr_z = *( (dtype_%(ten4)s*) PyArray_GETPTR4(%(ten4)s, n, s, ten4_2, ten4_3));
610                                    }
611                                }
612                            }
613                        }
614                    }
615        } // END NESTED SCOPE
616        """ % dict(ten4=inp[0], neib_shape=inp[1], neib_step=inp[2], z=out[0],
617                   fail=sub['fail'], mode=sub['params'])
618
619
620def images2neibs(ten4, neib_shape, neib_step=None, mode='valid'):
621    """
622    Function :func:`images2neibs <theano.tensor.nnet.neighbours.images2neibs>`
623    allows to apply a sliding window operation to a tensor containing
624    images or other two-dimensional objects.
625    The sliding window operation loops over points in input data and stores
626    a rectangular neighbourhood of each point.
627    It is possible to assign a step of selecting patches (parameter `neib_step`).
628
629    Parameters
630    ----------
631    ten4 : A 4d tensor-like
632        A 4-dimensional tensor which represents a list of lists of images.
633        It should have shape (list 1 dim, list 2 dim, row, col). The first
634        two dimensions can be useful to store different channels and batches.
635    neib_shape : A 1d tensor-like of 2 values
636        A tuple containing two values: height and width of the neighbourhood.
637        It should have shape (r,c) where r is the height of the neighborhood
638        in rows and c is the width of the neighborhood in columns.
639    neib_step : A 1d tensor-like of 2 values
640        (dr,dc) where dr is the number of rows to skip between patch and dc is
641        the number of columns. The parameter should be a tuple of two elements:
642        number of rows and number of columns to skip each iteration.
643        Basically, when the step is 1, the neighbourhood of every first element
644        is taken and every possible rectangular subset is returned.
645        By default it is equal to `neib_shape` in other words, the patches are
646        disjoint. When the step is greater than `neib_shape`, some elements are
647        omitted. When None, this is the same as neib_shape (patch are disjoint).
648    mode : {'valid', 'ignore_borders', 'wrap_centered', 'half'}
649        ``valid``
650            Requires an input that is a multiple of the
651            pooling factor (in each direction).
652        ``half``
653            Equivalent to 'valid' if we pre-pad with zeros the input on
654            each side by (neib_shape[0]//2, neib_shape[1]//2)
655        ``full``
656            Equivalent to 'valid' if we pre-pad with zeros the input on
657            each side by (neib_shape[0] - 1, neib_shape[1] - 1)
658        ``ignore_borders``
659            Same as valid, but will ignore the borders if the shape(s) of
660            the input is not a multiple of the pooling factor(s).
661        ``wrap_centered``
662            ?? TODO comment
663
664    Returns
665    -------
666    object
667        Reshapes the input as a 2D tensor where each row is an
668        pooling example. Pseudo-code of the output:
669
670          .. code-block:: python
671
672             idx = 0
673             for i in xrange(list 1 dim):
674                 for j in xrange(list 2 dim):
675                     for k in <image column coordinates>:
676                         for l in <image row coordinates>:
677                             output[idx,:]
678                                  = flattened version of ten4[i,j,l:l+r,k:k+c]
679                             idx += 1
680
681          .. note:: The operation isn't necessarily implemented internally with
682             these for loops, they're just the easiest way to describe the
683             output pattern.
684
685    Notes
686    -----
687    .. note::
688        Currently the step size should be chosen in the way that the
689        corresponding dimension :math:`i` (width or height) is equal
690        to :math:`n * step\_size_i + neib\_shape_i` for some :math:`n`.
691
692    Examples
693    --------
694
695    .. code-block:: python
696
697        # Defining variables
698        images = T.tensor4('images')
699        neibs = images2neibs(images, neib_shape=(5, 5))
700
701        # Constructing theano function
702        window_function = theano.function([images], neibs)
703
704        # Input tensor (one image 10x10)
705        im_val = np.arange(100.).reshape((1, 1, 10, 10))
706
707        # Function application
708        neibs_val = window_function(im_val)
709
710    .. note:: The underlying code will construct a 2D tensor of disjoint
711       patches 5x5. The output has shape 4x25.
712
713    """
714    return Images2Neibs(mode)(ten4, neib_shape, neib_step)
715
716
717def neibs2images(neibs, neib_shape, original_shape, mode='valid'):
718    """
719    Function :func:`neibs2images <theano.sandbox.neighbours.neibs2images>`
720    performs the inverse operation of
721    :func:`images2neibs <theano.sandbox.neigbours.neibs2images>`. It inputs
722    the output of :func:`images2neibs <theano.sandbox.neigbours.neibs2images>`
723    and reconstructs its input.
724
725    Parameters
726    ----------
727    neibs : 2d tensor
728        Like the one obtained by
729        :func:`images2neibs <theano.sandbox.neigbours.neibs2images>`.
730    neib_shape
731        `neib_shape` that was used in
732        :func:`images2neibs <theano.sandbox.neigbours.neibs2images>`.
733    original_shape
734        Original shape of the 4d tensor given to
735        :func:`images2neibs <theano.sandbox.neigbours.neibs2images>`
736
737    Returns
738    -------
739    object
740        Reconstructs the input of
741        :func:`images2neibs <theano.sandbox.neigbours.neibs2images>`,
742        a 4d tensor of shape `original_shape`.
743
744    Notes
745    -----
746    Currently, the function doesn't support tensors created with
747    `neib_step` different from default value. This means that it may be
748    impossible to compute the gradient of a variable gained by
749    :func:`images2neibs <theano.sandbox.neigbours.neibs2images>` w.r.t.
750    its inputs in this case, because it uses
751    :func:`images2neibs <theano.sandbox.neigbours.neibs2images>` for
752    gradient computation.
753
754    Examples
755    --------
756    Example, which uses a tensor gained in example for
757    :func:`images2neibs <theano.sandbox.neigbours.neibs2images>`:
758
759    .. code-block:: python
760
761        im_new = neibs2images(neibs, (5, 5), im_val.shape)
762        # Theano function definition
763        inv_window = theano.function([neibs], im_new)
764        # Function application
765        im_new_val = inv_window(neibs_val)
766
767    .. note:: The code will output the initial image array.
768
769    """
770    neibs = T.as_tensor_variable(neibs)
771    neib_shape = T.as_tensor_variable(neib_shape)
772    original_shape = T.as_tensor_variable(original_shape)
773
774    new_neib_shape = T.stack([original_shape[-1] // neib_shape[1],
775                              neib_shape[1]])
776    output_2d = images2neibs(neibs.dimshuffle('x', 'x', 0, 1),
777                             new_neib_shape, mode=mode)
778
779    if mode == 'ignore_borders':
780        # We use set_subtensor to accept original_shape we can't infer
781        # the shape and still raise error when it don't have the right
782        # shape.
783        valid_shape = original_shape
784        valid_shape = T.set_subtensor(
785            valid_shape[2],
786            (valid_shape[2] // neib_shape[0]) * neib_shape[0])
787        valid_shape = T.set_subtensor(
788            valid_shape[3],
789            (valid_shape[3] // neib_shape[1]) * neib_shape[1])
790        output_4d = output_2d.reshape(valid_shape, ndim=4)
791        # padding the borders with zeros
792        for d in [2, 3]:
793            pad_shape = list(output_4d.shape)
794            pad_shape[d] = original_shape[d] - valid_shape[d]
795            output_4d = T.concatenate([output_4d, T.zeros(pad_shape)], axis=d)
796    elif mode == 'valid':
797        # TODO: we do not implement all mode with this code.
798        # Add a check for the good cases.
799        output_4d = output_2d.reshape(original_shape, ndim=4)
800    else:
801        raise NotImplementedError("neibs2images do not support mode=%s" % mode)
802
803    return output_4d
804