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

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

importing new SZ files
  • Property mode set to 100644
Line 
1/**
2 *  @file double_compression.c
3 *  @author Sheng Di
4 *  @date April, 2016
5 *  @brief Compression Technique for double array
6 *  (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7 *      See COPYRIGHT in top-level directory.
8 */
9
10#include <stdio.h>
11#include <stdlib.h>
12#include <string.h>
13#include <unistd.h>
14#include "sz.h"
15#include "DynamicByteArray.h"
16#include "DynamicIntArray.h"
17#include "TightDataPointStorageD.h"
18#include "CompressElement.h"
19#include "dataCompression.h"
20
21int computeByteSizePerIntValue(long valueRangeSize)
22{
23        if(valueRangeSize<=256)
24                return 1;
25        else if(valueRangeSize<=65536)
26                return 2;
27        else if(valueRangeSize<=4294967296) //2^32
28                return 4;
29        else
30                return 8;
31}
32
33long computeRangeSize_int(void* oriData, int dataType, size_t size, int64_t* valueRangeSize)
34{
35        size_t i = 0;
36        long max = 0, min = 0;
37
38        if(dataType==SZ_UINT8)
39        {
40                unsigned char* data = (unsigned char*)oriData;
41                unsigned char data_; 
42                min = data[0], max = min;
43                computeMinMax(data);
44        }
45        else if(dataType == SZ_INT8)
46        {
47                char* data = (char*)oriData;
48                char data_;
49                min = data[0], max = min;
50                computeMinMax(data);
51        }
52        else if(dataType == SZ_UINT16)
53        {
54                unsigned short* data = (unsigned short*)oriData;
55                unsigned short data_; 
56                min = data[0], max = min;
57                computeMinMax(data);
58        }
59        else if(dataType == SZ_INT16)
60        { 
61                short* data = (short*)oriData;
62                short data_; 
63                min = data[0], max = min;
64                computeMinMax(data);
65        }
66        else if(dataType == SZ_UINT32)
67        {
68                unsigned int* data = (unsigned int*)oriData;
69                unsigned int data_; 
70                min = data[0], max = min;
71                computeMinMax(data);
72        }
73        else if(dataType == SZ_INT32)
74        {
75                int* data = (int*)oriData;
76                int data_; 
77                min = data[0], max = min;
78                computeMinMax(data);
79        }
80        else if(dataType == SZ_UINT64)
81        {
82                unsigned long* data = (unsigned long*)oriData;
83                unsigned long data_; 
84                min = data[0], max = min;
85                computeMinMax(data);
86        }
87        else if(dataType == SZ_INT64)
88        {
89                long* data = (long *)oriData;
90                long data_; 
91                min = data[0], max = min;
92                computeMinMax(data);
93        }
94
95        *valueRangeSize = max - min;
96        return min;     
97}
98
99float computeRangeSize_float(float* oriData, size_t size, float* valueRangeSize, float* medianValue)
100{
101        size_t i = 0;
102        float min = oriData[0];
103        float max = min;
104        for(i=1;i<size;i++)
105        {
106                float data = oriData[i];
107                if(min>data)
108                        min = data;
109                else if(max<data)
110                        max = data;
111        }
112
113        *valueRangeSize = max - min;
114        *medianValue = min + *valueRangeSize/2;
115        return min;
116}
117
118double computeRangeSize_double(double* oriData, size_t size, double* valueRangeSize, double* medianValue)
119{
120        size_t i = 0;
121        double min = oriData[0];
122        double max = min;
123        for(i=1;i<size;i++)
124        {
125                double data = oriData[i];
126                if(min>data)
127                        min = data;
128                else if(max<data)
129                        max = data;
130        }
131       
132        *valueRangeSize = max - min;
133        *medianValue = min + *valueRangeSize/2;
134        return min;
135}
136
137float computeRangeSize_float_subblock(float* oriData, float* valueRangeSize, float* medianValue,
138size_t r5, size_t r4, size_t r3, size_t r2, size_t r1,
139size_t s5, size_t s4, size_t s3, size_t s2, size_t s1,
140size_t e5, size_t e4, size_t e3, size_t e2, size_t e1)
141{
142        size_t i1, i2, i3, i4, i5;
143        size_t index_start = s5*(r4*r3*r2*r1) + s4*(r3*r2*r1) + s3*(r2*r1) + s2*r1 + s1;
144        float min = oriData[index_start];
145        float max = min;
146
147        for (i5 = s5; i5 <= e5; i5++)
148        for (i4 = s4; i4 <= e4; i4++)
149        for (i3 = s3; i3 <= e3; i3++)
150        for (i2 = s2; i2 <= e2; i2++)
151        for (i1 = s1; i1 <= e1; i1++)
152        {
153                size_t index = i5*(r4*r3*r2*r1) + i4*(r3*r2*r1) + i3*(r2*r1) + i2*r1 + i1;
154                float data = oriData[index];
155                if (min>data)
156                        min = data;
157                else if(max<data)
158                        max = data;
159        }
160
161        *valueRangeSize = max - min;
162        *medianValue = min + *valueRangeSize/2;
163        return min;
164}
165
166
167float computeRangeSize_double_subblock(double* oriData, double* valueRangeSize, double* medianValue,
168size_t r5, size_t r4, size_t r3, size_t r2, size_t r1,
169size_t s5, size_t s4, size_t s3, size_t s2, size_t s1,
170size_t e5, size_t e4, size_t e3, size_t e2, size_t e1)
171{
172        size_t i1, i2, i3, i4, i5;
173        size_t index_start = s5*(r4*r3*r2*r1) + s4*(r3*r2*r1) + s3*(r2*r1) + s2*r1 + s1;
174        double min = oriData[index_start];
175        double max = min;
176
177        for (i5 = s5; i5 <= e5; i5++)
178        for (i4 = s4; i4 <= e4; i4++)
179        for (i3 = s3; i3 <= e3; i3++)
180        for (i2 = s2; i2 <= e2; i2++)
181        for (i1 = s1; i1 <= e1; i1++)
182        {
183                size_t index = i5*(r4*r3*r2*r1) + i4*(r3*r2*r1) + i3*(r2*r1) + i2*r1 + i1;
184                double data = oriData[index];
185                if (min>data)
186                        min = data;
187                else if(max<data)
188                        max = data;
189        }
190
191        *valueRangeSize = max - min;
192        *medianValue = min + *valueRangeSize/2;
193        return min;
194}
195
196
197double min_d(double a, double b)
198{
199        if(a<b)
200                return a;
201        else
202                return b;
203}
204
205double max_d(double a, double b)
206{
207        if(a>b)
208                return a;
209        else
210                return b;
211}
212
213float min_f(float a, float b)
214{
215        if(a<b)
216                return a;
217        else
218                return b;
219}
220
221float max_f(float a, float b)
222{
223        if(a>b)
224                return a;
225        else
226                return b;
227}
228
229double getRealPrecision_double(double valueRangeSize, int errBoundMode, double absErrBound, double relBoundRatio, int *status)
230{
231        int state = SZ_SCES;
232        double precision = 0;
233        if(errBoundMode==ABS||errBoundMode==ABS_OR_PW_REL||errBoundMode==ABS_AND_PW_REL)
234                precision = absErrBound; 
235        else if(errBoundMode==REL||errBoundMode==REL_OR_PW_REL||errBoundMode==REL_AND_PW_REL)
236                precision = relBoundRatio*valueRangeSize;
237        else if(errBoundMode==ABS_AND_REL)
238                precision = min_d(absErrBound, relBoundRatio*valueRangeSize);
239        else if(errBoundMode==ABS_OR_REL)
240                precision = max_d(absErrBound, relBoundRatio*valueRangeSize);
241        else if(errBoundMode==PW_REL)
242                precision = 0;
243        else
244        {
245                printf("Error: error-bound-mode is incorrect!\n");
246                state = SZ_BERR;
247        }
248        *status = state;
249        return precision;
250}
251
252double getRealPrecision_float(float valueRangeSize, int errBoundMode, double absErrBound, double relBoundRatio, int *status)
253{
254        int state = SZ_SCES;
255        double precision = 0;
256        if(errBoundMode==ABS||errBoundMode==ABS_OR_PW_REL||errBoundMode==ABS_AND_PW_REL)
257                precision = absErrBound; 
258        else if(errBoundMode==REL||errBoundMode==REL_OR_PW_REL||errBoundMode==REL_AND_PW_REL)
259                precision = relBoundRatio*valueRangeSize;
260        else if(errBoundMode==ABS_AND_REL)
261                precision = min_f(absErrBound, relBoundRatio*valueRangeSize);
262        else if(errBoundMode==ABS_OR_REL)
263                precision = max_f(absErrBound, relBoundRatio*valueRangeSize);
264        else if(errBoundMode==PW_REL)
265                precision = 0;
266        else
267        {
268                printf("Error: error-bound-mode is incorrect!\n");
269                state = SZ_BERR;
270        }
271        *status = state;
272        return precision;
273}
274
275double getRealPrecision_int(long valueRangeSize, int errBoundMode, double absErrBound, double relBoundRatio, int *status)
276{
277        int state = SZ_SCES;
278        double precision = 0;
279        if(errBoundMode==ABS||errBoundMode==ABS_OR_PW_REL||errBoundMode==ABS_AND_PW_REL)
280                precision = absErrBound; 
281        else if(errBoundMode==REL||errBoundMode==REL_OR_PW_REL||errBoundMode==REL_AND_PW_REL)
282                precision = relBoundRatio*valueRangeSize;
283        else if(errBoundMode==ABS_AND_REL)
284                precision = min_f(absErrBound, relBoundRatio*valueRangeSize);
285        else if(errBoundMode==ABS_OR_REL)
286                precision = max_f(absErrBound, relBoundRatio*valueRangeSize);
287        else if(errBoundMode==PW_REL)
288                precision = -1;
289        else
290        {
291                printf("Error: error-bound-mode is incorrect!\n");
292                state = SZ_BERR;
293        }
294        *status = state;
295        return precision;
296}
297
298void symTransform_8bytes(unsigned char data[8])
299{
300        unsigned char tmp = data[0];
301        data[0] = data[7];
302        data[7] = tmp;
303
304        tmp = data[1];
305        data[1] = data[6];
306        data[6] = tmp;
307       
308        tmp = data[2];
309        data[2] = data[5];
310        data[5] = tmp;
311       
312        tmp = data[3];
313        data[3] = data[4];
314        data[4] = tmp;
315}
316
317inline void symTransform_2bytes(unsigned char data[2])
318{
319        unsigned char tmp = data[0];
320        data[0] = data[1];
321        data[1] = tmp;
322}
323
324inline void symTransform_4bytes(unsigned char data[4])
325{
326        unsigned char tmp = data[0];
327        data[0] = data[3];
328        data[3] = tmp;
329
330        tmp = data[1];
331        data[1] = data[2];
332        data[2] = tmp;
333}
334
335inline void compressInt8Value(int8_t tgtValue, int8_t minValue, int byteSize, unsigned char* bytes)
336{
337        uint8_t data = tgtValue - minValue;
338        memcpy(bytes, &data, byteSize); //byteSize==1
339}
340
341inline void compressInt16Value(int16_t tgtValue, int16_t minValue, int byteSize, unsigned char* bytes)
342{
343        uint16_t data = tgtValue - minValue;
344        unsigned char tmpBytes[2];
345        int16ToBytes_bigEndian(tmpBytes, data);
346        memcpy(bytes, tmpBytes + 2 - byteSize, byteSize);
347}
348
349inline void compressInt32Value(int32_t tgtValue, int32_t minValue, int byteSize, unsigned char* bytes)
350{
351        uint32_t data = tgtValue - minValue;
352        unsigned char tmpBytes[4];
353        int32ToBytes_bigEndian(tmpBytes, data);
354        memcpy(bytes, tmpBytes + 4 - byteSize, byteSize);
355}
356
357inline void compressInt64Value(int64_t tgtValue, int64_t minValue, int byteSize, unsigned char* bytes)
358{
359        uint64_t data = tgtValue - minValue;
360        unsigned char tmpBytes[8];
361        int64ToBytes_bigEndian(tmpBytes, data);
362        memcpy(bytes, tmpBytes + 8 - byteSize, byteSize);
363}
364
365inline void compressUInt8Value(uint8_t tgtValue, uint8_t minValue, int byteSize, unsigned char* bytes)
366{
367        uint8_t data = tgtValue - minValue;
368        memcpy(bytes, &data, byteSize); //byteSize==1
369}
370
371inline void compressUInt16Value(uint16_t tgtValue, uint16_t minValue, int byteSize, unsigned char* bytes)
372{
373        uint16_t data = tgtValue - minValue;
374        unsigned char tmpBytes[2];
375        int16ToBytes_bigEndian(tmpBytes, data);
376        memcpy(bytes, tmpBytes + 2 - byteSize, byteSize);
377}
378
379inline void compressUInt32Value(uint32_t tgtValue, uint32_t minValue, int byteSize, unsigned char* bytes)
380{
381        uint32_t data = tgtValue - minValue;
382        unsigned char tmpBytes[4];
383        int32ToBytes_bigEndian(tmpBytes, data);
384        memcpy(bytes, tmpBytes + 4 - byteSize, byteSize);
385}
386
387inline void compressUInt64Value(uint64_t tgtValue, uint64_t minValue, int byteSize, unsigned char* bytes)
388{
389        uint64_t data = tgtValue - minValue;
390        unsigned char tmpBytes[8];
391        int64ToBytes_bigEndian(tmpBytes, data);
392        memcpy(bytes, tmpBytes + 8 - byteSize, byteSize);
393}
394
395void compressSingleFloatValue(FloatValueCompressElement *vce, float tgtValue, float precision, float medianValue, 
396                int reqLength, int reqBytesLength, int resiBitsLength)
397{               
398        float normValue = tgtValue - medianValue;
399
400        lfloat lfBuf;
401        lfBuf.value = normValue;
402                       
403        int ignBytesLength = 32 - reqLength;
404        if(ignBytesLength<0)
405                ignBytesLength = 0;
406       
407        int tmp_int = lfBuf.ivalue;
408        intToBytes_bigEndian(vce->curBytes, tmp_int);
409               
410        lfBuf.ivalue = (lfBuf.ivalue >> ignBytesLength) << ignBytesLength;
411       
412        //float tmpValue = lfBuf.value;
413       
414        vce->data = lfBuf.value+medianValue;
415        vce->curValue = tmp_int;
416        vce->reqBytesLength = reqBytesLength;
417        vce->resiBitsLength = resiBitsLength;
418}
419
420void compressSingleDoubleValue(DoubleValueCompressElement *vce, double tgtValue, double precision, double medianValue, 
421                int reqLength, int reqBytesLength, int resiBitsLength)
422{               
423        double normValue = tgtValue - medianValue;
424
425        ldouble lfBuf;
426        lfBuf.value = normValue;
427                       
428        int ignBytesLength = 64 - reqLength;
429        if(ignBytesLength<0)
430                ignBytesLength = 0;
431
432        long tmp_long = lfBuf.lvalue;
433        longToBytes_bigEndian(vce->curBytes, tmp_long);
434                               
435        lfBuf.lvalue = (lfBuf.lvalue >> ignBytesLength)<<ignBytesLength;
436       
437        //double tmpValue = lfBuf.value;
438       
439        vce->data = lfBuf.value+medianValue;
440        vce->curValue = tmp_long;
441        vce->reqBytesLength = reqBytesLength;
442        vce->resiBitsLength = resiBitsLength;
443}
444
445int compIdenticalLeadingBytesCount_double(unsigned char* preBytes, unsigned char* curBytes)
446{
447        int i, n = 0;
448        for(i=0;i<8;i++)
449                if(preBytes[i]==curBytes[i])
450                        n++;
451                else
452                        break;
453        if(n>3) n = 3;
454        return n;
455}
456
457int compIdenticalLeadingBytesCount_float(unsigned char* preBytes, unsigned char* curBytes)
458{
459        int i, n = 0;
460        for(i=0;i<4;i++)
461                if(preBytes[i]==curBytes[i])
462                        n++;
463                else
464                        break;
465        if(n>3) n = 3;
466        return n;
467}
468
469//TODO double-check the correctness...
470void addExactData(DynamicByteArray *exactMidByteArray, DynamicIntArray *exactLeadNumArray, 
471                DynamicIntArray *resiBitArray, LossyCompressionElement *lce)
472{
473        int i;
474        int leadByteLength = lce->leadingZeroBytes;
475        addDIA_Data(exactLeadNumArray, leadByteLength);
476        unsigned char* intMidBytes = lce->integerMidBytes;
477        int integerMidBytesLength = lce->integerMidBytes_Length;
478        int resMidBitsLength = lce->resMidBitsLength;
479        if(intMidBytes!=NULL||resMidBitsLength!=0)
480        {
481                if(intMidBytes!=NULL)
482                        for(i = 0;i<integerMidBytesLength;i++)
483                                addDBA_Data(exactMidByteArray, intMidBytes[i]);
484                if(resMidBitsLength!=0)
485                        addDIA_Data(resiBitArray, lce->residualMidBits);
486        }
487}
488
489/**
490 * @deprecated
491 * @return: the length of the coefficient array.
492 * */
493int getPredictionCoefficients(int layers, int dimension, int **coeff_array, int *status)
494{
495        size_t size = 0;
496        switch(dimension)
497        {
498                case 1:
499                        switch(layers)
500                        {
501                                case 1:
502                                        *coeff_array = (int*)malloc(sizeof(int));
503                                        (*coeff_array)[0] = 1;
504                                        size = 1;
505                                        break;
506                                case 2:
507                                        *coeff_array = (int*)malloc(2*sizeof(int));
508                                        (*coeff_array)[0] = 2;
509                                        (*coeff_array)[1] = -1;
510                                        size = 2;
511                                        break;
512                                case 3:
513                                        *coeff_array = (int*)malloc(3*sizeof(int));
514                                        (*coeff_array)[0] = 3;
515                                        (*coeff_array)[1] = -3;
516                                        (*coeff_array)[2] = 1;
517                                        break;
518                        }       
519                        break;
520                case 2:
521                        switch(layers)
522                        {
523                                case 1:
524                               
525                                        break;
526                                case 2:
527                               
528                                        break;
529                                case 3:
530                               
531                                        break;
532                        }                               
533                        break;
534                case 3:
535                        switch(layers)
536                        {
537                                case 1:
538                               
539                                        break;
540                                case 2:
541                               
542                                        break;
543                                case 3:
544                               
545                                        break;
546                        }                       
547                        break;
548                default:
549                        printf("Error: dimension must be no greater than 3 in the current version.\n");
550                        *status = SZ_DERR;
551        }
552        *status = SZ_SCES;
553        return size;
554}
555
556int computeBlockEdgeSize_2D(int segmentSize)
557{
558        int i = 1;
559        for(i=1; i<segmentSize;i++)
560        {
561                if(i*i>segmentSize)
562                        break;
563        }
564        return i;
565        //return (int)(sqrt(segmentSize)+1);
566}
567
568int computeBlockEdgeSize_3D(int segmentSize)
569{
570        int i = 1;
571        for(i=1; i<segmentSize;i++)
572        {
573                if(i*i*i>segmentSize)
574                        break;
575        }
576        return i;       
577        //return (int)(pow(segmentSize, 1.0/3)+1);
578}
579
580//convert random-access version based bytes to output bytes
581int initRandomAccessBytes(unsigned char* raBytes)
582{
583        int k = 0, i = 0;
584        for (i = 0; i < 3; i++)//3
585                raBytes[k++] = versionNumber[i];
586        int sameByte = 0x80; //indicating this is random-access mode
587        if(exe_params->SZ_SIZE_TYPE==8)
588                sameByte = (unsigned char) (sameByte | 0x40); // 01000000, the 6th bit
589        sameByte = sameByte | (confparams_cpr->szMode << 1);
590
591        raBytes[k++] = sameByte;
592
593        convertSZParamsToBytes(confparams_cpr, &(raBytes[k]));
594        k = k + MetaDataByteLength;
595
596        return k;
597}
598
599//The following functions are float-precision version of dealing with the unpredictable data points
600int generateLossyCoefficients_float(float* oriData, double precision, size_t nbEle, int* reqBytesLength, int* resiBitsLength, float* medianValue, float* decData)
601{
602        float valueRangeSize;
603       
604        computeRangeSize_float(oriData, nbEle, &valueRangeSize, medianValue);
605        short radExpo = getExponent_float(valueRangeSize/2);
606       
607        int reqLength;
608        computeReqLength_float(precision, radExpo, &reqLength, medianValue);
609       
610        *reqBytesLength = reqLength/8;
611        *resiBitsLength = reqLength%8;
612       
613        size_t i = 0;
614        for(i = 0;i < nbEle;i++)
615        {
616                float normValue = oriData[i] - *medianValue;
617
618                lfloat lfBuf;
619                lfBuf.value = normValue;
620                               
621                int ignBytesLength = 32 - reqLength;
622                if(ignBytesLength<0)
623                        ignBytesLength = 0;
624                       
625                lfBuf.ivalue = (lfBuf.ivalue >> ignBytesLength) << ignBytesLength;
626               
627                //float tmpValue = lfBuf.value;
628               
629                decData[i] = lfBuf.value + *medianValue;
630        }
631        return reqLength;
632}       
633               
634/**
635 * @param float* oriData: inplace argument (input / output)
636 *
637 * */           
638int compressExactDataArray_float(float* oriData, double precision, size_t nbEle, unsigned char** leadArray, unsigned char** midArray, unsigned char** resiArray, 
639int reqLength, int reqBytesLength, int resiBitsLength, float medianValue)
640{
641        //allocate memory for coefficient compression arrays
642        DynamicIntArray *exactLeadNumArray;
643        new_DIA(&exactLeadNumArray, DynArrayInitLen);   
644        DynamicByteArray *exactMidByteArray;
645        new_DBA(&exactMidByteArray, DynArrayInitLen);
646        DynamicIntArray *resiBitArray;
647        new_DIA(&resiBitArray, DynArrayInitLen);
648        unsigned char preDataBytes[4] = {0,0,0,0};     
649
650        //allocate memory for vce and lce
651        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
652        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));       
653
654        size_t i = 0;
655        for(i = 0;i < nbEle;i++)
656        {
657                compressSingleFloatValue(vce, oriData[i], precision, medianValue, reqLength, reqBytesLength, resiBitsLength);
658                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
659                memcpy(preDataBytes,vce->curBytes,4);
660                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
661                oriData[i] = vce->data;
662        }
663        convertDIAtoInts(exactLeadNumArray, leadArray);
664        convertDBAtoBytes(exactMidByteArray,midArray);
665        convertDIAtoInts(resiBitArray, resiArray);
666
667        size_t midArraySize = exactMidByteArray->size;
668       
669        free(vce);
670        free(lce);
671       
672        free_DIA(exactLeadNumArray);
673        free_DBA(exactMidByteArray);
674        free_DIA(resiBitArray);
675       
676        return midArraySize;
677}
678
679void decompressExactDataArray_float(unsigned char* leadNum, unsigned char* exactMidBytes, unsigned char* residualMidBits, size_t nbEle, int reqLength, float medianValue, float** decData)
680{
681        *decData = (float*)malloc(nbEle*sizeof(float));
682        size_t i = 0, j = 0, k = 0, l = 0, p = 0, curByteIndex = 0;
683        float exactData = 0;
684        unsigned char preBytes[4] = {0,0,0,0};
685        unsigned char curBytes[4];
686        int resiBits; 
687        unsigned char leadingNum;               
688       
689        int reqBytesLength = reqLength/8;
690        int resiBitsLength = reqLength%8;
691       
692        for(i = 0; i<nbEle;i++)
693        {
694                // compute resiBits
695                resiBits = 0;
696                if (resiBitsLength != 0) {
697                        int kMod8 = k % 8;
698                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
699                        if (rightMovSteps > 0) {
700                                int code = getRightMovingCode(kMod8, resiBitsLength);
701                                resiBits = (residualMidBits[p] & code) >> rightMovSteps;
702                        } else if (rightMovSteps < 0) {
703                                int code1 = getLeftMovingCode(kMod8);
704                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
705                                int leftMovSteps = -rightMovSteps;
706                                rightMovSteps = 8 - leftMovSteps;
707                                resiBits = (residualMidBits[p] & code1) << leftMovSteps;
708                                p++;
709                                resiBits = resiBits
710                                                | ((residualMidBits[p] & code2) >> rightMovSteps);
711                        } else // rightMovSteps == 0
712                        {
713                                int code = getRightMovingCode(kMod8, resiBitsLength);
714                                resiBits = (residualMidBits[p] & code);
715                                p++;
716                        }
717                        k += resiBitsLength;
718                }
719
720                // recover the exact data       
721                memset(curBytes, 0, 4);
722                leadingNum = leadNum[l++];
723                memcpy(curBytes, preBytes, leadingNum);
724                for (j = leadingNum; j < reqBytesLength; j++)
725                        curBytes[j] = exactMidBytes[curByteIndex++];
726                if (resiBitsLength != 0) {
727                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
728                        curBytes[reqBytesLength] = resiByte;
729                }
730
731                exactData = bytesToFloat(curBytes);
732                (*decData)[i] = exactData + medianValue;
733                memcpy(preBytes,curBytes,4);
734        }       
735}
736
737//double-precision version of dealing with unpredictable data points in sz 2.0
738int generateLossyCoefficients_double(double* oriData, double precision, size_t nbEle, int* reqBytesLength, int* resiBitsLength, double* medianValue, double* decData)
739{
740        double valueRangeSize;
741       
742        computeRangeSize_double(oriData, nbEle, &valueRangeSize, medianValue);
743        short radExpo = getExponent_double(valueRangeSize/2);
744       
745        int reqLength;
746        computeReqLength_double(precision, radExpo, &reqLength, medianValue);
747       
748        *reqBytesLength = reqLength/8;
749        *resiBitsLength = reqLength%8;
750       
751        size_t i = 0;
752        for(i = 0;i < nbEle;i++)
753        {
754                double normValue = oriData[i] - *medianValue;
755
756                ldouble ldBuf;
757                ldBuf.value = normValue;
758                               
759                int ignBytesLength = 64 - reqLength;
760                if(ignBytesLength<0)
761                        ignBytesLength = 0;
762                       
763                ldBuf.lvalue = (ldBuf.lvalue >> ignBytesLength) << ignBytesLength;
764               
765                decData[i] = ldBuf.value + *medianValue;
766        }
767        return reqLength;
768}       
769               
770/**
771 * @param double* oriData: inplace argument (input / output)
772 *
773 * */           
774int compressExactDataArray_double(double* oriData, double precision, size_t nbEle, unsigned char** leadArray, unsigned char** midArray, unsigned char** resiArray, 
775int reqLength, int reqBytesLength, int resiBitsLength, double medianValue)
776{
777        //allocate memory for coefficient compression arrays
778        DynamicIntArray *exactLeadNumArray;
779        new_DIA(&exactLeadNumArray, DynArrayInitLen);   
780        DynamicByteArray *exactMidByteArray;
781        new_DBA(&exactMidByteArray, DynArrayInitLen);
782        DynamicIntArray *resiBitArray;
783        new_DIA(&resiBitArray, DynArrayInitLen);
784        unsigned char preDataBytes[8] = {0,0,0,0,0,0,0,0};     
785
786        //allocate memory for vce and lce
787        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
788        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));       
789
790        size_t i = 0;
791        for(i = 0;i < nbEle;i++)
792        {
793                compressSingleDoubleValue(vce, oriData[i], precision, medianValue, reqLength, reqBytesLength, resiBitsLength);
794                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
795                memcpy(preDataBytes,vce->curBytes,8);
796                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
797                oriData[i] = vce->data;
798        }
799        convertDIAtoInts(exactLeadNumArray, leadArray);
800        convertDBAtoBytes(exactMidByteArray,midArray);
801        convertDIAtoInts(resiBitArray, resiArray);
802
803        size_t midArraySize = exactMidByteArray->size;
804       
805        free(vce);
806        free(lce);
807       
808        free_DIA(exactLeadNumArray);
809        free_DBA(exactMidByteArray);
810        free_DIA(resiBitArray);
811       
812        return midArraySize;
813}
814
815void decompressExactDataArray_double(unsigned char* leadNum, unsigned char* exactMidBytes, unsigned char* residualMidBits, size_t nbEle, int reqLength, double medianValue, double** decData)
816{
817        *decData = (double*)malloc(nbEle*sizeof(double));
818        size_t i = 0, j = 0, k = 0, l = 0, p = 0, curByteIndex = 0;
819        double exactData = 0;
820        unsigned char preBytes[8] = {0,0,0,0,0,0,0,0};
821        unsigned char curBytes[8];
822        int resiBits; 
823        unsigned char leadingNum;               
824       
825        int reqBytesLength = reqLength/8;
826        int resiBitsLength = reqLength%8;
827       
828        for(i = 0; i<nbEle;i++)
829        {
830                // compute resiBits
831                resiBits = 0;
832                if (resiBitsLength != 0) {
833                        int kMod8 = k % 8;
834                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
835                        if (rightMovSteps > 0) {
836                                int code = getRightMovingCode(kMod8, resiBitsLength);
837                                resiBits = (residualMidBits[p] & code) >> rightMovSteps;
838                        } else if (rightMovSteps < 0) {
839                                int code1 = getLeftMovingCode(kMod8);
840                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
841                                int leftMovSteps = -rightMovSteps;
842                                rightMovSteps = 8 - leftMovSteps;
843                                resiBits = (residualMidBits[p] & code1) << leftMovSteps;
844                                p++;
845                                resiBits = resiBits
846                                                | ((residualMidBits[p] & code2) >> rightMovSteps);
847                        } else // rightMovSteps == 0
848                        {
849                                int code = getRightMovingCode(kMod8, resiBitsLength);
850                                resiBits = (residualMidBits[p] & code);
851                                p++;
852                        }
853                        k += resiBitsLength;
854                }
855
856                // recover the exact data       
857                memset(curBytes, 0, 8);
858                leadingNum = leadNum[l++];
859                memcpy(curBytes, preBytes, leadingNum);
860                for (j = leadingNum; j < reqBytesLength; j++)
861                        curBytes[j] = exactMidBytes[curByteIndex++];
862                if (resiBitsLength != 0) {
863                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
864                        curBytes[reqBytesLength] = resiByte;
865                }
866
867                exactData = bytesToDouble(curBytes);
868                (*decData)[i] = exactData + medianValue;
869                memcpy(preBytes,curBytes,8);
870        }
871}
Note: See TracBrowser for help on using the repository browser.