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

Revision 2c47b73, 6.4 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.c
3 *  @author Sheng Di and Dingwen Tao
4 *  @date Aug, 2016
5 *  @brief SZ_Init, Compression and Decompression functions
6 *  (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7 *      See COPYRIGHT in top-level directory.
8 */
9
10
11#include <stdio.h>
12#include <stdlib.h>
13#include <string.h>
14#include <unistd.h>
15#include <math.h>
16#include "sz.h"
17#include "CompressElement.h"
18#include "DynamicByteArray.h"
19#include "DynamicIntArray.h"
20#include "TightDataPointStorageF.h"
21#include "zlib.h"
22#include "rw.h"
23#include "sz_float_ts.h"
24
25unsigned int optimize_intervals_float_1D_ts(float *oriData, size_t dataLength, float* preData, double realPrecision)
26{       
27        size_t i = 0, radiusIndex;
28        float pred_value = 0, pred_err;
29        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
30        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
31        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
32        for(i=2;i<dataLength;i++)
33        {
34                if(i%confparams_cpr->sampleDistance==0)
35                {
36                        pred_value = preData[i];
37                        pred_err = fabs(pred_value - oriData[i]);
38                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
39                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
40                                radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
41                        intervals[radiusIndex]++;
42                }
43        }
44        //compute the appropriate number
45        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
46        size_t sum = 0;
47        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
48        {
49                sum += intervals[i];
50                if(sum>targetCount)
51                        break;
52        }
53        if(i>=confparams_cpr->maxRangeRadius)
54                i = confparams_cpr->maxRangeRadius-1;
55               
56        unsigned int accIntervals = 2*(i+1);
57        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
58       
59        if(powerOf2<32)
60                powerOf2 = 32;
61       
62        free(intervals);
63        return powerOf2;
64}
65
66TightDataPointStorageF* SZ_compress_float_1D_MDQ_ts(float *oriData, size_t dataLength, sz_multisteps* multisteps,
67double realPrecision, float valueRangeSize, float medianValue_f)
68{
69        float* preStepData = (float*)(multisteps->hist_data);
70
71        //store the decompressed data
72        float* decData = (float*)malloc(sizeof(float)*dataLength);
73        memset(decData, 0, sizeof(float)*dataLength);
74       
75        unsigned int quantization_intervals;
76        if(exe_params->optQuantMode==1)
77                quantization_intervals = optimize_intervals_float_1D_ts(oriData, dataLength, preStepData, realPrecision);
78        else
79                quantization_intervals = exe_params->intvCapacity;
80        updateQuantizationInfo(quantization_intervals); 
81
82        size_t i;
83        int reqLength;
84        float medianValue = medianValue_f;
85        short radExpo = getExponent_float(valueRangeSize/2);
86       
87        computeReqLength_float(realPrecision, radExpo, &reqLength, &medianValue);       
88
89        int* type = (int*) malloc(dataLength*sizeof(int));
90               
91        float* spaceFillingValue = oriData; //
92       
93        DynamicIntArray *exactLeadNumArray;
94        new_DIA(&exactLeadNumArray, DynArrayInitLen);
95       
96        DynamicByteArray *exactMidByteArray;
97        new_DBA(&exactMidByteArray, DynArrayInitLen);
98       
99        DynamicIntArray *resiBitArray;
100        new_DIA(&resiBitArray, DynArrayInitLen);
101       
102        unsigned char preDataBytes[4];
103        intToBytes_bigEndian(preDataBytes, 0);
104       
105        int reqBytesLength = reqLength/8;
106        int resiBitsLength = reqLength%8;
107
108        FloatValueCompressElement *vce = (FloatValueCompressElement*)malloc(sizeof(FloatValueCompressElement));
109        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
110                               
111        //add the first data   
112        type[0] = 0;
113        compressSingleFloatValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
114        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
115        memcpy(preDataBytes,vce->curBytes,4);
116        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
117        decData[0] = vce->data;
118               
119        //add the second data
120        type[1] = 0;
121        compressSingleFloatValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
122        updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
123        memcpy(preDataBytes,vce->curBytes,4);
124        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
125        decData[1] = vce->data; 
126       
127        int state = 0;
128        double checkRadius = 0;
129        float curData = 0;
130        float pred = 0;
131        float predAbsErr = 0;
132        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
133        double interval = 2*realPrecision;
134       
135        for(i=2;i<dataLength;i++)
136        {
137                curData = spaceFillingValue[i];
138                pred = preStepData[i];
139                predAbsErr = fabs(curData - pred);     
140                if(predAbsErr<=checkRadius)
141                {
142                        state = (predAbsErr/realPrecision+1)/2;
143                        if(curData>=pred)
144                        {
145                                type[i] = exe_params->intvRadius+state;
146                                pred = pred + state*interval;
147                        }
148                        else //curData<pred
149                        {
150                                type[i] = exe_params->intvRadius-state;
151                                pred = pred - state*interval;
152                        }
153                               
154                        //double-check the prediction error in case of machine-epsilon impact   
155                        if(fabs(curData-pred)>realPrecision)
156                        {       
157                                type[i] = 0;                           
158                                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
159                                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
160                                memcpy(preDataBytes,vce->curBytes,4);
161                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);         
162                                decData[i] = vce->data;
163                        }
164                        else
165                        {
166                                decData[i] = pred;
167                        }
168                       
169                        continue;
170                }
171               
172                //unpredictable data processing         
173                type[i] = 0;           
174                compressSingleFloatValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
175                updateLossyCompElement_Float(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
176                memcpy(preDataBytes,vce->curBytes,4);
177                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
178                decData[i] = vce->data;
179        }//end of for
180               
181        size_t exactDataNum = exactLeadNumArray->size;
182       
183        TightDataPointStorageF* tdps;
184                       
185        new_TightDataPointStorageF(&tdps, dataLength, exactDataNum, 
186                        type, exactMidByteArray->array, exactMidByteArray->size, 
187                        exactLeadNumArray->array, 
188                        resiBitArray->array, resiBitArray->size, 
189                        resiBitsLength,
190                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
191
192        //free memory
193        free_DIA(exactLeadNumArray);
194        free_DIA(resiBitArray);
195        free(type);     
196        free(vce);
197        free(lce);     
198        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
199               
200        memcpy(preStepData, decData, dataLength*sizeof(float)); //update the data
201        free(decData);
202       
203        return tdps;
204}
205
206
Note: See TracBrowser for help on using the repository browser.