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

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

importing new SZ files
  • Property mode set to 100644
Line 
1/**
2 *  @file szd_float_pwr.c
3 *  @author Sheng Di
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 <math.h>
14#include "TightDataPointStorageF.h"
15#include "CompressElement.h"
16#include "sz.h"
17#include "Huffman.h"
18#include "sz_float_pwr.h"
19#include "utility.h"
20//#include "rw.h"
21//
22#pragma GCC diagnostic push
23#pragma GCC diagnostic ignored "-Wchar-subscripts"
24
25
26void decompressDataSeries_float_1D_pwr(float** data, size_t dataSeriesLength, TightDataPointStorageF* tdps) 
27{
28        updateQuantizationInfo(tdps->intervals);
29        unsigned char tmpPrecBytes[4] = {0}; //used when needing to convert bytes to float values
30        unsigned char* bp = tdps->pwrErrBoundBytes;
31        size_t i, j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
32                                                                // in resiMidBits, p is to track the
33                                                                // byte_index of resiMidBits, l is for
34                                                                // leadNum
35        unsigned char* leadNum;
36        float interval = 0;// = (float)tdps->realPrecision*2;
37       
38        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
39
40        *data = (float*)malloc(sizeof(float)*dataSeriesLength);
41
42        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
43
44        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
45        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
46        SZ_ReleaseHuffman(huffmanTree); 
47
48        //sdi:Debug
49        //writeUShortData(type, dataSeriesLength, "decompressStateBytes.sb");
50
51        unsigned char preBytes[4];
52        unsigned char curBytes[4];
53       
54        memset(preBytes, 0, 4);
55
56        size_t curByteIndex = 0;
57        int reqLength = 0, reqBytesLength = 0, resiBitsLength = 0, resiBits = 0; 
58        unsigned char leadingNum;       
59        float medianValue, exactData, predValue = 0, realPrecision = 0;
60       
61        medianValue = tdps->medianValue;
62       
63        int type_, updateReqLength = 0;
64        for (i = 0; i < dataSeriesLength; i++) 
65        {
66                if(i%tdps->segment_size==0)
67                {
68                        tmpPrecBytes[0] = *(bp++);
69                        tmpPrecBytes[1] = *(bp++);
70                        tmpPrecBytes[2] = 0;
71                        tmpPrecBytes[3] = 0;
72                        realPrecision = bytesToFloat(tmpPrecBytes);
73                        interval = realPrecision*2;
74                        updateReqLength = 0;
75                }
76               
77                type_ = type[i];
78                switch (type_) {
79                case 0:
80                        // compute resiBits
81                        if(updateReqLength==0)
82                        {
83                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
84                                reqBytesLength = reqLength/8;
85                                resiBitsLength = reqLength%8;   
86                                updateReqLength = 1;   
87                        }
88                        resiBits = 0;
89                        if (resiBitsLength != 0) {
90                                int kMod8 = k % 8;
91                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
92                                if (rightMovSteps > 0) {
93                                        int code = getRightMovingCode(kMod8, resiBitsLength);
94                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
95                                } else if (rightMovSteps < 0) {
96                                        int code1 = getLeftMovingCode(kMod8);
97                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
98                                        int leftMovSteps = -rightMovSteps;
99                                        rightMovSteps = 8 - leftMovSteps;
100                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
101                                        p++;
102                                        resiBits = resiBits
103                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
104                                } else // rightMovSteps == 0
105                                {
106                                        int code = getRightMovingCode(kMod8, resiBitsLength);
107                                        resiBits = (tdps->residualMidBits[p] & code);
108                                        p++;
109                                }
110                                k += resiBitsLength;
111                        }
112
113                        // recover the exact data       
114                        memset(curBytes, 0, 4);
115                        leadingNum = leadNum[l++];
116                        memcpy(curBytes, preBytes, leadingNum);
117                        for (j = leadingNum; j < reqBytesLength; j++)
118                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
119                        if (resiBitsLength != 0) {
120                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
121                                curBytes[reqBytesLength] = resiByte;
122                        }
123                       
124                        exactData = bytesToFloat(curBytes);
125                        (*data)[i] = exactData + medianValue;
126                        memcpy(preBytes,curBytes,4);
127                        break;
128                default:
129                        //predValue = 2 * (*data)[i-1] - (*data)[i-2];
130                        predValue = (*data)[i-1];
131                        (*data)[i] = predValue + (type_-exe_params->intvRadius)*interval;
132                        break;
133                }
134                //printf("%.30G\n",(*data)[i]);
135        }
136        free(leadNum);
137        free(type);
138        return;
139}
140
141float* extractRealPrecision_2D_float(size_t R1, size_t R2, int blockSize, TightDataPointStorageF* tdps)
142{
143        size_t i,j,k=0, I;
144        unsigned char* bytes = tdps->pwrErrBoundBytes;
145        unsigned char tmpBytes[4] = {0};
146        float* result = (float*)malloc(sizeof(float)*R1*R2);
147        for(i=0;i<R1;i++)
148        {
149                I = i*R2;
150                for(j=0;j<R2;j++)
151                {
152                        tmpBytes[0] = bytes[k++];
153                        tmpBytes[1] = bytes[k++];
154                        result[I+j]=bytesToFloat(tmpBytes);
155                }
156        }
157        return result;
158}
159
160void decompressDataSeries_float_2D_pwr(float** data, size_t r1, size_t r2, TightDataPointStorageF* tdps) 
161{
162        updateQuantizationInfo(tdps->intervals);
163        //printf("tdps->intervals=%d, exe_params->intvRadius=%d\n", tdps->intervals, exe_params->intvRadius);
164       
165        size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
166        // in resiMidBits, p is to track the
167        // byte_index of resiMidBits, l is for
168        // leadNum
169        size_t dataSeriesLength = r1*r2;
170        //      printf ("%d %d\n", r1, r2);
171
172        unsigned char* leadNum;
173
174        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
175
176        *data = (float*)malloc(sizeof(float)*dataSeriesLength);
177
178        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
179
180        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
181        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
182        SZ_ReleaseHuffman(huffmanTree); 
183
184        unsigned char preBytes[4];
185        unsigned char curBytes[4];
186
187        memset(preBytes, 0, 4);
188
189        size_t curByteIndex = 0;
190        int reqLength, reqBytesLength, resiBitsLength, resiBits; 
191        unsigned char leadingNum;       
192        float medianValue, exactData, realPrecision;
193        int type_;     
194        float pred1D, pred2D;
195        size_t ii, jj, II = 0, JJ = 0, updateReqLength = 1;
196
197        int blockSize = computeBlockEdgeSize_2D(tdps->segment_size);
198        size_t R1 = 1+(r1-1)/blockSize;
199        size_t R2 = 1+(r2-1)/blockSize;         
200        float* pwrErrBound = extractRealPrecision_2D_float(R1, R2, blockSize, tdps);
201
202        realPrecision = pwrErrBound[0]; 
203        computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
204        reqBytesLength = reqLength/8;
205        resiBitsLength = reqLength%8;
206
207        /* Process Row-0, data 0 */
208
209        // compute resiBits
210        resiBits = 0;
211        if (resiBitsLength != 0) {
212                int kMod8 = k % 8;
213                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
214                if (rightMovSteps > 0) {
215                        int code = getRightMovingCode(kMod8, resiBitsLength);
216                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
217                } else if (rightMovSteps < 0) {
218                        int code1 = getLeftMovingCode(kMod8);
219                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
220                        int leftMovSteps = -rightMovSteps;
221                        rightMovSteps = 8 - leftMovSteps;
222                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
223                        p++;
224                        resiBits = resiBits
225                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
226                } else // rightMovSteps == 0
227                {
228                        int code = getRightMovingCode(kMod8, resiBitsLength);
229                        resiBits = (tdps->residualMidBits[p] & code);
230                        p++;
231                }
232                k += resiBitsLength;
233        }
234
235        // recover the exact data
236        memset(curBytes, 0, 4);
237        leadingNum = leadNum[l++];
238        memcpy(curBytes, preBytes, leadingNum);
239        for (j = leadingNum; j < reqBytesLength; j++)
240                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
241        if (resiBitsLength != 0) {
242                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
243                curBytes[reqBytesLength] = resiByte;
244        }
245
246        exactData = bytesToFloat(curBytes);
247        (*data)[0] = exactData + medianValue;
248        memcpy(preBytes,curBytes,4);
249
250        /* Process Row-0, data 1 */
251        type_ = type[1];
252        if (type_ != 0)
253        {
254                pred1D = (*data)[0];
255                (*data)[1] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
256        }
257        else
258        {
259                // compute resiBits             
260                resiBits = 0;
261                if (resiBitsLength != 0) {
262                        int kMod8 = k % 8;
263                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
264                        if (rightMovSteps > 0) {
265                                int code = getRightMovingCode(kMod8, resiBitsLength);
266                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
267                        } else if (rightMovSteps < 0) {
268                                int code1 = getLeftMovingCode(kMod8);
269                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
270                                int leftMovSteps = -rightMovSteps;
271                                rightMovSteps = 8 - leftMovSteps;
272                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
273                                p++;
274                                resiBits = resiBits
275                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
276                        } else // rightMovSteps == 0
277                        {
278                                int code = getRightMovingCode(kMod8, resiBitsLength);
279                                resiBits = (tdps->residualMidBits[p] & code);
280                                p++;
281                        }
282                        k += resiBitsLength;
283                }
284
285                // recover the exact data
286                memset(curBytes, 0, 4);
287                leadingNum = leadNum[l++];
288                memcpy(curBytes, preBytes, leadingNum);
289                for (j = leadingNum; j < reqBytesLength; j++)
290                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
291                if (resiBitsLength != 0) {
292                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
293                        curBytes[reqBytesLength] = resiByte;
294                }
295
296                exactData = bytesToFloat(curBytes);
297                (*data)[1] = exactData + medianValue;
298                memcpy(preBytes,curBytes,4);
299        }
300
301        /* Process Row-0, data 2 --> data r2-1 */
302        for (jj = 2; jj < r2; jj++)
303        {       
304                if(jj%blockSize==0)
305                {
306                        II = 0;
307                        JJ++;
308                        realPrecision = pwrErrBound[JJ];
309                        updateReqLength = 0;                   
310                }
311
312                type_ = type[jj];
313                if (type_ != 0)
314                {
315                        pred1D = 2*(*data)[jj-1] - (*data)[jj-2];                                               
316                        (*data)[jj] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
317                }
318                else
319                {
320                        // compute resiBits
321                        if(updateReqLength==0)
322                        {
323                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
324                                reqBytesLength = reqLength/8;
325                                resiBitsLength = reqLength%8;                           
326                                updateReqLength = 1;
327                        }
328                       
329                        resiBits = 0;
330                        if (resiBitsLength != 0) {
331                                int kMod8 = k % 8;
332                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
333                                if (rightMovSteps > 0) {
334                                        int code = getRightMovingCode(kMod8, resiBitsLength);
335                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
336                                } else if (rightMovSteps < 0) {
337                                        int code1 = getLeftMovingCode(kMod8);
338                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
339                                        int leftMovSteps = -rightMovSteps;
340                                        rightMovSteps = 8 - leftMovSteps;
341                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
342                                        p++;
343                                        resiBits = resiBits
344                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
345                                } else // rightMovSteps == 0
346                                {
347                                        int code = getRightMovingCode(kMod8, resiBitsLength);
348                                        resiBits = (tdps->residualMidBits[p] & code);
349                                        p++;
350                                }
351                                k += resiBitsLength;
352                        }
353
354                        // recover the exact data
355                        memset(curBytes, 0, 4);
356                        leadingNum = leadNum[l++];
357                        memcpy(curBytes, preBytes, leadingNum);
358                        for (j = leadingNum; j < reqBytesLength; j++)
359                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
360                        if (resiBitsLength != 0) {
361                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
362                                curBytes[reqBytesLength] = resiByte;
363                        }
364
365                        exactData = bytesToFloat(curBytes);
366                        (*data)[jj] = exactData + medianValue;
367                        memcpy(preBytes,curBytes,4);
368                }
369        }
370
371        size_t index;
372        /* Process Row-1 --> Row-r1-1 */
373        for (ii = 1; ii < r1; ii++)
374        {
375                /* Process row-ii data 0 */             
376                if(ii%blockSize==0)
377                        II++;
378                JJ = 0;
379                realPrecision = pwrErrBound[II*R2+JJ];                         
380                updateReqLength = 0;
381               
382                index = ii*r2;
383
384                type_ = type[index];
385                if (type_ != 0)
386                {
387                        pred1D = (*data)[index-r2];             
388                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
389                }
390                else
391                {
392                        // compute resiBits
393                        if(updateReqLength==0)
394                        {
395                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
396                                reqBytesLength = reqLength/8;
397                                resiBitsLength = reqLength%8;                           
398                                updateReqLength = 1;
399                        }
400                       
401                        resiBits = 0;
402                        if (resiBitsLength != 0) {
403                                int kMod8 = k % 8;
404                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
405                                if (rightMovSteps > 0) {
406                                        int code = getRightMovingCode(kMod8, resiBitsLength);
407                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
408                                } else if (rightMovSteps < 0) {
409                                        int code1 = getLeftMovingCode(kMod8);
410                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
411                                        int leftMovSteps = -rightMovSteps;
412                                        rightMovSteps = 8 - leftMovSteps;
413                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
414                                        p++;
415                                        resiBits = resiBits
416                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
417                                } else // rightMovSteps == 0
418                                {
419                                        int code = getRightMovingCode(kMod8, resiBitsLength);
420                                        resiBits = (tdps->residualMidBits[p] & code);
421                                        p++;
422                                }
423                                k += resiBitsLength;
424                        }
425
426                        // recover the exact data
427                        memset(curBytes, 0, 4);
428                        leadingNum = leadNum[l++];
429                        memcpy(curBytes, preBytes, leadingNum);
430                        for (j = leadingNum; j < reqBytesLength; j++)
431                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
432                        if (resiBitsLength != 0) {
433                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
434                                curBytes[reqBytesLength] = resiByte;
435                        }
436
437                        exactData = bytesToFloat(curBytes);
438                        (*data)[index] = exactData + medianValue;
439                        memcpy(preBytes,curBytes,4);
440                }
441
442                /* Process row-ii data 1 --> r2-1*/
443                for (jj = 1; jj < r2; jj++)
444                {
445                        index = ii*r2+jj;
446                       
447                        if(jj%blockSize==0)
448                                JJ++;
449                        realPrecision = pwrErrBound[II*R2+JJ];                 
450                        updateReqLength = 0;                   
451
452                        type_ = type[index];
453                        if (type_ != 0)
454                        {
455                                pred2D = (*data)[index-1] + (*data)[index-r2] - (*data)[index-r2-1];
456                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
457                        }
458                        else
459                        {
460                                // compute resiBits
461                                if(updateReqLength==0)
462                                {
463                                        computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
464                                        reqBytesLength = reqLength/8;
465                                        resiBitsLength = reqLength%8;                           
466                                        updateReqLength = 1;
467                                }                               
468                               
469                                resiBits = 0;
470                                if (resiBitsLength != 0) {
471                                        int kMod8 = k % 8;
472                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
473                                        if (rightMovSteps > 0) {
474                                                int code = getRightMovingCode(kMod8, resiBitsLength);
475                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
476                                        } else if (rightMovSteps < 0) {
477                                                int code1 = getLeftMovingCode(kMod8);
478                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
479                                                int leftMovSteps = -rightMovSteps;
480                                                rightMovSteps = 8 - leftMovSteps;
481                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
482                                                p++;
483                                                resiBits = resiBits
484                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
485                                        } else // rightMovSteps == 0
486                                        {
487                                                int code = getRightMovingCode(kMod8, resiBitsLength);
488                                                resiBits = (tdps->residualMidBits[p] & code);
489                                                p++;
490                                        }
491                                        k += resiBitsLength;
492                                }
493
494                                // recover the exact data
495                                memset(curBytes, 0, 4);
496                                leadingNum = leadNum[l++];
497                                memcpy(curBytes, preBytes, leadingNum);
498                                for (j = leadingNum; j < reqBytesLength; j++)
499                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
500                                if (resiBitsLength != 0) {
501                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
502                                        curBytes[reqBytesLength] = resiByte;
503                                }
504
505                                exactData = bytesToFloat(curBytes);
506                                (*data)[index] = exactData + medianValue;
507                                memcpy(preBytes,curBytes,4);
508                        }
509                }
510        }
511
512        free(pwrErrBound);
513        free(leadNum);
514        free(type);
515        return;
516}
517
518float* extractRealPrecision_3D_float(size_t R1, size_t R2, size_t R3, int blockSize, TightDataPointStorageF* tdps)
519{
520        size_t i,j,k=0, IR, JR, p = 0;
521        size_t R23 = R2*R3;
522        unsigned char* bytes = tdps->pwrErrBoundBytes;
523        unsigned char tmpBytes[4] = {0};
524        float* result = (float*)malloc(sizeof(float)*R1*R2*R3);
525        for(i=0;i<R1;i++)
526        {
527                IR = i*R23;
528                for(j=0;j<R2;j++)
529                {
530                        JR = j*R3;
531                        for(k=0;k<R3;k++)
532                        {
533                                tmpBytes[0] = bytes[p++];
534                                tmpBytes[1] = bytes[p++];
535                                result[IR+JR+k]=bytesToFloat(tmpBytes);                         
536                        }
537                }
538        }
539        return result;
540}
541
542void decompressDataSeries_float_3D_pwr(float** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageF* 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        unsigned char* leadNum;
553
554        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
555
556        *data = (float*)malloc(sizeof(float)*dataSeriesLength);
557        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
558
559        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
560        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
561        SZ_ReleaseHuffman(huffmanTree); 
562
563        unsigned char preBytes[4];
564        unsigned char curBytes[4];
565
566        memset(preBytes, 0, 4);
567        size_t curByteIndex = 0;
568        int reqLength, reqBytesLength, resiBitsLength, resiBits; 
569        unsigned char leadingNum;
570        float medianValue, exactData, realPrecision;
571        int type_;
572        float pred1D, pred2D, pred3D;
573        size_t ii, jj, kk, II = 0, JJ = 0, KK = 0, updateReqLength = 1;
574
575        int blockSize = computeBlockEdgeSize_3D(tdps->segment_size);
576        size_t R1 = 1+(r1-1)/blockSize;
577        size_t R2 = 1+(r2-1)/blockSize;         
578        size_t R3 = 1+(r3-1)/blockSize;
579        size_t R23 = R2*R3;
580        float* pwrErrBound = extractRealPrecision_3D_float(R1, R2, R3, blockSize, tdps);
581
582        realPrecision = pwrErrBound[0]; 
583        computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
584        reqBytesLength = reqLength/8;
585        resiBitsLength = reqLength%8;
586
587        ///////////////////////////     Process layer-0 ///////////////////////////
588        /* Process Row-0 data 0*/
589        // compute resiBits
590        resiBits = 0;
591        if (resiBitsLength != 0) {
592                int kMod8 = k % 8;
593                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
594                if (rightMovSteps > 0) {
595                        int code = getRightMovingCode(kMod8, resiBitsLength);
596                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
597                } else if (rightMovSteps < 0) {
598                        int code1 = getLeftMovingCode(kMod8);
599                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
600                        int leftMovSteps = -rightMovSteps;
601                        rightMovSteps = 8 - leftMovSteps;
602                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
603                        p++;
604                        resiBits = resiBits
605                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
606                } else // rightMovSteps == 0
607                {
608                        int code = getRightMovingCode(kMod8, resiBitsLength);
609                        resiBits = (tdps->residualMidBits[p] & code);
610                        p++;
611                }
612                k += resiBitsLength;
613        }
614
615        // recover the exact data
616        memset(curBytes, 0, 4);
617        leadingNum = leadNum[l++];
618        memcpy(curBytes, preBytes, leadingNum);
619        for (j = leadingNum; j < reqBytesLength; j++)
620                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
621        if (resiBitsLength != 0) {
622                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
623                curBytes[reqBytesLength] = resiByte;
624        }
625        exactData = bytesToFloat(curBytes);
626        (*data)[0] = exactData + medianValue;
627        memcpy(preBytes,curBytes,4);
628
629        /* Process Row-0, data 1 */
630        pred1D = (*data)[0];
631
632        type_ = type[1];
633        if (type_ != 0)
634        {
635                (*data)[1] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
636        }
637        else
638        {
639                // compute resiBits
640                resiBits = 0;
641                if (resiBitsLength != 0) {
642                        int kMod8 = k % 8;
643                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
644                        if (rightMovSteps > 0) {
645                                int code = getRightMovingCode(kMod8, resiBitsLength);
646                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
647                        } else if (rightMovSteps < 0) {
648                                int code1 = getLeftMovingCode(kMod8);
649                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
650                                int leftMovSteps = -rightMovSteps;
651                                rightMovSteps = 8 - leftMovSteps;
652                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
653                                p++;
654                                resiBits = resiBits
655                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
656                        } else // rightMovSteps == 0
657                        {
658                                int code = getRightMovingCode(kMod8, resiBitsLength);
659                                resiBits = (tdps->residualMidBits[p] & code);
660                                p++;
661                        }
662                        k += resiBitsLength;
663                }
664
665                // recover the exact data
666                memset(curBytes, 0, 4);
667                leadingNum = leadNum[l++];
668                memcpy(curBytes, preBytes, leadingNum);
669                for (j = leadingNum; j < reqBytesLength; j++)
670                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
671                if (resiBitsLength != 0) {
672                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
673                        curBytes[reqBytesLength] = resiByte;
674                }
675
676                exactData = bytesToFloat(curBytes);
677                (*data)[1] = exactData + medianValue;
678                memcpy(preBytes,curBytes,4);
679        }
680        /* Process Row-0, data 2 --> data r3-1 */
681        for (jj = 2; jj < r3; jj++)
682        {
683                if(jj%blockSize==0)
684                {
685                        KK = 0;//dimension 1 (top)
686                        II = 0;//dimension 2 (mid)
687                        JJ++;
688                        realPrecision = pwrErrBound[JJ];
689                        updateReqLength = 0;                   
690                }               
691                type_ = type[jj];
692                if (type_ != 0)
693                {
694                        pred1D = 2*(*data)[jj-1] - (*data)[jj-2];
695                        (*data)[jj] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
696                }
697                else
698                {
699                        // compute resiBits
700                        if(updateReqLength==0)
701                        {
702                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
703                                reqBytesLength = reqLength/8;
704                                resiBitsLength = reqLength%8;                           
705                                updateReqLength = 1;
706                        }
707
708                        resiBits = 0;
709                        if (resiBitsLength != 0) {
710                                int kMod8 = k % 8;
711                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
712                                if (rightMovSteps > 0) {
713                                        int code = getRightMovingCode(kMod8, resiBitsLength);
714                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
715                                } else if (rightMovSteps < 0) {
716                                        int code1 = getLeftMovingCode(kMod8);
717                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
718                                        int leftMovSteps = -rightMovSteps;
719                                        rightMovSteps = 8 - leftMovSteps;
720                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
721                                        p++;
722                                        resiBits = resiBits
723                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
724                                } else // rightMovSteps == 0
725                                {
726                                        int code = getRightMovingCode(kMod8, resiBitsLength);
727                                        resiBits = (tdps->residualMidBits[p] & code);
728                                        p++;
729                                }
730                                k += resiBitsLength;
731                        }
732
733                        // recover the exact data
734                        memset(curBytes, 0, 4);
735                        leadingNum = leadNum[l++];
736                        memcpy(curBytes, preBytes, leadingNum);
737                        for (j = leadingNum; j < reqBytesLength; j++)
738                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
739                        if (resiBitsLength != 0) {
740                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
741                                curBytes[reqBytesLength] = resiByte;
742                        }
743
744                        exactData = bytesToFloat(curBytes);
745                        (*data)[jj] = exactData + medianValue;
746                        memcpy(preBytes,curBytes,4);
747                }
748        }
749        size_t index;
750        /* Process Row-1 --> Row-r2-1 */
751        for (ii = 1; ii < r2; ii++)
752        {
753                /* Process row-ii data 0 */             
754                if(ii%blockSize==0)
755                        II++;           
756                JJ = 0;
757                realPrecision = pwrErrBound[II*R3+JJ];
758                updateReqLength = 0;           
759
760                index = ii*r3;
761               
762                type_ = type[index];
763                if (type_ != 0)
764                {
765                        pred1D = (*data)[index-r3];                     
766                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
767                }
768                else
769                {
770                        // compute resiBits
771                        if(updateReqLength==0)
772                        {
773                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
774                                reqBytesLength = reqLength/8;
775                                resiBitsLength = reqLength%8;                           
776                                updateReqLength = 1;
777                        }
778                       
779                        resiBits = 0;
780                        if (resiBitsLength != 0) {
781                                int kMod8 = k % 8;
782                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
783                                if (rightMovSteps > 0) {
784                                        int code = getRightMovingCode(kMod8, resiBitsLength);
785                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
786                                } else if (rightMovSteps < 0) {
787                                        int code1 = getLeftMovingCode(kMod8);
788                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
789                                        int leftMovSteps = -rightMovSteps;
790                                        rightMovSteps = 8 - leftMovSteps;
791                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
792                                        p++;
793                                        resiBits = resiBits
794                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
795                                } else // rightMovSteps == 0
796                                {
797                                        int code = getRightMovingCode(kMod8, resiBitsLength);
798                                        resiBits = (tdps->residualMidBits[p] & code);
799                                        p++;
800                                }
801                                k += resiBitsLength;
802                        }
803
804                        // recover the exact data
805                        memset(curBytes, 0, 4);
806                        leadingNum = leadNum[l++];
807                        memcpy(curBytes, preBytes, leadingNum);
808                        for (j = leadingNum; j < reqBytesLength; j++)
809                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
810                        if (resiBitsLength != 0) {
811                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
812                                curBytes[reqBytesLength] = resiByte;
813                        }
814
815                        exactData = bytesToFloat(curBytes);
816                        (*data)[index] = exactData + medianValue;
817                        memcpy(preBytes,curBytes,4);
818                }
819
820                /* Process row-ii data 1 --> r3-1*/
821                for (jj = 1; jj < r3; jj++)
822                {
823                        index = ii*r3+jj;
824
825                        if(jj%blockSize==0)
826                                JJ++;
827                        realPrecision = pwrErrBound[II*R3+JJ];                 
828                        updateReqLength = 0;                   
829                       
830                        type_ = type[index];
831                        if (type_ != 0)
832                        {
833                                pred2D = (*data)[index-1] + (*data)[index-r3] - (*data)[index-r3-1];                           
834                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
835                        }
836                        else
837                        {
838                                // compute resiBits
839                                if(updateReqLength==0)
840                                {
841                                        computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
842                                        reqBytesLength = reqLength/8;
843                                        resiBitsLength = reqLength%8;                           
844                                        updateReqLength = 1;
845                                }
846
847                                resiBits = 0;
848                                if (resiBitsLength != 0) {
849                                        int kMod8 = k % 8;
850                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
851                                        if (rightMovSteps > 0) {
852                                                int code = getRightMovingCode(kMod8, resiBitsLength);
853                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
854                                        } else if (rightMovSteps < 0) {
855                                                int code1 = getLeftMovingCode(kMod8);
856                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
857                                                int leftMovSteps = -rightMovSteps;
858                                                rightMovSteps = 8 - leftMovSteps;
859                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
860                                                p++;
861                                                resiBits = resiBits
862                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
863                                        } else // rightMovSteps == 0
864                                        {
865                                                int code = getRightMovingCode(kMod8, resiBitsLength);
866                                                resiBits = (tdps->residualMidBits[p] & code);
867                                                p++;
868                                        }
869                                        k += resiBitsLength;
870                                }
871
872                                // recover the exact data
873                                memset(curBytes, 0, 4);
874                                leadingNum = leadNum[l++];
875                                memcpy(curBytes, preBytes, leadingNum);
876                                for (j = leadingNum; j < reqBytesLength; j++)
877                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
878                                if (resiBitsLength != 0) {
879                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
880                                        curBytes[reqBytesLength] = resiByte;
881                                }
882
883                                exactData = bytesToFloat(curBytes);
884                                (*data)[index] = exactData + medianValue;
885                                memcpy(preBytes,curBytes,4);
886                        }
887                }
888        }
889
890        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
891
892        for (kk = 1; kk < r1; kk++)
893        {
894                /* Process Row-0 data 0*/
895                index = kk*r23;         
896                if(kk%blockSize==0)
897                        KK++;
898                II = 0;
899                JJ = 0;
900
901                realPrecision = pwrErrBound[KK*R23];                   
902                updateReqLength = 0;                   
903
904                type_ = type[index];
905                if (type_ != 0)
906                {
907                        pred1D = (*data)[index-r23];                   
908                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
909                }
910                else
911                {
912                        // compute resiBits
913                        if(updateReqLength==0)
914                        {
915                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
916                                reqBytesLength = reqLength/8;
917                                resiBitsLength = reqLength%8;                           
918                                updateReqLength = 1;
919                        }
920
921                        resiBits = 0;
922                        if (resiBitsLength != 0) {
923                                int kMod8 = k % 8;
924                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
925                                if (rightMovSteps > 0) {
926                                        int code = getRightMovingCode(kMod8, resiBitsLength);
927                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
928                                } else if (rightMovSteps < 0) {
929                                        int code1 = getLeftMovingCode(kMod8);
930                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
931                                        int leftMovSteps = -rightMovSteps;
932                                        rightMovSteps = 8 - leftMovSteps;
933                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
934                                        p++;
935                                        resiBits = resiBits
936                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
937                                } else // rightMovSteps == 0
938                                {
939                                        int code = getRightMovingCode(kMod8, resiBitsLength);
940                                        resiBits = (tdps->residualMidBits[p] & code);
941                                        p++;
942                                }
943                                k += resiBitsLength;
944                        }
945
946                        // recover the exact data
947                        memset(curBytes, 0, 4);
948                        leadingNum = leadNum[l++];
949                        memcpy(curBytes, preBytes, leadingNum);
950                        for (j = leadingNum; j < reqBytesLength; j++)
951                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
952                        if (resiBitsLength != 0) {
953                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
954                                curBytes[reqBytesLength] = resiByte;
955                        }
956
957                        exactData = bytesToFloat(curBytes);
958                        (*data)[index] = exactData + medianValue;
959                        memcpy(preBytes,curBytes,4);
960                }
961
962                /* Process Row-0 data 1 --> data r3-1 */
963                for (jj = 1; jj < r3; jj++)
964                {
965                        index = kk*r23+jj;
966
967                        if(jj%blockSize==0)
968                                JJ++;
969
970                        realPrecision = pwrErrBound[KK*R23+JJ];                 
971                        updateReqLength = 0;                   
972                       
973                        type_ = type[index];
974                        if (type_ != 0)
975                        {
976                                pred2D = (*data)[index-1] + (*data)[index-r23] - (*data)[index-r23-1];                 
977                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
978                        }
979                        else
980                        {
981                                // compute resiBits
982                                if(updateReqLength==0)
983                                {
984                                        computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
985                                        reqBytesLength = reqLength/8;
986                                        resiBitsLength = reqLength%8;                           
987                                        updateReqLength = 1;
988                                }
989                       
990                                resiBits = 0;
991                                if (resiBitsLength != 0) {
992                                        int kMod8 = k % 8;
993                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
994                                        if (rightMovSteps > 0) {
995                                                int code = getRightMovingCode(kMod8, resiBitsLength);
996                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
997                                        } else if (rightMovSteps < 0) {
998                                                int code1 = getLeftMovingCode(kMod8);
999                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1000                                                int leftMovSteps = -rightMovSteps;
1001                                                rightMovSteps = 8 - leftMovSteps;
1002                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1003                                                p++;
1004                                                resiBits = resiBits
1005                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1006                                        } else // rightMovSteps == 0
1007                                        {
1008                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1009                                                resiBits = (tdps->residualMidBits[p] & code);
1010                                                p++;
1011                                        }
1012                                        k += resiBitsLength;
1013                                }
1014
1015                                // recover the exact data
1016                                memset(curBytes, 0, 4);
1017                                leadingNum = leadNum[l++];
1018                                memcpy(curBytes, preBytes, leadingNum);
1019                                for (j = leadingNum; j < reqBytesLength; j++)
1020                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1021                                if (resiBitsLength != 0) {
1022                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1023                                        curBytes[reqBytesLength] = resiByte;
1024                                }
1025
1026                                exactData = bytesToFloat(curBytes);
1027                                (*data)[index] = exactData + medianValue;
1028                                memcpy(preBytes,curBytes,4);
1029                        }
1030                }
1031
1032                /* Process Row-1 --> Row-r2-1 */
1033                for (ii = 1; ii < r2; ii++)
1034                {
1035                        /* Process Row-i data 0 */
1036                        index = kk*r23 + ii*r3;
1037                       
1038                        if(ii%blockSize==0)
1039                                II++;
1040                        JJ = 0;
1041                       
1042                        realPrecision = pwrErrBound[KK*R23+II*R3];                     
1043                        updateReqLength = 0;                                           
1044
1045                        type_ = type[index];
1046                        if (type_ != 0)
1047                        {
1048                                pred2D = (*data)[index-r3] + (*data)[index-r23] - (*data)[index-r23-r3];                               
1049                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1050                        }
1051                        else
1052                        {
1053                                // compute resiBits
1054                                if(updateReqLength==0)
1055                                {
1056                                        computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
1057                                        reqBytesLength = reqLength/8;
1058                                        resiBitsLength = reqLength%8;                           
1059                                        updateReqLength = 1;
1060                                }
1061
1062                                resiBits = 0;
1063                                if (resiBitsLength != 0) {
1064                                        int kMod8 = k % 8;
1065                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1066                                        if (rightMovSteps > 0) {
1067                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1068                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1069                                        } else if (rightMovSteps < 0) {
1070                                                int code1 = getLeftMovingCode(kMod8);
1071                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1072                                                int leftMovSteps = -rightMovSteps;
1073                                                rightMovSteps = 8 - leftMovSteps;
1074                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1075                                                p++;
1076                                                resiBits = resiBits
1077                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1078                                        } else // rightMovSteps == 0
1079                                        {
1080                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1081                                                resiBits = (tdps->residualMidBits[p] & code);
1082                                                p++;
1083                                        }
1084                                        k += resiBitsLength;
1085                                }
1086
1087                                // recover the exact data
1088                                memset(curBytes, 0, 4);
1089                                leadingNum = leadNum[l++];
1090                                memcpy(curBytes, preBytes, leadingNum);
1091                                for (j = leadingNum; j < reqBytesLength; j++)
1092                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1093                                if (resiBitsLength != 0) {
1094                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1095                                        curBytes[reqBytesLength] = resiByte;
1096                                }
1097
1098                                exactData = bytesToFloat(curBytes);
1099                                (*data)[index] = exactData + medianValue;
1100                                memcpy(preBytes,curBytes,4);
1101                        }
1102
1103                        /* Process Row-i data 1 --> data r3-1 */
1104                        for (jj = 1; jj < r3; jj++)
1105                        {
1106                                index = kk*r23 + ii*r3 + jj;
1107                                if(jj%blockSize==0)
1108                                        JJ++;
1109
1110                                realPrecision = pwrErrBound[KK*R23+II*R3+JJ];                   
1111                                updateReqLength = 0;                           
1112
1113                                type_ = type[index];
1114                                if (type_ != 0)
1115                                {
1116                                        pred3D = (*data)[index-1] + (*data)[index-r3] + (*data)[index-r23]
1117                                        - (*data)[index-r3-1] - (*data)[index-r23-r3] - (*data)[index-r23-1] + (*data)[index-r23-r3-1];                                 
1118                                        (*data)[index] = pred3D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1119                                }
1120                                else
1121                                {
1122                                        // compute resiBits
1123                                        if(updateReqLength==0)
1124                                        {
1125                                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
1126                                                reqBytesLength = reqLength/8;
1127                                                resiBitsLength = reqLength%8;                           
1128                                                updateReqLength = 1;
1129                                        }
1130                               
1131                                        resiBits = 0;
1132                                        if (resiBitsLength != 0) {
1133                                                int kMod8 = k % 8;
1134                                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1135                                                if (rightMovSteps > 0) {
1136                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1137                                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1138                                                } else if (rightMovSteps < 0) {
1139                                                        int code1 = getLeftMovingCode(kMod8);
1140                                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1141                                                        int leftMovSteps = -rightMovSteps;
1142                                                        rightMovSteps = 8 - leftMovSteps;
1143                                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1144                                                        p++;
1145                                                        resiBits = resiBits
1146                                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1147                                                } else // rightMovSteps == 0
1148                                                {
1149                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1150                                                        resiBits = (tdps->residualMidBits[p] & code);
1151                                                        p++;
1152                                                }
1153                                                k += resiBitsLength;
1154                                        }
1155
1156                                        // recover the exact data
1157                                        memset(curBytes, 0, 4);
1158                                        leadingNum = leadNum[l++];
1159                                        memcpy(curBytes, preBytes, leadingNum);
1160                                        for (j = leadingNum; j < reqBytesLength; j++)
1161                                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1162                                        if (resiBitsLength != 0) {
1163                                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1164                                                curBytes[reqBytesLength] = resiByte;
1165                                        }
1166                                       
1167                                        exactData = bytesToFloat(curBytes);
1168                                        (*data)[index] = exactData + medianValue;
1169                                        memcpy(preBytes,curBytes,4);
1170                                }
1171                        }
1172                }
1173
1174        }
1175
1176        free(pwrErrBound);
1177        free(leadNum);
1178        free(type);
1179        return;
1180}
1181
1182void decompressDataSeries_float_1D_pwrgroup(float** data, size_t dataSeriesLength, TightDataPointStorageF* tdps) 
1183{
1184        float *posGroups, *negGroups, *groups;
1185        float pos_01_group, neg_01_group;
1186        int *posFlags, *negFlags;
1187       
1188        updateQuantizationInfo(tdps->intervals);
1189       
1190        unsigned char* leadNum;
1191        double interval;// = (float)tdps->realPrecision*2;
1192       
1193        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
1194
1195        *data = (float*)malloc(sizeof(float)*dataSeriesLength);
1196
1197        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
1198
1199        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
1200        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
1201        SZ_ReleaseHuffman(huffmanTree); 
1202       
1203        createRangeGroups_float(&posGroups, &negGroups, &posFlags, &negFlags);
1204       
1205        float realGroupPrecision;
1206        float realPrecision = tdps->realPrecision;
1207        char* groupID = decompressGroupIDArray(tdps->pwrErrBoundBytes, tdps->dataSeriesLength);
1208       
1209        //note that the groupID values here are [1,2,3,....,18] or [-1,-2,...,-18]
1210       
1211        double* groupErrorBounds = generateGroupErrBounds(confparams_dec->errorBoundMode, realPrecision, confparams_dec->pw_relBoundRatio);
1212        exe_params->intvRadius = generateGroupMaxIntervalCount(groupErrorBounds);
1213               
1214        size_t nbBins = (size_t)(1/confparams_dec->pw_relBoundRatio + 0.5);
1215        if(nbBins%2==1)
1216                nbBins++;
1217        exe_params->intvRadius = nbBins;
1218
1219        unsigned char preBytes[4];
1220        unsigned char curBytes[4];
1221       
1222        memset(preBytes, 0, 4);
1223
1224        size_t curByteIndex = 0;
1225        int reqLength, reqBytesLength = 0, resiBitsLength = 0, resiBits; 
1226        unsigned char leadingNum;       
1227        float medianValue, exactData, curValue, predValue;
1228       
1229        medianValue = tdps->medianValue;
1230       
1231        size_t i, j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
1232                                                        // in resiMidBits, p is to track the
1233                                                        // byte_index of resiMidBits, l is for
1234                                                        // leadNum
1235                                                       
1236        int type_, updateReqLength = 0;
1237        char rawGrpID = 0, indexGrpID = 0;
1238        for (i = 0; i < dataSeriesLength; i++) 
1239        {
1240                rawGrpID = groupID[i];
1241               
1242                if(rawGrpID >= 2)
1243                {
1244                        groups = posGroups;
1245                        indexGrpID = rawGrpID - 2;
1246                }
1247                else if(rawGrpID <= -2)
1248                {
1249                        groups = negGroups;
1250                        indexGrpID = -rawGrpID - 2;             }
1251                else if(rawGrpID == 1)
1252                {
1253                        groups = &pos_01_group;
1254                        indexGrpID = 0;
1255                }
1256                else //rawGrpID == -1
1257                {
1258                        groups = &neg_01_group;
1259                        indexGrpID = 0;                 
1260                }
1261               
1262                type_ = type[i];
1263                switch (type_) {
1264                case 0:
1265                        // compute resiBits
1266                        if(updateReqLength==0)
1267                        {
1268                                computeReqLength_float(realPrecision, tdps->radExpo, &reqLength, &medianValue);
1269                                reqBytesLength = reqLength/8;
1270                                resiBitsLength = reqLength%8;   
1271                                updateReqLength = 1;   
1272                        }
1273                        resiBits = 0;
1274                        if (resiBitsLength != 0) {
1275                                int kMod8 = k % 8;
1276                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1277                                if (rightMovSteps > 0) {
1278                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1279                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1280                                } else if (rightMovSteps < 0) {
1281                                        int code1 = getLeftMovingCode(kMod8);
1282                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1283                                        int leftMovSteps = -rightMovSteps;
1284                                        rightMovSteps = 8 - leftMovSteps;
1285                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1286                                        p++;
1287                                        resiBits = resiBits
1288                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1289                                } else // rightMovSteps == 0
1290                                {
1291                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1292                                        resiBits = (tdps->residualMidBits[p] & code);
1293                                        p++;
1294                                }
1295                                k += resiBitsLength;
1296                        }
1297
1298                        // recover the exact data       
1299                        memset(curBytes, 0, 4);
1300                        leadingNum = leadNum[l++];
1301                        memcpy(curBytes, preBytes, leadingNum);
1302                        for (j = leadingNum; j < reqBytesLength; j++)
1303                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1304                        if (resiBitsLength != 0) {
1305                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1306                                curBytes[reqBytesLength] = resiByte;
1307                        }
1308                       
1309                        exactData = bytesToFloat(curBytes);
1310                        exactData = exactData + medianValue;
1311                        (*data)[i] = exactData;
1312                        memcpy(preBytes,curBytes,4);
1313                       
1314                        groups[indexGrpID] = exactData;
1315                       
1316                        break;
1317                default:
1318                        predValue = groups[indexGrpID]; //Here, groups[indexGrpID] is the previous value.
1319                        realGroupPrecision = groupErrorBounds[indexGrpID];
1320                        interval = realGroupPrecision*2;               
1321                       
1322                        curValue = predValue + (type_-exe_params->intvRadius)*interval;
1323                       
1324                        //groupNum = computeGroupNum_float(curValue);
1325                       
1326                        if((curValue>0&&rawGrpID<0)||(curValue<0&&rawGrpID>0))
1327                                curValue = 0;
1328                        //else
1329                        //{
1330                        //      realGrpID = fabs(rawGrpID)-2;
1331                        //      if(groupNum<realGrpID)
1332                        //              curValue = rawGrpID>0?pow(2,realGrpID):-pow(2,realGrpID);
1333                        //      else if(groupNum>realGrpID)
1334                        //              curValue = rawGrpID>0?pow(2,groupNum):-pow(2,groupNum);                         
1335                        //}     
1336                               
1337                        (*data)[i] = curValue;
1338                        groups[indexGrpID] = curValue;
1339                        break;         
1340                }
1341        }       
1342       
1343        free(leadNum);
1344        free(type);
1345       
1346        free(posGroups);
1347        free(negGroups);
1348        free(posFlags);
1349        free(negFlags);
1350        free(groupErrorBounds);
1351        free(groupID);
1352}
1353
1354void decompressDataSeries_float_1D_pwr_pre_log(float** data, size_t dataSeriesLength, TightDataPointStorageF* tdps) {
1355
1356        decompressDataSeries_float_1D(data, dataSeriesLength, tdps);
1357        float threshold = tdps->minLogValue;
1358        if(tdps->pwrErrBoundBytes_size > 0){
1359                unsigned char * signs;
1360                sz_lossless_decompress(confparams_dec->losslessCompressor, tdps->pwrErrBoundBytes, tdps->pwrErrBoundBytes_size, &signs, dataSeriesLength);
1361                for(size_t i=0; i<dataSeriesLength; i++){
1362                        if((*data)[i] < threshold) (*data)[i] = 0;
1363                        else (*data)[i] = exp2((*data)[i]);
1364                        if(signs[i]) (*data)[i] = -((*data)[i]);
1365                }
1366                free(signs);
1367        }
1368        else{
1369                for(size_t i=0; i<dataSeriesLength; i++){
1370                        if((*data)[i] < threshold) (*data)[i] = 0;
1371                        else (*data)[i] = exp2((*data)[i]);
1372                }
1373        }
1374
1375}
1376
1377void decompressDataSeries_float_2D_pwr_pre_log(float** data, size_t r1, size_t r2, TightDataPointStorageF* tdps) {
1378
1379        size_t dataSeriesLength = r1 * r2;
1380        decompressDataSeries_float_2D(data, r1, r2, tdps);
1381        float threshold = tdps->minLogValue;
1382        if(tdps->pwrErrBoundBytes_size > 0){
1383                unsigned char * signs;
1384                sz_lossless_decompress(confparams_dec->losslessCompressor, tdps->pwrErrBoundBytes, tdps->pwrErrBoundBytes_size, &signs, dataSeriesLength);
1385                for(size_t i=0; i<dataSeriesLength; i++){
1386                        if((*data)[i] < threshold) (*data)[i] = 0;
1387                        else (*data)[i] = exp2((*data)[i]);
1388                        if(signs[i]) (*data)[i] = -((*data)[i]);
1389                }
1390                free(signs);
1391        }
1392        else{
1393                for(size_t i=0; i<dataSeriesLength; i++){
1394                        if((*data)[i] < threshold) (*data)[i] = 0;
1395                        else (*data)[i] = exp2((*data)[i]);
1396                }
1397        }
1398
1399}
1400
1401void decompressDataSeries_float_3D_pwr_pre_log(float** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageF* tdps) {
1402
1403        size_t dataSeriesLength = r1 * r2 * r3;
1404        decompressDataSeries_float_3D(data, r1, r2, r3, tdps);
1405        float threshold = tdps->minLogValue;
1406        if(tdps->pwrErrBoundBytes_size > 0){
1407                unsigned char * signs;
1408                sz_lossless_decompress(confparams_dec->losslessCompressor, tdps->pwrErrBoundBytes, tdps->pwrErrBoundBytes_size, &signs, dataSeriesLength);
1409                for(size_t i=0; i<dataSeriesLength; i++){
1410                        if((*data)[i] < threshold) (*data)[i] = 0;
1411                        else (*data)[i] = exp2((*data)[i]);
1412                        if(signs[i]) (*data)[i] = -((*data)[i]);
1413                }
1414                free(signs);
1415        }
1416        else{
1417                for(size_t i=0; i<dataSeriesLength; i++){
1418                        if((*data)[i] < threshold) (*data)[i] = 0;
1419                        else (*data)[i] = exp2((*data)[i]);
1420                }
1421        }
1422}
1423#pragma GCC diagnostic pop
Note: See TracBrowser for help on using the repository browser.