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

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

more work on adding SZ (latest version)

  • Property mode set to 100644
RevLine 
[2c47b73]1/**
2 *  @file sz_double_pwr.c
3 *  @author Sheng Di
4 *  @date Aug, 2016
5 *  @brief SZ_Init, Compression and Decompression functions
6 * This file contains the compression/decompression functions related to point-wise relative errors
7 *  (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
8 *      See COPYRIGHT in top-level directory.
9 */
10
11
12#include <stdio.h>
13#include <stdlib.h>
14#include <string.h>
15#include <unistd.h>
16#include <math.h>
17#include "sz.h"
18#include "CompressElement.h"
19#include "DynamicByteArray.h"
20#include "DynamicIntArray.h"
21#include "TightDataPointStorageD.h"
22#include "sz_double.h"
23#include "sz_double_pwr.h"
24#include "zlib.h"
25#include "rw.h"
26
27void compute_segment_precisions_double_1D(double *oriData, size_t dataLength, double* pwrErrBound, unsigned char* pwrErrBoundBytes, double globalPrecision)
28{
29        size_t i = 0, j = 0, k = 0;
30        double realPrecision = oriData[0]!=0?fabs(confparams_cpr->pw_relBoundRatio*oriData[0]):confparams_cpr->pw_relBoundRatio; 
31        double approxPrecision;
32        unsigned char realPrecBytes[8];
33        double curPrecision;
34        double curValue;
35        double sum = 0;
36        for(i=0;i<dataLength;i++)
37        {
38                curValue = oriData[i];
39                if(i%confparams_cpr->segment_size==0&&i>0)
40                {
41                        //get two first bytes of the realPrecision
42                        if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE)
43                        {
44                                realPrecision = sum/confparams_cpr->segment_size;
45                                sum = 0;                       
46                        }
47                        realPrecision *= confparams_cpr->pw_relBoundRatio;
48                        if(confparams_cpr->errorBoundMode==ABS_AND_PW_REL||confparams_cpr->errorBoundMode==REL_AND_PW_REL)
49                                realPrecision = realPrecision<globalPrecision?realPrecision:globalPrecision; 
50                        else if(confparams_cpr->errorBoundMode==ABS_OR_PW_REL||confparams_cpr->errorBoundMode==REL_OR_PW_REL)
51                                realPrecision = realPrecision<globalPrecision?globalPrecision:realPrecision;
52                               
53                        doubleToBytes(realPrecBytes, realPrecision);
54                        memset(&realPrecBytes[2], 0, 6);
55                        approxPrecision = bytesToDouble(realPrecBytes);
56                        //put the realPrecision in double* pwrErBound
57                        pwrErrBound[j++] = approxPrecision;
58                        //put the two bytes in pwrErrBoundBytes
59                        pwrErrBoundBytes[k++] = realPrecBytes[0];
60                        pwrErrBoundBytes[k++] = realPrecBytes[1];
61                       
62                        realPrecision = fabs(curValue);
63                }
64               
65                if(curValue!=0)
66                {
67                        curPrecision = fabs(curValue);
68                       
69                        switch(confparams_cpr->pwr_type)
70                        {
71                        case SZ_PWR_MIN_TYPE:
72                                if(realPrecision>curPrecision)
73                                        realPrecision = curPrecision;   
74                                break;
75                        case SZ_PWR_AVG_TYPE:
76                                sum += curPrecision;
77                                break;
78                        case SZ_PWR_MAX_TYPE:
79                                if(realPrecision<curPrecision)
80                                        realPrecision = curPrecision;                                   
81                                break;
82                        }
83                }
84        }
85        if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE)
86        {
87                int size = dataLength%confparams_cpr->segment_size==0?confparams_cpr->segment_size:dataLength%confparams_cpr->segment_size;
88                realPrecision = sum/size;               
89        }       
90        if(confparams_cpr->errorBoundMode==ABS_AND_PW_REL||confparams_cpr->errorBoundMode==REL_AND_PW_REL)
91                realPrecision = realPrecision<globalPrecision?realPrecision:globalPrecision; 
92        else if(confparams_cpr->errorBoundMode==ABS_OR_PW_REL||confparams_cpr->errorBoundMode==REL_OR_PW_REL)
93                realPrecision = realPrecision<globalPrecision?globalPrecision:realPrecision;   
94        doubleToBytes(realPrecBytes, realPrecision);
95        memset(&realPrecBytes[2], 0, 6);
96        approxPrecision = bytesToDouble(realPrecBytes);
97        //put the realPrecision in double* pwrErBound
98        pwrErrBound[j++] = approxPrecision;
99        //put the two bytes in pwrErrBoundBytes
100        pwrErrBoundBytes[k++] = realPrecBytes[0];
101        pwrErrBoundBytes[k++] = realPrecBytes[1];
102}
103
104unsigned int optimize_intervals_double_1D_pwr(double *oriData, size_t dataLength, double* pwrErrBound)
105{       
106        size_t i = 0, j = 0;
107        double realPrecision = pwrErrBound[j++];       
108        unsigned long radiusIndex;
109        double pred_value = 0, pred_err;
110        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
111        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
112        int totalSampleSize = dataLength/confparams_cpr->sampleDistance;
113        for(i=2;i<dataLength;i++)
114        {
115                if(i%confparams_cpr->segment_size==0)
116                        realPrecision = pwrErrBound[j++];
117                if(i%confparams_cpr->sampleDistance==0)
118                {
119                        //pred_value = 2*oriData[i-1] - oriData[i-2];
120                        pred_value = oriData[i-1];
121                        pred_err = fabs(pred_value - oriData[i]);
122                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
123                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
124                                radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
125                        intervals[radiusIndex]++;
126                }
127        }
128        //compute the appropriate number
129        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
130        size_t sum = 0;
131        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
132        {
133                sum += intervals[i];
134                if(sum>targetCount)
135                        break;
136        }
137        if(i>=confparams_cpr->maxRangeRadius)
138                i = confparams_cpr->maxRangeRadius-1;
139        unsigned int accIntervals = 2*(i+1);
140        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
141       
142        if(powerOf2<32)
143                powerOf2 = 32;
144       
145        free(intervals);
146        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
147        return powerOf2;
148}
149
150void compute_segment_precisions_double_2D(double *oriData, double* pwrErrBound, 
151size_t r1, size_t r2, size_t R2, size_t edgeSize, unsigned char* pwrErrBoundBytes, double Min, double Max, double globalPrecision)
152{
153        size_t i = 0, j = 0, k = 0, p = 0, index = 0, J = 0; //I=-1,J=-1 if they are needed
154        double realPrecision; 
155        double approxPrecision;
156        unsigned char realPrecBytes[8];
157        double curValue, curAbsValue;
158        double* statAbsValues = (double*)malloc(R2*sizeof(double));     
159       
160        double max = fabs(Min)<fabs(Max)?fabs(Max):fabs(Min); //get the max abs value.
161        double min = fabs(Min)<fabs(Max)?fabs(Min):fabs(Max);
162        for(i=0;i<R2;i++)
163        {
164                if(confparams_cpr->pwr_type == SZ_PWR_MIN_TYPE)
165                        statAbsValues[i] = max;
166                else if(confparams_cpr->pwr_type == SZ_PWR_MAX_TYPE)
167                        statAbsValues[i] = min;
168                else
169                        statAbsValues[i] = 0; //for SZ_PWR_AVG_TYPE
170        }
171        for(i=0;i<r1;i++)
172        {
173                for(j=0;j<r2;j++)
174                {
175                        index = i*r2+j;
176                        curValue = oriData[index];                             
177                        if(((i%edgeSize==edgeSize-1 || i==r1-1) &&j%edgeSize==0&&j>0) || (i%edgeSize==0&&j==0&&i>0))
178                        {
179                                if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE)
180                                {
181                                        int a = edgeSize, b = edgeSize;
182                                        if(j==0)
183                                        {
184                                                if(r2%edgeSize==0) 
185                                                        b = edgeSize;
186                                                else
187                                                        b = r2%edgeSize;
188                                        }
189                                        if(i==r1-1)
190                                        {
191                                                if(r1%edgeSize==0)
192                                                        a = edgeSize;
193                                                else
194                                                        a = r1%edgeSize;
195                                        }
196                                        realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J]/(a*b);
197                                }
198                                else
199                                        realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J];
200
201                                if(confparams_cpr->errorBoundMode==ABS_AND_PW_REL||confparams_cpr->errorBoundMode==REL_AND_PW_REL)
202                                        realPrecision = realPrecision<globalPrecision?realPrecision:globalPrecision; 
203                                else if(confparams_cpr->errorBoundMode==ABS_OR_PW_REL||confparams_cpr->errorBoundMode==REL_OR_PW_REL)
204                                        realPrecision = realPrecision<globalPrecision?globalPrecision:realPrecision;
205                                       
206                                doubleToBytes(realPrecBytes, realPrecision);
207                                memset(&realPrecBytes[2], 0, 6);
208                                approxPrecision = bytesToDouble(realPrecBytes);
209                                //put the realPrecision in double* pwrErBound
210                                pwrErrBound[p++] = approxPrecision;
211                                //put the two bytes in pwrErrBoundBytes
212                                pwrErrBoundBytes[k++] = realPrecBytes[0];
213                                pwrErrBoundBytes[k++] = realPrecBytes[1];       
214                               
215                                if(confparams_cpr->pwr_type == SZ_PWR_MIN_TYPE)
216                                        statAbsValues[J] = max;
217                                else if(confparams_cpr->pwr_type == SZ_PWR_MAX_TYPE)
218                                        statAbsValues[J] = min;
219                                else
220                                        statAbsValues[J] = 0; //for SZ_PWR_AVG_TYPE             
221                        }       
222                        if(j==0)
223                                J = 0;
224                        else if(j%edgeSize==0)
225                                J++;                   
226                        if(curValue!=0)
227                        {
228                                curAbsValue = fabs(curValue);
229                               
230                                switch(confparams_cpr->pwr_type)
231                                {
232                                case SZ_PWR_MIN_TYPE:
233                                        if(statAbsValues[J]>curAbsValue)
234                                                statAbsValues[J] = curAbsValue; 
235                                        break;
236                                case SZ_PWR_AVG_TYPE:
237                                        statAbsValues[J] += curAbsValue;
238                                        break;
239                                case SZ_PWR_MAX_TYPE:
240                                        if(statAbsValues[J]<curAbsValue)
241                                                statAbsValues[J] = curAbsValue;                                 
242                                        break;
243                                }
244                        }
245                }
246        }
247               
248        if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE)
249        {
250                int a = edgeSize, b = edgeSize;
251                if(r2%edgeSize==0) 
252                        b = edgeSize;
253                else
254                        b = r2%edgeSize;
255                if(r1%edgeSize==0)
256                        a = edgeSize;
257                else
258                        a = r1%edgeSize;
259                realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J]/(a*b);
260        }
261        else
262                realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J];             
263
264        if(confparams_cpr->errorBoundMode==ABS_AND_PW_REL||confparams_cpr->errorBoundMode==REL_AND_PW_REL)
265                realPrecision = realPrecision<globalPrecision?realPrecision:globalPrecision; 
266        else if(confparams_cpr->errorBoundMode==ABS_OR_PW_REL||confparams_cpr->errorBoundMode==REL_OR_PW_REL)
267                realPrecision = realPrecision<globalPrecision?globalPrecision:realPrecision;
268               
269        doubleToBytes(realPrecBytes, realPrecision);
270        realPrecBytes[2] = realPrecBytes[3] = 0;
271        approxPrecision = bytesToDouble(realPrecBytes);
272        //put the realPrecision in double* pwrErBound
273        pwrErrBound[p++] = approxPrecision;
274        //put the two bytes in pwrErrBoundBytes
275        pwrErrBoundBytes[k++] = realPrecBytes[0];
276        pwrErrBoundBytes[k++] = realPrecBytes[1];       
277       
278        free(statAbsValues);
279}
280
281unsigned int optimize_intervals_double_2D_pwr(double *oriData, size_t r1, size_t r2, size_t R2, size_t edgeSize, double* pwrErrBound)
282{       
283        size_t i = 0,j = 0, index, I=0, J=0;
284        double realPrecision = pwrErrBound[0]; 
285        unsigned long radiusIndex;
286        double pred_value = 0, pred_err;
287        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
288        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
289        size_t totalSampleSize = (r1-1)*(r2-1)/confparams_cpr->sampleDistance;
290        size_t ir2;
291        for(i=1;i<r1;i++)
292        {
293                ir2 = i*r2;
294                if(i%edgeSize==0)
295                {       
296                        I++;
297                        J = 0;
298                }
299                for(j=1;j<r2;j++)
300                {
301                        index = ir2+j;
302                        if(j%edgeSize==0)
303                                J++;
304                               
305                        if((i+j)%confparams_cpr->sampleDistance==0)
306                        {
307                                realPrecision = pwrErrBound[I*R2+J];
308                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
309                                pred_err = fabs(pred_value - oriData[index]);
310                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
311                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
312                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
313                                intervals[radiusIndex]++;
314                        }                       
315                }
316        }
317        //compute the appropriate number
318        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
319        size_t sum = 0;
320        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
321        {
322                sum += intervals[i];
323                if(sum>targetCount)
324                        break;
325        }
326        if(i>=confparams_cpr->maxRangeRadius)
327                i = confparams_cpr->maxRangeRadius-1;
328        unsigned int accIntervals = 2*(i+1);
329        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
330
331        if(powerOf2<32)
332                powerOf2 = 32;
333
334        free(intervals);
335        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
336        return powerOf2;
337}
338
339void compute_segment_precisions_double_3D(double *oriData, double* pwrErrBound, 
340size_t r1, size_t r2, size_t r3, size_t R2, size_t R3, size_t edgeSize, unsigned char* pwrErrBoundBytes, double Min, double Max, double globalPrecision)
341{
342        size_t i = 0, j = 0, k = 0, p = 0, q = 0, index = 0, J = 0, K = 0; //I=-1,J=-1 if they are needed
343        size_t r23 = r2*r3, ir, jr;
344        double realPrecision; 
345        double approxPrecision;
346        unsigned char realPrecBytes[8];
347        double curValue, curAbsValue;
348       
349        double** statAbsValues = create2DArray_double(R2, R3);
350        double max = fabs(Min)<fabs(Max)?fabs(Max):fabs(Min); //get the max abs value. 
351        double min = fabs(Min)<fabs(Max)?fabs(Min):fabs(Max);
352        for(i=0;i<R2;i++)
353                for(j=0;j<R3;j++)
354                {
355                        if(confparams_cpr->pwr_type == SZ_PWR_MIN_TYPE)
356                                statAbsValues[i][j] = max;
357                        else if(confparams_cpr->pwr_type == SZ_PWR_MAX_TYPE)
358                                statAbsValues[i][j] = min;
359                        else
360                                statAbsValues[i][j] = 0;
361                }
362        for(i=0;i<r1;i++)
363        {
364                ir = i*r23;             
365                if(i%edgeSize==0&&i>0)
366                {
367                        realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J][K];
368                        doubleToBytes(realPrecBytes, realPrecision);
369                        memset(&realPrecBytes[2], 0, 6);
370                        approxPrecision = bytesToDouble(realPrecBytes);
371                        //put the realPrecision in double* pwrErBound
372                        pwrErrBound[p++] = approxPrecision;
373                        //put the two bytes in pwrErrBoundBytes
374                        //printf("q=%d, i=%d, j=%d, k=%d\n",q,i,j,k);
375                        pwrErrBoundBytes[q++] = realPrecBytes[0];
376                        pwrErrBoundBytes[q++] = realPrecBytes[1];
377                        if(confparams_cpr->pwr_type == SZ_PWR_MIN_TYPE)
378                                statAbsValues[J][K] = max;
379                        else if(confparams_cpr->pwr_type == SZ_PWR_MAX_TYPE)
380                                statAbsValues[J][K] = min;
381                }               
382                for(j=0;j<r2;j++)
383                {
384                        jr = j*r3;
385                        if((i%edgeSize==edgeSize-1 || i == r1-1)&&j%edgeSize==0&&j>0)
386                        {
387                                realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J][K];
388                                doubleToBytes(realPrecBytes, realPrecision);
389                                memset(&realPrecBytes[2], 0, 6);
390                                approxPrecision = bytesToDouble(realPrecBytes);
391                                //put the realPrecision in double* pwrErBound
392                                pwrErrBound[p++] = approxPrecision;
393                                //put the two bytes in pwrErrBoundBytes
394                                //printf("q=%d, i=%d, j=%d, k=%d\n",q,i,j,k);
395                                pwrErrBoundBytes[q++] = realPrecBytes[0];
396                                pwrErrBoundBytes[q++] = realPrecBytes[1];
397                                if(confparams_cpr->pwr_type == SZ_PWR_MIN_TYPE)
398                                        statAbsValues[J][K] = max;
399                                else if(confparams_cpr->pwr_type == SZ_PWR_MAX_TYPE)
400                                        statAbsValues[J][K] = min;                     
401                        }
402                       
403                        if(j==0)
404                                J = 0;
405                        else if(j%edgeSize==0)
406                                J++;                                   
407                       
408                        for(k=0;k<r3;k++)
409                        {
410                                index = ir+jr+k;                               
411                                curValue = oriData[index];                             
412                                if((i%edgeSize==edgeSize-1 || i == r1-1)&&(j%edgeSize==edgeSize-1||j==r2-1)&&k%edgeSize==0&&k>0)
413                                {
414                                        realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J][K];
415                                        doubleToBytes(realPrecBytes, realPrecision);
416                                        memset(&realPrecBytes[2], 0, 6);
417                                        approxPrecision = bytesToDouble(realPrecBytes);
418                                        //put the realPrecision in double* pwrErBound
419                                        pwrErrBound[p++] = approxPrecision;
420                                        //put the two bytes in pwrErrBoundBytes
421                                        //printf("q=%d, i=%d, j=%d, k=%d\n",q,i,j,k);
422                                        pwrErrBoundBytes[q++] = realPrecBytes[0];
423                                        pwrErrBoundBytes[q++] = realPrecBytes[1];
424                                       
425                                        if(confparams_cpr->pwr_type == SZ_PWR_MIN_TYPE)
426                                                statAbsValues[J][K] = max;
427                                        else if(confparams_cpr->pwr_type == SZ_PWR_MAX_TYPE)
428                                                statAbsValues[J][K] = min;     
429                                }       
430
431                                if(k==0)
432                                        K = 0;
433                                else if(k%edgeSize==0)
434                                        K++;
435                                       
436                                if(curValue!=0)
437                                {
438                                        curAbsValue = fabs(curValue);
439                                        if(confparams_cpr->pwr_type == SZ_PWR_MIN_TYPE)
440                                        {
441                                                if(statAbsValues[J][K]>curAbsValue)
442                                                {
443                                                        statAbsValues[J][K] = curAbsValue;
444                                                }
445                                        }
446                                        else if(confparams_cpr->pwr_type == SZ_PWR_MAX_TYPE)
447                                        {
448                                                if(statAbsValues[J][K]<curAbsValue)
449                                                {
450                                                        statAbsValues[J][K] = curAbsValue;
451                                                }
452                                        }
453                                }
454                        }                       
455                }
456        }       
457       
458        realPrecision = confparams_cpr->pw_relBoundRatio*statAbsValues[J][K];
459        doubleToBytes(realPrecBytes, realPrecision);
460        memset(&realPrecBytes[2], 0, 6);
461        approxPrecision = bytesToDouble(realPrecBytes);
462        //put the realPrecision in double* pwrErBound
463        pwrErrBound[p++] = approxPrecision;
464        //put the two bytes in pwrErrBoundBytes
465        pwrErrBoundBytes[q++] = realPrecBytes[0];
466        pwrErrBoundBytes[q++] = realPrecBytes[1];
467       
468        free2DArray_double(statAbsValues, R2);
469}
470
471unsigned int optimize_intervals_double_3D_pwr(double *oriData, size_t r1, size_t r2, size_t r3, size_t R2, size_t R3, size_t edgeSize, double* pwrErrBound)
472{       
473        size_t i,j,k, ir,jr,index, I = 0,J=0,K=0;
474        double realPrecision = pwrErrBound[0];         
475        unsigned long radiusIndex;
476        size_t r23=r2*r3;
477        size_t R23 = R2*R3;
478        double pred_value = 0, pred_err;
479        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
480        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
481        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)/confparams_cpr->sampleDistance;
482        for(i=1;i<r1;i++)
483        {
484                ir = i*r23;
485                if(i%edgeSize==0)
486                {       
487                        I++;
488                        J = 0;
489                }
490                for(j=1;j<r2;j++)
491                {
492                        jr = j*r3;
493                        if(j%edgeSize==0)
494                        {       
495                                J++;
496                                K = 0;
497                        }                       
498                        for(k=1;k<r3;k++)
499                        {
500                                index = ir+jr+k;
501                                if(k%edgeSize==0)
502                                        K++;           
503                                if((i+j+k)%confparams_cpr->sampleDistance==0)
504                                {
505                                        realPrecision = pwrErrBound[I*R23+J*R2+K];                                     
506                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23] 
507                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
508                                        pred_err = fabs(pred_value - oriData[index]);
509                                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
510                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
511                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
512                                        intervals[radiusIndex]++;
513                                }
514                        }
515                }
516        }
517        //compute the appropriate number
518        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
519        size_t sum = 0;
520        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
521        {
522                sum += intervals[i];
523                if(sum>targetCount)
524                        break;
525        }
526        if(i>=confparams_cpr->maxRangeRadius)
527                i = confparams_cpr->maxRangeRadius-1;
528        unsigned int accIntervals = 2*(i+1);
529        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
530
531        if(powerOf2<32)
532                powerOf2 = 32;
533       
534        free(intervals);
535        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
536        return powerOf2;
537}
538
539void SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr(unsigned char** newByteData, double *oriData, double globalPrecision, 
540size_t dataLength, size_t *outSize, double min, double max)
541{
542        size_t pwrLength = dataLength%confparams_cpr->segment_size==0?dataLength/confparams_cpr->segment_size:dataLength/confparams_cpr->segment_size+1;
543        double* pwrErrBound = (double*)malloc(sizeof(double)*pwrLength);
544        size_t pwrErrBoundBytes_size = sizeof(unsigned char)*pwrLength*2;
545        unsigned char* pwrErrBoundBytes = (unsigned char*)malloc(pwrErrBoundBytes_size);
546       
547        compute_segment_precisions_double_1D(oriData, dataLength, pwrErrBound, pwrErrBoundBytes, globalPrecision);
548
549        unsigned int quantization_intervals;
550        if(exe_params->optQuantMode==1)
551        {
552                quantization_intervals = optimize_intervals_double_1D_pwr(oriData, dataLength, pwrErrBound);   
553                updateQuantizationInfo(quantization_intervals);
554        }
555        else
556                quantization_intervals = exe_params->intvCapacity;
557        size_t i = 0, j = 0;
558        int reqLength;
559        double realPrecision = pwrErrBound[j++];       
560        double medianValue = 0;
561        double radius = fabs(max)<fabs(min)?fabs(min):fabs(max);
562        short radExpo = getExponent_double(radius);
563       
564        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);     
565
566        int* type = (int*) malloc(dataLength*sizeof(int));
567        //type[dataLength]=0;
568               
569        double* spaceFillingValue = oriData; //
570       
571        DynamicByteArray *resiBitLengthArray;
572        new_DBA(&resiBitLengthArray, DynArrayInitLen);
573       
574        DynamicIntArray *exactLeadNumArray;
575        new_DIA(&exactLeadNumArray, DynArrayInitLen);
576       
577        DynamicByteArray *exactMidByteArray;
578        new_DBA(&exactMidByteArray, DynArrayInitLen);
579       
580        DynamicIntArray *resiBitArray;
581        new_DIA(&resiBitArray, DynArrayInitLen);
582       
583        type[0] = 0;
584       
585        unsigned char preDataBytes[8] = {0};
586        intToBytes_bigEndian(preDataBytes, 0);
587       
588        int reqBytesLength = reqLength/8;
589        int resiBitsLength = reqLength%8;
590        double last3CmprsData[3] = {0};
591
592        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
593        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
594                                               
595        //add the first data   
596        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
597        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
598        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
599        memcpy(preDataBytes,vce->curBytes,8);
600        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
601        listAdd_double(last3CmprsData, vce->data);
602        //printf("%.30G\n",last3CmprsData[0]); 
603               
604        //add the second data
605        type[1] = 0;
606        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);                 
607        compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
608        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
609        memcpy(preDataBytes,vce->curBytes,8);
610        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
611        listAdd_double(last3CmprsData, vce->data);
612        //printf("%.30G\n",last3CmprsData[0]); 
613       
614        int state;
615        double checkRadius;
616        double curData;
617        double pred;
618        double predAbsErr;
619        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
620        double interval = 2*realPrecision;
621        int updateReqLength = 0; //a marker: 1 means already updated
622       
623        for(i=2;i<dataLength;i++)
624        {
625                curData = spaceFillingValue[i];
626                if(i%confparams_cpr->segment_size==0)
627                {
628                        realPrecision = pwrErrBound[j++];
629                        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
630                        interval = 2*realPrecision;
631                        updateReqLength = 0;
632                }
633                //pred = 2*last3CmprsData[0] - last3CmprsData[1];
634                pred = last3CmprsData[0];
635                predAbsErr = fabs(curData - pred);     
636                if(predAbsErr<checkRadius)
637                {
638                        state = (predAbsErr/realPrecision+1)/2;
639                        if(curData>=pred)
640                        {
641                                type[i] = exe_params->intvRadius+state;
642                                pred = pred + state*interval;
643                        }
644                        else //curData<pred
645                        {
646                                type[i] = exe_params->intvRadius-state;
647                                pred = pred - state*interval;
648                        }
649                        listAdd_double(last3CmprsData, pred);                   
650                        continue;
651                }
652               
653                //unpredictable data processing         
654                if(updateReqLength==0)
655                {
656                        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);                             
657                        reqBytesLength = reqLength/8;
658                        resiBitsLength = reqLength%8;
659                        updateReqLength = 1;           
660                }
661               
662                type[i] = 0;
663                addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
664               
665                compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
666                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
667                memcpy(preDataBytes,vce->curBytes,8);
668                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
669
670                listAdd_double(last3CmprsData, vce->data);     
671        }//end of for
672               
673//      char* expSegmentsInBytes;
674//      int expSegmentsInBytes_size = convertESCToBytes(esc, &expSegmentsInBytes);
675        int exactDataNum = exactLeadNumArray->size;
676       
677        TightDataPointStorageD* tdps;
678                       
679        new_TightDataPointStorageD2(&tdps, dataLength, exactDataNum, 
680                        type, exactMidByteArray->array, exactMidByteArray->size, 
681                        exactLeadNumArray->array, 
682                        resiBitArray->array, resiBitArray->size, 
683                        resiBitLengthArray->array, resiBitLengthArray->size, 
684                        realPrecision, medianValue, (char)reqLength, quantization_intervals, pwrErrBoundBytes, pwrErrBoundBytes_size, radExpo);
685
686//sdi:Debug
687/*      int sum =0;
688        for(i=0;i<dataLength;i++)
689                if(type[i]==0) sum++;
690        printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);
691*/
692//      writeUShortData(type, dataLength, "compressStateBytes.sb");
693//      unsigned short type_[dataLength];
694//      SZ_Reset();
695//      decode_withTree(tdps->typeArray, tdps->typeArray_size, type_); 
696//      printf("tdps->typeArray_size=%d\n", tdps->typeArray_size);
697               
698        //free memory
699        free_DBA(resiBitLengthArray);
700        free_DIA(exactLeadNumArray);
701        free_DIA(resiBitArray);
702        free(type);
703               
704        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
705       
706        int doubleSize=sizeof(double);
707        if(*outSize>dataLength*doubleSize)
708        {
709                size_t k = 0, i;
710                tdps->isLossless = 1;
711                size_t totalByteLength = 3 + exe_params->SZ_SIZE_TYPE + 1 + doubleSize*dataLength;
712                *newByteData = (unsigned char*)malloc(totalByteLength);
713               
714                unsigned char dsLengthBytes[exe_params->SZ_SIZE_TYPE];
715                intToBytes_bigEndian(dsLengthBytes, dataLength);//4
716                for (i = 0; i < 3; i++)//3
717                        (*newByteData)[k++] = versionNumber[i];
718               
719                if(exe_params->SZ_SIZE_TYPE==4)
720                {
721                        (*newByteData)[k++] = 16;       //=00010000     
722                }
723                else 
724                {
725                        (*newByteData)[k++] = 80;
726                }
727                for (i = 0; i < exe_params->SZ_SIZE_TYPE; i++)//4 or 8
728                        (*newByteData)[k++] = dsLengthBytes[i];
729
730               
731                if(sysEndianType==BIG_ENDIAN_SYSTEM)
732                        memcpy((*newByteData)+4+exe_params->SZ_SIZE_TYPE, oriData, dataLength*doubleSize);
733                else
734                {
735                        unsigned char* p = (*newByteData)+4+exe_params->SZ_SIZE_TYPE;
736                        for(i=0;i<dataLength;i++,p+=doubleSize)
737                                doubleToBytes(p, oriData[i]);
738                }
739                *outSize = totalByteLength;
740        }
741       
742        free(pwrErrBound);
743       
744        free(vce);
745        free(lce);
746        free_TightDataPointStorageD(tdps);
747        free(exactMidByteArray);
748}
749
750
751/**
752 *
753 * Note: @r1 is high dimension
754 *               @r2 is low dimension
755 * */
756void SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr(unsigned char** newByteData, double *oriData, double globalPrecision, size_t r1, size_t r2,
757size_t *outSize, double min, double max)
758{
759        size_t dataLength=r1*r2;
760        int blockEdgeSize = computeBlockEdgeSize_2D(confparams_cpr->segment_size);
761        size_t R1 = 1+(r1-1)/blockEdgeSize;
762        size_t R2 = 1+(r2-1)/blockEdgeSize;
763        double* pwrErrBound = (double*)malloc(sizeof(double)*R1*R2);
764        size_t pwrErrBoundBytes_size = sizeof(unsigned char)*R1*R2*2;
765        unsigned char* pwrErrBoundBytes = (unsigned char*)malloc(pwrErrBoundBytes_size);
766       
767        compute_segment_precisions_double_2D(oriData, pwrErrBound, r1, r2, R2, blockEdgeSize, pwrErrBoundBytes, min, max, globalPrecision);
768        unsigned int quantization_intervals;
769        if(exe_params->optQuantMode==1)
770        {
771                quantization_intervals = optimize_intervals_double_2D_pwr(oriData, r1, r2, R2, blockEdgeSize, pwrErrBound);
772                updateQuantizationInfo(quantization_intervals);
773        }
774        else
775                quantization_intervals = exe_params->intvCapacity;     
776        //printf("quantization_intervals=%d\n",quantization_intervals);
777       
778        size_t i=0,j=0,I=0,J=0; 
779        int reqLength;
780        double realPrecision = pwrErrBound[I*R2+J];     
781        double pred1D, pred2D;
782        double diff = 0.0;
783        double itvNum = 0;
784        double *P0, *P1;
785       
786        P0 = (double*)malloc(r2*sizeof(double));
787        memset(P0, 0, r2*sizeof(double));
788        P1 = (double*)malloc(r2*sizeof(double));
789        memset(P1, 0, r2*sizeof(double));
790               
791        double medianValue = 0;
792        double radius = fabs(max)<fabs(min)?fabs(min):fabs(max);
793        short radExpo = getExponent_double(radius);
794        int updateReqLength = 1;
795       
796        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
797
798        int* type = (int*) malloc(dataLength*sizeof(int));
799        //type[dataLength]=0;
800               
801        double* spaceFillingValue = oriData; //
802       
803        DynamicByteArray *resiBitLengthArray;
804        new_DBA(&resiBitLengthArray, DynArrayInitLen);
805       
806        DynamicIntArray *exactLeadNumArray;
807        new_DIA(&exactLeadNumArray, DynArrayInitLen);
808       
809        DynamicByteArray *exactMidByteArray;
810        new_DBA(&exactMidByteArray, DynArrayInitLen);
811       
812        DynamicIntArray *resiBitArray;
813        new_DIA(&resiBitArray, DynArrayInitLen);
814       
815        type[0] = 0;
816       
817        unsigned char preDataBytes[8];
818        longToBytes_bigEndian(preDataBytes, 0);
819       
820        int reqBytesLength = reqLength/8;
821        int resiBitsLength = reqLength%8;
822
823        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
824        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
825                       
826        /* Process Row-0 data 0*/
827        type[0] = 0;
828        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
829        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
830        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
831        memcpy(preDataBytes,vce->curBytes,8);
832        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
833        P1[0] = vce->data;
834
835        /* Process Row-0 data 1*/
836        pred1D = P1[0];
837        diff = spaceFillingValue[1] - pred1D;
838
839        itvNum =  fabs(diff)/realPrecision + 1;
840
841        if (itvNum < exe_params->intvCapacity)
842        {
843                if (diff < 0) itvNum = -itvNum;
844                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
845                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
846        }
847        else
848        {               
849                type[1] = 0;
850
851                addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
852                compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
853                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
854                memcpy(preDataBytes,vce->curBytes,8);
855                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
856                P1[1] = vce->data;
857        }
858
859    /* Process Row-0 data 2 --> data r2-1 */
860        for (j = 2; j < r2; j++)
861        {
862                if(j%blockEdgeSize==0)
863                {
864                        J++;
865                        realPrecision = pwrErrBound[I*R2+J];
866                        updateReqLength = 0;
867                }
868
869                pred1D = 2*P1[j-1] - P1[j-2];
870                diff = spaceFillingValue[j] - pred1D;
871
872                itvNum = fabs(diff)/realPrecision + 1;
873
874                if (itvNum < exe_params->intvCapacity)
875                {
876                        if (diff < 0) itvNum = -itvNum;
877                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
878                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
879                }
880                else
881                {
882                        if(updateReqLength==0)
883                        {
884                                computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
885                                reqBytesLength = reqLength/8;
886                                resiBitsLength = reqLength%8;
887                                updateReqLength = 1;
888                        }
889
890                        type[j] = 0;
891
892                        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
893                        compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
894                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
895                        memcpy(preDataBytes,vce->curBytes,8);
896                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
897                        P1[j] = vce->data;
898                }
899        }
900
901        /* Process Row-1 --> Row-r1-1 */
902        size_t index;
903        for (i = 1; i < r1; i++)
904        {       
905                /* Process row-i data 0 */
906                index = i*r2;
907                J = 0;
908                if(i%blockEdgeSize==0)
909                        I++;
910                realPrecision = pwrErrBound[I*R2+J]; //J==0
911                updateReqLength = 0;
912               
913                pred1D = P1[0];
914                diff = spaceFillingValue[index] - pred1D;
915
916                itvNum = fabs(diff)/realPrecision + 1;
917
918                if (itvNum < exe_params->intvCapacity)
919                {
920                        if (diff < 0) itvNum = -itvNum;
921                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
922                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
923                }
924                else
925                {
926                        if(updateReqLength==0)
927                        {
928                                computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
929                                reqBytesLength = reqLength/8;
930                                resiBitsLength = reqLength%8;
931                                updateReqLength = 1;
932                        }
933                       
934                        type[index] = 0;
935
936                        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
937                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
938                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
939                        memcpy(preDataBytes,vce->curBytes,8);
940                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
941                        P0[0] = vce->data;
942                }
943                                                                       
944                /* Process row-i data 1 --> r2-1*/
945                for (j = 1; j < r2; j++)
946                {
947                        index = i*r2+j;
948                        if(j%blockEdgeSize==0)
949                        {
950                                J++;
951                                realPrecision = pwrErrBound[I*R2+J];
952                                updateReqLength = 0;
953                        }
954                        pred2D = P0[j-1] + P1[j] - P1[j-1];
955
956                        diff = spaceFillingValue[index] - pred2D;
957
958                        itvNum = fabs(diff)/realPrecision + 1;
959
960                        if (itvNum < exe_params->intvCapacity)
961                        {
962                                if (diff < 0) itvNum = -itvNum;
963                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
964                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
965                        }
966                        else
967                        {
968                                if(updateReqLength==0)
969                                {
970                                        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
971                                        reqBytesLength = reqLength/8;
972                                        resiBitsLength = reqLength%8;
973                                        updateReqLength = 1;
974                                }
975
976                                type[index] = 0;
977
978                                addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
979                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
980                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
981                                memcpy(preDataBytes,vce->curBytes,8);
982                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
983                                P0[j] = vce->data;
984                        }
985                }
986
987                double *Pt;
988                Pt = P1;
989                P1 = P0;
990                P0 = Pt;
991        }
992               
993        if(r2!=1)       
994                free(P0);
995        free(P1);
996        int exactDataNum = exactLeadNumArray->size;
997       
998        TightDataPointStorageD* tdps;
999                       
1000        new_TightDataPointStorageD2(&tdps, dataLength, exactDataNum, 
1001                        type, exactMidByteArray->array, exactMidByteArray->size, 
1002                        exactLeadNumArray->array, 
1003                        resiBitArray->array, resiBitArray->size, 
1004                        resiBitLengthArray->array, resiBitLengthArray->size, 
1005                        realPrecision, medianValue, (char)reqLength, quantization_intervals, pwrErrBoundBytes, pwrErrBoundBytes_size, radExpo);
1006       
1007        //free memory
1008        free_DBA(resiBitLengthArray);
1009        free_DIA(exactLeadNumArray);
1010        free_DIA(resiBitArray);
1011        free(type);
1012               
1013        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
1014
1015        free(pwrErrBound);
1016       
1017        free(vce);
1018        free(lce);
1019        free_TightDataPointStorageD(tdps);     
1020        free(exactMidByteArray);
1021}
1022
1023void SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(unsigned char** newByteData, double *oriData, double globalPrecision, 
1024size_t r1, size_t r2, size_t r3, size_t *outSize, double min, double max)
1025{
1026        size_t dataLength=r1*r2*r3;
1027       
1028        int blockEdgeSize = computeBlockEdgeSize_3D(confparams_cpr->segment_size);
1029        size_t R1 = 1+(r1-1)/blockEdgeSize;
1030        size_t R2 = 1+(r2-1)/blockEdgeSize;
1031        size_t R3 = 1+(r3-1)/blockEdgeSize;
1032        double* pwrErrBound = (double*)malloc(sizeof(double)*R1*R2*R3);
1033        size_t pwrErrBoundBytes_size = sizeof(unsigned char)*R1*R2*R3*2;
1034        unsigned char* pwrErrBoundBytes = (unsigned char*)malloc(pwrErrBoundBytes_size);       
1035       
1036        compute_segment_precisions_double_3D(oriData, pwrErrBound, r1, r2, r3, R2, R3, blockEdgeSize, pwrErrBoundBytes, min, max, globalPrecision);     
1037        unsigned int quantization_intervals;
1038        if(exe_params->optQuantMode==1)
1039        {
1040                quantization_intervals = optimize_intervals_double_3D_pwr(oriData, r1, r2, r3, R2, R3, blockEdgeSize, pwrErrBound);
1041                updateQuantizationInfo(quantization_intervals);
1042        }       
1043        else
1044                quantization_intervals = exe_params->intvCapacity;
1045        size_t i=0,j=0,k=0, I = 0, J = 0, K = 0;
1046        int reqLength;
1047        double realPrecision = pwrErrBound[0];         
1048        double pred1D, pred2D, pred3D;
1049        double diff = 0.0;
1050        double itvNum = 0;
1051        double *P0, *P1;
1052
1053        size_t r23 = r2*r3;
1054        size_t R23 = R2*R3;
1055        P0 = (double*)malloc(r23*sizeof(double));
1056        P1 = (double*)malloc(r23*sizeof(double));
1057        double radius = fabs(max)<fabs(min)?fabs(min):fabs(max);
1058        double medianValue = 0;
1059        short radExpo = getExponent_double(radius);
1060        int updateReqLength = 0;
1061        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);     
1062
1063        int* type = (int*) malloc(dataLength*sizeof(int));
1064        //type[dataLength]=0;
1065
1066        double* spaceFillingValue = oriData; //
1067       
1068        DynamicByteArray *resiBitLengthArray;
1069        new_DBA(&resiBitLengthArray, DynArrayInitLen);
1070
1071        DynamicIntArray *exactLeadNumArray;
1072        new_DIA(&exactLeadNumArray, DynArrayInitLen);
1073
1074        DynamicByteArray *exactMidByteArray;
1075        new_DBA(&exactMidByteArray, DynArrayInitLen);
1076
1077        DynamicIntArray *resiBitArray;
1078        new_DIA(&resiBitArray, DynArrayInitLen);
1079
1080        type[0] = 0;
1081
1082        unsigned char preDataBytes[8];
1083        longToBytes_bigEndian(preDataBytes, 0);
1084
1085        int reqBytesLength = reqLength/8;
1086        int resiBitsLength = reqLength%8;
1087
1088        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
1089        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
1090
1091
1092        ///////////////////////////     Process layer-0 ///////////////////////////
1093        /* Process Row-0 data 0*/
1094        type[0] = 0;
1095        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1096        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1097        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1098        memcpy(preDataBytes,vce->curBytes,8);
1099        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1100        P1[0] = vce->data;
1101
1102        /* Process Row-0 data 1*/
1103        pred1D = P1[0];
1104        diff = spaceFillingValue[1] - pred1D;
1105
1106        itvNum = fabs(diff)/realPrecision + 1;
1107
1108        if (itvNum < exe_params->intvCapacity)
1109        {
1110                if (diff < 0) itvNum = -itvNum;
1111                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
1112                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
1113        }
1114        else
1115        {
1116                if(updateReqLength==0)
1117                {
1118                        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1119                        reqBytesLength = reqLength/8;
1120                        resiBitsLength = reqLength%8;
1121                        updateReqLength = 1;
1122                }               
1123               
1124                type[1] = 0;
1125
1126                addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1127                compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1128                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1129                memcpy(preDataBytes,vce->curBytes,8);
1130                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1131                P1[1] = vce->data;
1132        }
1133
1134    /* Process Row-0 data 2 --> data r3-1 */
1135        for (j = 2; j < r3; j++)
1136        {
1137                if(j%blockEdgeSize==0)
1138                {
1139                        J++;
1140                        realPrecision = pwrErrBound[J];
1141                        updateReqLength = 0;
1142                }               
1143                pred1D = 2*P1[j-1] - P1[j-2];
1144                diff = spaceFillingValue[j] - pred1D;
1145
1146                itvNum = fabs(diff)/realPrecision + 1;
1147
1148                if (itvNum < exe_params->intvCapacity)
1149                {
1150                        if (diff < 0) itvNum = -itvNum;
1151                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
1152                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
1153                }
1154                else
1155                {
1156                        if(updateReqLength==0)
1157                        {
1158                                computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1159                                reqBytesLength = reqLength/8;
1160                                resiBitsLength = reqLength%8;
1161                                updateReqLength = 1;
1162                        }                       
1163
1164                        type[j] = 0;
1165
1166                        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1167                        compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1168                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1169                        memcpy(preDataBytes,vce->curBytes,8);
1170                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1171                        P1[j] = vce->data;
1172                }
1173        }
1174
1175        /* Process Row-1 --> Row-r2-1 */
1176        size_t index;
1177        K = 0;
1178        for (i = 1; i < r2; i++)
1179        {
1180                /* Process row-i data 0 */
1181                index = i*r3;   
1182
1183                J = 0;
1184                if(i%blockEdgeSize==0)
1185                        I++;
1186                realPrecision = pwrErrBound[I*R3+J]; //J==0
1187                updateReqLength = 0;
1188
1189                pred1D = P1[index-r3];
1190                diff = spaceFillingValue[index] - pred1D;
1191
1192                itvNum = fabs(diff)/realPrecision + 1;
1193
1194                if (itvNum < exe_params->intvCapacity)
1195                {
1196                        if (diff < 0) itvNum = -itvNum;
1197                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1198                        P1[index] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1199                }
1200                else
1201                {
1202                        if(updateReqLength==0)
1203                        {
1204                                computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1205                                reqBytesLength = reqLength/8;
1206                                resiBitsLength = reqLength%8;
1207                                updateReqLength = 1;
1208                        }               
1209                                               
1210                        type[index] = 0;
1211
1212                        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1213                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1214                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1215                        memcpy(preDataBytes,vce->curBytes,8);
1216                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1217                        P1[index] = vce->data;
1218                }
1219
1220                /* Process row-i data 1 --> data r3-1*/
1221                for (j = 1; j < r3; j++) //note that this j refers to fastest dimension (lowest order)
1222                {
1223                        index = i*r3+j;         
1224                        if(j%blockEdgeSize==0)
1225                        {
1226                                J++;
1227                                realPrecision = pwrErrBound[I*R3+J];
1228                                updateReqLength = 0;
1229                        }                       
1230               
1231                        pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];
1232
1233                        diff = spaceFillingValue[index] - pred2D;
1234
1235                        itvNum = fabs(diff)/realPrecision + 1;
1236
1237                        if (itvNum < exe_params->intvCapacity)
1238                        {
1239                                if (diff < 0) itvNum = -itvNum;
1240                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1241                                P1[index] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1242                        }
1243                        else
1244                        {
1245                                if(updateReqLength==0)
1246                                {
1247                                        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1248                                        reqBytesLength = reqLength/8;
1249                                        resiBitsLength = reqLength%8;
1250                                        updateReqLength = 1;
1251                                }                                               
1252                               
1253                                type[index] = 0;
1254
1255                                addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1256                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1257                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1258                                memcpy(preDataBytes,vce->curBytes,8);
1259                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1260                                P1[index] = vce->data;
1261                        }
1262                }
1263        }
1264
1265
1266        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
1267
1268        for (k = 1; k < r1; k++)
1269        {
1270                /* Process Row-0 data 0*/
1271                index = k*r23;                 
1272                I = 0;
1273                J = 0;
1274                if(k%blockEdgeSize==0)
1275                        K++;
1276                realPrecision = pwrErrBound[K*R23]; //J==0
1277                updateReqLength = 0;
1278               
1279                pred1D = P1[0];
1280                diff = spaceFillingValue[index] - pred1D;
1281
1282                itvNum = fabs(diff)/realPrecision + 1;
1283
1284                if (itvNum < exe_params->intvCapacity)
1285                {
1286                        if (diff < 0) itvNum = -itvNum;
1287                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1288                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1289                }
1290                else
1291                {
1292                        if(updateReqLength==0)
1293                        {
1294                                computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1295                                reqBytesLength = reqLength/8;
1296                                resiBitsLength = reqLength%8;
1297                                updateReqLength = 1;
1298                        }                                       
1299                       
1300                        type[index] = 0;
1301
1302                        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1303                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1304                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1305                        memcpy(preDataBytes,vce->curBytes,8);
1306                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1307                        P0[0] = vce->data;
1308                }
1309
1310
1311            /* Process Row-0 data 1 --> data r3-1 */
1312                for (j = 1; j < r3; j++)
1313                {
1314                        index = k*r23+j;       
1315
1316                        if(j%blockEdgeSize==0)
1317                        {
1318                                J++;
1319                                realPrecision = pwrErrBound[K*R23+J];
1320                                updateReqLength = 0;                   
1321                        }                                       
1322                        pred2D = P0[j-1] + P1[j] - P1[j-1];
1323                        diff = spaceFillingValue[index] - pred2D;
1324
1325                        itvNum = fabs(diff)/realPrecision + 1;
1326
1327                        if (itvNum < exe_params->intvCapacity)
1328                        {
1329                                if (diff < 0) itvNum = -itvNum;
1330                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1331                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1332                        }
1333                        else
1334                        {
1335                                if(updateReqLength==0)
1336                                {
1337                                        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1338                                        reqBytesLength = reqLength/8;
1339                                        resiBitsLength = reqLength%8;
1340                                        updateReqLength = 1;
1341                                }                                               
1342                               
1343                                type[index] = 0;
1344
1345                                addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1346                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1347                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1348                                memcpy(preDataBytes,vce->curBytes,8);
1349                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1350                                P0[j] = vce->data;
1351                        }
1352                }
1353
1354            /* Process Row-1 --> Row-r2-1 */
1355                size_t index2D;
1356                for (i = 1; i < r2; i++)
1357                {
1358                        /* Process Row-i data 0 */
1359                        index = k*r23 + i*r3;
1360
1361                        J = 0;
1362                        if(i%blockEdgeSize==0)
1363                                I++;
1364                        realPrecision = pwrErrBound[K*R23+I*R3+J]; //J==0
1365                        updateReqLength = 0;                   
1366                       
1367                        index2D = i*r3;         
1368                        pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
1369                        diff = spaceFillingValue[index] - pred2D;
1370
1371                        itvNum = fabs(diff)/realPrecision + 1;
1372
1373                        if (itvNum < exe_params->intvCapacity)
1374                        {
1375                                if (diff < 0) itvNum = -itvNum;
1376                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1377                                P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1378                        }
1379                        else
1380                        {
1381                                if(updateReqLength==0)
1382                                {
1383                                        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1384                                        reqBytesLength = reqLength/8;
1385                                        resiBitsLength = reqLength%8;
1386                                        updateReqLength = 1;
1387                                }                                               
1388                               
1389                                type[index] = 0;
1390
1391                                addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1392                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1393                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1394                                memcpy(preDataBytes,vce->curBytes,8);
1395                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1396                                P0[index2D] = vce->data;
1397                        }
1398
1399                        /* Process Row-i data 1 --> data r3-1 */
1400                        for (j = 1; j < r3; j++)
1401                        {
1402                                index = k*r23 + i*r3 + j;
1403                               
1404                                if(j%blockEdgeSize==0)
1405                                {
1406                                        J++;
1407                                        realPrecision = pwrErrBound[K*R23+I*R3+J];
1408                                        updateReqLength = 0;                   
1409                                }                                                       
1410                                index2D = i*r3 + j;
1411                                pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
1412                                diff = spaceFillingValue[index] - pred3D;
1413
1414                                itvNum = fabs(diff)/realPrecision + 1;
1415
1416                                if (itvNum < exe_params->intvCapacity)
1417                                {
1418                                        if (diff < 0) itvNum = -itvNum;
1419                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1420                                        P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1421                                }
1422                                else
1423                                {
1424                                        if(updateReqLength==0)
1425                                        {
1426                                                computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1427                                                reqBytesLength = reqLength/8;
1428                                                resiBitsLength = reqLength%8;
1429                                                updateReqLength = 1;
1430                                        }                                                       
1431                                       
1432                                        type[index] = 0;
1433
1434                                        addDBA_Data(resiBitLengthArray, (unsigned char)resiBitsLength);
1435                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1436                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1437                                        memcpy(preDataBytes,vce->curBytes,8);
1438                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1439                                        P0[index2D] = vce->data;
1440                                }
1441                        }
1442                }
1443
1444                double *Pt;
1445                Pt = P1;
1446                P1 = P0;
1447                P0 = Pt;
1448        }
1449        if(r23!=1)
1450                free(P0);
1451        free(P1);
1452        int exactDataNum = exactLeadNumArray->size;
1453
1454        TightDataPointStorageD* tdps;
1455
1456        new_TightDataPointStorageD2(&tdps, dataLength, exactDataNum,
1457                        type, exactMidByteArray->array, exactMidByteArray->size,
1458                        exactLeadNumArray->array,
1459                        resiBitArray->array, resiBitArray->size,
1460                        resiBitLengthArray->array, resiBitLengthArray->size, 
1461                        realPrecision, medianValue, (char)reqLength, quantization_intervals, pwrErrBoundBytes, pwrErrBoundBytes_size, radExpo);
1462
1463        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
1464
1465        //free memory
1466        free_DBA(resiBitLengthArray);
1467        free_DIA(exactLeadNumArray);
1468        free_DIA(resiBitArray);
1469        free(type);
1470
1471        free(pwrErrBound);
1472
1473        free(vce);
1474        free(lce);
1475        free_TightDataPointStorageD(tdps);
1476        free(exactMidByteArray);
1477}
1478
1479void createRangeGroups_double(double** posGroups, double** negGroups, int** posFlags, int** negFlags)
1480{
1481        size_t size = GROUP_COUNT*sizeof(double);
1482        size_t size2 = GROUP_COUNT*sizeof(int);
1483        *posGroups = (double*)malloc(size);
1484        *negGroups = (double*)malloc(size);
1485        *posFlags = (int*)malloc(size2);
1486        *negFlags = (int*)malloc(size2);
1487        memset(*posGroups, 0, size);
1488        memset(*negGroups, 0, size);
1489        memset(*posFlags, 0, size2);
1490        memset(*negFlags, 0, size2);
1491}
1492
1493void compressGroupIDArray_double(char* groupID, TightDataPointStorageD* tdps)
1494{
1495        size_t dataLength = tdps->dataSeriesLength;
1496        int* standGroupID = (int*)malloc(dataLength*sizeof(int));
1497
1498        size_t i;
1499        standGroupID[0] = groupID[0]+GROUP_COUNT; //plus an offset such that it would not be a negative number.
1500        char lastGroupIDValue = groupID[0], curGroupIDValue;
1501        int offset = 2*(GROUP_COUNT + 2);
1502        for(i=1; i<dataLength;i++)
1503        {
1504                curGroupIDValue = groupID[i];
1505                standGroupID[i] = (curGroupIDValue - lastGroupIDValue) + offset; 
1506                lastGroupIDValue = curGroupIDValue;
1507        }
1508       
1509        unsigned char* out = NULL;
1510        size_t outSize;
1511       
1512        HuffmanTree* huffmanTree = SZ_Reset();
1513        encode_withTree(huffmanTree, standGroupID, dataLength, &out, &outSize);
1514        SZ_ReleaseHuffman(huffmanTree);
1515       
1516        tdps->pwrErrBoundBytes = out; //groupIDArray
1517        tdps->pwrErrBoundBytes_size = outSize;
1518       
1519        free(standGroupID);
1520}
1521
1522TightDataPointStorageD* SZ_compress_double_1D_MDQ_pwrGroup(double* oriData, size_t dataLength, int errBoundMode, 
1523double absErrBound, double relBoundRatio, double pwrErrRatio, double valueRangeSize, double medianValue_f)
1524{
1525        size_t i;
1526        double *posGroups, *negGroups, *groups;
1527        double pos_01_group = 0, neg_01_group = 0; //[0,1] and [-1,0]
1528        int *posFlags, *negFlags, *flags;
1529        int pos_01_flag = 0, neg_01_flag = 0;
1530        createRangeGroups_double(&posGroups, &negGroups, &posFlags, &negFlags);
1531        size_t nbBins = (size_t)(1/pwrErrRatio);
1532        if(nbBins%2==1)
1533                nbBins++;
1534        exe_params->intvRadius = nbBins;
1535
1536        int reqLength, status;
1537        double medianValue = medianValue_f;
1538        double realPrecision = (double)getRealPrecision_double(valueRangeSize, errBoundMode, absErrBound, relBoundRatio, &status);
1539        if(realPrecision<0)
1540                realPrecision = pwrErrRatio;
1541        double realGroupPrecision; //precision (error) based on group ID
1542        getPrecisionReqLength_double(realPrecision);
1543        short radExpo = getExponent_double(valueRangeSize/2);
1544        short lastGroupNum = 0, groupNum, grpNum = 0;
1545       
1546        double* groupErrorBounds = generateGroupErrBounds(errBoundMode, realPrecision, pwrErrRatio);
1547        exe_params->intvRadius = generateGroupMaxIntervalCount(groupErrorBounds);
1548       
1549        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);
1550
1551        int* type = (int*) malloc(dataLength*sizeof(int));
1552        char *groupID = (char*) malloc(dataLength*sizeof(char));
1553        char *gp = groupID;
1554               
1555        double* spaceFillingValue = oriData; 
1556       
1557        DynamicIntArray *exactLeadNumArray;
1558        new_DIA(&exactLeadNumArray, DynArrayInitLen);
1559       
1560        DynamicByteArray *exactMidByteArray;
1561        new_DBA(&exactMidByteArray, DynArrayInitLen);
1562       
1563        DynamicIntArray *resiBitArray;
1564        new_DIA(&resiBitArray, DynArrayInitLen);
1565       
1566        unsigned char preDataBytes[8];
1567        intToBytes_bigEndian(preDataBytes, 0);
1568       
1569        int reqBytesLength = reqLength/8;
1570        int resiBitsLength = reqLength%8;
1571
1572        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
1573        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
1574                       
1575        int state;
1576        double curData, decValue;
1577        double pred;
1578        double predAbsErr;
1579        double interval = 0;
1580       
1581        //add the first data   
1582        type[0] = 0;
1583        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1584        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1585        memcpy(preDataBytes,vce->curBytes,8);
1586        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1587       
1588        curData = spaceFillingValue[0];
1589        groupNum = computeGroupNum_double(vce->data);
1590       
1591        if(curData > 0 && groupNum >= 0)
1592        {
1593                groups = posGroups;
1594                flags = posFlags;
1595                grpNum = groupNum;
1596        }
1597        else if(curData < 0 && groupNum >= 0)
1598        {
1599                groups = negGroups;
1600                flags = negFlags;
1601                grpNum = groupNum;
1602        }
1603        else if(curData >= 0 && groupNum == -1)
1604        {
1605                groups = &pos_01_group;
1606                flags = &pos_01_flag;
1607                grpNum = 0;
1608        }
1609        else //curData < 0 && groupNum == -1
1610        {
1611                groups = &neg_01_group;
1612                flags = &neg_01_flag;
1613                grpNum = 0;
1614        }
1615               
1616        listAdd_double_group(groups, flags, groupNum, spaceFillingValue[0], vce->data, gp);
1617        gp++;
1618       
1619        for(i=1;i<dataLength;i++)
1620        {
1621                curData = oriData[i];
1622                //printf("i=%d, posGroups[3]=%f, negGroups[3]=%f\n", i, posGroups[3], negGroups[3]);
1623               
1624                groupNum = computeGroupNum_double(curData);
1625               
1626                if(curData > 0 && groupNum >= 0)
1627                {
1628                        groups = posGroups;
1629                        flags = posFlags;
1630                        grpNum = groupNum;
1631                }
1632                else if(curData < 0 && groupNum >= 0)
1633                {
1634                        groups = negGroups;
1635                        flags = negFlags;
1636                        grpNum = groupNum;
1637                }
1638                else if(curData >= 0 && groupNum == -1)
1639                {
1640                        groups = &pos_01_group;
1641                        flags = &pos_01_flag;
1642                        grpNum = 0;
1643                }
1644                else //curData < 0 && groupNum == -1
1645                {
1646                        groups = &neg_01_group;
1647                        flags = &neg_01_flag;
1648                        grpNum = 0;
1649                }
1650
1651                if(groupNum>=GROUP_COUNT)
1652                {
1653                        type[i] = 0;
1654                        compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1655                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1656                        memcpy(preDataBytes,vce->curBytes,8);
1657                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1658                        listAdd_double_group(groups, flags, lastGroupNum, curData, vce->data, gp);      //set the group number to be last one in order to get the groupID array as smooth as possible.         
1659                }
1660                else if(flags[grpNum]==0) //the dec value may not be in the same group
1661                {       
1662                        type[i] = 0;
1663                        compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1664                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1665                        memcpy(preDataBytes,vce->curBytes,8);
1666                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1667                        //decGroupNum = computeGroupNum_double(vce->data);
1668                       
1669                        //if(decGroupNum < groupNum)
1670                        //      decValue = curData>0?pow(2, groupNum):-pow(2, groupNum);
1671                        //else if(decGroupNum > groupNum)
1672                        //      decValue = curData>0?pow(2, groupNum+1):-pow(2, groupNum+1);
1673                        //else
1674                        //      decValue = vce->data;
1675                       
1676                        decValue = vce->data;   
1677                        listAdd_double_group(groups, flags, groupNum, curData, decValue, gp);
1678                        lastGroupNum = curData>0?groupNum + 2: -(groupNum+2);
1679                }
1680                else //if flags[groupNum]==1, the dec value must be in the same group
1681                {
1682                        pred = groups[grpNum];
1683                        predAbsErr = fabs(curData - pred);
1684                        realGroupPrecision = groupErrorBounds[grpNum]; //compute real error bound
1685                        interval = realGroupPrecision*2;
1686                        state = (predAbsErr/realGroupPrecision+1)/2;
1687                        if(curData>=pred)
1688                        {
1689                                type[i] = exe_params->intvRadius+state;
1690                                decValue = pred + state*interval;
1691                        }
1692                        else //curData<pred
1693                        {
1694                                type[i] = exe_params->intvRadius-state;
1695                                decValue = pred - state*interval;
1696                        }
1697                        //decGroupNum = computeGroupNum_double(pred);
1698                       
1699                        if((decValue>0&&curData<0)||(decValue<0&&curData>=0))
1700                                decValue = 0;
1701                        //else
1702                        //{
1703                        //      if(decGroupNum < groupNum)
1704                        //              decValue = curData>0?pow(2, groupNum):-pow(2, groupNum);
1705                        //      else if(decGroupNum > groupNum)
1706                        //              decValue = curData>0?pow(2, groupNum+1):-pow(2, groupNum+1);
1707                        //      else
1708                        //              decValue = pred;                               
1709                        //}
1710                       
1711                        if(fabs(curData-decValue)>realGroupPrecision)
1712                        {       
1713                                type[i] = 0;
1714                                compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
1715                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
1716                                memcpy(preDataBytes,vce->curBytes,8);
1717                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
1718
1719                                decValue = vce->data;   
1720                        }
1721                       
1722                        listAdd_double_group(groups, flags, groupNum, curData, decValue, gp);                   
1723                        lastGroupNum = curData>=0?groupNum + 2: -(groupNum+2);                 
1724                }
1725                gp++;   
1726
1727        }
1728       
1729        int exactDataNum = exactLeadNumArray->size;
1730       
1731        TightDataPointStorageD* tdps;
1732                       
1733        //combineTypeAndGroupIDArray(nbBins, dataLength, &type, groupID);
1734
1735        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum, 
1736                        type, exactMidByteArray->array, exactMidByteArray->size, 
1737                        exactLeadNumArray->array, 
1738                        resiBitArray->array, resiBitArray->size, 
1739                        resiBitsLength, 
1740                        realPrecision, medianValue, (char)reqLength, nbBins, NULL, 0, radExpo); 
1741       
1742        compressGroupIDArray_double(groupID, tdps);
1743       
1744        free(posGroups);
1745        free(negGroups);
1746        free(posFlags);
1747        free(negFlags);
1748        free(groupID);
1749        free(groupErrorBounds);
1750       
1751        free_DIA(exactLeadNumArray);
1752        free_DIA(resiBitArray);
1753        free(type);     
1754        free(vce);
1755        free(lce);     
1756        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);     
1757       
1758        return tdps;
1759}
1760
1761void SZ_compress_args_double_NoCkRngeNoGzip_1D_pwrgroup(unsigned char** newByteData, double *oriData,
1762size_t dataLength, double absErrBound, double relBoundRatio, double pwrErrRatio, double valueRangeSize, double medianValue_f, size_t *outSize)
1763{
1764        TightDataPointStorageD* tdps = SZ_compress_double_1D_MDQ_pwrGroup(oriData, dataLength, confparams_cpr->errorBoundMode, 
1765        absErrBound, relBoundRatio, pwrErrRatio, 
1766        valueRangeSize, medianValue_f);
1767
1768        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
1769
1770        if(*outSize>dataLength*sizeof(double))
1771                SZ_compress_args_double_StoreOriData(oriData, dataLength+2, tdps, newByteData, outSize);
1772
1773        free_TightDataPointStorageD(tdps);
1774}
Note: See TracBrowser for help on using the repository browser.