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

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

more work on adding SZ (latest version)
  • Property mode set to 100644
Line 
1/**
2 *  @file sz_float.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 "TightDataPointStorageF.h"
21#include "sz_float.h"
22#include "sz_float_pwr.h"
23#include "szd_float.h"
24#include "szd_float_pwr.h"
25#include "zlib.h"
26#include "rw.h"
27#include "sz_float_ts.h"
28
29unsigned char* SZ_skip_compress_float(float* data, size_t dataLength, size_t* outSize)
30{
31        *outSize = dataLength*sizeof(float);
32        unsigned char* out = (unsigned char*)malloc(dataLength*sizeof(float));
33        memcpy(out, data, dataLength*sizeof(float));
34        return out;
35}
36unsigned int optimize_intervals_float_1D(float *oriData, size_t dataLength, double realPrecision)
37{       
38        size_t i = 0, radiusIndex;
39        float pred_value = 0, pred_err;
40        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
41        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
42        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
43        for(i=2;i<dataLength;i++)
44        {
45                if(i%confparams_cpr->sampleDistance==0)
46                {
47                        //pred_value = 2*oriData[i-1] - oriData[i-2];
48                        pred_value = oriData[i-1];
49                        pred_err = fabs(pred_value - oriData[i]);
50                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
51                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
52                                radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
53                        intervals[radiusIndex]++;
54                }
55        }
56        //compute the appropriate number
57        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
58        size_t sum = 0;
59        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
60        {
61                sum += intervals[i];
62                if(sum>targetCount)
63                        break;
64        }
65        if(i>=confparams_cpr->maxRangeRadius)
66                i = confparams_cpr->maxRangeRadius-1;
67               
68        unsigned int accIntervals = 2*(i+1);
69        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
70       
71        if(powerOf2<32)
72                powerOf2 = 32;
73       
74        free(intervals);
75        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
76        return powerOf2;
77}
78
79unsigned int optimize_intervals_float_2D(float *oriData, size_t r1, size_t r2, double realPrecision)
80{       
81        size_t i,j, index;
82        size_t radiusIndex;
83        float pred_value = 0, pred_err;
84        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
85        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
86        size_t totalSampleSize = (r1-1)*(r2-1)/confparams_cpr->sampleDistance;
87
88        //float max = oriData[0];
89        //float min = oriData[0];
90
91        for(i=1;i<r1;i++)
92        {
93                for(j=1;j<r2;j++)
94                {
95                        if((i+j)%confparams_cpr->sampleDistance==0)
96                        {
97                                index = i*r2+j;
98                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
99                                pred_err = fabs(pred_value - oriData[index]);
100                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
101                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
102                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
103                                intervals[radiusIndex]++;
104
105                        //      if (max < oriData[index]) max = oriData[index];
106                        //      if (min > oriData[index]) min = oriData[index];
107                        }                       
108                }
109        }
110        //compute the appropriate number
111        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
112        size_t sum = 0;
113        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
114        {
115                sum += intervals[i];
116                if(sum>targetCount)
117                        break;
118        }
119        if(i>=confparams_cpr->maxRangeRadius)
120                i = confparams_cpr->maxRangeRadius-1;
121        unsigned int accIntervals = 2*(i+1);
122        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
123
124        if(powerOf2<32)
125                powerOf2 = 32;
126
127        //      struct timeval costStart, costEnd;
128        //      double cost_est = 0;
129        //
130        //      gettimeofday(&costStart, NULL);
131        //
132        //      //compute estimate of bit-rate and distortion
133        //      double est_br = 0;
134        //      double est_psnr = 0;
135        //      double c1 = log2(targetCount)+1;
136        //      double c2 = -20.0*log10(realPrecision) + 20.0*log10(max-min) + 10.0*log10(3);
137        //
138        //      for (i = 0; i < powerOf2/2; i++)
139        //      {
140        //              int count = intervals[i];
141        //              if (count != 0)
142        //                      est_br += count*log2(count);
143        //              est_psnr += count;
144        //      }
145        //
146        //      //compute estimate of bit-rate
147        //      est_br -= c1*est_psnr;
148        //      est_br /= totalSampleSize;
149        //      est_br = -est_br;
150        //
151        //      //compute estimate of psnr
152        //      est_psnr /= totalSampleSize;
153        //      printf ("sum of P(i) = %lf\n", est_psnr);
154        //      est_psnr = -10.0*log10(est_psnr);
155        //      est_psnr += c2;
156        //
157        //      printf ("estimate bitrate = %.2f\n", est_br);
158        //      printf ("estimate psnr = %.2f\n",est_psnr);
159        //
160        //      gettimeofday(&costEnd, NULL);
161        //      cost_est = ((costEnd.tv_sec*1000000+costEnd.tv_usec)-(costStart.tv_sec*1000000+costStart.tv_usec))/1000000.0;
162        //
163        //      printf ("analysis time = %f\n", cost_est);
164
165        free(intervals);
166        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
167        return powerOf2;
168}
169
170unsigned int optimize_intervals_float_3D(float *oriData, size_t r1, size_t r2, size_t r3, double realPrecision)
171{       
172        size_t i,j,k, index;
173        size_t radiusIndex;
174        size_t r23=r2*r3;
175        float pred_value = 0, pred_err;
176        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
177        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
178        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)/confparams_cpr->sampleDistance;
179
180        //float max = oriData[0];
181        //float min = oriData[0];
182
183        for(i=1;i<r1;i++)
184        {
185                for(j=1;j<r2;j++)
186                {
187                        for(k=1;k<r3;k++)
188                        {                       
189                                if((i+j+k)%confparams_cpr->sampleDistance==0)
190                                {
191                                        index = i*r23+j*r3+k;
192                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23] 
193                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
194                                        pred_err = fabs(pred_value - oriData[index]);
195                                        radiusIndex = (pred_err/realPrecision+1)/2;
196                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
197                                        {
198                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
199                                                //printf("radiusIndex=%d\n", radiusIndex);
200                                        }
201                                        intervals[radiusIndex]++;
202
203                                        //      if (max < oriData[index]) max = oriData[index];
204                                        //      if (min > oriData[index]) min = oriData[index];
205                                }
206                        }
207                }
208        }
209        //compute the appropriate number
210        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
211        size_t sum = 0;
212        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
213        {
214                sum += intervals[i];
215                if(sum>targetCount)
216                        break;
217        }
218        if(i>=confparams_cpr->maxRangeRadius)
219                i = confparams_cpr->maxRangeRadius-1;
220        unsigned int accIntervals = 2*(i+1);
221        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
222
223        if(powerOf2<32)
224                powerOf2 = 32;
225       
226        //      struct timeval costStart, costEnd;
227        //      double cost_est = 0;
228        //
229        //      gettimeofday(&costStart, NULL);
230        //
231        //      //compute estimate of bit-rate and distortion
232        //      double est_br = 0;
233        //      double est_psnr = 0;
234        //      double c1 = log2(targetCount)+1;
235        //      double c2 = -20.0*log10(realPrecision) + 20.0*log10(max-min) + 10.0*log10(3);
236        //
237        //      for (i = 0; i < powerOf2/2; i++)
238        //      {
239        //              int count = intervals[i];
240        //              if (count != 0)
241        //                      est_br += count*log2(count);
242        //              est_psnr += count;
243        //      }
244        //
245        //      //compute estimate of bit-rate
246        //      est_br -= c1*est_psnr;
247        //      est_br /= totalSampleSize;
248        //      est_br = -est_br;
249        //
250        //      //compute estimate of psnr
251        //      est_psnr /= totalSampleSize;
252        //      printf ("sum of P(i) = %lf\n", est_psnr);
253        //      est_psnr = -10.0*log10(est_psnr);
254        //      est_psnr += c2;
255        //
256        //      printf ("estimate bitrate = %.2f\n", est_br);
257        //      printf ("estimate psnr = %.2f\n",est_psnr);
258        //
259        //      gettimeofday(&costEnd, NULL);
260        //      cost_est = ((costEnd.tv_sec*1000000+costEnd.tv_usec)-(costStart.tv_sec*1000000+costStart.tv_usec))/1000000.0;
261        //
262        //      printf ("analysis time = %f\n", cost_est);
263
264        free(intervals);
265        //printf("targetCount=%d, sum=%d, totalSampleSize=%d, ratio=%f, accIntervals=%d, powerOf2=%d\n", targetCount, sum, totalSampleSize, (double)sum/(double)totalSampleSize, accIntervals, powerOf2);
266        return powerOf2;
267}
268
269
270unsigned int optimize_intervals_float_4D(float *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision)
271{
272        size_t i,j,k,l, index;
273        size_t radiusIndex;
274        size_t r234=r2*r3*r4;
275        size_t r34=r3*r4;
276        float pred_value = 0, pred_err;
277        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
278        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
279        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)*(r4-1)/confparams_cpr->sampleDistance;
280        for(i=1;i<r1;i++)
281        {
282                for(j=1;j<r2;j++)
283                {
284                        for(k=1;k<r3;k++)
285                        {
286                                for (l=1;l<r4;l++)
287                                {
288                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
289                                        {
290                                                index = i*r234+j*r34+k*r4+l;
291                                                pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r34]
292                                                                - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
293                                                pred_err = fabs(pred_value - oriData[index]);
294                                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
295                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
296                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
297                                                intervals[radiusIndex]++;
298                                        }
299                                }
300                        }
301                }
302        }
303        //compute the appropriate number
304        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
305        size_t sum = 0;
306        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
307        {
308                sum += intervals[i];
309                if(sum>targetCount)
310                        break;
311        }
312        if(i>=confparams_cpr->maxRangeRadius)
313                i = confparams_cpr->maxRangeRadius-1;
314
315        unsigned int accIntervals = 2*(i+1);
316        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
317
318        if(powerOf2<32)
319                powerOf2 = 32;
320
321        free(intervals);
322        return powerOf2;
323}
324
325TightDataPointStorageF* SZ_compress_float_1D_MDQ(float *oriData, 
326size_t dataLength, double realPrecision, float valueRangeSize, float medianValue_f)
327{
328#ifdef HAVE_TIMECMPR   
329        float* decData = NULL;
330        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
331                decData = (float*)(multisteps->hist_data);
332#endif 
333       
334        unsigned int quantization_intervals;
335        if(exe_params->optQuantMode==1)
336                quantization_intervals = optimize_intervals_float_1D_opt(oriData, dataLength, realPrecision);
337        else
338                quantization_intervals = exe_params->intvCapacity;
339        updateQuantizationInfo(quantization_intervals); 
340
341        size_t i;
342        int reqLength;
343        float medianValue = medianValue_f;
344        short radExpo = getExponent_float(valueRangeSize/2);
345       
346        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);       
347
348        int* type = (int*) malloc(dataLength*sizeof(int));
349               
350        float* spaceFillingValue = oriData; //
351       
352        DynamicIntArray *exactLeadNumArray;
353        new_DIA(&exactLeadNumArray, DynArrayInitLen);
354       
355        DynamicByteArray *exactMidByteArray;
356        new_DBA(&exactMidByteArray, DynArrayInitLen);
357       
358        DynamicIntArray *resiBitArray;
359        new_DIA(&resiBitArray, DynArrayInitLen);
360       
361        unsigned char preDataBytes[4];
362        intToBytes_bigEndian(preDataBytes, 0);
363       
364        int reqBytesLength = reqLength/8;
365        int resiBitsLength = reqLength%8;
366        float last3CmprsData[3] = {0};
367
368        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
369        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
370                               
371        //add the first data   
372        type[0] = 0;
373        compressSingleFloatValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
374        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
375        memcpy(preDataBytes,vce->curBytes,4);
376        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
377        listAdd_float(last3CmprsData, vce->data);
378#ifdef HAVE_TIMECMPR   
379        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
380                decData[0] = vce->data;
381#endif         
382               
383        //add the second data
384        type[1] = 0;
385        compressSingleFloatValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
386        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
387        memcpy(preDataBytes,vce->curBytes,4);
388        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
389        listAdd_float(last3CmprsData, vce->data);
390#ifdef HAVE_TIMECMPR   
391        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
392                decData[1] = vce->data;
393#endif
394        int state;
395        double checkRadius;
396        float curData;
397        float pred;
398        float predAbsErr;
399        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
400        double interval = 2*realPrecision;
401       
402        for(i=2;i<dataLength;i++)
403        {       
404                curData = spaceFillingValue[i];
405                //pred = 2*last3CmprsData[0] - last3CmprsData[1];
406                pred = last3CmprsData[0];
407                predAbsErr = fabs(curData - pred);     
408                if(predAbsErr<=checkRadius)
409                {
410                        state = (predAbsErr/realPrecision+1)/2;
411                        if(curData>=pred)
412                        {
413                                type[i] = exe_params->intvRadius+state;
414                                pred = pred + state*interval;
415                        }
416                        else //curData<pred
417                        {
418                                type[i] = exe_params->intvRadius-state;
419                                pred = pred - state*interval;
420                        }
421                               
422                        //double-check the prediction error in case of machine-epsilon impact   
423                        if(fabs(curData-pred)>realPrecision)
424                        {       
425                                type[i] = 0;                           
426                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
427                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
428                                memcpy(preDataBytes,vce->curBytes,4);
429                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);         
430                               
431                                listAdd_float(last3CmprsData, vce->data);       
432#ifdef HAVE_TIMECMPR                                   
433                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
434                                        decData[i] = vce->data;
435#endif                                 
436                        }
437                        else
438                        {
439                                listAdd_float(last3CmprsData, pred);
440#ifdef HAVE_TIMECMPR                                   
441                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
442                                        decData[i] = pred;                     
443#endif 
444                        }       
445                        continue;
446                }
447               
448                //unpredictable data processing         
449                type[i] = 0;           
450                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
451                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
452                memcpy(preDataBytes,vce->curBytes,4);
453                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
454
455                listAdd_float(last3CmprsData, vce->data);
456#ifdef HAVE_TIMECMPR
457                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
458                        decData[i] = vce->data;
459#endif 
460               
461        }//end of for
462               
463//      char* expSegmentsInBytes;
464//      int expSegmentsInBytes_size = convertESCToBytes(esc, &expSegmentsInBytes);
465        size_t exactDataNum = exactLeadNumArray->size;
466       
467        TightDataPointStorageF* tdps;
468                       
469        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum, 
470                        type, exactMidByteArray->array, exactMidByteArray->size, 
471                        exactLeadNumArray->array, 
472                        resiBitArray->array, resiBitArray->size, 
473                        resiBitsLength,
474                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
475
476//sdi:Debug
477/*      int sum =0;
478        for(i=0;i<dataLength;i++)
479                if(type[i]==0) sum++;
480        printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);*/
481       
482        //free memory
483        free_DIA(exactLeadNumArray);
484        free_DIA(resiBitArray);
485        free(type);     
486        free(vce);
487        free(lce);     
488        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
489       
490        return tdps;
491}
492
493void SZ_compress_args_float_StoreOriData(float* oriData, size_t dataLength, TightDataPointStorageF* tdps, 
494unsigned char** newByteData, size_t *outSize)
495{
496        int floatSize=sizeof(float);   
497        size_t k = 0, i;
498        tdps->isLossless = 1;
499        size_t totalByteLength = 3 + MetaDataByteLength + exe_params->SZ_SIZE_TYPE + 1 + floatSize*dataLength;
500        *newByteData = (unsigned char*)malloc(totalByteLength);
501       
502        unsigned char dsLengthBytes[8];
503        for (i = 0; i < 3; i++)//3
504                (*newByteData)[k++] = versionNumber[i];
505
506        if(exe_params->SZ_SIZE_TYPE==4)//1
507                (*newByteData)[k++] = 16; //00010000
508        else
509                (*newByteData)[k++] = 80;       //01010000: 01000000 indicates the SZ_SIZE_TYPE=8
510       
511        convertSZParamsToBytes(confparams_cpr, &((*newByteData)[k]));
512        k = k + MetaDataByteLength;     
513       
514        sizeToBytes(dsLengthBytes,dataLength); //SZ_SIZE_TYPE: 4 or 8   
515        for (i = 0; i < exe_params->SZ_SIZE_TYPE; i++)
516                (*newByteData)[k++] = dsLengthBytes[i];
517               
518        if(sysEndianType==BIG_ENDIAN_SYSTEM)
519                memcpy((*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, oriData, dataLength*floatSize);
520        else
521        {
522                unsigned char* p = (*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
523                for(i=0;i<dataLength;i++,p+=floatSize)
524                        floatToBytes(p, oriData[i]);
525        }       
526        *outSize = totalByteLength;
527}
528
529char SZ_compress_args_float_NoCkRngeNoGzip_1D(unsigned char** newByteData, float *oriData, 
530size_t dataLength, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f)
531{               
532        char compressionType = 0;       
533        TightDataPointStorageF* tdps = NULL;   
534
535#ifdef HAVE_TIMECMPR
536        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
537        {
538                int timestep = sz_tsc->currentStep;
539                if(timestep % confparams_cpr->snapshotCmprStep != 0)
540                {
541                        tdps = SZ_compress_float_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_f);
542                        compressionType = 1; //time-series based compression
543                }
544                else
545                {       
546                        tdps = SZ_compress_float_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_f);
547                        compressionType = 0; //snapshot-based compression
548                        multisteps->lastSnapshotStep = timestep;
549                }               
550        }
551        else
552#endif
553                tdps = SZ_compress_float_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_f);     
554
555        convertTDPStoFlatBytes_float(tdps, newByteData, outSize);
556       
557        if(*outSize>dataLength*sizeof(float))
558                SZ_compress_args_float_StoreOriData(oriData, dataLength+2, tdps, newByteData, outSize);
559       
560        free_TightDataPointStorageF(tdps);
561        return compressionType;
562}
563
564TightDataPointStorageF* SZ_compress_float_2D_MDQ(float *oriData, size_t r1, size_t r2, double realPrecision, float valueRangeSize, float medianValue_f)
565{
566#ifdef HAVE_TIMECMPR
567        float* decData = NULL; 
568        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
569                decData = (float*)(multisteps->hist_data);
570#endif 
571       
572        unsigned int quantization_intervals;
573        if(exe_params->optQuantMode==1)
574        {
575                quantization_intervals = optimize_intervals_float_2D_opt(oriData, r1, r2, realPrecision);
576                updateQuantizationInfo(quantization_intervals);
577        }       
578        else
579                quantization_intervals = exe_params->intvCapacity;
580        size_t i,j; 
581        int reqLength;
582        float pred1D, pred2D;
583        float diff = 0.0;
584        double itvNum = 0;
585        float *P0, *P1;
586               
587        size_t dataLength = r1*r2;     
588       
589        P0 = (float*)malloc(r2*sizeof(float));
590        memset(P0, 0, r2*sizeof(float));
591        P1 = (float*)malloc(r2*sizeof(float));
592        memset(P1, 0, r2*sizeof(float));
593               
594        float medianValue = medianValue_f;
595        short radExpo = getExponent_float(valueRangeSize/2);
596        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);       
597
598        int* type = (int*) malloc(dataLength*sizeof(int));
599        //type[dataLength]=0;
600               
601        float* spaceFillingValue = oriData; //
602
603        DynamicIntArray *exactLeadNumArray;
604        new_DIA(&exactLeadNumArray, DynArrayInitLen);
605       
606        DynamicByteArray *exactMidByteArray;
607        new_DBA(&exactMidByteArray, DynArrayInitLen);
608       
609        DynamicIntArray *resiBitArray;
610        new_DIA(&resiBitArray, DynArrayInitLen);
611       
612        type[0] = 0;
613        unsigned char preDataBytes[4];
614        intToBytes_bigEndian(preDataBytes, 0);
615       
616        int reqBytesLength = reqLength/8;
617        int resiBitsLength = reqLength%8;
618
619        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
620        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
621                       
622        /* Process Row-0 data 0*/
623        type[0] = 0;
624        compressSingleFloatValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
625        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
626        memcpy(preDataBytes,vce->curBytes,4);
627        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
628        P1[0] = vce->data;
629#ifdef HAVE_TIMECMPR   
630        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
631                decData[0] = vce->data;
632#endif 
633
634        float curData;
635
636        /* Process Row-0 data 1*/
637        pred1D = P1[0];
638        curData = spaceFillingValue[1];
639        diff = curData - pred1D;
640
641        itvNum =  fabs(diff)/realPrecision + 1;
642
643        if (itvNum < exe_params->intvCapacity)
644        {
645                if (diff < 0) itvNum = -itvNum;
646                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
647                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;       
648
649                //ganrantee comporession error against the case of machine-epsilon
650                if(fabs(spaceFillingValue[1]-P1[1])>realPrecision)
651                {       
652                        type[1] = 0;                   
653                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
654                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
655                        memcpy(preDataBytes,vce->curBytes,4);
656                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
657                       
658                        P1[1] = vce->data;
659                }               
660        }
661        else
662        {
663                type[1] = 0;
664                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
665                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
666                memcpy(preDataBytes,vce->curBytes,4);
667                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
668                P1[1] = vce->data;
669        }
670#ifdef HAVE_TIMECMPR   
671        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
672                decData[1] = P1[1];
673#endif
674
675    /* Process Row-0 data 2 --> data r2-1 */
676        for (j = 2; j < r2; j++)
677        {
678                pred1D = 2*P1[j-1] - P1[j-2];
679                curData = spaceFillingValue[j];
680                diff = curData - pred1D;
681
682                itvNum = fabs(diff)/realPrecision + 1;
683
684                if (itvNum < exe_params->intvCapacity)
685                {
686                        if (diff < 0) itvNum = -itvNum;
687                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
688                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
689               
690                        //ganrantee comporession error against the case of machine-epsilon
691                        if(fabs(curData-P1[j])>realPrecision)
692                        {       
693                                type[j] = 0;                           
694                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
695                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
696                                memcpy(preDataBytes,vce->curBytes,4);
697                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
698                               
699                                P1[j] = vce->data;     
700                        }
701                }
702                else
703                {
704                        type[j] = 0;
705                        compressSingleFloatValue(vce,curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
706                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
707                        memcpy(preDataBytes,vce->curBytes,4);
708                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
709                        P1[j] = vce->data;
710                }
711#ifdef HAVE_TIMECMPR   
712                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
713                        decData[j] = P1[j];
714#endif         
715        }
716
717        /* Process Row-1 --> Row-r1-1 */
718        size_t index;
719        for (i = 1; i < r1; i++)
720        {       
721                /* Process row-i data 0 */
722                index = i*r2;
723                pred1D = P1[0];
724                curData = spaceFillingValue[index];
725                diff = curData - pred1D;
726
727                itvNum = fabs(diff)/realPrecision + 1;
728
729                if (itvNum < exe_params->intvCapacity)
730                {
731                        if (diff < 0) itvNum = -itvNum;
732                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
733                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
734
735                        //ganrantee comporession error against the case of machine-epsilon
736                        if(fabs(curData-P0[0])>realPrecision)
737                        {       
738                                type[index] = 0;                               
739                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
740                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
741                                memcpy(preDataBytes,vce->curBytes,4);
742                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
743                               
744                                P0[0] = vce->data;     
745                        }
746                }
747                else
748                {
749                        type[index] = 0;
750                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
751                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
752                        memcpy(preDataBytes,vce->curBytes,4);
753                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
754                        P0[0] = vce->data;
755                }
756#ifdef HAVE_TIMECMPR   
757                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
758                        decData[index] = P0[0];
759#endif
760                                                                       
761                /* Process row-i data 1 --> r2-1*/
762                for (j = 1; j < r2; j++)
763                {
764                        index = i*r2+j;
765                        pred2D = P0[j-1] + P1[j] - P1[j-1];
766
767                        curData = spaceFillingValue[index];
768                        diff = curData - pred2D;
769
770                        itvNum = fabs(diff)/realPrecision + 1;
771
772                        if (itvNum < exe_params->intvCapacity)
773                        {
774                                if (diff < 0) itvNum = -itvNum;
775                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
776                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
777                       
778                                //ganrantee comporession error against the case of machine-epsilon
779                                if(fabs(curData-P0[j])>realPrecision)
780                                {       
781                                        type[index] = 0;                                       
782                                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
783                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
784                                        memcpy(preDataBytes,vce->curBytes,4);
785                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
786                                       
787                                        P0[j] = vce->data;     
788                                }                       
789                        }
790                        else
791                        {
792                                type[index] = 0;
793                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
794                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
795                                memcpy(preDataBytes,vce->curBytes,4);
796                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
797                                P0[j] = vce->data;
798                        }
799#ifdef HAVE_TIMECMPR   
800                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
801                                decData[index] = P0[j];
802#endif                 
803                }
804
805                float *Pt;
806                Pt = P1;
807                P1 = P0;
808                P0 = Pt;
809        }
810       
811        if(r2!=1)
812                free(P0);
813        free(P1);                       
814        size_t exactDataNum = exactLeadNumArray->size;
815       
816        TightDataPointStorageF* tdps;
817                       
818        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum, 
819                        type, exactMidByteArray->array, exactMidByteArray->size, 
820                        exactLeadNumArray->array, 
821                        resiBitArray->array, resiBitArray->size, 
822                        resiBitsLength, 
823                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
824
825//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
826//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
827       
828//      for(i = 3800;i<3844;i++)
829//              printf("exactLeadNumArray->array[%d]=%d\n",i,exactLeadNumArray->array[i]);
830       
831        //free memory
832        free_DIA(exactLeadNumArray);
833        free_DIA(resiBitArray);
834        free(type);
835        free(vce);
836        free(lce);
837        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
838       
839        return tdps;   
840}
841
842/**
843 *
844 * Note: @r1 is high dimension
845 *               @r2 is low dimension
846 * */
847char SZ_compress_args_float_NoCkRngeNoGzip_2D(unsigned char** newByteData, float *oriData, size_t r1, size_t r2, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f)
848{       
849        size_t dataLength = r1*r2;
850        char compressionType = 0;       
851        TightDataPointStorageF* tdps = NULL; 
852
853#ifdef HAVE_TIMECMPR
854        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
855        {
856                int timestep = sz_tsc->currentStep;
857                if(timestep % confparams_cpr->snapshotCmprStep != 0)
858                {
859                        tdps = SZ_compress_float_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_f);
860                        compressionType = 1; //time-series based compression
861                }
862                else
863                {       
864                        tdps = SZ_compress_float_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_f);
865                        compressionType = 0; //snapshot-based compression
866                        multisteps->lastSnapshotStep = timestep;
867                }               
868        }
869        else
870#endif
871                tdps = SZ_compress_float_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_f); 
872
873        convertTDPStoFlatBytes_float(tdps, newByteData, outSize);
874
875        if(*outSize>dataLength*sizeof(float))
876                SZ_compress_args_float_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
877       
878        free_TightDataPointStorageF(tdps);     
879       
880        return compressionType;
881}
882
883TightDataPointStorageF* SZ_compress_float_3D_MDQ(float *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, float valueRangeSize, float medianValue_f)
884{
885#ifdef HAVE_TIMECMPR   
886        float* decData = NULL;
887        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
888                decData = (float*)(multisteps->hist_data);
889#endif         
890
891        unsigned int quantization_intervals;
892        if(exe_params->optQuantMode==1)
893        {
894                quantization_intervals = optimize_intervals_float_3D_opt(oriData, r1, r2, r3, realPrecision);
895                updateQuantizationInfo(quantization_intervals);
896        }       
897        else
898                quantization_intervals = exe_params->intvCapacity;
899        size_t i,j,k; 
900        int reqLength;
901        float pred1D, pred2D, pred3D;
902        float diff = 0.0;
903        double itvNum = 0;
904        float *P0, *P1;
905
906        size_t dataLength = r1*r2*r3;
907        size_t r23 = r2*r3;
908        P0 = (float*)malloc(r23*sizeof(float));
909        P1 = (float*)malloc(r23*sizeof(float));
910
911        float medianValue = medianValue_f;
912        short radExpo = getExponent_float(valueRangeSize/2);
913        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);       
914
915        int* type = (int*) malloc(dataLength*sizeof(int));
916
917        float* spaceFillingValue = oriData; //
918
919        DynamicIntArray *exactLeadNumArray;
920        new_DIA(&exactLeadNumArray, DynArrayInitLen);
921
922        DynamicByteArray *exactMidByteArray;
923        new_DBA(&exactMidByteArray, DynArrayInitLen);
924
925        DynamicIntArray *resiBitArray;
926        new_DIA(&resiBitArray, DynArrayInitLen);
927
928        unsigned char preDataBytes[4];
929        intToBytes_bigEndian(preDataBytes, 0);
930       
931        int reqBytesLength = reqLength/8;
932        int resiBitsLength = reqLength%8;
933
934        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
935        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
936
937
938        ///////////////////////////     Process layer-0 ///////////////////////////
939        /* Process Row-0 data 0*/
940        type[0] = 0;
941        compressSingleFloatValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
942        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
943        memcpy(preDataBytes,vce->curBytes,4);
944        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
945        P1[0] = vce->data;
946#ifdef HAVE_TIMECMPR   
947                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
948                        decData[0] = P1[0];
949#endif
950
951        float curData;
952
953        /* Process Row-0 data 1*/
954        pred1D = P1[0];
955        curData = spaceFillingValue[1];
956        diff = curData - pred1D;
957
958        itvNum = fabs(diff)/realPrecision + 1;
959
960        if (itvNum < exe_params->intvCapacity)
961        {
962                if (diff < 0) itvNum = -itvNum;
963                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
964                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
965               
966                //ganrantee comporession error against the case of machine-epsilon
967                if(fabs(curData-P1[1])>realPrecision)
968                {       
969                        type[1] = 0;                   
970                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
971                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
972                        memcpy(preDataBytes,vce->curBytes,4);
973                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
974                       
975                        P1[1] = vce->data;     
976                }                               
977        }
978        else
979        {
980                type[1] = 0;
981                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
982                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
983                memcpy(preDataBytes,vce->curBytes,4);
984                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
985                P1[1] = vce->data;
986        }
987#ifdef HAVE_TIMECMPR   
988        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
989                decData[1] = P1[1];
990#endif
991
992    /* Process Row-0 data 2 --> data r3-1 */
993        for (j = 2; j < r3; j++)
994        {
995                pred1D = 2*P1[j-1] - P1[j-2];
996                curData = spaceFillingValue[j];
997                diff = curData - pred1D;
998
999                itvNum = fabs(diff)/realPrecision + 1;
1000
1001                if (itvNum < exe_params->intvCapacity)
1002                {
1003                        if (diff < 0) itvNum = -itvNum;
1004                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
1005                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
1006                       
1007                        //ganrantee comporession error against the case of machine-epsilon
1008                        if(fabs(curData-P1[j])>realPrecision)
1009                        {       
1010                                type[j] = 0;                           
1011                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1012                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1013                                memcpy(preDataBytes,vce->curBytes,4);
1014                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
1015                               
1016                                P1[j] = vce->data;     
1017                        }                       
1018                }
1019                else
1020                {
1021                        type[j] = 0;
1022                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1023                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1024                        memcpy(preDataBytes,vce->curBytes,4);
1025                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1026                        P1[j] = vce->data;
1027                }
1028#ifdef HAVE_TIMECMPR   
1029                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1030                        decData[j] = P1[j];
1031#endif         
1032        }
1033
1034        /* Process Row-1 --> Row-r2-1 */
1035        size_t index;
1036        for (i = 1; i < r2; i++)
1037        {
1038                /* Process row-i data 0 */
1039                index = i*r3;   
1040                pred1D = P1[index-r3];
1041                curData = spaceFillingValue[index];
1042                diff = curData - pred1D;
1043
1044                itvNum = fabs(diff)/realPrecision + 1;
1045
1046                if (itvNum < exe_params->intvCapacity)
1047                {
1048                        if (diff < 0) itvNum = -itvNum;
1049                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1050                        P1[index] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1051                       
1052                        //ganrantee comporession error against the case of machine-epsilon
1053                        if(fabs(curData-P1[index])>realPrecision)
1054                        {       
1055                                type[index] = 0;                               
1056                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1057                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1058                                memcpy(preDataBytes,vce->curBytes,4);
1059                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
1060                               
1061                                P1[index] = vce->data; 
1062                        }                       
1063                }
1064                else
1065                {
1066                        type[index] = 0;
1067                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1068                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1069                        memcpy(preDataBytes,vce->curBytes,4);
1070                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1071                        P1[index] = vce->data;
1072                }
1073#ifdef HAVE_TIMECMPR   
1074                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1075                        decData[index] = P1[index];
1076#endif         
1077
1078                /* Process row-i data 1 --> data r3-1*/
1079                for (j = 1; j < r3; j++)
1080                {
1081                        index = i*r3+j;
1082                        pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];
1083
1084                        curData = spaceFillingValue[index];
1085                        diff = curData - pred2D;
1086
1087                        itvNum = fabs(diff)/realPrecision + 1;
1088
1089                        if (itvNum < exe_params->intvCapacity)
1090                        {
1091                                if (diff < 0) itvNum = -itvNum;
1092                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1093                                P1[index] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1094                               
1095                                //ganrantee comporession error against the case of machine-epsilon
1096                                if(fabs(curData-P1[index])>realPrecision)
1097                                {       
1098                                        type[index] = 0;                                       
1099                                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1100                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1101                                        memcpy(preDataBytes,vce->curBytes,4);
1102                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
1103                                       
1104                                        P1[index] = vce->data; 
1105                                }                               
1106                        }
1107                        else
1108                        {
1109                                type[index] = 0;
1110                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1111                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1112                                memcpy(preDataBytes,vce->curBytes,4);
1113                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1114                                P1[index] = vce->data;
1115                        }
1116#ifdef HAVE_TIMECMPR   
1117                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1118                                decData[index] = P1[index];
1119#endif                 
1120                }
1121        }
1122
1123
1124        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
1125
1126        for (k = 1; k < r1; k++)
1127        {
1128                /* Process Row-0 data 0*/
1129                index = k*r23;
1130                pred1D = P1[0];
1131                curData = spaceFillingValue[index];
1132                diff = curData - pred1D;
1133
1134                itvNum = fabs(diff)/realPrecision + 1;
1135
1136                if (itvNum < exe_params->intvCapacity)
1137                {
1138                        if (diff < 0) itvNum = -itvNum;
1139                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1140                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1141                       
1142                        //ganrantee comporession error against the case of machine-epsilon
1143                        if(fabs(curData-P0[0])>realPrecision)
1144                        {       
1145                                type[index] = 0;                               
1146                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1147                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1148                                memcpy(preDataBytes,vce->curBytes,4);
1149                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
1150                               
1151                                P0[0] = vce->data;     
1152                        }                       
1153                }
1154                else
1155                {
1156                        type[index] = 0;
1157                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1158                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1159                        memcpy(preDataBytes,vce->curBytes,4);
1160                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1161                        P0[0] = vce->data;
1162                }
1163#ifdef HAVE_TIMECMPR   
1164                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1165                        decData[index] = P0[0];
1166#endif
1167
1168            /* Process Row-0 data 1 --> data r3-1 */
1169                for (j = 1; j < r3; j++)
1170                {
1171                        //index = k*r2*r3+j;
1172                        index ++;
1173                        pred2D = P0[j-1] + P1[j] - P1[j-1];
1174                        curData = spaceFillingValue[index];
1175                        diff = spaceFillingValue[index] - pred2D;
1176
1177                        itvNum = fabs(diff)/realPrecision + 1;
1178
1179                        if (itvNum < exe_params->intvCapacity)
1180                        {
1181                                if (diff < 0) itvNum = -itvNum;
1182                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1183                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1184                                //ganrantee comporession error against the case of machine-epsilon
1185                                if(fabs(curData-P0[j])>realPrecision)
1186                                {       
1187                                        type[index] = 0;                                       
1188                                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1189                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1190                                        memcpy(preDataBytes,vce->curBytes,4);
1191                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
1192                                       
1193                                        P0[j] = vce->data;     
1194                                }
1195                        }
1196                        else
1197                        {
1198                                type[index] = 0;
1199                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1200                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1201                                memcpy(preDataBytes,vce->curBytes,4);
1202                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1203                                P0[j] = vce->data;
1204                        }
1205#ifdef HAVE_TIMECMPR   
1206                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1207                                decData[index] = P0[j];
1208#endif                 
1209                }
1210
1211            /* Process Row-1 --> Row-r2-1 */
1212                size_t index2D;
1213                for (i = 1; i < r2; i++)
1214                {
1215                        /* Process Row-i data 0 */
1216                        index = k*r23 + i*r3;
1217                        index2D = i*r3;         
1218                        pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
1219                        curData = spaceFillingValue[index];
1220                        diff = spaceFillingValue[index] - pred2D;
1221
1222                        itvNum = fabs(diff)/realPrecision + 1;
1223
1224                        if (itvNum < exe_params->intvCapacity)
1225                        {
1226                                if (diff < 0) itvNum = -itvNum;
1227                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1228                                P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1229                                //ganrantee comporession error against the case of machine-epsilon
1230                                if(fabs(curData-P0[index2D])>realPrecision)
1231                                {       
1232                                        type[index] = 0;                                       
1233                                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1234                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1235                                        memcpy(preDataBytes,vce->curBytes,4);
1236                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
1237                                       
1238                                        P0[index2D] = vce->data;       
1239                                }                               
1240                        }
1241                        else
1242                        {
1243                                type[index] = 0;
1244                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1245                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1246                                memcpy(preDataBytes,vce->curBytes,4);
1247                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1248                                P0[index2D] = vce->data;
1249                        }
1250#ifdef HAVE_TIMECMPR   
1251                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1252                                decData[index] = P0[index2D];
1253#endif                 
1254
1255                        /* Process Row-i data 1 --> data r3-1 */
1256                        for (j = 1; j < r3; j++)
1257                        {
1258//                              if(k==63&&i==43&&j==27)
1259//                                      printf("i=%d\n", i);
1260                                //index = k*r2*r3 + i*r3 + j;                   
1261                                index ++;
1262                                index2D = i*r3 + j;
1263                                pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
1264                                curData = spaceFillingValue[index];
1265                                diff = curData - pred3D;
1266
1267                                itvNum = fabs(diff)/realPrecision + 1;
1268
1269                                if (itvNum < exe_params->intvCapacity)
1270                                {
1271                                        if (diff < 0) itvNum = -itvNum;
1272                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1273                                        P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1274                                       
1275                                        //ganrantee comporession error against the case of machine-epsilon
1276                                        if(fabs(curData-P0[index2D])>realPrecision)
1277                                        {       
1278                                                type[index] = 0;                                               
1279                                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1280                                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1281                                                memcpy(preDataBytes,vce->curBytes,4);
1282                                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce); 
1283                                               
1284                                                P0[index2D] = vce->data;       
1285                                        }                                       
1286                                }
1287                                else
1288                                {
1289                                        type[index] = 0;
1290                                        compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1291                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1292                                        memcpy(preDataBytes,vce->curBytes,4);
1293                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1294                                        P0[index2D] = vce->data;
1295                                }
1296#ifdef HAVE_TIMECMPR   
1297                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1298                                        decData[index] = P0[index2D];
1299#endif                         
1300                        }
1301                }
1302
1303                float *Pt;
1304                Pt = P1;
1305                P1 = P0;
1306                P0 = Pt;
1307        }
1308        if(r23!=1)
1309                free(P0);
1310        free(P1);
1311        size_t exactDataNum = exactLeadNumArray->size;
1312
1313        TightDataPointStorageF* tdps;
1314
1315        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum,
1316                        type, exactMidByteArray->array, exactMidByteArray->size,
1317                        exactLeadNumArray->array,
1318                        resiBitArray->array, resiBitArray->size,
1319                        resiBitsLength, 
1320                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
1321
1322//sdi:Debug
1323/*      int sum =0;
1324        for(i=0;i<dataLength;i++)
1325                if(type[i]==0) sum++;
1326        printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);*/
1327
1328
1329//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
1330//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
1331
1332        //free memory
1333        free_DIA(exactLeadNumArray);
1334        free_DIA(resiBitArray);
1335        free(type);     
1336        free(vce);
1337        free(lce);
1338        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
1339       
1340        return tdps;   
1341}
1342
1343char SZ_compress_args_float_NoCkRngeNoGzip_3D(unsigned char** newByteData, float *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f)
1344{
1345        size_t dataLength = r1*r2*r3;
1346        char compressionType = 0;       
1347        TightDataPointStorageF* tdps = NULL; 
1348
1349#ifdef HAVE_TIMECMPR
1350        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1351        {
1352                int timestep = sz_tsc->currentStep;
1353                if(timestep % confparams_cpr->snapshotCmprStep != 0)
1354                {
1355                        tdps = SZ_compress_float_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_f);
1356                        compressionType = 1; //time-series based compression
1357                }
1358                else
1359                {       
1360                        tdps = SZ_compress_float_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_f);
1361                        compressionType = 0; //snapshot-based compression
1362                        multisteps->lastSnapshotStep = timestep;
1363                }               
1364        }
1365        else
1366#endif
1367                tdps = SZ_compress_float_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_f);
1368
1369
1370        convertTDPStoFlatBytes_float(tdps, newByteData, outSize);
1371
1372        if(*outSize>dataLength*sizeof(float))
1373                SZ_compress_args_float_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1374
1375        free_TightDataPointStorageF(tdps);
1376       
1377        return compressionType;
1378}
1379
1380
1381TightDataPointStorageF* SZ_compress_float_4D_MDQ(float *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, float valueRangeSize, float medianValue_f)
1382{       
1383        unsigned int quantization_intervals;
1384        if(exe_params->optQuantMode==1)
1385        {
1386                quantization_intervals = optimize_intervals_float_4D(oriData, r1, r2, r3, r4, realPrecision);
1387                updateQuantizationInfo(quantization_intervals);
1388        }
1389        else
1390                quantization_intervals = exe_params->intvCapacity;
1391
1392        size_t i,j,k; 
1393        int reqLength;
1394        float pred1D, pred2D, pred3D;
1395        float diff = 0.0;
1396        double itvNum = 0;
1397        float *P0, *P1;
1398
1399        size_t dataLength = r1*r2*r3*r4;
1400
1401        size_t r234 = r2*r3*r4;
1402        size_t r34 = r3*r4;
1403
1404        P0 = (float*)malloc(r34*sizeof(float));
1405        P1 = (float*)malloc(r34*sizeof(float));
1406
1407        float medianValue = medianValue_f;
1408        short radExpo = getExponent_float(valueRangeSize/2);
1409        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);
1410
1411        int* type = (int*) malloc(dataLength*sizeof(int));
1412
1413        float* spaceFillingValue = oriData; //
1414
1415        DynamicIntArray *exactLeadNumArray;
1416        new_DIA(&exactLeadNumArray, DynArrayInitLen);
1417
1418        DynamicByteArray *exactMidByteArray;
1419        new_DBA(&exactMidByteArray, DynArrayInitLen);
1420
1421        DynamicIntArray *resiBitArray;
1422        new_DIA(&resiBitArray, DynArrayInitLen);
1423
1424        unsigned char preDataBytes[4];
1425        intToBytes_bigEndian(preDataBytes, 0);
1426
1427        int reqBytesLength = reqLength/8;
1428        int resiBitsLength = reqLength%8;
1429
1430        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
1431        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
1432
1433
1434        size_t l;
1435        for (l = 0; l < r1; l++)
1436        {
1437
1438                ///////////////////////////     Process layer-0 ///////////////////////////
1439                /* Process Row-0 data 0*/
1440                size_t index = l*r234;
1441                size_t index2D = 0;
1442
1443                type[index] = 0;
1444                compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1445                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1446                memcpy(preDataBytes,vce->curBytes,4);
1447                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1448                P1[index2D] = vce->data;
1449
1450                /* Process Row-0 data 1*/
1451                index = l*r234+1;
1452                index2D = 1;
1453
1454                pred1D = P1[index2D-1];
1455                diff = spaceFillingValue[index] - pred1D;
1456
1457                itvNum = fabs(diff)/realPrecision + 1;
1458
1459                if (itvNum < exe_params->intvCapacity)
1460                {
1461                        if (diff < 0) itvNum = -itvNum;
1462                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1463                        P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1464                }
1465                else
1466                {
1467                        type[index] = 0;
1468                        compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1469                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1470                        memcpy(preDataBytes,vce->curBytes,4);
1471                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1472                        P1[index2D] = vce->data;
1473                }
1474
1475                /* Process Row-0 data 2 --> data r4-1 */
1476                for (j = 2; j < r4; j++)
1477                {
1478                        index = l*r234+j;
1479                        index2D = j;
1480
1481                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
1482                        diff = spaceFillingValue[index] - pred1D;
1483
1484                        itvNum = fabs(diff)/realPrecision + 1;
1485
1486                        if (itvNum < exe_params->intvCapacity)
1487                        {
1488                                if (diff < 0) itvNum = -itvNum;
1489                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1490                                P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1491                        }
1492                        else
1493                        {
1494                                type[index] = 0;
1495                                compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1496                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1497                                memcpy(preDataBytes,vce->curBytes,4);
1498                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1499                                P1[index2D] = vce->data;
1500                        }
1501                }
1502
1503                /* Process Row-1 --> Row-r3-1 */
1504                for (i = 1; i < r3; i++)
1505                {
1506                        /* Process row-i data 0 */
1507                        index = l*r234+i*r4;
1508                        index2D = i*r4;
1509
1510                        pred1D = P1[index2D-r4];
1511                        diff = spaceFillingValue[index] - pred1D;
1512
1513                        itvNum = fabs(diff)/realPrecision + 1;
1514
1515                        if (itvNum < exe_params->intvCapacity)
1516                        {
1517                                if (diff < 0) itvNum = -itvNum;
1518                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1519                                P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1520                        }
1521                        else
1522                        {
1523                                type[index] = 0;
1524                                compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1525                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1526                                memcpy(preDataBytes,vce->curBytes,4);
1527                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1528                                P1[index2D] = vce->data;
1529                        }
1530
1531                        /* Process row-i data 1 --> data r4-1*/
1532                        for (j = 1; j < r4; j++)
1533                        {
1534                                index = l*r234+i*r4+j;
1535                                index2D = i*r4+j;
1536
1537                                pred2D = P1[index2D-1] + P1[index2D-r4] - P1[index2D-r4-1];
1538
1539                                diff = spaceFillingValue[index] - pred2D;
1540
1541                                itvNum = fabs(diff)/realPrecision + 1;
1542
1543                                if (itvNum < exe_params->intvCapacity)
1544                                {
1545                                        if (diff < 0) itvNum = -itvNum;
1546                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1547                                        P1[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1548                                }
1549                                else
1550                                {
1551                                        type[index] = 0;
1552                                        compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1553                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1554                                        memcpy(preDataBytes,vce->curBytes,4);
1555                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1556                                        P1[index2D] = vce->data;
1557                                }
1558                        }
1559                }
1560
1561
1562                ///////////////////////////     Process layer-1 --> layer-r2-1 ///////////////////////////
1563
1564                for (k = 1; k < r2; k++)
1565                {
1566                        /* Process Row-0 data 0*/
1567                        index = l*r234+k*r34;
1568                        index2D = 0;
1569
1570                        pred1D = P1[index2D];
1571                        diff = spaceFillingValue[index] - pred1D;
1572
1573                        itvNum = fabs(diff)/realPrecision + 1;
1574
1575                        if (itvNum < exe_params->intvCapacity)
1576                        {
1577                                if (diff < 0) itvNum = -itvNum;
1578                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1579                                P0[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1580                        }
1581                        else
1582                        {
1583                                type[index] = 0;
1584                                compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1585                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1586                                memcpy(preDataBytes,vce->curBytes,4);
1587                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1588                                P0[index2D] = vce->data;
1589                        }
1590
1591                        /* Process Row-0 data 1 --> data r4-1 */
1592                        for (j = 1; j < r4; j++)
1593                        {
1594                                index = l*r234+k*r34+j;
1595                                index2D = j;
1596
1597                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
1598                                diff = spaceFillingValue[index] - pred2D;
1599
1600                                itvNum = fabs(diff)/realPrecision + 1;
1601
1602                                if (itvNum < exe_params->intvCapacity)
1603                                {
1604                                        if (diff < 0) itvNum = -itvNum;
1605                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1606                                        P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1607                                }
1608                                else
1609                                {
1610                                        type[index] = 0;
1611                                        compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1612                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1613                                        memcpy(preDataBytes,vce->curBytes,4);
1614                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1615                                        P0[index2D] = vce->data;
1616                                }
1617                        }
1618
1619                        /* Process Row-1 --> Row-r3-1 */
1620                        for (i = 1; i < r3; i++)
1621                        {
1622                                /* Process Row-i data 0 */
1623                                index = l*r234+k*r34+i*r4;
1624                                index2D = i*r4;
1625
1626                                pred2D = P0[index2D-r4] + P1[index2D] - P1[index2D-r4];
1627                                diff = spaceFillingValue[index] - pred2D;
1628
1629                                itvNum = fabs(diff)/realPrecision + 1;
1630
1631                                if (itvNum < exe_params->intvCapacity)
1632                                {
1633                                        if (diff < 0) itvNum = -itvNum;
1634                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1635                                        P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1636                                }
1637                                else
1638                                {
1639                                        type[index] = 0;
1640                                        compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1641                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1642                                        memcpy(preDataBytes,vce->curBytes,4);
1643                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1644                                        P0[index2D] = vce->data;
1645                                }
1646
1647                                /* Process Row-i data 1 --> data r4-1 */
1648                                for (j = 1; j < r4; j++)
1649                                {
1650                                        index = l*r234+k*r34+i*r4+j;
1651                                        index2D = i*r4+j;
1652
1653                                        pred3D = P0[index2D-1] + P0[index2D-r4]+ P1[index2D] - P0[index2D-r4-1] - P1[index2D-r4] - P1[index2D-1] + P1[index2D-r4-1];
1654                                        diff = spaceFillingValue[index] - pred3D;
1655
1656
1657                                        itvNum = fabs(diff)/realPrecision + 1;
1658
1659                                        if (itvNum < exe_params->intvCapacity)
1660                                        {
1661                                                if (diff < 0) itvNum = -itvNum;
1662                                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1663                                                P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1664                                        }
1665                                        else
1666                                        {
1667                                                type[index] = 0;
1668                                                compressSingleFloatValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1669                                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1670                                                memcpy(preDataBytes,vce->curBytes,4);
1671                                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1672                                                P0[index2D] = vce->data;
1673                                        }
1674                                }
1675                        }
1676
1677                        float *Pt;
1678                        Pt = P1;
1679                        P1 = P0;
1680                        P0 = Pt;
1681                }
1682        }
1683
1684        free(P0);
1685        free(P1);
1686        size_t exactDataNum = exactLeadNumArray->size;
1687
1688        TightDataPointStorageF* tdps;
1689
1690        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum,
1691                        type, exactMidByteArray->array, exactMidByteArray->size,
1692                        exactLeadNumArray->array,
1693                        resiBitArray->array, resiBitArray->size,
1694                        resiBitsLength,
1695                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
1696
1697        //free memory
1698        free_DIA(exactLeadNumArray);
1699        free_DIA(resiBitArray);
1700        free(type);
1701        free(vce);
1702        free(lce);
1703        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
1704
1705        return tdps;
1706}
1707
1708char SZ_compress_args_float_NoCkRngeNoGzip_4D(unsigned char** newByteData, float *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f)
1709{
1710        TightDataPointStorageF* tdps = SZ_compress_float_4D_MDQ(oriData, r1, r2, r3, r4, realPrecision, valueRangeSize, medianValue_f);
1711
1712        convertTDPStoFlatBytes_float(tdps, newByteData, outSize);
1713
1714        int dataLength = r1*r2*r3*r4;
1715        if(*outSize>dataLength*sizeof(float))
1716                SZ_compress_args_float_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1717
1718        free_TightDataPointStorageF(tdps);
1719       
1720        return 0;
1721}
1722
1723void SZ_compress_args_float_withinRange(unsigned char** newByteData, float *oriData, size_t dataLength, size_t *outSize)
1724{
1725        TightDataPointStorageF* tdps = (TightDataPointStorageF*) malloc(sizeof(TightDataPointStorageF));
1726        tdps->rtypeArray = NULL;
1727        tdps->typeArray = NULL; 
1728        tdps->leadNumArray = NULL;
1729        tdps->residualMidBits = NULL;
1730       
1731        tdps->allSameData = 1;
1732        tdps->dataSeriesLength = dataLength;
1733        tdps->exactMidBytes = (unsigned char*)malloc(sizeof(unsigned char)*4);
1734        tdps->pwrErrBoundBytes = NULL;
1735        tdps->isLossless = 0;
1736        float value = oriData[0];
1737        floatToBytes(tdps->exactMidBytes, value);
1738        tdps->exactMidBytes_size = 4;
1739       
1740        size_t tmpOutSize;
1741        //unsigned char *tmpByteData;
1742        convertTDPStoFlatBytes_float(tdps, newByteData, &tmpOutSize);
1743
1744        //*newByteData = (unsigned char*)malloc(sizeof(unsigned char)*12); //for floating-point data (1+3+4+4)
1745        //memcpy(*newByteData, tmpByteData, 12);
1746        *outSize = tmpOutSize; //8+SZ_SIZE_TYPE; //8==3+1+4(float_size)
1747        free_TightDataPointStorageF(tdps);     
1748}
1749
1750int SZ_compress_args_float_wRngeNoGzip(unsigned char** newByteData, float *oriData, 
1751size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1752int errBoundMode, double absErr_Bound, double relBoundRatio, double pwrErrRatio)
1753{
1754        int status = SZ_SCES;
1755        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1756        float valueRangeSize = 0, medianValue = 0;
1757       
1758        float min = computeRangeSize_float(oriData, dataLength, &valueRangeSize, &medianValue);
1759        float max = min+valueRangeSize;
1760        double realPrecision = getRealPrecision_float(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1761               
1762        if(valueRangeSize <= realPrecision)
1763        {
1764                SZ_compress_args_float_withinRange(newByteData, oriData, dataLength, outSize);
1765        }
1766        else
1767        {
1768//              SZ_compress_args_float_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize);
1769                if(r5==0&&r4==0&&r3==0&&r2==0)
1770                {
1771                        if(errBoundMode>=PW_REL)
1772                        {       
1773                                //SZ_compress_args_float_NoCkRngeNoGzip_1D_pwr(newByteData, oriData, realPrecision, r1, outSize, min, max);
1774                                SZ_compress_args_float_NoCkRngeNoGzip_1D_pwrgroup(newByteData, oriData, r1, absErr_Bound, relBoundRatio, pwrErrRatio, valueRangeSize, medianValue, outSize);
1775                        }
1776                        else
1777                                SZ_compress_args_float_NoCkRngeNoGzip_1D(newByteData, oriData, r1, realPrecision, outSize, valueRangeSize, medianValue);
1778                }
1779                else if(r5==0&&r4==0&&r3==0)
1780                {
1781                        if(errBoundMode>=PW_REL)
1782                                SZ_compress_args_float_NoCkRngeNoGzip_2D_pwr(newByteData, oriData, realPrecision, r2, r1, outSize, min, max);
1783                        else
1784                                SZ_compress_args_float_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1785                }
1786                else if(r5==0&&r4==0)
1787                {
1788                        if(errBoundMode>=PW_REL)
1789                                SZ_compress_args_float_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r3, r2, r1, outSize, min, max);
1790                        else
1791                                SZ_compress_args_float_NoCkRngeNoGzip_3D(newByteData, oriData, r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1792                }
1793                else if(r5==0)
1794                {
1795                        if(errBoundMode>=PW_REL)
1796                                SZ_compress_args_float_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r4*r3, r2, r1, outSize, min, max);
1797                        else
1798                                SZ_compress_args_float_NoCkRngeNoGzip_3D(newByteData, oriData, r4*r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
1799                }
1800        }
1801        return status;
1802}
1803
1804int SZ_compress_args_float(unsigned char** newByteData, float *oriData, 
1805size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1806int errBoundMode, double absErr_Bound, double relBoundRatio, double pwRelBoundRatio)
1807{
1808        confparams_cpr->errorBoundMode = errBoundMode;
1809        if(errBoundMode==PW_REL)
1810        {
1811                confparams_cpr->pw_relBoundRatio = pwRelBoundRatio;     
1812                //confparams_cpr->pwr_type = SZ_PWR_MIN_TYPE;
1813                if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE && r3 != 0 )
1814                {
1815                        printf("Error: Current version doesn't support 3D data compression with point-wise relative error bound being based on pwrType=AVG\n");
1816                        exit(0);
1817                        return SZ_NSCS;
1818                }
1819        }
1820        int status = SZ_SCES;
1821        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1822       
1823        if(dataLength <= MIN_NUM_OF_ELEMENTS)
1824        {
1825                *newByteData = SZ_skip_compress_float(oriData, dataLength, outSize);
1826                return status;
1827        }
1828       
1829        float valueRangeSize = 0, medianValue = 0;
1830       
1831        float min = computeRangeSize_float(oriData, dataLength, &valueRangeSize, &medianValue);
1832        float max = min+valueRangeSize;
1833        double realPrecision = 0; 
1834       
1835        if(confparams_cpr->errorBoundMode==PSNR)
1836        {
1837                confparams_cpr->errorBoundMode = ABS;
1838                realPrecision = confparams_cpr->absErrBound = computeABSErrBoundFromPSNR(confparams_cpr->psnr, (double)confparams_cpr->predThreshold, (double)valueRangeSize);
1839                //printf("realPrecision=%lf\n", realPrecision);
1840        }
1841        else
1842                realPrecision = getRealPrecision_float(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1843               
1844        if(valueRangeSize <= realPrecision)
1845        {
1846                SZ_compress_args_float_withinRange(newByteData, oriData, dataLength, outSize);
1847        }
1848        else
1849        {
1850                size_t tmpOutSize = 0;
1851                unsigned char* tmpByteData;
1852               
1853                if (r2==0)
1854                {
1855                        if(confparams_cpr->errorBoundMode>=PW_REL)
1856                        {
1857                                //SZ_compress_args_float_NoCkRngeNoGzip_1D_pwr(&tmpByteData, oriData, realPrecision, r1, &tmpOutSize, min, max);
1858                                SZ_compress_args_float_NoCkRngeNoGzip_1D_pwrgroup(&tmpByteData, oriData, r1, absErr_Bound, relBoundRatio, pwRelBoundRatio, 
1859                                valueRangeSize, medianValue, &tmpOutSize);
1860                        }
1861                        else
1862#ifdef HAVE_TIMECMPR
1863                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
1864                                        multisteps->compressionType = SZ_compress_args_float_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1865                                else
1866#endif                         
1867                                        SZ_compress_args_float_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1868                }
1869                else
1870                if (r3==0)
1871                {                       
1872                        if(confparams_cpr->errorBoundMode>=PW_REL)
1873                                SZ_compress_args_float_NoCkRngeNoGzip_2D_pwr(&tmpByteData, oriData, realPrecision, r2, r1, &tmpOutSize, min, max);
1874                        else
1875#ifdef HAVE_TIMECMPR
1876                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                           
1877                                        multisteps->compressionType = SZ_compress_args_float_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1878                                else
1879#endif
1880                                        SZ_compress_args_float_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1881                }
1882                else
1883                if (r4==0)
1884                {
1885                        if(confparams_cpr->errorBoundMode>=PW_REL)
1886                                SZ_compress_args_float_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r3, r2, r1, &tmpOutSize, min, max);
1887                        else
1888#ifdef HAVE_TIMECMPR
1889                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                           
1890                                        multisteps->compressionType = SZ_compress_args_float_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1891                                else
1892#endif
1893                                        SZ_compress_args_float_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1894                }
1895                else
1896                if (r5==0)
1897                {
1898                        if(confparams_cpr->errorBoundMode>=PW_REL)             
1899                                SZ_compress_args_float_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r4*r3, r2, r1, &tmpOutSize, min, max);
1900                                //ToDO
1901                                //SZ_compress_args_float_NoCkRngeNoGzip_4D_pwr(&tmpByteData, oriData, r4, r3, r2, r1, &tmpOutSize, min, max);
1902                        else
1903#ifdef HAVE_TIMECMPR
1904                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                           
1905                                        multisteps->compressionType = SZ_compress_args_float_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1906                                else
1907#endif
1908                                        SZ_compress_args_float_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
1909                }
1910                else
1911                {
1912                        printf("Error: doesn't support 5 dimensions for now.\n");
1913                        status = SZ_DERR; //dimension error
1914                }
1915                //Call Gzip to do the further compression.
1916                if(confparams_cpr->szMode==SZ_BEST_SPEED)
1917                {
1918                        *outSize = tmpOutSize;
1919                        *newByteData = tmpByteData;
1920                }
1921                else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION || confparams_cpr->szMode==SZ_TEMPORAL_COMPRESSION)
1922                {
1923                        *outSize = zlib_compress5(tmpByteData, tmpOutSize, newByteData, confparams_cpr->gzipMode);
1924                        free(tmpByteData);
1925                }
1926                else
1927                {
1928                        printf("Error: Wrong setting of confparams_cpr->szMode in the float compression.\n");
1929                        status = SZ_MERR; //mode error                 
1930                }
1931        }
1932       
1933        return status;
1934}
1935
1936
1937void computeReqLength_float(double realPrecision, short radExpo, int* reqLength, float* medianValue)
1938{
1939        short reqExpo = getPrecisionReqLength_double(realPrecision);
1940        *reqLength = 9+radExpo - reqExpo; //radExpo-reqExpo == reqMantiLength
1941        if(*reqLength<9)
1942                *reqLength = 9;
1943        if(*reqLength>32)
1944        {       
1945                *reqLength = 32;
1946                *medianValue = 0;
1947        }                       
1948}
1949
1950//TODO
1951int SZ_compress_args_float_subblock(unsigned char* compressedBytes, float *oriData,
1952size_t r5, size_t r4, size_t r3, size_t r2, size_t r1,
1953size_t s5, size_t s4, size_t s3, size_t s2, size_t s1,
1954size_t e5, size_t e4, size_t e3, size_t e2, size_t e1,
1955size_t *outSize, int errBoundMode, double absErr_Bound, double relBoundRatio)
1956{
1957        int status = SZ_SCES;
1958        float valueRangeSize = 0, medianValue = 0;
1959        computeRangeSize_float_subblock(oriData, &valueRangeSize, &medianValue, r5, r4, r3, r2, r1, s5, s4, s3, s2, s1, e5, e4, e3, e2, e1);
1960
1961        double realPrecision = getRealPrecision_float(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1962
1963        if(valueRangeSize <= realPrecision)
1964        {
1965                //TODO
1966                //SZ_compress_args_float_withinRange_subblock();
1967        }
1968        else
1969        {
1970                if (r2==0)
1971                {
1972                        if(errBoundMode>=PW_REL)
1973                        {
1974                                //TODO
1975                                //SZ_compress_args_float_NoCkRngeNoGzip_1D_pwr_subblock();
1976                                printf ("Current subblock version does not support point-wise relative error bound.\n");
1977                        }
1978                        else
1979                                SZ_compress_args_float_NoCkRnge_1D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r1, s1, e1);
1980                }
1981                else
1982                if (r3==0)
1983                {
1984                        //TODO
1985                        if(errBoundMode>=PW_REL)
1986                        {
1987                                //TODO
1988                                //SZ_compress_args_float_NoCkRngeNoGzip_2D_pwr_subblock();
1989                                printf ("Current subblock version does not support point-wise relative error bound.\n");
1990                        }
1991                        else
1992                                SZ_compress_args_float_NoCkRnge_2D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r2, r1, s2, s1, e2, e1);
1993                }
1994                else
1995                if (r4==0)
1996                {
1997                        if(errBoundMode>=PW_REL)
1998                        {
1999                                //TODO
2000                                //SZ_compress_args_float_NoCkRngeNoGzip_3D_pwr_subblock();
2001                                printf ("Current subblock version does not support point-wise relative error bound.\n");
2002                        }
2003                        else
2004                                SZ_compress_args_float_NoCkRnge_3D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r3, r2, r1, s3, s2, s1, e3, e2, e1);
2005                }
2006                else
2007                if (r5==0)
2008                {
2009                        if(errBoundMode>=PW_REL)
2010                        {
2011                                //TODO
2012                                //SZ_compress_args_float_NoCkRngeNoGzip_4D_pwr_subblock();
2013                                printf ("Current subblock version does not support point-wise relative error bound.\n");
2014                        }
2015                        else
2016                                SZ_compress_args_float_NoCkRnge_4D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);
2017                }
2018                else
2019                {
2020                        printf("Error: doesn't support 5 dimensions for now.\n");
2021                        status = SZ_DERR; //dimension error
2022                }
2023        }
2024        return status;
2025}
2026
2027void SZ_compress_args_float_NoCkRnge_1D_subblock(unsigned char* compressedBytes, float *oriData, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f,
2028size_t r1, size_t s1, size_t e1)
2029{
2030        TightDataPointStorageF* tdps = SZ_compress_float_1D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_f, r1, s1, e1);
2031
2032        if (confparams_cpr->szMode==SZ_BEST_SPEED)
2033                convertTDPStoFlatBytes_float_args(tdps, compressedBytes, outSize);
2034        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
2035        {
2036                unsigned char *tmpCompBytes;
2037                size_t tmpOutSize;
2038                convertTDPStoFlatBytes_float(tdps, &tmpCompBytes, &tmpOutSize);
2039                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
2040                free(tmpCompBytes);
2041        }
2042        else
2043        {
2044                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
2045        }
2046
2047        //TODO
2048//      if(*outSize>dataLength*sizeof(float))
2049//              SZ_compress_args_float_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
2050
2051        free_TightDataPointStorageF(tdps);
2052}
2053
2054void SZ_compress_args_float_NoCkRnge_2D_subblock(unsigned char* compressedBytes, float *oriData, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f,
2055size_t r2, size_t r1, size_t s2, size_t s1, size_t e2, size_t e1)
2056{
2057        TightDataPointStorageF* tdps = SZ_compress_float_2D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_f, r2, r1, s2, s1, e2, e1);
2058
2059        if (confparams_cpr->szMode==SZ_BEST_SPEED)
2060                convertTDPStoFlatBytes_float_args(tdps, compressedBytes, outSize);
2061        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
2062        {
2063                unsigned char *tmpCompBytes;
2064                size_t tmpOutSize;
2065                convertTDPStoFlatBytes_float(tdps, &tmpCompBytes, &tmpOutSize);
2066                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
2067                free(tmpCompBytes);
2068        }
2069        else
2070        {
2071                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
2072        }
2073
2074        //TODO
2075//      if(*outSize>dataLength*sizeof(float))
2076//              SZ_compress_args_float_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
2077
2078        free_TightDataPointStorageF(tdps);
2079}
2080
2081void SZ_compress_args_float_NoCkRnge_3D_subblock(unsigned char* compressedBytes, float *oriData, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f,
2082size_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)
2083{
2084        TightDataPointStorageF* tdps = SZ_compress_float_3D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_f, r3, r2, r1, s3, s2, s1, e3, e2, e1);
2085
2086        if (confparams_cpr->szMode==SZ_BEST_SPEED)
2087                convertTDPStoFlatBytes_float_args(tdps, compressedBytes, outSize);
2088        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
2089        {
2090                unsigned char *tmpCompBytes;
2091                size_t tmpOutSize;
2092                convertTDPStoFlatBytes_float(tdps, &tmpCompBytes, &tmpOutSize);
2093                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
2094                free(tmpCompBytes);
2095        }
2096        else
2097        {
2098                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
2099        }
2100
2101        //TODO
2102//      if(*outSize>dataLength*sizeof(float))
2103//              SZ_compress_args_float_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
2104
2105        free_TightDataPointStorageF(tdps);
2106}
2107
2108void SZ_compress_args_float_NoCkRnge_4D_subblock(unsigned char* compressedBytes, float *oriData, double realPrecision, size_t *outSize, float valueRangeSize, float medianValue_f,
2109size_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)
2110{
2111        TightDataPointStorageF* tdps = SZ_compress_float_4D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_f, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);
2112
2113        if (confparams_cpr->szMode==SZ_BEST_SPEED)
2114                convertTDPStoFlatBytes_float_args(tdps, compressedBytes, outSize);
2115        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
2116        {
2117                unsigned char *tmpCompBytes;
2118                size_t tmpOutSize;
2119                convertTDPStoFlatBytes_float(tdps, &tmpCompBytes, &tmpOutSize);
2120                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
2121                free(tmpCompBytes);
2122        }
2123        else
2124        {
2125                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
2126        }
2127
2128        //TODO
2129//      if(*outSize>dataLength*sizeof(float))
2130//              SZ_compress_args_float_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
2131
2132        free_TightDataPointStorageF(tdps);
2133
2134}
2135
2136unsigned int optimize_intervals_float_1D_subblock(float *oriData, double realPrecision, size_t r1, size_t s1, size_t e1)
2137{
2138        size_t dataLength = e1 - s1 + 1;
2139        oriData = oriData + s1;
2140
2141        size_t i = 0;
2142        unsigned long radiusIndex;
2143        float pred_value = 0, pred_err;
2144        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
2145        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
2146        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
2147        for(i=2;i<dataLength;i++)
2148        {
2149                if(i%confparams_cpr->sampleDistance==0)
2150                {
2151                        pred_value = 2*oriData[i-1] - oriData[i-2];
2152                        //pred_value = oriData[i-1];
2153                        pred_err = fabs(pred_value - oriData[i]);
2154                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
2155                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
2156                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
2157                        intervals[radiusIndex]++;
2158                }
2159        }
2160        //compute the appropriate number
2161        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
2162        size_t sum = 0;
2163        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2164        {
2165                sum += intervals[i];
2166                if(sum>targetCount)
2167                        break;
2168        }
2169        if(i>=confparams_cpr->maxRangeRadius)
2170                i = confparams_cpr->maxRangeRadius-1;
2171
2172        unsigned int accIntervals = 2*(i+1);
2173        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2174
2175        if(powerOf2<32)
2176                powerOf2 = 32;
2177
2178        free(intervals);
2179        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
2180        return powerOf2;
2181}
2182
2183unsigned int optimize_intervals_float_2D_subblock(float *oriData, double realPrecision, size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
2184{
2185        size_t R1 = e1 - s1 + 1;
2186        size_t R2 = e2 - s2 + 1;
2187
2188        size_t i,j, index;
2189        unsigned long radiusIndex;
2190        float pred_value = 0, pred_err;
2191        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
2192        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
2193        size_t totalSampleSize = R1*R2/confparams_cpr->sampleDistance;
2194        for(i=s1+1;i<=e1;i++)
2195        {
2196                for(j=s2+1;j<=e2;j++)
2197                {
2198                        if((i+j)%confparams_cpr->sampleDistance==0)
2199                        {
2200                                index = i*r2+j;
2201                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
2202                                pred_err = fabs(pred_value - oriData[index]);
2203                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
2204                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
2205                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
2206                                intervals[radiusIndex]++;
2207                        }
2208                }
2209        }
2210        //compute the appropriate number
2211        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
2212        size_t sum = 0;
2213        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2214        {
2215                sum += intervals[i];
2216                if(sum>targetCount)
2217                        break;
2218        }
2219        if(i>=confparams_cpr->maxRangeRadius)
2220                i = confparams_cpr->maxRangeRadius-1;
2221        unsigned int accIntervals = 2*(i+1);
2222        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2223
2224        if(powerOf2<32)
2225                powerOf2 = 32;
2226
2227        free(intervals);
2228        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
2229        return powerOf2;
2230}
2231
2232unsigned int optimize_intervals_float_3D_subblock(float *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)
2233{
2234        size_t R1 = e1 - s1 + 1;
2235        size_t R2 = e2 - s2 + 1;
2236        size_t R3 = e3 - s3 + 1;
2237
2238        size_t r23 = r2*r3;
2239
2240        size_t i,j,k, index;
2241        unsigned long radiusIndex;
2242        float pred_value = 0, pred_err;
2243        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
2244        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
2245        size_t totalSampleSize = R1*R2*R3/confparams_cpr->sampleDistance;
2246        for(i=s1+1;i<=e1;i++)
2247        {
2248                for(j=s2+1;j<=e2;j++)
2249                {
2250                        for(k=s3+1;k<=e3;k++)
2251                        {
2252                                if((i+j+k)%confparams_cpr->sampleDistance==0)
2253                                {
2254                                        index = i*r23+j*r3+k;
2255                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23]
2256                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
2257                                        pred_err = fabs(pred_value - oriData[index]);
2258                                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
2259                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
2260                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
2261                                        intervals[radiusIndex]++;
2262                                }
2263                        }
2264                }
2265        }
2266        //compute the appropriate number
2267        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
2268        size_t sum = 0;
2269        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2270        {
2271                sum += intervals[i];
2272                if(sum>targetCount)
2273                        break;
2274        }
2275        if(i>=confparams_cpr->maxRangeRadius)
2276                i = confparams_cpr->maxRangeRadius-1;
2277        unsigned int accIntervals = 2*(i+1);
2278        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2279
2280        if(powerOf2<32)
2281                powerOf2 = 32;
2282
2283        free(intervals);
2284        return powerOf2;
2285}
2286
2287unsigned int optimize_intervals_float_4D_subblock(float *oriData, double realPrecision,
2288size_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)
2289{
2290        size_t R1 = e1 - s1 + 1;
2291        size_t R2 = e2 - s2 + 1;
2292        size_t R3 = e3 - s3 + 1;
2293        size_t R4 = e4 - s4 + 1;
2294
2295        size_t r34 = r3*r4;
2296        size_t r234 = r2*r3*r4;
2297
2298        size_t i,j,k,l, index;
2299        unsigned long radiusIndex;
2300        float pred_value = 0, pred_err;
2301        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
2302        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
2303        size_t totalSampleSize = R1*R2*R3*R4/confparams_cpr->sampleDistance;
2304        for(i=s1+1;i<=e1;i++)
2305        {
2306                for(j=s2+1;j<=e2;j++)
2307                {
2308                        for(k=s3+1;k<=e3;k++)
2309                        {
2310                                for (l=s4+1;l<=e4;l++)
2311                                {
2312                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
2313                                        {
2314                                                index = i*r234+j*r34+k*r4+l;
2315                                                pred_value = oriData[index-1] + oriData[index-r4] + oriData[index-r34]
2316                                                                        - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
2317                                                pred_err = fabs(pred_value - oriData[index]);
2318                                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
2319                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
2320                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
2321                                                intervals[radiusIndex]++;
2322                                        }
2323                                }
2324                        }
2325                }
2326        }
2327        //compute the appropriate number
2328        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
2329        size_t sum = 0;
2330        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
2331        {
2332                sum += intervals[i];
2333                if(sum>targetCount)
2334                        break;
2335        }
2336        if(i>=confparams_cpr->maxRangeRadius)
2337                i = confparams_cpr->maxRangeRadius-1;
2338
2339        unsigned int accIntervals = 2*(i+1);
2340        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
2341
2342        if(powerOf2<32)
2343                powerOf2 = 32;
2344
2345        free(intervals);
2346        return powerOf2;
2347}
2348
2349TightDataPointStorageF* SZ_compress_float_1D_MDQ_subblock(float *oriData, double realPrecision, float valueRangeSize, float medianValue_f,
2350size_t r1, size_t s1, size_t e1)
2351{
2352        size_t dataLength = e1 - s1 + 1;
2353        unsigned int quantization_intervals;
2354        if(exe_params->optQuantMode==1)
2355                quantization_intervals = optimize_intervals_float_1D_subblock(oriData, realPrecision, r1, s1, e1);
2356        else
2357                quantization_intervals = exe_params->intvCapacity;
2358        updateQuantizationInfo(quantization_intervals);
2359
2360        size_t i; 
2361        int reqLength;
2362        float medianValue = medianValue_f;
2363        short radExpo = getExponent_float(valueRangeSize/2);
2364
2365        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);
2366
2367        int* type = (int*) malloc(dataLength*sizeof(int));
2368
2369        float* spaceFillingValue = oriData + s1;
2370
2371        DynamicIntArray *exactLeadNumArray;
2372        new_DIA(&exactLeadNumArray, DynArrayInitLen);
2373
2374        DynamicByteArray *exactMidByteArray;
2375        new_DBA(&exactMidByteArray, DynArrayInitLen);
2376
2377        DynamicIntArray *resiBitArray;
2378        new_DIA(&resiBitArray, DynArrayInitLen);
2379
2380        type[0] = 0;
2381
2382        unsigned char preDataBytes[4];
2383        intToBytes_bigEndian(preDataBytes, 0);
2384
2385        int reqBytesLength = reqLength/8;
2386        int resiBitsLength = reqLength%8;
2387        float last3CmprsData[3] = {0};
2388
2389        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
2390        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
2391
2392        //add the first data
2393        compressSingleFloatValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2394        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2395        memcpy(preDataBytes,vce->curBytes,4);
2396        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2397        listAdd_float(last3CmprsData, vce->data);
2398
2399        //add the second data
2400        type[1] = 0;
2401        compressSingleFloatValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2402        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2403        memcpy(preDataBytes,vce->curBytes,4);
2404        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2405        listAdd_float(last3CmprsData, vce->data);
2406
2407        int state;
2408        double checkRadius;
2409        float curData;
2410        float pred;
2411        float predAbsErr;
2412        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
2413        double interval = 2*realPrecision;
2414
2415        for(i=2;i<dataLength;i++)
2416        {
2417                curData = spaceFillingValue[i];
2418                pred = 2*last3CmprsData[0] - last3CmprsData[1];
2419                predAbsErr = fabs(curData - pred);
2420                if(predAbsErr<=checkRadius)
2421                {
2422                        state = (predAbsErr/realPrecision+1)/2;
2423                        if(curData>=pred)
2424                        {
2425                                type[i] = exe_params->intvRadius+state;
2426                                pred = pred + state*interval;
2427                        }
2428                        else
2429                        {
2430                                type[i] = exe_params->intvRadius-state;
2431                                pred = pred - state*interval;
2432                        }
2433
2434                        listAdd_float(last3CmprsData, pred);
2435                        continue;
2436                }
2437
2438                //unpredictable data processing
2439                type[i] = 0;
2440                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2441                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2442                memcpy(preDataBytes,vce->curBytes,4);
2443                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2444
2445                listAdd_float(last3CmprsData, vce->data);
2446        }
2447
2448        size_t exactDataNum = exactLeadNumArray->size;
2449
2450        TightDataPointStorageF* tdps;
2451
2452        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum,
2453                        type, exactMidByteArray->array, exactMidByteArray->size,
2454                        exactLeadNumArray->array,
2455                        resiBitArray->array, resiBitArray->size,
2456                        resiBitsLength,
2457                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2458
2459        //free memory
2460        free_DIA(exactLeadNumArray);
2461        free_DIA(resiBitArray);
2462        free(type);
2463        free(vce);
2464        free(lce);
2465        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2466
2467        return tdps;
2468}
2469
2470TightDataPointStorageF* SZ_compress_float_2D_MDQ_subblock(float *oriData, double realPrecision, float valueRangeSize, float medianValue_f,
2471size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
2472{
2473        unsigned int quantization_intervals;
2474        if(exe_params->optQuantMode==1)
2475        {
2476                quantization_intervals = optimize_intervals_float_2D_subblock(oriData, realPrecision, r1, r2, s1, s2, e1, e2);
2477                updateQuantizationInfo(quantization_intervals);
2478        }
2479        else
2480                quantization_intervals = exe_params->intvCapacity;
2481
2482        size_t i,j; 
2483        int reqLength;
2484        float pred1D, pred2D;
2485        float diff = 0.0;
2486        double itvNum = 0;
2487        float *P0, *P1;
2488
2489        size_t R1 = e1 - s1 + 1;
2490        size_t R2 = e2 - s2 + 1;
2491        size_t dataLength = R1*R2;
2492
2493        P0 = (float*)malloc(R2*sizeof(float));
2494        memset(P0, 0, R2*sizeof(float));
2495        P1 = (float*)malloc(R2*sizeof(float));
2496        memset(P1, 0, R2*sizeof(float));
2497
2498        float medianValue = medianValue_f;
2499        short radExpo = getExponent_float(valueRangeSize/2);
2500        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);
2501
2502        int* type = (int*) malloc(dataLength*sizeof(int));
2503
2504        float* spaceFillingValue = oriData; //
2505
2506        DynamicIntArray *exactLeadNumArray;
2507        new_DIA(&exactLeadNumArray, DynArrayInitLen);
2508
2509        DynamicByteArray *exactMidByteArray;
2510        new_DBA(&exactMidByteArray, DynArrayInitLen);
2511
2512        DynamicIntArray *resiBitArray;
2513        new_DIA(&resiBitArray, DynArrayInitLen);
2514
2515        unsigned char preDataBytes[4];
2516        intToBytes_bigEndian(preDataBytes, 0);
2517
2518        int reqBytesLength = reqLength/8;
2519        int resiBitsLength = reqLength%8;
2520
2521        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
2522        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
2523
2524        /* Process Row-s1 data s2*/
2525        size_t gIndex;
2526        size_t lIndex;
2527
2528        gIndex = s1*r2+s2;
2529        lIndex = 0;
2530
2531        type[lIndex] = 0;
2532        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2533        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2534        memcpy(preDataBytes,vce->curBytes,4);
2535        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2536        P1[0] = vce->data;
2537
2538        /* Process Row-s1 data s2+1*/
2539        gIndex = s1*r2+(s2+1);
2540        lIndex = 1;
2541
2542        pred1D = P1[0];
2543        diff = spaceFillingValue[gIndex] - pred1D;
2544
2545        itvNum =  fabs(diff)/realPrecision + 1;
2546
2547        if (itvNum < exe_params->intvCapacity)
2548        {
2549                if (diff < 0) itvNum = -itvNum;
2550                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2551                P1[1] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2552        }
2553        else
2554        {
2555                type[lIndex] = 0;
2556                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2557                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2558                memcpy(preDataBytes,vce->curBytes,4);
2559                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2560                P1[1] = vce->data;
2561        }
2562
2563    /* Process Row-s1 data s2+2 --> data e2 */
2564        for (j = 2; j < R2; j++)
2565        {
2566                gIndex = s1*r2+(s2+j);
2567                lIndex = j;
2568
2569                pred1D = 2*P1[j-1] - P1[j-2];
2570                diff = spaceFillingValue[gIndex] - pred1D;
2571
2572                itvNum = fabs(diff)/realPrecision + 1;
2573
2574                if (itvNum < exe_params->intvCapacity)
2575                {
2576                        if (diff < 0) itvNum = -itvNum;
2577                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2578                        P1[j] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2579                }
2580                else
2581                {
2582                        type[lIndex] = 0;
2583                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2584                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2585                        memcpy(preDataBytes,vce->curBytes,4);
2586                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2587                        P1[j] = vce->data;
2588                }
2589        }
2590
2591        /* Process Row-s1+1 --> Row-e1 */
2592        for (i = 1; i < R1; i++)
2593        {
2594                /* Process row-s1+i data s2 */
2595                gIndex = (s1+i)*r2+s2;
2596                lIndex = i*R2;
2597
2598                pred1D = P1[0];
2599                diff = spaceFillingValue[gIndex] - pred1D;
2600
2601                itvNum = fabs(diff)/realPrecision + 1;
2602
2603                if (itvNum < exe_params->intvCapacity)
2604                {
2605                        if (diff < 0) itvNum = -itvNum;
2606                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2607                        P0[0] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2608                }
2609                else
2610                {
2611                        type[lIndex] = 0;
2612                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2613                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2614                        memcpy(preDataBytes,vce->curBytes,4);
2615                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2616                        P0[0] = vce->data;
2617                }
2618
2619                /* Process row-s1+i data s2+1 --> e2 */
2620                for (j = 1; j < R2; j++)
2621                {
2622                        gIndex = (s1+i)*r2+(s2+j);
2623                        lIndex = i*R2+j;
2624
2625//                      printf ("global index = %d, local index = %d\n", gIndex, lIndex);
2626
2627                        pred2D = P0[j-1] + P1[j] - P1[j-1];
2628
2629                        diff = spaceFillingValue[gIndex] - pred2D;
2630
2631                        itvNum = fabs(diff)/realPrecision + 1;
2632
2633                        if (itvNum < exe_params->intvCapacity)
2634                        {
2635                                if (diff < 0) itvNum = -itvNum;
2636                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2637                                P0[j] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2638                        }
2639                        else
2640                        {
2641                                type[lIndex] = 0;
2642                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2643                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2644                                memcpy(preDataBytes,vce->curBytes,4);
2645                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2646                                P0[j] = vce->data;
2647                        }
2648                }
2649
2650                float *Pt;
2651                Pt = P1;
2652                P1 = P0;
2653                P0 = Pt;
2654        }
2655
2656        free(P0);
2657        free(P1);
2658        size_t exactDataNum = exactLeadNumArray->size;
2659
2660        TightDataPointStorageF* tdps;
2661
2662        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum,
2663                        type, exactMidByteArray->array, exactMidByteArray->size,
2664                        exactLeadNumArray->array,
2665                        resiBitArray->array, resiBitArray->size,
2666                        resiBitsLength,
2667                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
2668
2669        //free memory
2670        free_DIA(exactLeadNumArray);
2671        free_DIA(resiBitArray);
2672        free(type);
2673        free(vce);
2674        free(lce);
2675        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
2676
2677        return tdps;
2678}
2679
2680TightDataPointStorageF* SZ_compress_float_3D_MDQ_subblock(float *oriData, double realPrecision, float valueRangeSize, float medianValue_f,
2681size_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)
2682{
2683        unsigned int quantization_intervals;
2684        if(exe_params->optQuantMode==1)
2685        {
2686                quantization_intervals = optimize_intervals_float_3D_subblock(oriData, realPrecision, r1, r2, r3, s1, s2, s3, e1, e2, e3);
2687                updateQuantizationInfo(quantization_intervals);
2688        }
2689        else
2690                quantization_intervals = exe_params->intvCapacity;
2691
2692        size_t i,j,k; 
2693        int reqLength;
2694        float pred1D, pred2D, pred3D;
2695        float diff = 0.0;
2696        double itvNum = 0;
2697        float *P0, *P1;
2698
2699        size_t R1 = e1 - s1 + 1;
2700        size_t R2 = e2 - s2 + 1;
2701        size_t R3 = e3 - s3 + 1;
2702        size_t dataLength = R1*R2*R3;
2703
2704        size_t r23 = r2*r3;
2705        size_t R23 = R2*R3;
2706
2707        P0 = (float*)malloc(R23*sizeof(float));
2708        P1 = (float*)malloc(R23*sizeof(float));
2709
2710        float medianValue = medianValue_f;
2711        short radExpo = getExponent_float(valueRangeSize/2);
2712        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);
2713
2714        int* type = (int*) malloc(dataLength*sizeof(int));
2715        //type[dataLength]=0;
2716
2717        float* spaceFillingValue = oriData; //
2718
2719        DynamicIntArray *exactLeadNumArray;
2720        new_DIA(&exactLeadNumArray, DynArrayInitLen);
2721
2722        DynamicByteArray *exactMidByteArray;
2723        new_DBA(&exactMidByteArray, DynArrayInitLen);
2724
2725        DynamicIntArray *resiBitArray;
2726        new_DIA(&resiBitArray, DynArrayInitLen);
2727
2728        unsigned char preDataBytes[4];
2729        intToBytes_bigEndian(preDataBytes, 0);
2730
2731        int reqBytesLength = reqLength/8;
2732        int resiBitsLength = reqLength%8;
2733
2734        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
2735        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
2736
2737
2738        ///////////////////////////     Process layer-s1 ///////////////////////////
2739        /* Process Row-s2 data s3*/
2740        size_t gIndex;  //global index
2741        size_t lIndex;  //local index
2742        size_t index2D;         //local 2D index
2743
2744        gIndex = s1*r23+s2*r3+s3;
2745        lIndex = 0;
2746        index2D = 0;
2747
2748        type[lIndex] = 0;
2749        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2750        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2751        memcpy(preDataBytes,vce->curBytes,4);
2752        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2753        P1[index2D] = vce->data;
2754
2755        /* Process Row-s2 data s3+1*/
2756        gIndex = s1*r23+s2*r3+s3+1;
2757        lIndex = 1;
2758        index2D = 1;
2759
2760        pred1D = P1[index2D-1];
2761        diff = spaceFillingValue[gIndex] - pred1D;
2762
2763        itvNum = fabs(diff)/realPrecision + 1;
2764
2765        if (itvNum < exe_params->intvCapacity)
2766        {
2767                if (diff < 0) itvNum = -itvNum;
2768                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2769                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2770        }
2771        else
2772        {
2773                type[lIndex] = 0;
2774                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2775                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2776                memcpy(preDataBytes,vce->curBytes,4);
2777                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2778                P1[index2D] = vce->data;
2779        }
2780
2781    /* Process Row-s2 data s3+2 --> data e3 */
2782        for (j = 2; j < R3; j++)
2783        {
2784                gIndex = s1*r23+s2*r3+s3+j;
2785                lIndex = j;
2786                index2D = j;
2787
2788                pred1D = 2*P1[index2D-1] - P1[index2D-2];
2789                diff = spaceFillingValue[gIndex] - pred1D;
2790
2791                itvNum = fabs(diff)/realPrecision + 1;
2792
2793                if (itvNum < exe_params->intvCapacity)
2794                {
2795                        if (diff < 0) itvNum = -itvNum;
2796                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2797                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2798                }
2799                else
2800                {
2801                        type[lIndex] = 0;
2802                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2803                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2804                        memcpy(preDataBytes,vce->curBytes,4);
2805                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2806                        P1[index2D] = vce->data;
2807                }
2808        }
2809
2810        /* Process Row-s2+1 --> Row-e2 */
2811        for (i = 1; i < R2; i++)
2812        {
2813                /* Process row-s2+i data s3 */
2814                gIndex = s1*r23+(s2+i)*r3+s3;
2815                lIndex = i*R3;
2816                index2D = i*R3;
2817
2818                pred1D  = P1[index2D-R3];
2819                diff = spaceFillingValue[gIndex] - pred1D;
2820
2821                itvNum = fabs(diff)/realPrecision + 1;
2822
2823                if (itvNum < exe_params->intvCapacity)
2824                {
2825                        if (diff < 0) itvNum = -itvNum;
2826                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2827                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2828                }
2829                else
2830                {
2831                        type[lIndex] = 0;
2832                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2833                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2834                        memcpy(preDataBytes,vce->curBytes,4);
2835                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2836                        P1[index2D] = vce->data;
2837                }
2838
2839                /* Process row-s2+i data s3+1 --> data e3*/
2840                for (j = 1; j < R3; j++)
2841                {
2842                        gIndex = s1*r23+(s2+i)*r3+s3+j;
2843                        lIndex = i*R3+j;
2844                        index2D = i*R3+j;
2845
2846                        pred2D  = P1[index2D-1] + P1[index2D-R3] - P1[index2D-R3-1];
2847                        diff = spaceFillingValue[gIndex] - pred2D;
2848
2849                        itvNum = fabs(diff)/realPrecision + 1;
2850
2851                        if (itvNum < exe_params->intvCapacity)
2852                        {
2853                                if (diff < 0) itvNum = -itvNum;
2854                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2855                                P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2856                        }
2857                        else
2858                        {
2859                                type[lIndex] = 0;
2860                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2861                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2862                                memcpy(preDataBytes,vce->curBytes,4);
2863                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2864                                P1[index2D] = vce->data;
2865                        }
2866                }
2867        }
2868
2869
2870        ///////////////////////////     Process layer-s1+1 --> layer-e1 ///////////////////////////
2871
2872        for (k = 1; k < R1; k++)
2873        {
2874                /* Process Row-s2 data s3*/
2875                gIndex = (s1+k)*r23+s2*r3+s3;
2876                lIndex = k*R23;
2877                index2D = 0;
2878
2879                pred1D = P1[index2D];
2880                diff = spaceFillingValue[gIndex] - pred1D;
2881
2882                itvNum = fabs(diff)/realPrecision + 1;
2883
2884                if (itvNum < exe_params->intvCapacity)
2885                {
2886                        if (diff < 0) itvNum = -itvNum;
2887                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2888                        P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2889                }
2890                else
2891                {
2892                        type[lIndex] = 0;
2893                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2894                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2895                        memcpy(preDataBytes,vce->curBytes,4);
2896                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2897                        P0[index2D] = vce->data;
2898                }
2899
2900            /* Process Row-s2 data s3+1 --> data e3 */
2901                for (j = 1; j < R3; j++)
2902                {
2903                        gIndex = (s1+k)*r23+s2*r3+s3+j;
2904                        lIndex = k*R23+j;
2905                        index2D = j;
2906
2907                        pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
2908                        diff = spaceFillingValue[gIndex] - pred2D;
2909
2910                        itvNum = fabs(diff)/realPrecision + 1;
2911
2912                        if (itvNum < exe_params->intvCapacity)
2913                        {
2914                                if (diff < 0) itvNum = -itvNum;
2915                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2916                                P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2917                        }
2918                        else
2919                        {
2920                                type[lIndex] = 0;
2921                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2922                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2923                                memcpy(preDataBytes,vce->curBytes,4);
2924                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2925                                P0[index2D] = vce->data;
2926                        }
2927                }
2928
2929            /* Process Row-s2+1 --> Row-e2 */
2930                for (i = 1; i < R2; i++)
2931                {
2932                        /* Process Row-s2+i data s3 */
2933                        gIndex = (s1+k)*r23+(s2+i)*r3+s3;
2934                        lIndex = k*R23+i*R3;
2935                        index2D = i*R3;
2936
2937                        pred2D = P0[index2D-R3] + P1[index2D] - P1[index2D-R3];
2938                        diff = spaceFillingValue[gIndex] - pred2D;
2939
2940                        itvNum = fabs(diff)/realPrecision + 1;
2941
2942                        if (itvNum < exe_params->intvCapacity)
2943                        {
2944                                if (diff < 0) itvNum = -itvNum;
2945                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2946                                P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2947                        }
2948                        else
2949                        {
2950                                type[lIndex] = 0;
2951                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2952                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2953                                memcpy(preDataBytes,vce->curBytes,4);
2954                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2955                                P0[index2D] = vce->data;
2956                        }
2957
2958                        /* Process Row-s2+i data s3+1 --> data e3 */
2959                        for (j = 1; j < R3; j++)
2960                        {
2961                                gIndex = (s1+k)*r23+(s2+i)*r3+s3+j;
2962                                lIndex = k*R23+i*R3+j;
2963                                index2D = i*R3+j;
2964
2965//                              printf ("global index = %d, local index = %d\n", gIndex, lIndex);
2966
2967                                pred3D = P0[index2D-1] + P0[index2D-R3]+ P1[index2D] - P0[index2D-R3-1] - P1[index2D-R3] - P1[index2D-1] + P1[index2D-R3-1];
2968                                diff = spaceFillingValue[gIndex] - pred3D;
2969
2970                                itvNum = fabs(diff)/realPrecision + 1;
2971
2972                                if (itvNum < exe_params->intvCapacity)
2973                                {
2974                                        if (diff < 0) itvNum = -itvNum;
2975                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
2976                                        P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
2977                                }
2978                                else
2979                                {
2980                                        type[lIndex] = 0;
2981                                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
2982                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
2983                                        memcpy(preDataBytes,vce->curBytes,4);
2984                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
2985                                        P0[index2D] = vce->data;
2986                                }
2987                        }
2988                }
2989
2990                float *Pt;
2991                Pt = P1;
2992                P1 = P0;
2993                P0 = Pt;
2994        }
2995
2996        free(P0);
2997        free(P1);
2998        size_t exactDataNum = exactLeadNumArray->size;
2999
3000        TightDataPointStorageF* tdps;
3001
3002        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum,
3003                        type, exactMidByteArray->array, exactMidByteArray->size,
3004                        exactLeadNumArray->array,
3005                        resiBitArray->array, resiBitArray->size,
3006                        resiBitsLength,
3007                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
3008
3009        //free memory
3010        free_DIA(exactLeadNumArray);
3011        free_DIA(resiBitArray);
3012        free(type);
3013        free(vce);
3014        free(lce);
3015        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
3016
3017        return tdps;
3018}
3019
3020TightDataPointStorageF* SZ_compress_float_4D_MDQ_subblock(float *oriData, double realPrecision, float valueRangeSize, float medianValue_f,
3021size_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)
3022{
3023        unsigned int quantization_intervals;
3024        if(exe_params->optQuantMode==1)
3025        {
3026                quantization_intervals = optimize_intervals_float_4D_subblock(oriData, realPrecision, r1, r2, r3, r4, s1, s2, s3, s4, e1, e2, e3, e4);
3027                updateQuantizationInfo(quantization_intervals);
3028        }
3029        else
3030                quantization_intervals = exe_params->intvCapacity;
3031
3032        size_t i,j,k; 
3033        int reqLength;
3034        float pred1D, pred2D, pred3D;
3035        float diff = 0.0;
3036        double itvNum = 0;
3037        float *P0, *P1;
3038
3039        size_t R1 = e1 - s1 + 1;
3040        size_t R2 = e2 - s2 + 1;
3041        size_t R3 = e3 - s3 + 1;
3042        size_t R4 = e4 - s4 + 1;
3043
3044        size_t dataLength = R1*R2*R3*R4;
3045
3046        size_t r34 = r3*r4;
3047        size_t r234 = r2*r3*r4;
3048        size_t R34 = R3*R4;
3049        size_t R234 = R2*R3*R4;
3050
3051        P0 = (float*)malloc(R34*sizeof(float));
3052        P1 = (float*)malloc(R34*sizeof(float));
3053
3054        float medianValue = medianValue_f;
3055        short radExpo = getExponent_float(valueRangeSize/2);
3056        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);
3057
3058        int* type = (int*) malloc(dataLength*sizeof(int));
3059
3060        float* spaceFillingValue = oriData; //
3061
3062        DynamicIntArray *exactLeadNumArray;
3063        new_DIA(&exactLeadNumArray, DynArrayInitLen);
3064
3065        DynamicByteArray *exactMidByteArray;
3066        new_DBA(&exactMidByteArray, DynArrayInitLen);
3067
3068        DynamicIntArray *resiBitArray;
3069        new_DIA(&resiBitArray, DynArrayInitLen);
3070
3071        unsigned char preDataBytes[4];
3072        intToBytes_bigEndian(preDataBytes, 0);
3073
3074        int reqBytesLength = reqLength/8;
3075        int resiBitsLength = reqLength%8;
3076
3077        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
3078        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
3079
3080
3081        size_t l;
3082        for (l = 0; l < R1; l++)
3083        {
3084
3085                ///////////////////////////     Process layer-s2 ///////////////////////////
3086                /* Process Row-s3 data s4*/
3087                size_t gIndex;  //global index
3088                size_t lIndex;  //local index
3089                size_t index2D;         //local 2D index
3090
3091                gIndex = (s1+l)*r234+s2*r34+s3*r4+s4;
3092                lIndex = l*R234;
3093                index2D = 0;
3094
3095                type[lIndex] = 0;
3096                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3097                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3098                memcpy(preDataBytes,vce->curBytes,4);
3099                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3100                P1[index2D] = vce->data;
3101
3102                /* Process Row-s3 data s4+1*/
3103                gIndex = (s1+l)*r234+s2*r34+s3*r4+s4+1;
3104                lIndex = l*R234+1;
3105                index2D = 1;
3106
3107                pred1D = P1[index2D-1];
3108                diff = spaceFillingValue[gIndex] - pred1D;
3109
3110                itvNum = fabs(diff)/realPrecision + 1;
3111
3112                if (itvNum < exe_params->intvCapacity)
3113                {
3114                        if (diff < 0) itvNum = -itvNum;
3115                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3116                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3117                }
3118                else
3119                {
3120                        type[lIndex] = 0;
3121                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3122                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3123                        memcpy(preDataBytes,vce->curBytes,4);
3124                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3125                        P1[index2D] = vce->data;
3126                }
3127
3128                /* Process Row-s3 data s4+2 --> data e4 */
3129                for (j = 2; j < R4; j++)
3130                {
3131                        gIndex = (s1+l)*r234+s2*r34+s3*r4+s4+j;
3132                        lIndex = l*R234+j;
3133                        index2D = j;
3134
3135                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
3136                        diff = spaceFillingValue[gIndex] - pred1D;
3137
3138                        itvNum = fabs(diff)/realPrecision + 1;
3139
3140                        if (itvNum < exe_params->intvCapacity)
3141                        {
3142                                if (diff < 0) itvNum = -itvNum;
3143                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3144                                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3145                        }
3146                        else
3147                        {
3148                                type[lIndex] = 0;
3149                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3150                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3151                                memcpy(preDataBytes,vce->curBytes,4);
3152                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3153                                P1[index2D] = vce->data;
3154                        }
3155                }
3156
3157                /* Process Row-s3+1 --> Row-e3 */
3158                for (i = 1; i < R3; i++)
3159                {
3160                        /* Process row-s2+i data s3 */
3161                        gIndex = (s1+l)*r234+s2*r34+(s3+i)*r4+s4;
3162                        lIndex = l*R234+i*R4;
3163                        index2D = i*R4;
3164
3165                        pred1D  = P1[index2D-R4];
3166                        diff = spaceFillingValue[gIndex] - pred1D;
3167
3168                        itvNum = fabs(diff)/realPrecision + 1;
3169
3170                        if (itvNum < exe_params->intvCapacity)
3171                        {
3172                                if (diff < 0) itvNum = -itvNum;
3173                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3174                                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3175                        }
3176                        else
3177                        {
3178                                type[lIndex] = 0;
3179                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3180                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3181                                memcpy(preDataBytes,vce->curBytes,4);
3182                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3183                                P1[index2D] = vce->data;
3184                        }
3185
3186                        /* Process row-s3+i data s4+1 --> data e4*/
3187                        for (j = 1; j < R4; j++)
3188                        {
3189                                gIndex = (s1+l)*r234+s2*r34+(s3+i)*r4+s4+j;
3190                                lIndex = l*R234+i*R4+j;
3191                                index2D = i*R4+j;
3192
3193                                pred2D  = P1[index2D-1] + P1[index2D-R4] - P1[index2D-R4-1];
3194                                diff = spaceFillingValue[gIndex] - pred2D;
3195
3196                                itvNum = fabs(diff)/realPrecision + 1;
3197
3198                                if (itvNum < exe_params->intvCapacity)
3199                                {
3200                                        if (diff < 0) itvNum = -itvNum;
3201                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3202                                        P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3203                                }
3204                                else
3205                                {
3206                                        type[lIndex] = 0;
3207                                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3208                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3209                                        memcpy(preDataBytes,vce->curBytes,4);
3210                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3211                                        P1[index2D] = vce->data;
3212                                }
3213                        }
3214                }
3215
3216
3217                ///////////////////////////     Process layer-s2+1 --> layer-e2 ///////////////////////////
3218
3219                for (k = 1; k < R2; k++)
3220                {
3221                        /* Process Row-s3 data s4*/
3222                        gIndex = (s1+l)*r234+(s2+k)*r34+s3*r4+s4;
3223                        lIndex = l*R234+k*R34;
3224                        index2D = 0;
3225
3226                        pred1D = P1[index2D];
3227                        diff = spaceFillingValue[gIndex] - pred1D;
3228
3229                        itvNum = fabs(diff)/realPrecision + 1;
3230
3231                        if (itvNum < exe_params->intvCapacity)
3232                        {
3233                                if (diff < 0) itvNum = -itvNum;
3234                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3235                                P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3236                        }
3237                        else
3238                        {
3239                                type[lIndex] = 0;
3240                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3241                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3242                                memcpy(preDataBytes,vce->curBytes,4);
3243                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3244                                P0[index2D] = vce->data;
3245                        }
3246
3247                        /* Process Row-s3 data s4+1 --> data e4 */
3248                        for (j = 1; j < R4; j++)
3249                        {
3250                                gIndex = (s1+l)*r234+(s2+k)*r34+s3*r4+s4+j;
3251                                lIndex = l*R234+k*R34+j;
3252                                index2D = j;
3253
3254                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
3255                                diff = spaceFillingValue[gIndex] - pred2D;
3256
3257                                itvNum = fabs(diff)/realPrecision + 1;
3258
3259                                if (itvNum < exe_params->intvCapacity)
3260                                {
3261                                        if (diff < 0) itvNum = -itvNum;
3262                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3263                                        P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3264                                }
3265                                else
3266                                {
3267                                        type[lIndex] = 0;
3268                                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3269                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3270                                        memcpy(preDataBytes,vce->curBytes,4);
3271                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3272                                        P0[index2D] = vce->data;
3273                                }
3274                        }
3275
3276                        /* Process Row-s3+1 --> Row-e3 */
3277                        for (i = 1; i < R3; i++)
3278                        {
3279                                /* Process Row-s3+i data s4 */
3280                                gIndex = (s1+l)*r234+(s2+k)*r34+(s3+i)*r4+s4;
3281                                lIndex = l*R234+k*R34+i*R4;
3282                                index2D = i*R4;
3283
3284                                pred2D = P0[index2D-R4] + P1[index2D] - P1[index2D-R4];
3285                                diff = spaceFillingValue[gIndex] - pred2D;
3286
3287                                itvNum = fabs(diff)/realPrecision + 1;
3288
3289                                if (itvNum < exe_params->intvCapacity)
3290                                {
3291                                        if (diff < 0) itvNum = -itvNum;
3292                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3293                                        P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3294                                }
3295                                else
3296                                {
3297                                        type[lIndex] = 0;
3298                                        compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3299                                        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3300                                        memcpy(preDataBytes,vce->curBytes,4);
3301                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3302                                        P0[index2D] = vce->data;
3303                                }
3304
3305                                /* Process Row-s3+i data s4+1 --> data e4 */
3306                                for (j = 1; j < R4; j++)
3307                                {
3308                                        gIndex = (s1+l)*r234+(s2+k)*r34+(s3+i)*r4+s4+j;
3309                                        lIndex = l*R234+k*R34+i*R4+j;
3310                                        index2D = i*R4+j;
3311
3312//                                      printf ("global index = %d, local index = %d\n", gIndex, lIndex);
3313
3314                                        pred3D = P0[index2D-1] + P0[index2D-R4]+ P1[index2D] - P0[index2D-R4-1] - P1[index2D-R4] - P1[index2D-1] + P1[index2D-R4-1];
3315                                        diff = spaceFillingValue[gIndex] - pred3D;
3316
3317                                        itvNum = fabs(diff)/realPrecision + 1;
3318
3319                                        if (itvNum < exe_params->intvCapacity)
3320                                        {
3321                                                if (diff < 0) itvNum = -itvNum;
3322                                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
3323                                                P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
3324                                        }
3325                                        else
3326                                        {
3327                                                type[lIndex] = 0;
3328                                                compressSingleFloatValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
3329                                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
3330                                                memcpy(preDataBytes,vce->curBytes,4);
3331                                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
3332                                                P0[index2D] = vce->data;
3333                                        }
3334                                }
3335                        }
3336
3337                        float *Pt;
3338                        Pt = P1;
3339                        P1 = P0;
3340                        P0 = Pt;
3341                }
3342
3343        }
3344
3345        free(P0);
3346        free(P1);
3347        size_t exactDataNum = exactLeadNumArray->size;
3348
3349        TightDataPointStorageF* tdps;
3350
3351        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum,
3352                        type, exactMidByteArray->array, exactMidByteArray->size,
3353                        exactLeadNumArray->array,
3354                        resiBitArray->array, resiBitArray->size,
3355                        resiBitsLength,
3356                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
3357
3358        //free memory
3359        free_DIA(exactLeadNumArray);
3360        free_DIA(resiBitArray);
3361        free(type);
3362        free(vce);
3363        free(lce);
3364        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
3365
3366        return tdps;
3367}
3368
3369unsigned int optimize_intervals_float_3D_opt(float *oriData, size_t r1, size_t r2, size_t r3, double realPrecision)
3370{       
3371        size_t i;
3372        size_t radiusIndex;
3373        size_t r23=r2*r3;
3374        float pred_value = 0, pred_err;
3375        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3376        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3377        size_t totalSampleSize = 0;//(r1-1)*(r2-1)*(r3-1)/confparams_cpr->sampleDistance;
3378
3379        size_t offset_count = confparams_cpr->sampleDistance - 2; // count r3 offset
3380        size_t offset_count_2;
3381        float * data_pos = oriData + r23 + r3 + offset_count;
3382        size_t n1_count = 1, n2_count = 1; // count i,j sum
3383        size_t len = r1 * r2 * r3;
3384        while(data_pos - oriData < len){
3385                totalSampleSize++;
3386                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];
3387                pred_err = fabs(pred_value - *data_pos);
3388                radiusIndex = (pred_err/realPrecision+1)/2;
3389                if(radiusIndex>=confparams_cpr->maxRangeRadius)
3390                {
3391                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
3392                        //printf("radiusIndex=%d\n", radiusIndex);
3393                }
3394                intervals[radiusIndex]++;
3395                // printf("TEST: %ld, i: %ld\tj: %ld\tk: %ld\n", data_pos - oriData);
3396                // fflush(stdout);
3397                offset_count += confparams_cpr->sampleDistance;
3398                if(offset_count >= r3){
3399                        n2_count ++;
3400                        if(n2_count == r2){
3401                                n1_count ++;
3402                                n2_count = 1;
3403                                data_pos += r3;
3404                        }
3405                        offset_count_2 = (n1_count + n2_count) % confparams_cpr->sampleDistance;
3406                        data_pos += (r3 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
3407                        offset_count = (confparams_cpr->sampleDistance - offset_count_2);
3408                        if(offset_count == 0) offset_count ++;
3409                }
3410                else data_pos += confparams_cpr->sampleDistance;
3411        }       
3412        // printf("sample_count: %ld\n", sample_count);
3413        // fflush(stdout);
3414        // if(*max_freq < 0.15) *max_freq *= 2;
3415        //compute the appropriate number
3416        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3417        size_t sum = 0;
3418        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3419        {
3420                sum += intervals[i];
3421                if(sum>targetCount)
3422                        break;
3423        }
3424        if(i>=confparams_cpr->maxRangeRadius)
3425                i = confparams_cpr->maxRangeRadius-1;
3426        unsigned int accIntervals = 2*(i+1);
3427        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3428
3429        if(powerOf2<32)
3430                powerOf2 = 32;
3431        free(intervals);
3432        //printf("targetCount=%d, sum=%d, totalSampleSize=%d, ratio=%f, accIntervals=%d, powerOf2=%d\n", targetCount, sum, totalSampleSize, (double)sum/(double)totalSampleSize, accIntervals, powerOf2);
3433        return powerOf2;
3434}
3435
3436size_t SZ_compress_float_3D_MDQ_RA_block(float * block_ori_data, float * mean, size_t dim_0, size_t dim_1, size_t dim_2, size_t block_dim_0, size_t block_dim_1, size_t block_dim_2, double realPrecision, float * P0, float * P1, int * type, float * unpredictable_data){
3437
3438        size_t dim0_offset = dim_1 * dim_2;
3439        size_t dim1_offset = dim_2;
3440
3441        // data_pos = block_ori_data;
3442        // for(size_t i=0; i<block_dim_0; i++){
3443        //      for(size_t j=0; j<block_dim_1; j++){
3444        //              for(size_t k=0; k<block_dim_2; k++){
3445        //                      sum += *data_pos;
3446        //                      data_pos ++;
3447        //              }
3448        //              data_pos += dim1_offset - block_dim_2;
3449        //      }
3450        //      data_pos += dim0_offset - block_dim_1 * dim1_offset;
3451        // }
3452        // size_t num_elements = block_dim_0 * block_dim_1 * block_dim_2;
3453        // if(num_elements > 0) mean[0] = sum / num_elements;
3454        // else mean[0] = 0.0;
3455        mean[0] = block_ori_data[0];
3456
3457        size_t unpredictable_count = 0;
3458        size_t r1, r2, r3;
3459        r1 = block_dim_0;
3460        r2 = block_dim_1;
3461        r3 = block_dim_2;
3462
3463        float * cur_data_pos = block_ori_data;
3464        float curData;
3465        float pred1D, pred2D, pred3D;
3466        double itvNum;
3467        double diff;
3468        size_t i, j, k;
3469        size_t r23 = r2*r3;
3470        // Process Row-0 data 0
3471        pred1D = mean[0];
3472        curData = *cur_data_pos;
3473        diff = curData - pred1D;
3474        itvNum = fabs(diff)/realPrecision + 1;
3475        if (itvNum < exe_params->intvCapacity){
3476                if (diff < 0) itvNum = -itvNum;
3477                type[0] = (int) (itvNum/2) + exe_params->intvRadius;
3478                P1[0] = pred1D + 2 * (type[0] - exe_params->intvRadius) * realPrecision;
3479                //ganrantee comporession error against the case of machine-epsilon
3480                if(fabs(curData-P1[0])>realPrecision){ 
3481                        type[0] = 0;
3482                        P1[0] = curData;
3483                        unpredictable_data[unpredictable_count ++] = curData;
3484                }               
3485        }
3486        else{
3487                type[0] = 0;
3488                P1[0] = curData;
3489                unpredictable_data[unpredictable_count ++] = curData;
3490        }
3491
3492        /* Process Row-0 data 1*/
3493        pred1D = P1[0];
3494        curData = cur_data_pos[1];
3495        diff = curData - pred1D;
3496        itvNum = fabs(diff)/realPrecision + 1;
3497        if (itvNum < exe_params->intvCapacity){
3498                if (diff < 0) itvNum = -itvNum;
3499                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
3500                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
3501                //ganrantee comporession error against the case of machine-epsilon
3502                if(fabs(curData-P1[1])>realPrecision){ 
3503                        type[1] = 0;
3504                        P1[1] = curData;       
3505                        unpredictable_data[unpredictable_count ++] = curData;
3506                }               
3507        }
3508        else{
3509                type[1] = 0;
3510                P1[1] = curData;
3511                unpredictable_data[unpredictable_count ++] = curData;
3512        }
3513    /* Process Row-0 data 2 --> data r3-1 */
3514        for (j = 2; j < r3; j++){
3515                pred1D = 2*P1[j-1] - P1[j-2];
3516                curData = cur_data_pos[j];
3517                diff = curData - pred1D;
3518                itvNum = fabs(diff)/realPrecision + 1;
3519                if (itvNum < exe_params->intvCapacity){
3520                        if (diff < 0) itvNum = -itvNum;
3521                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
3522                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
3523                        //ganrantee comporession error against the case of machine-epsilon
3524                        if(fabs(curData-P1[j])>realPrecision){ 
3525                                type[j] = 0;
3526                                P1[j] = curData;       
3527                                unpredictable_data[unpredictable_count ++] = curData;
3528                        }                       
3529                }
3530                else{
3531                        type[j] = 0;
3532                        P1[j] = curData;
3533                        unpredictable_data[unpredictable_count ++] = curData;
3534                }
3535        }
3536        cur_data_pos += dim1_offset;
3537
3538        /* Process Row-1 --> Row-r2-1 */
3539        size_t index;
3540        for (i = 1; i < r2; i++)
3541        {
3542                /* Process row-i data 0 */
3543                index = i*r3;   
3544                pred1D = P1[index-r3];
3545                curData = *cur_data_pos;
3546                diff = curData - pred1D;
3547
3548                itvNum = fabs(diff)/realPrecision + 1;
3549
3550                if (itvNum < exe_params->intvCapacity)
3551                {
3552                        if (diff < 0) itvNum = -itvNum;
3553                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
3554                        P1[index] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
3555                       
3556                        //ganrantee comporession error against the case of machine-epsilon
3557                        if(fabs(curData-P1[index])>realPrecision)
3558                        {       
3559                                type[index] = 0;
3560                                P1[index] = curData;   
3561                                unpredictable_data[unpredictable_count ++] = curData;
3562                        }                       
3563                }
3564                else
3565                {
3566                        type[index] = 0;
3567                        P1[index] = curData;
3568                        unpredictable_data[unpredictable_count ++] = curData;
3569                }
3570
3571                /* Process row-i data 1 --> data r3-1*/
3572                for (j = 1; j < r3; j++)
3573                {
3574                        index = i*r3+j;
3575                        pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];
3576
3577                        curData = cur_data_pos[j];
3578                        diff = curData - pred2D;
3579
3580                        itvNum = fabs(diff)/realPrecision + 1;
3581
3582                        if (itvNum < exe_params->intvCapacity)
3583                        {
3584                                if (diff < 0) itvNum = -itvNum;
3585                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
3586                                P1[index] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
3587                               
3588                                //ganrantee comporession error against the case of machine-epsilon
3589                                if(fabs(curData-P1[index])>realPrecision)
3590                                {       
3591                                        type[index] = 0;
3592                                        P1[index] = curData;   
3593                                        unpredictable_data[unpredictable_count ++] = curData;
3594                                }                               
3595                        }
3596                        else
3597                        {
3598                                type[index] = 0;
3599                                P1[index] = curData;
3600                                unpredictable_data[unpredictable_count ++] = curData;
3601                        }
3602                }
3603                cur_data_pos += dim1_offset;
3604        }
3605        cur_data_pos += dim0_offset - r2 * dim1_offset;
3606
3607        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
3608
3609        for (k = 1; k < r1; k++)
3610        {
3611                /* Process Row-0 data 0*/
3612                index = k*r23;
3613                pred1D = P1[0];
3614                curData = *cur_data_pos;
3615                diff = curData - pred1D;
3616                itvNum = fabs(diff)/realPrecision + 1;
3617                if (itvNum < exe_params->intvCapacity)
3618                {
3619                        if (diff < 0) itvNum = -itvNum;
3620                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
3621                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
3622                        //ganrantee comporession error against the case of machine-epsilon
3623                        if(fabs(curData-P0[0])>realPrecision)
3624                        {       
3625                                type[index] = 0;
3626                                P0[0] = curData;       
3627                                unpredictable_data[unpredictable_count ++] = curData;
3628                        }                       
3629                }
3630                else
3631                {
3632                        type[index] = 0;
3633                        P0[0] = curData;
3634                        unpredictable_data[unpredictable_count ++] = curData;
3635                }
3636            /* Process Row-0 data 1 --> data r3-1 */
3637                for (j = 1; j < r3; j++)
3638                {
3639                        //index = k*r2*r3+j;
3640                        index ++;
3641                        pred2D = P0[j-1] + P1[j] - P1[j-1];
3642                        curData = cur_data_pos[j];
3643                        diff = curData - pred2D;
3644                        itvNum = fabs(diff)/realPrecision + 1;
3645                        if (itvNum < exe_params->intvCapacity)
3646                        {
3647                                if (diff < 0) itvNum = -itvNum;
3648                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
3649                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
3650                                //ganrantee comporession error against the case of machine-epsilon
3651                                if(fabs(curData-P0[j])>realPrecision)
3652                                {       
3653                                        type[index] = 0;
3654                                        P0[j] = curData;       
3655                                        unpredictable_data[unpredictable_count ++] = curData;
3656                                }
3657                        }
3658                        else
3659                        {
3660                                type[index] = 0;
3661                                P0[j] = curData;
3662                                unpredictable_data[unpredictable_count ++] = curData;
3663                        }
3664                }
3665
3666                cur_data_pos += dim1_offset;
3667            /* Process Row-1 --> Row-r2-1 */
3668                size_t index2D;
3669                for (i = 1; i < r2; i++)
3670                {
3671                        /* Process Row-i data 0 */
3672                        index = k*r23 + i*r3;
3673                        index2D = i*r3;         
3674                        pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
3675                        curData = *cur_data_pos;
3676                        diff = curData - pred2D;
3677
3678                        itvNum = fabs(diff)/realPrecision + 1;
3679
3680                        if (itvNum < exe_params->intvCapacity)
3681                        {
3682                                if (diff < 0) itvNum = -itvNum;
3683                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
3684                                P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
3685                                //ganrantee comporession error against the case of machine-epsilon
3686                                if(fabs(curData-P0[index2D])>realPrecision)
3687                                {       
3688                                        type[index] = 0;
3689                                        P0[index2D] = curData; 
3690                                        unpredictable_data[unpredictable_count ++] = curData;
3691                                }                               
3692                        }
3693                        else
3694                        {
3695                                type[index] = 0;
3696                                P0[index2D] = curData;
3697                                unpredictable_data[unpredictable_count ++] = curData;
3698                        }
3699
3700                        /* Process Row-i data 1 --> data r3-1 */
3701                        for (j = 1; j < r3; j++)
3702                        {
3703                                //index = k*r2*r3 + i*r3 + j;                   
3704                                index ++;
3705                                index2D = i*r3 + j;
3706                                pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
3707                                curData = cur_data_pos[j];
3708                                diff = curData - pred3D;
3709
3710                                itvNum = fabs(diff)/realPrecision + 1;
3711
3712                                if (itvNum < exe_params->intvCapacity)
3713                                {
3714                                        if (diff < 0) itvNum = -itvNum;
3715                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
3716                                        P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
3717                                       
3718                                        //ganrantee comporession error against the case of machine-epsilon
3719                                        if(fabs(curData-P0[index2D])>realPrecision)
3720                                        {       
3721                                                type[index] = 0;
3722                                                P0[index2D] = curData; 
3723                                                unpredictable_data[unpredictable_count ++] = curData;
3724                                        }                                       
3725                                }
3726                                else
3727                                {
3728                                        type[index] = 0;
3729                                        P0[index2D] = curData;
3730                                        unpredictable_data[unpredictable_count ++] = curData;
3731                                }
3732                        }
3733                        cur_data_pos += dim1_offset;
3734                }
3735                cur_data_pos += dim0_offset - r2 * dim1_offset;
3736                float *Pt;
3737                Pt = P1;
3738                P1 = P0;
3739                P0 = Pt;
3740        }
3741
3742        return unpredictable_count;
3743}
3744
3745unsigned int optimize_intervals_float_2D_opt(float *oriData, size_t r1, size_t r2, double realPrecision)
3746{       
3747        size_t i;
3748        size_t radiusIndex;
3749        float pred_value = 0, pred_err;
3750        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3751        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3752        size_t totalSampleSize = 0;//(r1-1)*(r2-1)/confparams_cpr->sampleDistance;
3753
3754        //float max = oriData[0];
3755        //float min = oriData[0];
3756
3757        size_t offset_count = confparams_cpr->sampleDistance - 1; // count r2 offset
3758        size_t offset_count_2;
3759        float * data_pos = oriData + r2 + offset_count;
3760        size_t n1_count = 1; // count i sum
3761        size_t len = r1 * r2;
3762        while(data_pos - oriData < len){
3763                totalSampleSize++;
3764                pred_value = data_pos[-1] + data_pos[-r2] - data_pos[-r2-1];
3765                pred_err = fabs(pred_value - *data_pos);
3766                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
3767                if(radiusIndex>=confparams_cpr->maxRangeRadius)
3768                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
3769                intervals[radiusIndex]++;
3770
3771                offset_count += confparams_cpr->sampleDistance;
3772                if(offset_count >= r2){
3773                        n1_count ++;
3774                        offset_count_2 = n1_count % confparams_cpr->sampleDistance;
3775                        data_pos += (r2 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
3776                        offset_count = (confparams_cpr->sampleDistance - offset_count_2);
3777                        if(offset_count == 0) offset_count ++;
3778                }
3779                else data_pos += confparams_cpr->sampleDistance;
3780        }
3781
3782        //compute the appropriate number
3783        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3784        size_t sum = 0;
3785        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3786        {
3787                sum += intervals[i];
3788                if(sum>targetCount)
3789                        break;
3790        }
3791        if(i>=confparams_cpr->maxRangeRadius)
3792                i = confparams_cpr->maxRangeRadius-1;
3793        unsigned int accIntervals = 2*(i+1);
3794        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3795
3796        if(powerOf2<32)
3797                powerOf2 = 32;
3798
3799        free(intervals);
3800        return powerOf2;
3801}
3802
3803unsigned int optimize_intervals_float_1D_opt(float *oriData, size_t dataLength, double realPrecision)
3804{       
3805        size_t i = 0, radiusIndex;
3806        float pred_value = 0, pred_err;
3807        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
3808        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
3809        size_t totalSampleSize = 0;//dataLength/confparams_cpr->sampleDistance;
3810
3811        float * data_pos = oriData + 2;
3812        while(data_pos - oriData < dataLength){
3813                totalSampleSize++;
3814                //pred_value = 2*data_pos[-1] - data_pos[-2];
3815                pred_value = data_pos[-1];
3816                pred_err = fabs(pred_value - *data_pos);
3817                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
3818                if(radiusIndex>=confparams_cpr->maxRangeRadius)
3819                        radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
3820                intervals[radiusIndex]++;
3821
3822                data_pos += confparams_cpr->sampleDistance;
3823        }
3824        //compute the appropriate number
3825        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
3826        size_t sum = 0;
3827        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
3828        {
3829                sum += intervals[i];
3830                if(sum>targetCount)
3831                        break;
3832        }
3833        if(i>=confparams_cpr->maxRangeRadius)
3834                i = confparams_cpr->maxRangeRadius-1;
3835               
3836        unsigned int accIntervals = 2*(i+1);
3837        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
3838       
3839        if(powerOf2<32)
3840                powerOf2 = 32;
3841       
3842        free(intervals);
3843        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
3844        return powerOf2;
3845}
3846
3847size_t SZ_compress_float_1D_MDQ_RA_block(float * block_ori_data, float * mean, size_t dim_0, size_t block_dim_0, double realPrecision, int * type, float * unpredictable_data){
3848
3849        mean[0] = block_ori_data[0];
3850        unsigned short unpredictable_count = 0;
3851
3852        float curData;
3853        double itvNum;
3854        double diff;
3855        float last_over_thres = mean[0];
3856        float pred1D;
3857        size_t type_index = 0;
3858        float * data_pos = block_ori_data;
3859        for(size_t i=0; i<block_dim_0; i++){
3860                curData = *data_pos;
3861
3862                pred1D = last_over_thres;
3863                diff = curData - pred1D;
3864                itvNum = fabs(diff)/realPrecision + 1;
3865                if (itvNum < exe_params->intvCapacity){
3866                        if (diff < 0) itvNum = -itvNum;
3867                        type[type_index] = (int) (itvNum/2) + exe_params->intvRadius;   
3868                        last_over_thres = pred1D + 2 * (type[type_index] - exe_params->intvRadius) * realPrecision;
3869                        if(fabs(curData-last_over_thres)>realPrecision){
3870                                type[type_index] = 0;
3871                                last_over_thres = curData;
3872                                unpredictable_data[unpredictable_count ++] = curData;
3873                        }
3874
3875                }
3876                else{
3877                        type[type_index] = 0;
3878                        unpredictable_data[unpredictable_count ++] = curData;
3879                        last_over_thres = curData;
3880                }
3881                type_index ++;
3882                data_pos ++;
3883        }
3884        return unpredictable_count;
3885
3886}
3887
3888size_t SZ_compress_float_2D_MDQ_RA_block(float * block_ori_data, float * mean, size_t dim_0, size_t dim_1, size_t block_dim_0, size_t block_dim_1, double realPrecision, float * P0, float * P1, int * type, float * unpredictable_data){
3889
3890        size_t dim0_offset = dim_1;
3891        mean[0] = block_ori_data[0];
3892
3893        size_t unpredictable_count = 0;
3894        size_t r1, r2;
3895        r1 = block_dim_0;
3896        r2 = block_dim_1;
3897
3898        float * cur_data_pos = block_ori_data;
3899        float curData;
3900        float pred1D, pred2D;
3901        double itvNum;
3902        double diff;
3903        size_t i, j;
3904        /* Process Row-0 data 0*/
3905        curData = *cur_data_pos;
3906        pred1D = mean[0];
3907        diff = curData - pred1D;
3908        itvNum = fabs(diff)/realPrecision + 1;
3909        if (itvNum < exe_params->intvCapacity){
3910                if (diff < 0) itvNum = -itvNum;
3911                type[0] = (int) (itvNum/2) + exe_params->intvRadius;
3912                P1[0] = pred1D + 2 * (type[0] - exe_params->intvRadius) * realPrecision;
3913                //ganrantee comporession error against the case of machine-epsilon
3914                if(fabs(curData-P1[0])>realPrecision){ 
3915                        type[0] = 0;
3916                        P1[0] = curData;
3917                        unpredictable_data[unpredictable_count ++] = curData;
3918                }               
3919        }
3920        else{
3921                type[0] = 0;
3922                P1[0] = curData;
3923                unpredictable_data[unpredictable_count ++] = curData;
3924        }
3925
3926        /* Process Row-0 data 1*/
3927        curData = cur_data_pos[1];
3928        pred1D = P1[0];
3929        diff = curData - pred1D;
3930        itvNum = fabs(diff)/realPrecision + 1;
3931        if (itvNum < exe_params->intvCapacity){
3932                if (diff < 0) itvNum = -itvNum;
3933                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
3934                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
3935                //ganrantee comporession error against the case of machine-epsilon
3936                if(fabs(curData-P1[1])>realPrecision){ 
3937                        type[1] = 0;
3938                        P1[1] = curData;       
3939                        unpredictable_data[unpredictable_count ++] = curData;
3940                }               
3941        }
3942        else{
3943                type[1] = 0;
3944                P1[1] = curData;
3945                unpredictable_data[unpredictable_count ++] = curData;
3946        }
3947
3948    /* Process Row-0 data 2 --> data r2-1 */
3949        for (j = 2; j < r2; j++)
3950        {
3951                curData = cur_data_pos[j];
3952                pred1D = 2*P1[j-1] - P1[j-2];
3953                diff = curData - pred1D;
3954                itvNum = fabs(diff)/realPrecision + 1;
3955                if (itvNum < exe_params->intvCapacity){
3956                        if (diff < 0) itvNum = -itvNum;
3957                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
3958                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
3959                        //ganrantee comporession error against the case of machine-epsilon
3960                        if(fabs(curData-P1[j])>realPrecision){ 
3961                                type[j] = 0;
3962                                P1[j] = curData;       
3963                                unpredictable_data[unpredictable_count ++] = curData;
3964                        }                       
3965                }
3966                else{
3967                        type[j] = 0;
3968                        P1[j] = curData;
3969                        unpredictable_data[unpredictable_count ++] = curData;
3970                }
3971        }
3972        cur_data_pos += dim0_offset;
3973        /* Process Row-1 --> Row-r1-1 */
3974        size_t index;
3975        for (i = 1; i < r1; i++)
3976        {       
3977                /* Process row-i data 0 */
3978                index = i*r2;
3979                curData = *cur_data_pos;
3980                pred1D = P1[0];
3981                diff = curData - pred1D;
3982                itvNum = fabs(diff)/realPrecision + 1;
3983                if (itvNum < exe_params->intvCapacity){
3984                        if (diff < 0) itvNum = -itvNum;
3985                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
3986                        P0[0] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
3987                        //ganrantee comporession error against the case of machine-epsilon
3988                        if(fabs(curData-P0[0])>realPrecision){ 
3989                                type[index] = 0;
3990                                P0[0] = curData;       
3991                                unpredictable_data[unpredictable_count ++] = curData;
3992                        }                       
3993                }
3994                else{
3995                        type[index] = 0;
3996                        P0[0] = curData;
3997                        unpredictable_data[unpredictable_count ++] = curData;
3998                }
3999                                                                       
4000                /* Process row-i data 1 --> r2-1*/
4001                for (j = 1; j < r2; j++)
4002                {
4003                        index = i*r2+j;
4004                        curData = cur_data_pos[j];
4005                        pred2D = P0[j-1] + P1[j] - P1[j-1];
4006                        diff = curData - pred2D;
4007                        itvNum = fabs(diff)/realPrecision + 1;
4008                        if (itvNum < exe_params->intvCapacity)
4009                        {
4010                                if (diff < 0) itvNum = -itvNum;
4011                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
4012                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
4013                               
4014                                //ganrantee comporession error against the case of machine-epsilon
4015                                if(fabs(curData-P0[j])>realPrecision)
4016                                {       
4017                                        type[index] = 0;
4018                                        P0[j] = curData;       
4019                                        unpredictable_data[unpredictable_count ++] = curData;
4020                                }                               
4021                        }
4022                        else
4023                        {
4024                                type[index] = 0;
4025                                P0[j] = curData;
4026                                unpredictable_data[unpredictable_count ++] = curData;
4027                        }
4028                }
4029                cur_data_pos += dim0_offset;
4030
4031                float *Pt;
4032                Pt = P1;
4033                P1 = P0;
4034                P0 = Pt;
4035        }
4036        return unpredictable_count;
4037}
4038
Note: See TracBrowser for help on using the repository browser.