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

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

more work on adding SZ (latest version)
  • 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                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                unsigned 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}
Note: See TracBrowser for help on using the repository browser.