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