source: thirdparty/SZ/sz/src/szd_double.c @ 9ee2ce3

Revision 9ee2ce3, 153.6 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)

importing new SZ files

  • Property mode set to 100644
RevLine 
[2c47b73]1/**
2 *  @file szd_double.c
3 *  @author Sheng Di and Dingwen Tao
4 *  @date Aug, 2016
5 *  @brief
6 *  (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7 *      See COPYRIGHT in top-level directory.
8 */
9
10#include <stdlib.h>
11#include <stdio.h>
12#include <string.h>
13#include "szd_double.h"
14#include "TightDataPointStorageD.h"
15#include "sz.h"
16#include "Huffman.h"
17#include "szd_double_pwr.h"
18#include "szd_double_ts.h"
[9ee2ce3]19#include "utility.h"
[2c47b73]20
21int SZ_decompress_args_double(double** newData, size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, unsigned char* cmpBytes, size_t cmpSize)
22{
23        int status = SZ_SCES;
24        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
25       
26        //unsigned char* tmpBytes;
27        size_t targetUncompressSize = dataLength <<3; //i.e., *8
28        //tmpSize must be "much" smaller than dataLength
29        size_t i, tmpSize = 12+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
30        unsigned char* szTmpBytes;
31        if(cmpSize!=12+4+MetaDataByteLength && cmpSize!=12+8+MetaDataByteLength)
32        {
[9ee2ce3]33                confparams_dec->losslessCompressor = is_lossless_compressed_data(cmpBytes, cmpSize);
[2c47b73]34                if(confparams_dec->szMode!=SZ_TEMPORAL_COMPRESSION)
35                {
[9ee2ce3]36                        if(confparams_dec->losslessCompressor!=-1)
[2c47b73]37                                confparams_dec->szMode = SZ_BEST_COMPRESSION;
38                        else
39                                confparams_dec->szMode = SZ_BEST_SPEED;                 
40                }
41                if(confparams_dec->szMode==SZ_BEST_SPEED)
42                {
43                        tmpSize = cmpSize;
44                        szTmpBytes = cmpBytes; 
45                }       
46                else if(confparams_dec->szMode==SZ_BEST_COMPRESSION || confparams_dec->szMode==SZ_DEFAULT_COMPRESSION || confparams_dec->szMode==SZ_TEMPORAL_COMPRESSION)
47                {
48                        if(targetUncompressSize<MIN_ZLIB_DEC_ALLOMEM_BYTES) //Considering the minimum size
49                                targetUncompressSize = MIN_ZLIB_DEC_ALLOMEM_BYTES;                     
[9ee2ce3]50                        tmpSize = sz_lossless_decompress(confparams_dec->losslessCompressor, cmpBytes, (unsigned long)cmpSize, &szTmpBytes, (unsigned long)targetUncompressSize+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE);                 
[2c47b73]51                        //szTmpBytes = (unsigned char*)malloc(sizeof(unsigned char)*tmpSize);
52                        //memcpy(szTmpBytes, tmpBytes, tmpSize);
53                        //free(tmpBytes); //release useless memory             
54                }
55                else
56                {
57                        printf("Wrong value of confparams_dec->szMode in the double compressed bytes.\n");
58                        status = SZ_MERR;
59                        return status;
60                }       
61        }
62        else
63                szTmpBytes = cmpBytes;
64        //TODO: convert szTmpBytes to double array.
65        TightDataPointStorageD* tdps;
66        int errBoundMode = new_TightDataPointStorageD_fromFlatBytes(&tdps, szTmpBytes, tmpSize);
67
68        int dim = computeDimension(r5,r4,r3,r2,r1);
69        int doubleSize = sizeof(double);
70        if(tdps->isLossless)
71        {
72                *newData = (double*)malloc(doubleSize*dataLength);
73                if(sysEndianType==BIG_ENDIAN_SYSTEM)
74                {
75                        memcpy(*newData, szTmpBytes+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, dataLength*doubleSize);
76                }
77                else
78                {
79                        unsigned char* p = szTmpBytes+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
80                        for(i=0;i<dataLength;i++,p+=doubleSize)
81                                (*newData)[i] = bytesToDouble(p);
82                }               
83        }
[9ee2ce3]84        else 
[2c47b73]85        {
[9ee2ce3]86                if(tdps->raBytes_size > 0) //v2.0
87                {
88                        if (dim == 1)
89                                getSnapshotData_double_1D(newData,r1,tdps, errBoundMode);
90                        else if(dim == 2)
91                                decompressDataSeries_double_2D_nonblocked_with_blocked_regression(newData, r2, r1, tdps->raBytes);
92                        else if(dim == 3)
93                                decompressDataSeries_double_3D_nonblocked_with_blocked_regression(newData, r3, r2, r1, tdps->raBytes);
94                        else if(dim == 4)
95                                decompressDataSeries_double_3D_nonblocked_with_blocked_regression(newData, r4*r3, r2, r1, tdps->raBytes);
96                        else
97                        {
98                                printf("Error: currently support only at most 4 dimensions!\n");
99                                status = SZ_DERR;
100                        }       
101                }
102                else //1.4.13
103                {
104                        if (dim == 1)
105                                getSnapshotData_double_1D(newData,r1,tdps, errBoundMode);
106                        else
107                        if (dim == 2)
108                                getSnapshotData_double_2D(newData,r2,r1,tdps, errBoundMode);
109                        else
110                        if (dim == 3)
111                                getSnapshotData_double_3D(newData,r3,r2,r1,tdps, errBoundMode);
112                        else
113                        if (dim == 4)
114                                getSnapshotData_double_4D(newData,r4,r3,r2,r1,tdps, errBoundMode);                     
115                        else
116                        {
117                                printf("Error: currently support only at most 4 dimensions!\n");
118                                status = SZ_DERR;
119                        }                       
120                }
121        }       
122
[2c47b73]123        free_TightDataPointStorageD2(tdps);
124        if(confparams_dec->szMode!=SZ_BEST_SPEED && cmpSize!=12+MetaDataByteLength+exe_params->SZ_SIZE_TYPE)
125                free(szTmpBytes);       
126        return status;
127}
128
129void decompressDataSeries_double_1D(double** data, size_t dataSeriesLength, TightDataPointStorageD* tdps) 
130{
131        updateQuantizationInfo(tdps->intervals);
132        size_t i, j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
133                                                                // in resiMidBits, p is to track the
134                                                                // byte_index of resiMidBits, l is for
135                                                                // leadNum
136        unsigned char* leadNum;
137        double interval = tdps->realPrecision*2;
138       
139        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
140        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
141
142        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
143
144        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
145        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
146        SZ_ReleaseHuffman(huffmanTree); 
147       
148        unsigned char preBytes[8];
149        unsigned char curBytes[8];
150       
151        memset(preBytes, 0, 8);
152
153        size_t curByteIndex = 0;
154        int reqBytesLength, resiBitsLength, resiBits; 
155        unsigned char leadingNum;       
156        double medianValue, exactData, predValue;
157       
158        reqBytesLength = tdps->reqLength/8;
159        resiBitsLength = tdps->reqLength%8;
160        medianValue = tdps->medianValue;
161       
162        int type_;
163        for (i = 0; i < dataSeriesLength; i++) {
164                type_ = type[i];
165                switch (type_) {
166                case 0:
167                        // compute resiBits
168                        resiBits = 0;
169                        if (resiBitsLength != 0) {
170                                int kMod8 = k % 8;
171                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
172                                if (rightMovSteps > 0) {
173                                        int code = getRightMovingCode(kMod8, resiBitsLength);
174                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
175                                } else if (rightMovSteps < 0) {
176                                        int code1 = getLeftMovingCode(kMod8);
177                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
178                                        int leftMovSteps = -rightMovSteps;
179                                        rightMovSteps = 8 - leftMovSteps;
180                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
181                                        p++;
182                                        resiBits = resiBits
183                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
184                                } else // rightMovSteps == 0
185                                {
186                                        int code = getRightMovingCode(kMod8, resiBitsLength);
187                                        resiBits = (tdps->residualMidBits[p] & code);
188                                        p++;
189                                }
190                                k += resiBitsLength;
191                        }
192
193                        // recover the exact data
194                        memset(curBytes, 0, 8);
195                        leadingNum = leadNum[l++];
196                        memcpy(curBytes, preBytes, leadingNum);
197                        for (j = leadingNum; j < reqBytesLength; j++)
198                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
199                        if (resiBitsLength != 0) {
200                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
201                                curBytes[reqBytesLength] = resiByte;
202                        }
203                       
204                        exactData = bytesToDouble(curBytes);
205                        (*data)[i] = exactData + medianValue;
206                        memcpy(preBytes,curBytes,8);
207                        break;
208                default:
209                        //predValue = 2 * (*data)[i-1] - (*data)[i-2];
210                        predValue = (*data)[i-1];
211                        (*data)[i] = predValue + (type_-exe_params->intvRadius)*interval;
212                        break;
213                }
214                //printf("%.30G\n",(*data)[i]);
215        }
216       
217#ifdef HAVE_TIMECMPR   
218        if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
219                memcpy(multisteps->hist_data, (*data), dataSeriesLength*sizeof(double));
220#endif 
221       
222        free(leadNum);
223        free(type);
224        return;
225}
226
227void decompressDataSeries_double_2D(double** data, size_t r1, size_t r2, TightDataPointStorageD* tdps) 
228{
229        updateQuantizationInfo(tdps->intervals);
230        //printf("tdps->intervals=%d, exe_params->intvRadius=%d\n", tdps->intervals, exe_params->intvRadius);
231       
232        size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
233        // in resiMidBits, p is to track the
234        // byte_index of resiMidBits, l is for
235        // leadNum
236        size_t dataSeriesLength = r1*r2;
237        //      printf ("%d %d\n", r1, r2);
238
239        unsigned char* leadNum;
240        double realPrecision = tdps->realPrecision;
241
242        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
243
244        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
245
246        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
247
248        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
249        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
250        SZ_ReleaseHuffman(huffmanTree); 
251
252        unsigned char preBytes[8];
253        unsigned char curBytes[8];
254
255        memset(preBytes, 0, 8);
256
257        size_t curByteIndex = 0;
258        int reqBytesLength, resiBitsLength, resiBits; 
259        unsigned char leadingNum;       
260        double medianValue, exactData;
261        int type_;
262
263        reqBytesLength = tdps->reqLength/8;
264        resiBitsLength = tdps->reqLength%8;
265        medianValue = tdps->medianValue;
266
267        double pred1D, pred2D;
268        size_t ii, jj;
269
270        /* Process Row-0, data 0 */
271
272        // compute resiBits
273        resiBits = 0;
274        if (resiBitsLength != 0) {
275                int kMod8 = k % 8;
276                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
277                if (rightMovSteps > 0) {
278                        int code = getRightMovingCode(kMod8, resiBitsLength);
279                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
280                } else if (rightMovSteps < 0) {
281                        int code1 = getLeftMovingCode(kMod8);
282                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
283                        int leftMovSteps = -rightMovSteps;
284                        rightMovSteps = 8 - leftMovSteps;
285                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
286                        p++;
287                        resiBits = resiBits
288                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
289                } else // rightMovSteps == 0
290                {
291                        int code = getRightMovingCode(kMod8, resiBitsLength);
292                        resiBits = (tdps->residualMidBits[p] & code);
293                        p++;
294                }
295                k += resiBitsLength;
296        }
297
298        // recover the exact data
299        memset(curBytes, 0, 8);
300        leadingNum = leadNum[l++];
301        memcpy(curBytes, preBytes, leadingNum);
302        for (j = leadingNum; j < reqBytesLength; j++)
303                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
304        if (resiBitsLength != 0) {
305                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
306                curBytes[reqBytesLength] = resiByte;
307        }
308
309        exactData = bytesToDouble(curBytes);
310        (*data)[0] = exactData + medianValue;
311        memcpy(preBytes,curBytes,8);
312
313        /* Process Row-0, data 1 */
314        type_ = type[1]; 
315        if (type_ != 0)
316        {
317                pred1D = (*data)[0];
318                (*data)[1] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
319        }
320        else
321        {
322                // compute resiBits
323                resiBits = 0;
324                if (resiBitsLength != 0) {
325                        int kMod8 = k % 8;
326                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
327                        if (rightMovSteps > 0) {
328                                int code = getRightMovingCode(kMod8, resiBitsLength);
329                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
330                        } else if (rightMovSteps < 0) {
331                                int code1 = getLeftMovingCode(kMod8);
332                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
333                                int leftMovSteps = -rightMovSteps;
334                                rightMovSteps = 8 - leftMovSteps;
335                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
336                                p++;
337                                resiBits = resiBits
338                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
339                        } else // rightMovSteps == 0
340                        {
341                                int code = getRightMovingCode(kMod8, resiBitsLength);
342                                resiBits = (tdps->residualMidBits[p] & code);
343                                p++;
344                        }
345                        k += resiBitsLength;
346                }
347
348                // recover the exact data
349                memset(curBytes, 0, 8);
350                leadingNum = leadNum[l++];
351                memcpy(curBytes, preBytes, leadingNum);
352                for (j = leadingNum; j < reqBytesLength; j++)
353                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
354                if (resiBitsLength != 0) {
355                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
356                        curBytes[reqBytesLength] = resiByte;
357                }
358               
359                exactData = bytesToDouble(curBytes);
360                (*data)[1] = exactData + medianValue;
361                memcpy(preBytes,curBytes,8);
362        }
363
364        /* Process Row-0, data 2 --> data r2-1 */
365        for (jj = 2; jj < r2; jj++)
366        {
367                type_ = type[jj];
368                if (type_ != 0)
369                {
370                        pred1D = 2*(*data)[jj-1] - (*data)[jj-2];                       
371                        (*data)[jj] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
372                }
373                else
374                {
375                        // compute resiBits
376                        resiBits = 0;
377                        if (resiBitsLength != 0) {
378                                int kMod8 = k % 8;
379                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
380                                if (rightMovSteps > 0) {
381                                        int code = getRightMovingCode(kMod8, resiBitsLength);
382                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
383                                } else if (rightMovSteps < 0) {
384                                        int code1 = getLeftMovingCode(kMod8);
385                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
386                                        int leftMovSteps = -rightMovSteps;
387                                        rightMovSteps = 8 - leftMovSteps;
388                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
389                                        p++;
390                                        resiBits = resiBits
391                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
392                                } else // rightMovSteps == 0
393                                {
394                                        int code = getRightMovingCode(kMod8, resiBitsLength);
395                                        resiBits = (tdps->residualMidBits[p] & code);
396                                        p++;
397                                }
398                                k += resiBitsLength;
399                        }
400
401                        // recover the exact data
402                        memset(curBytes, 0, 8);
403                        leadingNum = leadNum[l++];
404                        memcpy(curBytes, preBytes, leadingNum);
405                        for (j = leadingNum; j < reqBytesLength; j++)
406                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
407                        if (resiBitsLength != 0) {
408                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
409                                curBytes[reqBytesLength] = resiByte;
410                        }
411
412                        exactData = bytesToDouble(curBytes);
413                        (*data)[jj] = exactData + medianValue;
414                        memcpy(preBytes,curBytes,8);
415                }
416        }
417
418        size_t index;
419        /* Process Row-1 --> Row-r1-1 */
420        for (ii = 1; ii < r1; ii++)
421        {
422                /* Process row-ii data 0 */
423                index = ii*r2;
424
425                type_ = type[index];
426                if (type_ != 0)
427                {
428                        pred1D = (*data)[index-r2];
429                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
430                }
431                else
432                {
433                        // compute resiBits
434                        resiBits = 0;
435                        if (resiBitsLength != 0) {
436                                int kMod8 = k % 8;
437                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
438                                if (rightMovSteps > 0) {
439                                        int code = getRightMovingCode(kMod8, resiBitsLength);
440                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
441                                } else if (rightMovSteps < 0) {
442                                        int code1 = getLeftMovingCode(kMod8);
443                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
444                                        int leftMovSteps = -rightMovSteps;
445                                        rightMovSteps = 8 - leftMovSteps;
446                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
447                                        p++;
448                                        resiBits = resiBits
449                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
450                                } else // rightMovSteps == 0
451                                {
452                                        int code = getRightMovingCode(kMod8, resiBitsLength);
453                                        resiBits = (tdps->residualMidBits[p] & code);
454                                        p++;
455                                }
456                                k += resiBitsLength;
457                        }
458
459                        // recover the exact data
460                        memset(curBytes, 0, 8);
461                        leadingNum = leadNum[l++];
462                        memcpy(curBytes, preBytes, leadingNum);
463                        for (j = leadingNum; j < reqBytesLength; j++)
464                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
465                        if (resiBitsLength != 0) {
466                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
467                                curBytes[reqBytesLength] = resiByte;
468                        }
469
470                        exactData = bytesToDouble(curBytes);
471                        (*data)[index] = exactData + medianValue;
472                        memcpy(preBytes,curBytes,8);
473                }
474
475                /* Process row-ii data 1 --> r2-1*/
476                for (jj = 1; jj < r2; jj++)
477                {
478                        index = ii*r2+jj;
479                        pred2D = (*data)[index-1] + (*data)[index-r2] - (*data)[index-r2-1];
480
481                        type_ = type[index];
482                        if (type_ != 0)
483                        {
484                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
485                        }
486                        else
487                        {
488                                // compute resiBits
489                                resiBits = 0;
490                                if (resiBitsLength != 0) {
491                                        int kMod8 = k % 8;
492                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
493                                        if (rightMovSteps > 0) {
494                                                int code = getRightMovingCode(kMod8, resiBitsLength);
495                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
496                                        } else if (rightMovSteps < 0) {
497                                                int code1 = getLeftMovingCode(kMod8);
498                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
499                                                int leftMovSteps = -rightMovSteps;
500                                                rightMovSteps = 8 - leftMovSteps;
501                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
502                                                p++;
503                                                resiBits = resiBits
504                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
505                                        } else // rightMovSteps == 0
506                                        {
507                                                int code = getRightMovingCode(kMod8, resiBitsLength);
508                                                resiBits = (tdps->residualMidBits[p] & code);
509                                                p++;
510                                        }
511                                        k += resiBitsLength;
512                                }
513
514                                // recover the exact data
515                                memset(curBytes, 0, 8);
516                                leadingNum = leadNum[l++];
517                                memcpy(curBytes, preBytes, leadingNum);
518                                for (j = leadingNum; j < reqBytesLength; j++)
519                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
520                                if (resiBitsLength != 0) {
521                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
522                                        curBytes[reqBytesLength] = resiByte;
523                                }
524
525                                exactData = bytesToDouble(curBytes);
526                                (*data)[index] = exactData + medianValue;
527                                memcpy(preBytes,curBytes,8);
528                        }
529                }
530        }
531
532#ifdef HAVE_TIMECMPR   
533        if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
534                memcpy(multisteps->hist_data, (*data), dataSeriesLength*sizeof(double));
535#endif 
536
537        free(leadNum);
538        free(type);
539        return;
540}
541
542void decompressDataSeries_double_3D(double** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageD* tdps) 
543{
544        updateQuantizationInfo(tdps->intervals);
545        size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
546        // in resiMidBits, p is to track the
547        // byte_index of resiMidBits, l is for
548        // leadNum
549        size_t dataSeriesLength = r1*r2*r3;
550        size_t r23 = r2*r3;
551//      printf ("%d %d %d\n", r1, r2, r3);
552
553        unsigned char* leadNum;
554        double realPrecision = tdps->realPrecision;
555
556        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
557
558        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
559
560        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
561
562        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
563        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
564        SZ_ReleaseHuffman(huffmanTree); 
565
566        unsigned char preBytes[8];
567        unsigned char curBytes[8];
568
569        memset(preBytes, 0, 8);
570
571        size_t curByteIndex = 0;
572        int reqBytesLength, resiBitsLength, resiBits;
573        unsigned char leadingNum;
574        double medianValue, exactData;
575        int type_;
576
577        reqBytesLength = tdps->reqLength/8;
578        resiBitsLength = tdps->reqLength%8;
579        medianValue = tdps->medianValue;
580
581        double pred1D, pred2D, pred3D;
582        size_t ii, jj, kk;
583
584        ///////////////////////////     Process layer-0 ///////////////////////////
585        /* Process Row-0 data 0*/
586        // compute resiBits
587        resiBits = 0;
588        if (resiBitsLength != 0) {
589                int kMod8 = k % 8;
590                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
591                if (rightMovSteps > 0) {
592                        int code = getRightMovingCode(kMod8, resiBitsLength);
593                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
594                } else if (rightMovSteps < 0) {
595                        int code1 = getLeftMovingCode(kMod8);
596                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
597                        int leftMovSteps = -rightMovSteps;
598                        rightMovSteps = 8 - leftMovSteps;
599                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
600                        p++;
601                        resiBits = resiBits
602                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
603                } else // rightMovSteps == 0
604                {
605                        int code = getRightMovingCode(kMod8, resiBitsLength);
606                        resiBits = (tdps->residualMidBits[p] & code);
607                        p++;
608                }
609                k += resiBitsLength;
610        }
611
612        // recover the exact data
613        memset(curBytes, 0, 8);
614        leadingNum = leadNum[l++];
615        memcpy(curBytes, preBytes, leadingNum);
616        for (j = leadingNum; j < reqBytesLength; j++)
617                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
618        if (resiBitsLength != 0) {
619                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
620                curBytes[reqBytesLength] = resiByte;
621        }
622
623        exactData = bytesToDouble(curBytes);
624        (*data)[0] = exactData + medianValue;
625        memcpy(preBytes,curBytes,8);
626
627        /* Process Row-0, data 1 */
628        pred1D = (*data)[0];
629
630        type_ = type[1];
631        if (type_ != 0)
632        {
633                (*data)[1] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
634        }
635        else
636        {
637                // compute resiBits
638                resiBits = 0;
639                if (resiBitsLength != 0) {
640                        int kMod8 = k % 8;
641                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
642                        if (rightMovSteps > 0) {
643                                int code = getRightMovingCode(kMod8, resiBitsLength);
644                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
645                        } else if (rightMovSteps < 0) {
646                                int code1 = getLeftMovingCode(kMod8);
647                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
648                                int leftMovSteps = -rightMovSteps;
649                                rightMovSteps = 8 - leftMovSteps;
650                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
651                                p++;
652                                resiBits = resiBits
653                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
654                        } else // rightMovSteps == 0
655                        {
656                                int code = getRightMovingCode(kMod8, resiBitsLength);
657                                resiBits = (tdps->residualMidBits[p] & code);
658                                p++;
659                        }
660                        k += resiBitsLength;
661                }
662
663                // recover the exact data
664                memset(curBytes, 0, 8);
665                leadingNum = leadNum[l++];
666                memcpy(curBytes, preBytes, leadingNum);
667                for (j = leadingNum; j < reqBytesLength; j++)
668                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
669                if (resiBitsLength != 0) {
670                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
671                        curBytes[reqBytesLength] = resiByte;
672                }
673
674                exactData = bytesToDouble(curBytes);
675                (*data)[1] = exactData + medianValue;
676                memcpy(preBytes,curBytes,8);
677        }
678
679        /* Process Row-0, data 2 --> data r3-1 */
680        for (jj = 2; jj < r3; jj++)
681        {
682                pred1D = 2*(*data)[jj-1] - (*data)[jj-2];
683
684                type_ = type[jj];
685                if (type_ != 0)
686                {
687                        (*data)[jj] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
688                }
689                else
690                {
691                        // compute resiBits
692                        resiBits = 0;
693                        if (resiBitsLength != 0) {
694                                int kMod8 = k % 8;
695                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
696                                if (rightMovSteps > 0) {
697                                        int code = getRightMovingCode(kMod8, resiBitsLength);
698                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
699                                } else if (rightMovSteps < 0) {
700                                        int code1 = getLeftMovingCode(kMod8);
701                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
702                                        int leftMovSteps = -rightMovSteps;
703                                        rightMovSteps = 8 - leftMovSteps;
704                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
705                                        p++;
706                                        resiBits = resiBits
707                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
708                                } else // rightMovSteps == 0
709                                {
710                                        int code = getRightMovingCode(kMod8, resiBitsLength);
711                                        resiBits = (tdps->residualMidBits[p] & code);
712                                        p++;
713                                }
714                                k += resiBitsLength;
715                        }
716
717                        // recover the exact data
718                        memset(curBytes, 0, 8);
719                        leadingNum = leadNum[l++];
720                        memcpy(curBytes, preBytes, leadingNum);
721                        for (j = leadingNum; j < reqBytesLength; j++)
722                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
723                        if (resiBitsLength != 0) {
724                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
725                                curBytes[reqBytesLength] = resiByte;
726                        }
727
728                        exactData = bytesToDouble(curBytes);
729                        (*data)[jj] = exactData + medianValue;
730                        memcpy(preBytes,curBytes,8);
731                }
732        }
733
734        size_t index;
735        /* Process Row-1 --> Row-r2-1 */
736        for (ii = 1; ii < r2; ii++)
737        {
738                /* Process row-ii data 0 */
739                index = ii*r3;
740                pred1D = (*data)[index-r3];
741
742                type_ = type[index];
743                if (type_ != 0)
744                {
745                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
746                }
747                else
748                {
749                        // compute resiBits
750                        resiBits = 0;
751                        if (resiBitsLength != 0) {
752                                int kMod8 = k % 8;
753                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
754                                if (rightMovSteps > 0) {
755                                        int code = getRightMovingCode(kMod8, resiBitsLength);
756                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
757                                } else if (rightMovSteps < 0) {
758                                        int code1 = getLeftMovingCode(kMod8);
759                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
760                                        int leftMovSteps = -rightMovSteps;
761                                        rightMovSteps = 8 - leftMovSteps;
762                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
763                                        p++;
764                                        resiBits = resiBits
765                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
766                                } else // rightMovSteps == 0
767                                {
768                                        int code = getRightMovingCode(kMod8, resiBitsLength);
769                                        resiBits = (tdps->residualMidBits[p] & code);
770                                        p++;
771                                }
772                                k += resiBitsLength;
773                        }
774
775                        // recover the exact data
776                        memset(curBytes, 0, 8);
777                        leadingNum = leadNum[l++];
778                        memcpy(curBytes, preBytes, leadingNum);
779                        for (j = leadingNum; j < reqBytesLength; j++)
780                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
781                        if (resiBitsLength != 0) {
782                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
783                                curBytes[reqBytesLength] = resiByte;
784                        }
785
786                        exactData = bytesToDouble(curBytes);
787                        (*data)[index] = exactData + medianValue;
788                        memcpy(preBytes,curBytes,8);
789                }
790
791                /* Process row-ii data 1 --> r3-1*/
792                for (jj = 1; jj < r3; jj++)
793                {
794                        index = ii*r3+jj;
795                        pred2D = (*data)[index-1] + (*data)[index-r3] - (*data)[index-r3-1];
796
797                        type_ = type[index];
798                        if (type_ != 0)
799                        {
800                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
801                        }
802                        else
803                        {
804                                // compute resiBits
805                                resiBits = 0;
806                                if (resiBitsLength != 0) {
807                                        int kMod8 = k % 8;
808                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
809                                        if (rightMovSteps > 0) {
810                                                int code = getRightMovingCode(kMod8, resiBitsLength);
811                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
812                                        } else if (rightMovSteps < 0) {
813                                                int code1 = getLeftMovingCode(kMod8);
814                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
815                                                int leftMovSteps = -rightMovSteps;
816                                                rightMovSteps = 8 - leftMovSteps;
817                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
818                                                p++;
819                                                resiBits = resiBits
820                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
821                                        } else // rightMovSteps == 0
822                                        {
823                                                int code = getRightMovingCode(kMod8, resiBitsLength);
824                                                resiBits = (tdps->residualMidBits[p] & code);
825                                                p++;
826                                        }
827                                        k += resiBitsLength;
828                                }
829
830                                // recover the exact data
831                                memset(curBytes, 0, 8);
832                                leadingNum = leadNum[l++];
833                                memcpy(curBytes, preBytes, leadingNum);
834                                for (j = leadingNum; j < reqBytesLength; j++)
835                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
836                                if (resiBitsLength != 0) {
837                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
838                                        curBytes[reqBytesLength] = resiByte;
839                                }
840
841                                exactData = bytesToDouble(curBytes);
842                                (*data)[index] = exactData + medianValue;
843                                memcpy(preBytes,curBytes,8);
844                        }
845                }
846        }
847
848        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
849
850        for (kk = 1; kk < r1; kk++)
851        {
852                /* Process Row-0 data 0*/
853                index = kk*r23;
854                pred1D = (*data)[index-r23];
855
856                type_ = type[index];
857                if (type_ != 0)
858                {
859                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
860                }
861                else
862                {
863                        // compute resiBits
864                        resiBits = 0;
865                        if (resiBitsLength != 0) {
866                                int kMod8 = k % 8;
867                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
868                                if (rightMovSteps > 0) {
869                                        int code = getRightMovingCode(kMod8, resiBitsLength);
870                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
871                                } else if (rightMovSteps < 0) {
872                                        int code1 = getLeftMovingCode(kMod8);
873                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
874                                        int leftMovSteps = -rightMovSteps;
875                                        rightMovSteps = 8 - leftMovSteps;
876                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
877                                        p++;
878                                        resiBits = resiBits
879                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
880                                } else // rightMovSteps == 0
881                                {
882                                        int code = getRightMovingCode(kMod8, resiBitsLength);
883                                        resiBits = (tdps->residualMidBits[p] & code);
884                                        p++;
885                                }
886                                k += resiBitsLength;
887                        }
888
889                        // recover the exact data
890                        memset(curBytes, 0, 8);
891                        leadingNum = leadNum[l++];
892                        memcpy(curBytes, preBytes, leadingNum);
893                        for (j = leadingNum; j < reqBytesLength; j++)
894                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
895                        if (resiBitsLength != 0) {
896                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
897                                curBytes[reqBytesLength] = resiByte;
898                        }
899
900                        exactData = bytesToDouble(curBytes);
901                        (*data)[index] = exactData + medianValue;
902                        memcpy(preBytes,curBytes,8);
903                }
904
905                /* Process Row-0 data 1 --> data r3-1 */
906                for (jj = 1; jj < r3; jj++)
907                {
908                        index = kk*r23+jj;
909                        pred2D = (*data)[index-1] + (*data)[index-r23] - (*data)[index-r23-1];
910
911                        type_ = type[index];
912                        if (type_ != 0)
913                        {
914                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
915                        }
916                        else
917                        {
918                                // compute resiBits
919                                resiBits = 0;
920                                if (resiBitsLength != 0) {
921                                        int kMod8 = k % 8;
922                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
923                                        if (rightMovSteps > 0) {
924                                                int code = getRightMovingCode(kMod8, resiBitsLength);
925                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
926                                        } else if (rightMovSteps < 0) {
927                                                int code1 = getLeftMovingCode(kMod8);
928                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
929                                                int leftMovSteps = -rightMovSteps;
930                                                rightMovSteps = 8 - leftMovSteps;
931                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
932                                                p++;
933                                                resiBits = resiBits
934                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
935                                        } else // rightMovSteps == 0
936                                        {
937                                                int code = getRightMovingCode(kMod8, resiBitsLength);
938                                                resiBits = (tdps->residualMidBits[p] & code);
939                                                p++;
940                                        }
941                                        k += resiBitsLength;
942                                }
943
944                                // recover the exact data
945                                memset(curBytes, 0, 8);
946                                leadingNum = leadNum[l++];
947                                memcpy(curBytes, preBytes, leadingNum);
948                                for (j = leadingNum; j < reqBytesLength; j++)
949                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
950                                if (resiBitsLength != 0) {
951                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
952                                        curBytes[reqBytesLength] = resiByte;
953                                }
954
955                                exactData = bytesToDouble(curBytes);
956                                (*data)[index] = exactData + medianValue;
957                                memcpy(preBytes,curBytes,8);
958                        }
959                }
960
961                /* Process Row-1 --> Row-r2-1 */
962                for (ii = 1; ii < r2; ii++)
963                {
964                        /* Process Row-i data 0 */
965                        index = kk*r23 + ii*r3;
966                        pred2D = (*data)[index-r3] + (*data)[index-r23] - (*data)[index-r23-r3];
967
968                        type_ = type[index];
969                        if (type_ != 0)
970                        {
971                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
972                        }
973                        else
974                        {
975                                // compute resiBits
976                                resiBits = 0;
977                                if (resiBitsLength != 0) {
978                                        int kMod8 = k % 8;
979                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
980                                        if (rightMovSteps > 0) {
981                                                int code = getRightMovingCode(kMod8, resiBitsLength);
982                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
983                                        } else if (rightMovSteps < 0) {
984                                                int code1 = getLeftMovingCode(kMod8);
985                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
986                                                int leftMovSteps = -rightMovSteps;
987                                                rightMovSteps = 8 - leftMovSteps;
988                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
989                                                p++;
990                                                resiBits = resiBits
991                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
992                                        } else // rightMovSteps == 0
993                                        {
994                                                int code = getRightMovingCode(kMod8, resiBitsLength);
995                                                resiBits = (tdps->residualMidBits[p] & code);
996                                                p++;
997                                        }
998                                        k += resiBitsLength;
999                                }
1000
1001                                // recover the exact data
1002                                memset(curBytes, 0, 8);
1003                                leadingNum = leadNum[l++];
1004                                memcpy(curBytes, preBytes, leadingNum);
1005                                for (j = leadingNum; j < reqBytesLength; j++)
1006                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1007                                if (resiBitsLength != 0) {
1008                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1009                                        curBytes[reqBytesLength] = resiByte;
1010                                }
1011
1012                                exactData = bytesToDouble(curBytes);
1013                                (*data)[index] = exactData + medianValue;
1014                                memcpy(preBytes,curBytes,8);
1015                        }
1016
1017                        /* Process Row-i data 1 --> data r3-1 */
1018                        for (jj = 1; jj < r3; jj++)
1019                        {
1020                                index = kk*r23 + ii*r3 + jj;
1021                                pred3D = (*data)[index-1] + (*data)[index-r3] + (*data)[index-r23]
1022                                        - (*data)[index-r3-1] - (*data)[index-r23-r3] - (*data)[index-r23-1] + (*data)[index-r23-r3-1];
1023
1024                                type_ = type[index];
1025                                if (type_ != 0)
1026                                {
1027                                        (*data)[index] = pred3D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1028                                }
1029                                else
1030                                {
1031                                        // compute resiBits
1032                                        resiBits = 0;
1033                                        if (resiBitsLength != 0) {
1034                                                int kMod8 = k % 8;
1035                                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1036                                                if (rightMovSteps > 0) {
1037                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1038                                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1039                                                } else if (rightMovSteps < 0) {
1040                                                        int code1 = getLeftMovingCode(kMod8);
1041                                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1042                                                        int leftMovSteps = -rightMovSteps;
1043                                                        rightMovSteps = 8 - leftMovSteps;
1044                                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1045                                                        p++;
1046                                                        resiBits = resiBits
1047                                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1048                                                } else // rightMovSteps == 0
1049                                                {
1050                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1051                                                        resiBits = (tdps->residualMidBits[p] & code);
1052                                                        p++;
1053                                                }
1054                                                k += resiBitsLength;
1055                                        }
1056
1057                                        // recover the exact data
1058                                        memset(curBytes, 0, 8);
1059                                        leadingNum = leadNum[l++];
1060                                        memcpy(curBytes, preBytes, leadingNum);
1061                                        for (j = leadingNum; j < reqBytesLength; j++)
1062                                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1063                                        if (resiBitsLength != 0) {
1064                                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1065                                                curBytes[reqBytesLength] = resiByte;
1066                                        }
1067
1068                                        exactData = bytesToDouble(curBytes);
1069                                        (*data)[index] = exactData + medianValue;
1070                                        memcpy(preBytes,curBytes,8);
1071                                }
1072                        }
1073                }
1074        }
1075
1076#ifdef HAVE_TIMECMPR   
1077        if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1078                memcpy(multisteps->hist_data, (*data), dataSeriesLength*sizeof(double));
1079#endif 
1080
1081        free(leadNum);
1082        free(type);
1083        return;
1084}
1085
1086void decompressDataSeries_double_4D(double** data, size_t r1, size_t r2, size_t r3, size_t r4, TightDataPointStorageD* tdps)
1087{
1088        updateQuantizationInfo(tdps->intervals);
1089        size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
1090        // in resiMidBits, p is to track the
1091        // byte_index of resiMidBits, l is for
1092        // leadNum
1093        size_t dataSeriesLength = r1*r2*r3*r4;
1094        size_t r234 = r2*r3*r4;
1095        size_t r34 = r3*r4;
1096//      printf ("%d %d %d\n", r1, r2, r3, r4);
1097
1098        unsigned char* leadNum;
1099        double realPrecision = tdps->realPrecision;
1100
1101        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
1102
1103        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1104
1105        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
1106
1107        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
1108        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
1109        SZ_ReleaseHuffman(huffmanTree); 
1110
1111        unsigned char preBytes[8];
1112        unsigned char curBytes[8];
1113
1114        memset(preBytes, 0, 8);
1115
1116        size_t curByteIndex = 0;
1117        int reqBytesLength, resiBitsLength, resiBits;
1118        unsigned char leadingNum;
1119        double medianValue, exactData;
1120        int type_;
1121
1122        reqBytesLength = tdps->reqLength/8;
1123        resiBitsLength = tdps->reqLength%8;
1124        medianValue = tdps->medianValue;
1125
1126        double pred1D, pred2D, pred3D;
1127        size_t ii, jj, kk, ll;
1128        size_t index;
1129
1130        for (ll = 0; ll < r1; ll++)
1131        {
1132
1133                ///////////////////////////     Process layer-0 ///////////////////////////
1134                /* Process Row-0 data 0*/
1135                index = ll*r234;
1136
1137                // compute resiBits
1138                resiBits = 0;
1139                if (resiBitsLength != 0) {
1140                        int kMod8 = k % 8;
1141                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1142                        if (rightMovSteps > 0) {
1143                                int code = getRightMovingCode(kMod8, resiBitsLength);
1144                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1145                        } else if (rightMovSteps < 0) {
1146                                int code1 = getLeftMovingCode(kMod8);
1147                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1148                                int leftMovSteps = -rightMovSteps;
1149                                rightMovSteps = 8 - leftMovSteps;
1150                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1151                                p++;
1152                                resiBits = resiBits
1153                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1154                        } else // rightMovSteps == 0
1155                        {
1156                                int code = getRightMovingCode(kMod8, resiBitsLength);
1157                                resiBits = (tdps->residualMidBits[p] & code);
1158                                p++;
1159                        }
1160                        k += resiBitsLength;
1161                }
1162
1163                // recover the exact data
1164                memset(curBytes, 0, 8);
1165                leadingNum = leadNum[l++];
1166                memcpy(curBytes, preBytes, leadingNum);
1167                for (j = leadingNum; j < reqBytesLength; j++)
1168                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1169                if (resiBitsLength != 0) {
1170                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1171                        curBytes[reqBytesLength] = resiByte;
1172                }
1173
1174                exactData = bytesToDouble(curBytes);
1175                (*data)[index] = exactData + medianValue;
1176                memcpy(preBytes,curBytes,8);
1177
1178                /* Process Row-0, data 1 */
1179                index = ll*r234+1;
1180
1181                pred1D = (*data)[index-1];
1182
1183                type_ = type[index];
1184                if (type_ != 0)
1185                {
1186                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1187                }
1188                else
1189                {
1190                        // compute resiBits
1191                        resiBits = 0;
1192                        if (resiBitsLength != 0) {
1193                                int kMod8 = k % 8;
1194                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1195                                if (rightMovSteps > 0) {
1196                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1197                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1198                                } else if (rightMovSteps < 0) {
1199                                        int code1 = getLeftMovingCode(kMod8);
1200                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1201                                        int leftMovSteps = -rightMovSteps;
1202                                        rightMovSteps = 8 - leftMovSteps;
1203                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1204                                        p++;
1205                                        resiBits = resiBits
1206                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1207                                } else // rightMovSteps == 0
1208                                {
1209                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1210                                        resiBits = (tdps->residualMidBits[p] & code);
1211                                        p++;
1212                                }
1213                                k += resiBitsLength;
1214                        }
1215
1216                        // recover the exact data
1217                        memset(curBytes, 0, 8);
1218                        leadingNum = leadNum[l++];
1219                        memcpy(curBytes, preBytes, leadingNum);
1220                        for (j = leadingNum; j < reqBytesLength; j++)
1221                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1222                        if (resiBitsLength != 0) {
1223                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1224                                curBytes[reqBytesLength] = resiByte;
1225                        }
1226
1227                        exactData = bytesToDouble(curBytes);
1228                        (*data)[index] = exactData + medianValue;
1229                        memcpy(preBytes,curBytes,8);
1230                }
1231
1232                /* Process Row-0, data 2 --> data r4-1 */
1233                for (jj = 2; jj < r4; jj++)
1234                {
1235                        index = ll*r234+jj;
1236
1237                        pred1D = 2*(*data)[index-1] - (*data)[index-2];
1238
1239                        type_ = type[index];
1240                        if (type_ != 0)
1241                        {
1242                                (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1243                        }
1244                        else
1245                        {
1246                                // compute resiBits
1247                                resiBits = 0;
1248                                if (resiBitsLength != 0) {
1249                                        int kMod8 = k % 8;
1250                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1251                                        if (rightMovSteps > 0) {
1252                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1253                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1254                                        } else if (rightMovSteps < 0) {
1255                                                int code1 = getLeftMovingCode(kMod8);
1256                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1257                                                int leftMovSteps = -rightMovSteps;
1258                                                rightMovSteps = 8 - leftMovSteps;
1259                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1260                                                p++;
1261                                                resiBits = resiBits
1262                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1263                                        } else // rightMovSteps == 0
1264                                        {
1265                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1266                                                resiBits = (tdps->residualMidBits[p] & code);
1267                                                p++;
1268                                        }
1269                                        k += resiBitsLength;
1270                                }
1271
1272                                // recover the exact data
1273                                memset(curBytes, 0, 8);
1274                                leadingNum = leadNum[l++];
1275                                memcpy(curBytes, preBytes, leadingNum);
1276                                for (j = leadingNum; j < reqBytesLength; j++)
1277                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1278                                if (resiBitsLength != 0) {
1279                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1280                                        curBytes[reqBytesLength] = resiByte;
1281                                }
1282
1283                                exactData = bytesToDouble(curBytes);
1284                                (*data)[index] = exactData + medianValue;
1285                                memcpy(preBytes,curBytes,8);
1286                        }
1287                }
1288
1289                /* Process Row-1 --> Row-r3-1 */
1290                for (ii = 1; ii < r3; ii++)
1291                {
1292                        /* Process row-ii data 0 */
1293                        index = ll*r234+ii*r4;
1294
1295                        pred1D = (*data)[index-r4];
1296
1297                        type_ = type[index];
1298                        if (type_ != 0)
1299                        {
1300                                (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1301                        }
1302                        else
1303                        {
1304                                // compute resiBits
1305                                resiBits = 0;
1306                                if (resiBitsLength != 0) {
1307                                        int kMod8 = k % 8;
1308                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1309                                        if (rightMovSteps > 0) {
1310                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1311                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1312                                        } else if (rightMovSteps < 0) {
1313                                                int code1 = getLeftMovingCode(kMod8);
1314                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1315                                                int leftMovSteps = -rightMovSteps;
1316                                                rightMovSteps = 8 - leftMovSteps;
1317                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1318                                                p++;
1319                                                resiBits = resiBits
1320                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1321                                        } else // rightMovSteps == 0
1322                                        {
1323                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1324                                                resiBits = (tdps->residualMidBits[p] & code);
1325                                                p++;
1326                                        }
1327                                        k += resiBitsLength;
1328                                }
1329
1330                                // recover the exact data
1331                                memset(curBytes, 0, 8);
1332                                leadingNum = leadNum[l++];
1333                                memcpy(curBytes, preBytes, leadingNum);
1334                                for (j = leadingNum; j < reqBytesLength; j++)
1335                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1336                                if (resiBitsLength != 0) {
1337                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1338                                        curBytes[reqBytesLength] = resiByte;
1339                                }
1340
1341                                exactData = bytesToDouble(curBytes);
1342                                (*data)[index] = exactData + medianValue;
1343                                memcpy(preBytes,curBytes,8);
1344                        }
1345
1346                        /* Process row-ii data 1 --> r4-1*/
1347                        for (jj = 1; jj < r4; jj++)
1348                        {
1349                                index = ll*r234+ii*r4+jj;
1350
1351                                pred2D = (*data)[index-1] + (*data)[index-r4] - (*data)[index-r4-1];
1352
1353                                type_ = type[index];
1354                                if (type_ != 0)
1355                                {
1356                                        (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1357                                }
1358                                else
1359                                {
1360                                        // compute resiBits
1361                                        resiBits = 0;
1362                                        if (resiBitsLength != 0) {
1363                                                int kMod8 = k % 8;
1364                                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1365                                                if (rightMovSteps > 0) {
1366                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1367                                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1368                                                } else if (rightMovSteps < 0) {
1369                                                        int code1 = getLeftMovingCode(kMod8);
1370                                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1371                                                        int leftMovSteps = -rightMovSteps;
1372                                                        rightMovSteps = 8 - leftMovSteps;
1373                                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1374                                                        p++;
1375                                                        resiBits = resiBits
1376                                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1377                                                } else // rightMovSteps == 0
1378                                                {
1379                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1380                                                        resiBits = (tdps->residualMidBits[p] & code);
1381                                                        p++;
1382                                                }
1383                                                k += resiBitsLength;
1384                                        }
1385
1386                                        // recover the exact data
1387                                        memset(curBytes, 0, 8);
1388                                        leadingNum = leadNum[l++];
1389                                        memcpy(curBytes, preBytes, leadingNum);
1390                                        for (j = leadingNum; j < reqBytesLength; j++)
1391                                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1392                                        if (resiBitsLength != 0) {
1393                                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1394                                                curBytes[reqBytesLength] = resiByte;
1395                                        }
1396
1397                                        exactData = bytesToDouble(curBytes);
1398                                        (*data)[index] = exactData + medianValue;
1399                                        memcpy(preBytes,curBytes,8);
1400                                }
1401                        }
1402                }
1403
1404                ///////////////////////////     Process layer-1 --> layer-r2-1 ///////////////////////////
1405
1406                for (kk = 1; kk < r2; kk++)
1407                {
1408                        /* Process Row-0 data 0*/
1409                        index = ll*r234+kk*r34;
1410
1411                        pred1D = (*data)[index-r34];
1412
1413                        type_ = type[index];
1414                        if (type_ != 0)
1415                        {
1416                                (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1417                        }
1418                        else
1419                        {
1420                                // compute resiBits
1421                                resiBits = 0;
1422                                if (resiBitsLength != 0) {
1423                                        int kMod8 = k % 8;
1424                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1425                                        if (rightMovSteps > 0) {
1426                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1427                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1428                                        } else if (rightMovSteps < 0) {
1429                                                int code1 = getLeftMovingCode(kMod8);
1430                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1431                                                int leftMovSteps = -rightMovSteps;
1432                                                rightMovSteps = 8 - leftMovSteps;
1433                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1434                                                p++;
1435                                                resiBits = resiBits
1436                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1437                                        } else // rightMovSteps == 0
1438                                        {
1439                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1440                                                resiBits = (tdps->residualMidBits[p] & code);
1441                                                p++;
1442                                        }
1443                                        k += resiBitsLength;
1444                                }
1445
1446                                // recover the exact data
1447                                memset(curBytes, 0, 8);
1448                                leadingNum = leadNum[l++];
1449                                memcpy(curBytes, preBytes, leadingNum);
1450                                for (j = leadingNum; j < reqBytesLength; j++)
1451                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1452                                if (resiBitsLength != 0) {
1453                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1454                                        curBytes[reqBytesLength] = resiByte;
1455                                }
1456
1457                                exactData = bytesToDouble(curBytes);
1458                                (*data)[index] = exactData + medianValue;
1459                                memcpy(preBytes,curBytes,8);
1460                        }
1461
1462                        /* Process Row-0 data 1 --> data r4-1 */
1463                        for (jj = 1; jj < r4; jj++)
1464                        {
1465                                index = ll*r234+kk*r34+jj;
1466
1467                                pred2D = (*data)[index-1] + (*data)[index-r34] - (*data)[index-r34-1];
1468
1469                                type_ = type[index];
1470                                if (type_ != 0)
1471                                {
1472                                        (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1473                                }
1474                                else
1475                                {
1476                                        // compute resiBits
1477                                        resiBits = 0;
1478                                        if (resiBitsLength != 0) {
1479                                                int kMod8 = k % 8;
1480                                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1481                                                if (rightMovSteps > 0) {
1482                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1483                                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1484                                                } else if (rightMovSteps < 0) {
1485                                                        int code1 = getLeftMovingCode(kMod8);
1486                                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1487                                                        int leftMovSteps = -rightMovSteps;
1488                                                        rightMovSteps = 8 - leftMovSteps;
1489                                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1490                                                        p++;
1491                                                        resiBits = resiBits
1492                                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1493                                                } else // rightMovSteps == 0
1494                                                {
1495                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1496                                                        resiBits = (tdps->residualMidBits[p] & code);
1497                                                        p++;
1498                                                }
1499                                                k += resiBitsLength;
1500                                        }
1501
1502                                        // recover the exact data
1503                                        memset(curBytes, 0, 8);
1504                                        leadingNum = leadNum[l++];
1505                                        memcpy(curBytes, preBytes, leadingNum);
1506                                        for (j = leadingNum; j < reqBytesLength; j++)
1507                                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1508                                        if (resiBitsLength != 0) {
1509                                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1510                                                curBytes[reqBytesLength] = resiByte;
1511                                        }
1512
1513                                        exactData = bytesToDouble(curBytes);
1514                                        (*data)[index] = exactData + medianValue;
1515                                        memcpy(preBytes,curBytes,8);
1516                                }
1517                        }
1518
1519                        /* Process Row-1 --> Row-r3-1 */
1520                        for (ii = 1; ii < r3; ii++)
1521                        {
1522                                /* Process Row-i data 0 */
1523                                index = ll*r234+kk*r34+ii*r4;
1524
1525                                pred2D = (*data)[index-r4] + (*data)[index-r34] - (*data)[index-r34-r4];
1526
1527                                type_ = type[index];
1528                                if (type_ != 0)
1529                                {
1530                                        (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1531                                }
1532                                else
1533                                {
1534                                        // compute resiBits
1535                                        resiBits = 0;
1536                                        if (resiBitsLength != 0) {
1537                                                int kMod8 = k % 8;
1538                                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1539                                                if (rightMovSteps > 0) {
1540                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1541                                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1542                                                } else if (rightMovSteps < 0) {
1543                                                        int code1 = getLeftMovingCode(kMod8);
1544                                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1545                                                        int leftMovSteps = -rightMovSteps;
1546                                                        rightMovSteps = 8 - leftMovSteps;
1547                                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1548                                                        p++;
1549                                                        resiBits = resiBits
1550                                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1551                                                } else // rightMovSteps == 0
1552                                                {
1553                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1554                                                        resiBits = (tdps->residualMidBits[p] & code);
1555                                                        p++;
1556                                                }
1557                                                k += resiBitsLength;
1558                                        }
1559
1560                                        // recover the exact data
1561                                        memset(curBytes, 0, 8);
1562                                        leadingNum = leadNum[l++];
1563                                        memcpy(curBytes, preBytes, leadingNum);
1564                                        for (j = leadingNum; j < reqBytesLength; j++)
1565                                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1566                                        if (resiBitsLength != 0) {
1567                                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1568                                                curBytes[reqBytesLength] = resiByte;
1569                                        }
1570
1571                                        exactData = bytesToDouble(curBytes);
1572                                        (*data)[index] = exactData + medianValue;
1573                                        memcpy(preBytes,curBytes,8);
1574                                }
1575
1576                                /* Process Row-i data 1 --> data r4-1 */
1577                                for (jj = 1; jj < r4; jj++)
1578                                {
1579                                        index = ll*r234+kk*r34+ii*r4+jj;
1580
1581                                        pred3D = (*data)[index-1] + (*data)[index-r4] + (*data)[index-r34]
1582                                                        - (*data)[index-r4-1] - (*data)[index-r34-r4] - (*data)[index-r34-1] + (*data)[index-r34-r4-1];
1583
1584                                        type_ = type[index];
1585                                        if (type_ != 0)
1586                                        {
1587                                                (*data)[index] = pred3D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1588                                        }
1589                                        else
1590                                        {
1591                                                // compute resiBits
1592                                                resiBits = 0;
1593                                                if (resiBitsLength != 0) {
1594                                                        int kMod8 = k % 8;
1595                                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1596                                                        if (rightMovSteps > 0) {
1597                                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1598                                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1599                                                        } else if (rightMovSteps < 0) {
1600                                                                int code1 = getLeftMovingCode(kMod8);
1601                                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1602                                                                int leftMovSteps = -rightMovSteps;
1603                                                                rightMovSteps = 8 - leftMovSteps;
1604                                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1605                                                                p++;
1606                                                                resiBits = resiBits
1607                                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1608                                                        } else // rightMovSteps == 0
1609                                                        {
1610                                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1611                                                                resiBits = (tdps->residualMidBits[p] & code);
1612                                                                p++;
1613                                                        }
1614                                                        k += resiBitsLength;
1615                                                }
1616
1617                                                // recover the exact data
1618                                                memset(curBytes, 0, 8);
1619                                                leadingNum = leadNum[l++];
1620                                                memcpy(curBytes, preBytes, leadingNum);
1621                                                for (j = leadingNum; j < reqBytesLength; j++)
1622                                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1623                                                if (resiBitsLength != 0) {
1624                                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1625                                                        curBytes[reqBytesLength] = resiByte;
1626                                                }
1627
1628                                                exactData = bytesToDouble(curBytes);
1629                                                (*data)[index] = exactData + medianValue;
1630                                                memcpy(preBytes,curBytes,8);
1631                                        }
1632                                }
1633                        }
1634                }
1635        }
1636
1637//I didn't implement time-based compression for 4D actually.
1638//#ifdef HAVE_TIMECMPR 
1639//      if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1640//              memcpy(multisteps->hist_data, (*data), dataSeriesLength*sizeof(double));
1641//#endif       
1642
1643        free(leadNum);
1644        free(type);
1645        return;
1646}
1647
1648void getSnapshotData_double_1D(double** data, size_t dataSeriesLength, TightDataPointStorageD* tdps, int errBoundMode) 
1649{
1650        size_t i;
1651        if (tdps->allSameData) {
1652                double value = bytesToDouble(tdps->exactMidBytes);
1653                *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1654                for (i = 0; i < dataSeriesLength; i++)
1655                        (*data)[i] = value;
1656        } else {
1657                if (tdps->rtypeArray == NULL) {
1658                        if(errBoundMode < PW_REL)
1659                        {
1660#ifdef HAVE_TIMECMPR                           
1661                                if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1662                                {
1663                                        if(multisteps->compressionType == 0) //snapshot
1664                                                decompressDataSeries_double_1D(data, dataSeriesLength, tdps);
1665                                        else
1666                                                decompressDataSeries_double_1D_ts(data, dataSeriesLength, multisteps, tdps);                                   
1667                                }
1668                                else
1669#endif                                                         
1670                                        decompressDataSeries_double_1D(data, dataSeriesLength, tdps);
1671                        }
1672                        else 
1673                        {
[9ee2ce3]1674                                decompressDataSeries_double_1D_pwr_pre_log(data, dataSeriesLength, tdps);
1675                                //decompressDataSeries_double_1D_pwrgroup(data, dataSeriesLength, tdps);
[2c47b73]1676                        }
1677                        return;
1678                } else {
1679                        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1680                        // insert the reserved values
1681                        //int[] rtypes = TypeManager.convertByteArray2IntArray_fast_1b(
1682                        //              dataSeriesLength, rtypeArray);
1683                        int* rtypes;
1684                        int validLength = computeBitNumRequired(dataSeriesLength);
1685                        decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1686                        size_t count = 0;
1687                        for (i = 0; i < dataSeriesLength; i++) {
1688                                if (rtypes[i] == 1)
1689                                        (*data)[i] = tdps->reservedValue;
1690                                else
1691                                        count++;
1692                        }
1693                        // get the decompressed data
1694                        double* decmpData;
1695                        if(errBoundMode < PW_REL)
1696                                decompressDataSeries_double_1D(&decmpData, dataSeriesLength, tdps);
1697                        else 
[9ee2ce3]1698                                //decompressDataSeries_double_1D_pwr(&decmpData, dataSeriesLength, tdps);
1699                                decompressDataSeries_double_1D_pwr_pre_log(&decmpData, dataSeriesLength, tdps);
[2c47b73]1700                        // insert the decompressed data
1701                        size_t k = 0;
1702                        for (i = 0; i < dataSeriesLength; i++) {
1703                                if (rtypes[i] == 0) {
1704                                        (*data)[i] = decmpData[k++];
1705                                }
1706                        }
1707                        free(decmpData);
1708                        free(rtypes);
1709                }
1710        }
1711}
1712
1713void getSnapshotData_double_2D(double** data, size_t r1, size_t r2, TightDataPointStorageD* tdps, int errBoundMode) 
1714{
1715        size_t i;
1716        size_t dataSeriesLength = r1*r2;
1717        if (tdps->allSameData) {
1718                double value = bytesToDouble(tdps->exactMidBytes);
1719                *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1720                for (i = 0; i < dataSeriesLength; i++)
1721                        (*data)[i] = value;
1722        } else {
1723                if (tdps->rtypeArray == NULL) {
1724                        if(errBoundMode < PW_REL)
1725                        {
1726#ifdef HAVE_TIMECMPR                           
1727                                if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1728                                {
1729                                        if(multisteps->compressionType == 0) //snapshot
1730                                                decompressDataSeries_double_2D(data, r1, r2, tdps);
1731                                        else
1732                                                decompressDataSeries_double_1D_ts(data, dataSeriesLength, multisteps, tdps);                                   
1733                                }
1734                                else
1735#endif                                         
1736                                        decompressDataSeries_double_2D(data, r1, r2, tdps);
1737                        }
1738                        else 
[9ee2ce3]1739                                //decompressDataSeries_double_2D_pwr(data, r1, r2, tdps);
1740                                decompressDataSeries_double_2D_pwr_pre_log(data, r1, r2, tdps);
[2c47b73]1741                        return;
1742                } else {
1743                        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1744                        // insert the reserved values
1745                        //int[] rtypes = TypeManager.convertByteArray2IntArray_fast_1b(
1746                        //              dataSeriesLength, rtypeArray);
1747                        int* rtypes;
1748                        int validLength = computeBitNumRequired(dataSeriesLength);
1749                        decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1750                        size_t count = 0;
1751                        for (i = 0; i < dataSeriesLength; i++) {
1752                                if (rtypes[i] == 1)
1753                                        (*data)[i] = tdps->reservedValue;
1754                                else
1755                                        count++;
1756                        }
1757                        // get the decompressed data
1758                        double* decmpData;
1759                        if(errBoundMode < PW_REL)
1760                                decompressDataSeries_double_2D(&decmpData, r1, r2, tdps);
1761                        else 
[9ee2ce3]1762                                //decompressDataSeries_double_2D_pwr(&decmpData, r1, r2, tdps);
1763                                decompressDataSeries_double_2D_pwr_pre_log(&decmpData, r1, r2, tdps);
[2c47b73]1764                        // insert the decompressed data
1765                        size_t k = 0;
1766                        for (i = 0; i < dataSeriesLength; i++) {
1767                                if (rtypes[i] == 0) {
1768                                        (*data)[i] = decmpData[k++];
1769                                }
1770                        }
1771                        free(decmpData);
1772                        free(rtypes);
1773                }
1774        }
1775}
1776
1777void getSnapshotData_double_3D(double** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageD* tdps, int errBoundMode) 
1778{
1779        size_t i;
1780        size_t dataSeriesLength = r1*r2*r3;
1781        if (tdps->allSameData) {
1782                double value = bytesToDouble(tdps->exactMidBytes);
1783                *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1784                for (i = 0; i < dataSeriesLength; i++)
1785                        (*data)[i] = value;
1786        } else {
1787                if (tdps->rtypeArray == NULL) {
1788                        if(errBoundMode < PW_REL)
1789                        {
1790#ifdef HAVE_TIMECMPR                           
1791                                if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1792                                {
1793                                        if(multisteps->compressionType == 0) //snapshot
1794                                                decompressDataSeries_double_3D(data, r1, r2, r3, tdps);
1795                                        else
1796                                                decompressDataSeries_double_1D_ts(data, dataSeriesLength, multisteps, tdps);                                   
1797                                }
1798                                else
1799#endif                                         
1800                                        decompressDataSeries_double_3D(data, r1, r2, r3, tdps);
1801                        }
1802                        else 
[9ee2ce3]1803                                //decompressDataSeries_double_3D_pwr(data, r1, r2, r3, tdps);
1804                                decompressDataSeries_double_3D_pwr_pre_log(data, r1, r2, r3, tdps);
[2c47b73]1805                        return;
1806                } else {
1807                        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1808                        // insert the reserved values
1809                        //int[] rtypes = TypeManager.convertByteArray2IntArray_fast_1b(
1810                        //              dataSeriesLength, rtypeArray);
1811                        int* rtypes;
1812                        int validLength = computeBitNumRequired(dataSeriesLength);
1813                        decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1814                        size_t count = 0;
1815                        for (i = 0; i < dataSeriesLength; i++) {
1816                                if (rtypes[i] == 1)
1817                                        (*data)[i] = tdps->reservedValue;
1818                                else
1819                                        count++;
1820                        }
1821                        // get the decompressed data
1822                        double* decmpData;
1823                        if(errBoundMode < PW_REL)
1824                                decompressDataSeries_double_3D(&decmpData, r1, r2, r3, tdps);
1825                        else 
[9ee2ce3]1826                                //decompressDataSeries_double_3D_pwr(&decmpData, r1, r2, r3, tdps);                     
1827                                decompressDataSeries_double_3D_pwr_pre_log(&decmpData, r1, r2, r3, tdps);                       
[2c47b73]1828                        // insert the decompressed data
1829                        size_t k = 0;
1830                        for (i = 0; i < dataSeriesLength; i++) {
1831                                if (rtypes[i] == 0) {
1832                                        (*data)[i] = decmpData[k++];
1833                                }
1834                        }
1835                        free(decmpData);
1836                        free(rtypes);
1837                }
1838        }
1839}
1840
1841void getSnapshotData_double_4D(double** data, size_t r1, size_t r2, size_t r3, size_t r4, TightDataPointStorageD* tdps, int errBoundMode)
1842{
1843        size_t i;
1844        size_t dataSeriesLength = r1*r2*r3*r4;
1845        if (tdps->allSameData) {
1846                double value = bytesToDouble(tdps->exactMidBytes);
1847                *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1848                for (i = 0; i < dataSeriesLength; i++)
1849                        (*data)[i] = value;
1850        } else {
1851                if (tdps->rtypeArray == NULL) {
1852                        if(errBoundMode < PW_REL)
1853                        {
1854#ifdef HAVE_TIMECMPR                                   
1855                                if(confparams_dec->szMode == SZ_TEMPORAL_COMPRESSION)
1856                                {
1857                                        if(multisteps->compressionType == 0)
1858                                                decompressDataSeries_double_4D(data, r1, r2, r3, r4, tdps);
1859                                        else
1860                                                decompressDataSeries_double_1D_ts(data, r1*r2*r3*r4, multisteps, tdps);                                 
1861                                }
1862                                else
1863#endif                         
1864                                        decompressDataSeries_double_4D(data, r1, r2, r3, r4, tdps);
1865                        }
1866                        else 
1867                        {
[9ee2ce3]1868                                //decompressDataSeries_double_3D_pwr(data, r1*r2, r3, r4, tdps);
1869                                decompressDataSeries_double_3D_pwr_pre_log(data, r1*r2, r3, r4, tdps);
[2c47b73]1870                                //ToDO
1871                                //decompressDataSeries_double_4D_pwr(data, r1, r2, r3, r4, tdps);
1872                        }                                       
1873                        return;
1874                } else {
1875                        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1876                        int* rtypes;
1877                        int validLength = computeBitNumRequired(dataSeriesLength);
1878                        decompressBitArraybySimpleLZ77(&rtypes, tdps->rtypeArray, tdps->rtypeArray_size, dataSeriesLength, validLength);
1879                        size_t count = 0;
1880                        for (i = 0; i < dataSeriesLength; i++) {
1881                                if (rtypes[i] == 1)
1882                                        (*data)[i] = tdps->reservedValue;
1883                                else
1884                                        count++;
1885                        }
1886                        // get the decompressed data
1887                        double* decmpData;
1888                        if(errBoundMode < PW_REL)
1889                                decompressDataSeries_double_4D(&decmpData, r1, r2, r3, r4, tdps);
1890                        else
[9ee2ce3]1891                                //decompressDataSeries_double_3D_pwr(&decmpData, r1*r2, r3, r4, tdps);
1892                                decompressDataSeries_double_3D_pwr_pre_log(&decmpData, r1*r2, r3, r4, tdps);
[2c47b73]1893                                //ToDo
1894                                //decompressDataSeries_double_4D_pwr(&decmpData, r1, r2, r3, r4, tdps);
1895                        // insert the decompressed data
1896                        size_t k = 0;
1897                        for (i = 0; i < dataSeriesLength; i++) {
1898                                if (rtypes[i] == 0) {
1899                                        (*data)[i] = decmpData[k++];
1900                                }
1901                        }
1902                        free(decmpData);
1903                        free(rtypes);
1904                }
1905        }
1906}
[9ee2ce3]1907
1908void decompressDataSeries_double_2D_nonblocked_with_blocked_regression(double** data, size_t r1, size_t r2, unsigned char* comp_data){
1909
1910        size_t dim0_offset = r2;
1911        size_t num_elements = r1 * r2;
1912
1913        *data = (double*)malloc(sizeof(double)*num_elements);
1914
1915        unsigned char * comp_data_pos = comp_data;
1916
1917        size_t block_size = bytesToInt_bigEndian(comp_data_pos);
1918        comp_data_pos += sizeof(int);
1919        // calculate block dims
1920        size_t num_x, num_y;
1921        SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r1, num_x, block_size);
1922        SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r2, num_y, block_size);
1923
1924        size_t split_index_x, split_index_y;
1925        size_t early_blockcount_x, early_blockcount_y;
1926        size_t late_blockcount_x, late_blockcount_y;
1927        SZ_COMPUTE_BLOCKCOUNT(r1, num_x, split_index_x, early_blockcount_x, late_blockcount_x);
1928        SZ_COMPUTE_BLOCKCOUNT(r2, num_y, split_index_y, early_blockcount_y, late_blockcount_y);
1929
1930        size_t num_blocks = num_x * num_y;
1931
1932        double realPrecision = bytesToDouble(comp_data_pos);
1933        comp_data_pos += sizeof(double);
1934        unsigned int intervals = bytesToInt_bigEndian(comp_data_pos);
1935        comp_data_pos += sizeof(int);
1936
1937        updateQuantizationInfo(intervals);
1938
1939        unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
1940        comp_data_pos += sizeof(int);
1941
1942        int stateNum = 2*intervals;
1943        HuffmanTree* huffmanTree = createHuffmanTree(stateNum);
1944       
1945        int nodeCount = bytesToInt_bigEndian(comp_data_pos);
1946       
1947        node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree,comp_data_pos+sizeof(int), nodeCount);
1948        comp_data_pos += sizeof(int) + tree_size;
1949
1950        double mean;
1951        unsigned char use_mean;
1952        memcpy(&use_mean, comp_data_pos, sizeof(unsigned char));
1953        comp_data_pos += sizeof(unsigned char);
1954        memcpy(&mean, comp_data_pos, sizeof(double));
1955        comp_data_pos += sizeof(double);
1956        size_t reg_count = 0;
1957
1958        unsigned char * indicator;
1959        size_t indicator_bitlength = (num_blocks - 1)/8 + 1;
1960        convertByteArray2IntArray_fast_1b(num_blocks, comp_data_pos, indicator_bitlength, &indicator);
1961        comp_data_pos += indicator_bitlength;
1962        for(size_t i=0; i<num_blocks; i++){
1963                if(!indicator[i]) reg_count ++;
1964        }
1965        //printf("reg_count: %ld\n", reg_count);
1966
1967        int coeff_intvRadius[3];
1968        int * coeff_result_type = (int *) malloc(num_blocks*3*sizeof(int));
1969        int * coeff_type[3];
1970        double precision[3];
1971        double * coeff_unpred_data[3];
1972        if(reg_count > 0){
1973                for(int i=0; i<3; i++){
1974                        precision[i] = bytesToDouble(comp_data_pos);
1975                        comp_data_pos += sizeof(double);
1976                        coeff_intvRadius[i] = bytesToInt_bigEndian(comp_data_pos);
1977                        comp_data_pos += sizeof(int);
1978                        unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
1979                        comp_data_pos += sizeof(int);
1980                        int stateNum = 2*coeff_intvRadius[i]*2;
1981                        HuffmanTree* huffmanTree = createHuffmanTree(stateNum); 
1982                        int nodeCount = bytesToInt_bigEndian(comp_data_pos);
1983                        node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree, comp_data_pos+sizeof(int), nodeCount);
1984                        comp_data_pos += sizeof(int) + tree_size;
1985
1986                        coeff_type[i] = coeff_result_type + i * num_blocks;
1987                        size_t typeArray_size = bytesToSize(comp_data_pos);
1988                        decode(comp_data_pos + sizeof(size_t), reg_count, root, coeff_type[i]);
1989                        comp_data_pos += sizeof(size_t) + typeArray_size;
1990                        int coeff_unpred_count = bytesToInt_bigEndian(comp_data_pos);
1991                        comp_data_pos += sizeof(int);
1992                        coeff_unpred_data[i] = (double *) comp_data_pos;
1993                        comp_data_pos += coeff_unpred_count * sizeof(double);
1994                        SZ_ReleaseHuffman(huffmanTree);
1995                }
1996        }
1997        double last_coefficients[3] = {0.0};
1998        int coeff_unpred_data_count[3] = {0};
1999        int coeff_index = 0;
2000        updateQuantizationInfo(intervals);
2001
2002        size_t total_unpred;
2003        memcpy(&total_unpred, comp_data_pos, sizeof(size_t));
2004        comp_data_pos += sizeof(size_t);
2005        double * unpred_data = (double *) comp_data_pos;
2006        comp_data_pos += total_unpred * sizeof(double);
2007
2008        int * result_type = (int *) malloc(num_elements * sizeof(int));
2009        decode(comp_data_pos, num_elements, root, result_type);
2010        SZ_ReleaseHuffman(huffmanTree);
2011       
2012        int intvRadius = exe_params->intvRadius;
2013       
2014        int * type;
2015
2016        double * data_pos = *data;
2017        size_t offset_x, offset_y;
2018        size_t current_blockcount_x, current_blockcount_y;
2019        size_t cur_unpred_count;
2020
2021        unsigned char * indicator_pos = indicator;
2022        if(use_mean){
2023                type = result_type;
2024                for(size_t i=0; i<num_x; i++){
2025                        for(size_t j=0; j<num_y; j++){
2026                                offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
2027                                offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2028                                data_pos = *data + offset_x * dim0_offset + offset_y;
2029
2030                                current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
2031                                current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2032
2033                                size_t current_block_elements = current_blockcount_x * current_blockcount_y;
2034                                if(*indicator_pos){
2035                                        // decompress by SZ
2036
2037                                        double * block_data_pos = data_pos;
2038                                        double pred;
2039                                        size_t index = 0;
2040                                        int type_;
2041                                        // d11 is current data
2042                                        size_t unpredictable_count = 0;
2043                                        double d00, d01, d10;
2044                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
2045                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
2046                                                        type_ = type[index];
2047                                                        if(type_ == intvRadius){
2048                                                                *block_data_pos = mean;
2049                                                        }
2050                                                        else if(type_ == 0){
2051                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2052                                                        }
2053                                                        else{
2054                                                                d00 = d01 = d10 = 1;
2055                                                                if(i == 0 && ii == 0){
2056                                                                        d00 = d01 = 0;
2057                                                                }
2058                                                                if(j == 0 && jj == 0){
2059                                                                        d00 = d10 = 0;
2060                                                                }
2061                                                                if(d00){
2062                                                                        d00 = block_data_pos[- dim0_offset - 1];
2063                                                                }
2064                                                                if(d01){
2065                                                                        d01 = block_data_pos[- dim0_offset];
2066                                                                }
2067                                                                if(d10){
2068                                                                        d10 = block_data_pos[- 1];
2069                                                                }
2070                                                                if(type_ < intvRadius) type_ += 1;
2071                                                                pred = d10 + d01 - d00;
2072                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2073                                                        }
2074                                                        index ++;
2075                                                        block_data_pos ++;
2076                                                }
2077                                                block_data_pos += dim0_offset - current_blockcount_y;
2078                                        }
2079                                        cur_unpred_count = unpredictable_count;
2080                                }
2081                                else{
2082                                        // decompress by regression
2083                                        {
2084                                                //restore regression coefficients
2085                                                double pred;
2086                                                int type_;
2087                                                for(int e=0; e<3; e++){
2088                                                        type_ = coeff_type[e][coeff_index];
2089                                                        if (type_ != 0){
2090                                                                pred = last_coefficients[e];
2091                                                                last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2092                                                        }
2093                                                        else{
2094                                                                last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2095                                                                coeff_unpred_data_count[e] ++;
2096                                                        }
2097                                                }
2098                                                coeff_index ++;
2099                                        }
2100                                        {
2101                                                double * block_data_pos = data_pos;
2102                                                double pred;
2103                                                int type_;
2104                                                size_t index = 0;
2105                                                size_t unpredictable_count = 0;
2106                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
2107                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
2108                                                                type_ = type[index];
2109                                                                if (type_ != 0){
2110                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2];
2111                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2112                                                                }
2113                                                                else{
2114                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
2115                                                                }
2116
2117                                                                index ++;       
2118                                                                block_data_pos ++;
2119                                                        }
2120                                                        block_data_pos += dim0_offset - current_blockcount_y;
2121                                                }
2122                                                cur_unpred_count = unpredictable_count;
2123                                        }
2124                                }
2125
2126                                type += current_block_elements;
2127                                indicator_pos ++;
2128                                unpred_data += cur_unpred_count;
2129                        }
2130                }
2131        }
2132        else{
2133                type = result_type;
2134                for(size_t i=0; i<num_x; i++){
2135                        for(size_t j=0; j<num_y; j++){
2136                                offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
2137                                offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2138                                data_pos = *data + offset_x * dim0_offset + offset_y;
2139
2140                                current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
2141                                current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2142
2143                                size_t current_block_elements = current_blockcount_x * current_blockcount_y;
2144                                if(*indicator_pos){
2145                                        // decompress by SZ
2146                                       
2147                                        double * block_data_pos = data_pos;
2148                                        double pred;
2149                                        size_t index = 0;
2150                                        int type_;
2151                                        // d11 is current data
2152                                        size_t unpredictable_count = 0;
2153                                        double d00, d01, d10;
2154                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
2155                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
2156                                                        type_ = type[index];
2157                                                        if(type_ == 0){
2158                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2159                                                        }
2160                                                        else{
2161                                                                d00 = d01 = d10 = 1;
2162                                                                if(i == 0 && ii == 0){
2163                                                                        d00 = d01 = 0;
2164                                                                }
2165                                                                if(j == 0 && jj == 0){
2166                                                                        d00 = d10 = 0;
2167                                                                }
2168                                                                if(d00){
2169                                                                        d00 = block_data_pos[- dim0_offset - 1];
2170                                                                }
2171                                                                if(d01){
2172                                                                        d01 = block_data_pos[- dim0_offset];
2173                                                                }
2174                                                                if(d10){
2175                                                                        d10 = block_data_pos[- 1];
2176                                                                }
2177                                                                pred = d10 + d01 - d00;
2178                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2179                                                        }
2180                                                        index ++;
2181                                                        block_data_pos ++;
2182                                                }
2183                                                block_data_pos += dim0_offset - current_blockcount_y;
2184                                        }
2185                                        cur_unpred_count = unpredictable_count;
2186                                }
2187                                else{
2188                                        // decompress by regression
2189                                        {
2190                                                //restore regression coefficients
2191                                                double pred;
2192                                                int type_;
2193                                                for(int e=0; e<3; e++){
2194                                                        type_ = coeff_type[e][coeff_index];
2195                                                        if (type_ != 0){
2196                                                                pred = last_coefficients[e];
2197                                                                last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2198                                                        }
2199                                                        else{
2200                                                                last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2201                                                                coeff_unpred_data_count[e] ++;
2202                                                        }
2203                                                }
2204                                                coeff_index ++;
2205                                        }
2206                                        {
2207                                                double * block_data_pos = data_pos;
2208                                                double pred;
2209                                                int type_;
2210                                                size_t index = 0;
2211                                                size_t unpredictable_count = 0;
2212                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
2213                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
2214                                                                type_ = type[index];
2215                                                                if (type_ != 0){
2216                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2];
2217                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2218                                                                }
2219                                                                else{
2220                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
2221                                                                }
2222                                                                index ++;       
2223                                                                block_data_pos ++;
2224                                                        }
2225                                                        block_data_pos += dim0_offset - current_blockcount_y;
2226                                                }
2227                                                cur_unpred_count = unpredictable_count;
2228                                        }
2229                                }
2230
2231                                type += current_block_elements;
2232                                indicator_pos ++;
2233                                unpred_data += cur_unpred_count;
2234                        }
2235                }
2236        }
2237        free(coeff_result_type);
2238
2239        free(indicator);
2240        free(result_type);
2241}
2242
2243
2244void decompressDataSeries_double_3D_nonblocked_with_blocked_regression(double** data, size_t r1, size_t r2, size_t r3, unsigned char* comp_data){
2245
2246        size_t dim0_offset = r2 * r3;
2247        size_t dim1_offset = r3;
2248        size_t num_elements = r1 * r2 * r3;
2249
2250        *data = (double*)malloc(sizeof(double)*num_elements);
2251
2252        unsigned char * comp_data_pos = comp_data;
2253
2254        size_t block_size = bytesToInt_bigEndian(comp_data_pos);
2255        comp_data_pos += sizeof(int);
2256        // calculate block dims
2257        size_t num_x, num_y, num_z;
2258        SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r1, num_x, block_size);
2259        SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r2, num_y, block_size);
2260        SZ_COMPUTE_3D_NUMBER_OF_BLOCKS(r3, num_z, block_size);
2261
2262        size_t split_index_x, split_index_y, split_index_z;
2263        size_t early_blockcount_x, early_blockcount_y, early_blockcount_z;
2264        size_t late_blockcount_x, late_blockcount_y, late_blockcount_z;
2265        SZ_COMPUTE_BLOCKCOUNT(r1, num_x, split_index_x, early_blockcount_x, late_blockcount_x);
2266        SZ_COMPUTE_BLOCKCOUNT(r2, num_y, split_index_y, early_blockcount_y, late_blockcount_y);
2267        SZ_COMPUTE_BLOCKCOUNT(r3, num_z, split_index_z, early_blockcount_z, late_blockcount_z);
2268
2269        size_t num_blocks = num_x * num_y * num_z;
2270
2271        double realPrecision = bytesToDouble(comp_data_pos);
2272        comp_data_pos += sizeof(double);
2273        unsigned int intervals = bytesToInt_bigEndian(comp_data_pos);
2274        comp_data_pos += sizeof(int);
2275
2276        updateQuantizationInfo(intervals);
2277
2278        unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
2279        comp_data_pos += sizeof(int);
2280       
2281        int stateNum = 2*intervals;
2282        HuffmanTree* huffmanTree = createHuffmanTree(stateNum); 
2283       
2284        int nodeCount = bytesToInt_bigEndian(comp_data_pos);
2285        node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree,comp_data_pos+4, nodeCount);
2286        comp_data_pos += sizeof(int) + tree_size;
2287
2288        double mean;
2289        unsigned char use_mean;
2290        memcpy(&use_mean, comp_data_pos, sizeof(unsigned char));
2291        comp_data_pos += sizeof(unsigned char);
2292        memcpy(&mean, comp_data_pos, sizeof(double));
2293        comp_data_pos += sizeof(double);
2294        size_t reg_count = 0;
2295
2296        unsigned char * indicator;
2297        size_t indicator_bitlength = (num_blocks - 1)/8 + 1;
2298        convertByteArray2IntArray_fast_1b(num_blocks, comp_data_pos, indicator_bitlength, &indicator);
2299        comp_data_pos += indicator_bitlength;
2300        for(size_t i=0; i<num_blocks; i++){
2301                if(!indicator[i]) reg_count ++;
2302        }
2303
2304        int coeff_intvRadius[4];
2305        int * coeff_result_type = (int *) malloc(num_blocks*4*sizeof(int));
2306        int * coeff_type[4];
2307        double precision[4];
2308        double * coeff_unpred_data[4];
2309        if(reg_count > 0){
2310                for(int i=0; i<4; i++){
2311                        precision[i] = bytesToDouble(comp_data_pos);
2312                        comp_data_pos += sizeof(double);
2313                        coeff_intvRadius[i] = bytesToInt_bigEndian(comp_data_pos);
2314                        comp_data_pos += sizeof(int);
2315                        unsigned int tree_size = bytesToInt_bigEndian(comp_data_pos);
2316                        comp_data_pos += sizeof(int);
2317                        int stateNum = 2*coeff_intvRadius[i]*2;
2318                        HuffmanTree* huffmanTree = createHuffmanTree(stateNum); 
2319                        int nodeCount = bytesToInt_bigEndian(comp_data_pos);
2320                        node root = reconstruct_HuffTree_from_bytes_anyStates(huffmanTree, comp_data_pos+4, nodeCount);
2321                        comp_data_pos += sizeof(int) + tree_size;
2322
2323                        coeff_type[i] = coeff_result_type + i * num_blocks;
2324                        size_t typeArray_size = bytesToSize(comp_data_pos);
2325                        decode(comp_data_pos + sizeof(size_t), reg_count, root, coeff_type[i]);
2326                        comp_data_pos += sizeof(size_t) + typeArray_size;
2327                        int coeff_unpred_count = bytesToInt_bigEndian(comp_data_pos);
2328                        comp_data_pos += sizeof(int);
2329                        coeff_unpred_data[i] = (double *) comp_data_pos;
2330                        comp_data_pos += coeff_unpred_count * sizeof(double);
2331                        SZ_ReleaseHuffman(huffmanTree);
2332                }
2333        }
2334        double last_coefficients[4] = {0.0};
2335        int coeff_unpred_data_count[4] = {0};
2336        int coeff_index = 0;
2337        updateQuantizationInfo(intervals);
2338
2339        size_t total_unpred;
2340        memcpy(&total_unpred, comp_data_pos, sizeof(size_t));
2341        comp_data_pos += sizeof(size_t);
2342        double * unpred_data = (double *) comp_data_pos;
2343        comp_data_pos += total_unpred * sizeof(double);
2344
2345        int * result_type = (int *) malloc(num_elements * sizeof(int));
2346        decode(comp_data_pos, num_elements, root, result_type);
2347        SZ_ReleaseHuffman(huffmanTree);
2348       
2349        int intvRadius = exe_params->intvRadius;
2350       
2351        int * type;
2352        double * data_pos = *data;
2353        size_t offset_x, offset_y, offset_z;
2354        size_t current_blockcount_x, current_blockcount_y, current_blockcount_z;
2355        size_t cur_unpred_count;
2356        unsigned char * indicator_pos = indicator;
2357        if(use_mean){
2358                // type = result_type;
2359
2360                // for(size_t i=0; i<num_x; i++){
2361                //      for(size_t j=0; j<num_y; j++){
2362                //              for(size_t k=0; k<num_z; k++){
2363                //                      offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
2364                //                      offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2365                //                      offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
2366                //                      data_pos = *data + offset_x * dim0_offset + offset_y * dim1_offset + offset_z;
2367
2368                //                      current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
2369                //                      current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2370                //                      current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
2371
2372                //                      // type_offset = offset_x * dim0_offset +  offset_y * current_blockcount_x * dim1_offset + offset_z * current_blockcount_x * current_blockcount_y;
2373                //                      // type = result_type + type_offset;
2374                //                      size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2375                //                      // index = i * num_y * num_z + j * num_z + k;
2376
2377                //                      // printf("i j k: %ld %ld %ld\toffset: %ld %ld %ld\tindicator: %ld\n", i, j, k, offset_x, offset_y, offset_z, indicator[index]);
2378                //                      if(*indicator_pos){
2379                //                              // decompress by SZ
2380                //                              // cur_unpred_count = decompressDataSeries_double_3D_blocked_nonblock_pred(data_pos, r1, r2, r3, current_blockcount_x, current_blockcount_y, current_blockcount_z, i, j, k, realPrecision, type, unpred_data);
2381                //                              double * block_data_pos = data_pos;
2382                //                              double pred;
2383                //                              size_t index = 0;
2384                //                              int type_;
2385                //                              // d111 is current data
2386                //                              size_t unpredictable_count = 0;
2387                //                              double d000, d001, d010, d011, d100, d101, d110;
2388                //                              for(size_t ii=0; ii<current_blockcount_x; ii++){
2389                //                                      for(size_t jj=0; jj<current_blockcount_y; jj++){
2390                //                                              for(size_t kk=0; kk<current_blockcount_z; kk++){
2391                //                                                      type_ = type[index];
2392                //                                                      if(type_ == intvRadius){
2393                //                                                              *block_data_pos = mean;
2394                //                                                      }
2395                //                                                      else if(type_ == 0){
2396                //                                                              *block_data_pos = unpred_data[unpredictable_count ++];
2397                //                                                      }
2398                //                                                      else{
2399                //                                                              d000 = d001 = d010 = d011 = d100 = d101 = d110 = 1;
2400                //                                                              if(i == 0 && ii == 0){
2401                //                                                                      d000 = d001 = d010 = d011 = 0;
2402                //                                                              }
2403                //                                                              if(j == 0 && jj == 0){
2404                //                                                                      d000 = d001 = d100 = d101 = 0;
2405                //                                                              }
2406                //                                                              if(k == 0 && kk == 0){
2407                //                                                                      d000 = d010 = d100 = d110 = 0;
2408                //                                                              }
2409                //                                                              if(d000){
2410                //                                                                      d000 = block_data_pos[- dim0_offset - dim1_offset - 1];
2411                //                                                              }
2412                //                                                              if(d001){
2413                //                                                                      d001 = block_data_pos[- dim0_offset - dim1_offset];
2414                //                                                              }
2415                //                                                              if(d010){
2416                //                                                                      d010 = block_data_pos[- dim0_offset - 1];
2417                //                                                              }
2418                //                                                              if(d011){
2419                //                                                                      d011 = block_data_pos[- dim0_offset];
2420                //                                                              }
2421                //                                                              if(d100){
2422                //                                                                      d100 = block_data_pos[- dim1_offset - 1];
2423                //                                                              }
2424                //                                                              if(d101){
2425                //                                                                      d101 = block_data_pos[- dim1_offset];
2426                //                                                              }
2427                //                                                              if(d110){
2428                //                                                                      d110 = block_data_pos[- 1];
2429                //                                                              }
2430                //                                                              if(type_ < intvRadius) type_ += 1;
2431                //                                                              pred = d110 + d101 + d011 - d100 - d010 - d001 + d000;
2432                //                                                              *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2433                //                                                      }
2434                //                                                      index ++;
2435                //                                                      block_data_pos ++;
2436                //                                              }
2437                //                                              block_data_pos += dim1_offset - current_blockcount_z;
2438                //                                      }
2439                //                                      block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2440                //                              }
2441                //                              cur_unpred_count = unpredictable_count;
2442                //                      }
2443                //                      else{
2444                //                              // decompress by regression
2445                //                              {
2446                //                                      //restore regression coefficients
2447                //                                      double pred;
2448                //                                      int type_;
2449                //                                      for(int e=0; e<4; e++){
2450                //                                              // if(i == 0 && j == 0 && k == 19){
2451                //                                              //      printf("~\n");
2452                //                                              // }
2453                //                                              type_ = coeff_type[e][coeff_index];
2454                //                                              if (type_ != 0){
2455                //                                                      pred = last_coefficients[e];
2456                //                                                      last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2457                //                                              }
2458                //                                              else{
2459                //                                                      last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2460                //                                                      coeff_unpred_data_count[e] ++;
2461                //                                              }
2462                //                                              if(fabs(last_coefficients[e]) > 10000){
2463                //                                                      printf("%d %d %d-%d: pred %.4f type %d precision %.4g last_coefficients %.4g\n", i, j, k, e, pred, type_, precision[e], last_coefficients[e]);
2464                //                                                      exit(0);
2465                //                                              }
2466                //                                      }
2467                //                                      coeff_index ++;
2468                //                              }
2469                //                              {
2470                //                                      double * block_data_pos = data_pos;
2471                //                                      double pred;
2472                //                                      int type_;
2473                //                                      size_t index = 0;
2474                //                                      size_t unpredictable_count = 0;
2475                //                                      for(size_t ii=0; ii<current_blockcount_x; ii++){
2476                //                                              for(size_t jj=0; jj<current_blockcount_y; jj++){
2477                //                                                      for(size_t kk=0; kk<current_blockcount_z; kk++){
2478                //                                                              if(block_data_pos - (*data) == 19470788){
2479                //                                                                      printf("dec stop\n");
2480                //                                                              }
2481
2482                //                                                              type_ = type[index];
2483                //                                                              if (type_ != 0){
2484                //                                                                      pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
2485                //                                                                      *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2486                //                                                              }
2487                //                                                              else{
2488                //                                                                      *block_data_pos = unpred_data[unpredictable_count ++];
2489                //                                                              }
2490                //                                                              index ++;       
2491                //                                                              block_data_pos ++;
2492                //                                                      }
2493                //                                                      block_data_pos += dim1_offset - current_blockcount_z;
2494                //                                              }
2495                //                                              block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2496                //                                      }
2497                //                                      cur_unpred_count = unpredictable_count;
2498                //                              }
2499                //                      }
2500
2501                //                      type += current_block_elements;
2502                //                      indicator_pos ++;
2503                //                      unpred_data += cur_unpred_count;
2504                //                      // decomp_unpred += cur_unpred_count;
2505                //                      // printf("block comp done, data_offset from %ld to %ld: diff %ld\n", *data, data_pos, data_pos - *data);
2506                //                      // fflush(stdout);
2507                //              }
2508                //      }
2509                // }
2510
2511                type = result_type;
2512                // i == 0
2513                {
2514                        // j == 0
2515                        {
2516                                // k == 0
2517                                {
2518                                        data_pos = *data;
2519
2520                                        current_blockcount_x = early_blockcount_x;
2521                                        current_blockcount_y = early_blockcount_y;
2522                                        current_blockcount_z = early_blockcount_z;
2523                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2524                                        if(*indicator_pos){
2525                                                // decompress by SZ
2526                                                double * block_data_pos = data_pos;
2527                                                double pred;
2528                                                size_t index = 0;
2529                                                int type_;
2530                                                size_t unpredictable_count = 0;
2531                                                // ii == 0
2532                                                {
2533                                                        // jj == 0
2534                                                        {
2535                                                                {
2536                                                                        // kk == 0
2537                                                                        type_ = type[index];
2538                                                                        if(type_ == intvRadius){
2539                                                                                *block_data_pos = mean;
2540                                                                        }
2541                                                                        else if(type_ == 0){
2542                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2543                                                                        }
2544                                                                        else{
2545                                                                                if(type_ < intvRadius) type_ += 1;
2546                                                                                pred = 0;
2547                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2548                                                                        }
2549                                                                        index ++;
2550                                                                        block_data_pos ++;
2551                                                                }
2552                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
2553                                                                        type_ = type[index];
2554                                                                        if(type_ == intvRadius){
2555                                                                                *block_data_pos = mean;
2556                                                                        }
2557                                                                        else if(type_ == 0){
2558                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2559                                                                        }
2560                                                                        else{
2561                                                                                if(type_ < intvRadius) type_ += 1;
2562                                                                                pred = block_data_pos[- 1];
2563                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2564                                                                        }
2565                                                                        index ++;
2566                                                                        block_data_pos ++;
2567                                                                }
2568                                                                block_data_pos += dim1_offset - current_blockcount_z;
2569                                                        }
2570                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
2571                                                                {
2572                                                                        // kk == 0
2573                                                                        type_ = type[index];
2574                                                                        if(type_ == intvRadius){
2575                                                                                *block_data_pos = mean;
2576                                                                        }
2577                                                                        else if(type_ == 0){
2578                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2579                                                                        }
2580                                                                        else{
2581                                                                                if(type_ < intvRadius) type_ += 1;
2582                                                                                pred = block_data_pos[- dim1_offset];
2583                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2584                                                                        }
2585                                                                        index ++;
2586                                                                        block_data_pos ++;
2587                                                                }
2588                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
2589                                                                        type_ = type[index];
2590                                                                        if(type_ == intvRadius){
2591                                                                                *block_data_pos = mean;
2592                                                                        }
2593                                                                        else if(type_ == 0){
2594                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2595                                                                        }
2596                                                                        else{
2597                                                                                if(type_ < intvRadius) type_ += 1;
2598                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
2599                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2600                                                                        }
2601                                                                        index ++;
2602                                                                        block_data_pos ++;
2603                                                                }
2604                                                                block_data_pos += dim1_offset - current_blockcount_z;
2605                                                        }
2606                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;                                             
2607                                                }
2608                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
2609                                                        // jj == 0
2610                                                        {
2611                                                                {
2612                                                                        // kk == 0
2613                                                                        type_ = type[index];
2614                                                                        if(type_ == intvRadius){
2615                                                                                *block_data_pos = mean;
2616                                                                        }
2617                                                                        else if(type_ == 0){
2618                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2619                                                                        }
2620                                                                        else{
2621                                                                                if(type_ < intvRadius) type_ += 1;
2622                                                                                pred = block_data_pos[- dim0_offset];
2623                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2624                                                                        }
2625                                                                        index ++;
2626                                                                        block_data_pos ++;
2627                                                                }
2628                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
2629                                                                        type_ = type[index];
2630                                                                        if(type_ == intvRadius){
2631                                                                                *block_data_pos = mean;
2632                                                                        }
2633                                                                        else if(type_ == 0){
2634                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2635                                                                        }
2636                                                                        else{
2637                                                                                if(type_ < intvRadius) type_ += 1;
2638                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
2639                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2640                                                                        }
2641                                                                        index ++;
2642                                                                        block_data_pos ++;
2643                                                                }
2644                                                                block_data_pos += dim1_offset - current_blockcount_z;
2645                                                        }
2646                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
2647                                                                {
2648                                                                        // kk == 0
2649                                                                        type_ = type[index];
2650                                                                        if(type_ == intvRadius){
2651                                                                                *block_data_pos = mean;
2652                                                                        }
2653                                                                        else if(type_ == 0){
2654                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2655                                                                        }
2656                                                                        else{
2657                                                                                if(type_ < intvRadius) type_ += 1;
2658                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
2659                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2660                                                                        }
2661                                                                        index ++;
2662                                                                        block_data_pos ++;
2663                                                                }
2664                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
2665                                                                        type_ = type[index];
2666                                                                        if(type_ == intvRadius){
2667                                                                                *block_data_pos = mean;
2668                                                                        }
2669                                                                        else if(type_ == 0){
2670                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2671                                                                        }
2672                                                                        else{
2673                                                                                if(type_ < intvRadius) type_ += 1;
2674                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
2675                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2676                                                                        }
2677                                                                        index ++;
2678                                                                        block_data_pos ++;
2679                                                                }
2680                                                                block_data_pos += dim1_offset - current_blockcount_z;
2681                                                        }
2682                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2683                                                }
2684                                                cur_unpred_count = unpredictable_count;
2685                                        }
2686                                        else{
2687                                                // decompress by regression
2688                                                {
2689                                                        //restore regression coefficients
2690                                                        double pred;
2691                                                        int type_;
2692                                                        for(int e=0; e<4; e++){
2693                                                                type_ = coeff_type[e][coeff_index];
2694                                                                if (type_ != 0){
2695                                                                        pred = last_coefficients[e];
2696                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2697                                                                }
2698                                                                else{
2699                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2700                                                                        coeff_unpred_data_count[e] ++;
2701                                                                }
2702                                                        }
2703                                                        coeff_index ++;
2704                                                }
2705                                                {
2706                                                        double * block_data_pos = data_pos;
2707                                                        double pred;
2708                                                        int type_;
2709                                                        size_t index = 0;
2710                                                        size_t unpredictable_count = 0;
2711                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
2712                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
2713                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
2714                                                                                type_ = type[index];
2715                                                                                if (type_ != 0){
2716                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
2717                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2718                                                                                }
2719                                                                                else{
2720                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
2721                                                                                }
2722                                                                                index ++;       
2723                                                                                block_data_pos ++;
2724                                                                        }
2725                                                                        block_data_pos += dim1_offset - current_blockcount_z;
2726                                                                }
2727                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2728                                                        }
2729                                                        cur_unpred_count = unpredictable_count;
2730                                                }
2731                                        }
2732                                        indicator_pos ++;
2733                                        type += current_block_elements;
2734                                        unpred_data += cur_unpred_count;
2735                                } // end k == 0
2736                                // i == 0 j == 0 k != 0
2737                                for(size_t k=1; k<num_z; k++){
2738                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
2739                                        data_pos = *data + offset_z;
2740
2741                                        current_blockcount_x = early_blockcount_x;
2742                                        current_blockcount_y = early_blockcount_y;
2743                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
2744
2745                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2746                                        if(*indicator_pos){
2747                                                // decompress by SZ
2748                                                double * block_data_pos = data_pos;
2749                                                double pred;
2750                                                size_t index = 0;
2751                                                int type_;
2752                                                size_t unpredictable_count = 0;
2753                                                // ii == 0
2754                                                {
2755                                                        // jj == 0
2756                                                        {
2757                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
2758                                                                        type_ = type[index];
2759                                                                        if(type_ == intvRadius){
2760                                                                                *block_data_pos = mean;
2761                                                                        }
2762                                                                        else if(type_ == 0){
2763                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2764                                                                        }
2765                                                                        else{
2766                                                                                if(type_ < intvRadius) type_ += 1;
2767                                                                                pred = block_data_pos[- 1];
2768                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2769                                                                        }
2770                                                                        index ++;
2771                                                                        block_data_pos ++;
2772                                                                }
2773                                                                block_data_pos += dim1_offset - current_blockcount_z;
2774                                                        }
2775                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
2776                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
2777                                                                        type_ = type[index];
2778                                                                        if(type_ == intvRadius){
2779                                                                                *block_data_pos = mean;
2780                                                                        }
2781                                                                        else if(type_ == 0){
2782                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2783                                                                        }
2784                                                                        else{
2785                                                                                if(type_ < intvRadius) type_ += 1;
2786                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
2787                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2788                                                                        }
2789                                                                        index ++;
2790                                                                        block_data_pos ++;
2791                                                                }
2792                                                                block_data_pos += dim1_offset - current_blockcount_z;
2793                                                        }
2794                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2795                                                }
2796                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
2797                                                        // jj == 0
2798                                                        {
2799                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
2800                                                                        type_ = type[index];
2801                                                                        if(type_ == intvRadius){
2802                                                                                *block_data_pos = mean;
2803                                                                        }
2804                                                                        else if(type_ == 0){
2805                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2806                                                                        }
2807                                                                        else{
2808                                                                                if(type_ < intvRadius) type_ += 1;
2809                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
2810                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2811                                                                        }
2812                                                                        index ++;
2813                                                                        block_data_pos ++;
2814                                                                }
2815                                                                block_data_pos += dim1_offset - current_blockcount_z;
2816                                                        }
2817                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
2818                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
2819                                                                        type_ = type[index];
2820                                                                        if(type_ == intvRadius){
2821                                                                                *block_data_pos = mean;
2822                                                                        }
2823                                                                        else if(type_ == 0){
2824                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2825                                                                        }
2826                                                                        else{
2827                                                                                if(type_ < intvRadius) type_ += 1;
2828                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
2829                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2830                                                                        }
2831                                                                        index ++;
2832                                                                        block_data_pos ++;
2833                                                                }
2834                                                                block_data_pos += dim1_offset - current_blockcount_z;
2835                                                        }
2836                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2837                                                }
2838                                                cur_unpred_count = unpredictable_count;
2839                                        }
2840                                        else{
2841                                                // decompress by regression
2842                                                {
2843                                                        //restore regression coefficients
2844                                                        double pred;
2845                                                        int type_;
2846                                                        for(int e=0; e<4; e++){
2847                                                                type_ = coeff_type[e][coeff_index];
2848                                                                if (type_ != 0){
2849                                                                        pred = last_coefficients[e];
2850                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
2851                                                                }
2852                                                                else{
2853                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
2854                                                                        coeff_unpred_data_count[e] ++;
2855                                                                }
2856                                                        }
2857                                                        coeff_index ++;
2858                                                }
2859                                                {
2860                                                        double * block_data_pos = data_pos;
2861                                                        double pred;
2862                                                        int type_;
2863                                                        size_t index = 0;
2864                                                        size_t unpredictable_count = 0;
2865                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
2866                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
2867                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
2868                                                                                type_ = type[index];
2869                                                                                if (type_ != 0){
2870                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
2871                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2872                                                                                }
2873                                                                                else{
2874                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
2875                                                                                }
2876                                                                                index ++;       
2877                                                                                block_data_pos ++;
2878                                                                        }
2879                                                                        block_data_pos += dim1_offset - current_blockcount_z;
2880                                                                }
2881                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2882                                                        }
2883                                                        cur_unpred_count = unpredictable_count;
2884                                                }
2885                                        }
2886                                        indicator_pos ++;
2887                                        type += current_block_elements;
2888                                        unpred_data += cur_unpred_count;
2889                                }
2890                        }// end j==0
2891                        for(size_t j=1; j<num_y; j++){
2892                                // k == 0
2893                                {
2894                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
2895                                        data_pos = *data + offset_y * dim1_offset;
2896
2897                                        current_blockcount_x = early_blockcount_x;
2898                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
2899                                        current_blockcount_z = early_blockcount_z;
2900                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
2901                                        if(*indicator_pos){
2902                                                // decompress by SZ
2903                                                double * block_data_pos = data_pos;
2904                                                double pred;
2905                                                size_t index = 0;
2906                                                int type_;
2907                                                size_t unpredictable_count = 0;
2908                                                // ii == 0
2909                                                {
2910                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
2911                                                                {
2912                                                                        // kk == 0
2913                                                                        type_ = type[index];
2914                                                                        if(type_ == intvRadius){
2915                                                                                *block_data_pos = mean;
2916                                                                        }
2917                                                                        else if(type_ == 0){
2918                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2919                                                                        }
2920                                                                        else{
2921                                                                                if(type_ < intvRadius) type_ += 1;
2922                                                                                pred = block_data_pos[- dim1_offset];
2923                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2924                                                                        }
2925                                                                        index ++;
2926                                                                        block_data_pos ++;
2927                                                                }
2928                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
2929                                                                        type_ = type[index];
2930                                                                        if(type_ == intvRadius){
2931                                                                                *block_data_pos = mean;
2932                                                                        }
2933                                                                        else if(type_ == 0){
2934                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2935                                                                        }
2936                                                                        else{
2937                                                                                if(type_ < intvRadius) type_ += 1;
2938                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
2939                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2940                                                                        }
2941                                                                        index ++;
2942                                                                        block_data_pos ++;
2943                                                                }
2944                                                                block_data_pos += dim1_offset - current_blockcount_z;
2945                                                        }
2946                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2947                                                }
2948                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
2949                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
2950                                                                {
2951                                                                        // kk == 0
2952                                                                        type_ = type[index];
2953                                                                        if(type_ == intvRadius){
2954                                                                                *block_data_pos = mean;
2955                                                                        }
2956                                                                        else if(type_ == 0){
2957                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2958                                                                        }
2959                                                                        else{
2960                                                                                if(type_ < intvRadius) type_ += 1;
2961                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
2962                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2963                                                                        }
2964                                                                        index ++;
2965                                                                        block_data_pos ++;
2966                                                                }
2967                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
2968                                                                        type_ = type[index];
2969                                                                        if(type_ == intvRadius){
2970                                                                                *block_data_pos = mean;
2971                                                                        }
2972                                                                        else if(type_ == 0){
2973                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
2974                                                                        }
2975                                                                        else{
2976                                                                                if(type_ < intvRadius) type_ += 1;
2977                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
2978                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
2979                                                                        }
2980                                                                        index ++;
2981                                                                        block_data_pos ++;
2982                                                                }
2983                                                                block_data_pos += dim1_offset - current_blockcount_z;
2984                                                        }
2985                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
2986                                                }
2987                                                cur_unpred_count = unpredictable_count;
2988                                        }
2989                                        else{
2990                                                // decompress by regression
2991                                                {
2992                                                        //restore regression coefficients
2993                                                        double pred;
2994                                                        int type_;
2995                                                        for(int e=0; e<4; e++){
2996                                                                type_ = coeff_type[e][coeff_index];
2997                                                                if (type_ != 0){
2998                                                                        pred = last_coefficients[e];
2999                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3000                                                                }
3001                                                                else{
3002                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3003                                                                        coeff_unpred_data_count[e] ++;
3004                                                                }
3005                                                        }
3006                                                        coeff_index ++;
3007                                                }
3008                                                {
3009                                                        double * block_data_pos = data_pos;
3010                                                        double pred;
3011                                                        int type_;
3012                                                        size_t index = 0;
3013                                                        size_t unpredictable_count = 0;
3014                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3015                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3016                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3017                                                                                type_ = type[index];
3018                                                                                if (type_ != 0){
3019                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3020                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3021                                                                                }
3022                                                                                else{
3023                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3024                                                                                }
3025                                                                                index ++;       
3026                                                                                block_data_pos ++;
3027                                                                        }
3028                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3029                                                                }
3030                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3031                                                        }
3032                                                        cur_unpred_count = unpredictable_count;
3033                                                }
3034                                        }
3035                                        indicator_pos ++;
3036                                        type += current_block_elements;
3037                                        unpred_data += cur_unpred_count;
3038                                } // end k == 0
3039                                for(size_t k=1; k<num_z; k++){
3040                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3041                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3042                                        data_pos = *data + offset_y * dim1_offset + offset_z;
3043
3044                                        current_blockcount_x = early_blockcount_x;
3045                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3046                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3047
3048                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3049                                        if(*indicator_pos){
3050                                                // decompress by SZ
3051                                                double * block_data_pos = data_pos;
3052                                                double pred;
3053                                                size_t index = 0;
3054                                                int type_;
3055                                                size_t unpredictable_count = 0;
3056                                                // ii == 0
3057                                                {
3058                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
3059                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3060                                                                        type_ = type[index];
3061                                                                        if(type_ == intvRadius){
3062                                                                                *block_data_pos = mean;
3063                                                                        }
3064                                                                        else if(type_ == 0){
3065                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3066                                                                        }
3067                                                                        else{
3068                                                                                if(type_ < intvRadius) type_ += 1;
3069                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3070                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3071                                                                        }
3072                                                                        index ++;
3073                                                                        block_data_pos ++;
3074                                                                }
3075                                                                block_data_pos += dim1_offset - current_blockcount_z;
3076                                                        }
3077                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3078                                                }
3079                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
3080                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
3081                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3082                                                                        type_ = type[index];
3083                                                                        if(type_ == intvRadius){
3084                                                                                *block_data_pos = mean;
3085                                                                        }
3086                                                                        else if(type_ == 0){
3087                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3088                                                                        }
3089                                                                        else{
3090                                                                                if(type_ < intvRadius) type_ += 1;
3091                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3092                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3093                                                                        }
3094                                                                        index ++;
3095                                                                        block_data_pos ++;
3096                                                                }
3097                                                                block_data_pos += dim1_offset - current_blockcount_z;
3098                                                        }
3099                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3100                                                }
3101                                                cur_unpred_count = unpredictable_count;
3102                                        }
3103                                        else{
3104                                                // decompress by regression
3105                                                {
3106                                                        //restore regression coefficients
3107                                                        double pred;
3108                                                        int type_;
3109                                                        for(int e=0; e<4; e++){
3110                                                                type_ = coeff_type[e][coeff_index];
3111                                                                if (type_ != 0){
3112                                                                        pred = last_coefficients[e];
3113                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3114                                                                }
3115                                                                else{
3116                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3117                                                                        coeff_unpred_data_count[e] ++;
3118                                                                }
3119                                                        }
3120                                                        coeff_index ++;
3121                                                }
3122                                                {
3123                                                        double * block_data_pos = data_pos;
3124                                                        double pred;
3125                                                        int type_;
3126                                                        size_t index = 0;
3127                                                        size_t unpredictable_count = 0;
3128                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3129                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3130                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3131                                                                                type_ = type[index];
3132                                                                                if (type_ != 0){
3133                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3134                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3135                                                                                }
3136                                                                                else{
3137                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3138                                                                                }
3139                                                                                index ++;       
3140                                                                                block_data_pos ++;
3141                                                                        }
3142                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3143                                                                }
3144                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3145                                                        }
3146                                                        cur_unpred_count = unpredictable_count;
3147                                                }
3148                                        }
3149                                        indicator_pos ++;
3150                                        type += current_block_elements;
3151                                        unpred_data += cur_unpred_count;
3152                                }
3153                        }
3154                } // end i==0
3155                for(size_t i=1; i<num_x; i++){
3156                        // j == 0
3157                        {
3158                                // k == 0
3159                                {
3160                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3161                                        data_pos = *data + offset_x * dim0_offset;
3162
3163                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3164                                        current_blockcount_y = early_blockcount_y;
3165                                        current_blockcount_z = early_blockcount_z;
3166                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3167                                        if(*indicator_pos){
3168                                                // decompress by SZ
3169                                                double * block_data_pos = data_pos;
3170                                                double pred;
3171                                                size_t index = 0;
3172                                                int type_;
3173                                                size_t unpredictable_count = 0;
3174                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
3175                                                        // jj == 0
3176                                                        {
3177                                                                {
3178                                                                        // kk == 0
3179                                                                        type_ = type[index];
3180                                                                        if(type_ == intvRadius){
3181                                                                                *block_data_pos = mean;
3182                                                                        }
3183                                                                        else if(type_ == 0){
3184                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3185                                                                        }
3186                                                                        else{
3187                                                                                if(type_ < intvRadius) type_ += 1;
3188                                                                                pred = block_data_pos[- dim0_offset];
3189                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3190                                                                        }
3191                                                                        index ++;
3192                                                                        block_data_pos ++;
3193                                                                }
3194                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3195                                                                        type_ = type[index];
3196                                                                        if(type_ == intvRadius){
3197                                                                                *block_data_pos = mean;
3198                                                                        }
3199                                                                        else if(type_ == 0){
3200                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3201                                                                        }
3202                                                                        else{
3203                                                                                if(type_ < intvRadius) type_ += 1;
3204                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3205                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3206                                                                        }
3207                                                                        index ++;
3208                                                                        block_data_pos ++;
3209                                                                }
3210                                                                block_data_pos += dim1_offset - current_blockcount_z;
3211                                                        }
3212                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
3213                                                                {
3214                                                                        // kk == 0
3215                                                                        type_ = type[index];
3216                                                                        if(type_ == intvRadius){
3217                                                                                *block_data_pos = mean;
3218                                                                        }
3219                                                                        else if(type_ == 0){
3220                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3221                                                                        }
3222                                                                        else{
3223                                                                                if(type_ < intvRadius) type_ += 1;
3224                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
3225                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3226                                                                        }
3227                                                                        index ++;
3228                                                                        block_data_pos ++;
3229                                                                }
3230                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3231                                                                        type_ = type[index];
3232                                                                        if(type_ == intvRadius){
3233                                                                                *block_data_pos = mean;
3234                                                                        }
3235                                                                        else if(type_ == 0){
3236                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3237                                                                        }
3238                                                                        else{
3239                                                                                if(type_ < intvRadius) type_ += 1;
3240                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3241                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3242                                                                        }
3243                                                                        index ++;
3244                                                                        block_data_pos ++;
3245                                                                }
3246                                                                block_data_pos += dim1_offset - current_blockcount_z;
3247                                                        }
3248                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3249                                                }
3250                                                cur_unpred_count = unpredictable_count;
3251                                        }
3252                                        else{
3253                                                // decompress by regression
3254                                                {
3255                                                        //restore regression coefficients
3256                                                        double pred;
3257                                                        int type_;
3258                                                        for(int e=0; e<4; e++){
3259                                                                type_ = coeff_type[e][coeff_index];
3260                                                                if (type_ != 0){
3261                                                                        pred = last_coefficients[e];
3262                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3263                                                                }
3264                                                                else{
3265                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3266                                                                        coeff_unpred_data_count[e] ++;
3267                                                                }
3268                                                        }
3269                                                        coeff_index ++;
3270                                                }
3271                                                {
3272                                                        double * block_data_pos = data_pos;
3273                                                        double pred;
3274                                                        int type_;
3275                                                        size_t index = 0;
3276                                                        size_t unpredictable_count = 0;
3277                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3278                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3279                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3280                                                                                type_ = type[index];
3281                                                                                if (type_ != 0){
3282                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3283                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3284                                                                                }
3285                                                                                else{
3286                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3287                                                                                }
3288                                                                                index ++;       
3289                                                                                block_data_pos ++;
3290                                                                        }
3291                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3292                                                                }
3293                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3294                                                        }
3295                                                        cur_unpred_count = unpredictable_count;
3296                                                }
3297                                        }
3298                                        indicator_pos ++;
3299                                        type += current_block_elements;
3300                                        unpred_data += cur_unpred_count;
3301                                } // end k == 0
3302                                for(size_t k=1; k<num_z; k++){
3303                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3304                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3305                                        data_pos = *data + offset_x * dim0_offset + offset_z;
3306
3307                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3308                                        current_blockcount_y = early_blockcount_y;
3309                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3310                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3311                                        if(*indicator_pos){
3312                                                // decompress by SZ
3313                                                double * block_data_pos = data_pos;
3314                                                double pred;
3315                                                size_t index = 0;
3316                                                int type_;
3317                                                size_t unpredictable_count = 0;
3318                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
3319                                                        // jj == 0
3320                                                        {
3321                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3322                                                                        type_ = type[index];
3323                                                                        if(type_ == intvRadius){
3324                                                                                *block_data_pos = mean;
3325                                                                        }
3326                                                                        else if(type_ == 0){
3327                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3328                                                                        }
3329                                                                        else{
3330                                                                                if(type_ < intvRadius) type_ += 1;
3331                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3332                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3333                                                                        }
3334                                                                        index ++;
3335                                                                        block_data_pos ++;
3336                                                                }
3337                                                                block_data_pos += dim1_offset - current_blockcount_z;
3338                                                        }
3339                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
3340                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3341                                                                        type_ = type[index];
3342                                                                        if(type_ == intvRadius){
3343                                                                                *block_data_pos = mean;
3344                                                                        }
3345                                                                        else if(type_ == 0){
3346                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3347                                                                        }
3348                                                                        else{
3349                                                                                if(type_ < intvRadius) type_ += 1;
3350                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3351                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3352                                                                        }
3353                                                                        index ++;
3354                                                                        block_data_pos ++;
3355                                                                }
3356                                                                block_data_pos += dim1_offset - current_blockcount_z;
3357                                                        }
3358                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3359                                                }
3360                                                cur_unpred_count = unpredictable_count;
3361                                        }
3362                                        else{
3363                                                // decompress by regression
3364                                                {
3365                                                        //restore regression coefficients
3366                                                        double pred;
3367                                                        int type_;
3368                                                        for(int e=0; e<4; e++){
3369                                                                type_ = coeff_type[e][coeff_index];
3370                                                                if (type_ != 0){
3371                                                                        pred = last_coefficients[e];
3372                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3373                                                                }
3374                                                                else{
3375                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3376                                                                        coeff_unpred_data_count[e] ++;
3377                                                                }
3378                                                        }
3379                                                        coeff_index ++;
3380                                                }
3381                                                {
3382                                                        double * block_data_pos = data_pos;
3383                                                        double pred;
3384                                                        int type_;
3385                                                        size_t index = 0;
3386                                                        size_t unpredictable_count = 0;
3387                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3388                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3389                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3390                                                                                type_ = type[index];
3391                                                                                if (type_ != 0){
3392                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3393                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3394                                                                                }
3395                                                                                else{
3396                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3397                                                                                }
3398                                                                                index ++;       
3399                                                                                block_data_pos ++;
3400                                                                        }
3401                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3402                                                                }
3403                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3404                                                        }
3405                                                        cur_unpred_count = unpredictable_count;
3406                                                }
3407                                        }
3408                                        indicator_pos ++;
3409                                        type += current_block_elements;
3410                                        unpred_data += cur_unpred_count;
3411                                }
3412                        }// end j = 0
3413                        for(size_t j=1; j<num_y; j++){
3414                                // k == 0
3415                                {
3416                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3417                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3418                                        data_pos = *data + offset_x * dim0_offset + offset_y * dim1_offset;
3419
3420                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3421                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3422                                        current_blockcount_z = early_blockcount_z;
3423                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3424                                        if(*indicator_pos){
3425                                                // decompress by SZ
3426                                                double * block_data_pos = data_pos;
3427                                                double pred;
3428                                                size_t index = 0;
3429                                                int type_;
3430                                                size_t unpredictable_count = 0;
3431                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
3432                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
3433                                                                {
3434                                                                        // kk == 0
3435                                                                        type_ = type[index];
3436                                                                        if(type_ == intvRadius){
3437                                                                                *block_data_pos = mean;
3438                                                                        }
3439                                                                        else if(type_ == 0){
3440                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3441                                                                        }
3442                                                                        else{
3443                                                                                if(type_ < intvRadius) type_ += 1;
3444                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
3445                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3446                                                                        }
3447                                                                        index ++;
3448                                                                        block_data_pos ++;
3449                                                                }
3450                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3451                                                                        type_ = type[index];
3452                                                                        if(type_ == intvRadius){
3453                                                                                *block_data_pos = mean;
3454                                                                        }
3455                                                                        else if(type_ == 0){
3456                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3457                                                                        }
3458                                                                        else{
3459                                                                                if(type_ < intvRadius) type_ += 1;
3460                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3461                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3462                                                                        }
3463                                                                        index ++;
3464                                                                        block_data_pos ++;
3465                                                                }
3466                                                                block_data_pos += dim1_offset - current_blockcount_z;
3467                                                        }
3468                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3469                                                }
3470                                                cur_unpred_count = unpredictable_count;
3471                                        }
3472                                        else{
3473                                                // decompress by regression
3474                                                {
3475                                                        //restore regression coefficients
3476                                                        double pred;
3477                                                        int type_;
3478                                                        for(int e=0; e<4; e++){
3479                                                                type_ = coeff_type[e][coeff_index];
3480                                                                if (type_ != 0){
3481                                                                        pred = last_coefficients[e];
3482                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3483                                                                }
3484                                                                else{
3485                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3486                                                                        coeff_unpred_data_count[e] ++;
3487                                                                }
3488                                                        }
3489                                                        coeff_index ++;
3490                                                }
3491                                                {
3492                                                        double * block_data_pos = data_pos;
3493                                                        double pred;
3494                                                        int type_;
3495                                                        size_t index = 0;
3496                                                        size_t unpredictable_count = 0;
3497                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3498                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3499                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3500                                                                                type_ = type[index];
3501                                                                                if (type_ != 0){
3502                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3503                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3504                                                                                }
3505                                                                                else{
3506                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3507                                                                                }
3508                                                                                index ++;       
3509                                                                                block_data_pos ++;
3510                                                                        }
3511                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3512                                                                }
3513                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3514                                                        }
3515                                                        cur_unpred_count = unpredictable_count;
3516                                                }
3517                                        }
3518                                        indicator_pos ++;
3519                                        type += current_block_elements;
3520                                        unpred_data += cur_unpred_count;
3521                                } // end k == 0
3522                                for(size_t k=1; k<num_z; k++){
3523                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
3524                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3525                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3526                                        data_pos = *data + offset_x * dim0_offset + offset_y * dim1_offset + offset_z;
3527
3528                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
3529                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3530                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3531
3532                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3533                                        if(*indicator_pos){
3534                                                // decompress by SZ
3535                                                double * block_data_pos = data_pos;
3536                                                double pred;
3537                                                size_t index = 0;
3538                                                int type_;
3539                                                size_t unpredictable_count = 0;
3540                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
3541                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
3542                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3543                                                                        type_ = type[index];
3544                                                                        if(type_ == intvRadius){
3545                                                                                *block_data_pos = mean;
3546                                                                        }
3547                                                                        else if(type_ == 0){
3548                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3549                                                                        }
3550                                                                        else{
3551                                                                                if(type_ < intvRadius) type_ += 1;
3552                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3553                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3554                                                                        }
3555                                                                        index ++;
3556                                                                        block_data_pos ++;
3557                                                                }
3558                                                                block_data_pos += dim1_offset - current_blockcount_z;
3559                                                        }
3560                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3561                                                }
3562                                                cur_unpred_count = unpredictable_count;
3563                                        }
3564                                        else{
3565                                                // decompress by regression
3566                                                {
3567                                                        //restore regression coefficients
3568                                                        double pred;
3569                                                        int type_;
3570                                                        for(int e=0; e<4; e++){
3571                                                                type_ = coeff_type[e][coeff_index];
3572                                                                if (type_ != 0){
3573                                                                        pred = last_coefficients[e];
3574                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3575                                                                }
3576                                                                else{
3577                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3578                                                                        coeff_unpred_data_count[e] ++;
3579                                                                }
3580                                                        }
3581                                                        coeff_index ++;
3582                                                }
3583                                                {
3584                                                        double * block_data_pos = data_pos;
3585                                                        double pred;
3586                                                        int type_;
3587                                                        size_t index = 0;
3588                                                        size_t unpredictable_count = 0;
3589                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3590                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3591                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3592                                                                                type_ = type[index];
3593                                                                                if (type_ != 0){
3594                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3595                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3596                                                                                }
3597                                                                                else{
3598                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3599                                                                                }
3600                                                                                index ++;       
3601                                                                                block_data_pos ++;
3602                                                                        }
3603                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3604                                                                }
3605                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3606                                                        }
3607                                                        cur_unpred_count = unpredictable_count;
3608                                                }
3609                                        }
3610                                        indicator_pos ++;
3611                                        type += current_block_elements;
3612                                        unpred_data += cur_unpred_count;
3613                                }
3614                        }
3615                }
3616        }
3617        else{
3618                type = result_type;
3619                // i == 0
3620                {
3621                        // j == 0
3622                        {
3623                                // k == 0
3624                                {
3625                                        data_pos = *data;
3626
3627                                        current_blockcount_x = early_blockcount_x;
3628                                        current_blockcount_y = early_blockcount_y;
3629                                        current_blockcount_z = early_blockcount_z;
3630                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3631                                        if(*indicator_pos){
3632                                                // decompress by SZ
3633                                                double * block_data_pos = data_pos;
3634                                                double pred;
3635                                                size_t index = 0;
3636                                                int type_;
3637                                                size_t unpredictable_count = 0;
3638                                                // ii == 0
3639                                                {
3640                                                        // jj == 0
3641                                                        {
3642                                                                {
3643                                                                        // kk == 0
3644                                                                        type_ = type[index];
3645                                                                        if(type_ == 0){
3646                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3647                                                                        }
3648                                                                        else{
3649                                                                                pred = 0;
3650                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3651                                                                        }
3652                                                                        index ++;
3653                                                                        block_data_pos ++;
3654                                                                }
3655                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3656                                                                        type_ = type[index];
3657                                                                        if(type_ == 0){
3658                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3659                                                                        }
3660                                                                        else{
3661                                                                                pred = block_data_pos[- 1];
3662                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3663                                                                        }
3664                                                                        index ++;
3665                                                                        block_data_pos ++;
3666                                                                }
3667                                                                block_data_pos += dim1_offset - current_blockcount_z;
3668                                                        }
3669                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
3670                                                                {
3671                                                                        // kk == 0
3672                                                                        type_ = type[index];
3673                                                                        if(type_ == 0){
3674                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3675                                                                        }
3676                                                                        else{
3677                                                                                pred = block_data_pos[- dim1_offset];
3678                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3679                                                                        }
3680                                                                        index ++;
3681                                                                        block_data_pos ++;
3682                                                                }
3683                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3684                                                                        type_ = type[index];
3685                                                                        if(type_ == 0){
3686                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3687                                                                        }
3688                                                                        else{
3689                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3690                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3691                                                                        }
3692                                                                        index ++;
3693                                                                        block_data_pos ++;
3694                                                                }
3695                                                                block_data_pos += dim1_offset - current_blockcount_z;
3696                                                        }
3697                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;                                             
3698                                                }
3699                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
3700                                                        // jj == 0
3701                                                        {
3702                                                                {
3703                                                                        // kk == 0
3704                                                                        type_ = type[index];
3705                                                                        if(type_ == 0){
3706                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3707                                                                        }
3708                                                                        else{
3709                                                                                pred = block_data_pos[- dim0_offset];
3710                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3711                                                                        }
3712                                                                        index ++;
3713                                                                        block_data_pos ++;
3714                                                                }
3715                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3716                                                                        type_ = type[index];
3717                                                                        if(type_ == 0){
3718                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3719                                                                        }
3720                                                                        else{
3721                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3722                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3723                                                                        }
3724                                                                        index ++;
3725                                                                        block_data_pos ++;
3726                                                                }
3727                                                                block_data_pos += dim1_offset - current_blockcount_z;
3728                                                        }
3729                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
3730                                                                {
3731                                                                        // kk == 0
3732                                                                        type_ = type[index];
3733                                                                        if(type_ == 0){
3734                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3735                                                                        }
3736                                                                        else{
3737                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
3738                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3739                                                                        }
3740                                                                        index ++;
3741                                                                        block_data_pos ++;
3742                                                                }
3743                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3744                                                                        type_ = type[index];
3745                                                                        if(type_ == 0){
3746                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3747                                                                        }
3748                                                                        else{
3749                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3750                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3751                                                                        }
3752                                                                        index ++;
3753                                                                        block_data_pos ++;
3754                                                                }
3755                                                                block_data_pos += dim1_offset - current_blockcount_z;
3756                                                        }
3757                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3758                                                }
3759                                                cur_unpred_count = unpredictable_count;
3760                                        }
3761                                        else{
3762                                                // decompress by regression
3763                                                {
3764                                                        //restore regression coefficients
3765                                                        double pred;
3766                                                        int type_;
3767                                                        for(int e=0; e<4; e++){
3768                                                                type_ = coeff_type[e][coeff_index];
3769                                                                if (type_ != 0){
3770                                                                        pred = last_coefficients[e];
3771                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3772                                                                }
3773                                                                else{
3774                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3775                                                                        coeff_unpred_data_count[e] ++;
3776                                                                }
3777                                                        }
3778                                                        coeff_index ++;
3779                                                }
3780                                                {
3781                                                        double * block_data_pos = data_pos;
3782                                                        double pred;
3783                                                        int type_;
3784                                                        size_t index = 0;
3785                                                        size_t unpredictable_count = 0;
3786                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3787                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3788                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3789                                                                                type_ = type[index];
3790                                                                                if (type_ != 0){
3791                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3792                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3793                                                                                }
3794                                                                                else{
3795                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3796                                                                                }
3797                                                                                index ++;       
3798                                                                                block_data_pos ++;
3799                                                                        }
3800                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3801                                                                }
3802                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3803                                                        }
3804                                                        cur_unpred_count = unpredictable_count;
3805                                                }
3806                                        }
3807                                        indicator_pos ++;
3808                                        type += current_block_elements;
3809                                        unpred_data += cur_unpred_count;
3810                                } // end k == 0
3811                                // i == 0 j == 0 k != 0
3812                                for(size_t k=1; k<num_z; k++){
3813                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
3814                                        data_pos = *data + offset_z;
3815
3816                                        current_blockcount_x = early_blockcount_x;
3817                                        current_blockcount_y = early_blockcount_y;
3818                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
3819
3820                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3821                                        if(*indicator_pos){
3822                                                // decompress by SZ
3823                                                double * block_data_pos = data_pos;
3824                                                double pred;
3825                                                size_t index = 0;
3826                                                int type_;
3827                                                size_t unpredictable_count = 0;
3828                                                // ii == 0
3829                                                {
3830                                                        // jj == 0
3831                                                        {
3832                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3833                                                                        type_ = type[index];
3834                                                                        if(type_ == 0){
3835                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3836                                                                        }
3837                                                                        else{
3838                                                                                pred = block_data_pos[- 1];
3839                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3840                                                                        }
3841                                                                        index ++;
3842                                                                        block_data_pos ++;
3843                                                                }
3844                                                                block_data_pos += dim1_offset - current_blockcount_z;
3845                                                        }
3846                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
3847                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3848                                                                        type_ = type[index];
3849                                                                        if(type_ == 0){
3850                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3851                                                                        }
3852                                                                        else{
3853                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3854                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3855                                                                        }
3856                                                                        index ++;
3857                                                                        block_data_pos ++;
3858                                                                }
3859                                                                block_data_pos += dim1_offset - current_blockcount_z;
3860                                                        }
3861                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3862                                                }
3863                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
3864                                                        // jj == 0
3865                                                        {
3866                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3867                                                                        type_ = type[index];
3868                                                                        if(type_ == 0){
3869                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3870                                                                        }
3871                                                                        else{
3872                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
3873                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3874                                                                        }
3875                                                                        index ++;
3876                                                                        block_data_pos ++;
3877                                                                }
3878                                                                block_data_pos += dim1_offset - current_blockcount_z;
3879                                                        }
3880                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
3881                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
3882                                                                        type_ = type[index];
3883                                                                        if(type_ == 0){
3884                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3885                                                                        }
3886                                                                        else{
3887                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
3888                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3889                                                                        }
3890                                                                        index ++;
3891                                                                        block_data_pos ++;
3892                                                                }
3893                                                                block_data_pos += dim1_offset - current_blockcount_z;
3894                                                        }
3895                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3896                                                }
3897                                                cur_unpred_count = unpredictable_count;
3898                                        }
3899                                        else{
3900                                                // decompress by regression
3901                                                {
3902                                                        //restore regression coefficients
3903                                                        double pred;
3904                                                        int type_;
3905                                                        for(int e=0; e<4; e++){
3906                                                                type_ = coeff_type[e][coeff_index];
3907                                                                if (type_ != 0){
3908                                                                        pred = last_coefficients[e];
3909                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
3910                                                                }
3911                                                                else{
3912                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
3913                                                                        coeff_unpred_data_count[e] ++;
3914                                                                }
3915                                                        }
3916                                                        coeff_index ++;
3917                                                }
3918                                                {
3919                                                        double * block_data_pos = data_pos;
3920                                                        double pred;
3921                                                        int type_;
3922                                                        size_t index = 0;
3923                                                        size_t unpredictable_count = 0;
3924                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
3925                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
3926                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
3927                                                                                type_ = type[index];
3928                                                                                if (type_ != 0){
3929                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
3930                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3931                                                                                }
3932                                                                                else{
3933                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
3934                                                                                }
3935                                                                                index ++;       
3936                                                                                block_data_pos ++;
3937                                                                        }
3938                                                                        block_data_pos += dim1_offset - current_blockcount_z;
3939                                                                }
3940                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3941                                                        }
3942                                                        cur_unpred_count = unpredictable_count;
3943                                                }
3944                                        }
3945                                        indicator_pos ++;
3946                                        type += current_block_elements;
3947                                        unpred_data += cur_unpred_count;
3948                                }
3949                        }// end j==0
3950                        for(size_t j=1; j<num_y; j++){
3951                                // k == 0
3952                                {
3953                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
3954                                        data_pos = *data + offset_y * dim1_offset;
3955
3956                                        current_blockcount_x = early_blockcount_x;
3957                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
3958                                        current_blockcount_z = early_blockcount_z;
3959                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
3960                                        if(*indicator_pos){
3961                                                // decompress by SZ
3962                                                double * block_data_pos = data_pos;
3963                                                double pred;
3964                                                size_t index = 0;
3965                                                int type_;
3966                                                size_t unpredictable_count = 0;
3967                                                // ii == 0
3968                                                {
3969                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
3970                                                                {
3971                                                                        // kk == 0
3972                                                                        type_ = type[index];
3973                                                                        if(type_ == 0){
3974                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3975                                                                        }
3976                                                                        else{
3977                                                                                pred = block_data_pos[- dim1_offset];
3978                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3979                                                                        }
3980                                                                        index ++;
3981                                                                        block_data_pos ++;
3982                                                                }
3983                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
3984                                                                        type_ = type[index];
3985                                                                        if(type_ == 0){
3986                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
3987                                                                        }
3988                                                                        else{
3989                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
3990                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
3991                                                                        }
3992                                                                        index ++;
3993                                                                        block_data_pos ++;
3994                                                                }
3995                                                                block_data_pos += dim1_offset - current_blockcount_z;
3996                                                        }
3997                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
3998                                                }
3999                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
4000                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
4001                                                                {
4002                                                                        // kk == 0
4003                                                                        type_ = type[index];
4004                                                                        if(type_ == 0){
4005                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4006                                                                        }
4007                                                                        else{
4008                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
4009                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4010                                                                        }
4011                                                                        index ++;
4012                                                                        block_data_pos ++;
4013                                                                }
4014                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
4015                                                                        type_ = type[index];
4016                                                                        if(type_ == 0){
4017                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4018                                                                        }
4019                                                                        else{
4020                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4021                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4022                                                                        }
4023                                                                        index ++;
4024                                                                        block_data_pos ++;
4025                                                                }
4026                                                                block_data_pos += dim1_offset - current_blockcount_z;
4027                                                        }
4028                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4029                                                }
4030                                                cur_unpred_count = unpredictable_count;
4031                                        }
4032                                        else{
4033                                                // decompress by regression
4034                                                {
4035                                                        //restore regression coefficients
4036                                                        double pred;
4037                                                        int type_;
4038                                                        for(int e=0; e<4; e++){
4039                                                                type_ = coeff_type[e][coeff_index];
4040                                                                if (type_ != 0){
4041                                                                        pred = last_coefficients[e];
4042                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4043                                                                }
4044                                                                else{
4045                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4046                                                                        coeff_unpred_data_count[e] ++;
4047                                                                }
4048                                                        }
4049                                                        coeff_index ++;
4050                                                }
4051                                                {
4052                                                        double * block_data_pos = data_pos;
4053                                                        double pred;
4054                                                        int type_;
4055                                                        size_t index = 0;
4056                                                        size_t unpredictable_count = 0;
4057                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
4058                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
4059                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
4060                                                                                type_ = type[index];
4061                                                                                if (type_ != 0){
4062                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4063                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4064                                                                                }
4065                                                                                else{
4066                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
4067                                                                                }
4068                                                                                index ++;       
4069                                                                                block_data_pos ++;
4070                                                                        }
4071                                                                        block_data_pos += dim1_offset - current_blockcount_z;
4072                                                                }
4073                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4074                                                        }
4075                                                        cur_unpred_count = unpredictable_count;
4076                                                }
4077                                        }
4078                                        indicator_pos ++;
4079                                        type += current_block_elements;
4080                                        unpred_data += cur_unpred_count;
4081                                } // end k == 0
4082                                for(size_t k=1; k<num_z; k++){
4083                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
4084                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
4085                                        data_pos = *data + offset_y * dim1_offset + offset_z;
4086
4087                                        current_blockcount_x = early_blockcount_x;
4088                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
4089                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
4090
4091                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4092                                        if(*indicator_pos){
4093                                                // decompress by SZ
4094                                                double * block_data_pos = data_pos;
4095                                                double pred;
4096                                                size_t index = 0;
4097                                                int type_;
4098                                                size_t unpredictable_count = 0;
4099                                                // ii == 0
4100                                                {
4101                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
4102                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
4103                                                                        type_ = type[index];
4104                                                                        if(type_ == 0){
4105                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4106                                                                        }
4107                                                                        else{
4108                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] - block_data_pos[- dim1_offset - 1];
4109                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4110                                                                        }
4111                                                                        index ++;
4112                                                                        block_data_pos ++;
4113                                                                }
4114                                                                block_data_pos += dim1_offset - current_blockcount_z;
4115                                                        }
4116                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4117                                                }
4118                                                for(size_t ii=1; ii<current_blockcount_x; ii++){
4119                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
4120                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
4121                                                                        type_ = type[index];
4122                                                                        if(type_ == 0){
4123                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4124                                                                        }
4125                                                                        else{
4126                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4127                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4128                                                                        }
4129                                                                        index ++;
4130                                                                        block_data_pos ++;
4131                                                                }
4132                                                                block_data_pos += dim1_offset - current_blockcount_z;
4133                                                        }
4134                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4135                                                }
4136                                                cur_unpred_count = unpredictable_count;
4137                                        }
4138                                        else{
4139                                                // decompress by regression
4140                                                {
4141                                                        //restore regression coefficients
4142                                                        double pred;
4143                                                        int type_;
4144                                                        for(int e=0; e<4; e++){
4145                                                                type_ = coeff_type[e][coeff_index];
4146                                                                if (type_ != 0){
4147                                                                        pred = last_coefficients[e];
4148                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4149                                                                }
4150                                                                else{
4151                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4152                                                                        coeff_unpred_data_count[e] ++;
4153                                                                }
4154                                                        }
4155                                                        coeff_index ++;
4156                                                }
4157                                                {
4158                                                        double * block_data_pos = data_pos;
4159                                                        double pred;
4160                                                        int type_;
4161                                                        size_t index = 0;
4162                                                        size_t unpredictable_count = 0;
4163                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
4164                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
4165                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
4166                                                                                type_ = type[index];
4167                                                                                if (type_ != 0){
4168                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4169                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4170                                                                                }
4171                                                                                else{
4172                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
4173                                                                                }
4174                                                                                index ++;       
4175                                                                                block_data_pos ++;
4176                                                                        }
4177                                                                        block_data_pos += dim1_offset - current_blockcount_z;
4178                                                                }
4179                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4180                                                        }
4181                                                        cur_unpred_count = unpredictable_count;
4182                                                }
4183                                        }
4184                                        indicator_pos ++;
4185                                        type += current_block_elements;
4186                                        unpred_data += cur_unpred_count;
4187                                }
4188                        }
4189                } // end i==0
4190                for(size_t i=1; i<num_x; i++){
4191                        // j == 0
4192                        {
4193                                // k == 0
4194                                {
4195                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4196                                        data_pos = *data + offset_x * dim0_offset;
4197
4198                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4199                                        current_blockcount_y = early_blockcount_y;
4200                                        current_blockcount_z = early_blockcount_z;
4201                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4202                                        if(*indicator_pos){
4203                                                // decompress by SZ
4204                                                double * block_data_pos = data_pos;
4205                                                double pred;
4206                                                size_t index = 0;
4207                                                int type_;
4208                                                size_t unpredictable_count = 0;
4209                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
4210                                                        // jj == 0
4211                                                        {
4212                                                                {
4213                                                                        // kk == 0
4214                                                                        type_ = type[index];
4215                                                                        if(type_ == 0){
4216                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4217                                                                        }
4218                                                                        else{
4219                                                                                pred = block_data_pos[- dim0_offset];
4220                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4221                                                                        }
4222                                                                        index ++;
4223                                                                        block_data_pos ++;
4224                                                                }
4225                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
4226                                                                        type_ = type[index];
4227                                                                        if(type_ == 0){
4228                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4229                                                                        }
4230                                                                        else{
4231                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
4232                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4233                                                                        }
4234                                                                        index ++;
4235                                                                        block_data_pos ++;
4236                                                                }
4237                                                                block_data_pos += dim1_offset - current_blockcount_z;
4238                                                        }
4239                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
4240                                                                {
4241                                                                        // kk == 0
4242                                                                        type_ = type[index];
4243                                                                        if(type_ == 0){
4244                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4245                                                                        }
4246                                                                        else{
4247                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
4248                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4249                                                                        }
4250                                                                        index ++;
4251                                                                        block_data_pos ++;
4252                                                                }
4253                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
4254                                                                        type_ = type[index];
4255                                                                        if(type_ == 0){
4256                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4257                                                                        }
4258                                                                        else{
4259                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4260                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4261                                                                        }
4262                                                                        index ++;
4263                                                                        block_data_pos ++;
4264                                                                }
4265                                                                block_data_pos += dim1_offset - current_blockcount_z;
4266                                                        }
4267                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4268                                                }
4269                                                cur_unpred_count = unpredictable_count;
4270                                        }
4271                                        else{
4272                                                // decompress by regression
4273                                                {
4274                                                        //restore regression coefficients
4275                                                        double pred;
4276                                                        int type_;
4277                                                        for(int e=0; e<4; e++){
4278                                                                type_ = coeff_type[e][coeff_index];
4279                                                                if (type_ != 0){
4280                                                                        pred = last_coefficients[e];
4281                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4282                                                                }
4283                                                                else{
4284                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4285                                                                        coeff_unpred_data_count[e] ++;
4286                                                                }
4287                                                        }
4288                                                        coeff_index ++;
4289                                                }
4290                                                {
4291                                                        double * block_data_pos = data_pos;
4292                                                        double pred;
4293                                                        int type_;
4294                                                        size_t index = 0;
4295                                                        size_t unpredictable_count = 0;
4296                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
4297                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
4298                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
4299                                                                                type_ = type[index];
4300                                                                                if (type_ != 0){
4301                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4302                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4303                                                                                }
4304                                                                                else{
4305                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
4306                                                                                }
4307                                                                                index ++;       
4308                                                                                block_data_pos ++;
4309                                                                        }
4310                                                                        block_data_pos += dim1_offset - current_blockcount_z;
4311                                                                }
4312                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4313                                                        }
4314                                                        cur_unpred_count = unpredictable_count;
4315                                                }
4316                                        }
4317                                        indicator_pos ++;
4318                                        type += current_block_elements;
4319                                        unpred_data += cur_unpred_count;
4320                                } // end k == 0
4321                                for(size_t k=1; k<num_z; k++){
4322                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4323                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
4324                                        data_pos = *data + offset_x * dim0_offset + offset_z;
4325
4326                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4327                                        current_blockcount_y = early_blockcount_y;
4328                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
4329
4330                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4331                                        if(*indicator_pos){
4332                                                // decompress by SZ
4333                                                double * block_data_pos = data_pos;
4334                                                double pred;
4335                                                size_t index = 0;
4336                                                int type_;
4337                                                size_t unpredictable_count = 0;
4338                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
4339                                                        // jj == 0
4340                                                        {
4341                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
4342                                                                        type_ = type[index];
4343                                                                        if(type_ == 0){
4344                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4345                                                                        }
4346                                                                        else{
4347                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - 1];
4348                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4349                                                                        }
4350                                                                        index ++;
4351                                                                        block_data_pos ++;
4352                                                                }
4353                                                                block_data_pos += dim1_offset - current_blockcount_z;
4354                                                        }
4355                                                        for(size_t jj=1; jj<current_blockcount_y; jj++){
4356                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
4357                                                                        type_ = type[index];
4358                                                                        if(type_ == 0){
4359                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4360                                                                        }
4361                                                                        else{
4362                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4363                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4364                                                                        }
4365                                                                        index ++;
4366                                                                        block_data_pos ++;
4367                                                                }
4368                                                                block_data_pos += dim1_offset - current_blockcount_z;
4369                                                        }
4370                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4371                                                }
4372                                                cur_unpred_count = unpredictable_count;
4373                                        }
4374                                        else{
4375                                                // decompress by regression
4376                                                {
4377                                                        //restore regression coefficients
4378                                                        double pred;
4379                                                        int type_;
4380                                                        for(int e=0; e<4; e++){
4381                                                                type_ = coeff_type[e][coeff_index];
4382                                                                if (type_ != 0){
4383                                                                        pred = last_coefficients[e];
4384                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4385                                                                }
4386                                                                else{
4387                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4388                                                                        coeff_unpred_data_count[e] ++;
4389                                                                }
4390                                                        }
4391                                                        coeff_index ++;
4392                                                }
4393                                                {
4394                                                        double * block_data_pos = data_pos;
4395                                                        double pred;
4396                                                        int type_;
4397                                                        size_t index = 0;
4398                                                        size_t unpredictable_count = 0;
4399                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
4400                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
4401                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
4402                                                                                type_ = type[index];
4403                                                                                if (type_ != 0){
4404                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4405                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4406                                                                                }
4407                                                                                else{
4408                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
4409                                                                                }
4410                                                                                index ++;       
4411                                                                                block_data_pos ++;
4412                                                                        }
4413                                                                        block_data_pos += dim1_offset - current_blockcount_z;
4414                                                                }
4415                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4416                                                        }
4417                                                        cur_unpred_count = unpredictable_count;
4418                                                }
4419                                        }
4420                                        indicator_pos ++;
4421                                        type += current_block_elements;
4422                                        unpred_data += cur_unpred_count;
4423                                }
4424                        }// end j = 0
4425                        for(size_t j=1; j<num_y; j++){
4426                                // k == 0
4427                                {
4428                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4429                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
4430                                        data_pos = *data + offset_x * dim0_offset + offset_y * dim1_offset;
4431
4432                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4433                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
4434                                        current_blockcount_z = early_blockcount_z;
4435                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4436                                        if(*indicator_pos){
4437                                                // decompress by SZ
4438                                                double * block_data_pos = data_pos;
4439                                                double pred;
4440                                                size_t index = 0;
4441                                                int type_;
4442                                                size_t unpredictable_count = 0;
4443                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
4444                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
4445                                                                {
4446                                                                        // kk == 0
4447                                                                        type_ = type[index];
4448                                                                        if(type_ == 0){
4449                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4450                                                                        }
4451                                                                        else{
4452                                                                                pred = block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim0_offset - dim1_offset];
4453                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4454                                                                        }
4455                                                                        index ++;
4456                                                                        block_data_pos ++;
4457                                                                }
4458                                                                for(size_t kk=1; kk<current_blockcount_z; kk++){
4459                                                                        type_ = type[index];
4460                                                                        if(type_ == 0){
4461                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4462                                                                        }
4463                                                                        else{
4464                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4465                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4466                                                                        }
4467                                                                        index ++;
4468                                                                        block_data_pos ++;
4469                                                                }
4470                                                                block_data_pos += dim1_offset - current_blockcount_z;
4471                                                        }
4472                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4473                                                }
4474                                                cur_unpred_count = unpredictable_count;
4475                                        }
4476                                        else{
4477                                                // decompress by regression
4478                                                {
4479                                                        //restore regression coefficients
4480                                                        double pred;
4481                                                        int type_;
4482                                                        for(int e=0; e<4; e++){
4483                                                                type_ = coeff_type[e][coeff_index];
4484                                                                if (type_ != 0){
4485                                                                        pred = last_coefficients[e];
4486                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4487                                                                }
4488                                                                else{
4489                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4490                                                                        coeff_unpred_data_count[e] ++;
4491                                                                }
4492                                                        }
4493                                                        coeff_index ++;
4494                                                }
4495                                                {
4496                                                        double * block_data_pos = data_pos;
4497                                                        double pred;
4498                                                        int type_;
4499                                                        size_t index = 0;
4500                                                        size_t unpredictable_count = 0;
4501                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
4502                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
4503                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
4504                                                                                type_ = type[index];
4505                                                                                if (type_ != 0){
4506                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4507                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4508                                                                                }
4509                                                                                else{
4510                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
4511                                                                                }
4512                                                                                index ++;       
4513                                                                                block_data_pos ++;
4514                                                                        }
4515                                                                        block_data_pos += dim1_offset - current_blockcount_z;
4516                                                                }
4517                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4518                                                        }
4519                                                        cur_unpred_count = unpredictable_count;
4520                                                }
4521                                        }
4522                                        indicator_pos ++;
4523                                        type += current_block_elements;
4524                                        unpred_data += cur_unpred_count;
4525                                } // end k == 0
4526                                for(size_t k=1; k<num_z; k++){
4527                                        offset_x = (i < split_index_x) ? i * early_blockcount_x : i * late_blockcount_x + split_index_x;
4528                                        offset_y = (j < split_index_y) ? j * early_blockcount_y : j * late_blockcount_y + split_index_y;
4529                                        offset_z = (k < split_index_z) ? k * early_blockcount_z : k * late_blockcount_z + split_index_z;
4530                                        data_pos = *data + offset_x * dim0_offset + offset_y * dim1_offset + offset_z;
4531
4532                                        current_blockcount_x = (i < split_index_x) ? early_blockcount_x : late_blockcount_x;
4533                                        current_blockcount_y = (j < split_index_y) ? early_blockcount_y : late_blockcount_y;
4534                                        current_blockcount_z = (k < split_index_z) ? early_blockcount_z : late_blockcount_z;
4535
4536                                        size_t current_block_elements = current_blockcount_x * current_blockcount_y * current_blockcount_z;
4537                                        if(*indicator_pos){
4538                                                // decompress by SZ
4539                                                double * block_data_pos = data_pos;
4540                                                double pred;
4541                                                size_t index = 0;
4542                                                int type_;
4543                                                size_t unpredictable_count = 0;
4544                                                for(size_t ii=0; ii<current_blockcount_x; ii++){
4545                                                        for(size_t jj=0; jj<current_blockcount_y; jj++){
4546                                                                for(size_t kk=0; kk<current_blockcount_z; kk++){
4547                                                                        type_ = type[index];
4548                                                                        if(type_ == 0){
4549                                                                                *block_data_pos = unpred_data[unpredictable_count ++];
4550                                                                        }
4551                                                                        else{
4552                                                                                pred = block_data_pos[- 1] + block_data_pos[- dim1_offset] + block_data_pos[- dim0_offset] - block_data_pos[- dim1_offset - 1] - block_data_pos[- dim0_offset - 1] - block_data_pos[- dim0_offset - dim1_offset] + block_data_pos[- dim0_offset - dim1_offset - 1];
4553                                                                                *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4554                                                                        }
4555                                                                        index ++;
4556                                                                        block_data_pos ++;
4557                                                                }
4558                                                                block_data_pos += dim1_offset - current_blockcount_z;
4559                                                        }
4560                                                        block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4561                                                }
4562                                                cur_unpred_count = unpredictable_count;
4563                                        }
4564                                        else{
4565                                                // decompress by regression
4566                                                {
4567                                                        //restore regression coefficients
4568                                                        double pred;
4569                                                        int type_;
4570                                                        for(int e=0; e<4; e++){
4571                                                                type_ = coeff_type[e][coeff_index];
4572                                                                if (type_ != 0){
4573                                                                        pred = last_coefficients[e];
4574                                                                        last_coefficients[e] = pred + 2 * (type_ - coeff_intvRadius[e]) * precision[e];
4575                                                                }
4576                                                                else{
4577                                                                        last_coefficients[e] = coeff_unpred_data[e][coeff_unpred_data_count[e]];
4578                                                                        coeff_unpred_data_count[e] ++;
4579                                                                }
4580                                                        }
4581                                                        coeff_index ++;
4582                                                }
4583                                                {
4584                                                        double * block_data_pos = data_pos;
4585                                                        double pred;
4586                                                        int type_;
4587                                                        size_t index = 0;
4588                                                        size_t unpredictable_count = 0;
4589                                                        for(size_t ii=0; ii<current_blockcount_x; ii++){
4590                                                                for(size_t jj=0; jj<current_blockcount_y; jj++){
4591                                                                        for(size_t kk=0; kk<current_blockcount_z; kk++){
4592                                                                                type_ = type[index];
4593                                                                                if (type_ != 0){
4594                                                                                        pred = last_coefficients[0] * ii + last_coefficients[1] * jj + last_coefficients[2] * kk + last_coefficients[3];
4595                                                                                        *block_data_pos = pred + 2 * (type_ - intvRadius) * realPrecision;
4596                                                                                }
4597                                                                                else{
4598                                                                                        *block_data_pos = unpred_data[unpredictable_count ++];
4599                                                                                }
4600                                                                                index ++;       
4601                                                                                block_data_pos ++;
4602                                                                        }
4603                                                                        block_data_pos += dim1_offset - current_blockcount_z;
4604                                                                }
4605                                                                block_data_pos += dim0_offset - current_blockcount_y * dim1_offset;
4606                                                        }
4607                                                        cur_unpred_count = unpredictable_count;
4608                                                }
4609                                        }
4610                                        indicator_pos ++;
4611                                        type += current_block_elements;
4612                                        unpred_data += cur_unpred_count;
4613                                }
4614                        }
4615                }
4616        }
4617
4618        free(coeff_result_type);
4619
4620        free(indicator);
4621        free(result_type);
4622}
Note: See TracBrowser for help on using the repository browser.