1 /***************************************************************************
2                           smoothPolyDnans.hpp -- smooth for n>2 Dims
3                              -------------------
4     begin                : May 30 2017
5     copyright            : (C) 2017 by G. Duvert
6     email                : m_schellens@users.sourceforge.net
7  ***************************************************************************/
8 
9 /***************************************************************************
10  *                                                                         *
11  *   This program is free software; you can redistribute it and/or modify  *
12  *   it under the terms of the GNU General Public License as published by  *
13  *   the Free Software Foundation; either version 2 of the License, or     *
14  *   (at your option) any later version.                                   *
15  *                                                                         *
16  ***************************************************************************/
17 
18 // to be included from datatypes.cpp
19 #ifdef INCLUDE_SMOOTH_POLYD
20 
21 SMOOTH_Ty* src=srcIn;
22 SMOOTH_Ty* dest=destIn;
23 SMOOTH_Ty* store;
24 
25 SizeT srcDim[MAXRANK];
26 for (int i = 0; i < rank; ++i) srcDim[i] = datainDim[i];
27 
28 SizeT nEl=1;
29 for (int i=0; i< rank; ++i) nEl*=srcDim[i];
30 
31 DUInt perm[rank]; //turntable transposition index list
32 for (int i = 0; i < rank; ++i) perm[i] = ((i + 1) % rank); //[1,2,...,0]
33 
34 SizeT destStride[ MAXRANK + 1];
35 
36 SizeT destDim[ MAXRANK]; //use destDim mainly to turn resDim each time
37 for (int i = 0; i < rank; ++i) destDim[i] = srcDim[i];
38 
39 // successively apply smooth 1d and write perm-transposed values in res, then exchange res and data for next iteration.
40 // the trick is to update srcDim and resStride each time.
41 for (int r = 0; r < rank; ++r) {
42  //accelerator: compute destStride from a perm-uted srcDim:
43  destStride[0] = 1;
44  destStride[1] = srcDim[perm[0]];
45  int m = 1;
46  for (; m < rank; ++m) destStride[m + 1] = destStride[m] * srcDim[perm[m]];
47  for (; m < MAXRANK; ++m) destStride[m + 1] = destStride[rank];
48 
49  SizeT dimx = srcDim[0]; //which has been permutated
50  SizeT dimy = nEl / dimx;
51  SizeT w = width[r] / 2;
52  if (w == 0) {//fast transpose
53 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
54   {
55 #pragma omp for nowait
56    for (SizeT i = 0; i < nEl; ++i) dest[transposed1Index(i, srcDim, destStride, rank)] = src[i];
57   }
58  } else { //smooth & transpose
59 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
60   {
61 #pragma omp for nowait
62 
63    for (SizeT j = 0; j < dimy; ++j) {
64     //initiate mean of Width first values:
65     DDouble z;
66     DDouble n = 0;
67     DDouble mean = 0;
68     for (SizeT i = j * dimx; i < j * dimx + (2 * w + 1); ++i) {
69      DDouble v = src[i];
70      n += 1.0;
71      z = 1. / n;
72      mean = (1. - z) * mean + z * v;
73     }
74 
75 #if defined(USE_EDGE)
76     //start: use mean1
77     DDouble mean1 = mean;
78     DDouble n1 = n;
79     for (SizeT i = 0; i < w; ++i) {
80      dest[transposed1Index(w - i + j * dimx, srcDim, destStride, rank)] = mean1;
81      DDouble v = src[2 * w - i + j * dimx];
82      mean1 -= z*v;
83 
84 #if defined (EDGE_WRAP)
85      v = src[dimx - i - 1 + j * dimx];
86 #elif defined (EDGE_TRUNCATE)
87      v = src[0 + j * dimx];
88 #elif defined (EDGE_MIRROR)
89      v = src[i + j * dimx];
90 #elif defined (EDGE_ZERO)
91      v = 0;
92 #endif
93       mean1 += z*v;
94     }
95     dest[transposed1Index(0 + j*dimx, srcDim, destStride, rank)] = mean1;
96 #else //just transpose
97     for (SizeT i = j * dimx; i < j * dimx + w; ++i) dest[transposed1Index(i, srcDim, destStride, rank)] = src[i];
98 #endif
99 
100     //center
101     for (SizeT i = w + j * dimx, ibef = j * dimx, iend = 2 * w + 1 + j * dimx; i < dimx - w - 1 + j * dimx; ++i, ++ibef, ++iend) {
102      dest[transposed1Index(i, srcDim, destStride, rank)] = mean;
103      DDouble v = src[ibef];
104      mean -= z*v;
105      v = src[iend];
106      mean += z*v;
107     }
108     dest[transposed1Index(dimx - 1 - w + j*dimx, srcDim, destStride, rank)] = mean;
109 
110 #if defined(USE_EDGE)
111     //end: use mean
112     for (SizeT i = dimx - 1 - w; i < dimx - 1; ++i) {
113      dest[transposed1Index(i + j*dimx, srcDim, destStride, rank)] = mean;
114      DDouble v = src[i - w + j * dimx];
115      mean -= z*v;
116 
117 #if defined (EDGE_WRAP)
118      v = src[i - dimx + w + 1 + j * dimx];
119 #elif defined (EDGE_TRUNCATE)
120      v = src[dimx - 1 + j * dimx];
121 #elif defined (EDGE_MIRROR)
122      v = src[2 * dimx - i - w - 2 + j * dimx];
123 #elif defined (EDGE_ZERO)
124      v = 0;
125 #endif
126       mean += z*v;
127     }
128     dest[transposed1Index(dimx - 1 + j*dimx, srcDim, destStride, rank)] = mean;
129 #else //just transpose
130     for (SizeT i = (dimx - w + j * dimx); i < (dimx + j * dimx); ++i) dest[transposed1Index(i, srcDim, destStride, rank)] = src[i];
131 #endif
132    }
133   }
134  }
135  //pseudo-dim of src is now rotated by 1
136  SizeT tempSize[MAXRANK];
137  for (int i = 0; i < rank; ++i) tempSize[i] = srcDim[i];
138  for (int i = 0; i < rank; ++i) srcDim[i] = tempSize[perm[i]];
139  // exchange dest and src
140   store=src;
141   src=dest;
142   dest=store; //now 'src' is the last computed array, whatever srcIn or destIn...
143 }
144 //if rank is even, must copy back src to destIn (since only the pointers were exchanged)
145 if (rank%2 == 0) memcpy(destIn, src, nEl * sizeof (src[0]));
146 #endif
147