source: thirdparty/SZ/sz/src/szd_float_pwr.c @ 2c47b73

Revision 2c47b73, 39.5 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)

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