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

Revision 2c47b73, 40.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 sz_uint16.c
3 *  @author Sheng Di
4 *  @date Aug, 2017
5 *  @brief sz_uint16, Compression and Decompression functions
6 *  (C) 2017 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7 *      See COPYRIGHT in top-level directory.
8 */
9
10
11#include <stdio.h>
12#include <stdlib.h>
13#include <string.h>
14#include <unistd.h>
15#include <math.h>
16#include "sz.h"
17#include "CompressElement.h"
18#include "DynamicByteArray.h"
19#include "DynamicIntArray.h"
20#include "zlib.h"
21#include "rw.h"
22#include "TightDataPointStorageI.h"
23#include "sz_uint16.h"
24
25unsigned int optimize_intervals_uint16_1D(uint16_t *oriData, size_t dataLength, double realPrecision)
26{       
27        size_t i = 0, radiusIndex;
28        int64_t pred_value = 0, pred_err;
29        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
30        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
31        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
32        for(i=2;i<dataLength;i++)
33        {
34                if(i%confparams_cpr->sampleDistance==0)
35                {
36                        //pred_value = 2*oriData[i-1] - oriData[i-2];
37                        pred_value = oriData[i-1];
38                        pred_err = llabs(pred_value - oriData[i]);
39                        radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
40                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
41                                radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
42                        intervals[radiusIndex]++;
43                }
44        }
45        //compute the appropriate number
46        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
47        size_t sum = 0;
48        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
49        {
50                sum += intervals[i];
51                if(sum>targetCount)
52                        break;
53        }
54        if(i>=confparams_cpr->maxRangeRadius)
55                i = confparams_cpr->maxRangeRadius-1;
56               
57        unsigned int accIntervals = 2*(i+1);
58        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
59       
60        if(powerOf2<32)
61                powerOf2 = 32;
62       
63        free(intervals);
64        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
65        return powerOf2;
66}
67
68unsigned int optimize_intervals_uint16_2D(uint16_t *oriData, size_t r1, size_t r2, double realPrecision)
69{       
70        size_t i,j, index;
71        size_t radiusIndex;
72        int64_t pred_value = 0, pred_err;
73        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
74        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
75        size_t totalSampleSize = (r1-1)*(r2-1)/confparams_cpr->sampleDistance;
76        for(i=1;i<r1;i++)
77        {
78                for(j=1;j<r2;j++)
79                {
80                        if((i+j)%confparams_cpr->sampleDistance==0)
81                        {
82                                index = i*r2+j;
83                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
84                                pred_err = llabs(pred_value - oriData[index]);
85                                radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
86                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
87                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
88                                intervals[radiusIndex]++;
89                        }                       
90                }
91        }
92        //compute the appropriate number
93        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
94        size_t sum = 0;
95        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
96        {
97                sum += intervals[i];
98                if(sum>targetCount)
99                        break;
100        }
101        if(i>=confparams_cpr->maxRangeRadius)
102                i = confparams_cpr->maxRangeRadius-1;
103        unsigned int accIntervals = 2*(i+1);
104        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
105
106        if(powerOf2<32)
107                powerOf2 = 32;
108
109        free(intervals);
110        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
111        return powerOf2;
112}
113
114unsigned int optimize_intervals_uint16_3D(uint16_t *oriData, size_t r1, size_t r2, size_t r3, double realPrecision)
115{       
116        size_t i,j,k, index;
117        size_t radiusIndex;
118        size_t r23=r2*r3;
119        int64_t pred_value = 0, pred_err;
120        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
121        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
122        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)/confparams_cpr->sampleDistance;
123        for(i=1;i<r1;i++)
124        {
125                for(j=1;j<r2;j++)
126                {
127                        for(k=1;k<r3;k++)
128                        {                       
129                                if((i+j+k)%confparams_cpr->sampleDistance==0)
130                                {
131                                        index = i*r23+j*r3+k;
132                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23] 
133                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
134                                        pred_err = llabs(pred_value - oriData[index]);
135                                        radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
136                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
137                                        {
138                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
139                                                //printf("radiusIndex=%d\n", radiusIndex);
140                                        }
141                                        intervals[radiusIndex]++;
142                                }
143                        }
144                }
145        }
146        //compute the appropriate number
147        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
148        size_t sum = 0;
149        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
150        {
151                sum += intervals[i];
152                if(sum>targetCount)
153                        break;
154        }
155        if(i>=confparams_cpr->maxRangeRadius)
156                i = confparams_cpr->maxRangeRadius-1;
157        unsigned int accIntervals = 2*(i+1);
158        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
159
160        if(powerOf2<32)
161                powerOf2 = 32;
162       
163        free(intervals);
164        //printf("targetCount=%d, sum=%d, totalSampleSize=%d, ratio=%f, accIntervals=%d, powerOf2=%d\n", targetCount, sum, totalSampleSize, (double)sum/(double)totalSampleSize, accIntervals, powerOf2);
165        return powerOf2;
166}
167
168
169unsigned int optimize_intervals_uint16_4D(uint16_t *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision)
170{
171        size_t i,j,k,l, index;
172        size_t radiusIndex;
173        size_t r234=r2*r3*r4;
174        size_t r34=r3*r4;
175        int64_t pred_value = 0, pred_err;
176        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
177        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
178        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)*(r4-1)/confparams_cpr->sampleDistance;
179        for(i=1;i<r1;i++)
180        {
181                for(j=1;j<r2;j++)
182                {
183                        for(k=1;k<r3;k++)
184                        {
185                                for (l=1;l<r4;l++)
186                                {
187                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
188                                        {
189                                                index = i*r234+j*r34+k*r4+l;
190                                                pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r34]
191                                                                - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
192                                                pred_err = llabs(pred_value - oriData[index]);
193                                                radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
194                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
195                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
196                                                intervals[radiusIndex]++;
197                                        }
198                                }
199                        }
200                }
201        }
202        //compute the appropriate number
203        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
204        size_t sum = 0;
205        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
206        {
207                sum += intervals[i];
208                if(sum>targetCount)
209                        break;
210        }
211        if(i>=confparams_cpr->maxRangeRadius)
212                i = confparams_cpr->maxRangeRadius-1;
213
214        unsigned int accIntervals = 2*(i+1);
215        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
216
217        if(powerOf2<32)
218                powerOf2 = 32;
219
220        free(intervals);
221        return powerOf2;
222}
223
224TightDataPointStorageI* SZ_compress_uint16_1D_MDQ(uint16_t *oriData, size_t dataLength, double realPrecision, int64_t valueRangeSize, int64_t minValue)
225{
226        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
227        int byteSize = computeByteSizePerIntValue(valueRangeSize);
228        unsigned int quantization_intervals;
229        if(exe_params->optQuantMode==1)
230                quantization_intervals = optimize_intervals_uint16_1D(oriData, dataLength, realPrecision);
231        else
232                quantization_intervals = exe_params->intvCapacity;
233        updateQuantizationInfo(quantization_intervals); 
234        size_t i;
235
236        int* type = (int*) malloc(dataLength*sizeof(int));
237               
238        uint16_t* spaceFillingValue = oriData; //
239       
240        DynamicByteArray *exactDataByteArray;
241        new_DBA(&exactDataByteArray, DynArrayInitLen);
242               
243        int64_t last3CmprsData[3] = {0,0,0};
244                               
245        //add the first data   
246        type[0] = 0;
247        compressUInt16Value(spaceFillingValue[0], minValue, byteSize, bytes);
248        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
249        listAdd_int(last3CmprsData, spaceFillingValue[0]);
250               
251        type[1] = 0;
252        compressUInt16Value(spaceFillingValue[1], minValue, byteSize, bytes);
253        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
254        listAdd_int(last3CmprsData, spaceFillingValue[1]);
255        //printf("%.30G\n",last3CmprsData[0]); 
256       
257        int state;
258        double checkRadius = (exe_params->intvCapacity-1)*realPrecision;
259        int64_t curData;
260        int64_t pred, predAbsErr;
261        double interval = 2*realPrecision;
262       
263        for(i=2;i<dataLength;i++)
264        {
265                curData = spaceFillingValue[i];
266                //pred = 2*last3CmprsData[0] - last3CmprsData[1];
267                pred = last3CmprsData[0];
268                predAbsErr = llabs(curData - pred);     
269                if(predAbsErr<=checkRadius)
270                {
271                        state = (predAbsErr/realPrecision+1)/2;
272                        if(curData>=pred)
273                        {
274                                type[i] = exe_params->intvRadius+state;
275                                pred = pred + state*interval;
276                        }
277                        else //curData<pred
278                        {
279                                type[i] = exe_params->intvRadius-state;
280                                pred = pred - state*interval;
281                        }
282                        if(pred>SZ_UINT16_MAX) pred = SZ_UINT16_MAX;
283                        if(pred<SZ_UINT16_MIN) pred = SZ_UINT16_MIN;                   
284                        listAdd_int(last3CmprsData, pred);                                     
285                        continue;
286                }
287               
288                //unpredictable data processing         
289                type[i] = 0;
290                compressUInt16Value(curData, minValue, byteSize, bytes);
291                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
292                listAdd_int(last3CmprsData, curData);
293        }//end of for
294               
295        size_t exactDataNum = exactDataByteArray->size / byteSize;
296       
297        TightDataPointStorageI* tdps;   
298                       
299        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
300                        type, exactDataByteArray->array, exactDataByteArray->size, 
301                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
302
303//sdi:Debug
304/*      int sum =0;
305        for(i=0;i<dataLength;i++)
306                if(type[i]==0) sum++;
307        printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);*/
308       
309        //free memory
310        free(type);     
311        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
312       
313        return tdps;
314}
315
316void SZ_compress_args_uint16_StoreOriData(uint16_t* oriData, size_t dataLength, TightDataPointStorageI* tdps, 
317unsigned char** newByteData, size_t *outSize)
318{
319        int intSize=sizeof(uint16_t);   
320        size_t k = 0, i;
321        tdps->isLossless = 1;
322        size_t totalByteLength = 3 + MetaDataByteLength + exe_params->SZ_SIZE_TYPE + 1 + intSize*dataLength;
323        *newByteData = (unsigned char*)malloc(totalByteLength);
324       
325        unsigned char dsLengthBytes[8];
326        for (i = 0; i < 3; i++)//3
327                (*newByteData)[k++] = versionNumber[i]; 
328
329        if(exe_params->SZ_SIZE_TYPE==4)//1
330                (*newByteData)[k++] = 16; //00010000
331        else
332                (*newByteData)[k++] = 80;       //01010000: 01000000 indicates the SZ_SIZE_TYPE=8
333       
334        convertSZParamsToBytes(confparams_cpr, &((*newByteData)[k]));
335        k = k + MetaDataByteLength;     
336       
337        sizeToBytes(dsLengthBytes,dataLength); //SZ_SIZE_TYPE: 4 or 8   
338        for (i = 0; i < exe_params->SZ_SIZE_TYPE; i++)
339                (*newByteData)[k++] = dsLengthBytes[i];
340               
341        if(sysEndianType==BIG_ENDIAN_SYSTEM)
342                memcpy((*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, oriData, dataLength*intSize);
343        else
344        {
345                unsigned char* p = (*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
346                for(i=0;i<dataLength;i++,p+=intSize)
347                        int16ToBytes_bigEndian(p, oriData[i]);
348        }       
349        *outSize = totalByteLength;
350}
351
352void SZ_compress_args_uint16_NoCkRngeNoGzip_1D(unsigned char** newByteData, uint16_t *oriData, 
353size_t dataLength, double realPrecision, size_t *outSize, int64_t valueRangeSize, uint16_t minValue)
354{
355        TightDataPointStorageI* tdps = SZ_compress_uint16_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, minValue);
356        //TODO: return bytes....
357        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
358        if(*outSize > dataLength*sizeof(uint16_t))
359                SZ_compress_args_uint16_StoreOriData(oriData, dataLength+2, tdps, newByteData, outSize);
360        free_TightDataPointStorageI(tdps);
361}
362
363TightDataPointStorageI* SZ_compress_uint16_2D_MDQ(uint16_t *oriData, size_t r1, size_t r2, double realPrecision, int64_t valueRangeSize, int64_t minValue)
364{
365        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
366        int byteSize = computeByteSizePerIntValue(valueRangeSize);
367       
368        unsigned int quantization_intervals;
369        if(exe_params->optQuantMode==1)
370        {
371                quantization_intervals = optimize_intervals_uint16_2D(oriData, r1, r2, realPrecision);
372                updateQuantizationInfo(quantization_intervals);
373        }       
374        else
375                quantization_intervals = exe_params->intvCapacity;
376        size_t i,j; 
377        int64_t pred1D, pred2D, curValue, tmp;
378        int diff = 0.0;
379        double itvNum = 0;
380        uint16_t *P0, *P1;
381               
382        size_t dataLength = r1*r2;     
383       
384        P0 = (uint16_t*)malloc(r2*sizeof(uint16_t));
385        memset(P0, 0, r2*sizeof(uint16_t));
386        P1 = (uint16_t*)malloc(r2*sizeof(uint16_t));
387        memset(P1, 0, r2*sizeof(uint16_t));
388               
389        int* type = (int*) malloc(dataLength*sizeof(int));
390        //type[dataLength]=0;
391               
392        uint16_t* spaceFillingValue = oriData; //
393       
394        DynamicByteArray *exactDataByteArray;
395        new_DBA(&exactDataByteArray, DynArrayInitLen); 
396
397        type[0] = 0;
398        curValue = P1[0] = spaceFillingValue[0];
399        compressUInt16Value(curValue, minValue, byteSize, bytes);
400        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
401
402        /* Process Row-0 data 1*/
403        pred1D = P1[0];
404        diff = spaceFillingValue[1] - pred1D;
405
406        itvNum =  llabs(diff)/realPrecision + 1;
407
408        if (itvNum < exe_params->intvCapacity)
409        {
410                if (diff < 0) itvNum = -itvNum;
411                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
412                tmp = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
413                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
414                        P1[1] = tmp;
415                else if(tmp < SZ_UINT16_MIN)
416                        P1[1] = SZ_UINT16_MIN;
417                else
418                        P1[1] = SZ_UINT16_MAX;
419        }
420        else
421        {
422                type[1] = 0;
423                curValue = P1[1] = spaceFillingValue[1];
424                compressUInt16Value(curValue, minValue, byteSize, bytes);
425                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
426        }
427
428    /* Process Row-0 data 2 --> data r2-1 */
429        for (j = 2; j < r2; j++)
430        {
431                pred1D = 2*P1[j-1] - P1[j-2];
432                diff = spaceFillingValue[j] - pred1D;
433
434                itvNum = llabs(diff)/realPrecision + 1;
435
436                if (itvNum < exe_params->intvCapacity)
437                {
438                        if (diff < 0) itvNum = -itvNum;
439                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
440                        tmp = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
441                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
442                                P1[j] = tmp;
443                        else if(tmp < SZ_UINT16_MIN)
444                                P1[j] = SZ_UINT16_MIN;
445                        else
446                                P1[j] = SZ_UINT16_MAX;                 
447                }
448                else
449                {
450                        type[j] = 0;
451                        curValue = P1[j] = spaceFillingValue[j];
452                        compressUInt16Value(curValue, minValue, byteSize, bytes);
453                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
454                }
455        }
456
457        /* Process Row-1 --> Row-r1-1 */
458        size_t index;
459        for (i = 1; i < r1; i++)
460        {       
461                /* Process row-i data 0 */
462                index = i*r2;
463                pred1D = P1[0];
464                diff = spaceFillingValue[index] - pred1D;
465
466                itvNum = llabs(diff)/realPrecision + 1;
467
468                if (itvNum < exe_params->intvCapacity)
469                {
470                        if (diff < 0) itvNum = -itvNum;
471                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
472                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
473                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
474                                P0[0] = tmp;
475                        else if(tmp < SZ_UINT16_MIN)
476                                P0[0] = SZ_UINT16_MIN;
477                        else
478                                P0[0] = SZ_UINT16_MAX;                 
479                }
480                else
481                {
482                        type[index] = 0;
483                        curValue = P0[0] = spaceFillingValue[index];
484                        compressUInt16Value(curValue, minValue, byteSize, bytes);
485                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
486                }
487                                                                       
488                /* Process row-i data 1 --> r2-1*/
489                for (j = 1; j < r2; j++)
490                {
491                        index = i*r2+j;
492                        pred2D = P0[j-1] + P1[j] - P1[j-1];
493
494                        diff = spaceFillingValue[index] - pred2D;
495
496                        itvNum = llabs(diff)/realPrecision + 1;
497
498                        if (itvNum < exe_params->intvCapacity)
499                        {
500                                if (diff < 0) itvNum = -itvNum;
501                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
502                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
503                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
504                                        P0[j] = tmp;
505                                else if(tmp < SZ_UINT16_MIN)
506                                        P0[j] = SZ_UINT16_MIN;
507                                else
508                                        P0[j] = SZ_UINT16_MAX;                                         
509                        }
510                        else
511                        {
512                                type[index] = 0;
513                                curValue = P0[j] = spaceFillingValue[index];
514                                compressUInt16Value(curValue, minValue, byteSize, bytes);
515                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
516                        }
517                }
518
519                uint16_t *Pt;
520                Pt = P1;
521                P1 = P0;
522                P0 = Pt;
523        }
524       
525        if(r2!=1)
526                free(P0);
527        free(P1);                       
528       
529        size_t exactDataNum = exactDataByteArray->size;
530       
531        TightDataPointStorageI* tdps;   
532                       
533        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
534                        type, exactDataByteArray->array, exactDataByteArray->size, 
535                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
536                       
537        //free memory
538        free(type);     
539        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
540       
541        return tdps;   
542}
543
544/**
545 *
546 * Note: @r1 is high dimension
547 *               @r2 is low dimension
548 * */
549void SZ_compress_args_uint16_NoCkRngeNoGzip_2D(unsigned char** newByteData, uint16_t *oriData, size_t r1, size_t r2, double realPrecision, size_t *outSize, 
550int64_t valueRangeSize, uint16_t minValue)
551{
552        TightDataPointStorageI* tdps = SZ_compress_uint16_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, minValue);
553
554        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
555
556        size_t dataLength = r1*r2;
557        if(*outSize>dataLength*sizeof(uint16_t))
558                SZ_compress_args_uint16_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
559       
560        free_TightDataPointStorageI(tdps);     
561}
562
563TightDataPointStorageI* SZ_compress_uint16_3D_MDQ(uint16_t *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, int64_t valueRangeSize, int64_t minValue)
564{
565        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
566        int byteSize = computeByteSizePerIntValue(valueRangeSize);
567       
568        unsigned int quantization_intervals;
569        if(exe_params->optQuantMode==1)
570        {
571                quantization_intervals = optimize_intervals_uint16_3D(oriData, r1, r2, r3, realPrecision);
572                updateQuantizationInfo(quantization_intervals);
573        }       
574        else
575                quantization_intervals = exe_params->intvCapacity;
576        size_t i,j,k; 
577        int64_t pred1D, pred2D, pred3D, curValue, tmp;
578        int diff = 0.0;
579        double itvNum = 0;
580        uint16_t *P0, *P1;
581               
582        size_t dataLength = r1*r2*r3;           
583
584        size_t r23 = r2*r3;
585        P0 = (uint16_t*)malloc(r23*sizeof(uint16_t));
586        P1 = (uint16_t*)malloc(r23*sizeof(uint16_t));
587
588        int* type = (int*) malloc(dataLength*sizeof(int));
589
590        uint16_t* spaceFillingValue = oriData; //
591       
592        DynamicByteArray *exactDataByteArray;
593        new_DBA(&exactDataByteArray, DynArrayInitLen); 
594
595        type[0] = 0;
596        P1[0] = spaceFillingValue[0];
597        compressUInt16Value(spaceFillingValue[0], minValue, byteSize, bytes);
598        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
599
600        /* Process Row-0 data 1*/
601        pred1D = P1[0];
602        diff = spaceFillingValue[1] - pred1D;
603
604        itvNum = llabs(diff)/realPrecision + 1;
605
606        if (itvNum < exe_params->intvCapacity)
607        {
608                if (diff < 0) itvNum = -itvNum;
609                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
610                tmp = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
611                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
612                        P1[1] = tmp;
613                else if(tmp < SZ_UINT16_MIN)
614                        P1[1] = SZ_UINT16_MIN;
615                else
616                        P1[1] = SZ_UINT16_MAX;         
617        }
618        else
619        {
620                type[1] = 0;
621                curValue = P1[1] = spaceFillingValue[1];
622                compressUInt16Value(curValue, minValue, byteSize, bytes);
623                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
624        }
625
626    /* Process Row-0 data 2 --> data r3-1 */
627        for (j = 2; j < r3; j++)
628        {
629                pred1D = 2*P1[j-1] - P1[j-2];
630                diff = spaceFillingValue[j] - pred1D;
631
632                itvNum = llabs(diff)/realPrecision + 1;
633
634                if (itvNum < exe_params->intvCapacity)
635                {
636                        if (diff < 0) itvNum = -itvNum;
637                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
638                        tmp = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
639                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
640                                P1[j] = tmp;
641                        else if(tmp < SZ_UINT16_MIN)
642                                P1[j] = SZ_UINT16_MIN;
643                        else
644                                P1[j] = SZ_UINT16_MAX;                 
645                }
646                else
647                {
648                        type[j] = 0;
649                        curValue = P1[j] = spaceFillingValue[j];
650                        compressUInt16Value(curValue, minValue, byteSize, bytes);
651                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
652                }
653        }
654
655        /* Process Row-1 --> Row-r2-1 */
656        size_t index;
657        for (i = 1; i < r2; i++)
658        {
659                /* Process row-i data 0 */
660                index = i*r3;   
661                pred1D = P1[index-r3];
662                diff = spaceFillingValue[index] - pred1D;
663
664                itvNum = llabs(diff)/realPrecision + 1;
665
666                if (itvNum < exe_params->intvCapacity)
667                {
668                        if (diff < 0) itvNum = -itvNum;
669                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
670                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
671                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
672                                P1[index] = tmp;
673                        else if(tmp < SZ_UINT16_MIN)
674                                P1[index] = SZ_UINT16_MIN;
675                        else
676                                P1[index] = SZ_UINT16_MAX;                     
677                }
678                else
679                {
680                        type[index] = 0;
681                        curValue = P1[index] = spaceFillingValue[index];
682                        compressUInt16Value(curValue, minValue, byteSize, bytes);
683                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
684                }
685
686                /* Process row-i data 1 --> data r3-1*/
687                for (j = 1; j < r3; j++)
688                {
689                        index = i*r3+j;
690                        pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];
691
692                        diff = spaceFillingValue[index] - pred2D;
693
694                        itvNum = llabs(diff)/realPrecision + 1;
695
696                        if (itvNum < exe_params->intvCapacity)
697                        {
698                                if (diff < 0) itvNum = -itvNum;
699                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
700                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
701                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
702                                        P1[index] = tmp;
703                                else if(tmp < SZ_UINT16_MIN)
704                                        P1[index] = SZ_UINT16_MIN;
705                                else
706                                        P1[index] = SZ_UINT16_MAX;                             
707                        }
708                        else
709                        {
710                                type[index] = 0;
711                                curValue = P1[index] = spaceFillingValue[index];
712                                compressUInt16Value(curValue, minValue, byteSize, bytes);
713                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
714                        }
715                }
716        }
717
718
719        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
720
721        for (k = 1; k < r1; k++)
722        {
723                /* Process Row-0 data 0*/
724                index = k*r23;
725                pred1D = P1[0];
726                diff = spaceFillingValue[index] - pred1D;
727
728                itvNum = llabs(diff)/realPrecision + 1;
729
730                if (itvNum < exe_params->intvCapacity)
731                {
732                        if (diff < 0) itvNum = -itvNum;
733                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
734                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
735                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
736                                P0[0] = tmp;
737                        else if(tmp < SZ_UINT16_MIN)
738                                P0[0] = SZ_UINT16_MIN;
739                        else
740                                P0[0] = SZ_UINT16_MAX;
741                }
742                else
743                {
744                        type[index] = 0;
745                        curValue = P0[0] = spaceFillingValue[index];
746                        compressUInt16Value(curValue, minValue, byteSize, bytes);
747                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
748                }
749
750
751            /* Process Row-0 data 1 --> data r3-1 */
752                for (j = 1; j < r3; j++)
753                {
754                        //index = k*r2*r3+j;
755                        index ++;
756                        pred2D = P0[j-1] + P1[j] - P1[j-1];
757                        diff = spaceFillingValue[index] - pred2D;
758
759                        itvNum = llabs(diff)/realPrecision + 1;
760
761                        if (itvNum < exe_params->intvCapacity)
762                        {
763                                if (diff < 0) itvNum = -itvNum;
764                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
765                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
766                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
767                                        P0[j] = tmp;
768                                else if(tmp < SZ_UINT16_MIN)
769                                        P0[j] = SZ_UINT16_MIN;
770                                else
771                                        P0[j] = SZ_UINT16_MAX;                         
772                        }
773                        else
774                        {
775                                type[index] = 0;
776                                curValue = P0[j] = spaceFillingValue[index];
777                                compressUInt16Value(curValue, minValue, byteSize, bytes);
778                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
779                        }
780                }
781
782            /* Process Row-1 --> Row-r2-1 */
783                size_t index2D;
784                for (i = 1; i < r2; i++)
785                {
786                        /* Process Row-i data 0 */
787                        index = k*r23 + i*r3;
788                        index2D = i*r3;         
789                        pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
790                        diff = spaceFillingValue[index] - pred2D;
791
792                        itvNum = llabs(diff)/realPrecision + 1;
793
794                        if (itvNum < exe_params->intvCapacity)
795                        {
796                                if (diff < 0) itvNum = -itvNum;
797                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
798                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
799                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
800                                        P0[index2D] = tmp;
801                                else if(tmp < SZ_UINT16_MIN)
802                                        P0[index2D] = SZ_UINT16_MIN;
803                                else
804                                        P0[index2D] = SZ_UINT16_MAX;
805                        }
806                        else
807                        {
808                                type[index] = 0;
809                                curValue = P0[index2D] = spaceFillingValue[index];
810                                compressUInt16Value(curValue, minValue, byteSize, bytes);
811                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
812                        }
813
814                        /* Process Row-i data 1 --> data r3-1 */
815                        for (j = 1; j < r3; j++)
816                        {
817//                              if(k==63&&i==43&&j==27)
818//                                      printf("i=%d\n", i);
819                                //index = k*r2*r3 + i*r3 + j;                   
820                                index ++;
821                                index2D = i*r3 + j;
822                                pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
823                                diff = spaceFillingValue[index] - pred3D;
824
825                                itvNum = llabs(diff)/realPrecision + 1;
826
827                                if (itvNum < exe_params->intvCapacity)
828                                {
829                                        if (diff < 0) itvNum = -itvNum;
830                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
831                                        tmp = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
832                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
833                                                P0[index2D] = tmp;
834                                        else if(tmp < SZ_UINT16_MIN)
835                                                P0[index2D] = SZ_UINT16_MIN;
836                                        else
837                                                P0[index2D] = SZ_UINT16_MAX;
838                                }
839                                else
840                                {
841                                        type[index] = 0;
842                                        curValue = P0[index2D] = spaceFillingValue[index];
843                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
844                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
845                                }
846                        }
847                }
848
849                uint16_t *Pt;
850                Pt = P1;
851                P1 = P0;
852                P0 = Pt;
853        }
854        if(r23!=1)
855                free(P0);
856        free(P1);
857
858        size_t exactDataNum = exactDataByteArray->size;
859       
860        TightDataPointStorageI* tdps;   
861                       
862        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
863                        type, exactDataByteArray->array, exactDataByteArray->size, 
864                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
865                       
866        //free memory
867        free(type);     
868        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
869       
870        return tdps;   
871}
872
873
874void SZ_compress_args_uint16_NoCkRngeNoGzip_3D(unsigned char** newByteData, uint16_t *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, size_t *outSize, 
875int64_t valueRangeSize, int64_t minValue)
876{       
877        TightDataPointStorageI* tdps = SZ_compress_uint16_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, minValue);
878
879        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
880
881        size_t dataLength = r1*r2*r3;
882        if(*outSize>dataLength*sizeof(uint16_t))
883                SZ_compress_args_uint16_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
884       
885        free_TightDataPointStorageI(tdps);     
886}
887
888
889TightDataPointStorageI* SZ_compress_uint16_4D_MDQ(uint16_t *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, int64_t valueRangeSize, int64_t minValue)
890{
891        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
892        int byteSize = computeByteSizePerIntValue(valueRangeSize);
893       
894        unsigned int quantization_intervals;
895        if(exe_params->optQuantMode==1)
896        {
897                quantization_intervals = optimize_intervals_uint16_4D(oriData, r1, r2, r3, r4, realPrecision);
898                updateQuantizationInfo(quantization_intervals);
899        }       
900        else
901                quantization_intervals = exe_params->intvCapacity;
902        size_t i,j,k; 
903        int64_t pred1D, pred2D, pred3D, curValue, tmp;
904        int diff = 0.0;
905        double itvNum = 0;
906        uint16_t *P0, *P1;
907               
908        size_t dataLength = r1*r2*r3*r4;               
909
910        size_t r234 = r2*r3*r4;
911        size_t r34 = r3*r4;
912
913        P0 = (uint16_t*)malloc(r34*sizeof(uint16_t));
914        P1 = (uint16_t*)malloc(r34*sizeof(uint16_t));
915       
916        int* type = (int*) malloc(dataLength*sizeof(int));
917
918        uint16_t* spaceFillingValue = oriData; //
919       
920        DynamicByteArray *exactDataByteArray;
921        new_DBA(&exactDataByteArray, DynArrayInitLen); 
922
923        size_t l;
924        for (l = 0; l < r1; l++)
925        {
926
927                ///////////////////////////     Process layer-0 ///////////////////////////
928                /* Process Row-0 data 0*/
929                size_t index = l*r234;
930                size_t index2D = 0;
931
932                type[index] = 0;
933                curValue = P1[index2D] = spaceFillingValue[index];
934                compressUInt16Value(curValue, minValue, byteSize, bytes);
935                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
936
937                /* Process Row-0 data 1*/
938                index = l*r234+1;
939                index2D = 1;
940
941                pred1D = P1[index2D-1];
942                diff = curValue - pred1D;
943
944                itvNum = llabs(diff)/realPrecision + 1;
945
946                if (itvNum < exe_params->intvCapacity)
947                {
948                        if (diff < 0) itvNum = -itvNum;
949                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
950                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
951                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
952                                P1[index2D] = tmp;
953                        else if(tmp < SZ_UINT16_MIN)
954                                P1[index2D] = SZ_UINT16_MIN;
955                        else
956                                P1[index2D] = SZ_UINT16_MAX;                   
957                }
958                else
959                {
960                        type[index] = 0;
961
962                        curValue = P1[index2D] = spaceFillingValue[0];
963                        compressUInt16Value(curValue, minValue, byteSize, bytes);
964                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
965                }
966
967                /* Process Row-0 data 2 --> data r4-1 */
968                for (j = 2; j < r4; j++)
969                {
970                        index = l*r234+j;
971                        index2D = j;
972
973                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
974                        diff = spaceFillingValue[index] - pred1D;
975
976                        itvNum = llabs(diff)/realPrecision + 1;
977
978                        if (itvNum < exe_params->intvCapacity)
979                        {
980                                if (diff < 0) itvNum = -itvNum;
981                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
982                                tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
983                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
984                                        P1[index2D] = tmp;
985                                else if(tmp < SZ_UINT16_MIN)
986                                        P1[index2D] = SZ_UINT16_MIN;
987                                else
988                                        P1[index2D] = SZ_UINT16_MAX;                                   
989                        }
990                        else
991                        {
992                                type[index] = 0;
993
994                                curValue = P1[index2D] = spaceFillingValue[0];
995                                compressUInt16Value(curValue, minValue, byteSize, bytes);
996                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
997                        }
998                }
999
1000                /* Process Row-1 --> Row-r3-1 */
1001                for (i = 1; i < r3; i++)
1002                {
1003                        /* Process row-i data 0 */
1004                        index = l*r234+i*r4;
1005                        index2D = i*r4;
1006
1007                        pred1D = P1[index2D-r4];
1008                        diff = spaceFillingValue[index] - pred1D;
1009
1010                        itvNum = llabs(diff)/realPrecision + 1;
1011
1012                        if (itvNum < exe_params->intvCapacity)
1013                        {
1014                                if (diff < 0) itvNum = -itvNum;
1015                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1016                                tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1017                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1018                                        P1[index2D] = tmp;
1019                                else if(tmp < SZ_UINT16_MIN)
1020                                        P1[index2D] = SZ_UINT16_MIN;
1021                                else
1022                                        P1[index2D] = SZ_UINT16_MAX;                                   
1023                        }
1024                        else
1025                        {
1026                                type[index] = 0;
1027
1028                                curValue = P1[index2D] = spaceFillingValue[0];
1029                                compressUInt16Value(curValue, minValue, byteSize, bytes);
1030                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1031                        }
1032
1033                        /* Process row-i data 1 --> data r4-1*/
1034                        for (j = 1; j < r4; j++)
1035                        {
1036                                index = l*r234+i*r4+j;
1037                                index2D = i*r4+j;
1038
1039                                pred2D = P1[index2D-1] + P1[index2D-r4] - P1[index2D-r4-1];
1040
1041                                diff = spaceFillingValue[index] - pred2D;
1042
1043                                itvNum = llabs(diff)/realPrecision + 1;
1044
1045                                if (itvNum < exe_params->intvCapacity)
1046                                {
1047                                        if (diff < 0) itvNum = -itvNum;
1048                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1049                                        tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1050                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1051                                                P1[index2D] = tmp;
1052                                        else if(tmp < SZ_UINT16_MIN)
1053                                                P1[index2D] = SZ_UINT16_MIN;
1054                                        else
1055                                                P1[index2D] = SZ_UINT16_MAX;                                           
1056                                }
1057                                else
1058                                {
1059                                        type[index] = 0;
1060
1061                                        curValue = P1[index2D] = spaceFillingValue[0];
1062                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
1063                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1064                                }
1065                        }
1066                }
1067
1068
1069                ///////////////////////////     Process layer-1 --> layer-r2-1 ///////////////////////////
1070
1071                for (k = 1; k < r2; k++)
1072                {
1073                        /* Process Row-0 data 0*/
1074                        index = l*r234+k*r34;
1075                        index2D = 0;
1076
1077                        pred1D = P1[index2D];
1078                        diff = spaceFillingValue[index] - pred1D;
1079
1080                        itvNum = llabs(diff)/realPrecision + 1;
1081
1082                        if (itvNum < exe_params->intvCapacity)
1083                        {
1084                                if (diff < 0) itvNum = -itvNum;
1085                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1086                                tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1087                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1088                                        P0[index2D] = tmp;
1089                                else if(tmp < SZ_UINT16_MIN)
1090                                        P0[index2D] = SZ_UINT16_MIN;
1091                                else
1092                                        P0[index2D] = SZ_UINT16_MAX;                                   
1093                        }
1094                        else
1095                        {
1096                                type[index] = 0;
1097
1098                                curValue = P0[index2D] = spaceFillingValue[0];
1099                                compressUInt16Value(curValue, minValue, byteSize, bytes);
1100                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1101                        }
1102
1103                        /* Process Row-0 data 1 --> data r4-1 */
1104                        for (j = 1; j < r4; j++)
1105                        {
1106                                index = l*r234+k*r34+j;
1107                                index2D = j;
1108
1109                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
1110                                diff = spaceFillingValue[index] - pred2D;
1111
1112                                itvNum = llabs(diff)/realPrecision + 1;
1113
1114                                if (itvNum < exe_params->intvCapacity)
1115                                {
1116                                        if (diff < 0) itvNum = -itvNum;
1117                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1118                                        tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1119                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1120                                                P0[index2D] = tmp;
1121                                        else if(tmp < SZ_UINT16_MIN)
1122                                                P0[index2D] = SZ_UINT16_MIN;
1123                                        else
1124                                                P0[index2D] = SZ_UINT16_MAX;                                           
1125                                }
1126                                else
1127                                {
1128                                        type[index] = 0;
1129
1130                                        curValue = P0[index2D] = spaceFillingValue[0];
1131                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
1132                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1133                                }
1134                        }
1135
1136                        /* Process Row-1 --> Row-r3-1 */
1137                        for (i = 1; i < r3; i++)
1138                        {
1139                                /* Process Row-i data 0 */
1140                                index = l*r234+k*r34+i*r4;
1141                                index2D = i*r4;
1142
1143                                pred2D = P0[index2D-r4] + P1[index2D] - P1[index2D-r4];
1144                                diff = spaceFillingValue[index] - pred2D;
1145
1146                                itvNum = llabs(diff)/realPrecision + 1;
1147
1148                                if (itvNum < exe_params->intvCapacity)
1149                                {
1150                                        if (diff < 0) itvNum = -itvNum;
1151                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1152                                        tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1153                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1154                                                P0[index2D] = tmp;
1155                                        else if(tmp < SZ_UINT16_MIN)
1156                                                P0[index2D] = SZ_UINT16_MIN;
1157                                        else
1158                                                P0[index2D] = SZ_UINT16_MAX;                                           
1159                                }
1160                                else
1161                                {
1162                                        type[index] = 0;
1163
1164                                        curValue = P0[index2D] = spaceFillingValue[0];
1165                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
1166                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1167                                }
1168
1169                                /* Process Row-i data 1 --> data r4-1 */
1170                                for (j = 1; j < r4; j++)
1171                                {
1172                                        index = l*r234+k*r34+i*r4+j;
1173                                        index2D = i*r4+j;
1174
1175                                        pred3D = P0[index2D-1] + P0[index2D-r4]+ P1[index2D] - P0[index2D-r4-1] - P1[index2D-r4] - P1[index2D-1] + P1[index2D-r4-1];
1176                                        diff = spaceFillingValue[index] - pred3D;
1177
1178
1179                                        itvNum = llabs(diff)/realPrecision + 1;
1180
1181                                        if (itvNum < exe_params->intvCapacity)
1182                                        {
1183                                                if (diff < 0) itvNum = -itvNum;
1184                                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1185                                                tmp = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1186                                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1187                                                        P0[index2D] = tmp;
1188                                                else if(tmp < SZ_UINT16_MIN)
1189                                                        P0[index2D] = SZ_UINT16_MIN;
1190                                                else
1191                                                        P0[index2D] = SZ_UINT16_MAX;                                                   
1192                                        }
1193                                        else
1194                                        {
1195                                                type[index] = 0;
1196
1197                                                curValue = P0[index2D] = spaceFillingValue[0];
1198                                                compressUInt16Value(curValue, minValue, byteSize, bytes);
1199                                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1200                                        }
1201                                }
1202                        }
1203
1204                        uint16_t *Pt;
1205                        Pt = P1;
1206                        P1 = P0;
1207                        P0 = Pt;
1208                }
1209        }
1210
1211        free(P0);
1212        free(P1);
1213
1214        size_t exactDataNum = exactDataByteArray->size;
1215       
1216        TightDataPointStorageI* tdps;   
1217                       
1218        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
1219                        type, exactDataByteArray->array, exactDataByteArray->size, 
1220                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
1221                       
1222        //free memory
1223        free(type);     
1224        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
1225       
1226        return tdps;   
1227}
1228
1229void SZ_compress_args_uint16_NoCkRngeNoGzip_4D(unsigned char** newByteData, uint16_t *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, 
1230size_t *outSize, int64_t valueRangeSize, int64_t minValue)
1231{
1232        TightDataPointStorageI* tdps = SZ_compress_uint16_4D_MDQ(oriData, r1, r2, r3, r4, realPrecision, valueRangeSize, minValue);
1233
1234        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
1235
1236        size_t dataLength = r1*r2*r3*r4;
1237        if(*outSize>dataLength*sizeof(uint16_t))
1238                SZ_compress_args_uint16_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1239
1240        free_TightDataPointStorageI(tdps);
1241}
1242
1243void SZ_compress_args_uint16_withinRange(unsigned char** newByteData, uint16_t *oriData, size_t dataLength, size_t *outSize)
1244{
1245        TightDataPointStorageI* tdps = (TightDataPointStorageI*) malloc(sizeof(TightDataPointStorageI));
1246        tdps->typeArray = NULL; 
1247       
1248        tdps->allSameData = 1;
1249        tdps->dataSeriesLength = dataLength;
1250        tdps->exactDataBytes = (unsigned char*)malloc(sizeof(unsigned char)*2);
1251        tdps->isLossless = 0;
1252        //tdps->exactByteSize = 4;
1253        tdps->exactDataNum = 1;
1254        tdps->exactDataBytes_size = 2;
1255       
1256        uint16_t value = oriData[0];
1257        int16ToBytes_bigEndian(tdps->exactDataBytes, value);
1258       
1259        size_t tmpOutSize;
1260        convertTDPStoFlatBytes_int(tdps, newByteData, &tmpOutSize);
1261
1262        *outSize = tmpOutSize;//3+1+sizeof(uint16_t)+SZ_SIZE_TYPE; //8==3+1+4(uint16_size)
1263        free_TightDataPointStorageI(tdps);     
1264}
1265
1266int SZ_compress_args_uint16_wRngeNoGzip(unsigned char** newByteData, uint16_t *oriData, 
1267size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1268int errBoundMode, double absErr_Bound, double relBoundRatio)
1269{
1270        int status = SZ_SCES;
1271        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1272        int64_t valueRangeSize = 0;
1273       
1274        uint16_t minValue = computeRangeSize_int(oriData, SZ_UINT16, dataLength, &valueRangeSize);
1275        double realPrecision = getRealPrecision_int(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1276               
1277        if(valueRangeSize <= realPrecision)
1278        {
1279                SZ_compress_args_uint16_withinRange(newByteData, oriData, dataLength, outSize);
1280        }
1281        else
1282        {
1283//              SZ_compress_args_uint16_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize);
1284                if(r5==0&&r4==0&&r3==0&&r2==0)
1285                {
1286                        SZ_compress_args_uint16_NoCkRngeNoGzip_1D(newByteData, oriData, r1, realPrecision, outSize, valueRangeSize, minValue);
1287                }
1288                else if(r5==0&&r4==0&&r3==0)
1289                {
1290                        SZ_compress_args_uint16_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize, valueRangeSize, minValue);
1291                }
1292                else if(r5==0&&r4==0)
1293                {
1294                        SZ_compress_args_uint16_NoCkRngeNoGzip_3D(newByteData, oriData, r3, r2, r1, realPrecision, outSize, valueRangeSize, minValue);
1295                }
1296                else if(r5==0)
1297                {
1298                        SZ_compress_args_uint16_NoCkRngeNoGzip_3D(newByteData, oriData, r4*r3, r2, r1, realPrecision, outSize, valueRangeSize, minValue);
1299                }
1300        }
1301        return status;
1302}
1303
1304int SZ_compress_args_uint16(unsigned char** newByteData, uint16_t *oriData, 
1305size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1306int errBoundMode, double absErr_Bound, double relBoundRatio)
1307{
1308        confparams_cpr->errorBoundMode = errBoundMode;
1309       
1310        if(errBoundMode>=PW_REL)
1311        {
1312                printf("Error: Current SZ version doesn't support integer data compression with point-wise relative error bound being based on pwrType=AVG\n");
1313                exit(0);
1314                return SZ_NSCS;
1315        }
1316        int status = SZ_SCES;
1317        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1318        int64_t valueRangeSize = 0;
1319
1320        uint16_t minValue = (uint16_t)computeRangeSize_int(oriData, SZ_UINT16, dataLength, &valueRangeSize);
1321        double realPrecision = 0; 
1322       
1323        if(confparams_cpr->errorBoundMode==PSNR)
1324        {
1325                confparams_cpr->errorBoundMode = ABS;
1326                realPrecision = confparams_cpr->absErrBound = computeABSErrBoundFromPSNR(confparams_cpr->psnr, (double)confparams_cpr->predThreshold, (double)valueRangeSize);
1327                //printf("realPrecision=%lf\n", realPrecision);
1328        }
1329        else
1330                realPrecision = getRealPrecision_int(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1331
1332        if(valueRangeSize <= realPrecision)
1333        {
1334                SZ_compress_args_uint16_withinRange(newByteData, oriData, dataLength, outSize);
1335        }
1336        else
1337        {
1338                size_t tmpOutSize = 0;
1339                unsigned char* tmpByteData;
1340                if (r2==0)
1341                {
1342                        SZ_compress_args_uint16_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1343                }
1344                else
1345                if (r3==0)
1346                {
1347                        SZ_compress_args_uint16_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1348                }
1349                else
1350                if (r4==0)
1351                {
1352                        SZ_compress_args_uint16_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1353                }
1354                else
1355                if (r5==0)
1356                {
1357                        SZ_compress_args_uint16_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1358                }
1359                else
1360                {
1361                        printf("Error: doesn't support 5 dimensions for now.\n");
1362                        status = SZ_DERR; //dimension error
1363                }
1364                //Call Gzip to do the further compression.
1365                if(confparams_cpr->szMode==SZ_BEST_SPEED)
1366                {
1367                        *outSize = tmpOutSize;
1368                        *newByteData = tmpByteData;
1369                }
1370                else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1371                {
1372                        *outSize = zlib_compress5(tmpByteData, tmpOutSize, newByteData, confparams_cpr->gzipMode);
1373                        free(tmpByteData);
1374                }
1375                else
1376                {
1377                        printf("Error: Wrong setting of confparams_cpr->szMode in the uint16_t compression.\n");
1378                        status = SZ_MERR; //mode error                 
1379                }
1380        }
1381       
1382        return status;
1383}
Note: See TracBrowser for help on using the repository browser.