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

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

more work on adding SZ (latest version)

  • Property mode set to 100644
RevLine 
[2c47b73]1/**
2 *  @file sz_double.c
3 *  @author Sheng Di and Dingwen Tao
4 *  @date Aug, 2016
5 *  @brief SZ_Init, Compression and Decompression functions
6 *  (C) 2016 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 "TightDataPointStorageD.h"
21#include "sz_double.h"
22#include "sz_double_pwr.h"
23#include "szd_double.h"
24#include "szd_double_pwr.h"
25#include "zlib.h"
26#include "rw.h"
27#include "sz_double_ts.h"
28
29unsigned char* SZ_skip_compress_double(double* data, size_t dataLength, size_t* outSize)
30{
31        *outSize = dataLength*sizeof(double);
32        unsigned char* out = (unsigned char*)malloc(dataLength*sizeof(double));
33        memcpy(out, data, dataLength*sizeof(double));
34        return out;
35}
36
37void computeReqLength_double(double realPrecision, short radExpo, int* reqLength, double* medianValue)
38{
39        short reqExpo = getPrecisionReqLength_double(realPrecision);
40        *reqLength = 12+radExpo - reqExpo; //radExpo-reqExpo == reqMantiLength
41        if(*reqLength<12)
42                *reqLength = 12;
43        if(*reqLength>64)
44        {
45                *reqLength = 64;
46                *medianValue = 0;
47        }
48}
49
50unsigned int optimize_intervals_double_1D(double *oriData, size_t dataLength, double realPrecision)
51{       
52        size_t i = 0, radiusIndex;
53        double pred_value = 0, pred_err;
54        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
55        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
56        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
57        for(i=2;i<dataLength;i++)
58        {
59                if(i%confparams_cpr->sampleDistance==0)
60                {
61                        //pred_value = 2*oriData[i-1] - oriData[i-2];
62                        pred_value = oriData[i-1];
63                        pred_err = fabs(pred_value - oriData[i]);
64                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
65                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
66                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
67                        intervals[radiusIndex]++;
68                }
69        }
70        //compute the appropriate number
71        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
72        size_t sum = 0;
73        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
74        {
75                sum += intervals[i];
76                if(sum>targetCount)
77                        break;
78        }
79
80        if(i>=confparams_cpr->maxRangeRadius)
81                i = confparams_cpr->maxRangeRadius-1;
82        unsigned int accIntervals = 2*(i+1);
83        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
84
85        if(powerOf2<32)
86                powerOf2 = 32;
87
88        free(intervals);
89        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
90        return powerOf2;
91}
92
93unsigned int optimize_intervals_double_2D(double *oriData, size_t r1, size_t r2, double realPrecision)
94{       
95        size_t i,j, index;
96        size_t radiusIndex;
97        double pred_value = 0, pred_err;
98        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
99        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
100        size_t totalSampleSize = (r1-1)*(r2-1)/confparams_cpr->sampleDistance;
101        for(i=1;i<r1;i++)
102        {
103                for(j=1;j<r2;j++)
104                {
105                        if((i+j)%confparams_cpr->sampleDistance==0)
106                        {
107                                index = i*r2+j;
108                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
109                                pred_err = fabs(pred_value - oriData[index]);
110                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
111                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
112                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
113                                intervals[radiusIndex]++;
114                        }                       
115                }
116        }
117        //compute the appropriate number
118        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
119        size_t sum = 0;
120        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
121        {
122                sum += intervals[i];
123                if(sum>targetCount)
124                        break;
125        }
126        if(i>=confparams_cpr->maxRangeRadius)
127                i = confparams_cpr->maxRangeRadius-1;   
128        unsigned int accIntervals = 2*(i+1);
129        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
130        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
131
132        if(powerOf2<32)
133                powerOf2 = 32;
134
135        free(intervals);
136        return powerOf2;
137}
138
139unsigned int optimize_intervals_double_3D(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision)
140{       
141        size_t i,j,k, index;
142        size_t radiusIndex;
143        size_t r23=r2*r3;
144        double pred_value = 0, pred_err;
145        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
146        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
147        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)/confparams_cpr->sampleDistance;
148        for(i=1;i<r1;i++)
149        {
150                for(j=1;j<r2;j++)
151                {
152                        for(k=1;k<r3;k++)
153                        {
154                                if((i+j+k)%confparams_cpr->sampleDistance==0)
155                                {
156                                        index = i*r23+j*r3+k;
157                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23] 
158                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
159                                        pred_err = fabs(pred_value - oriData[index]);
160                                        radiusIndex = (pred_err/realPrecision+1)/2;
161                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
162                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
163                                        intervals[radiusIndex]++;
164                                }                               
165                        }
166                       
167                }
168        }
169        //compute the appropriate number
170        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
171        size_t sum = 0;
172        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
173        {
174                sum += intervals[i];
175                if(sum>targetCount)
176                        break;
177        }
178        if(i>=confparams_cpr->maxRangeRadius)
179                i = confparams_cpr->maxRangeRadius-1;
180               
181        unsigned int accIntervals = 2*(i+1);
182        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
183
184        if(powerOf2<32)
185                powerOf2 = 32;
186
187        free(intervals);
188        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
189        return powerOf2;
190}
191
192unsigned int optimize_intervals_double_4D(double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision)
193{
194        size_t i,j,k,l, index;
195        size_t radiusIndex;
196        size_t r234=r2*r3*r4;
197        size_t r34=r3*r4;
198        double pred_value = 0, pred_err;
199        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
200        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
201        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)*(r4-1)/confparams_cpr->sampleDistance;
202        for(i=1;i<r1;i++)
203        {
204                for(j=1;j<r2;j++)
205                {
206                        for(k=1;k<r3;k++)
207                        {
208                                for (l=1;l<r4;l++)
209                                {
210                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
211                                        {
212                                                index = i*r234+j*r34+k*r4+l;
213                                                pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r34]
214                                                                - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
215                                                pred_err = fabs(pred_value - oriData[index]);
216                                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
217                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
218                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
219                                                intervals[radiusIndex]++;
220                                        }
221                                }
222                        }
223                }
224        }
225        //compute the appropriate number
226        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
227        size_t sum = 0;
228        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
229        {
230                sum += intervals[i];
231                if(sum>targetCount)
232                        break;
233        }
234        if(i>=confparams_cpr->maxRangeRadius)
235                i = confparams_cpr->maxRangeRadius-1;
236
237        unsigned int accIntervals = 2*(i+1);
238        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
239
240        if(powerOf2<32)
241                powerOf2 = 32;
242
243        free(intervals);
244        return powerOf2;
245}
246
247TightDataPointStorageD* SZ_compress_double_1D_MDQ(double *oriData, 
248size_t dataLength, double realPrecision, double valueRangeSize, double medianValue_d)
249{
250#ifdef HAVE_TIMECMPR
251        double* decData = NULL; 
252        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
253                decData = (double*)(multisteps->hist_data);
254#endif 
255       
256        unsigned int quantization_intervals;
257        if(exe_params->optQuantMode==1)
258                quantization_intervals = optimize_intervals_double_1D_opt(oriData, dataLength, realPrecision);
259        else
260                quantization_intervals = exe_params->intvCapacity;
261        updateQuantizationInfo(quantization_intervals); 
262
263        size_t i;
264        int reqLength;
265        double medianValue = medianValue_d;
266        short radExpo = getExponent_double(valueRangeSize/2);
267
268        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);     
269
270        int* type = (int*) malloc(dataLength*sizeof(int));
271               
272        double* spaceFillingValue = oriData; //
273       
274        DynamicIntArray *exactLeadNumArray;
275        new_DIA(&exactLeadNumArray, DynArrayInitLen);
276       
277        DynamicByteArray *exactMidByteArray;
278        new_DBA(&exactMidByteArray, DynArrayInitLen);
279       
280        DynamicIntArray *resiBitArray;
281        new_DIA(&resiBitArray, DynArrayInitLen);
282
283        unsigned char preDataBytes[8];
284        longToBytes_bigEndian(preDataBytes, 0);
285       
286        int reqBytesLength = reqLength/8;
287        int resiBitsLength = reqLength%8;
288        double last3CmprsData[3] = {0};
289
290        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
291        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));                       
292                               
293        //add the first data   
294        type[0] = 0;
295        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
296        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
297        memcpy(preDataBytes,vce->curBytes,8);
298        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
299        listAdd_double(last3CmprsData, vce->data);
300#ifdef HAVE_TIMECMPR   
301        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
302                decData[0] = vce->data;
303#endif         
304               
305        //add the second data
306        type[1] = 0;
307        compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
308        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
309        memcpy(preDataBytes,vce->curBytes,8);
310        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
311        listAdd_double(last3CmprsData, vce->data);
312#ifdef HAVE_TIMECMPR   
313        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
314                decData[1] = vce->data;
315#endif
316        int state;
317        double checkRadius;
318        double curData;
319        double pred;
320        double predAbsErr;
321        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
322        double interval = 2*realPrecision;
323
324        for(i=2;i<dataLength;i++)
325        {                               
326                //printf("%.30G\n",last3CmprsData[0]);
327                curData = spaceFillingValue[i];
328                //pred = 2*last3CmprsData[0] - last3CmprsData[1];
329                pred = last3CmprsData[0];
330                predAbsErr = fabs(curData - pred);     
331                if(predAbsErr<=checkRadius)
332                {
333                        state = (predAbsErr/realPrecision+1)/2;
334                        if(curData>=pred)
335                        {
336                                type[i] = exe_params->intvRadius+state;
337                                pred = pred + state*interval;
338                        }
339                        else //curData<pred
340                        {
341                                type[i] = exe_params->intvRadius-state;
342                                pred = pred - state*interval;
343                        }
344                        listAdd_double(last3CmprsData, pred);
345#ifdef HAVE_TIMECMPR                                   
346                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
347                                decData[i] = pred;                     
348#endif 
349                        continue;
350                }
351               
352                //unpredictable data processing
353                type[i] = 0;           
354                compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
355                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
356                memcpy(preDataBytes,vce->curBytes,8);
357                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
358                                                       
359                listAdd_double(last3CmprsData, vce->data);
360#ifdef HAVE_TIMECMPR
361                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
362                        decData[i] = vce->data;
363#endif 
364               
365        }//end of for
366               
367        int exactDataNum = exactLeadNumArray->size;
368       
369        TightDataPointStorageD* tdps;
370                       
371        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum, 
372                        type, exactMidByteArray->array, exactMidByteArray->size, 
373                        exactLeadNumArray->array, 
374                        resiBitArray->array, resiBitArray->size, 
375                        resiBitsLength, 
376                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
377       
378//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
379//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
380       
381        //free memory
382        free_DIA(exactLeadNumArray);
383        free_DIA(resiBitArray);
384        free(type);
385        free(vce);
386        free(lce);     
387        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);     
388       
389        return tdps;   
390}
391
392void SZ_compress_args_double_StoreOriData(double* oriData, size_t dataLength, TightDataPointStorageD* tdps, 
393unsigned char** newByteData, size_t *outSize)
394{
395        int doubleSize = sizeof(double);
396        size_t k = 0, i;
397        tdps->isLossless = 1;
398        size_t totalByteLength = 3 + MetaDataByteLength + exe_params->SZ_SIZE_TYPE + 1 + doubleSize*dataLength;
399        *newByteData = (unsigned char*)malloc(totalByteLength);
400       
401        unsigned char dsLengthBytes[8];
402        for (i = 0; i < 3; i++)//3
403                (*newByteData)[k++] = versionNumber[i];
404       
405        if(exe_params->SZ_SIZE_TYPE==4)//1
406                (*newByteData)[k++] = 16; //00010000
407        else
408                (*newByteData)[k++] = 80;       //01010000: 01000000 indicates the SZ_SIZE_TYPE=8
409
410        convertSZParamsToBytes(confparams_cpr, &((*newByteData)[k]));
411        k = k + MetaDataByteLength;
412
413        sizeToBytes(dsLengthBytes,dataLength);
414        for (i = 0; i < exe_params->SZ_SIZE_TYPE; i++)//ST: 4 or 8
415                (*newByteData)[k++] = dsLengthBytes[i];
416
417        if(sysEndianType==BIG_ENDIAN_SYSTEM)
418                memcpy((*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, oriData, dataLength*doubleSize);
419        else
420        {
421                unsigned char* p = (*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
422                for(i=0;i<dataLength;i++,p+=doubleSize)
423                        doubleToBytes(p, oriData[i]);
424        }
425        *outSize = totalByteLength;
426}
427
428
429char SZ_compress_args_double_NoCkRngeNoGzip_1D(unsigned char** newByteData, double *oriData, 
430size_t dataLength, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
431{
432        char compressionType = 0;       
433        TightDataPointStorageD* tdps = NULL;   
434#ifdef HAVE_TIMECMPR
435        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
436        {
437                int timestep = sz_tsc->currentStep;
438                if(timestep % confparams_cpr->snapshotCmprStep != 0)
439                {
440                        tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
441                        compressionType = 1; //time-series based compression
442                }
443                else
444                {       
445                        tdps = SZ_compress_double_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_d);
446                        compressionType = 0; //snapshot-based compression
447                        multisteps->lastSnapshotStep = timestep;
448                }               
449        }
450        else
451#endif
452                tdps = SZ_compress_double_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_d);                   
453       
454        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
455       
456        if(*outSize>dataLength*sizeof(double))
457                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
458       
459        free_TightDataPointStorageD(tdps);     
460        return compressionType;
461}
462
463TightDataPointStorageD* SZ_compress_double_2D_MDQ(double *oriData, size_t r1, size_t r2, double realPrecision, double valueRangeSize, double medianValue_d)
464{
465#ifdef HAVE_TIMECMPR   
466        double* decData = NULL;
467        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
468                decData = (double*)(multisteps->hist_data);
469#endif 
470       
471        unsigned int quantization_intervals;
472        if(exe_params->optQuantMode==1)
473        {
474                quantization_intervals = optimize_intervals_double_2D_opt(oriData, r1, r2, realPrecision);
475                updateQuantizationInfo(quantization_intervals);
476        }
477        else
478                quantization_intervals = exe_params->intvCapacity;     
479        size_t i,j; 
480        int reqLength;
481        double pred1D, pred2D;
482        double diff = 0.0;
483        double itvNum = 0;
484        double *P0, *P1;
485               
486        size_t dataLength = r1*r2;     
487       
488        P0 = (double*)malloc(r2*sizeof(double));
489        memset(P0, 0, r2*sizeof(double));
490        P1 = (double*)malloc(r2*sizeof(double));
491        memset(P1, 0, r2*sizeof(double));
492               
493        double medianValue = medianValue_d;
494        short radExpo = getExponent_double(valueRangeSize/2);
495        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);     
496
497        int* type = (int*) malloc(dataLength*sizeof(int));
498        //type[dataLength]=0;
499               
500        double* spaceFillingValue = oriData; //
501       
502        DynamicIntArray *exactLeadNumArray;
503        new_DIA(&exactLeadNumArray, DynArrayInitLen);
504       
505        DynamicByteArray *exactMidByteArray;
506        new_DBA(&exactMidByteArray, DynArrayInitLen);
507       
508        DynamicIntArray *resiBitArray;
509        new_DIA(&resiBitArray, DynArrayInitLen);
510       
511        type[0] = 0;
512       
513        unsigned char preDataBytes[8];
514        longToBytes_bigEndian(preDataBytes, 0);
515       
516        int reqBytesLength = reqLength/8;
517        int resiBitsLength = reqLength%8;
518
519        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
520        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
521                       
522        /* Process Row-0 data 0*/
523        type[0] = 0;
524        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
525        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
526        memcpy(preDataBytes,vce->curBytes,8);
527        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
528        P1[0] = vce->data;
529#ifdef HAVE_TIMECMPR   
530        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
531                decData[0] = vce->data;
532#endif 
533
534        /* Process Row-0 data 1*/
535        pred1D = P1[0];
536        diff = spaceFillingValue[1] - pred1D;
537
538        itvNum =  fabs(diff)/realPrecision + 1;
539
540        if (itvNum < exe_params->intvCapacity)
541        {
542                if (diff < 0) itvNum = -itvNum;
543                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
544                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
545        }
546        else
547        {
548                type[1] = 0;
549                compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
550                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
551                memcpy(preDataBytes,vce->curBytes,8);
552                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
553                P1[1] = vce->data;
554        }
555#ifdef HAVE_TIMECMPR   
556        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
557                decData[1] = P1[1];
558#endif
559
560    /* Process Row-0 data 2 --> data r2-1 */
561        for (j = 2; j < r2; j++)
562        {
563                pred1D = 2*P1[j-1] - P1[j-2];
564                diff = spaceFillingValue[j] - pred1D;
565
566                itvNum = fabs(diff)/realPrecision + 1;
567
568                if (itvNum < exe_params->intvCapacity)
569                {
570                        if (diff < 0) itvNum = -itvNum;
571                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
572                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
573                }
574                else
575                {
576                        type[j] = 0;
577                        compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
578                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
579                        memcpy(preDataBytes,vce->curBytes,8);
580                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
581                        P1[j] = vce->data;
582                }
583#ifdef HAVE_TIMECMPR   
584                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
585                        decData[j] = P1[j];
586#endif         
587        }
588
589        /* Process Row-1 --> Row-r1-1 */
590        size_t index;
591        for (i = 1; i < r1; i++)
592        {       
593                /* Process row-i data 0 */
594                index = i*r2;
595                pred1D = P1[0];
596                diff = spaceFillingValue[index] - pred1D;
597
598                itvNum = fabs(diff)/realPrecision + 1;
599
600                if (itvNum < exe_params->intvCapacity)
601                {
602                        if (diff < 0) itvNum = -itvNum;
603                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
604                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
605                }
606                else
607                {
608                        type[index] = 0;
609                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
610                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
611                        memcpy(preDataBytes,vce->curBytes,8);
612                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
613                        P0[0] = vce->data;
614                }
615#ifdef HAVE_TIMECMPR   
616                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
617                        decData[index] = P0[0];
618#endif
619                                                                       
620                /* Process row-i data 1 --> r2-1*/
621                for (j = 1; j < r2; j++)
622                {
623                        index = i*r2+j;
624                        pred2D = P0[j-1] + P1[j] - P1[j-1];
625
626                        diff = spaceFillingValue[index] - pred2D;
627
628                        itvNum = fabs(diff)/realPrecision + 1;
629
630                        if (itvNum < exe_params->intvCapacity)
631                        {
632                                if (diff < 0) itvNum = -itvNum;
633                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
634                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
635                        }
636                        else
637                        {
638                                type[index] = 0;
639                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
640                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
641                                memcpy(preDataBytes,vce->curBytes,8);
642                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
643                                P0[j] = vce->data;
644                        }
645#ifdef HAVE_TIMECMPR   
646                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
647                                decData[index] = P0[j];
648#endif                 
649                }
650
651                double *Pt;
652                Pt = P1;
653                P1 = P0;
654                P0 = Pt;
655        }
656               
657        if(r2!=1)       
658                free(P0);
659        free(P1);
660        size_t exactDataNum = exactLeadNumArray->size;
661       
662        TightDataPointStorageD* tdps;
663                       
664        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum, 
665                        type, exactMidByteArray->array, exactMidByteArray->size, 
666                        exactLeadNumArray->array, 
667                        resiBitArray->array, resiBitArray->size, 
668                        resiBitsLength, 
669                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
670
671/*      int sum =0;
672        for(i=0;i<dataLength;i++)
673                if(type[i]==0) sum++;
674        printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);
675
676        for(i=0;i<dataLength;i++)
677                printf("%d ", type[i]);
678        printf("\n");*/
679
680//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
681//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
682       
683//      for(i = 3800;i<3844;i++)
684//              printf("exactLeadNumArray->array[%d]=%d\n",i,exactLeadNumArray->array[i]);
685       
686        //free memory
687        free_DIA(exactLeadNumArray);
688        free_DIA(resiBitArray);
689        free(type);     
690        free(vce);
691        free(lce);     
692        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
693       
694        return tdps;
695}
696
697/**
698 *
699 * Note: @r1 is high dimension
700 *               @r2 is low dimension
701 * */
702char SZ_compress_args_double_NoCkRngeNoGzip_2D(unsigned char** newByteData, double *oriData, size_t r1, size_t r2, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
703{
704        size_t dataLength = r1*r2;
705        char compressionType = 0;       
706        TightDataPointStorageD* tdps = NULL;   
707#ifdef HAVE_TIMECMPR
708        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
709        {
710                int timestep = sz_tsc->currentStep;
711                if(timestep % confparams_cpr->snapshotCmprStep != 0)
712                {
713                        tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
714                        compressionType = 1; //time-series based compression
715                }
716                else
717                {       
718                        tdps = SZ_compress_double_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_d);
719                        compressionType = 0; //snapshot-based compression
720                        multisteps->lastSnapshotStep = timestep;
721                }               
722        }
723        else
724#endif
725                tdps = SZ_compress_double_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_d);       
726       
727        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
728       
729        if(*outSize>dataLength*sizeof(double))
730                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize); 
731       
732        free_TightDataPointStorageD(tdps);
733        return compressionType;
734}
735
736TightDataPointStorageD* SZ_compress_double_3D_MDQ(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, double valueRangeSize, double medianValue_d)
737{
738#ifdef HAVE_TIMECMPR
739        double* decData = NULL;
740        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
741                decData = (double*)(multisteps->hist_data);
742#endif         
743
744        unsigned int quantization_intervals;
745        if(exe_params->optQuantMode==1)
746        {
747                quantization_intervals = optimize_intervals_double_3D_opt(oriData, r1, r2, r3, realPrecision);
748                updateQuantizationInfo(quantization_intervals);
749        }       
750        else
751                quantization_intervals = exe_params->intvCapacity;
752        size_t i,j,k; 
753        int reqLength;
754        double pred1D, pred2D, pred3D;
755        double diff = 0.0;
756        double itvNum = 0;
757        double *P0, *P1;
758
759        size_t dataLength = r1*r2*r3;
760
761        size_t r23 = r2*r3;
762
763        P0 = (double*)malloc(r23*sizeof(double));
764        P1 = (double*)malloc(r23*sizeof(double));
765
766        double medianValue = medianValue_d;
767        short radExpo = getExponent_double(valueRangeSize/2);
768        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);     
769
770        int* type = (int*) malloc(dataLength*sizeof(int));
771        //type[dataLength]=0;
772
773        double* spaceFillingValue = oriData; //
774
775        DynamicIntArray *exactLeadNumArray;
776        new_DIA(&exactLeadNumArray, DynArrayInitLen);
777
778        DynamicByteArray *exactMidByteArray;
779        new_DBA(&exactMidByteArray, DynArrayInitLen);
780
781        DynamicIntArray *resiBitArray;
782        new_DIA(&resiBitArray, DynArrayInitLen);
783
784        type[0] = 0;
785
786        unsigned char preDataBytes[8];
787        longToBytes_bigEndian(preDataBytes, 0);
788
789        int reqBytesLength = reqLength/8;
790        int resiBitsLength = reqLength%8;
791
792        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
793        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
794
795
796        ///////////////////////////     Process layer-0 ///////////////////////////
797        /* Process Row-0 data 0*/
798        type[0] = 0;
799        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
800        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
801        memcpy(preDataBytes,vce->curBytes,8);
802        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
803        P1[0] = vce->data;
804#ifdef HAVE_TIMECMPR   
805                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
806                        decData[0] = P1[0];
807#endif
808
809        /* Process Row-0 data 1*/
810        pred1D = P1[0];
811        diff = spaceFillingValue[1] - pred1D;
812
813        itvNum = fabs(diff)/realPrecision + 1;
814
815        if (itvNum < exe_params->intvCapacity)
816        {
817                if (diff < 0) itvNum = -itvNum;
818                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
819                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
820        }
821        else
822        {
823                type[1] = 0;
824                compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
825                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
826                memcpy(preDataBytes,vce->curBytes,8);
827                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
828                P1[1] = vce->data;
829        }
830#ifdef HAVE_TIMECMPR   
831        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
832                decData[1] = P1[1];
833#endif
834
835    /* Process Row-0 data 2 --> data r3-1 */
836        for (j = 2; j < r3; j++)
837        {
838                pred1D = 2*P1[j-1] - P1[j-2];
839                diff = spaceFillingValue[j] - pred1D;
840
841                itvNum = fabs(diff)/realPrecision + 1;
842
843                if (itvNum < exe_params->intvCapacity)
844                {
845                        if (diff < 0) itvNum = -itvNum;
846                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
847                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
848                }
849                else
850                {
851                        type[j] = 0;
852                        compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
853                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
854                        memcpy(preDataBytes,vce->curBytes,8);
855                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
856                        P1[j] = vce->data;
857                }
858#ifdef HAVE_TIMECMPR   
859                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
860                        decData[j] = P1[j];
861#endif         
862        }
863
864        /* Process Row-1 --> Row-r2-1 */
865        size_t index;
866        for (i = 1; i < r2; i++)
867        {
868                /* Process row-i data 0 */
869                index = i*r3;
870                pred1D = P1[index-r3];
871                diff = spaceFillingValue[index] - pred1D;
872
873                itvNum = fabs(diff)/realPrecision + 1;
874
875                if (itvNum < exe_params->intvCapacity)
876                {
877                        if (diff < 0) itvNum = -itvNum;
878                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
879                        P1[index] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
880                }
881                else
882                {
883                        type[index] = 0;
884                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
885                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
886                        memcpy(preDataBytes,vce->curBytes,8);
887                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
888                        P1[index] = vce->data;
889                }
890#ifdef HAVE_TIMECMPR   
891                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
892                        decData[index] = P1[index];
893#endif         
894
895                /* Process row-i data 1 --> data r3-1*/
896                for (j = 1; j < r3; j++)
897                {
898                        index = i*r3+j;
899                        pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];
900
901                        diff = spaceFillingValue[index] - pred2D;
902
903                        itvNum = fabs(diff)/realPrecision + 1;
904
905                        if (itvNum < exe_params->intvCapacity)
906                        {
907                                if (diff < 0) itvNum = -itvNum;
908                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
909                                P1[index] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
910                        }
911                        else
912                        {
913                                type[index] = 0;
914                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
915                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
916                                memcpy(preDataBytes,vce->curBytes,8);
917                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
918                                P1[index] = vce->data;
919                        }
920#ifdef HAVE_TIMECMPR   
921                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
922                                decData[index] = P1[index];
923#endif                 
924                }
925        }
926
927
928        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
929
930        for (k = 1; k < r1; k++)
931        {
932                /* Process Row-0 data 0*/
933                index = k*r23;
934                pred1D = P1[0];
935                diff = spaceFillingValue[index] - pred1D;
936
937                itvNum = fabs(diff)/realPrecision + 1;
938
939                if (itvNum < exe_params->intvCapacity)
940                {
941                        if (diff < 0) itvNum = -itvNum;
942                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
943                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
944                }
945                else
946                {
947                        type[index] = 0;
948                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
949                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
950                        memcpy(preDataBytes,vce->curBytes,8);
951                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
952                        P0[0] = vce->data;
953                }
954#ifdef HAVE_TIMECMPR   
955                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
956                        decData[index] = P0[0];
957#endif
958
959            /* Process Row-0 data 1 --> data r3-1 */
960                for (j = 1; j < r3; j++)
961                {
962                        //index = k*r2*r3+j;
963                        index ++;
964                        pred2D = P0[j-1] + P1[j] - P1[j-1];
965                        diff = spaceFillingValue[index] - pred2D;
966
967                        itvNum = fabs(diff)/realPrecision + 1;
968
969                        if (itvNum < exe_params->intvCapacity)
970                        {
971                                if (diff < 0) itvNum = -itvNum;
972                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
973                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
974                        }
975                        else
976                        {
977                                type[index] = 0;
978                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
979                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
980                                memcpy(preDataBytes,vce->curBytes,8);
981                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
982                                P0[j] = vce->data;
983                        }
984#ifdef HAVE_TIMECMPR   
985                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
986                                decData[index] = P0[j];
987#endif                 
988                }
989
990            /* Process Row-1 --> Row-r2-1 */
991                size_t index2D;
992                for (i = 1; i < r2; i++)
993                {
994                        /* Process Row-i data 0 */
995                        index = k*r23 + i*r3;
996                        index2D = i*r3;
997                        pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
998                        diff = spaceFillingValue[index] - pred2D;
999
1000                        itvNum = fabs(diff)/realPrecision + 1;
1001
1002                        if (itvNum < exe_params->intvCapacity)
1003                        {
1004                                if (diff < 0) itvNum = -itvNum;
1005                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1006                                P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1007                        }
1008                        else
1009                        {
1010                                type[index] = 0;
1011                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1012                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1013                                memcpy(preDataBytes,vce->curBytes,8);
1014                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1015                                P0[index2D] = vce->data;
1016                        }
1017#ifdef HAVE_TIMECMPR   
1018                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1019                                decData[index] = P0[index2D];
1020#endif                 
1021
1022                        /* Process Row-i data 1 --> data r3-1 */
1023                        for (j = 1; j < r3; j++)
1024                        {
1025                                //index = k*r2*r3 + i*r3 + j;
1026                                index ++;
1027                                index2D = i*r3 + j;
1028                                pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
1029                                diff = spaceFillingValue[index] - pred3D;
1030
1031                                itvNum = fabs(diff)/realPrecision + 1;
1032
1033                                if (itvNum < exe_params->intvCapacity)
1034                                {
1035                                        if (diff < 0) itvNum = -itvNum;
1036                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1037                                        P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1038                                }
1039                                else
1040                                {
1041                                        type[index] = 0;
1042                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1043                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1044                                        memcpy(preDataBytes,vce->curBytes,8);
1045                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1046                                        P0[index2D] = vce->data;
1047                                }
1048#ifdef HAVE_TIMECMPR   
1049                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1050                                        decData[index] = P0[index2D];
1051#endif                         
1052                        }
1053                }
1054
1055                double *Pt;
1056                Pt = P1;
1057                P1 = P0;
1058                P0 = Pt;
1059        }
1060        if(r23!=1)
1061                free(P0);
1062        free(P1);
1063        size_t exactDataNum = exactLeadNumArray->size;
1064
1065        TightDataPointStorageD* tdps;
1066
1067        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
1068                        type, exactMidByteArray->array, exactMidByteArray->size,
1069                        exactLeadNumArray->array,
1070                        resiBitArray->array, resiBitArray->size,
1071                        resiBitsLength, 
1072                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
1073
1074//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
1075//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
1076
1077//      for(i = 3800;i<3844;i++)
1078//              printf("exactLeadNumArray->array[%d]=%d\n",i,exactLeadNumArray->array[i]);
1079
1080        //free memory
1081        free_DIA(exactLeadNumArray);
1082        free_DIA(resiBitArray);
1083        free(type);
1084        free(vce);
1085        free(lce);     
1086        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);     
1087       
1088        return tdps;   
1089}
1090
1091
1092char SZ_compress_args_double_NoCkRngeNoGzip_3D(unsigned char** newByteData, double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
1093{
1094        size_t dataLength = r1*r2*r3;
1095        char compressionType = 0;       
1096        TightDataPointStorageD* tdps = NULL;   
1097#ifdef HAVE_TIMECMPR
1098        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1099        {
1100                int timestep = sz_tsc->currentStep;
1101                if(timestep % confparams_cpr->snapshotCmprStep != 0)
1102                {
1103                        tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
1104                        compressionType = 1; //time-series based compression
1105                }
1106                else
1107                {       
1108                        tdps = SZ_compress_double_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_d);
1109                        compressionType = 0; //snapshot-based compression
1110                        multisteps->lastSnapshotStep = timestep;
1111                }               
1112        }
1113        else
1114#endif
1115                tdps = SZ_compress_double_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_d);           
1116
1117        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
1118
1119        if(*outSize>dataLength*sizeof(double))
1120                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1121
1122        free_TightDataPointStorageD(tdps);
1123        return compressionType;
1124}
1125
1126TightDataPointStorageD* SZ_compress_double_4D_MDQ(double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, double valueRangeSize, double medianValue_d)
1127{
1128        unsigned int quantization_intervals;
1129        if(exe_params->optQuantMode==1)
1130        {
1131                quantization_intervals = optimize_intervals_double_4D(oriData, r1, r2, r3, r4, realPrecision);
1132                updateQuantizationInfo(quantization_intervals);
1133        }
1134        else
1135                quantization_intervals = exe_params->intvCapacity;
1136
1137        size_t i,j,k; 
1138        int reqLength;
1139        double pred1D, pred2D, pred3D;
1140        double diff = 0.0;
1141        double itvNum = 0;
1142        double *P0, *P1;
1143
1144        size_t dataLength = r1*r2*r3*r4;
1145
1146        size_t r234 = r2*r3*r4;
1147        size_t r34 = r3*r4;
1148
1149        P0 = (double*)malloc(r34*sizeof(double));
1150        P1 = (double*)malloc(r34*sizeof(double));
1151
1152        double medianValue = medianValue_d;
1153        short radExpo = getExponent_double(valueRangeSize/2);
1154        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1155
1156        int* type = (int*) malloc(dataLength*sizeof(int));
1157
1158        double* spaceFillingValue = oriData; //
1159
1160        DynamicIntArray *exactLeadNumArray;
1161        new_DIA(&exactLeadNumArray, DynArrayInitLen);
1162
1163        DynamicByteArray *exactMidByteArray;
1164        new_DBA(&exactMidByteArray, DynArrayInitLen);
1165
1166        DynamicIntArray *resiBitArray;
1167        new_DIA(&resiBitArray, DynArrayInitLen);
1168
1169        unsigned char preDataBytes[8];
1170        longToBytes_bigEndian(preDataBytes, 0);
1171
1172        int reqBytesLength = reqLength/8;
1173        int resiBitsLength = reqLength%8;
1174
1175        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
1176        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
1177
1178
1179        size_t l;
1180        for (l = 0; l < r1; l++)
1181        {
1182
1183                ///////////////////////////     Process layer-0 ///////////////////////////
1184                /* Process Row-0 data 0*/
1185                size_t index = l*r234;
1186                size_t index2D = 0;
1187
1188                type[index] = 0;
1189                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1190                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1191                memcpy(preDataBytes,vce->curBytes,8);
1192                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1193                P1[index2D] = vce->data;
1194
1195                /* Process Row-0 data 1*/
1196                index = l*r234+1;
1197                index2D = 1;
1198
1199                pred1D = P1[index2D-1];
1200                diff = spaceFillingValue[index] - pred1D;
1201
1202                itvNum = fabs(diff)/realPrecision + 1;
1203
1204                if (itvNum < exe_params->intvCapacity)
1205                {
1206                        if (diff < 0) itvNum = -itvNum;
1207                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1208                        P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1209                }
1210                else
1211                {
1212                        type[index] = 0;
1213                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1214                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1215                        memcpy(preDataBytes,vce->curBytes,8);
1216                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1217                        P1[index2D] = vce->data;
1218                }
1219
1220                /* Process Row-0 data 2 --> data r4-1 */
1221                for (j = 2; j < r4; j++)
1222                {
1223                        index = l*r234+j;
1224                        index2D = j;
1225
1226                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
1227                        diff = spaceFillingValue[index] - pred1D;
1228
1229                        itvNum = fabs(diff)/realPrecision + 1;
1230
1231                        if (itvNum < exe_params->intvCapacity)
1232                        {
1233                                if (diff < 0) itvNum = -itvNum;
1234                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1235                                P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1236                        }
1237                        else
1238                        {
1239                                type[index] = 0;
1240                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1241                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1242                                memcpy(preDataBytes,vce->curBytes,8);
1243                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1244                                P1[index2D] = vce->data;
1245                        }
1246                }
1247
1248                /* Process Row-1 --> Row-r3-1 */
1249                for (i = 1; i < r3; i++)
1250                {
1251                        /* Process row-i data 0 */
1252                        index = l*r234+i*r4;
1253                        index2D = i*r4;
1254
1255                        pred1D = P1[index2D-r4];
1256                        diff = spaceFillingValue[index] - pred1D;
1257
1258                        itvNum = fabs(diff)/realPrecision + 1;
1259
1260                        if (itvNum < exe_params->intvCapacity)
1261                        {
1262                                if (diff < 0) itvNum = -itvNum;
1263                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1264                                P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1265                        }
1266                        else
1267                        {
1268                                type[index] = 0;
1269                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1270                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1271                                memcpy(preDataBytes,vce->curBytes,8);
1272                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1273                                P1[index2D] = vce->data;
1274                        }
1275
1276                        /* Process row-i data 1 --> data r4-1*/
1277                        for (j = 1; j < r4; j++)
1278                        {
1279                                index = l*r234+i*r4+j;
1280                                index2D = i*r4+j;
1281
1282                                pred2D = P1[index2D-1] + P1[index2D-r4] - P1[index2D-r4-1];
1283
1284                                diff = spaceFillingValue[index] - pred2D;
1285
1286                                itvNum = fabs(diff)/realPrecision + 1;
1287
1288                                if (itvNum < exe_params->intvCapacity)
1289                                {
1290                                        if (diff < 0) itvNum = -itvNum;
1291                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1292                                        P1[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1293                                }
1294                                else
1295                                {
1296                                        type[index] = 0;
1297                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1298                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1299                                        memcpy(preDataBytes,vce->curBytes,8);
1300                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1301                                        P1[index2D] = vce->data;
1302                                }
1303                        }
1304                }
1305
1306
1307                ///////////////////////////     Process layer-1 --> layer-r2-1 ///////////////////////////
1308
1309                for (k = 1; k < r2; k++)
1310                {
1311                        /* Process Row-0 data 0*/
1312                        index = l*r234+k*r34;
1313                        index2D = 0;
1314
1315                        pred1D = P1[index2D];
1316                        diff = spaceFillingValue[index] - pred1D;
1317
1318                        itvNum = fabs(diff)/realPrecision + 1;
1319
1320                        if (itvNum < exe_params->intvCapacity)
1321                        {
1322                                if (diff < 0) itvNum = -itvNum;
1323                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1324                                P0[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1325                        }
1326                        else
1327                        {
1328                                type[index] = 0;
1329                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1330                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1331                                memcpy(preDataBytes,vce->curBytes,8);
1332                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1333                                P0[index2D] = vce->data;
1334                        }
1335
1336
1337                        /* Process Row-0 data 1 --> data r4-1 */
1338                        for (j = 1; j < r4; j++)
1339                        {
1340                                index = l*r234+k*r34+j;
1341                                index2D = j;
1342
1343                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
1344                                diff = spaceFillingValue[index] - pred2D;
1345
1346                                itvNum = fabs(diff)/realPrecision + 1;
1347
1348                                if (itvNum < exe_params->intvCapacity)
1349                                {
1350                                        if (diff < 0) itvNum = -itvNum;
1351                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1352                                        P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1353                                }
1354                                else
1355                                {
1356                                        type[index] = 0;
1357                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1358                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1359                                        memcpy(preDataBytes,vce->curBytes,8);
1360                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1361                                        P0[index2D] = vce->data;
1362                                }
1363                        }
1364
1365                        /* Process Row-1 --> Row-r3-1 */
1366                        for (i = 1; i < r3; i++)
1367                        {
1368                                /* Process Row-i data 0 */
1369                                index = l*r234+k*r34+i*r4;
1370                                index2D = i*r4;
1371
1372                                pred2D = P0[index2D-r4] + P1[index2D] - P1[index2D-r4];
1373                                diff = spaceFillingValue[index] - pred2D;
1374
1375                                itvNum = fabs(diff)/realPrecision + 1;
1376
1377                                if (itvNum < exe_params->intvCapacity)
1378                                {
1379                                        if (diff < 0) itvNum = -itvNum;
1380                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1381                                        P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1382                                }
1383                                else
1384                                {
1385                                        type[index] = 0;
1386                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1387                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1388                                        memcpy(preDataBytes,vce->curBytes,8);
1389                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1390                                        P0[index2D] = vce->data;
1391                                }
1392
1393                                /* Process Row-i data 1 --> data r4-1 */
1394                                for (j = 1; j < r4; j++)
1395                                {
1396                                        index = l*r234+k*r34+i*r4+j;
1397                                        index2D = i*r4+j;
1398
1399                                        pred3D = P0[index2D-1] + P0[index2D-r4]+ P1[index2D] - P0[index2D-r4-1] - P1[index2D-r4] - P1[index2D-1] + P1[index2D-r4-1];
1400                                        diff = spaceFillingValue[index] - pred3D;
1401
1402
1403                                        itvNum = fabs(diff)/realPrecision + 1;
1404
1405                                        if (itvNum < exe_params->intvCapacity)
1406                                        {
1407                                                if (diff < 0) itvNum = -itvNum;
1408                                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1409                                                P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1410                                        }
1411                                        else
1412                                        {
1413                                                type[index] = 0;
1414                                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1415                                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1416                                                memcpy(preDataBytes,vce->curBytes,8);
1417                                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1418                                                P0[index2D] = vce->data;
1419                                        }
1420                                }
1421                        }
1422
1423                        double *Pt;
1424                        Pt = P1;
1425                        P1 = P0;
1426                        P0 = Pt;
1427                }
1428        }
1429
1430        free(P0);
1431        free(P1);
1432        size_t exactDataNum = exactLeadNumArray->size;
1433
1434        TightDataPointStorageD* tdps;
1435
1436        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
1437                        type, exactMidByteArray->array, exactMidByteArray->size,
1438                        exactLeadNumArray->array,
1439                        resiBitArray->array, resiBitArray->size,
1440                        resiBitsLength,
1441                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
1442
1443        //free memory
1444        free_DIA(exactLeadNumArray);
1445        free_DIA(resiBitArray);
1446        free(type);
1447        free(vce);
1448        free(lce);
1449        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
1450
1451        return tdps;
1452}
1453
1454
1455char SZ_compress_args_double_NoCkRngeNoGzip_4D(unsigned char** newByteData, double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
1456{
1457        TightDataPointStorageD* tdps = SZ_compress_double_4D_MDQ(oriData, r1, r2, r3, r4, realPrecision, valueRangeSize, medianValue_d);
1458
1459        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
1460
1461        size_t dataLength = r1*r2*r3*r4;
1462        if(*outSize>dataLength*sizeof(double))
1463                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1464
1465        free_TightDataPointStorageD(tdps);
1466        return 0;
1467}
1468
1469void SZ_compress_args_double_withinRange(unsigned char** newByteData, double *oriData, size_t dataLength, size_t *outSize)
1470{
1471        TightDataPointStorageD* tdps = (TightDataPointStorageD*) malloc(sizeof(TightDataPointStorageD));
1472        tdps->rtypeArray = NULL;
1473        tdps->typeArray = NULL;
1474        tdps->leadNumArray = NULL;
1475        tdps->residualMidBits = NULL;
1476       
1477        tdps->allSameData = 1;
1478        tdps->dataSeriesLength = dataLength;
1479        tdps->exactMidBytes = (unsigned char*)malloc(sizeof(unsigned char)*8);
1480        tdps->pwrErrBoundBytes = NULL;
1481        tdps->isLossless = 0;
1482        double value = oriData[0];
1483        doubleToBytes(tdps->exactMidBytes, value);
1484        tdps->exactMidBytes_size = 8;
1485       
1486        size_t tmpOutSize;
1487        //unsigned char *tmpByteData;
1488        convertTDPStoFlatBytes_double(tdps, newByteData, &tmpOutSize);
1489        //convertTDPStoFlatBytes_double(tdps, &tmpByteData, &tmpOutSize);
1490
1491        //*newByteData = (unsigned char*)malloc(sizeof(unsigned char)*16); //for floating-point data (1+3+4+4)
1492        //memcpy(*newByteData, tmpByteData, 16);
1493        *outSize = tmpOutSize;//12==3+1+8(double_size)+MetaDataByteLength
1494        free_TightDataPointStorageD(tdps);     
1495}
1496
1497int SZ_compress_args_double_wRngeNoGzip(unsigned char** newByteData, double *oriData, 
1498size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1499int errBoundMode, double absErr_Bound, double relBoundRatio, double pwrErrRatio)
1500{
1501        int status = SZ_SCES;
1502        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1503        double valueRangeSize = 0, medianValue = 0;
1504       
1505        double min = computeRangeSize_double(oriData, dataLength, &valueRangeSize, &medianValue);
1506        double max = min+valueRangeSize;
1507        double realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1508               
1509        if(valueRangeSize <= realPrecision)
1510        {
1511                SZ_compress_args_double_withinRange(newByteData, oriData, dataLength, outSize);
1512        }
1513        else
1514        {
1515                if(r5==0&&r4==0&&r3==0&&r2==0)
1516                {
1517                        if(errBoundMode>=PW_REL)
1518                        {
1519                                //SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr(newByteData, oriData, realPrecision, r1, outSize, min, max);
1520                                SZ_compress_args_double_NoCkRngeNoGzip_1D_pwrgroup(newByteData, oriData, r1, absErr_Bound, relBoundRatio, pwrErrRatio, valueRangeSize, medianValue, outSize);                           
1521                        }
1522                        else
1523                                SZ_compress_args_double_NoCkRngeNoGzip_1D(newByteData, oriData, r1, realPrecision, outSize, valueRangeSize, medianValue);
1524                }
1525                else if(r5==0&&r4==0&&r3==0)
1526                {
1527                        if(errBoundMode>=PW_REL)
1528                                SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr(newByteData, oriData, realPrecision, r2, r1, outSize, min, max);
1529                        else
1530                                SZ_compress_args_double_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1531                }
1532                else if(r5==0&&r4==0)
1533                {
1534                        if(errBoundMode>=PW_REL)
1535                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r3, r2, r1, outSize, min, max);
1536                        else
1537                                SZ_compress_args_double_NoCkRngeNoGzip_3D(newByteData, oriData, r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1538                }
1539                else if(r5==0)
1540                {
1541                        if(errBoundMode>=PW_REL)
1542                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r4*r3, r2, r1, outSize, min, max);
1543                        else
1544                                SZ_compress_args_double_NoCkRngeNoGzip_3D(newByteData, oriData, r4*r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1545                }
1546        }
1547        return status;
1548}
1549
1550int SZ_compress_args_double(unsigned char** newByteData, double *oriData, 
1551size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1552int errBoundMode, double absErr_Bound, double relBoundRatio, double pwRelBoundRatio)
1553{
1554        confparams_cpr->errorBoundMode = errBoundMode;
1555        if(errBoundMode==PW_REL)
1556        {
1557                confparams_cpr->pw_relBoundRatio = pwRelBoundRatio;     
1558                //confparams_cpr->pwr_type = SZ_PWR_MIN_TYPE;
1559                if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE && r3 != 0 )
1560                {
1561                        printf("Error: Current version doesn't support 3D data compression with point-wise relative error bound being based on pwrType=AVG\n");
1562                        exit(0);
1563                        return SZ_NSCS;
1564                }
1565        }                               
1566               
1567        int status = SZ_SCES;
1568        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1569       
1570        if(dataLength <= MIN_NUM_OF_ELEMENTS)
1571        {
1572                *newByteData = SZ_skip_compress_double(oriData, dataLength, outSize);
1573                return status;
1574        }
1575       
1576        double valueRangeSize = 0, medianValue = 0;
1577       
1578        double min = computeRangeSize_double(oriData, dataLength, &valueRangeSize, &medianValue);
1579        double max = min+valueRangeSize;
1580
1581        double realPrecision = 0; 
1582       
1583        if(confparams_cpr->errorBoundMode==PSNR)
1584        {
1585                confparams_cpr->errorBoundMode = ABS;
1586                realPrecision = confparams_cpr->absErrBound = computeABSErrBoundFromPSNR(confparams_cpr->psnr, (double)confparams_cpr->predThreshold, valueRangeSize);
1587        }
1588        else
1589                realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1590               
1591        if(valueRangeSize <= realPrecision)
1592        {
1593                SZ_compress_args_double_withinRange(newByteData, oriData, dataLength, outSize);
1594        }
1595        else
1596        {
1597                size_t tmpOutSize = 0;
1598                unsigned char* tmpByteData;
1599                if (r2==0)
1600                {
1601                        if(confparams_cpr->errorBoundMode>=PW_REL)
1602                        {
1603                                //SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr(&tmpByteData, oriData, realPrecision, r1, &tmpOutSize, min, max);
1604                                SZ_compress_args_double_NoCkRngeNoGzip_1D_pwrgroup(&tmpByteData, oriData, r1, absErr_Bound, relBoundRatio, pwRelBoundRatio, 
1605                                valueRangeSize, medianValue, &tmpOutSize);
1606                        }
1607                        else
1608#ifdef HAVE_TIMECMPR
1609                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                   
1610                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1611                                else
1612#endif
1613                                        SZ_compress_args_double_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1614                }
1615                else
1616                if (r3==0)
1617                {
1618                        if(confparams_cpr->errorBoundMode>=PW_REL)
1619                                SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr(&tmpByteData, oriData, realPrecision, r2, r1, &tmpOutSize, min, max);
1620                        else
1621#ifdef HAVE_TIMECMPR
1622                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                   
1623                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1624                                else
1625#endif
1626                                        SZ_compress_args_double_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1627                }
1628                else
1629                if (r4==0)
1630                {
1631                        if(confparams_cpr->errorBoundMode>=PW_REL)
1632                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r3, r2, r1, &tmpOutSize, min, max);
1633                        else
1634#ifdef HAVE_TIMECMPR
1635                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1636                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1637                                else
1638#endif
1639                                        SZ_compress_args_double_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1640                }
1641                else
1642                if (r5==0)
1643                {
1644                        if(confparams_cpr->errorBoundMode>=PW_REL)
1645                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r4*r3, r2, r1, &tmpOutSize, min, max);
1646                        else
1647#ifdef HAVE_TIMECMPR
1648                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                   
1649                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1650                                else
1651#endif
1652                                        SZ_compress_args_double_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1653                }
1654                else
1655                {
1656                        printf("Error: doesn't support 5 dimensions for now.\n");
1657                        status = SZ_DERR;
1658                }
1659                               
1660                //Call Gzip to do the further compression.
1661                if(confparams_cpr->szMode==SZ_BEST_SPEED)
1662                {
1663                        *outSize = tmpOutSize;
1664                        *newByteData = tmpByteData;                     
1665                }
1666                else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1667                {
1668                        *outSize = zlib_compress5(tmpByteData, tmpOutSize, newByteData, confparams_cpr->gzipMode);
1669                        free(tmpByteData);
1670                }
1671                else
1672                {
1673                        printf("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1674                        status = SZ_MERR;       
1675                }
1676        }
1677
1678        return status;
1679}
1680
1681//TODO
1682int SZ_compress_args_double_subblock(unsigned char* compressedBytes, double *oriData,
1683size_t r5, size_t r4, size_t r3, size_t r2, size_t r1,
1684size_t s5, size_t s4, size_t s3, size_t s2, size_t s1,
1685size_t e5, size_t e4, size_t e3, size_t e2, size_t e1,
1686size_t *outSize, int errBoundMode, double absErr_Bound, double relBoundRatio)
1687{
1688        int status = SZ_SCES;
1689        double valueRangeSize = 0, medianValue = 0;
1690        computeRangeSize_double_subblock(oriData, &valueRangeSize, &medianValue, r5, r4, r3, r2, r1, s5, s4, s3, s2, s1, e5, e4, e3, e2, e1);
1691
1692        double realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1693
1694        if(valueRangeSize <= realPrecision)
1695        {
1696                //TODO
1697                //SZ_compress_args_double_withinRange_subblock();
1698        }
1699        else
1700        {
1701                if (r2==0)
1702                {
1703                        //TODO
1704                        if(errBoundMode==PW_REL)
1705                        {
1706                                //TODO
1707                                //SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr_subblock();
1708                                printf ("Current subblock version does not support point-wise relative error bound.\n");
1709                        }
1710                        else
1711                                SZ_compress_args_double_NoCkRnge_1D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r1, s1, e1);
1712                }
1713                else
1714                if (r3==0)
1715                {
1716                        if(errBoundMode==PW_REL)
1717                        {
1718                                //TODO
1719                                //SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr_subblock();
1720                                printf ("Current subblock version does not support point-wise relative error bound.\n");
1721                        }
1722                        else
1723                                SZ_compress_args_double_NoCkRnge_2D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r2, r1, s2, s1, e2, e1);
1724                }
1725                else
1726                if (r4==0)
1727                {
1728                        if(errBoundMode==PW_REL)
1729                        {
1730                                //TODO
1731                                //SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr_subblock();
1732                                printf ("Current subblock version does not support point-wise relative error bound.\n");
1733                        }
1734                        else
1735                                SZ_compress_args_double_NoCkRnge_3D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r3, r2, r1, s3, s2, s1, e3, e2, e1);
1736                }
1737                else
1738                if (r5==0)
1739                {
1740                        if(errBoundMode==PW_REL)
1741                        {
1742                                //TODO
1743                                //SZ_compress_args_double_NoCkRngeNoGzip_4D_pwr_subblock();
1744                                printf ("Current subblock version does not support point-wise relative error bound.\n");
1745                        }
1746                        else
1747                                SZ_compress_args_double_NoCkRnge_4D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);
1748                }
1749                else
1750                {
1751                        printf("Error: doesn't support 5 dimensions for now.\n");
1752                        status = SZ_DERR; //dimension error
1753                }
1754        }
1755        return status;
1756}
1757
1758void SZ_compress_args_double_NoCkRnge_1D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
1759size_t r1, size_t s1, size_t e1)
1760{
1761        TightDataPointStorageD* tdps = SZ_compress_double_1D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r1, s1, e1);
1762
1763        if (confparams_cpr->szMode==SZ_BEST_SPEED)
1764                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1765        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1766        {
1767                unsigned char *tmpCompBytes;
1768                size_t tmpOutSize;
1769                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1770                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1771                free(tmpCompBytes);
1772        }
1773        else
1774        {
1775                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1776        }
1777
1778        //TODO
1779//      if(*outSize>dataLength*sizeof(double))
1780//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1781
1782        free_TightDataPointStorageD(tdps);
1783}
1784
1785void SZ_compress_args_double_NoCkRnge_2D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
1786size_t r2, size_t r1, size_t s2, size_t s1, size_t e2, size_t e1)
1787{
1788        TightDataPointStorageD* tdps = SZ_compress_double_2D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r2, r1, s2, s1, e2, e1);
1789
1790        if (confparams_cpr->szMode==SZ_BEST_SPEED)
1791                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1792        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1793        {
1794                unsigned char *tmpCompBytes;
1795                size_t tmpOutSize;
1796                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1797                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1798                free(tmpCompBytes);
1799        }
1800        else
1801        {
1802                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1803        }
1804
1805        //TODO
1806//      if(*outSize>dataLength*sizeof(double))
1807//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1808
1809        free_TightDataPointStorageD(tdps);
1810}
1811
1812void SZ_compress_args_double_NoCkRnge_3D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
1813size_t r3, size_t r2, size_t r1, size_t s3, size_t s2, size_t s1, size_t e3, size_t e2, size_t e1)
1814{
1815        TightDataPointStorageD* tdps = SZ_compress_double_3D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r3, r2, r1, s3, s2, s1, e3, e2, e1);
1816
1817        if (confparams_cpr->szMode==SZ_BEST_SPEED)
1818                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1819        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1820        {
1821                unsigned char *tmpCompBytes;
1822                size_t tmpOutSize;
1823                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1824                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1825                free(tmpCompBytes);
1826        }
1827        else
1828        {
1829                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1830        }
1831
1832        //TODO
1833//      if(*outSize>dataLength*sizeof(double))
1834//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1835
1836        free_TightDataPointStorageD(tdps);
1837}
1838
1839void SZ_compress_args_double_NoCkRnge_4D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
1840size_t r4, size_t r3, size_t r2, size_t r1, size_t s4, size_t s3, size_t s2, size_t s1, size_t e4, size_t e3, size_t e2, size_t e1)
1841{
1842        TightDataPointStorageD* tdps = SZ_compress_double_4D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);
1843
1844        if (confparams_cpr->szMode==SZ_BEST_SPEED)
1845                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
1846        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1847        {
1848                unsigned char *tmpCompBytes;
1849                size_t tmpOutSize;
1850                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
1851                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
1852                free(tmpCompBytes);
1853        }
1854        else
1855        {
1856                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
1857        }
1858
1859        //TODO
1860//      if(*outSize>dataLength*sizeof(double))
1861//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1862
1863        free_TightDataPointStorageD(tdps);
1864}
1865
1866
1867unsigned int optimize_intervals_double_1D_subblock(double *oriData, double realPrecision, size_t r1, size_t s1, size_t e1)
1868{
1869        size_t dataLength = e1 - s1 + 1;
1870        oriData = oriData + s1;
1871
1872        size_t i = 0;
1873        unsigned long radiusIndex;
1874        double pred_value = 0, pred_err;
1875        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
1876        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
1877        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
1878        for(i=2;i<dataLength;i++)
1879        {
1880                if(i%confparams_cpr->sampleDistance==0)
1881                {
1882                        pred_value = 2*oriData[i-1] - oriData[i-2];
1883                        //pred_value = oriData[i-1];
1884                        pred_err = fabs(pred_value - oriData[i]);
1885                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
1886                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
1887                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
1888                        intervals[radiusIndex]++;
1889                }
1890        }
1891        //compute the appropriate number
1892        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
1893        size_t sum = 0;
1894        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
1895        {
1896                sum += intervals[i];
1897                if(sum>targetCount)
1898                        break;
1899        }
1900
1901        if(i>=confparams_cpr->maxRangeRadius)
1902                i = confparams_cpr->maxRangeRadius-1;
1903        unsigned int accIntervals = 2*(i+1);
1904        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
1905
1906        if(powerOf2<32)
1907                powerOf2 = 32;
1908
1909        free(intervals);
1910        return powerOf2;
1911}
1912
1913unsigned int optimize_intervals_double_2D_subblock(double *oriData, double realPrecision, size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
1914{
1915        size_t R1 = e1 - s1 + 1;
1916        size_t R2 = e2 - s2 + 1;
1917
1918        size_t i,j, index;
1919        unsigned long radiusIndex;
1920        double pred_value = 0, pred_err;
1921        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
1922        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
1923        size_t totalSampleSize = R1*R2/confparams_cpr->sampleDistance;
1924        for(i=s1+1;i<=e1;i++)
1925        {
1926                for(j=s2+1;j<=e2;j++)
1927                {
1928                        if((i+j)%confparams_cpr->sampleDistance==0)
1929                        {
1930                                index = i*r2+j;
1931                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
1932                                pred_err = fabs(pred_value - oriData[index]);
1933                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
1934                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
1935                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
1936                                intervals[radiusIndex]++;
1937                        }
1938                }
1939        }
1940        //compute the appropriate number
1941        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
1942        size_t sum = 0;
1943        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
1944        {
1945                sum += intervals[i];
1946                if(sum>targetCount)
1947                        break;
1948        }
1949        if(i>=confparams_cpr->maxRangeRadius)
1950                i = confparams_cpr->maxRangeRadius-1;
1951        unsigned int accIntervals = 2*(i+1);
1952        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
1953
1954        if(powerOf2<32)
1955                powerOf2 = 32;
1956
1957        free(intervals);
1958        return powerOf2;
1959}
1960
1961unsigned int optimize_intervals_double_3D_subblock(double *oriData, double realPrecision, size_t r1, size_t r2, size_t r3, size_t s1, size_t s2, size_t s3, size_t e1, size_t e2, size_t e3)
1962{
1963        size_t R1 = e1 - s1 + 1;
1964        size_t R2 = e2 - s2 + 1;
1965        size_t R3 = e3 - s3 + 1;
1966
1967        size_t r23 = r2*r3;
1968
1969        size_t i,j,k, index;
1970        unsigned long radiusIndex;
1971        double pred_value = 0, pred_err;
1972        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
1973        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
1974        size_t totalSampleSize = R1*R2*R3/confparams_cpr->sampleDistance;
1975        for(i=s1+1;i<=e1;i++)
1976        {
1977                for(j=s2+1;j<=e2;j++)
1978                {
1979                        for(k=s3+1;k<=e3;k++)
1980                        {
1981                                if((i+j+k)%confparams_cpr->sampleDistance==0)
1982                                {
1983                                        index = i*r23+j*r3+k;
1984                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23]
1985                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
1986                                        pred_err = fabs(pred_value - oriData[index]);
1987                                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
1988                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
1989                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
1990                                        intervals[radiusIndex]++;
1991                                }
1992                        }
1993
1994                }
1995        }
1996        //compute the appropriate number
1997        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
1998        size_t sum = 0;
1999        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2000        {
2001                sum += intervals[i];
2002                if(sum>targetCount)
2003                        break;
2004        }
2005        if(i>=confparams_cpr->maxRangeRadius)
2006                i = confparams_cpr->maxRangeRadius-1;
2007
2008        unsigned int accIntervals = 2*(i+1);
2009        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2010
2011        if(powerOf2<32)
2012                powerOf2 = 32;
2013
2014        free(intervals);
2015        return powerOf2;
2016}
2017
2018unsigned int optimize_intervals_double_4D_subblock(double *oriData, double realPrecision,
2019size_t r1, size_t r2, size_t r3, size_t r4, size_t s1, size_t s2, size_t s3, size_t s4, size_t e1, size_t e2, size_t e3, size_t e4)
2020{
2021        size_t R1 = e1 - s1 + 1;
2022        size_t R2 = e2 - s2 + 1;
2023        size_t R3 = e3 - s3 + 1;
2024        size_t R4 = e4 - s4 + 1;
2025
2026        size_t r34 = r3*r4;
2027        size_t r234 = r2*r3*r4;
2028
2029        size_t i,j,k,l, index;
2030        unsigned long radiusIndex;
2031        double pred_value = 0, pred_err;
2032        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
2033        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
2034        size_t totalSampleSize = R1*R2*R3*R4/confparams_cpr->sampleDistance;
2035        for(i=s1+1;i<=e1;i++)
2036        {
2037                for(j=s2+1;j<=e2;j++)
2038                {
2039                        for(k=s3+1;k<=e3;k++)
2040                        {
2041                                for(l=s4+1;l<=e4;l++)
2042                                {
2043                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
2044                                        {
2045                                                index = i*r234+j*r34+k*r4+l;
2046                                                pred_value = oriData[index-1] + oriData[index-r4] + oriData[index-r34]
2047                                                                - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
2048                                                pred_err = fabs(pred_value - oriData[index]);
2049                                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
2050                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
2051                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
2052                                                intervals[radiusIndex]++;
2053                                        }
2054                                }
2055                        }
2056
2057                }
2058        }
2059        //compute the appropriate number
2060        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
2061        size_t sum = 0;
2062        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2063        {
2064                sum += intervals[i];
2065                if(sum>targetCount)
2066                        break;
2067        }
2068        if(i>=confparams_cpr->maxRangeRadius)
2069                i = confparams_cpr->maxRangeRadius-1;
2070
2071        unsigned int accIntervals = 2*(i+1);
2072        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2073
2074        if(powerOf2<32)
2075                powerOf2 = 32;
2076
2077        free(intervals);
2078        return powerOf2;
2079}
2080
2081TightDataPointStorageD* SZ_compress_double_1D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2082size_t r1, size_t s1, size_t e1)
2083{
2084        size_t dataLength = e1 - s1 + 1;
2085
2086        unsigned int quantization_intervals;
2087        if(exe_params->optQuantMode==1)
2088                quantization_intervals = optimize_intervals_double_1D_subblock(oriData, realPrecision, r1, s1, e1);
2089        else
2090                quantization_intervals = exe_params->intvCapacity;
2091        updateQuantizationInfo(quantization_intervals);
2092
2093        size_t i; 
2094        int reqLength;
2095        double medianValue = medianValue_d;
2096        short radExpo = getExponent_double(valueRangeSize/2);
2097
2098        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2099
2100        int* type = (int*) malloc(dataLength*sizeof(int));
2101
2102        double* spaceFillingValue = oriData + s1; //
2103
2104        DynamicIntArray *exactLeadNumArray;
2105        new_DIA(&exactLeadNumArray, DynArrayInitLen);
2106
2107        DynamicByteArray *exactMidByteArray;
2108        new_DBA(&exactMidByteArray, DynArrayInitLen);
2109
2110        DynamicIntArray *resiBitArray;
2111        new_DIA(&resiBitArray, DynArrayInitLen);
2112
2113        type[0] = 0;
2114
2115        unsigned char preDataBytes[8];
2116        longToBytes_bigEndian(preDataBytes, 0);
2117
2118        int reqBytesLength = reqLength/8;
2119        int resiBitsLength = reqLength%8;
2120        double last3CmprsData[3] = {0};
2121
2122        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
2123        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
2124
2125        //add the first data
2126        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2127        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2128        memcpy(preDataBytes,vce->curBytes,8);
2129        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2130        listAdd_double(last3CmprsData, vce->data);
2131
2132        //add the second data
2133        type[1] = 0;
2134        compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2135        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2136        memcpy(preDataBytes,vce->curBytes,8);
2137        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2138        listAdd_double(last3CmprsData, vce->data);
2139
2140        int state;
2141        double checkRadius;
2142        double curData;
2143        double pred;
2144        double predAbsErr;
2145        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
2146        double interval = 2*realPrecision;
2147
2148        for(i=2;i<dataLength;i++)
2149        {
2150                //printf("%.30G\n",last3CmprsData[0]);
2151                curData = spaceFillingValue[i];
2152                pred = 2*last3CmprsData[0] - last3CmprsData[1];
2153                //pred = last3CmprsData[0];
2154                predAbsErr = fabs(curData - pred);
2155                if(predAbsErr<=checkRadius)
2156                {
2157                        state = (predAbsErr/realPrecision+1)/2;
2158                        if(curData>=pred)
2159                        {
2160                                type[i] = exe_params->intvRadius+state;
2161                                pred = pred + state*interval;
2162                        }
2163                        else //curData<pred
2164                        {
2165                                type[i] = exe_params->intvRadius-state;
2166                                pred = pred - state*interval;
2167                        }
2168                        listAdd_double(last3CmprsData, pred);
2169                        continue;
2170                }
2171
2172                //unpredictable data processing
2173                type[i] = 0;
2174                compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2175                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2176                memcpy(preDataBytes,vce->curBytes,8);
2177                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2178
2179                listAdd_double(last3CmprsData, vce->data);
2180        }//end of for
2181
2182        size_t exactDataNum = exactLeadNumArray->size;
2183
2184        TightDataPointStorageD* tdps;
2185
2186        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
2187                        type, exactMidByteArray->array, exactMidByteArray->size,
2188                        exactLeadNumArray->array,
2189                        resiBitArray->array, resiBitArray->size,
2190                        resiBitsLength,
2191                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2192
2193        //free memory
2194        free_DIA(exactLeadNumArray);
2195        free_DIA(resiBitArray);
2196        free(type);
2197        free(vce);
2198        free(lce);
2199        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2200
2201        return tdps;
2202}
2203
2204
2205TightDataPointStorageD* SZ_compress_double_2D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2206size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
2207{
2208        unsigned int quantization_intervals;
2209        if(exe_params->optQuantMode==1)
2210        {
2211                quantization_intervals = optimize_intervals_double_2D_subblock(oriData, realPrecision, r1, r2, s1, s2, e1, e2);
2212                updateQuantizationInfo(quantization_intervals);
2213        }
2214        else
2215                quantization_intervals = exe_params->intvCapacity;
2216
2217        size_t i,j; 
2218        int reqLength;
2219        double pred1D, pred2D;
2220        double diff = 0.0;
2221        double itvNum = 0;
2222        double *P0, *P1;
2223
2224        size_t R1 = e1 - s1 + 1;
2225        size_t R2 = e2 - s2 + 1;
2226        size_t dataLength = R1*R2;
2227
2228        P0 = (double*)malloc(R2*sizeof(double));
2229        memset(P0, 0, R2*sizeof(double));
2230        P1 = (double*)malloc(R2*sizeof(double));
2231        memset(P1, 0, R2*sizeof(double));
2232
2233        double medianValue = medianValue_d;
2234        short radExpo = getExponent_double(valueRangeSize/2);
2235        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2236
2237        int* type = (int*) malloc(dataLength*sizeof(int));
2238
2239        double* spaceFillingValue = oriData; //
2240
2241        DynamicIntArray *exactLeadNumArray;
2242        new_DIA(&exactLeadNumArray, DynArrayInitLen);
2243
2244        DynamicByteArray *exactMidByteArray;
2245        new_DBA(&exactMidByteArray, DynArrayInitLen);
2246
2247        DynamicIntArray *resiBitArray;
2248        new_DIA(&resiBitArray, DynArrayInitLen);
2249
2250        unsigned char preDataBytes[8];
2251        longToBytes_bigEndian(preDataBytes, 0);
2252
2253        int reqBytesLength = reqLength/8;
2254        int resiBitsLength = reqLength%8;
2255
2256        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
2257        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
2258
2259        /* Process Row-s1 data s2*/
2260        size_t gIndex;
2261        size_t lIndex;
2262
2263        gIndex = s1*r2+s2;
2264        lIndex = 0;
2265
2266        type[lIndex] = 0;
2267        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2268        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2269        memcpy(preDataBytes,vce->curBytes,8);
2270        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2271        P1[0] = vce->data;
2272
2273        /* Process Row-s1 data s2+1*/
2274        gIndex = s1*r2+(s2+1);
2275        lIndex = 1;
2276
2277        pred1D = P1[0];
2278        diff = spaceFillingValue[gIndex] - pred1D;
2279
2280        itvNum =  fabs(diff)/realPrecision + 1;
2281
2282        if (itvNum < exe_params->intvCapacity)
2283        {
2284                if (diff < 0) itvNum = -itvNum;
2285                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2286                P1[1] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2287        }
2288        else
2289        {
2290                type[lIndex] = 0;
2291                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2292                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2293                memcpy(preDataBytes,vce->curBytes,8);
2294                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2295                P1[1] = vce->data;
2296        }
2297
2298    /* Process Row-s1 data s2+2 --> data e2 */
2299        for (j = 2; j < R2; j++)
2300        {
2301                gIndex = s1*r2+(s2+j);
2302                lIndex = j;
2303
2304                pred1D = 2*P1[j-1] - P1[j-2];
2305                diff = spaceFillingValue[gIndex] - pred1D;
2306
2307                itvNum = fabs(diff)/realPrecision + 1;
2308
2309                if (itvNum < exe_params->intvCapacity)
2310                {
2311                        if (diff < 0) itvNum = -itvNum;
2312                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2313                        P1[j] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2314                }
2315                else
2316                {
2317                        type[lIndex] = 0;
2318                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2319                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2320                        memcpy(preDataBytes,vce->curBytes,8);
2321                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2322                        P1[j] = vce->data;
2323                }
2324        }
2325
2326        /* Process Row-s1+1 --> Row-e1 */
2327        for (i = 1; i < R1; i++)
2328        {
2329                /* Process row-s1+i data s2 */
2330                gIndex = (s1+i)*r2+s2;
2331                lIndex = i*R2;
2332
2333                pred1D = P1[0];
2334                diff = spaceFillingValue[gIndex] - pred1D;
2335
2336                itvNum = fabs(diff)/realPrecision + 1;
2337
2338                if (itvNum < exe_params->intvCapacity)
2339                {
2340                        if (diff < 0) itvNum = -itvNum;
2341                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2342                        P0[0] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2343                }
2344                else
2345                {
2346                        type[lIndex] = 0;
2347                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2348                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2349                        memcpy(preDataBytes,vce->curBytes,8);
2350                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2351                        P0[0] = vce->data;
2352                }
2353
2354                /* Process row-s1+i data s2+1 --> e2 */
2355                for (j = 1; j < R2; j++)
2356                {
2357                        gIndex = (s1+i)*r2+(s2+j);
2358                        lIndex = i*R2+j;
2359
2360                        pred2D = P0[j-1] + P1[j] - P1[j-1];
2361                        diff = spaceFillingValue[gIndex] - pred2D;
2362
2363                        itvNum = fabs(diff)/realPrecision + 1;
2364
2365                        if (itvNum < exe_params->intvCapacity)
2366                        {
2367                                if (diff < 0) itvNum = -itvNum;
2368                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2369                                P0[j] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2370                        }
2371                        else
2372                        {
2373                                type[lIndex] = 0;
2374                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2375                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2376                                memcpy(preDataBytes,vce->curBytes,8);
2377                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2378                                P0[j] = vce->data;
2379                        }
2380                }
2381
2382                double *Pt;
2383                Pt = P1;
2384                P1 = P0;
2385                P0 = Pt;
2386        }
2387
2388        free(P0);
2389        free(P1);
2390        size_t exactDataNum = exactLeadNumArray->size;
2391
2392        TightDataPointStorageD* tdps;
2393
2394        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
2395                        type, exactMidByteArray->array, exactMidByteArray->size,
2396                        exactLeadNumArray->array,
2397                        resiBitArray->array, resiBitArray->size,
2398                        resiBitsLength,
2399                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2400
2401        //free memory
2402        free_DIA(exactLeadNumArray);
2403        free_DIA(resiBitArray);
2404        free(type);
2405        free(vce);
2406        free(lce);
2407        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2408
2409        return tdps;
2410}
2411
2412TightDataPointStorageD* SZ_compress_double_3D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2413size_t r1, size_t r2, size_t r3, size_t s1, size_t s2, size_t s3, size_t e1, size_t e2, size_t e3)
2414{
2415        unsigned int quantization_intervals;
2416        if(exe_params->optQuantMode==1)
2417        {
2418                quantization_intervals = optimize_intervals_double_3D_subblock(oriData, realPrecision, r1, r2, r3, s1, s2, s3, e1, e2, e3);
2419                updateQuantizationInfo(quantization_intervals);
2420        }
2421        else
2422                quantization_intervals = exe_params->intvCapacity;
2423
2424        size_t i,j,k; 
2425        int reqLength;
2426        double pred1D, pred2D, pred3D;
2427        double diff = 0.0;
2428        double itvNum = 0;
2429        double *P0, *P1;
2430
2431        size_t R1 = e1 - s1 + 1;
2432        size_t R2 = e2 - s2 + 1;
2433        size_t R3 = e3 - s3 + 1;
2434        size_t dataLength = R1*R2*R3;
2435
2436        size_t r23 = r2*r3;
2437        size_t R23 = R2*R3;
2438
2439        P0 = (double*)malloc(R23*sizeof(double));
2440        P1 = (double*)malloc(R23*sizeof(double));
2441
2442        double medianValue = medianValue_d;
2443        short radExpo = getExponent_double(valueRangeSize/2);
2444        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2445
2446        int* type = (int*) malloc(dataLength*sizeof(int));
2447
2448        double* spaceFillingValue = oriData; //
2449
2450        DynamicIntArray *exactLeadNumArray;
2451        new_DIA(&exactLeadNumArray, DynArrayInitLen);
2452
2453        DynamicByteArray *exactMidByteArray;
2454        new_DBA(&exactMidByteArray, DynArrayInitLen);
2455
2456        DynamicIntArray *resiBitArray;
2457        new_DIA(&resiBitArray, DynArrayInitLen);
2458
2459        unsigned char preDataBytes[8];
2460        longToBytes_bigEndian(preDataBytes, 0);
2461
2462        int reqBytesLength = reqLength/8;
2463        int resiBitsLength = reqLength%8;
2464
2465        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
2466        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
2467
2468
2469        ///////////////////////////     Process layer-s1 ///////////////////////////
2470        /* Process Row-s2 data s3*/
2471        size_t gIndex;  //global index
2472        size_t lIndex;  //local index
2473        size_t index2D;         //local 2D index
2474
2475        gIndex = s1*r23+s2*r3+s3;
2476        lIndex = 0;
2477        index2D = 0;
2478
2479        type[lIndex] = 0;
2480        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2481        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2482        memcpy(preDataBytes,vce->curBytes,8);
2483        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2484        P1[index2D] = vce->data;
2485
2486        /* Process Row-s2 data s3+1*/
2487        gIndex = s1*r23+s2*r3+s3+1;
2488        lIndex = 1;
2489        index2D = 1;
2490
2491        pred1D = P1[index2D-1];
2492        diff = spaceFillingValue[gIndex] - pred1D;
2493
2494        itvNum = fabs(diff)/realPrecision + 1;
2495
2496        if (itvNum < exe_params->intvCapacity)
2497        {
2498                if (diff < 0) itvNum = -itvNum;
2499                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2500                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2501        }
2502        else
2503        {
2504                type[lIndex] = 0;
2505                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2506                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2507                memcpy(preDataBytes,vce->curBytes,8);
2508                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2509                P1[index2D] = vce->data;
2510        }
2511
2512    /* Process Row-s2 data s3+2 --> data e3 */
2513        for (j = 2; j < R3; j++)
2514        {
2515                gIndex = s1*r23+s2*r3+s3+j;
2516                lIndex = j;
2517                index2D = j;
2518
2519                pred1D = 2*P1[index2D-1] - P1[index2D-2];
2520                diff = spaceFillingValue[gIndex] - pred1D;
2521
2522                itvNum = fabs(diff)/realPrecision + 1;
2523
2524                if (itvNum < exe_params->intvCapacity)
2525                {
2526                        if (diff < 0) itvNum = -itvNum;
2527                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2528                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2529                }
2530                else
2531                {
2532                        type[lIndex] = 0;
2533                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2534                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2535                        memcpy(preDataBytes,vce->curBytes,8);
2536                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2537                        P1[index2D] = vce->data;
2538                }
2539        }
2540
2541        /* Process Row-s2+1 --> Row-e2 */
2542        for (i = 1; i < R2; i++)
2543        {
2544                /* Process row-s2+i data s3 */
2545                gIndex = s1*r23+(s2+i)*r3+s3;
2546                lIndex = i*R3;
2547                index2D = i*R3;
2548
2549                pred1D  = P1[index2D-R3];
2550                diff    = spaceFillingValue[gIndex] - pred1D;
2551
2552                itvNum  = fabs(diff)/realPrecision + 1;
2553
2554                if (itvNum < exe_params->intvCapacity)
2555                {
2556                        if (diff < 0) itvNum = -itvNum;
2557                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2558                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2559                }
2560                else
2561                {
2562                        type[lIndex] = 0;
2563                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2564                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2565                        memcpy(preDataBytes,vce->curBytes,8);
2566                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2567                        P1[index2D] = vce->data;
2568                }
2569
2570                /* Process row-s2+i data s3+1 --> data e3*/
2571                for (j = 1; j < R3; j++)
2572                {
2573                        gIndex = s1*r23+(s2+i)*r3+s3+j;
2574                        lIndex = i*R3+j;
2575                        index2D = i*R3+j;
2576
2577                        pred2D  = P1[index2D-1] + P1[index2D-R3] - P1[index2D-R3-1];
2578                        diff = spaceFillingValue[gIndex] - pred2D;
2579
2580                        itvNum = fabs(diff)/realPrecision + 1;
2581
2582                        if (itvNum < exe_params->intvCapacity)
2583                        {
2584                                if (diff < 0) itvNum = -itvNum;
2585                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2586                                P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2587                        }
2588                        else
2589                        {
2590                                type[lIndex] = 0;
2591                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2592                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2593                                memcpy(preDataBytes,vce->curBytes,8);
2594                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2595                                P1[index2D] = vce->data;
2596                        }
2597                }
2598        }
2599
2600
2601        ///////////////////////////     Process layer-s1+1 --> layer-e1 ///////////////////////////
2602
2603        for (k = 1; k < R1; k++)
2604        {
2605                /* Process Row-s2 data s3*/
2606                gIndex = (s1+k)*r23+s2*r3+s3;
2607                lIndex = k*R23;
2608                index2D = 0;
2609
2610                pred1D = P1[index2D];
2611                diff = spaceFillingValue[gIndex] - pred1D;
2612
2613                itvNum = fabs(diff)/realPrecision + 1;
2614
2615                if (itvNum < exe_params->intvCapacity)
2616                {
2617                        if (diff < 0) itvNum = -itvNum;
2618                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2619                        P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2620                }
2621                else
2622                {
2623                        type[lIndex] = 0;
2624                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2625                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2626                        memcpy(preDataBytes,vce->curBytes,8);
2627                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2628                        P0[index2D] = vce->data;
2629                }
2630
2631
2632            /* Process Row-s2 data s3+1 --> data e3 */
2633                for (j = 1; j < R3; j++)
2634                {
2635                        gIndex = (s1+k)*r23+s2*r3+s3+j;
2636                        lIndex = k*R23+j;
2637                        index2D = j;
2638
2639                        pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
2640                        diff = spaceFillingValue[gIndex] - pred2D;
2641
2642                        itvNum = fabs(diff)/realPrecision + 1;
2643
2644                        if (itvNum < exe_params->intvCapacity)
2645                        {
2646                                if (diff < 0) itvNum = -itvNum;
2647                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2648                                P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2649                        }
2650                        else
2651                        {
2652                                type[lIndex] = 0;
2653                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2654                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2655                                memcpy(preDataBytes,vce->curBytes,8);
2656                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2657                                P0[index2D] = vce->data;
2658                        }
2659                }
2660
2661            /* Process Row-s2+1 --> Row-e2 */
2662                for (i = 1; i < R2; i++)
2663                {
2664                        /* Process Row-s2+i data s3 */
2665                        gIndex = (s1+k)*r23+(s2+i)*r3+s3;
2666                        lIndex = k*R23+i*R3;
2667                        index2D = i*R3;
2668
2669                        pred2D = P0[index2D-R3] + P1[index2D] - P1[index2D-R3];
2670                        diff = spaceFillingValue[gIndex] - pred2D;
2671
2672                        itvNum = fabs(diff)/realPrecision + 1;
2673
2674                        if (itvNum < exe_params->intvCapacity)
2675                        {
2676                                if (diff < 0) itvNum = -itvNum;
2677                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2678                                P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2679                        }
2680                        else
2681                        {
2682                                type[lIndex] = 0;
2683                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2684                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2685                                memcpy(preDataBytes,vce->curBytes,8);
2686                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2687                                P0[index2D] = vce->data;
2688                        }
2689
2690                        /* Process Row-s2+i data s3+1 --> data e3 */
2691                        for (j = 1; j < R3; j++)
2692                        {
2693                                gIndex = (s1+k)*r23+(s2+i)*r3+s3+j;
2694                                lIndex = k*R23+i*R3+j;
2695                                index2D = i*R3+j;
2696
2697                                pred3D = P0[index2D-1] + P0[index2D-R3]+ P1[index2D] - P0[index2D-R3-1] - P1[index2D-R3] - P1[index2D-1] + P1[index2D-R3-1];
2698                                diff = spaceFillingValue[gIndex] - pred3D;
2699
2700                                itvNum = fabs(diff)/realPrecision + 1;
2701
2702                                if (itvNum < exe_params->intvCapacity)
2703                                {
2704                                        if (diff < 0) itvNum = -itvNum;
2705                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2706                                        P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2707                                }
2708                                else
2709                                {
2710                                        type[lIndex] = 0;
2711                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2712                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2713                                        memcpy(preDataBytes,vce->curBytes,8);
2714                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2715                                        P0[index2D] = vce->data;
2716                                }
2717                        }
2718                }
2719
2720                double *Pt;
2721                Pt = P1;
2722                P1 = P0;
2723                P0 = Pt;
2724        }
2725
2726        free(P0);
2727        free(P1);
2728        size_t exactDataNum = exactLeadNumArray->size;
2729
2730        TightDataPointStorageD* tdps;
2731
2732        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
2733                        type, exactMidByteArray->array, exactMidByteArray->size,
2734                        exactLeadNumArray->array,
2735                        resiBitArray->array, resiBitArray->size,
2736                        resiBitsLength,
2737                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2738
2739        //free memory
2740        free_DIA(exactLeadNumArray);
2741        free_DIA(resiBitArray);
2742        free(type);
2743        free(vce);
2744        free(lce);
2745        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2746
2747        return tdps;
2748}
2749
2750TightDataPointStorageD* SZ_compress_double_4D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
2751size_t r1, size_t r2, size_t r3, size_t r4, size_t s1, size_t s2, size_t s3, size_t s4, size_t e1, size_t e2, size_t e3, size_t e4)
2752{
2753        unsigned int quantization_intervals;
2754        if(exe_params->optQuantMode==1)
2755        {
2756                quantization_intervals = optimize_intervals_double_4D_subblock(oriData, realPrecision, r1, r2, r3, r4, s1, s2, s3, s4, e1, e2, e3, e4);
2757                updateQuantizationInfo(quantization_intervals);
2758        }
2759        else
2760                quantization_intervals = exe_params->intvCapacity;
2761
2762        size_t i,j,k; 
2763        int reqLength;
2764        double pred1D, pred2D, pred3D;
2765        double diff = 0.0;
2766        double itvNum = 0;
2767        double *P0, *P1;
2768
2769        size_t R1 = e1 - s1 + 1;
2770        size_t R2 = e2 - s2 + 1;
2771        size_t R3 = e3 - s3 + 1;
2772        size_t R4 = e4 - s4 + 1;
2773
2774        size_t dataLength = R1*R2*R3*R4;
2775
2776        size_t r34 = r3*r4;
2777        size_t r234 = r2*r3*r4;
2778        size_t R34 = R3*R4;
2779        size_t R234 = R2*R3*R4;
2780
2781        P0 = (double*)malloc(R34*sizeof(double));
2782        P1 = (double*)malloc(R34*sizeof(double));
2783
2784        double medianValue = medianValue_d;
2785        short radExpo = getExponent_double(valueRangeSize/2);
2786        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
2787
2788        int* type = (int*) malloc(dataLength*sizeof(int));
2789
2790        double* spaceFillingValue = oriData; //
2791
2792        DynamicIntArray *exactLeadNumArray;
2793        new_DIA(&exactLeadNumArray, DynArrayInitLen);
2794
2795        DynamicByteArray *exactMidByteArray;
2796        new_DBA(&exactMidByteArray, DynArrayInitLen);
2797
2798        DynamicIntArray *resiBitArray;
2799        new_DIA(&resiBitArray, DynArrayInitLen);
2800
2801        unsigned char preDataBytes[8];
2802        longToBytes_bigEndian(preDataBytes, 0);
2803
2804        int reqBytesLength = reqLength/8;
2805        int resiBitsLength = reqLength%8;
2806
2807        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
2808        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
2809
2810        size_t l;
2811        for (l = 0; l < R1; l++)
2812        {
2813
2814                ///////////////////////////     Process layer-s2 ///////////////////////////
2815                /* Process Row-s3 data s4*/
2816                size_t gIndex;  //global index
2817                size_t lIndex;  //local index
2818                size_t index2D;         //local 2D index
2819
2820                gIndex = (s1+l)*r234+s2*r34+s3*r4+s4;
2821                lIndex = l*R234;
2822                index2D = 0;
2823
2824                type[lIndex] = 0;
2825                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2826                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2827                memcpy(preDataBytes,vce->curBytes,8);
2828                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2829                P1[index2D] = vce->data;
2830
2831                /* Process Row-s3 data s4+1*/
2832                gIndex = (s1+l)*r234+s2*r34+s3*r4+s4+1;
2833                lIndex = l*R234+1;
2834                index2D = 1;
2835
2836                pred1D = P1[index2D-1];
2837                diff = spaceFillingValue[gIndex] - pred1D;
2838
2839                itvNum = fabs(diff)/realPrecision + 1;
2840
2841                if (itvNum < exe_params->intvCapacity)
2842                {
2843                        if (diff < 0) itvNum = -itvNum;
2844                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2845                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2846                }
2847                else
2848                {
2849                        type[lIndex] = 0;
2850                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2851                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2852                        memcpy(preDataBytes,vce->curBytes,8);
2853                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2854                        P1[index2D] = vce->data;
2855                }
2856
2857                /* Process Row-s3 data s4+2 --> data e4 */
2858                for (j = 2; j < R4; j++)
2859                {
2860                        gIndex = (s1+l)*r234+s2*r34+s3*r4+s4+j;
2861                        lIndex = l*R234+j;
2862                        index2D = j;
2863
2864                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
2865                        diff = spaceFillingValue[gIndex] - pred1D;
2866
2867                        itvNum = fabs(diff)/realPrecision + 1;
2868
2869                        if (itvNum < exe_params->intvCapacity)
2870                        {
2871                                if (diff < 0) itvNum = -itvNum;
2872                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2873                                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2874                        }
2875                        else
2876                        {
2877                                type[lIndex] = 0;
2878                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2879                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2880                                memcpy(preDataBytes,vce->curBytes,8);
2881                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2882                                P1[index2D] = vce->data;
2883                        }
2884                }
2885
2886                /* Process Row-s3+1 --> Row-e3 */
2887                for (i = 1; i < R3; i++)
2888                {
2889                        /* Process row-s2+i data s3 */
2890                        gIndex = (s1+l)*r234+s2*r34+(s3+i)*r4+s4;
2891                        lIndex = l*R234+i*R4;
2892                        index2D = i*R4;
2893
2894                        pred1D  = P1[index2D-R4];
2895                        diff    = spaceFillingValue[gIndex] - pred1D;
2896
2897                        itvNum  = fabs(diff)/realPrecision + 1;
2898
2899                        if (itvNum < exe_params->intvCapacity)
2900                        {
2901                                if (diff < 0) itvNum = -itvNum;
2902                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2903                                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2904                        }
2905                        else
2906                        {
2907                                type[lIndex] = 0;
2908                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2909                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2910                                memcpy(preDataBytes,vce->curBytes,8);
2911                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2912                                P1[index2D] = vce->data;
2913                        }
2914
2915                        /* Process row-s3+i data s4+1 --> data e4*/
2916                        for (j = 1; j < R4; j++)
2917                        {
2918                                gIndex = (s1+l)*r234+s2*r34+(s3+i)*r4+s4+j;
2919                                lIndex = l*R234+i*R4+j;
2920                                index2D = i*R4+j;
2921
2922                                pred2D  = P1[index2D-1] + P1[index2D-R4] - P1[index2D-R4-1];
2923                                diff = spaceFillingValue[gIndex] - pred2D;
2924
2925                                itvNum = fabs(diff)/realPrecision + 1;
2926
2927                                if (itvNum < exe_params->intvCapacity)
2928                                {
2929                                        if (diff < 0) itvNum = -itvNum;
2930                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2931                                        P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2932                                }
2933                                else
2934                                {
2935                                        type[lIndex] = 0;
2936                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2937                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2938                                        memcpy(preDataBytes,vce->curBytes,8);
2939                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2940                                        P1[index2D] = vce->data;
2941                                }
2942                        }
2943                }
2944
2945
2946                ///////////////////////////     Process layer-s2+1 --> layer-e2 ///////////////////////////
2947
2948                for (k = 1; k < R2; k++)
2949                {
2950                        /* Process Row-s3 data s4*/
2951                        gIndex = (s1+l)*r234+(s2+k)*r34+s3*r4+s4;
2952                        lIndex = l*R234+k*R34;
2953                        index2D = 0;
2954
2955                        pred1D = P1[index2D];
2956                        diff = spaceFillingValue[gIndex] - pred1D;
2957
2958                        itvNum = fabs(diff)/realPrecision + 1;
2959
2960                        if (itvNum < exe_params->intvCapacity)
2961                        {
2962                                if (diff < 0) itvNum = -itvNum;
2963                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2964                                P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2965                        }
2966                        else
2967                        {
2968                                type[lIndex] = 0;
2969                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2970                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2971                                memcpy(preDataBytes,vce->curBytes,8);
2972                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2973                                P0[index2D] = vce->data;
2974                        }
2975
2976
2977                        /* Process Row-s3 data s4+1 --> data e4 */
2978                        for (j = 1; j < R4; j++)
2979                        {
2980                                gIndex = (s1+l)*r234+(s2+k)*r34+s3*r4+s4+j;
2981                                lIndex = l*R234+k*R34+j;
2982                                index2D = j;
2983
2984                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
2985                                diff = spaceFillingValue[gIndex] - pred2D;
2986
2987                                itvNum = fabs(diff)/realPrecision + 1;
2988
2989                                if (itvNum < exe_params->intvCapacity)
2990                                {
2991                                        if (diff < 0) itvNum = -itvNum;
2992                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2993                                        P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2994                                }
2995                                else
2996                                {
2997                                        type[lIndex] = 0;
2998                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2999                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3000                                        memcpy(preDataBytes,vce->curBytes,8);
3001                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3002                                        P0[index2D] = vce->data;
3003                                }
3004                        }
3005
3006                        /* Process Row-s3+1 --> Row-e3 */
3007                        for (i = 1; i < R3; i++)
3008                        {
3009                                /* Process Row-s3+i data s4 */
3010                                gIndex = (s1+l)*r234+(s2+k)*r34+(s3+i)*r4+s4;
3011                                lIndex = l*R234+k*R34+i*R4;
3012                                index2D = i*R4;
3013
3014                                pred2D = P0[index2D-R4] + P1[index2D] - P1[index2D-R4];
3015                                diff = spaceFillingValue[gIndex] - pred2D;
3016
3017                                itvNum = fabs(diff)/realPrecision + 1;
3018
3019                                if (itvNum < exe_params->intvCapacity)
3020                                {
3021                                        if (diff < 0) itvNum = -itvNum;
3022                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3023                                        P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3024                                }
3025                                else
3026                                {
3027                                        type[lIndex] = 0;
3028                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3029                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3030                                        memcpy(preDataBytes,vce->curBytes,8);
3031                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3032                                        P0[index2D] = vce->data;
3033                                }
3034
3035                                /* Process Row-s3+i data s4+1 --> data e4 */
3036                                for (j = 1; j < R4; j++)
3037                                {
3038                                        gIndex = (s1+l)*r234+(s2+k)*r34+(s3+i)*r4+s4+j;
3039                                        lIndex = l*R234+k*R34+i*R4+j;
3040                                        index2D = i*R4+j;
3041
3042//                                      printf ("global index = %d, local index = %d\n", gIndex, lIndex);
3043
3044                                        pred3D = P0[index2D-1] + P0[index2D-R4]+ P1[index2D] - P0[index2D-R4-1] - P1[index2D-R4] - P1[index2D-1] + P1[index2D-R4-1];
3045                                        diff = spaceFillingValue[gIndex] - pred3D;
3046
3047                                        itvNum = fabs(diff)/realPrecision + 1;
3048
3049                                        if (itvNum < exe_params->intvCapacity)
3050                                        {
3051                                                if (diff < 0) itvNum = -itvNum;
3052                                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3053                                                P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3054                                        }
3055                                        else
3056                                        {
3057                                                type[lIndex] = 0;
3058                                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3059                                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3060                                                memcpy(preDataBytes,vce->curBytes,8);
3061                                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3062                                                P0[index2D] = vce->data;
3063                                        }
3064                                }
3065                        }
3066
3067                        double *Pt;
3068                        Pt = P1;
3069                        P1 = P0;
3070                        P0 = Pt;
3071                }
3072        }
3073
3074        free(P0);
3075        free(P1);
3076        size_t exactDataNum = exactLeadNumArray->size;
3077
3078        TightDataPointStorageD* tdps;
3079
3080        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
3081                        type, exactMidByteArray->array, exactMidByteArray->size,
3082                        exactLeadNumArray->array,
3083                        resiBitArray->array, resiBitArray->size,
3084                        resiBitsLength,
3085                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
3086
3087        //free memory
3088        free_DIA(exactLeadNumArray);
3089        free_DIA(resiBitArray);
3090        free(type);
3091        free(vce);
3092        free(lce);
3093        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
3094
3095        return tdps;
3096}
3097
3098/**
3099 *
3100 * This is a fast implementation for optimize_intervals_double_3D()
3101 * */
3102unsigned int optimize_intervals_double_3D_opt(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision){ 
3103        size_t i;
3104        size_t radiusIndex;
3105        size_t r23=r2*r3;
3106        double pred_value = 0, pred_err;
3107        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3108        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3109        size_t totalSampleSize = 0;
3110
3111        size_t offset_count = confparams_cpr->sampleDistance - 2; // count r3 offset
3112        size_t offset_count_2;
3113        double * data_pos = oriData + r23 + r3 + offset_count;
3114        size_t n1_count = 1, n2_count = 1; // count i,j sum
3115        size_t len = r1 * r2 * r3;
3116        while(data_pos - oriData < len){
3117                totalSampleSize++;
3118                pred_value = data_pos[-1] + data_pos[-r3] + data_pos[-r23] - data_pos[-1-r23] - data_pos[-r3-1] - data_pos[-r3-r23] + data_pos[-r3-r23-1];
3119                pred_err = fabs(pred_value - *data_pos);
3120                radiusIndex = (pred_err/realPrecision+1)/2;
3121                if(radiusIndex>=confparams_cpr->maxRangeRadius)
3122                {
3123                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
3124                        //printf("radiusIndex=%d\n", radiusIndex);
3125                }
3126                intervals[radiusIndex]++;
3127                // printf("TEST: %ld, i: %ld\tj: %ld\tk: %ld\n", data_pos - oriData);
3128                // fflush(stdout);
3129                offset_count += confparams_cpr->sampleDistance;
3130                if(offset_count >= r3){
3131                        n2_count ++;
3132                        if(n2_count == r2){
3133                                n1_count ++;
3134                                n2_count = 1;
3135                                data_pos += r3;
3136                        }
3137                        offset_count_2 = (n1_count + n2_count) % confparams_cpr->sampleDistance;
3138                        data_pos += (r3 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
3139                        offset_count = (confparams_cpr->sampleDistance - offset_count_2);
3140                        if(offset_count == 0) offset_count ++;
3141                }
3142                else data_pos += confparams_cpr->sampleDistance;
3143        }       
3144        // printf("sample_count: %ld\n", sample_count);
3145        // fflush(stdout);
3146        // if(*max_freq < 0.15) *max_freq *= 2;
3147        //compute the appropriate number
3148        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3149        size_t sum = 0;
3150        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3151        {
3152                sum += intervals[i];
3153                if(sum>targetCount)
3154                        break;
3155        }
3156        if(i>=confparams_cpr->maxRangeRadius)
3157                i = confparams_cpr->maxRangeRadius-1;
3158        unsigned int accIntervals = 2*(i+1);
3159        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3160
3161        if(powerOf2<32)
3162                powerOf2 = 32;
3163        free(intervals);
3164        //printf("targetCount=%d, sum=%d, totalSampleSize=%d, ratio=%f, accIntervals=%d, powerOf2=%d\n", targetCount, sum, totalSampleSize, (double)sum/(double)totalSampleSize, accIntervals, powerOf2);
3165        return powerOf2;
3166}
3167
3168unsigned int optimize_intervals_double_2D_opt(double *oriData, size_t r1, size_t r2, double realPrecision)
3169{       
3170        size_t i;
3171        size_t radiusIndex;
3172        double pred_value = 0, pred_err;
3173        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3174        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3175        size_t totalSampleSize = 0;//(r1-1)*(r2-1)/confparams_cpr->sampleDistance;
3176
3177        size_t offset_count = confparams_cpr->sampleDistance - 1; // count r2 offset
3178        size_t offset_count_2;
3179        double * data_pos = oriData + r2 + offset_count;
3180        size_t n1_count = 1; // count i sum
3181        size_t len = r1 * r2;
3182        while(data_pos - oriData < len){
3183                totalSampleSize++;
3184                pred_value = data_pos[-1] + data_pos[-r2] - data_pos[-r2-1];
3185                pred_err = fabs(pred_value - *data_pos);
3186                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
3187                if(radiusIndex>=confparams_cpr->maxRangeRadius)
3188                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
3189                intervals[radiusIndex]++;
3190
3191                offset_count += confparams_cpr->sampleDistance;
3192                if(offset_count >= r2){
3193                        n1_count ++;
3194                        offset_count_2 = n1_count % confparams_cpr->sampleDistance;
3195                        data_pos += (r2 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
3196                        offset_count = (confparams_cpr->sampleDistance - offset_count_2);
3197                        if(offset_count == 0) offset_count ++;
3198                }
3199                else data_pos += confparams_cpr->sampleDistance;
3200        }
3201
3202        //compute the appropriate number
3203        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3204        size_t sum = 0;
3205        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3206        {
3207                sum += intervals[i];
3208                if(sum>targetCount)
3209                        break;
3210        }
3211        if(i>=confparams_cpr->maxRangeRadius)
3212                i = confparams_cpr->maxRangeRadius-1;
3213        unsigned int accIntervals = 2*(i+1);
3214        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3215
3216        if(powerOf2<32)
3217                powerOf2 = 32;
3218
3219        free(intervals);
3220        return powerOf2;
3221}
3222
3223unsigned int optimize_intervals_double_1D_opt(double *oriData, size_t dataLength, double realPrecision)
3224{       
3225        size_t i = 0, radiusIndex;
3226        double pred_value = 0, pred_err;
3227        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3228        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3229        size_t totalSampleSize = 0;//dataLength/confparams_cpr->sampleDistance;
3230
3231        double * data_pos = oriData + 2;
3232        while(data_pos - oriData < dataLength){
3233                totalSampleSize++;
3234                //pred_value = 2*data_pos[-1] - data_pos[-2];
3235                pred_value = data_pos[-1];
3236                pred_err = fabs(pred_value - *data_pos);
3237                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
3238                if(radiusIndex>=confparams_cpr->maxRangeRadius)
3239                        radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
3240                intervals[radiusIndex]++;
3241
3242                data_pos += confparams_cpr->sampleDistance;
3243        }
3244        //compute the appropriate number
3245        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3246        size_t sum = 0;
3247        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3248        {
3249                sum += intervals[i];
3250                if(sum>targetCount)
3251                        break;
3252        }
3253        if(i>=confparams_cpr->maxRangeRadius)
3254                i = confparams_cpr->maxRangeRadius-1;
3255               
3256        unsigned int accIntervals = 2*(i+1);
3257        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3258       
3259        if(powerOf2<32)
3260                powerOf2 = 32;
3261       
3262        free(intervals);
3263        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
3264        return powerOf2;
3265}
Note: See TracBrowser for help on using the repository browser.