1 : #include <stdio.h>
2 : #include <stdlib.h>
3 : #include <string.h>
4 : #include <math.h>
5 : #include "tdlpack.h"
6 : #include "myerror.h"
7 : #include "meta.h"
8 : #include "memendian.h"
9 : #include "fileendian.h"
10 : #include "myassert.h"
11 : #include "myutil.h"
12 : #include "clock.h"
13 : #ifdef MEMWATCH
14 : #include "memwatch.h"
15 : #endif
16 :
17 : #define GRIB_UNSIGN_INT3(a,b,c) ((a<<16)+(b<<8)+c)
18 : #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b)
19 : #define GRIB_SIGN_INT3(a,b,c) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 127) << 16)+(b<<8)+c))
20 : #define GRIB_SIGN_INT2(a,b) ((1-(int) ((unsigned) (a & 0x80) >> 6)) * (int) (((a & 0x7f) << 8) + b))
21 :
22 : /* *INDENT-OFF* */
23 : TDLP_TableType TDLP_B_Table[5] = {
24 : /* 0 */ {0, "Continuous field"},
25 : /* 1 */ {1, "Point Binary - cumulative from above"},
26 : /* 2 */ {2, "Point Binary - cumulative from below"},
27 : /* 3 */ {3, "Point Binary - discrete"},
28 : /* 4 */ {5, "Grid Binary - "},
29 : };
30 :
31 : TDLP_TableType TDLP_DD_Table[9] = {
32 : /* 0 */ {0, "Independent of model"},
33 : /* 1 */ {6, "NGM model"},
34 : /* 2 */ {7, "Eta model"},
35 : /* 3 */ {8, "AVN model"},
36 : /* 4 */ {9, "MRF model"},
37 : /* 5 */ {79, "NDFD forecast"},
38 : /* 6 */ {80, "MOS AEV forecast"},
39 : /* 7 */ {81, "Local AEV Firecasts"},
40 : /* 8 */ {82, "Obs matching AEV Forecasts"},
41 : };
42 :
43 : TDLP_TableType TDLP_V_Table[4] = {
44 : /* 0 */ {0, "No vertical processing"},
45 : /* 1 */ {1, "Difference Levels (UUUU - LLLL)"},
46 : /* 2 */ {2, "Sum Levels (UUUU + LLLL)"},
47 : /* 3 */ {3, "Mean Levels (UUUU + LLLL) / 2."},
48 : };
49 :
50 : TDLP_TableType TDLP_T_Table[3] = {
51 : /* 0 */ {0, "No nolinear tranform"},
52 : /* 1 */ {1, "Square transform"},
53 : /* 2 */ {2, "Square root tranform"},
54 : };
55 :
56 : TDLP_TableType TDLP_Oper_Table[9] = {
57 : /* 0 */ {0, "No time operator"},
58 : /* 1 */ {1, "Mean (Var 1, Var 2)"},
59 : /* 2 */ {2, "Difference (Var 1 - Var 2)"},
60 : /* 3 */ {3, "Maximum (Var 1, Var 2)"},
61 : /* 4 */ {4, "Minimum (Var 1, Var 2)"},
62 : /* 5 */ {5, "Mean (Var 1, Var 2)"},
63 : /* 6 */ {6, "Difference (Var 2 - Var 1)"},
64 : /* 7 */ {7, "Maximum (Var 1, Var 2)"},
65 : /* 8 */ {8, "Minimum (Var 1, Var 2)"},
66 : };
67 :
68 : TDLP_TableType TDLP_I_Table[4] = {
69 : /* 0 */ {0, "No interpolation"},
70 : /* 1 */ {1, "Bi-quadratic interpolation"},
71 : /* 2 */ {2, "Bi-linear interpolation"},
72 : /* 3 */ {3, "Special interpolation for QPF"},
73 : };
74 :
75 : TDLP_TableType TDLP_S_Table[6] = {
76 : /* 0 */ {0, "No smoothing"},
77 : /* 1 */ {1, "5-point smoothing"},
78 : /* 2 */ {2, "9-point smoothing"},
79 : /* 3 */ {3, "25-point smoothing"},
80 : /* 4 */ {4, "81-point smoothing"},
81 : /* 5 */ {5, "168-point smoothing"},
82 : };
83 : /* *INDENT-ON* */
84 :
85 : typedef struct {
86 : sInt4 min; /* Minimum value in a group. Can cause problems if
87 : * this is unsigned. */
88 : uChar bit; /* # of bits in a group. 2^31 is the largest # of
89 : * bits that can be used to hold # of bits in a
90 : * group. However, the # of bits for a number can't
91 : * be > 64, and is probably < 32, so bit is < 6 and
92 : * probably < 5. */
93 : uInt4 num; /* number of values in the group. May need this to be
94 : * signed. */
95 : sInt4 max; /* Max value for the group */
96 : uInt4 start; /* index in Data where group starts. */
97 : uChar f_trySplit; /* Flag, if we have tried to split this group. */
98 : uChar f_tryShift; /* Flag, if we have tried to shift this group. */
99 : } TDLGroupType;
100 :
101 : /*****************************************************************************
102 : * ReadTDLPSect1() --
103 : *
104 : * Arthur Taylor / MDL
105 : *
106 : * PURPOSE
107 : * Parses the TDLP "Product Definition Section" or section 1, filling out
108 : * the pdsMeta data structure.
109 : *
110 : * ARGUMENTS
111 : * pds = The compressed part of the message dealing with "PDS". (Input)
112 : * gribLen = The total length of the TDLP message. (Input)
113 : * curLoc = Current location in the TDLP message. (Output)
114 : * pdsMeta = The filled out pdsMeta data structure. (Output)
115 : * f_gds = boolean if there is a Grid Definition Section. (Output)
116 : * f_bms = boolean if there is a Bitmap Section. (Output)
117 : * DSF = Decimal Scale Factor for unpacking the data. (Output)
118 : * BSF = Binary Scale Factor for unpacking the data. (Output)
119 : *
120 : * FILES/DATABASES: None
121 : *
122 : * RETURNS: int (could use errSprintf())
123 : * 0 = OK
124 : * -1 = gribLen is too small.
125 : *
126 : * HISTORY
127 : * 10/2004 Arthur Taylor (MDL): Created.
128 : *
129 : * NOTES
130 : *****************************************************************************
131 : */
132 0 : static int ReadTDLPSect1 (uChar *pds, sInt4 tdlpLen, sInt4 *curLoc,
133 : pdsTDLPType * pdsMeta, char *f_gds, char *f_bms,
134 : short int *DSF, short int *BSF)
135 : {
136 : char sectLen; /* Length of section. */
137 : sInt4 li_temp; /* Temporary variable. */
138 : int W, XXXX, YY; /* Helps with Threshold calculation. */
139 : int year, t_year; /* The reference year, and a consistency test */
140 : uChar month, t_month; /* The reference month, and a consistency test */
141 : uChar day, t_day; /* The reference day, and a consistency test */
142 : uChar hour, t_hour; /* The reference hour, and a consistency test */
143 : uChar min; /* The reference minute */
144 : uShort2 project_hr; /* The projection in hours. */
145 : sInt4 tau; /* Used to cross check project_hr */
146 : int lenPL; /* Length of the Plain Language descriptor. */
147 :
148 0 : sectLen = *(pds++);
149 0 : *curLoc += sectLen;
150 0 : if (*curLoc > tdlpLen) {
151 0 : errSprintf ("Ran out of data in PDS (TDLP Section 1)\n");
152 0 : return -1;
153 : }
154 0 : myAssert (sectLen <= 71);
155 0 : if (sectLen < 39) {
156 0 : errSprintf ("TDLP Section 1 is too small.\n");
157 0 : return -1;
158 : }
159 0 : *f_bms = (GRIB2BIT_7 & *pds) ? 1 : 0;
160 0 : *f_gds = (GRIB2BIT_8 & *pds) ? 1 : 0;
161 0 : pds++;
162 0 : year = GRIB_UNSIGN_INT2 (*pds, pds[1]);
163 0 : pds += 2;
164 0 : month = *(pds++);
165 0 : day = *(pds++);
166 0 : hour = *(pds++);
167 0 : min = *(pds++);
168 0 : MEMCPY_BIG (&li_temp, pds, sizeof (sInt4));
169 0 : pds += 4;
170 0 : t_year = li_temp / 1000000L;
171 0 : li_temp -= t_year * 1000000L;
172 0 : t_month = li_temp / 10000L;
173 0 : li_temp -= t_month * 10000L;
174 0 : t_day = li_temp / 100;
175 0 : t_hour = li_temp - t_day * 100;
176 0 : if ((t_year != year) || (t_month != month) || (t_day != day) ||
177 : (t_hour != hour)) {
178 0 : errSprintf ("Error Inconsistant Times in ReadTDLPSect1.\n");
179 0 : return -1;
180 : }
181 0 : if (ParseTime (&(pdsMeta->refTime), year, month, day, hour, min, 0) != 0) {
182 0 : preErrSprintf ("Error In call to ParseTime in ReadTDLPSect1.\n");
183 0 : return -1;
184 : }
185 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
186 0 : pds += 4;
187 0 : pdsMeta->ID1 = li_temp;
188 0 : pdsMeta->CCC = li_temp / 1000000L;
189 0 : li_temp -= pdsMeta->CCC * 1000000L;
190 0 : pdsMeta->FFF = li_temp / 1000;
191 0 : li_temp -= pdsMeta->FFF * 1000;
192 0 : pdsMeta->B = li_temp / 100;
193 0 : pdsMeta->DD = li_temp - pdsMeta->B * 100;
194 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
195 0 : pds += 4;
196 0 : pdsMeta->ID2 = li_temp;
197 0 : pdsMeta->V = li_temp / 100000000L;
198 0 : li_temp -= pdsMeta->V * 100000000L;
199 0 : pdsMeta->LLLL = li_temp / 10000;
200 0 : pdsMeta->UUUU = li_temp - pdsMeta->LLLL * 10000;
201 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
202 0 : pds += 4;
203 0 : pdsMeta->ID3 = li_temp;
204 0 : pdsMeta->T = li_temp / 100000000L;
205 0 : li_temp -= pdsMeta->T * 100000000L;
206 0 : pdsMeta->RR = li_temp / 1000000L;
207 0 : li_temp -= pdsMeta->RR * 1000000L;
208 0 : pdsMeta->Oper = li_temp / 100000L;
209 0 : li_temp -= pdsMeta->Oper * 100000L;
210 0 : pdsMeta->HH = li_temp / 1000;
211 0 : pdsMeta->ttt = li_temp - pdsMeta->HH * 1000;
212 0 : MEMCPY_BIG (&(li_temp), pds, sizeof (sInt4));
213 0 : pds += 4;
214 0 : pdsMeta->ID4 = li_temp;
215 0 : W = li_temp / 1000000000L;
216 0 : li_temp -= W * 1000000000L;
217 0 : XXXX = li_temp / 100000L;
218 0 : li_temp -= XXXX * 100000L;
219 0 : if (W) {
220 0 : XXXX = -1 * XXXX;
221 : }
222 0 : YY = li_temp / 1000L;
223 0 : li_temp -= YY * 1000L;
224 0 : if (YY >= 50) {
225 0 : YY = -1 * (YY - 50);
226 : }
227 0 : pdsMeta->thresh = (XXXX / 10000.) * pow (10.0, YY);
228 0 : pdsMeta->I = li_temp / 100;
229 0 : li_temp -= pdsMeta->I * 100L;
230 0 : pdsMeta->S = li_temp / 10;
231 0 : pdsMeta->G = li_temp - pdsMeta->S * 10;
232 0 : project_hr = GRIB_UNSIGN_INT2 (*pds, pds[1]);
233 0 : tau = pdsMeta->ID3 - ((pdsMeta->ID3 / 1000) * 1000);
234 0 : if (tau != project_hr) {
235 : printf ("Warning: Inconsistant Projections in hours in "
236 0 : "ReadTDLPSect1 (%d vs %d)\n", tau, project_hr);
237 : /*
238 : errSprintf ("Warning: Inconsistant Projections in hours in "
239 : "ReadTDLPSect1 (%ld vs %d)\n", tau, project_hr);
240 : */
241 0 : project_hr = tau;
242 : /* return -1; */
243 : }
244 0 : pds += 2;
245 0 : pdsMeta->project = project_hr * 3600 + (*(pds++)) * 60;
246 0 : pdsMeta->procNum = (*(pds++));
247 0 : pdsMeta->seqNum = (*(pds++));
248 0 : *DSF = (*pds > 128) ? 128 - (*(pds++)) : (*(pds++));
249 0 : *BSF = (*pds > 128) ? 128 - (*(pds++)) : (*(pds++));
250 0 : if ((*pds != 0) || (pds[1] != 0) || (pds[2] != 0)) {
251 0 : errSprintf ("Error Reserved was not set to 0 in ReadTDLPSect1.\n");
252 0 : return -1;
253 : }
254 0 : pds += 3;
255 0 : lenPL = (*(pds++));
256 0 : if (sectLen - lenPL != 39) {
257 : errSprintf ("Error sectLen(%d) - lenPL(%d) != 39 in ReadTDLPSect1.\n",
258 0 : sectLen, lenPL);
259 0 : return -1;
260 : }
261 0 : if (lenPL > 32) {
262 0 : lenPL = 32;
263 : }
264 0 : strncpy (pdsMeta->Descriptor, (char *) pds, lenPL);
265 0 : pdsMeta->Descriptor[lenPL] = '\0';
266 0 : strTrim (pdsMeta->Descriptor);
267 0 : return 0;
268 : }
269 :
270 : /*****************************************************************************
271 : * TDLP_TableLookUp() --
272 : *
273 : * Arthur Taylor / MDL
274 : *
275 : * PURPOSE
276 : * To look up TDL Ids information in a given table.
277 : *
278 : * ARGUMENTS
279 : * table = The Table to look up the Id in. (Input)
280 : * tableLen = The length of the Table. (Input)
281 : * index = The index into the Table. (Input)
282 : *
283 : * FILES/DATABASES: None
284 : *
285 : * RETURNS: char *
286 : * Returns the meta data that is associated with index in the table.
287 : *
288 : * HISTORY
289 : * 10/2004 Arthur Taylor (MDL): Created.
290 : *
291 : * NOTES
292 : *****************************************************************************
293 : */
294 0 : static const char *TDLP_TableLookUp(TDLP_TableType * table, int tableLen,
295 : int index)
296 : {
297 : int i; /* Loop counter. */
298 :
299 0 : for (i = 0; i < tableLen; i++) {
300 0 : if (table[i].index == index)
301 0 : return (table[i].data);
302 : }
303 0 : return "Unknown";
304 : }
305 :
306 : /*****************************************************************************
307 : * PrintPDS_TDLP() --
308 : *
309 : * Arthur Taylor / MDL
310 : *
311 : * PURPOSE
312 : * To generate the message for the Product Definition Sections of the TDLP
313 : * Message.
314 : *
315 : * ARGUMENTS
316 : * pds = The TDLP Product Definition Section to print. (Input)
317 : *
318 : * FILES/DATABASES: None
319 : *
320 : * RETURNS: void
321 : *
322 : * HISTORY
323 : * 10/2004 Arthur Taylor (MDL): Created.
324 : *
325 : * NOTES
326 : *****************************************************************************
327 : */
328 0 : void PrintPDS_TDLP (pdsTDLPType * pds)
329 : {
330 : char buffer[25]; /* Stores format of pds1->refTime. */
331 :
332 : /*
333 : strftime (buffer, 25, "%m/%d/%Y %H:%M:%S UTC", gmtime (&(pds->refTime)));
334 : */
335 0 : Clock_Print (buffer, 25, pds->refTime, "%m/%d/%Y %H:%M:%S UTC", 0);
336 :
337 0 : Print ("PDS-TDLP", "Reference Time", Prt_S, buffer);
338 0 : Print ("PDS-TDLP", "Plain Language", Prt_S, pds->Descriptor);
339 0 : sprintf (buffer, "%09d", pds->ID1);
340 0 : Print ("PDS-TDLP", "ID 1", Prt_S, buffer);
341 0 : sprintf (buffer, "%09d", pds->ID2);
342 0 : Print ("PDS-TDLP", "ID 2", Prt_S, buffer);
343 0 : sprintf (buffer, "%09d", pds->ID3);
344 0 : Print ("PDS-TDLP", "ID 3", Prt_S, buffer);
345 0 : Print ("PDS-TDLP", "ID 4", Prt_D, pds->ID4);
346 0 : Print ("PDS-TDLP", "Model or Process Number", Prt_D, pds->procNum);
347 0 : Print ("PDS-TDLP", "Sequence Number", Prt_D, pds->seqNum);
348 :
349 0 : sprintf (buffer, "%03d", pds->CCC);
350 0 : Print ("PDS-TDLP", "ID1-CCC", Prt_S, buffer);
351 0 : sprintf (buffer, "%03d", pds->FFF);
352 0 : Print ("PDS-TDLP", "ID1-FFF", Prt_S, buffer);
353 : Print ("PDS-TDLP", "ID1-B", Prt_DS, pds->B,
354 0 : TDLP_TableLookUp (TDLP_B_Table, sizeof (TDLP_B_Table), pds->B));
355 0 : sprintf (buffer, "%02d", pds->DD);
356 : Print ("PDS-TDLP", "ID1-DD", Prt_SS, buffer,
357 0 : TDLP_TableLookUp (TDLP_DD_Table, sizeof (TDLP_DD_Table), pds->DD));
358 :
359 : Print ("PDS-TDLP", "ID2-V", Prt_DS, pds->V,
360 0 : TDLP_TableLookUp (TDLP_V_Table, sizeof (TDLP_V_Table), pds->V));
361 0 : sprintf (buffer, "%04d", pds->LLLL);
362 0 : Print ("PDS-TDLP", "ID2-LLLL", Prt_S, buffer);
363 0 : sprintf (buffer, "%04d", pds->UUUU);
364 0 : Print ("PDS-TDLP", "ID2-UUUU", Prt_S, buffer);
365 :
366 0 : if (pds->Oper != 0) {
367 : Print ("PDS-TDLP", "ID3-T", Prt_DS, pds->T,
368 0 : TDLP_TableLookUp (TDLP_T_Table, sizeof (TDLP_T_Table), pds->T));
369 0 : sprintf (buffer, "%02d", pds->RR);
370 : Print ("PDS-TDLP", "ID3-RR", Prt_SS, buffer,
371 0 : "Run time offset in hours");
372 : Print ("PDS-TDLP", "ID3-Oper", Prt_DS, pds->Oper,
373 : TDLP_TableLookUp (TDLP_Oper_Table, sizeof (TDLP_Oper_Table),
374 0 : pds->Oper));
375 0 : sprintf (buffer, "%02d", pds->HH);
376 : Print ("PDS-TDLP", "ID3-HH", Prt_SS, buffer,
377 0 : "Number of hours between variables");
378 : } else {
379 : Print ("PDS-TDLP", "ID3-Oper", Prt_DS, pds->Oper,
380 : TDLP_TableLookUp (TDLP_Oper_Table, sizeof (TDLP_Oper_Table),
381 0 : pds->Oper));
382 : }
383 0 : sprintf (buffer, "%03d", pds->ttt);
384 0 : Print ("PDS-TDLP", "ID3-ttt", Prt_SS, buffer, "Forecast Projection");
385 :
386 0 : Print ("PDS-TDLP", "ID4-thresh", Prt_F, pds->thresh);
387 : Print ("PDS-TDLP", "ID4-I", Prt_DS, pds->I,
388 0 : TDLP_TableLookUp (TDLP_I_Table, sizeof (TDLP_I_Table), pds->I));
389 : Print ("PDS-TDLP", "ID4-S", Prt_DS, pds->S,
390 0 : TDLP_TableLookUp (TDLP_S_Table, sizeof (TDLP_S_Table), pds->S));
391 0 : Print ("PDS-TDLP", "ID4-G", Prt_D, pds->G);
392 0 : }
393 :
394 : /*****************************************************************************
395 : * TDLP_ElemSurfUnit() --
396 : *
397 : * Arthur Taylor / MDL
398 : *
399 : * PURPOSE
400 : * Deal with element name, unit, and comment for TDLP items. Also deals
401 : * with the surface information stored in ID2.
402 : *
403 : * ARGUMENTS
404 : * pds = The TDLP Product Definition Section to parse. (Input)
405 : * element = The resulting element name. (Output)
406 : * unitName = The resulting unit name. (Output)
407 : * comment = The resulting comment. (Output)
408 : * shortFstLevel = The resulting short forecast level. (Output)
409 : * longFstLevel = The resulting long forecast level. (Output)
410 : *
411 : * FILES/DATABASES: None
412 : *
413 : * RETURNS: void
414 : *
415 : * HISTORY
416 : * 10/2004 Arthur Taylor (MDL): Created.
417 : *
418 : * NOTES
419 : *****************************************************************************
420 : */
421 0 : static void TDLP_ElemSurfUnit (pdsTDLPType * pds, char **element,
422 : char **unitName, char **comment,
423 : char **shortFstLevel, char **longFstLevel)
424 : {
425 : char *ptr; /* Help guess unitName, and clean the elemName. */
426 : char *ptr2; /* Help guess unitName, and clean the elemName. */
427 :
428 0 : myAssert (*element == NULL);
429 0 : myAssert (*unitName == NULL);
430 0 : myAssert (*comment == NULL);
431 0 : myAssert (*shortFstLevel == NULL);
432 0 : myAssert (*longFstLevel == NULL);
433 :
434 0 : *element = (char *) malloc (1 + strlen (pds->Descriptor) * sizeof (char));
435 0 : strcpy (*element, pds->Descriptor);
436 0 : (*element)[strlen (pds->Descriptor)] = '\0';
437 :
438 0 : ptr = strchr (*element, '(');
439 0 : if (ptr != NULL) {
440 0 : ptr2 = strchr (ptr, ')');
441 0 : *ptr2 = '\0';
442 0 : if (strcmp (ptr + 1, "unofficial id") == 0) {
443 0 : *unitName = (char *) malloc ((1 + 3) * sizeof (char));
444 0 : strcpy (*unitName, "[-]");
445 : } else {
446 0 : reallocSprintf (unitName, "[%s]", ptr + 1);
447 : }
448 : /* Trim any parens from element. */
449 0 : *ptr = '\0';
450 0 : strTrimRight (*element, ' ');
451 : } else {
452 0 : *unitName = (char *) malloc ((1 + 3) * sizeof (char));
453 0 : strcpy (*unitName, "[-]");
454 : }
455 0 : ptr = *element;
456 0 : while (*ptr != '\0') {
457 0 : if (*ptr == ' ') {
458 0 : *ptr = '-';
459 : }
460 0 : ptr++;
461 : }
462 0 : strCompact (*element, '-');
463 :
464 : reallocSprintf (comment, "%09ld-%09ld-%09ld-%ld %s", pds->ID1,
465 0 : pds->ID2, pds->ID3, pds->ID4, *unitName);
466 0 : reallocSprintf (shortFstLevel, "%09ld", pds->ID2);
467 0 : reallocSprintf (longFstLevel, "%09ld", pds->ID2);
468 0 : }
469 :
470 : /*****************************************************************************
471 : * TDLP_Inventory() --
472 : *
473 : * Arthur Taylor / MDL
474 : *
475 : * PURPOSE
476 : * Parses the TDLP "Product Definition Section" for enough information to
477 : * fill out the inventory data structure so we can do a simple inventory on
478 : * the file in a similar way to how we did it for GRIB1 and GRIB2.
479 : *
480 : * ARGUMENTS
481 : * fp = An opened TDLP file already at the correct message. (Input)
482 : * tdlpLen = The total length of the TDLP message. (Input)
483 : * inv = The inventory data structure that we need to fill. (Output)
484 : *
485 : * FILES/DATABASES: None
486 : *
487 : * RETURNS: int (could use errSprintf())
488 : * 0 = OK
489 : * -1 = tdlpLen is too small.
490 : *
491 : * HISTORY
492 : * 10/2004 Arthur Taylor (MDL): Created.
493 : *
494 : * NOTES
495 : * Speed improvements...
496 : * 1) pds doen't need to be allocated each time.
497 : * 2) Not all data is needed, do something like TDLP_RefTime
498 : * 3) TDLP_ElemSurfUnit may be slow?
499 : *****************************************************************************
500 : */
501 0 : int TDLP_Inventory (DataSource &fp, sInt4 tdlpLen, inventoryType *inv)
502 : {
503 : sInt4 curLoc; /* Where we are in the current TDLP message. */
504 : int i_temp;
505 : uChar sectLen; /* Length of section. */
506 : uChar *pds; /* The part of the message dealing with the PDS. */
507 : pdsTDLPType pdsMeta; /* The pds parsed into a usable data structure. */
508 : char f_gds; /* flag if there is a GDS section. */
509 : char f_bms; /* flag if there is a BMS section. */
510 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
511 : short int BSF; /* Binary Scale Factor for unpacking the data. */
512 :
513 0 : curLoc = 8;
514 0 : if ((i_temp = fp.DataSourceFgetc()) == EOF) {
515 0 : errSprintf ("Ran out of file in PDS (TDLP_Inventory).\n");
516 0 : return -1;
517 : }
518 0 : sectLen = (uChar) i_temp;
519 0 : curLoc += sectLen;
520 0 : if (curLoc > tdlpLen) {
521 0 : errSprintf ("Ran out of data in PDS (TDLP_Inventory)\n");
522 0 : return -1;
523 : }
524 0 : pds = (uChar *) malloc (sectLen * sizeof (uChar));
525 0 : *pds = sectLen;
526 0 : if (fp.DataSourceFread (pds + 1, sizeof (char), sectLen - 1) + 1 != sectLen) {
527 0 : errSprintf ("Ran out of file.\n");
528 0 : free (pds);
529 0 : return -1;
530 : }
531 :
532 0 : if (ReadTDLPSect1 (pds, tdlpLen, &curLoc, &pdsMeta, &f_gds, &f_bms,
533 : &DSF, &BSF) != 0) {
534 0 : preErrSprintf ("Inside TDLP_Inventory\n");
535 0 : free (pds);
536 0 : return -1;
537 : }
538 0 : free (pds);
539 :
540 0 : inv->element = NULL;
541 0 : inv->unitName = NULL;
542 0 : inv->comment = NULL;
543 0 : free (inv->shortFstLevel);
544 0 : inv->shortFstLevel = NULL;
545 0 : free (inv->longFstLevel);
546 0 : inv->longFstLevel = NULL;
547 : TDLP_ElemSurfUnit (&pdsMeta, &(inv->element), &(inv->unitName),
548 : &(inv->comment), &(inv->shortFstLevel),
549 0 : &(inv->longFstLevel));
550 :
551 0 : inv->refTime = pdsMeta.refTime;
552 0 : inv->validTime = pdsMeta.refTime + pdsMeta.project;
553 0 : inv->foreSec = pdsMeta.project;
554 0 : return 0;
555 : }
556 :
557 : /*****************************************************************************
558 : * TDLP_RefTime() --
559 : *
560 : * Arthur Taylor / MDL
561 : *
562 : * PURPOSE
563 : * Return just the reference time of a TDLP message.
564 : *
565 : * ARGUMENTS
566 : * fp = An opened TDLP file already at the correct message. (Input)
567 : * tdlpLen = The total length of the TDLP message. (Input)
568 : * refTime = The reference time of interest. (Output)
569 : *
570 : * FILES/DATABASES: None
571 : *
572 : * RETURNS: int (could use errSprintf())
573 : * 0 = OK
574 : * -1 = tdlpLen is too small.
575 : *
576 : * HISTORY
577 : * 10/2004 Arthur Taylor (MDL): Created.
578 : *
579 : * NOTES
580 : *****************************************************************************
581 : */
582 0 : int TDLP_RefTime (DataSource &fp, sInt4 tdlpLen, double *refTime)
583 : {
584 : int sectLen; /* Length of section. */
585 : sInt4 curLoc; /* Where we are in the current TDLP message. */
586 : int c_temp; /* Temporary variable for use with fgetc */
587 : short int si_temp; /* Temporary variable. */
588 : int year, t_year; /* The reference year, and a consistency test */
589 : uChar month, t_month; /* The reference month, and a consistency test */
590 : uChar day, t_day; /* The reference day, and a consistency test */
591 : uChar hour, t_hour; /* The reference hour, and a consistency test */
592 : uChar min; /* The reference minute */
593 : sInt4 li_temp; /* Temporary variable. */
594 :
595 0 : if ((sectLen = fp.DataSourceFgetc ()) == EOF)
596 0 : goto error;
597 0 : curLoc = 8 + sectLen;
598 0 : if (curLoc > tdlpLen) {
599 0 : errSprintf ("Ran out of data in PDS (TDLP_RefTime)\n");
600 0 : return -1;
601 : }
602 0 : myAssert (sectLen <= 71);
603 0 : if (sectLen < 39) {
604 0 : errSprintf ("TDLP Section 1 is too small.\n");
605 0 : return -1;
606 : }
607 :
608 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
609 0 : goto error;
610 0 : if (FREAD_BIG (&si_temp, sizeof (short int), 1, fp) != 1)
611 0 : goto error;
612 0 : year = si_temp;
613 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
614 0 : goto error;
615 0 : month = c_temp;
616 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
617 0 : goto error;
618 0 : day = c_temp;
619 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
620 0 : goto error;
621 0 : hour = c_temp;
622 0 : if ((c_temp = fp.DataSourceFgetc()) == EOF)
623 0 : goto error;
624 0 : min = c_temp;
625 :
626 0 : if (FREAD_BIG (&li_temp, sizeof (sInt4), 1, fp) != 1)
627 0 : goto error;
628 0 : t_year = li_temp / 1000000L;
629 0 : li_temp -= t_year * 1000000L;
630 0 : t_month = li_temp / 10000L;
631 0 : li_temp -= t_month * 10000L;
632 0 : t_day = li_temp / 100;
633 0 : t_hour = li_temp - t_day * 100;
634 :
635 0 : if ((t_year != year) || (t_month != month) || (t_day != day) ||
636 : (t_hour != hour)) {
637 0 : errSprintf ("Error Inconsistant Times in TDLP_RefTime.\n");
638 0 : return -1;
639 : }
640 0 : if (ParseTime (refTime, year, month, day, hour, min, 0) != 0) {
641 0 : preErrSprintf ("Error In call to ParseTime in TDLP_RefTime.\n");
642 0 : return -1;
643 : }
644 :
645 : /* Get to the end of the TDLP message. */
646 : /* (inventory.c : GRIB2Inventory), is responsible for this. */
647 : /* fseek (fp, gribLen - sectLen, SEEK_CUR); */
648 0 : return 0;
649 : error:
650 0 : errSprintf ("Ran out of file in PDS (TDLP_RefTime).\n");
651 0 : return -1;
652 : }
653 :
654 : /*****************************************************************************
655 : * ReadTDLPSect2() --
656 : *
657 : * Arthur Taylor / MDL
658 : *
659 : * PURPOSE
660 : * Parses the TDLP "Grid Definition Section" or section 2, filling out
661 : * the gdsMeta data structure.
662 : *
663 : * ARGUMENTS
664 : * gds = The compressed part of the message dealing with "GDS". (Input)
665 : * tdlpLen = The total length of the TDLP message. (Input)
666 : * curLoc = Current location in the TDLP message. (Output)
667 : * gdsMeta = The filled out gdsMeta data structure. (Output)
668 : *
669 : * FILES/DATABASES: None
670 : *
671 : * RETURNS: int (could use errSprintf())
672 : * 0 = OK
673 : * -1 = tdlpLen is too small.
674 : * -2 = unexpected values in gds.
675 : *
676 : * HISTORY
677 : * 10/2004 Arthur Taylor (MDL): Created.
678 : *
679 : * NOTES
680 : *****************************************************************************
681 : */
682 0 : static int ReadTDLPSect2 (uChar *gds, sInt4 tdlpLen, sInt4 *curLoc,
683 : gdsType *gdsMeta)
684 : {
685 : char sectLen; /* Length of section. */
686 : int gridType; /* Which type of grid. (see enumerated types). */
687 : sInt4 li_temp; /* Temporary variable. */
688 :
689 0 : sectLen = *(gds++);
690 0 : *curLoc += sectLen;
691 0 : if (*curLoc > tdlpLen) {
692 0 : errSprintf ("Ran out of data in GDS (TDLP Section 2)\n");
693 0 : return -1;
694 : }
695 0 : myAssert (sectLen == 28);
696 :
697 0 : gridType = *(gds++);
698 0 : gdsMeta->Nx = GRIB_UNSIGN_INT2 (*gds, gds[1]);
699 0 : gds += 2;
700 0 : gdsMeta->Ny = GRIB_UNSIGN_INT2 (*gds, gds[1]);
701 0 : gds += 2;
702 0 : gdsMeta->lat1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
703 0 : gds += 3;
704 0 : gdsMeta->lon1 = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
705 0 : gdsMeta->lon1 = 360 - gdsMeta->lon1;
706 0 : if (gdsMeta->lon1 < 0) {
707 0 : gdsMeta->lon1 += 360;
708 : }
709 0 : if (gdsMeta->lon1 > 360) {
710 0 : gdsMeta->lon1 -= 360;
711 : }
712 0 : gds += 3;
713 0 : gdsMeta->orientLon = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
714 0 : gdsMeta->orientLon = 360 - gdsMeta->orientLon;
715 0 : if (gdsMeta->orientLon < 0) {
716 0 : gdsMeta->orientLon += 360;
717 : }
718 0 : if (gdsMeta->orientLon > 360) {
719 0 : gdsMeta->orientLon -= 360;
720 : }
721 0 : gds += 3;
722 0 : MEMCPY_BIG (&li_temp, gds, sizeof (sInt4));
723 0 : gdsMeta->Dx = li_temp / 1000.0;
724 0 : gds += 4;
725 0 : gdsMeta->meshLat = GRIB_SIGN_INT3 (*gds, gds[1], gds[2]) / 10000.0;
726 0 : gds += 3;
727 0 : if ((*gds != 0) || (gds[1] != 0) || (gds[2] != 0) || (gds[3] != 0) ||
728 0 : (gds[4] != 0) || (gds[5] != 0)) {
729 0 : errSprintf ("Error Reserved was not set to 0 in ReadTDLPSect2.\n");
730 0 : return -1;
731 : }
732 :
733 0 : gdsMeta->numPts = gdsMeta->Nx * gdsMeta->Ny;
734 0 : gdsMeta->f_sphere = 1;
735 0 : gdsMeta->majEarth = 6371.2;
736 0 : gdsMeta->minEarth = 6371.2;
737 0 : gdsMeta->Dy = gdsMeta->Dx;
738 0 : gdsMeta->resFlag = 0;
739 0 : gdsMeta->center = 0;
740 0 : gdsMeta->scan = 64;
741 0 : gdsMeta->lat2 = 0;
742 0 : gdsMeta->lon2 = 0;
743 0 : gdsMeta->scaleLat1 = gdsMeta->meshLat;
744 0 : gdsMeta->scaleLat2 = gdsMeta->meshLat;
745 0 : gdsMeta->southLat = 0;
746 0 : gdsMeta->southLon = 0;
747 0 : switch (gridType) {
748 : case TDLP_POLAR:
749 0 : gdsMeta->projType = GS3_POLAR;
750 : /* 4/24/2006 Added the following. */
751 0 : gdsMeta->scaleLat1 = 90;
752 0 : gdsMeta->scaleLat2 = 90;
753 0 : break;
754 : case TDLP_LAMBERT:
755 0 : gdsMeta->projType = GS3_LAMBERT;
756 0 : break;
757 : case TDLP_MERCATOR:
758 0 : gdsMeta->projType = GS3_MERCATOR;
759 : /* 4/24/2006 Added the following. */
760 0 : gdsMeta->scaleLat1 = 0;
761 0 : gdsMeta->scaleLat2 = 0;
762 0 : break;
763 : default:
764 0 : errSprintf ("Grid projection number is %d\n", gridType);
765 0 : return -2;
766 : }
767 0 : return 0;
768 : }
769 :
770 : /*****************************************************************************
771 : * ReadTDLPSect3() --
772 : *
773 : * Arthur Taylor / MDL
774 : *
775 : * PURPOSE
776 : * Parses the TDLP "Bit Map Section" or section 3, filling out the bitmap
777 : * as needed.
778 : *
779 : * ARGUMENTS
780 : * bms = The compressed part of the message dealing with "BMS". (Input)
781 : * tdlpLen = The total length of the TDLP message. (Input)
782 : * curLoc = Current location in the TDLP message. (Output)
783 : * bitmap = The extracted bitmap. (Output)
784 : * NxNy = The total size of the grid. (Input)
785 : *
786 : * FILES/DATABASES: None
787 : *
788 : * RETURNS: int (could use errSprintf())
789 : * 0 = OK
790 : * -1 = tdlpLen is too small.
791 : * -2 = unexpected values in bms.
792 : *
793 : * HISTORY
794 : * 10/2004 Arthur Taylor (MDL): Created.
795 : *
796 : * NOTES
797 : *****************************************************************************
798 : */
799 0 : static int ReadTDLPSect3 (uChar *bms, sInt4 tdlpLen, sInt4 *curLoc,
800 : uChar *bitmap, sInt4 NxNy)
801 : {
802 0 : errSprintf ("Bitmap data is Not Supported\n");
803 0 : return -1;
804 : }
805 :
806 : /*****************************************************************************
807 : * ReadTDLPSect4() --
808 : *
809 : * Arthur Taylor / MDL
810 : *
811 : * PURPOSE
812 : * Unpacks the "Binary Data Section" or section 4.
813 : *
814 : * ARGUMENTS
815 : * bds = The compressed part of the message dealing with "BDS". (Input)
816 : * tdlpLen = The total length of the TDLP message. (Input)
817 : * curLoc = Current location in the TDLP message. (Output)
818 : * DSF = Decimal Scale Factor for unpacking the data. (Input)
819 : * BSF = Binary Scale Factor for unpacking the data. (Input)
820 : * data = The extracted grid. (Output)
821 : * meta = The meta data associated with the grid (Input/Output)
822 : * unitM = The M unit conversion value in equation y = Mx + B. (Input)
823 : * unitB = The B unit conversion value in equation y = Mx + B. (Input)
824 : *
825 : * FILES/DATABASES: None
826 : *
827 : * RETURNS: int (could use errSprintf())
828 : * 0 = OK
829 : * -1 = gribLen is too small.
830 : * -2 = unexpected values in bds.
831 : *
832 : * HISTORY
833 : * 4/2003 Arthur Taylor (MDL/RSIS): Created
834 : * 3/2004 AAT: Switched {# Pts * (# Bits in a Group) +
835 : * # of unused bits != # of available bits} to a warning from an
836 : * error.
837 : * 2/2005 AAT: Found bug: memBitRead grp[i].bit was sizeof sInt4 instead
838 : * of uChar.
839 : * 2/2005 AAT: Second order diff, no miss value bug (lastData - 1) should
840 : * be lastData.
841 : * 2/2005 AAT: Added test to see if the number of bits needed matches the
842 : * section length.
843 : *
844 : * NOTES
845 : * 1) See metaparse.c : ParseGrid()
846 : *****************************************************************************
847 : */
848 0 : static int ReadTDLPSect4 (uChar *bds, sInt4 tdlpLen, sInt4 *curLoc,
849 : short int DSF, short int BSF, double *data,
850 : grib_MetaData *meta, double unitM, double unitB)
851 : {
852 : uInt4 sectLen; /* Length in bytes of the current section. */
853 : uChar f_notGridPnt; /* Not Grid point data? */
854 : uChar f_complexPack; /* Complex packing? */
855 : uChar f_sndOrder; /* Second order differencing? */
856 : uChar f_primMiss; /* Primary missing value? */
857 : uChar f_secMiss; /* Secondary missing value? */
858 : uInt4 numPack; /* Number of points packed. */
859 : sInt4 li_temp; /* Temporary variable. */
860 : uInt4 uli_temp; /* Temporary variable. */
861 : uChar bufLoc; /* Keeps track of where to start getting more data
862 : * out of the packed data stream. */
863 : uChar f_negative; /* used to help with signs of numbers. */
864 : size_t numUsed; /* How many bytes were used in a given call to
865 : * memBitRead. */
866 0 : sInt4 origVal = 0; /* Original value. */
867 : uChar mbit; /* # of bits for abs (first first order difference) */
868 0 : sInt4 fstDiff = 0; /* First first order difference. */
869 0 : sInt4 diff = 0; /* general first order difference. */
870 : uChar nbit; /* # of bits for abs (overall min value) */
871 : sInt4 minVal; /* Minimum value. */
872 : size_t LX; /* Number of groups. */
873 : uChar ibit; /* # of bits for group min values. */
874 : uChar jbit; /* # of bits for # of bits for group. */
875 : uChar kbit; /* # of bits for # values in a group. */
876 : TDLGroupType *grp; /* Holds the info about each group. */
877 : size_t i, j; /* Loop counters. */
878 : uInt4 t_numPack; /* Used to total number of values in a group to check
879 : * the numPack value. */
880 : uInt4 t_numBits; /* Used to total number of bits used in the groups. */
881 : uInt4 t_numBytes; /* Used to total number of bytes used to compare to
882 : * sectLen. */
883 : sInt4 maxVal; /* The max value in a group. */
884 : uInt4 dataCnt; /* How many values (miss or othewise) we have read. */
885 : uInt4 lastData; /* Index to last actual data. */
886 : uInt4 numVal; /* # of actual (non-missing values) we have. */
887 : double scale; /* Amount to scale values by. */
888 : uInt4 dataInd; /* Index into data for this value (used to switch
889 : * from a11..a1n,a2n..a21 to normal grid of
890 : * a11..a1n,a21..a2n. */
891 : uChar f_missing; /* Used to help with primary missing values, and the
892 : * 0 bit possibility. */
893 : #ifdef DEBUG
894 0 : sInt4 t_UK1 = 0; /* Used to test theories about un defined values. */
895 0 : sInt4 t_UK2 = 0; /* Used to test theories about un defined values. */
896 : #endif
897 :
898 0 : sectLen = GRIB_UNSIGN_INT3 (*bds, bds[1], bds[2]);
899 0 : *curLoc += sectLen;
900 0 : if (*curLoc > tdlpLen) {
901 0 : errSprintf ("Ran out of data in BDS (TDLP Section 4)\n");
902 0 : return -1;
903 : }
904 0 : bds += 3;
905 0 : t_numBytes = 3;
906 0 : f_notGridPnt = (GRIB2BIT_4 & *bds) ? 1 : 0;
907 0 : f_complexPack = (GRIB2BIT_5 & *bds) ? 1 : 0;
908 0 : f_sndOrder = (GRIB2BIT_6 & *bds) ? 1 : 0;
909 0 : f_primMiss = (GRIB2BIT_7 & *bds) ? 1 : 0;
910 0 : f_secMiss = (GRIB2BIT_8 & *bds) ? 1 : 0;
911 :
912 0 : if (f_secMiss && (!f_primMiss)) {
913 0 : errSprintf ("Secondary missing value without a primary!\n");
914 0 : return -1;
915 : }
916 0 : if (f_complexPack) {
917 0 : if (!f_sndOrder) {
918 0 : meta->gridAttrib.packType = GS5_CMPLX;
919 : } else {
920 0 : meta->gridAttrib.packType = GS5_CMPLXSEC;
921 : }
922 : } else {
923 0 : errSprintf ("Simple pack is not supported at this time.\n");
924 0 : return -1;
925 : }
926 0 : bds++;
927 0 : MEMCPY_BIG (&numPack, bds, sizeof (sInt4));
928 0 : bds += 4;
929 0 : t_numBytes += 5;
930 0 : if (!f_notGridPnt) {
931 0 : if (numPack != meta->gds.numPts) {
932 : errSprintf ("Number packed %d != number of points %d\n", numPack,
933 0 : meta->gds.numPts);
934 0 : return -1;
935 : }
936 : }
937 0 : meta->gridAttrib.DSF = DSF;
938 0 : meta->gridAttrib.ESF = BSF;
939 0 : meta->gridAttrib.fieldType = 0;
940 0 : if (!f_primMiss) {
941 0 : meta->gridAttrib.f_miss = 0;
942 0 : meta->gridAttrib.missPri = 0;
943 0 : meta->gridAttrib.missSec = 0;
944 : } else {
945 0 : MEMCPY_BIG (&li_temp, bds, sizeof (sInt4));
946 0 : bds += 4;
947 0 : t_numBytes += 4;
948 0 : meta->gridAttrib.missPri = li_temp / 10000.0;
949 0 : if (!f_secMiss) {
950 0 : meta->gridAttrib.f_miss = 1;
951 0 : meta->gridAttrib.missSec = 0;
952 : } else {
953 0 : MEMCPY_BIG (&li_temp, bds, sizeof (sInt4));
954 0 : bds += 4;
955 0 : t_numBytes += 4;
956 0 : meta->gridAttrib.missSec = li_temp / 10000.0;
957 0 : meta->gridAttrib.f_miss = 2;
958 : }
959 : }
960 : /* Init the buffer location. */
961 0 : bufLoc = 8;
962 : /* The origValue and fstDiff are only present if sndOrder packed. */
963 0 : if (f_sndOrder) {
964 : memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc,
965 0 : &numUsed);
966 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, 31, &bufLoc, &numUsed);
967 0 : myAssert (numUsed == 4);
968 0 : bds += numUsed;
969 0 : t_numBytes += numUsed;
970 0 : origVal = (f_negative) ? -1 * uli_temp : uli_temp;
971 0 : memBitRead (&mbit, sizeof (mbit), bds, 5, &bufLoc, &numUsed);
972 : memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc,
973 0 : &numUsed);
974 0 : myAssert (numUsed == 0);
975 0 : myAssert ((mbit > 0) && (mbit < 32));
976 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, mbit, &bufLoc, &numUsed);
977 0 : bds += numUsed;
978 0 : t_numBytes += numUsed;
979 0 : fstDiff = (f_negative) ? -1 * uli_temp : uli_temp;
980 : }
981 0 : memBitRead (&nbit, sizeof (nbit), bds, 5, &bufLoc, &numUsed);
982 0 : bds += numUsed;
983 0 : t_numBytes += numUsed;
984 0 : memBitRead (&f_negative, sizeof (f_negative), bds, 1, &bufLoc, &numUsed);
985 0 : bds += numUsed;
986 0 : t_numBytes += numUsed;
987 0 : myAssert ((nbit > 0) && (nbit < 32));
988 0 : memBitRead (&uli_temp, sizeof (sInt4), bds, nbit, &bufLoc, &numUsed);
989 0 : bds += numUsed;
990 0 : t_numBytes += numUsed;
991 0 : minVal = (f_negative) ? -1 * uli_temp : uli_temp;
992 0 : memBitRead (&LX, sizeof (LX), bds, 16, &bufLoc, &numUsed);
993 0 : bds += numUsed;
994 0 : t_numBytes += numUsed;
995 0 : grp = (TDLGroupType *) malloc (LX * sizeof (TDLGroupType));
996 0 : memBitRead (&ibit, sizeof (ibit), bds, 5, &bufLoc, &numUsed);
997 0 : bds += numUsed;
998 0 : t_numBytes += numUsed;
999 0 : memBitRead (&jbit, sizeof (jbit), bds, 5, &bufLoc, &numUsed);
1000 : /* Following assert is because it is the # of bits of # of bits. Which
1001 : * means that # of bits of value that has a max of 64. */
1002 0 : myAssert (jbit < 6);
1003 0 : bds += numUsed;
1004 0 : t_numBytes += numUsed;
1005 0 : memBitRead (&kbit, sizeof (kbit), bds, 5, &bufLoc, &numUsed);
1006 0 : bds += numUsed;
1007 0 : t_numBytes += numUsed;
1008 0 : myAssert (ibit < 33);
1009 0 : for (i = 0; i < LX; i++) {
1010 0 : if (ibit == 0) {
1011 0 : grp[i].min = 0;
1012 : } else {
1013 0 : memBitRead (&(grp[i].min), sizeof (sInt4), bds, ibit, &bufLoc,
1014 0 : &numUsed);
1015 0 : bds += numUsed;
1016 0 : t_numBytes += numUsed;
1017 : }
1018 : }
1019 0 : myAssert (jbit < 8);
1020 0 : for (i = 0; i < LX; i++) {
1021 0 : if (jbit == 0) {
1022 0 : grp[i].bit = 0;
1023 : } else {
1024 0 : myAssert (jbit <= sizeof (uChar) * 8);
1025 0 : memBitRead (&(grp[i].bit), sizeof (uChar), bds, jbit, &bufLoc,
1026 0 : &numUsed);
1027 0 : bds += numUsed;
1028 0 : t_numBytes += numUsed;
1029 : }
1030 0 : myAssert (grp[i].bit < 32);
1031 : }
1032 0 : myAssert (kbit < 33);
1033 0 : t_numPack = 0;
1034 0 : t_numBits = 0;
1035 0 : for (i = 0; i < LX; i++) {
1036 0 : if (kbit == 0) {
1037 0 : grp[i].num = 0;
1038 : } else {
1039 0 : memBitRead (&(grp[i].num), sizeof (sInt4), bds, kbit, &bufLoc,
1040 0 : &numUsed);
1041 0 : bds += numUsed;
1042 0 : t_numBytes += numUsed;
1043 : }
1044 0 : t_numPack += grp[i].num;
1045 0 : t_numBits += grp[i].num * grp[i].bit;
1046 : }
1047 0 : if (t_numPack != numPack) {
1048 : errSprintf ("Number packed %d != number of values in groups %d\n",
1049 0 : numPack, t_numPack);
1050 0 : free (grp);
1051 0 : return -1;
1052 : }
1053 0 : if ((t_numBytes + ceil (t_numBits / 8.)) > sectLen) {
1054 : errSprintf ("# bytes in groups %ld (%ld + %ld / 8) > sectLen %ld\n",
1055 : (sInt4) (t_numBytes + ceil (t_numBits / 8.)),
1056 0 : t_numBytes, t_numBits, sectLen);
1057 0 : free (grp);
1058 0 : return -1;
1059 : }
1060 0 : dataCnt = 0;
1061 0 : dataInd = 0;
1062 :
1063 : #ifdef DEBUG
1064 0 : printf ("nbit %d, ibit %d, jbit %d, kbit %d\n", nbit, ibit, jbit, kbit);
1065 0 : if ((t_numBytes + ceil (t_numBits / 8.)) != sectLen) {
1066 : printf ("Caution: # bytes in groups %d (%d + %d / 8) != "
1067 : "sectLen %d\n", (sInt4) (t_numBytes + ceil (t_numBits / 8.)),
1068 0 : t_numBytes, t_numBits, sectLen);
1069 : }
1070 : #endif
1071 :
1072 : /* Binary scale factor in TDLP has reverse sign from GRIB definition. */
1073 0 : scale = pow (10.0, -1 * DSF) * pow (2.0, -1 * BSF);
1074 :
1075 0 : meta->gridAttrib.f_maxmin = 0;
1076 : /* Work with Second order complex packed data. */
1077 0 : if (f_sndOrder) {
1078 : /* *INDENT-OFF* */
1079 : /* The algorithm appears to be:
1080 : * Data: a1 a2 a3 a4 a5 ...
1081 : * 1st diff: 0 b2 b3 b4 b5 ...
1082 : * 2nd diff: UK1 UK2 c3 c4 c5 ...
1083 : * We already know a1 and b2, and unpack a stream of UK1 UK2 c3 c4
1084 : * The problem is that UK1, UK2 is undefined. Originally I thought
1085 : * this was 0, or c3, but it appears that if b2 != 0, then
1086 : * UK2 = c3 + 2 b2, and UK1 = c3 + 1 * b2, otherwise it appears that
1087 : * UK1 == UK2, and typically UK1 == c3 (but not always). */
1088 : /* *INDENT-ON* */
1089 0 : myAssert (numPack >= 2);
1090 0 : if (f_secMiss) {
1091 0 : numVal = 0;
1092 0 : lastData = 0;
1093 0 : for (i = 0; i < LX; i++) {
1094 0 : maxVal = (1 << grp[i].bit) - 1;
1095 0 : for (j = 0; j < grp[i].num; j++) {
1096 : /* signed int. */
1097 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1098 0 : &bufLoc, &numUsed);
1099 0 : bds += numUsed;
1100 0 : if (li_temp == maxVal) {
1101 0 : data[dataInd] = meta->gridAttrib.missPri;
1102 0 : } else if (li_temp == (maxVal - 1)) {
1103 0 : data[dataInd] = meta->gridAttrib.missSec;
1104 : } else {
1105 0 : if (numVal > 1) {
1106 : #ifdef DEBUG
1107 0 : if (numVal == 2) {
1108 0 : if (fstDiff != 0) {
1109 : /*
1110 : myAssert (t_UK1 == li_temp + fstDiff);
1111 : myAssert (t_UK2 == li_temp + 2 * fstDiff);
1112 : */
1113 : } else {
1114 0 : myAssert (t_UK1 == t_UK2);
1115 : }
1116 : }
1117 : #endif
1118 0 : diff += (li_temp + grp[i].min + minVal);
1119 0 : data[dataInd] = data[lastData] + diff * scale;
1120 0 : lastData = dataInd;
1121 0 : } else if (numVal == 1) {
1122 0 : data[dataInd] = (origVal + fstDiff) * scale;
1123 0 : lastData = dataInd;
1124 0 : diff = fstDiff;
1125 : #ifdef DEBUG
1126 0 : t_UK2 = li_temp;
1127 : #endif
1128 : } else {
1129 0 : data[dataInd] = origVal * scale;
1130 : #ifdef DEBUG
1131 0 : t_UK1 = li_temp;
1132 : #endif
1133 : }
1134 0 : numVal++;
1135 0 : if (!meta->gridAttrib.f_maxmin) {
1136 0 : meta->gridAttrib.min = data[dataInd];
1137 0 : meta->gridAttrib.max = data[dataInd];
1138 0 : meta->gridAttrib.f_maxmin = 1;
1139 : } else {
1140 0 : if (data[dataInd] < meta->gridAttrib.min) {
1141 0 : meta->gridAttrib.min = data[dataInd];
1142 : }
1143 0 : if (data[dataInd] > meta->gridAttrib.max) {
1144 0 : meta->gridAttrib.max = data[dataInd];
1145 : }
1146 : }
1147 : }
1148 0 : dataCnt++;
1149 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1150 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1151 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1152 0 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1153 : }
1154 : }
1155 0 : } else if (f_primMiss) {
1156 0 : numVal = 0;
1157 0 : lastData = 0;
1158 0 : for (i = 0; i < LX; i++) {
1159 0 : maxVal = (1 << grp[i].bit) - 1;
1160 0 : for (j = 0; j < grp[i].num; j++) {
1161 : /* signed int. */
1162 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1163 0 : &bufLoc, &numUsed);
1164 0 : bds += numUsed;
1165 0 : f_missing = 0;
1166 0 : if (li_temp == maxVal) {
1167 0 : data[dataInd] = meta->gridAttrib.missPri;
1168 : /* In the case of grp[i].bit == 0, if grp[i].min == 0, then
1169 : * it is the missing value, otherwise regular value. Only
1170 : * need to be concerned for primary missing values. */
1171 0 : f_missing = 1;
1172 0 : if ((grp[i].bit == 0) && (grp[i].min != 0)) {
1173 : #ifdef DEBUG
1174 0 : printf ("This doesn't happen often.\n");
1175 0 : printf ("%d %d %d\n", (int) i, grp[i].bit, grp[i].min);
1176 : #endif
1177 0 : myAssert (1 == 2);
1178 0 : f_missing = 0;
1179 : }
1180 : }
1181 0 : if (!f_missing) {
1182 0 : if (numVal > 1) {
1183 : #ifdef DEBUG
1184 0 : if (numVal == 2) {
1185 0 : if (fstDiff != 0) {
1186 : /*
1187 : myAssert (t_UK1 == li_temp + fstDiff);
1188 : myAssert (t_UK2 == li_temp + 2 * fstDiff);
1189 : */
1190 : } else {
1191 0 : myAssert (t_UK1 == t_UK2);
1192 : }
1193 : }
1194 : #endif
1195 0 : diff += (li_temp + grp[i].min + minVal);
1196 0 : data[dataInd] = data[lastData] + diff * scale;
1197 0 : lastData = dataInd;
1198 0 : } else if (numVal == 1) {
1199 0 : data[dataInd] = (origVal + fstDiff) * scale;
1200 0 : lastData = dataInd;
1201 0 : diff = fstDiff;
1202 : #ifdef DEBUG
1203 0 : t_UK2 = li_temp;
1204 : #endif
1205 : } else {
1206 0 : data[dataInd] = origVal * scale;
1207 : #ifdef DEBUG
1208 0 : t_UK1 = li_temp;
1209 : #endif
1210 : }
1211 0 : numVal++;
1212 0 : if (!meta->gridAttrib.f_maxmin) {
1213 0 : meta->gridAttrib.min = data[dataInd];
1214 0 : meta->gridAttrib.max = data[dataInd];
1215 0 : meta->gridAttrib.f_maxmin = 1;
1216 : } else {
1217 0 : if (data[dataInd] < meta->gridAttrib.min) {
1218 0 : meta->gridAttrib.min = data[dataInd];
1219 : }
1220 0 : if (data[dataInd] > meta->gridAttrib.max) {
1221 0 : meta->gridAttrib.max = data[dataInd];
1222 : }
1223 : }
1224 : }
1225 0 : dataCnt++;
1226 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1227 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1228 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1229 0 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1230 : }
1231 : }
1232 : } else {
1233 0 : lastData = 0;
1234 0 : for (i = 0; i < LX; i++) {
1235 0 : for (j = 0; j < grp[i].num; j++) {
1236 : /* signed int. */
1237 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1238 0 : &bufLoc, &numUsed);
1239 0 : bds += numUsed;
1240 0 : if (dataCnt > 1) {
1241 : #ifdef DEBUG
1242 0 : if (dataCnt == 2) {
1243 0 : if (fstDiff != 0) {
1244 : /*
1245 : myAssert (t_UK1 == li_temp + fstDiff);
1246 : myAssert (t_UK2 == li_temp + 2 * fstDiff);
1247 : */
1248 : } else {
1249 0 : myAssert (t_UK1 == t_UK2);
1250 : }
1251 : }
1252 : #endif
1253 0 : diff += (li_temp + grp[i].min + minVal);
1254 0 : data[dataInd] = data[lastData] + diff * scale;
1255 0 : lastData = dataInd;
1256 0 : } else if (dataCnt == 1) {
1257 0 : data[dataInd] = (origVal + fstDiff) * scale;
1258 0 : lastData = dataInd;
1259 0 : diff = fstDiff;
1260 : #ifdef DEBUG
1261 0 : t_UK2 = li_temp;
1262 : #endif
1263 : } else {
1264 0 : data[dataInd] = origVal * scale;
1265 : #ifdef DEBUG
1266 0 : t_UK1 = li_temp;
1267 : #endif
1268 : }
1269 : #ifdef DEBUG
1270 : /*
1271 : if (i >= 4153) {
1272 : */
1273 : /*
1274 : if ((data[dataInd] > 100) || (data[dataInd] < -100)) {
1275 : */
1276 : /*
1277 : if ((diff > 50) || (diff < -50)) {
1278 : printf ("li_temp :: %ld, diff = %ld\n", li_temp, diff);
1279 : printf ("data[dataInd] :: %f\n", data[dataInd]);
1280 : printf ("Group # %d element %d, grp[i].min %ld, "
1281 : "grp[i].bit %d, minVal %ld\n", i, j, grp[i].min,
1282 : grp[i].bit, minVal);
1283 : }
1284 : */
1285 : #endif
1286 0 : if (!meta->gridAttrib.f_maxmin) {
1287 0 : meta->gridAttrib.min = data[dataInd];
1288 0 : meta->gridAttrib.max = data[dataInd];
1289 0 : meta->gridAttrib.f_maxmin = 1;
1290 : } else {
1291 0 : if (data[dataInd] < meta->gridAttrib.min) {
1292 0 : meta->gridAttrib.min = data[dataInd];
1293 : }
1294 0 : if (data[dataInd] > meta->gridAttrib.max) {
1295 0 : meta->gridAttrib.max = data[dataInd];
1296 : }
1297 : }
1298 0 : dataCnt++;
1299 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1300 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1301 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1302 0 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1303 : }
1304 : }
1305 0 : numVal = dataCnt;
1306 : }
1307 :
1308 : /* Work with regular complex packed data. */
1309 : } else {
1310 : #ifdef DEBUG
1311 : /*
1312 : printf ("Work with regular complex packed data\n");
1313 : */
1314 : #endif
1315 0 : if (f_secMiss) {
1316 0 : numVal = 0;
1317 0 : for (i = 0; i < LX; i++) {
1318 0 : maxVal = (1 << grp[i].bit) - 1;
1319 0 : for (j = 0; j < grp[i].num; j++) {
1320 : /* signed int. */
1321 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1322 0 : &bufLoc, &numUsed);
1323 0 : bds += numUsed;
1324 0 : if (li_temp == maxVal) {
1325 0 : data[dataInd] = meta->gridAttrib.missPri;
1326 0 : } else if (li_temp == (maxVal - 1)) {
1327 0 : data[dataInd] = meta->gridAttrib.missSec;
1328 : } else {
1329 0 : data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1330 0 : numVal++;
1331 0 : if (!meta->gridAttrib.f_maxmin) {
1332 0 : meta->gridAttrib.min = data[dataInd];
1333 0 : meta->gridAttrib.max = data[dataInd];
1334 0 : meta->gridAttrib.f_maxmin = 1;
1335 : } else {
1336 0 : if (data[dataInd] < meta->gridAttrib.min) {
1337 0 : meta->gridAttrib.min = data[dataInd];
1338 : }
1339 0 : if (data[dataInd] > meta->gridAttrib.max) {
1340 0 : meta->gridAttrib.max = data[dataInd];
1341 : }
1342 : }
1343 : }
1344 0 : dataCnt++;
1345 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1346 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1347 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1348 0 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1349 : }
1350 : }
1351 0 : } else if (f_primMiss) {
1352 : #ifdef DEBUG
1353 : /*
1354 : printf ("Work with primary missing data\n");
1355 : */
1356 : #endif
1357 0 : numVal = 0;
1358 0 : for (i = 0; i < LX; i++) {
1359 0 : maxVal = (1 << grp[i].bit) - 1;
1360 0 : for (j = 0; j < grp[i].num; j++) {
1361 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1362 0 : &bufLoc, &numUsed);
1363 0 : bds += numUsed;
1364 0 : f_missing = 0;
1365 0 : if (li_temp == maxVal) {
1366 0 : data[dataInd] = meta->gridAttrib.missPri;
1367 : /* In the case of grp[i].bit == 0, if grp[i].min == 0, then
1368 : * it is the missing value, otherwise regular value. Only
1369 : * need to be concerned for primary missing values. */
1370 0 : f_missing = 1;
1371 0 : if ((grp[i].bit == 0) && (grp[i].min != 0)) {
1372 : #ifdef DEBUG
1373 0 : printf ("This doesn't happen often.\n");
1374 0 : printf ("%d %d %d\n", (int) i, grp[i].bit, grp[i].min);
1375 0 : myAssert (1 == 2);
1376 : #endif
1377 0 : f_missing = 0;
1378 : }
1379 : }
1380 0 : if (!f_missing) {
1381 0 : data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1382 0 : numVal++;
1383 0 : if (!meta->gridAttrib.f_maxmin) {
1384 0 : meta->gridAttrib.min = data[dataInd];
1385 0 : meta->gridAttrib.max = data[dataInd];
1386 0 : meta->gridAttrib.f_maxmin = 1;
1387 : } else {
1388 0 : if (data[dataInd] < meta->gridAttrib.min) {
1389 0 : meta->gridAttrib.min = data[dataInd];
1390 : }
1391 0 : if (data[dataInd] > meta->gridAttrib.max) {
1392 0 : meta->gridAttrib.max = data[dataInd];
1393 : }
1394 : }
1395 : }
1396 0 : dataCnt++;
1397 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1398 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1399 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1400 0 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1401 : }
1402 : }
1403 : } else {
1404 0 : for (i = 0; i < LX; i++) {
1405 0 : for (j = 0; j < grp[i].num; j++) {
1406 0 : memBitRead (&(li_temp), sizeof (sInt4), bds, grp[i].bit,
1407 0 : &bufLoc, &numUsed);
1408 0 : bds += numUsed;
1409 0 : data[dataInd] = (li_temp + grp[i].min + minVal) * scale;
1410 0 : if (!meta->gridAttrib.f_maxmin) {
1411 0 : meta->gridAttrib.min = data[dataInd];
1412 0 : meta->gridAttrib.max = data[dataInd];
1413 0 : meta->gridAttrib.f_maxmin = 1;
1414 : } else {
1415 0 : if (data[dataInd] < meta->gridAttrib.min) {
1416 0 : meta->gridAttrib.min = data[dataInd];
1417 : }
1418 0 : if (data[dataInd] > meta->gridAttrib.max) {
1419 0 : meta->gridAttrib.max = data[dataInd];
1420 : }
1421 : }
1422 0 : dataCnt++;
1423 : dataInd = ((dataCnt / meta->gds.Nx) % 2) ?
1424 : (2 * (dataCnt / meta->gds.Nx) + 1) *
1425 0 : meta->gds.Nx - dataCnt - 1 : dataCnt;
1426 0 : myAssert ((dataInd < numPack) || (dataCnt == numPack));
1427 : }
1428 : }
1429 0 : numVal = dataCnt;
1430 : }
1431 : }
1432 0 : meta->gridAttrib.numMiss = dataCnt - numVal;
1433 0 : meta->gridAttrib.refVal = minVal * scale;
1434 :
1435 0 : free (grp);
1436 0 : return 0;
1437 : }
1438 :
1439 : /*****************************************************************************
1440 : * ReadTDLPRecord() --
1441 : *
1442 : * Arthur Taylor / MDL
1443 : *
1444 : * PURPOSE
1445 : * Reads in a TDLP message, and parses the data into various data
1446 : * structures, for use with other code.
1447 : *
1448 : * ARGUMENTS
1449 : * fp = An opened TDLP file already at the correct message. (Input)
1450 : * TDLP_Data = The read in TDLP data. (Output)
1451 : * tdlp_DataLen = Size of TDLP_Data. (Output)
1452 : * meta = A filled in meta structure (Output)
1453 : * IS = The structure containing all the arrays that the
1454 : * unpacker uses (Output)
1455 : * sect0 = Already read in section 0 data. (Input)
1456 : * tdlpLen = Length of the TDLP message. (Input)
1457 : * majEarth = Used to override the TDLP major axis of earth. (Input)
1458 : * minEarth = Used to override the TDLP minor axis of earth. (Input)
1459 : *
1460 : * FILES/DATABASES:
1461 : * An already opened file pointing to the desired TDLP message.
1462 : *
1463 : * RETURNS: int (could use errSprintf())
1464 : * 0 = OK
1465 : * -1 = Problems reading in the PDS.
1466 : * -2 = Problems reading in the GDS.
1467 : * -3 = Problems reading in the BMS.
1468 : * -4 = Problems reading in the BDS.
1469 : * -5 = Problems reading the closing section.
1470 : *
1471 : * HISTORY
1472 : * 10/2004 Arthur Taylor (MDL): Created
1473 : *
1474 : * NOTES
1475 : *****************************************************************************
1476 : */
1477 0 : int ReadTDLPRecord (DataSource &fp, double **TDLP_Data, uInt4 *tdlp_DataLen,
1478 : grib_MetaData *meta, IS_dataType *IS,
1479 : sInt4 sect0[SECT0LEN_WORD], uInt4 tdlpLen,
1480 : double majEarth, double minEarth)
1481 : {
1482 : sInt4 nd5; /* Size of TDLP message rounded up to the nearest *
1483 : * sInt4. */
1484 : uChar *c_ipack; /* A char ptr to the message stored in IS->ipack */
1485 : sInt4 curLoc; /* Current location in the GRIB message. */
1486 : char f_gds; /* flag if there is a GDS section. */
1487 : char f_bms; /* flag if there is a BMS section. */
1488 : short int DSF; /* Decimal Scale Factor for unpacking the data. */
1489 : short int BSF; /* Binary Scale Factor for unpacking the data. */
1490 : double *tdlp_Data; /* A pointer to TDLP_Data for ease of manipulation. */
1491 0 : double unitM = 1; /* M in y = Mx + B, for unit conversion. */
1492 0 : double unitB = 0; /* B in y = Mx + B, for unit conversion. */
1493 : sInt4 li_temp; /* Used to make sure section 5 is 7777. */
1494 : size_t pad; /* Number of bytes to pad the message to get to the
1495 : * correct byte boundary. */
1496 : char buffer[24]; /* Read the trailing bytes in the TDLPack record. */
1497 : uChar *bitmap; /* Would contain bitmap data if it was supported. */
1498 :
1499 : /* Make room for entire message, and read it in. */
1500 : /* nd5 needs to be tdlpLen in (sInt4) units rounded up. */
1501 0 : nd5 = (tdlpLen + 3) / 4;
1502 0 : if (nd5 > IS->ipackLen) {
1503 0 : IS->ipackLen = nd5;
1504 : IS->ipack = (sInt4 *) realloc ((void *) (IS->ipack),
1505 0 : (IS->ipackLen) * sizeof (sInt4));
1506 : }
1507 0 : c_ipack = (uChar *) IS->ipack;
1508 : /* Init last sInt4 to 0, to make sure that the padded bytes are 0. */
1509 0 : IS->ipack[nd5 - 1] = 0;
1510 : /* Init first 2 sInt4 to sect0. */
1511 0 : memcpy (c_ipack, sect0, SECT0LEN_WORD * 2);
1512 : /* Read in the rest of the message. */
1513 0 : if (fp.DataSourceFread(c_ipack + SECT0LEN_WORD * 2, sizeof (char),
1514 0 : (tdlpLen - SECT0LEN_WORD * 2)) + SECT0LEN_WORD * 2 != tdlpLen) {
1515 0 : errSprintf ("Ran out of file\n");
1516 0 : return -1;
1517 : }
1518 :
1519 : /* Preceeding was in degrib2, next part is specific to TDLP. */
1520 0 : curLoc = 8;
1521 0 : if (ReadTDLPSect1 (c_ipack + curLoc, tdlpLen, &curLoc, &(meta->pdsTdlp),
1522 : &f_gds, &f_bms, &DSF, &BSF) != 0) {
1523 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1524 0 : return -1;
1525 : }
1526 :
1527 : /* Figure out some basic stuff about the grid. */
1528 0 : free (meta->element);
1529 0 : meta->element = NULL;
1530 0 : free (meta->unitName);
1531 0 : meta->unitName = NULL;
1532 0 : free (meta->comment);
1533 0 : meta->comment = NULL;
1534 0 : free (meta->shortFstLevel);
1535 0 : meta->shortFstLevel = NULL;
1536 0 : free (meta->longFstLevel);
1537 0 : meta->longFstLevel = NULL;
1538 : TDLP_ElemSurfUnit (&(meta->pdsTdlp), &(meta->element), &(meta->unitName),
1539 : &(meta->comment), &(meta->shortFstLevel),
1540 0 : &(meta->longFstLevel));
1541 0 : meta->center = 7; /* US NWS, NCEP */
1542 0 : meta->subcenter = 14; /* NWS Meteorological Development Laboratory */
1543 :
1544 : /* strftime (meta->refTime, 20, "%Y%m%d%H%M",
1545 : gmtime (&(meta->pdsTdlp.refTime)));
1546 : */
1547 0 : Clock_Print (meta->refTime, 20, meta->pdsTdlp.refTime, "%Y%m%d%H%M", 0);
1548 :
1549 : /*
1550 : validTime = meta->pdsTdlp.refTime + meta->pdsTdlp.project;
1551 : strftime (meta->validTime, 20, "%Y%m%d%H%M", gmtime (&(validTime)));
1552 : */
1553 : Clock_Print (meta->validTime, 20, meta->pdsTdlp.refTime +
1554 0 : meta->pdsTdlp.project, "%Y%m%d%H%M", 0);
1555 :
1556 0 : meta->deltTime = meta->pdsTdlp.project;
1557 :
1558 : /* Get the Grid Definition Section. */
1559 0 : if (f_gds) {
1560 0 : if (ReadTDLPSect2 (c_ipack + curLoc, tdlpLen, &curLoc,
1561 : &(meta->gds)) != 0) {
1562 0 : preErrSprintf ("Inside ReadGrib1Record\n");
1563 0 : return -2;
1564 : }
1565 : } else {
1566 0 : errSprintf ("Don't know how to handle vector data yet.\n");
1567 0 : return -2;
1568 : }
1569 :
1570 : /* Allow over ride of the earth radii. */
1571 0 : if ((majEarth > 6300) && (majEarth < 6400)) {
1572 0 : if ((minEarth > 6300) && (minEarth < 6400)) {
1573 0 : meta->gds.f_sphere = 0;
1574 0 : meta->gds.majEarth = majEarth;
1575 0 : meta->gds.minEarth = minEarth;
1576 0 : if (majEarth == minEarth) {
1577 0 : meta->gds.f_sphere = 1;
1578 : }
1579 : } else {
1580 0 : meta->gds.f_sphere = 1;
1581 0 : meta->gds.majEarth = majEarth;
1582 0 : meta->gds.minEarth = majEarth;
1583 : }
1584 : }
1585 :
1586 : /* Allocate memory for the grid. */
1587 0 : if (meta->gds.numPts > *tdlp_DataLen) {
1588 0 : *tdlp_DataLen = meta->gds.numPts;
1589 : *TDLP_Data = (double *) realloc ((void *) (*TDLP_Data),
1590 0 : (*tdlp_DataLen) * sizeof (double));
1591 : }
1592 0 : tdlp_Data = *TDLP_Data;
1593 :
1594 : /* Get the Bit Map Section. */
1595 0 : if (f_bms) {
1596 : /* errSprintf ("Bitmap data is Not Supported\n");*/
1597 : /* Need to allocate bitmap when this is implemented. */
1598 0 : bitmap = NULL;
1599 : ReadTDLPSect3 (c_ipack + curLoc, tdlpLen, &curLoc, bitmap,
1600 0 : meta->gds.numPts);
1601 0 : return -1;
1602 : }
1603 :
1604 : /* Read the GRID. */
1605 0 : if (ReadTDLPSect4 (c_ipack + curLoc, tdlpLen, &curLoc, DSF, BSF,
1606 : tdlp_Data, meta, unitM, unitB) != 0) {
1607 0 : preErrSprintf ("Inside ReadTDLPRecord\n");
1608 0 : return -4;
1609 : }
1610 :
1611 : /* Read section 5. If it is "7777" == 926365495 we are done. */
1612 0 : memcpy (&li_temp, c_ipack + curLoc, 4);
1613 0 : if (li_temp != 926365495L) {
1614 0 : errSprintf ("Did not find the end of the message.\n");
1615 0 : return -5;
1616 : }
1617 0 : curLoc += 4;
1618 : /* Read the trailing part of the message. */
1619 : /* TDLPack uses 4 bytes for FORTRAN record size, then another 8 bytes for
1620 : * the size of the record (so FORTRAN can see it), then the data rounded
1621 : * up to an 8 byte boundary, then a trailing 4 bytes for a final FORTRAN
1622 : * record size. However it only stores in the gribLen the non-rounded
1623 : * amount, so we need to take care of the rounding, and the trailing 4
1624 : * bytes here. */
1625 0 : pad = ((sInt4) (ceil (tdlpLen / 8.0))) * 8 - tdlpLen + 4;
1626 0 : if (fp.DataSourceFread(buffer, sizeof (char), pad) != pad) {
1627 0 : errSprintf ("Ran out of file\n");
1628 0 : return -1;
1629 : }
1630 0 : return 0;
1631 : }
1632 :
1633 : /*****************************************************************************
1634 : * TDL_ScaleData() --
1635 : *
1636 : * Arthur Taylor / MDL
1637 : *
1638 : * PURPOSE
1639 : * Deal with scaling while excluding the primary and secondary missing
1640 : * values. After this, dst should contain scaled data + primary or secondary
1641 : * missing values
1642 : * "tdlpack library"::pack2d.f line 257 or search for:
1643 : "the above statement"
1644 : *
1645 : * ARGUMENTS
1646 : * Src = The original data. (Input)
1647 : * Dst = The scaled data. (Output)
1648 : * numData = The number of elements in data. (Input)
1649 : * DSF = Decimal Scale Factor for scaling the data. (Input)
1650 : * BSF = Binary Scale Factor for scaling the data. (Input)
1651 : * f_primMiss = Flag saying if we have a primary missing value (In/Out)
1652 : * primMiss = primary missing value. (In/Out)
1653 : * f_secMiss = Flag saying if we have a secondary missing value (In/Out)
1654 : * secMiss = secondary missing value. (In/Out)
1655 : * f_min = Flag saying if we have the minimum value. (Output)
1656 : * min = minimum scaled value in the grid. (Output)
1657 : *
1658 : * FILES/DATABASES: None
1659 : *
1660 : * RETURNS: void
1661 : *
1662 : * HISTORY
1663 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1664 : *
1665 : * NOTES
1666 : *****************************************************************************
1667 : */
1668 : #define SCALE_MISSING 10000L
1669 0 : static void TDL_ScaleData (double *Src, sInt4 *Dst, sInt4 numData,
1670 : int DSF, int BSF, char *f_primMiss,
1671 : double *primMiss, char *f_secMiss,
1672 : double *secMiss, char *f_min, sInt4 *min)
1673 : {
1674 : sInt4 cnt;
1675 0 : double *src = Src;
1676 0 : sInt4 *dst = Dst;
1677 0 : double scale = pow (10.0, -1 * DSF) * pow (2.0, -1 * BSF);
1678 0 : char f_actualPrim = 0;
1679 0 : char f_actualSec = 0;
1680 0 : sInt4 li_primMiss = (sInt4) (*primMiss * SCALE_MISSING + .5);
1681 0 : sInt4 li_secMiss = (sInt4) (*secMiss * SCALE_MISSING + .5);
1682 :
1683 0 : *f_min = 0;
1684 0 : for (cnt = 0; cnt < numData; cnt++) {
1685 0 : if (((*f_primMiss) || (*f_secMiss)) && (*src == *primMiss)) {
1686 0 : *(dst++) = li_primMiss;
1687 0 : src++;
1688 0 : f_actualPrim = 1;
1689 0 : } else if ((*f_secMiss) && (*src == *secMiss)) {
1690 0 : *(dst++) = li_secMiss;
1691 0 : src++;
1692 0 : f_actualSec = 1;
1693 : } else {
1694 0 : *(dst) = (long int) (floor ((*(src++) / scale) + .5));
1695 : /* Check if scaled value == primary missing value. */
1696 0 : if (((*f_primMiss) || (*f_secMiss)) && (*dst == li_primMiss)) {
1697 0 : *dst = *dst - 1;
1698 : }
1699 : /* Check if scaled value == secondary missing value. */
1700 0 : if ((*f_secMiss) && (*dst == li_secMiss)) {
1701 0 : *dst = *dst - 1;
1702 : /* Check if adjustment caused scaled value == primary missing. */
1703 0 : if (*dst == li_primMiss) {
1704 0 : *dst = *dst - 1;
1705 : }
1706 : }
1707 0 : if (!(*f_min)) {
1708 0 : *min = *dst;
1709 0 : *f_min = 1;
1710 0 : } else if (*min > *dst) {
1711 0 : *min = *dst;
1712 : }
1713 0 : dst++;
1714 : }
1715 : }
1716 0 : if ((*f_secMiss) && (!f_actualSec)) {
1717 0 : *f_secMiss = 0;
1718 : }
1719 0 : if (((*f_secMiss) || (*f_primMiss)) && (!f_actualPrim)) {
1720 0 : *f_primMiss = 0;
1721 : /* Check consistency. */
1722 0 : if (*f_secMiss) {
1723 0 : *f_secMiss = 0;
1724 0 : *f_primMiss = 1;
1725 0 : *primMiss = *secMiss;
1726 : }
1727 : }
1728 0 : }
1729 :
1730 : /*****************************************************************************
1731 : * TDL_ReorderGrid() --
1732 : *
1733 : * Arthur Taylor / MDL
1734 : *
1735 : * PURPOSE
1736 : * Loop through the data, so that
1737 : * data is: "a1,1 ... a1,n a2,n ... a2,1 ..."
1738 : * instead of: "a1,1 ... a1,n a2,1 ... a2,n ..."
1739 : *
1740 : * ARGUMENTS
1741 : * Src = The data. (Input/Output)
1742 : * NX = The number of X values. (Input)
1743 : * NY = The number of Y values. (Input)
1744 : *
1745 : * FILES/DATABASES: None
1746 : *
1747 : * RETURNS: void
1748 : *
1749 : * HISTORY
1750 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1751 : *
1752 : * NOTES
1753 : *****************************************************************************
1754 : */
1755 0 : static void TDL_ReorderGrid (sInt4 *Src, short int NX, short int NY)
1756 : {
1757 : int i, j;
1758 : sInt4 *src1, *src2;
1759 : sInt4 li_temp;
1760 :
1761 0 : for (j = 1; j < NY; j += 2) {
1762 0 : src1 = Src + j * NX;
1763 0 : src2 = Src + (j + 1) * NX - 1;
1764 0 : for (i = 0; i < (NX / 2); i++) {
1765 0 : li_temp = *src1;
1766 0 : *(src1++) = *src2;
1767 0 : *(src2--) = li_temp;
1768 : }
1769 : }
1770 0 : }
1771 :
1772 : /*****************************************************************************
1773 : * TDL_GetSecDiff() --
1774 : *
1775 : * Arthur Taylor / MDL
1776 : *
1777 : * PURPOSE
1778 : * Get the second order difference where we have special values for missing,
1779 : * and for actual data we have the following scheme.
1780 : * Data: a1 a2 a3 a4 a5 ...
1781 : * 1st diff: 0 b2 b3 b4 b5 ...
1782 : * 2nd diff: UK1 UK2 c3 c4 c5 ...
1783 : * where UK1 = c3 + b2, and UK2 = c3 + 2 * b2. Note: The choice of UK1, and
1784 : * UK2 doesn't matter because of the following FORTRAN unpacking code:
1785 : * IWORK(1)=IFIRST
1786 : * IWORK(2)=IWORK(1)+IFOD
1787 : * ISUM=IFOD
1788 : * DO 385 K=3,IS4(3)
1789 : * ISUM=IWORK(K)+ISUM
1790 : * IWORK(K)=IWORK(K-1)+ISUM
1791 : * 385 CONTINUE
1792 : * So ISUM is a function of IWORK(3), not IWORK(1).
1793 : *
1794 : * ARGUMENTS
1795 : * Data = The data. (Input)
1796 : * numData = The number of elements in data. (Input)
1797 : * SecDiff = The secondary differences of the data. (Output)
1798 : * f_primMiss = Flag saying if we have a primary missing value (Input)
1799 : * li_primMiss = Scaled primary missing value. (Input)
1800 : * a1 = First non-missing value in the field. (Output)
1801 : * b2 = First non-missing value in the 1st order delta field (Out)
1802 : * min = The minimum value (Input).
1803 : *
1804 : * FILES/DATABASES: None
1805 : *
1806 : * RETURNS: int
1807 : * 0 = Success
1808 : * 1 = Couldn't find second differences (don't use).
1809 : *
1810 : * HISTORY
1811 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1812 : *
1813 : * NOTES
1814 : *****************************************************************************
1815 : */
1816 0 : static int TDL_GetSecDiff (sInt4 *Data, int numData, sInt4 *SecDiff,
1817 : char f_primMiss, sInt4 li_primMiss,
1818 : sInt4 *a1, sInt4 *b2, sInt4 *min)
1819 : {
1820 : int i;
1821 0 : char f_min = 0;
1822 0 : sInt4 last = 0, before_last = 0;
1823 0 : int a1Index = -1;
1824 0 : int a2Index = -1;
1825 :
1826 0 : if (numData < 3) {
1827 0 : return 1;
1828 : }
1829 0 : if (f_primMiss) {
1830 0 : for (i = 0; i < numData; i++) {
1831 0 : if (Data[i] == li_primMiss) {
1832 0 : SecDiff[i] = li_primMiss;
1833 0 : } else if (a1Index == -1) {
1834 0 : a1Index = i;
1835 0 : *a1 = Data[a1Index];
1836 0 : } else if (a2Index == -1) {
1837 0 : a2Index = i;
1838 0 : *b2 = Data[a2Index] - Data[a1Index];
1839 0 : before_last = Data[a1Index];
1840 0 : last = Data[a2Index];
1841 : } else {
1842 0 : SecDiff[i] = Data[i] - 2 * last + before_last;
1843 0 : before_last = last;
1844 0 : last = Data[i];
1845 0 : if (!f_min) {
1846 : /* Set the UK1, UK2 values. */
1847 0 : *min = SecDiff[i];
1848 0 : f_min = 1;
1849 0 : SecDiff[a1Index] = SecDiff[i] + *b2;
1850 0 : SecDiff[a2Index] = SecDiff[i] + 2 * (*b2);
1851 0 : } else if (*min > SecDiff[i]) {
1852 0 : *min = SecDiff[i];
1853 : }
1854 : }
1855 : }
1856 0 : if (!f_min) {
1857 0 : return 1;
1858 : }
1859 : } else {
1860 0 : *a1 = Data[0];
1861 0 : *b2 = Data[1] - Data[0];
1862 0 : for (i = 3; i < numData; i++) {
1863 0 : SecDiff[i] = Data[i] - 2 * Data[i - 1] - Data[i - 2];
1864 0 : if (i == 3) {
1865 0 : *min = SecDiff[i];
1866 : /* Set the UK1, UK2 values. */
1867 0 : SecDiff[0] = SecDiff[i] + *b2;
1868 0 : SecDiff[1] = SecDiff[i] + 2 * (*b2);
1869 0 : } else if (*min > SecDiff[i]) {
1870 0 : *min = SecDiff[i];
1871 : }
1872 : }
1873 : }
1874 0 : return 0;
1875 : }
1876 :
1877 : /*****************************************************************************
1878 : * TDL_UseSecDiff_Prim() --
1879 : *
1880 : * Arthur Taylor / MDL
1881 : *
1882 : * PURPOSE
1883 : * Checks if the average range of 2nd order differences < average range of
1884 : * 0 order differnces, to determine if we should use second order differences
1885 : * This deals with the case when we have primary missing values.
1886 : *
1887 : * ARGUMENTS
1888 : * Data = The data. (Input)
1889 : * numData = The number of elements in data. (Input)
1890 : * SecDiff = The secondary differences of the data. (Input)
1891 : * li_primMiss = Scaled primary missing value. (Input)
1892 : * minGroup = The minimum group size. (Input)
1893 : *
1894 : * FILES/DATABASES: None
1895 : *
1896 : * RETURNS: int
1897 : * 0 = Don't use 2nd order differences.
1898 : * 1 = Use 2nd order differences.
1899 : *
1900 : * HISTORY
1901 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
1902 : *
1903 : * NOTES
1904 : *****************************************************************************
1905 : */
1906 0 : static int TDL_UseSecDiff_Prim (sInt4 *Data, sInt4 numData,
1907 : sInt4 *SecDiff, sInt4 li_primMiss,
1908 : int minGroup)
1909 : {
1910 : int i, locCnt;
1911 : int range0, range2;
1912 : int tot0, tot2;
1913 : char f_min;
1914 0 : sInt4 min = 0, max = 0;
1915 :
1916 0 : locCnt = 0;
1917 0 : range0 = 0;
1918 0 : tot0 = 0;
1919 0 : f_min = 0;
1920 : /* Compute scores for no differences */
1921 0 : for (i = 0; i < numData; i++) {
1922 0 : if (Data[i] != li_primMiss) {
1923 0 : if (!f_min) {
1924 0 : min = Data[i];
1925 0 : max = Data[i];
1926 0 : f_min = 1;
1927 : } else {
1928 0 : if (min > Data[i])
1929 0 : min = Data[i];
1930 0 : if (max < Data[i])
1931 0 : max = Data[i];
1932 : }
1933 : }
1934 0 : locCnt++;
1935 : /* Fake a "group" by using the minimum group size. */
1936 0 : if (locCnt == minGroup) {
1937 0 : if (f_min) {
1938 0 : range0 += (max - min);
1939 0 : tot0++;
1940 0 : f_min = 0;
1941 : }
1942 0 : locCnt = 0;
1943 : }
1944 : }
1945 0 : if (locCnt != 0) {
1946 0 : range0 += (max - min);
1947 0 : tot0++;
1948 : }
1949 :
1950 0 : locCnt = 0;
1951 0 : range2 = 0;
1952 0 : tot2 = 0;
1953 0 : f_min = 0;
1954 : /* Compute scores for second order differences */
1955 0 : for (i = 0; i < numData; i++) {
1956 0 : if (SecDiff[i] != li_primMiss) {
1957 0 : if (!f_min) {
1958 0 : min = SecDiff[i];
1959 0 : max = SecDiff[i];
1960 0 : f_min = 1;
1961 : } else {
1962 0 : if (min > SecDiff[i])
1963 0 : min = SecDiff[i];
1964 0 : if (max < SecDiff[i])
1965 0 : max = SecDiff[i];
1966 : }
1967 : }
1968 0 : locCnt++;
1969 : /* Fake a "group" by using the minimum group size. */
1970 0 : if (locCnt == minGroup) {
1971 0 : if (f_min) {
1972 0 : range2 += (max - min);
1973 0 : tot2++;
1974 0 : f_min = 0;
1975 : }
1976 0 : locCnt = 0;
1977 : }
1978 : }
1979 0 : if (locCnt != 0) {
1980 0 : range2 += (max - min);
1981 0 : tot2++;
1982 : }
1983 :
1984 : /* Compare average group size of no differencing to second order. */
1985 0 : if ((range0 / (tot0 + 0.0)) <= (range2 / (tot2 + 0.0))) {
1986 0 : return 0;
1987 : } else {
1988 0 : return 1;
1989 : }
1990 : }
1991 :
1992 : /*****************************************************************************
1993 : * TDL_UseSecDiff() --
1994 : *
1995 : * Arthur Taylor / MDL
1996 : *
1997 : * PURPOSE
1998 : * Checks if the average range of 2nd order differences < average range of
1999 : * 0 order differnces, to determine if we should use second order differences
2000 : *
2001 : * ARGUMENTS
2002 : * Data = The data. (Input)
2003 : * numData = The number of elements in data. (Input)
2004 : * SecDiff = The secondary differences of the data. (Input)
2005 : * minGroup = The minimum group size. (Input)
2006 : *
2007 : * FILES/DATABASES: None
2008 : *
2009 : * RETURNS: int
2010 : * 0 = Don't use 2nd order differences.
2011 : * 1 = Use 2nd order differences.
2012 : *
2013 : * HISTORY
2014 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2015 : *
2016 : * NOTES
2017 : *****************************************************************************
2018 : */
2019 0 : static int TDL_UseSecDiff (sInt4 *Data, sInt4 numData,
2020 : sInt4 *SecDiff, int minGroup)
2021 : {
2022 : int i, locCnt;
2023 : int range0, range2;
2024 : int tot0, tot2;
2025 0 : sInt4 min = 0, max = 0;
2026 :
2027 0 : locCnt = 0;
2028 0 : range0 = 0;
2029 0 : tot0 = 0;
2030 : /* Compute scores for no differences */
2031 0 : for (i = 0; i < numData; i++) {
2032 0 : if (locCnt == 0) {
2033 0 : min = Data[i];
2034 0 : max = Data[i];
2035 : } else {
2036 0 : if (min > Data[i])
2037 0 : min = Data[i];
2038 0 : if (max < Data[i])
2039 0 : max = Data[i];
2040 : }
2041 0 : locCnt++;
2042 : /* Fake a "group" by using the minimum group size. */
2043 0 : if (locCnt == minGroup) {
2044 0 : range0 += (max - min);
2045 0 : tot0++;
2046 0 : locCnt = 0;
2047 : }
2048 : }
2049 0 : if (locCnt != 0) {
2050 0 : range0 += (max - min);
2051 0 : tot0++;
2052 : }
2053 :
2054 0 : locCnt = 0;
2055 0 : range2 = 0;
2056 0 : tot2 = 0;
2057 : /* Compute scores for second order differences */
2058 0 : for (i = 0; i < numData; i++) {
2059 0 : if (locCnt == 0) {
2060 0 : min = SecDiff[i];
2061 0 : max = SecDiff[i];
2062 : } else {
2063 0 : if (min > SecDiff[i])
2064 0 : min = SecDiff[i];
2065 0 : if (max < SecDiff[i])
2066 0 : max = SecDiff[i];
2067 : }
2068 0 : locCnt++;
2069 : /* Fake a "group" by using the minimum group size. */
2070 0 : if (locCnt == minGroup) {
2071 0 : range2 += (max - min);
2072 0 : tot2++;
2073 0 : locCnt = 0;
2074 : }
2075 : }
2076 0 : if (locCnt != 0) {
2077 0 : range2 += (max - min);
2078 0 : tot2++;
2079 : }
2080 :
2081 : /* Compare average group size of no differencing to second order. */
2082 0 : if ((range0 / (tot0 + 0.0)) <= (range2 / (tot2 + 0.0))) {
2083 0 : return 0;
2084 : } else {
2085 0 : return 1;
2086 : }
2087 : }
2088 :
2089 : /*****************************************************************************
2090 : * power() --
2091 : *
2092 : * Arthur Taylor / MDL
2093 : *
2094 : * PURPOSE
2095 : * Calculate the number of bits required to store a given positive number.
2096 : *
2097 : * ARGUMENTS
2098 : * val = The number to store (Input)
2099 : * extra = number of slots to allocate for prim/sec missing values (Input)
2100 : *
2101 : * FILES/DATABASES: None
2102 : *
2103 : * RETURNS: int
2104 : * The number of bits needed to store this number.
2105 : *
2106 : * HISTORY
2107 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2108 : *
2109 : * NOTES
2110 : *****************************************************************************
2111 : */
2112 0 : static int power (uInt4 val, int extra)
2113 : {
2114 : int i;
2115 :
2116 0 : val += extra;
2117 0 : if (val == 0) {
2118 0 : return 1;
2119 : }
2120 0 : for (i = 0; val != 0; i++) {
2121 0 : val = val >> 1;
2122 : }
2123 0 : return i;
2124 : }
2125 :
2126 : /*****************************************************************************
2127 : * findMaxMin2() --
2128 : *
2129 : * Arthur Taylor / MDL
2130 : *
2131 : * PURPOSE
2132 : * Find the min/max value between start/stop index values in the data
2133 : * Assuming primary and secondary missing values.
2134 : *
2135 : * ARGUMENTS
2136 : * Data = The data. (Input)
2137 : * start = The starting index in data (Input)
2138 : * stop = The stoping index in data (Input)
2139 : * li_primMiss = scaled primary missing value (Input)
2140 : * li_secMiss = scaled secondary missing value (Input)
2141 : * min = The min value found (Output)
2142 : * max = The max value found (Output)
2143 : *
2144 : * FILES/DATABASES: None
2145 : *
2146 : * RETURNS: char
2147 : * Flag if min/max are valid.
2148 : *
2149 : * HISTORY
2150 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2151 : *
2152 : * NOTES
2153 : *****************************************************************************
2154 : */
2155 0 : static char findMaxMin2 (sInt4 *Data, int start, int stop,
2156 : sInt4 li_primMiss, sInt4 li_secMiss,
2157 : sInt4 *min, sInt4 *max)
2158 : {
2159 0 : char f_min = 0; /* Flag if we found the max/min values */
2160 : int i; /* Loop counter. */
2161 :
2162 0 : *max = *min = Data[start];
2163 0 : for (i = start; i < stop; i++) {
2164 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2165 0 : if (!f_min) {
2166 0 : *max = Data[i];
2167 0 : *min = Data[i];
2168 0 : f_min = 1;
2169 : } else {
2170 0 : if (*max < Data[i]) {
2171 0 : *max = Data[i];
2172 0 : } else if (*min > Data[i]) {
2173 0 : *min = Data[i];
2174 : }
2175 : }
2176 : }
2177 : }
2178 0 : return f_min;
2179 : }
2180 :
2181 : /*****************************************************************************
2182 : * findMaxMin1() --
2183 : *
2184 : * Arthur Taylor / MDL
2185 : *
2186 : * PURPOSE
2187 : * Find the min/max value between start/stop index values in the data
2188 : * Assuming primary missing values.
2189 : *
2190 : * ARGUMENTS
2191 : * Data = The data. (Input)
2192 : * start = The starting index in data (Input)
2193 : * stop = The stoping index in data (Input)
2194 : * li_primMiss = scaled primary missing value (Input)
2195 : * min = The min value found (Output)
2196 : * max = The max value found (Output)
2197 : *
2198 : * FILES/DATABASES: None
2199 : *
2200 : * RETURNS: char
2201 : * Flag if min/max are valid.
2202 : *
2203 : * HISTORY
2204 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2205 : *
2206 : * NOTES
2207 : *****************************************************************************
2208 : */
2209 0 : static char findMaxMin1 (sInt4 *Data, int start, int stop,
2210 : sInt4 li_primMiss, sInt4 *min, sInt4 *max)
2211 : {
2212 0 : char f_min = 0; /* Flag if we found the max/min values */
2213 : int i; /* Loop counter. */
2214 :
2215 0 : *max = *min = Data[start];
2216 0 : for (i = start; i < stop; i++) {
2217 0 : if (Data[i] != li_primMiss) {
2218 0 : if (!f_min) {
2219 0 : *max = Data[i];
2220 0 : *min = Data[i];
2221 0 : f_min = 1;
2222 : } else {
2223 0 : if (*max < Data[i]) {
2224 0 : *max = Data[i];
2225 0 : } else if (*min > Data[i]) {
2226 0 : *min = Data[i];
2227 : }
2228 : }
2229 : }
2230 : }
2231 0 : return f_min;
2232 : }
2233 :
2234 : /*****************************************************************************
2235 : * findMaxMin0() --
2236 : *
2237 : * Arthur Taylor / MDL
2238 : *
2239 : * PURPOSE
2240 : * Find the min/max value between start/stop index values in the data
2241 : * Assuming no missing values.
2242 : *
2243 : * ARGUMENTS
2244 : * Data = The data. (Input)
2245 : * start = The starting index in data (Input)
2246 : * stop = The stoping index in data (Input)
2247 : * min = The min value found (Output)
2248 : * max = The max value found (Output)
2249 : *
2250 : * FILES/DATABASES: None
2251 : *
2252 : * RETURNS: void
2253 : *
2254 : * HISTORY
2255 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
2256 : *
2257 : * NOTES
2258 : *****************************************************************************
2259 : */
2260 0 : static void findMaxMin0 (sInt4 *Data, int start, int stop, sInt4 *min,
2261 : sInt4 *max)
2262 : {
2263 : int i; /* Loop counter. */
2264 :
2265 0 : *max = *min = Data[start];
2266 0 : for (i = start + 1; i < stop; i++) {
2267 0 : if (*max < Data[i]) {
2268 0 : *max = Data[i];
2269 0 : } else if (*min > Data[i]) {
2270 0 : *min = Data[i];
2271 : }
2272 : }
2273 0 : }
2274 :
2275 : /*****************************************************************************
2276 : * findGroup2() --
2277 : *
2278 : * Arthur Taylor / MDL
2279 : *
2280 : * PURPOSE
2281 : * Find "split" so that the numbers between start and split are within
2282 : * "range" of each other... stops if it reaches "stop".
2283 : * Assumes primary and secondary missing values.
2284 : *
2285 : * ARGUMENTS
2286 : * Data = The data. (Input)
2287 : * start = The starting index in data (Input)
2288 : * stop = The stoping index in data (Input)
2289 : * li_primMiss = scaled primary missing value (Input)
2290 : * li_secMiss = scaled secondary missing value (Input)
2291 : * range = The range to use (Input)
2292 : * split = The first index that is out of the range (Output)
2293 : * min = The min value for the group. (Output)
2294 : * max = The max value for the group. (Output)
2295 : *
2296 : * FILES/DATABASES: None
2297 : *
2298 : * RETURNS: void
2299 : *
2300 : * HISTORY
2301 : * 1/2005 Arthur Taylor (MDL): Created
2302 : *
2303 : * NOTES
2304 : *****************************************************************************
2305 : */
2306 0 : static void findGroup2 (sInt4 *Data, int start, int stop,
2307 : sInt4 li_primMiss, sInt4 li_secMiss,
2308 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2309 : {
2310 0 : char f_min = 0; /* Flag if we found the max/min values */
2311 : int i; /* Loop counter. */
2312 :
2313 0 : *min = *max = 0;
2314 0 : for (i = start; i < stop; i++) {
2315 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2316 0 : if (!f_min) {
2317 0 : *max = *min = Data[i];
2318 0 : f_min = 1;
2319 : } else {
2320 0 : if (*max < Data[i]) {
2321 0 : if ((Data[i] - *min) > range) {
2322 0 : *split = i;
2323 0 : return;
2324 : }
2325 0 : *max = Data[i];
2326 0 : } else if (*min > Data[i]) {
2327 0 : if ((*max - Data[i]) > range) {
2328 0 : *split = i;
2329 0 : return;
2330 : }
2331 0 : *min = Data[i];
2332 : }
2333 : }
2334 : }
2335 : }
2336 0 : *split = stop;
2337 : }
2338 :
2339 : /*****************************************************************************
2340 : * findGroup1() --
2341 : *
2342 : * Arthur Taylor / MDL
2343 : *
2344 : * PURPOSE
2345 : * Find "split" so that the numbers between start and split are within
2346 : * "range" of each other... stops if it reaches "stop".
2347 : * Assumes primary missing values.
2348 : *
2349 : * ARGUMENTS
2350 : * Data = The data. (Input)
2351 : * start = The starting index in data (Input)
2352 : * stop = The stoping index in data (Input)
2353 : * li_primMiss = scaled primary missing value (Input)
2354 : * range = The range to use (Input)
2355 : * split = The first index that is out of the range (Output)
2356 : * min = The min value for the group. (Output)
2357 : * max = The max value for the group. (Output)
2358 : *
2359 : * FILES/DATABASES: None
2360 : *
2361 : * RETURNS: void
2362 : *
2363 : * HISTORY
2364 : * 1/2005 Arthur Taylor (MDL): Created
2365 : *
2366 : * NOTES
2367 : *****************************************************************************
2368 : */
2369 0 : static void findGroup1 (sInt4 *Data, int start, int stop,
2370 : sInt4 li_primMiss, sInt4 range, int *split,
2371 : sInt4 *min, sInt4 *max)
2372 : {
2373 0 : char f_min = 0; /* Flag if we found the max/min values */
2374 : int i; /* Loop counter. */
2375 :
2376 0 : *min = *max = 0;
2377 0 : for (i = start; i < stop; i++) {
2378 0 : if (Data[i] != li_primMiss) {
2379 0 : if (!f_min) {
2380 0 : *max = *min = Data[i];
2381 0 : f_min = 1;
2382 : } else {
2383 0 : if (*max < Data[i]) {
2384 0 : if ((Data[i] - *min) > range) {
2385 0 : *split = i;
2386 0 : return;
2387 : }
2388 0 : *max = Data[i];
2389 0 : } else if (*min > Data[i]) {
2390 0 : if ((*max - Data[i]) > range) {
2391 0 : *split = i;
2392 0 : return;
2393 : }
2394 0 : *min = Data[i];
2395 : }
2396 : }
2397 : }
2398 : }
2399 0 : *split = stop;
2400 : }
2401 :
2402 : /*****************************************************************************
2403 : * findGroup0() --
2404 : *
2405 : * Arthur Taylor / MDL
2406 : *
2407 : * PURPOSE
2408 : * Find "split" so that the numbers between start and split are within
2409 : * "range" of each other... stops if it reaches "stop".
2410 : * Assumes no missing values.
2411 : *
2412 : * ARGUMENTS
2413 : * Data = The data. (Input)
2414 : * start = The starting index in data (Input)
2415 : * stop = The stoping index in data (Input)
2416 : * range = The range to use (Input)
2417 : * split = The first index that is out of the range (Output)
2418 : * min = The min value for the group. (Output)
2419 : * max = The max value for the group. (Output)
2420 : *
2421 : * FILES/DATABASES: None
2422 : *
2423 : * RETURNS: void
2424 : *
2425 : * HISTORY
2426 : * 1/2005 Arthur Taylor (MDL): Created
2427 : *
2428 : * NOTES
2429 : *****************************************************************************
2430 : */
2431 0 : static void findGroup0 (sInt4 *Data, int start, int stop,
2432 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2433 : {
2434 : int i; /* Loop counter. */
2435 :
2436 0 : *max = *min = Data[0];
2437 0 : for (i = start + 1; i < stop; i++) {
2438 0 : if (*max < Data[i]) {
2439 0 : if ((Data[i] - *min) > range) {
2440 0 : *split = i;
2441 0 : return;
2442 : }
2443 0 : *max = Data[i];
2444 0 : } else if (*min > Data[i]) {
2445 0 : if ((*max - Data[i]) > range) {
2446 0 : *split = i;
2447 0 : return;
2448 : }
2449 0 : *min = Data[i];
2450 : }
2451 : }
2452 0 : *split = stop;
2453 : }
2454 :
2455 : /*****************************************************************************
2456 : * findGroupRev2() --
2457 : *
2458 : * Arthur Taylor / MDL
2459 : *
2460 : * PURPOSE
2461 : * Find "split" so that the numbers between split and stop are within
2462 : * "range" of each other... stops if it reaches "start".
2463 : * Assumes primary and secondary missing values.
2464 : *
2465 : * ARGUMENTS
2466 : * Data = The data. (Input)
2467 : * start = The starting index in data (Input)
2468 : * stop = The stoping index in data (Input)
2469 : * li_primMiss = scaled primary missing value (Input)
2470 : * li_secMiss = scaled secondary missing value (Input)
2471 : * range = The range to use (Input)
2472 : * split = The first index that is still in the range (Output)
2473 : * min = The min value for the group. (Output)
2474 : * max = The max value for the group. (Output)
2475 : *
2476 : * FILES/DATABASES: None
2477 : *
2478 : * RETURNS: void
2479 : *
2480 : * HISTORY
2481 : * 1/2005 Arthur Taylor (MDL): Created
2482 : *
2483 : * NOTES
2484 : *****************************************************************************
2485 : */
2486 0 : static void findGroupRev2 (sInt4 *Data, int start, int stop,
2487 : sInt4 li_primMiss, sInt4 li_secMiss,
2488 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2489 : {
2490 0 : char f_min = 0; /* Flag if we found the max/min values */
2491 : int i; /* Loop counter. */
2492 :
2493 0 : *min = *max = 0;
2494 0 : for (i = stop - 1; i >= start; i--) {
2495 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
2496 0 : if (!f_min) {
2497 0 : *max = *min = Data[i];
2498 0 : f_min = 1;
2499 : } else {
2500 0 : if (*max < Data[i]) {
2501 0 : if ((Data[i] - *min) > range) {
2502 0 : *split = i + 1;
2503 0 : return;
2504 : }
2505 0 : *max = Data[i];
2506 0 : } else if (*min > Data[i]) {
2507 0 : if ((*max - Data[i]) > range) {
2508 0 : *split = i + 1;
2509 0 : return;
2510 : }
2511 0 : *min = Data[i];
2512 : }
2513 : }
2514 : }
2515 : }
2516 0 : *split = start;
2517 : }
2518 :
2519 : /*****************************************************************************
2520 : * findGroupRev1() --
2521 : *
2522 : * Arthur Taylor / MDL
2523 : *
2524 : * PURPOSE
2525 : * Find "split" so that the numbers between split and stop are within
2526 : * "range" of each other... stops if it reaches "start".
2527 : * Assumes primary missing values.
2528 : *
2529 : * ARGUMENTS
2530 : * Data = The data. (Input)
2531 : * start = The starting index in data (Input)
2532 : * stop = The stoping index in data (Input)
2533 : * li_primMiss = scaled primary missing value (Input)
2534 : * range = The range to use (Input)
2535 : * split = The first index that is still in the range (Output)
2536 : * min = The min value for the group. (Output)
2537 : * max = The max value for the group. (Output)
2538 : *
2539 : * FILES/DATABASES: None
2540 : *
2541 : * RETURNS: void
2542 : *
2543 : * HISTORY
2544 : * 1/2005 Arthur Taylor (MDL): Created
2545 : *
2546 : * NOTES
2547 : *****************************************************************************
2548 : */
2549 0 : static void findGroupRev1 (sInt4 *Data, int start, int stop,
2550 : sInt4 li_primMiss, sInt4 range, int *split,
2551 : sInt4 *min, sInt4 *max)
2552 : {
2553 0 : char f_min = 0; /* Flag if we found the max/min values */
2554 : int i; /* Loop counter. */
2555 :
2556 0 : *min = *max = 0;
2557 0 : for (i = stop - 1; i >= start; i--) {
2558 0 : if (Data[i] != li_primMiss) {
2559 0 : if (!f_min) {
2560 0 : *max = *min = Data[i];
2561 0 : f_min = 1;
2562 : } else {
2563 0 : if (*max < Data[i]) {
2564 0 : if ((Data[i] - *min) > range) {
2565 0 : *split = i + 1;
2566 0 : return;
2567 : }
2568 0 : *max = Data[i];
2569 0 : } else if (*min > Data[i]) {
2570 0 : if ((*max - Data[i]) > range) {
2571 0 : *split = i + 1;
2572 0 : return;
2573 : }
2574 0 : *min = Data[i];
2575 : }
2576 : }
2577 : }
2578 : }
2579 0 : *split = start;
2580 : }
2581 :
2582 : /*****************************************************************************
2583 : * findGroupRev0() --
2584 : *
2585 : * Arthur Taylor / MDL
2586 : *
2587 : * PURPOSE
2588 : * Find "split" so that the numbers between split and stop are within
2589 : * "range" of each other... stops if it reaches "start".
2590 : * Assumes no missing values.
2591 : *
2592 : * ARGUMENTS
2593 : * Data = The data. (Input)
2594 : * start = The starting index in data (Input)
2595 : * stop = The stoping index in data (Input)
2596 : * li_primMiss = scaled primary missing value (Input)
2597 : * range = The range to use (Input)
2598 : * split = The first index that is still in the range (Output)
2599 : * min = The min value for the group. (Output)
2600 : * max = The max value for the group. (Output)
2601 : *
2602 : * FILES/DATABASES: None
2603 : *
2604 : * RETURNS: void
2605 : *
2606 : * HISTORY
2607 : * 1/2005 Arthur Taylor (MDL): Created
2608 : *
2609 : * NOTES
2610 : *****************************************************************************
2611 : */
2612 0 : static void findGroupRev0 (sInt4 *Data, int start, int stop,
2613 : sInt4 range, int *split, sInt4 *min, sInt4 *max)
2614 : {
2615 : int i; /* Loop counter. */
2616 :
2617 0 : *max = *min = Data[stop - 1];
2618 0 : for (i = stop - 2; i >= start; i--) {
2619 0 : if (*max < Data[i]) {
2620 0 : if ((Data[i] - *min) > range) {
2621 0 : *split = i + 1;
2622 0 : return;
2623 : }
2624 0 : *max = Data[i];
2625 0 : } else if (*min > Data[i]) {
2626 0 : if ((*max - Data[i]) > range) {
2627 0 : *split = i + 1;
2628 0 : return;
2629 : }
2630 0 : *min = Data[i];
2631 : }
2632 : }
2633 0 : *split = start;
2634 : }
2635 :
2636 : /*****************************************************************************
2637 : * shiftGroup2() --
2638 : *
2639 : * Arthur Taylor / MDL
2640 : *
2641 : * PURPOSE
2642 : * Find "split" so that the numbers between split and start1 are all inside
2643 : * the range defined by max, min and bit. It allows max and min to change,
2644 : * as long as it doesn't exceed the range defined by bit.
2645 : * This is very similar to findGroupRev?() but here we already know
2646 : * information about the min/max values, and are just expanding the group a
2647 : * little, while the other one knew nothing about the group, and just wanted
2648 : * a group of a given range.
2649 : * Assumes primary and secondary missing values.
2650 : *
2651 : * ARGUMENTS
2652 : * Data = The data. (Input)
2653 : * start1 = The starting index in data (Input)
2654 : * start2 = The starting index of the earlier group (ie don't go to any
2655 : * earlier indicies than this. (Input)
2656 : * li_primMiss = scaled primary missing value (Input)
2657 : * li_secMiss = scaled secondary missing value (Input)
2658 : * bit = The range we are allowed to store this in. (Input)
2659 : * min = The min value for the group. (Input/Output)
2660 : * max = The max value for the group. (Input/Output)
2661 : * split = The first index that is still in the range (Output)
2662 : *
2663 : * FILES/DATABASES: None
2664 : *
2665 : * RETURNS: void
2666 : *
2667 : * HISTORY
2668 : * 1/2005 Arthur Taylor (MDL): Created
2669 : *
2670 : * NOTES
2671 : *****************************************************************************
2672 : */
2673 0 : static void shiftGroup2 (sInt4 *Data, int start1, int start2,
2674 : sInt4 li_primMiss, sInt4 li_secMiss, int bit,
2675 : sInt4 *min, sInt4 *max, size_t *split)
2676 : {
2677 : int i; /* Loop counter. */
2678 : int range; /* The range defined by bit. */
2679 :
2680 0 : range = (int) (pow (2.0, bit) - 1) - 1;
2681 0 : myAssert (start2 <= start1);
2682 0 : for (i = start1; i >= start2; i--) {
2683 0 : if ((Data[i] != li_primMiss) && (Data[i] != li_secMiss)) {
2684 0 : if (Data[i] > *max) {
2685 0 : if ((Data[i] - *min) <= range) {
2686 0 : *max = Data[i];
2687 : } else {
2688 0 : *split = i + 1;
2689 0 : return;
2690 : }
2691 0 : } else if (Data[i] < *min) {
2692 0 : if ((*max - Data[i]) <= range) {
2693 0 : *min = Data[i];
2694 : } else {
2695 0 : *split = i + 1;
2696 0 : return;
2697 : }
2698 : }
2699 : }
2700 : }
2701 0 : *split = start2;
2702 : }
2703 :
2704 : /*****************************************************************************
2705 : * shiftGroup1() --
2706 : *
2707 : * Arthur Taylor / MDL
2708 : *
2709 : * PURPOSE
2710 : * Find "split" so that the numbers between split and start1 are all inside
2711 : * the range defined by max, min and bit. It allows max and min to change,
2712 : * as long as it doesn't exceed the range defined by bit.
2713 : * This is very similar to findGroupRev?() but here we already know
2714 : * information about the min/max values, and are just expanding the group a
2715 : * little, while the other one knew nothing about the group, and just wanted
2716 : * a group of a given range.
2717 : * Assumes primary missing values.
2718 : *
2719 : * ARGUMENTS
2720 : * Data = The data. (Input)
2721 : * start1 = The starting index in data (Input)
2722 : * start2 = The starting index of the earlier group (ie don't go to any
2723 : * earlier indicies than this. (Input)
2724 : * li_primMiss = scaled primary missing value (Input)
2725 : * bit = The range we are allowed to store this in. (Input)
2726 : * min = The min value for the group. (Input/Output)
2727 : * max = The max value for the group. (Input/Output)
2728 : * split = The first index that is still in the range (Output)
2729 : *
2730 : * FILES/DATABASES: None
2731 : *
2732 : * RETURNS: void
2733 : *
2734 : * HISTORY
2735 : * 1/2005 Arthur Taylor (MDL): Created
2736 : *
2737 : * NOTES
2738 : *****************************************************************************
2739 : */
2740 0 : static void shiftGroup1 (sInt4 *Data, int start1, int start2,
2741 : sInt4 li_primMiss, int bit,
2742 : sInt4 *min, sInt4 *max, size_t *split)
2743 : {
2744 : int i; /* Loop counter. */
2745 : int range; /* The range defined by bit. */
2746 :
2747 0 : range = (int) (pow (2.0, bit) - 1) - 1;
2748 0 : myAssert (start2 <= start1);
2749 0 : for (i = start1; i >= start2; i--) {
2750 0 : if (Data[i] != li_primMiss) {
2751 0 : if (Data[i] > *max) {
2752 0 : if ((Data[i] - *min) <= range) {
2753 0 : *max = Data[i];
2754 : } else {
2755 0 : *split = i + 1;
2756 0 : return;
2757 : }
2758 0 : } else if (Data[i] < *min) {
2759 0 : if ((*max - Data[i]) <= range) {
2760 0 : *min = Data[i];
2761 : } else {
2762 0 : *split = i + 1;
2763 0 : return;
2764 : }
2765 : }
2766 : }
2767 : }
2768 0 : *split = start2;
2769 : }
2770 :
2771 : /*****************************************************************************
2772 : * shiftGroup0() --
2773 : *
2774 : * Arthur Taylor / MDL
2775 : *
2776 : * PURPOSE
2777 : * Find "split" so that the numbers between split and start1 are all inside
2778 : * the range defined by max, min and bit. It allows max and min to change,
2779 : * as long as it doesn't exceed the range defined by bit.
2780 : * This is very similar to findGroupRev?() but here we already know
2781 : * information about the min/max values, and are just expanding the group a
2782 : * little, while the other one knew nothing about the group, and just wanted
2783 : * a group of a given range.
2784 : * Assumes no missing values.
2785 : *
2786 : * ARGUMENTS
2787 : * Data = The data. (Input)
2788 : * start1 = The starting index in data (Input)
2789 : * start2 = The starting index of the earlier group (ie don't go to any
2790 : * earlier indicies than this. (Input)
2791 : * bit = The range we are allowed to store this in. (Input)
2792 : * min = The min value for the group. (Input/Output)
2793 : * max = The max value for the group. (Input/Output)
2794 : * split = The first index that is still in the range (Output)
2795 : *
2796 : * FILES/DATABASES: None
2797 : *
2798 : * RETURNS: void
2799 : *
2800 : * HISTORY
2801 : * 1/2005 Arthur Taylor (MDL): Created
2802 : *
2803 : * NOTES
2804 : *****************************************************************************
2805 : */
2806 0 : static void shiftGroup0 (sInt4 *Data, int start1, int start2, int bit,
2807 : sInt4 *min, sInt4 *max, size_t *split)
2808 : {
2809 : int i; /* Loop counter. */
2810 : int range; /* The range defined by bit. */
2811 :
2812 0 : range = (int) (pow (2.0, bit) - 1) - 0;
2813 0 : myAssert (start2 <= start1);
2814 0 : for (i = start1; i >= start2; i--) {
2815 0 : if (Data[i] > *max) {
2816 0 : if ((Data[i] - *min) <= range) {
2817 0 : *max = Data[i];
2818 : } else {
2819 0 : *split = i + 1;
2820 0 : return;
2821 : }
2822 0 : } else if (Data[i] < *min) {
2823 0 : if ((*max - Data[i]) <= range) {
2824 0 : *min = Data[i];
2825 : } else {
2826 0 : *split = i + 1;
2827 0 : return;
2828 : }
2829 : }
2830 : }
2831 0 : *split = start2;
2832 : }
2833 :
2834 : /*****************************************************************************
2835 : * doSplit() --
2836 : *
2837 : * Arthur Taylor / MDL
2838 : *
2839 : * PURPOSE
2840 : * Reduce the "bit range", and create groups that grab as much as they can
2841 : * to the right. Then reduce those groups if they improve the score.
2842 : *
2843 : * ARGUMENTS
2844 : * Data = The data. (Input)
2845 : * numData = The number of elements in data. (Input)
2846 : * G = The group to split. (Input)
2847 : * lclGroup = The resulting groups (Output)
2848 : * numLclGroup = The number of resulting groups. (Output)
2849 : * f_primMiss = Flag if we have a primary missing value (Input)
2850 : * li_primMiss = scaled primary missing value (Input)
2851 : * f_secMiss = Flag if we have a secondary missing value (Input)
2852 : * li_secMiss = scaled secondary missing value (Input)
2853 : * xFactor = Estimate of cost (in bits) of a group. (Input)
2854 : *
2855 : * FILES/DATABASES: None
2856 : *
2857 : * RETURNS: void
2858 : *
2859 : * HISTORY
2860 : * 1/2005 Arthur Taylor (MDL): Created.
2861 : *
2862 : * NOTES
2863 : *****************************************************************************
2864 : */
2865 0 : static void doSplit (sInt4 *Data, int numData, TDLGroupType * G,
2866 : TDLGroupType ** lclGroup, int *numLclGroup,
2867 : char f_primMiss, sInt4 li_primMiss,
2868 : char f_secMiss, sInt4 li_secMiss, int xFactor)
2869 : {
2870 : int start; /* Where to start the current group. */
2871 : int range; /* The range to make the groups. */
2872 : int final; /* One more than the last index in the group G. */
2873 : int split; /* Where to split the group. */
2874 : TDLGroupType G1; /* The current group to add. */
2875 : TDLGroupType G2; /* The group if we evaporated the previous group. */
2876 : int evaporate; /* How many groups we have "evaporated". */
2877 : int i; /* Loop counter to help "evaporate" groups. */
2878 : sInt4 scoreA; /* The original score for 2 groups */
2879 : sInt4 scoreB; /* The new score (having evaporated a group) */
2880 : int GroupLen; /* Actual alloc'ed group len. */
2881 :
2882 : /* *INDENT-OFF* */
2883 : /* The (pow (2, ..) -1) is because 2^n - 1 is max range of a group.
2884 : * Example n = 1, 2^1 -1 = 1 range of (1,0) is 1.
2885 : * Example n = 2, 2^2 -1 = 3 range of (3,2,1,0) is 3.
2886 : * The G.bit - 1 is because we are trying to reduce the range. */
2887 : /* *INDENT-ON* */
2888 0 : range = (int) (pow (2.0, G->bit - 1) - 1) - (f_secMiss + f_primMiss);
2889 0 : split = G->start;
2890 0 : start = G->start;
2891 0 : final = G->start + G->num;
2892 0 : myAssert (final <= numData);
2893 0 : *numLclGroup = 0;
2894 0 : GroupLen = 1;
2895 0 : *lclGroup = (TDLGroupType *) malloc (GroupLen * sizeof (TDLGroupType));
2896 0 : while (split < final) {
2897 0 : if (f_secMiss) {
2898 : findGroup2 (Data, start, final, li_primMiss, li_secMiss, range,
2899 0 : &split, &(G1.min), &(G1.max));
2900 0 : } else if (f_primMiss) {
2901 : findGroup1 (Data, start, final, li_primMiss, range, &split,
2902 0 : &(G1.min), &(G1.max));
2903 : } else {
2904 0 : findGroup0 (Data, start, final, range, &split, &(G1.min), &(G1.max));
2905 : }
2906 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
2907 0 : f_secMiss + f_primMiss);
2908 0 : G1.num = split - start;
2909 0 : G1.start = start;
2910 0 : G1.f_trySplit = 1;
2911 0 : G1.f_tryShift = 1;
2912 : /* Test if we should add to previous group, or create a new group. */
2913 0 : if (*numLclGroup == 0) {
2914 0 : *numLclGroup = 1;
2915 0 : (*lclGroup)[0] = G1;
2916 : } else {
2917 0 : G2.start = (*lclGroup)[*numLclGroup - 1].start;
2918 0 : G2.num = (*lclGroup)[*numLclGroup - 1].num + G1.num;
2919 0 : G2.min = ((*lclGroup)[*numLclGroup - 1].min < G1.min) ?
2920 0 : (*lclGroup)[*numLclGroup - 1].min : G1.min;
2921 0 : G2.max = ((*lclGroup)[*numLclGroup - 1].max > G1.max) ?
2922 0 : (*lclGroup)[*numLclGroup - 1].max : G1.max;
2923 : G2.bit = (char) power ((uInt4) (G2.max - G2.min),
2924 0 : f_secMiss + f_primMiss);
2925 0 : G2.f_trySplit = 1;
2926 0 : G2.f_tryShift = 1;
2927 0 : scoreA = ((*lclGroup)[*numLclGroup - 1].bit *
2928 0 : (*lclGroup)[*numLclGroup - 1].num) + xFactor;
2929 0 : scoreA += G1.bit * G1.num + xFactor;
2930 0 : scoreB = G2.bit * G2.num + xFactor;
2931 0 : if (scoreB < scoreA) {
2932 0 : (*lclGroup)[*numLclGroup - 1] = G2;
2933 : /* See if we can evaporate any of the old groups */
2934 0 : evaporate = 0;
2935 0 : for (i = *numLclGroup - 1; i > 0; i--) {
2936 0 : G1.start = (*lclGroup)[i - 1].start;
2937 0 : G1.num = (*lclGroup)[i].num + (*lclGroup)[i - 1].num;
2938 0 : G1.min = ((*lclGroup)[i].min < (*lclGroup)[i - 1].min) ?
2939 0 : (*lclGroup)[i].min : (*lclGroup)[i - 1].min;
2940 0 : G1.max = ((*lclGroup)[i].max > (*lclGroup)[i - 1].max) ?
2941 0 : (*lclGroup)[i].max : (*lclGroup)[i - 1].max;
2942 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
2943 0 : f_secMiss + f_primMiss);
2944 0 : G1.f_trySplit = 1;
2945 0 : G1.f_tryShift = 1;
2946 0 : scoreA = (*lclGroup)[i].bit * (*lclGroup)[i].num + xFactor;
2947 0 : scoreA += ((*lclGroup)[i - 1].bit * (*lclGroup)[i - 1].num +
2948 0 : xFactor);
2949 0 : scoreB = G1.bit * G1.num + xFactor;
2950 0 : if (scoreB < scoreA) {
2951 0 : evaporate++;
2952 0 : (*lclGroup)[i - 1] = G1;
2953 : } else {
2954 0 : break;
2955 : }
2956 : }
2957 0 : if (evaporate != 0) {
2958 0 : *numLclGroup = *numLclGroup - evaporate;
2959 : }
2960 : } else {
2961 0 : *numLclGroup = *numLclGroup + 1;
2962 0 : if (*numLclGroup > GroupLen) {
2963 0 : GroupLen = *numLclGroup;
2964 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
2965 : GroupLen *
2966 0 : sizeof (TDLGroupType));
2967 : }
2968 0 : (*lclGroup)[*numLclGroup - 1] = G1;
2969 : }
2970 : }
2971 0 : start = split;
2972 : }
2973 0 : }
2974 :
2975 : /*****************************************************************************
2976 : * doSplitRight() --
2977 : *
2978 : * Arthur Taylor / MDL
2979 : *
2980 : * PURPOSE
2981 : * Break into two groups right has range n - 1, left has range n.
2982 : *
2983 : * ARGUMENTS
2984 : * Data = The data. (Input)
2985 : * numData = The number of elements in data. (Input)
2986 : * G = The group to split. (Input)
2987 : * G1 = The right most group (Output)
2988 : * G2 = The remainder. (Output)
2989 : * f_primMiss = Flag if we have a primary missing value (Input)
2990 : * li_primMiss = scaled primary missing value (Input)
2991 : * f_secMiss = Flag if we have a secondary missing value (Input)
2992 : * li_secMiss = scaled secondary missing value (Input)
2993 : *
2994 : * FILES/DATABASES: None
2995 : *
2996 : * RETURNS: void
2997 : *
2998 : * HISTORY
2999 : * 1/2005 Arthur Taylor (MDL): Created.
3000 : *
3001 : * NOTES
3002 : *****************************************************************************
3003 : */
3004 0 : static void doSplitRight (sInt4 *Data, int numData, TDLGroupType * G,
3005 : TDLGroupType * G1, TDLGroupType * G2,
3006 : char f_primMiss, sInt4 li_primMiss,
3007 : char f_secMiss, sInt4 li_secMiss)
3008 : {
3009 : int range; /* The range to make the right most group. */
3010 : int final; /* One more than the last index in the group. */
3011 : int split; /* Where to split the group. */
3012 :
3013 : /* *INDENT-OFF* */
3014 : /* The (pow (2, ..) -1) is because 2^n - 1 is max range of a group.
3015 : * Example n = 1, 2^1 -1 = 1 range of (1,0) is 1.
3016 : * Example n = 2, 2^2 -1 = 3 range of (3,2,1,0) is 3.
3017 : * The G.bit - 1 is because we are trying to reduce the range. */
3018 : /* *INDENT-ON* */
3019 0 : range = (int) (pow (2.0, G->bit - 1) - 1) - (f_secMiss + f_primMiss);
3020 0 : final = G->start + G->num;
3021 0 : split = final;
3022 0 : myAssert (final <= numData);
3023 :
3024 0 : if (f_secMiss) {
3025 : findGroupRev2 (Data, G->start, final, li_primMiss, li_secMiss, range,
3026 0 : &split, &(G1->min), &(G1->max));
3027 : findMaxMin2 (Data, G->start, split, li_primMiss, li_secMiss,
3028 0 : &(G2->min), &(G2->max));
3029 0 : } else if (f_primMiss) {
3030 : findGroupRev1 (Data, G->start, final, li_primMiss, range, &split,
3031 0 : &(G1->min), &(G1->max));
3032 : findMaxMin1 (Data, G->start, split, li_primMiss, &(G2->min),
3033 0 : &(G2->max));
3034 : } else {
3035 : findGroupRev0 (Data, G->start, final, range, &split, &(G1->min),
3036 0 : &(G1->max));
3037 0 : findMaxMin0 (Data, G->start, split, &(G2->min), &(G2->max));
3038 : }
3039 :
3040 : G1->bit = (char) power ((uInt4) (G1->max - G1->min),
3041 0 : f_secMiss + f_primMiss);
3042 : G2->bit = (char) power ((uInt4) (G2->max - G2->min),
3043 0 : f_secMiss + f_primMiss);
3044 0 : G1->start = split;
3045 0 : G2->start = G->start;
3046 0 : G1->num = final - split;
3047 0 : G2->num = split - G->start;
3048 0 : G1->f_trySplit = 1;
3049 0 : G1->f_tryShift = 1;
3050 0 : G2->f_trySplit = 1;
3051 0 : G2->f_tryShift = 1;
3052 0 : }
3053 :
3054 : /*****************************************************************************
3055 : * ComputeGroupSize() --
3056 : *
3057 : * Arthur Taylor / MDL
3058 : *
3059 : * PURPOSE
3060 : * Compute the number of bits needed for the various elements of the groups
3061 : * as well as the total number of bits needed.
3062 : *
3063 : * ARGUMENTS
3064 : * group = Groups (Input)
3065 : * numGroup = Number of groups (Input)
3066 : * ibit = Number of bits needed for the minimum values of each group.
3067 : * Find max absolute value of group mins. Note: all group mins
3068 : * are positive (Output)
3069 : * jbit = Number of bits needed for the number of bits for each group.
3070 : * Find max absolute value of number of bits. (Output)
3071 : * kbit = Number of bits needed for the number of values for each group.
3072 : * Find max absolute value of number of values. (Output)
3073 : *
3074 : * FILES/DATABASES: None
3075 : *
3076 : * RETURNS: sInt4
3077 : * number of bits needed by the groups
3078 : *
3079 : * HISTORY
3080 : * 12/2004 Arthur Taylor (MDL): Updated from "tdlpack.c" in "C" tdlpack code
3081 : *
3082 : * NOTES
3083 : *****************************************************************************
3084 : */
3085 0 : static sInt4 ComputeGroupSize (TDLGroupType * group, int numGroup,
3086 : size_t *ibit, size_t *jbit, size_t *kbit)
3087 : {
3088 : int i; /* loop counter. */
3089 0 : sInt4 ans = 0; /* The number of bits needed. */
3090 0 : sInt4 maxMin = 0; /* The largest min value in the groups */
3091 0 : uChar maxBit = 0; /* The largest needed bits in the groups */
3092 0 : uInt4 maxNum = 0; /* The largest number of values in the groups. */
3093 :
3094 0 : for (i = 0; i < numGroup; i++) {
3095 0 : ans += group[i].bit * group[i].num;
3096 0 : if (group[i].min > maxMin) {
3097 0 : maxMin = group[i].min;
3098 : }
3099 0 : if (group[i].bit > maxBit) {
3100 0 : maxBit = group[i].bit;
3101 : }
3102 0 : if (group[i].num > maxNum) {
3103 0 : maxNum = group[i].num;
3104 : }
3105 : }
3106 : /* This only works for pos numbers... */
3107 0 : for (i = 0; (maxMin != 0); i++) {
3108 0 : maxMin = maxMin >> 1;
3109 : }
3110 : /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3111 0 : *ibit = i;
3112 : /* This only works for pos numbers... */
3113 0 : for (i = 0; (maxBit != 0); i++) {
3114 0 : maxBit = maxBit >> 1;
3115 : }
3116 : /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3117 0 : *jbit = i;
3118 : /* This only works for pos numbers... */
3119 0 : for (i = 0; (maxNum != 0); i++) {
3120 0 : maxNum = maxNum >> 1;
3121 : }
3122 : /* Allow 0 bits for min. Assumes that decoder allows 0 bits */
3123 0 : *kbit = i;
3124 0 : ans += ((*ibit) + (*jbit) + (*kbit)) * numGroup;
3125 0 : return ans;
3126 : }
3127 :
3128 : /*****************************************************************************
3129 : * splitGroup() --
3130 : *
3131 : * Arthur Taylor / MDL
3132 : *
3133 : * PURPOSE
3134 : * Tries to reduce (split) each group by 1 bit. It does this by:
3135 : * A) reduce the "bit range", and create groups that grab as much as they
3136 : * can to the right. Then reduce those groups if they improve the score.
3137 : * B) reduce the bit range and grab the left most group only, leaving the
3138 : * rest unchanged.
3139 : * C) reduce the bit range and grab the right most group only, leaving the
3140 : * rest unchanged.
3141 : *
3142 : * ARGUMENTS
3143 : * Data = The data. (Input)
3144 : * numData = The number of elements in data. (Input)
3145 : * group = The groups using the "best minGroup. (Input)
3146 : * numGroup = Number of groups (Input)
3147 : * lclGroup = The local copy of the groups (Output)
3148 : * numLclGroup = Number of local groups (Output)
3149 : * f_primMiss = Flag if we have a primary missing value (Input)
3150 : * li_primMiss = scaled primary missing value (Input)
3151 : * f_secMiss = Flag if we have a secondary missing value (Input)
3152 : * li_secMiss = scaled secondary missing value (Input)
3153 : * xFactor = Estimate of cost (in bits) of a group. (Input)
3154 : *
3155 : * FILES/DATABASES: None
3156 : *
3157 : * RETURNS: void
3158 : *
3159 : * HISTORY
3160 : * 1/2005 Arthur Taylor (MDL): Created.
3161 : *
3162 : * NOTES
3163 : *****************************************************************************
3164 : */
3165 0 : static int splitGroup (sInt4 *Data, int numData, TDLGroupType * group,
3166 : int numGroup, TDLGroupType ** lclGroup,
3167 : int *numLclGroup, char f_primMiss,
3168 : sInt4 li_primMiss, char f_secMiss,
3169 : sInt4 li_secMiss, size_t xFactor)
3170 : {
3171 : uInt4 minBit; /* The fewest # of bits, with no subdivision. */
3172 : TDLGroupType *subGroup; /* The subgroups that we tried splitting the
3173 : * primary group into. */
3174 : int numSubGroup; /* The number of groups in subGroup. */
3175 : sInt4 A_max; /* Max value of a given group. */
3176 : sInt4 A_min; /* Min value of a given group. */
3177 : sInt4 scoreA; /* The original score for a given group */
3178 : sInt4 scoreB; /* The new score */
3179 0 : int f_adjust = 0; /* Flag if group has changed. */
3180 : int f_keep; /* Flag to keep the subgroup instead of original */
3181 : int i; /* Loop counters */
3182 : int sub; /* loop counter over the sub group. */
3183 : int lclIndex; /* Used to help copy data from subGroup to answer. */
3184 : int GroupLen; /* Actual alloc'ed group len. */
3185 : int extra; /* Used to reduce number of allocs. */
3186 :
3187 : /* Figure out how few bits a group can have without being able to further
3188 : * divide it. */
3189 0 : if (f_secMiss) {
3190 : /* 11 = primMiss 10 = secMiss 01, 00 = data. */
3191 0 : minBit = 2;
3192 0 : } else if (f_primMiss) {
3193 : /* 1 = primMiss 0 = data. */
3194 : /* might try minBit = 1 here. */
3195 0 : minBit = 1;
3196 : } else {
3197 : /* 1, 0 = data. */
3198 0 : minBit = 1;
3199 : }
3200 :
3201 0 : *numLclGroup = 0;
3202 0 : *lclGroup = (TDLGroupType *) malloc (numGroup * sizeof (TDLGroupType));
3203 0 : GroupLen = numGroup;
3204 0 : extra = 0;
3205 0 : for (i = 0; i < numGroup; i++) {
3206 : /* Check if we have already tried to split this group, or it has too
3207 : * few members, or it doesn't have enough bits to split. If so, skip
3208 : * this group. */
3209 0 : if ((group[i].f_trySplit) && (group[i].num > xFactor) &&
3210 0 : (group[i].bit > minBit)) {
3211 0 : f_keep = 0;
3212 : doSplit (Data, numData, &(group[i]), &subGroup, &numSubGroup,
3213 0 : f_primMiss, li_primMiss, f_secMiss, li_secMiss, xFactor);
3214 0 : if (numSubGroup != 1) {
3215 0 : scoreA = group[i].bit * group[i].num + xFactor;
3216 0 : scoreB = 0;
3217 0 : for (sub = 0; sub < numSubGroup; sub++) {
3218 0 : scoreB += subGroup[sub].bit * subGroup[sub].num + xFactor;
3219 : }
3220 0 : if (scoreB < scoreA) {
3221 0 : f_keep = 1;
3222 0 : } else if (numSubGroup > 2) {
3223 : /* We can do "doSplitLeft" (which is breaking it into 2 groups
3224 : * the first having range n - 1, the second having range n,
3225 : * using what we know from doSplit. */
3226 0 : subGroup[1].num = group[i].num - subGroup[0].num;
3227 0 : if (f_secMiss) {
3228 0 : findMaxMin2 (Data, subGroup[1].start,
3229 0 : subGroup[1].start + subGroup[1].num,
3230 0 : li_primMiss, li_secMiss, &A_min, &A_max);
3231 0 : } else if (f_primMiss) {
3232 0 : findMaxMin1 (Data, subGroup[1].start,
3233 0 : subGroup[1].start + subGroup[1].num,
3234 0 : li_primMiss, &A_min, &A_max);
3235 : } else {
3236 0 : findMaxMin0 (Data, subGroup[1].start,
3237 0 : subGroup[1].start + subGroup[1].num,
3238 0 : &A_min, &A_max);
3239 : }
3240 0 : subGroup[1].min = A_min;
3241 0 : subGroup[1].max = A_max;
3242 0 : subGroup[1].bit = power ((uInt4) (A_max - A_min),
3243 0 : f_secMiss + f_primMiss);
3244 0 : subGroup[1].f_trySplit = 1;
3245 0 : subGroup[1].f_tryShift = 1;
3246 0 : numSubGroup = 2;
3247 0 : scoreB = subGroup[0].bit * subGroup[0].num + xFactor;
3248 0 : scoreB += subGroup[1].bit * subGroup[1].num + xFactor;
3249 0 : if (scoreB < scoreA) {
3250 0 : f_keep = 1;
3251 : }
3252 : }
3253 : }
3254 0 : if (!f_keep) {
3255 0 : if (numSubGroup == 1) {
3256 : subGroup = (TDLGroupType *) realloc (subGroup,
3257 : 2 *
3258 0 : sizeof (TDLGroupType));
3259 : }
3260 0 : numSubGroup = 2;
3261 : doSplitRight (Data, numData, &(group[i]), &(subGroup[1]),
3262 : &(subGroup[0]), f_primMiss, li_primMiss, f_secMiss,
3263 0 : li_secMiss);
3264 0 : scoreA = group[i].bit * group[i].num + xFactor;
3265 0 : scoreB = subGroup[0].bit * subGroup[0].num + xFactor;
3266 0 : scoreB += subGroup[1].bit * subGroup[1].num + xFactor;
3267 0 : if (scoreB < scoreA) {
3268 0 : f_keep = 1;
3269 : }
3270 : }
3271 0 : if (f_keep) {
3272 0 : lclIndex = *numLclGroup;
3273 0 : *numLclGroup = *numLclGroup + numSubGroup;
3274 0 : if (*numLclGroup > GroupLen) {
3275 0 : GroupLen += extra;
3276 0 : extra = 0;
3277 0 : if (*numLclGroup > GroupLen) {
3278 0 : GroupLen = *numLclGroup;
3279 : }
3280 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3281 : GroupLen *
3282 0 : sizeof (TDLGroupType));
3283 : } else {
3284 0 : extra += numSubGroup - 1;
3285 : }
3286 : memcpy ((*lclGroup) + lclIndex, subGroup,
3287 0 : numSubGroup * sizeof (TDLGroupType));
3288 0 : f_adjust = 1;
3289 : } else {
3290 0 : *numLclGroup = *numLclGroup + 1;
3291 0 : if (*numLclGroup > GroupLen) {
3292 0 : GroupLen += extra;
3293 0 : extra = 0;
3294 0 : if (*numLclGroup > GroupLen) {
3295 0 : GroupLen = *numLclGroup;
3296 : }
3297 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3298 : GroupLen *
3299 0 : sizeof (TDLGroupType));
3300 : }
3301 0 : (*lclGroup)[*numLclGroup - 1] = group[i];
3302 0 : (*lclGroup)[*numLclGroup - 1].f_trySplit = 0;
3303 : }
3304 0 : free (subGroup);
3305 0 : subGroup = NULL;
3306 : } else {
3307 0 : *numLclGroup = *numLclGroup + 1;
3308 0 : if (*numLclGroup > GroupLen) {
3309 0 : GroupLen += extra;
3310 0 : extra = 0;
3311 0 : if (*numLclGroup > GroupLen) {
3312 0 : GroupLen = *numLclGroup;
3313 : }
3314 : *lclGroup = (TDLGroupType *) realloc ((void *) *lclGroup,
3315 : GroupLen *
3316 0 : sizeof (TDLGroupType));
3317 : }
3318 0 : (*lclGroup)[*numLclGroup - 1] = group[i];
3319 0 : (*lclGroup)[*numLclGroup - 1].f_trySplit = 0;
3320 : }
3321 : }
3322 0 : myAssert (GroupLen == *numLclGroup);
3323 0 : return f_adjust;
3324 : }
3325 :
3326 : /*****************************************************************************
3327 : * shiftGroup() --
3328 : *
3329 : * Arthur Taylor / MDL
3330 : *
3331 : * PURPOSE
3332 : * Tries to shift / join the groups together. It does this by first
3333 : * calculating if a group should still exist. If it should, then it has
3334 : * each group grab as much as it can to the left without increasing its "bit
3335 : * range".
3336 : *
3337 : * ARGUMENTS
3338 : * Data = The data. (Input)
3339 : * numData = The number of elements in data. (Input)
3340 : * group = The resulting groups. (Output)
3341 : * numGroup = Number of groups (Output)
3342 : * f_primMiss = Flag if we have a primary missing value (Input)
3343 : * li_primMiss = scaled primary missing value (Input)
3344 : * f_secMiss = Flag if we have a secondary missing value (Input)
3345 : * li_secMiss = scaled secondary missing value (Input)
3346 : * xFactor = Estimate of cost (in bits) of a group. (Input)
3347 : *
3348 : * FILES/DATABASES: None
3349 : *
3350 : * RETURNS: void
3351 : *
3352 : * HISTORY
3353 : * 1/2005 Arthur Taylor (MDL): Created.
3354 : *
3355 : * NOTES
3356 : *****************************************************************************
3357 : */
3358 0 : static void shiftGroup (sInt4 *Data, int numData, TDLGroupType ** Group,
3359 : size_t *NumGroup, char f_primMiss, sInt4 li_primMiss,
3360 : char f_secMiss, sInt4 li_secMiss, int xFactor)
3361 : {
3362 0 : TDLGroupType *group = (*Group); /* Local pointer to Group. */
3363 0 : int numGroup = (*NumGroup); /* # elements in group. */
3364 : int i, j; /* loop counters. */
3365 : sInt4 A_max; /* Max value of a given group. */
3366 : sInt4 A_min; /* Min value of a given group. */
3367 : size_t begin; /* New start to the group. */
3368 : sInt4 scoreA; /* The original score for group i and i - 1 */
3369 : sInt4 scoreB; /* The new score for group i and i - 1 */
3370 : TDLGroupType G1; /* The "new" group[i - 1]. */
3371 : TDLGroupType G2; /* The "new" group[i]. */
3372 0 : int evaporate = 0; /* number of groups that "evaporated" */
3373 : int index; /* index while getting rid of "evaporated groups". */
3374 :
3375 0 : for (i = numGroup - 1; i > 0; i--) {
3376 0 : myAssert (group[i].num > 0);
3377 : /* See if we can evaporate the group n - 1 */
3378 0 : G1.start = group[i - 1].start;
3379 0 : G1.num = group[i].num + group[i - 1].num;
3380 0 : G1.min = (group[i].min < group[i - 1].min) ?
3381 0 : group[i].min : group[i - 1].min;
3382 0 : G1.max = (group[i].max > group[i - 1].max) ?
3383 0 : group[i].max : group[i - 1].max;
3384 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
3385 0 : f_secMiss + f_primMiss);
3386 0 : G1.f_trySplit = 1;
3387 0 : G1.f_tryShift = 1;
3388 0 : scoreA = group[i].bit * group[i].num + xFactor;
3389 0 : scoreA += group[i - 1].bit * group[i - 1].num + xFactor;
3390 0 : scoreB = G1.bit * G1.num + xFactor;
3391 0 : if (scoreB < scoreA) {
3392 : /* One of the groups evaporated. */
3393 0 : evaporate++;
3394 0 : group[i - 1] = G1;
3395 0 : group[i].num = 0;
3396 : /* See if that affects any of the previous groups. */
3397 0 : for (j = i + 1; j < numGroup; j++) {
3398 0 : if (group[j].num != 0) {
3399 0 : G1.start = G1.start;
3400 0 : G1.num = group[i - 1].num + group[j].num;
3401 0 : G1.min = (group[i - 1].min < group[j].min) ?
3402 0 : group[i - 1].min : group[j].min;
3403 0 : G1.max = (group[i - 1].max > group[j].max) ?
3404 0 : group[i - 1].max : group[j].max;
3405 : G1.bit = (char) power ((uInt4) (G1.max - G1.min),
3406 0 : f_secMiss + f_primMiss);
3407 0 : G1.f_trySplit = 1;
3408 0 : G1.f_tryShift = 1;
3409 0 : scoreA = group[i - 1].bit * group[i - 1].num + xFactor;
3410 0 : scoreA += group[j].bit * group[j].num + xFactor;
3411 0 : scoreB = G1.bit * G1.num + xFactor;
3412 0 : if (scoreB < scoreA) {
3413 0 : evaporate++;
3414 0 : group[i - 1] = G1;
3415 0 : group[j].num = 0;
3416 : } else {
3417 0 : break;
3418 : }
3419 : }
3420 : }
3421 0 : } else if (group[i].f_tryShift) {
3422 : /* Group did not evaporate, so do the "grabby" algorithm. */
3423 0 : if ((group[i].bit != 0) && (group[i - 1].bit >= group[i].bit)) {
3424 0 : if (f_secMiss) {
3425 0 : A_max = group[i].max;
3426 0 : A_min = group[i].min;
3427 0 : shiftGroup2 (Data, group[i].start - 1, group[i - 1].start,
3428 0 : li_primMiss, li_secMiss, group[i].bit, &A_min,
3429 0 : &A_max, &begin);
3430 0 : } else if (f_primMiss) {
3431 0 : A_max = group[i].max;
3432 0 : A_min = group[i].min;
3433 0 : shiftGroup1 (Data, group[i].start - 1, group[i - 1].start,
3434 0 : li_primMiss, group[i].bit, &A_min, &A_max,
3435 0 : &begin);
3436 : } else {
3437 0 : A_max = group[i].max;
3438 0 : A_min = group[i].min;
3439 0 : shiftGroup0 (Data, group[i].start - 1, group[i - 1].start,
3440 0 : group[i].bit, &A_min, &A_max, &begin);
3441 : }
3442 0 : if (begin != group[i].start) {
3443 : /* Re-Calculate min/max of group[i - 1], since it could have
3444 : * moved to group[i] */
3445 0 : G1 = group[i - 1];
3446 0 : G2 = group[i];
3447 0 : G2.min = A_min;
3448 0 : G2.max = A_max;
3449 0 : G1.num -= (group[i].start - begin);
3450 0 : if (f_secMiss) {
3451 : findMaxMin2 (Data, G1.start, G1.start + G1.num,
3452 0 : li_primMiss, li_secMiss, &A_min, &A_max);
3453 0 : } else if (f_primMiss) {
3454 : findMaxMin1 (Data, G1.start, G1.start + G1.num,
3455 0 : li_primMiss, &A_min, &A_max);
3456 : } else {
3457 : findMaxMin0 (Data, G1.start, G1.start + G1.num,
3458 0 : &A_min, &A_max);
3459 : }
3460 0 : if ((A_min != G1.min) || (A_max != G1.max)) {
3461 0 : G1.min = A_min;
3462 0 : G1.max = A_max;
3463 : G1.bit = (char) power ((uInt4) (A_max - A_min),
3464 0 : f_secMiss + f_primMiss);
3465 0 : G1.f_trySplit = 1;
3466 0 : G1.f_tryShift = 1;
3467 : }
3468 0 : G2.num += (group[i].start - begin);
3469 0 : G2.start = begin;
3470 0 : G2.f_trySplit = 1;
3471 0 : G2.f_tryShift = 1;
3472 0 : scoreA = group[i].bit * group[i].num + xFactor;
3473 0 : scoreA += group[i - 1].bit * group[i - 1].num + xFactor;
3474 0 : scoreB = G2.bit * G2.num + xFactor;
3475 0 : if (G1.num != 0) {
3476 0 : scoreB += G1.bit * G1.num + xFactor;
3477 : }
3478 : #ifdef DEBUG
3479 0 : if (scoreB > scoreA) {
3480 0 : printf ("Made score worse!\n");
3481 : }
3482 : #endif
3483 0 : if (begin == group[i - 1].start) {
3484 : /* Grabby algorithm evaporated a group. */
3485 0 : evaporate++;
3486 0 : myAssert (G1.num == 0);
3487 : /* Switch the evaporating group to other side so we have
3488 : * potential to continue evaporating the next group down. */
3489 0 : group[i - 1] = G2;
3490 0 : group[i] = G1;
3491 : } else {
3492 0 : group[i - 1] = G1;
3493 0 : group[i] = G2;
3494 : }
3495 : } else {
3496 0 : group[i].f_tryShift = 0;
3497 : }
3498 : }
3499 : }
3500 : }
3501 : /* Loop through the grid removing the evaporated groups. */
3502 0 : if (evaporate != 0) {
3503 0 : index = 0;
3504 0 : for (i = 0; i < numGroup; i++) {
3505 0 : if (group[i].num != 0) {
3506 0 : group[index] = group[i];
3507 0 : index++;
3508 : }
3509 : }
3510 0 : *NumGroup = numGroup - evaporate;
3511 : *Group = (TDLGroupType *) realloc ((void *) (*Group),
3512 0 : *NumGroup * sizeof (TDLGroupType));
3513 : }
3514 0 : }
3515 :
3516 : /*****************************************************************************
3517 : * GroupIt() --
3518 : *
3519 : * Arthur Taylor / MDL
3520 : *
3521 : * PURPOSE
3522 : * Attempts to find groups for packing the data. It starts by preparing
3523 : * the data, by removing the overall min value.
3524 : *
3525 : * Next it Creates any 0 bit groups (primary missing: missings are all 0
3526 : * bit groups, const values are 0 bit groups. No missing: const values are
3527 : * 0 bit groups.)
3528 : *
3529 : * Next it tries to reduce (split) each group by 1 bit. It does this by:
3530 : * A) reduce the "bit range", and create groups that grab as much as they
3531 : * can to the right. Then reduce those groups if they improve the score.
3532 : * B) reduce the bit range and grab the left most group only, leaving the
3533 : * rest unchanged.
3534 : * C) reduce the bit range and grab the right most group only, leaving the
3535 : * rest unchanged.
3536 : *
3537 : * Next it tries to shift / join those groups together. It does this by
3538 : * first calculating if a group should still exist. If it should, then it
3539 : * has each group grab as much as it can to the left without increasing its
3540 : * "bit range".
3541 : *
3542 : * ARGUMENTS
3543 : * OverallMin = The overall min value in the data. (Input)
3544 : * Data = The data. (Input)
3545 : * numData = The number of elements in data. (Input)
3546 : * group = The resulting groups. (Output)
3547 : * numGroup = Number of groups (Output)
3548 : * f_primMiss = Flag if we have a primary missing value (Input)
3549 : * li_primMiss = scaled primary missing value (Input)
3550 : * f_secMiss = Flag if we have a secondary missing value (Input)
3551 : * li_secMiss = scaled secondary missing value (Input)
3552 : * groupSize = How many bytes the groups and data will take. (Output)
3553 : * ibit = # of bits for largest minimum value in groups (Output)
3554 : * jbit = # of bits for largest # bits in groups (Output)
3555 : * kbit = # of bits for largest # values in groups (Output)
3556 : *
3557 : * FILES/DATABASES: None
3558 : *
3559 : * RETURNS: void
3560 : *
3561 : * HISTORY
3562 : * 1/2005 Arthur Taylor (MDL): Created.
3563 : *
3564 : * NOTES
3565 : * 1) Have not implemented const 0 bit groups for prim miss or no miss.
3566 : * 2) MAX_GROUP_LEN (found experimentally) is useful to keep cost of groups
3567 : * down.
3568 : *****************************************************************************
3569 : */
3570 : #define MAX_GROUP_LEN 255
3571 0 : static void GroupIt (sInt4 OverallMin, sInt4 *Data, size_t numData,
3572 : TDLGroupType ** group, size_t *numGroup, char f_primMiss,
3573 : sInt4 li_primMiss, char f_secMiss,
3574 : sInt4 li_secMiss, sInt4 *groupSize, size_t *ibit,
3575 : size_t *jbit, size_t *kbit)
3576 : {
3577 : sInt4 A_max; /* Max value of a given group. */
3578 : sInt4 A_min; /* Min value of a given group. */
3579 : TDLGroupType G; /* Used to init the groups. */
3580 : int f_adjust; /* Flag if we have changed the groups. */
3581 : TDLGroupType *lclGroup; /* A temporary copy of the groups. */
3582 : int numLclGroup; /* # of groups in lclGroup. */
3583 : size_t xFactor; /* Estimate of cost (in bits) of a group. */
3584 : size_t i; /* loop counter. */
3585 :
3586 : /* Subtract the Overall Min Value. */
3587 0 : if (OverallMin != 0) {
3588 0 : if (f_secMiss) {
3589 0 : for (i = 0; i < numData; i++) {
3590 0 : if ((Data[i] != li_secMiss) && (Data[i] != li_primMiss)) {
3591 0 : Data[i] -= OverallMin;
3592 : /* Check if we accidently adjusted to prim or sec, if so add
3593 : * 1. */
3594 0 : if ((Data[i] == li_secMiss) || (Data[i] == li_primMiss)) {
3595 0 : myAssert (1 == 2);
3596 0 : Data[i]++;
3597 0 : if ((Data[i] == li_secMiss) || (Data[i] == li_primMiss)) {
3598 0 : myAssert (1 == 2);
3599 0 : Data[i]++;
3600 : }
3601 : }
3602 : }
3603 : }
3604 0 : } else if (f_primMiss) {
3605 0 : for (i = 0; i < numData; i++) {
3606 0 : if (Data[i] != li_primMiss) {
3607 0 : Data[i] -= OverallMin;
3608 : /* Check if we accidently adjusted to prim or sec, if so add
3609 : * 1. */
3610 0 : if (Data[i] == li_primMiss) {
3611 0 : myAssert (1 == 2);
3612 0 : Data[i]++;
3613 : }
3614 : }
3615 : }
3616 : } else {
3617 0 : for (i = 0; i < numData; i++) {
3618 0 : Data[i] -= OverallMin;
3619 : }
3620 : }
3621 : }
3622 :
3623 0 : myAssert ((f_secMiss == 0) || (f_secMiss == 1));
3624 0 : myAssert ((f_primMiss == 0) || (f_primMiss == 1));
3625 :
3626 : /* Create zero groups. */
3627 0 : *numGroup = 0;
3628 0 : *group = NULL;
3629 0 : if (f_primMiss) {
3630 0 : G.min = Data[0];
3631 0 : G.max = Data[0];
3632 0 : G.num = 1;
3633 0 : G.start = 0;
3634 0 : for (i = 1; i < numData; i++) {
3635 0 : if (G.min == li_primMiss) {
3636 0 : if (Data[i] == li_primMiss) {
3637 0 : G.num++;
3638 0 : if (G.num == (MAX_GROUP_LEN + 1)) {
3639 : /* Close a missing group */
3640 0 : G.f_trySplit = 0;
3641 0 : G.f_tryShift = 1;
3642 0 : G.bit = 0;
3643 0 : G.min = 0;
3644 0 : G.max = 0;
3645 0 : G.num = MAX_GROUP_LEN;
3646 0 : (*numGroup)++;
3647 : *group = (TDLGroupType *) realloc ((void *) *group,
3648 : *numGroup *
3649 0 : sizeof (TDLGroupType));
3650 0 : (*group)[(*numGroup) - 1] = G;
3651 : /* Init a missing group. */
3652 0 : G.min = Data[i];
3653 0 : G.max = Data[i];
3654 0 : G.num = 1;
3655 0 : G.start = i;
3656 : }
3657 : } else {
3658 : /* Close a missing group */
3659 0 : G.f_trySplit = 0;
3660 0 : G.f_tryShift = 1;
3661 0 : G.bit = 0;
3662 0 : G.min = 0;
3663 0 : G.max = 0;
3664 0 : (*numGroup)++;
3665 : *group = (TDLGroupType *) realloc ((void *) *group,
3666 : *numGroup *
3667 0 : sizeof (TDLGroupType));
3668 0 : (*group)[(*numGroup) - 1] = G;
3669 : /* Init a non-missing group. */
3670 0 : G.min = Data[i];
3671 0 : G.max = Data[i];
3672 0 : G.num = 1;
3673 0 : G.start = i;
3674 : }
3675 : } else {
3676 0 : if (Data[i] == li_primMiss) {
3677 : /* Close a non-missing group */
3678 0 : G.f_trySplit = 1;
3679 0 : G.f_tryShift = 1;
3680 : G.bit = (char) power ((uInt4) (G.max - G.min),
3681 0 : f_secMiss + f_primMiss);
3682 0 : myAssert (G.bit != 0);
3683 0 : if ((G.min == 0) && (G.bit == 0) && (f_primMiss == 1)) {
3684 : printf ("Warning: potential confusion between const value "
3685 0 : "and prim-missing.\n");
3686 0 : G.bit = 1;
3687 : }
3688 0 : (*numGroup)++;
3689 : *group = (TDLGroupType *) realloc ((void *) *group,
3690 : *numGroup *
3691 0 : sizeof (TDLGroupType));
3692 0 : (*group)[(*numGroup) - 1] = G;
3693 : /* Init a missing group. */
3694 0 : G.min = Data[i];
3695 0 : G.max = Data[i];
3696 0 : G.num = 1;
3697 0 : G.start = i;
3698 : } else {
3699 0 : if (G.min > Data[i]) {
3700 0 : G.min = Data[i];
3701 0 : } else if (G.max < Data[i]) {
3702 0 : G.max = Data[i];
3703 : }
3704 0 : G.num++;
3705 : }
3706 : }
3707 : }
3708 0 : if (G.min == li_primMiss) {
3709 : /* Close a missing group */
3710 0 : G.f_trySplit = 0;
3711 0 : G.f_tryShift = 1;
3712 0 : G.bit = 0;
3713 0 : G.min = 0;
3714 0 : G.max = 0;
3715 0 : (*numGroup)++;
3716 : *group = (TDLGroupType *) realloc ((void *) *group,
3717 : *numGroup *
3718 0 : sizeof (TDLGroupType));
3719 0 : (*group)[(*numGroup) - 1] = G;
3720 : } else {
3721 : /* Close a non-missing group */
3722 0 : G.f_trySplit = 1;
3723 0 : G.f_tryShift = 1;
3724 : G.bit = (char) power ((uInt4) (G.max - G.min),
3725 0 : f_secMiss + f_primMiss);
3726 0 : myAssert (G.bit != 0);
3727 0 : if ((G.min == 0) && (G.bit == 0) && (f_primMiss == 1)) {
3728 : printf ("Warning: potential confusion between const value and "
3729 0 : "prim-missing.\n");
3730 0 : G.bit = 1;
3731 : }
3732 0 : (*numGroup)++;
3733 : *group = (TDLGroupType *) realloc ((void *) *group,
3734 : *numGroup *
3735 0 : sizeof (TDLGroupType));
3736 0 : (*group)[(*numGroup) - 1] = G;
3737 : }
3738 : } else {
3739 : /* Already handled the f_primMiss case */
3740 0 : if (f_secMiss) {
3741 : findMaxMin2 (Data, 0, numData, li_primMiss, li_secMiss, &A_min,
3742 0 : &A_max);
3743 : } else {
3744 0 : findMaxMin0 (Data, 0, numData, &A_min, &A_max);
3745 : }
3746 0 : G.start = 0;
3747 0 : G.num = numData;
3748 0 : G.min = A_min;
3749 0 : G.max = A_max;
3750 0 : G.bit = (char) power ((uInt4) (A_max - A_min), f_secMiss + f_primMiss);
3751 0 : G.f_trySplit = 1;
3752 0 : G.f_tryShift = 1;
3753 0 : *numGroup = 1;
3754 0 : *group = (TDLGroupType *) malloc (sizeof (TDLGroupType));
3755 0 : (*group)[0] = G;
3756 : }
3757 :
3758 0 : lclGroup = NULL;
3759 0 : numLclGroup = 0;
3760 0 : *groupSize = ComputeGroupSize (*group, *numGroup, ibit, jbit, kbit);
3761 0 : xFactor = *ibit + *jbit + *kbit;
3762 : #ifdef DEBUG
3763 : /*
3764 : printf ("NumGroup = %d: Bytes = %ld: XFactor %d\n", *numGroup,
3765 : (*groupSize / 8) + 1, xFactor);
3766 : */
3767 : #endif
3768 :
3769 0 : f_adjust = 1;
3770 0 : while (f_adjust) {
3771 : f_adjust = splitGroup (Data, numData, *group, *numGroup, &lclGroup,
3772 : &numLclGroup, f_primMiss, li_primMiss,
3773 0 : f_secMiss, li_secMiss, xFactor);
3774 0 : free (*group);
3775 0 : *group = lclGroup;
3776 0 : *numGroup = numLclGroup;
3777 :
3778 0 : if (f_adjust) {
3779 : shiftGroup (Data, numData, group, numGroup, f_primMiss,
3780 0 : li_primMiss, f_secMiss, li_secMiss, xFactor);
3781 0 : *groupSize = ComputeGroupSize (*group, *numGroup, ibit, jbit, kbit);
3782 0 : if (xFactor != *ibit + *jbit + *kbit) {
3783 0 : for (i = 0; i < *numGroup; i++) {
3784 0 : if (((*group)[i].num > *ibit + *jbit + *kbit) &&
3785 0 : ((*group)[i].num <= xFactor)) {
3786 0 : (*group)[i].f_trySplit = 1;
3787 : }
3788 : }
3789 : }
3790 0 : xFactor = *ibit + *jbit + *kbit;
3791 : #ifdef DEBUG
3792 : /*
3793 : printf ("NumGroup = %d: Bytes = %ld: XFactor %d\n", *numGroup,
3794 : (*groupSize / 8) + 1, xFactor);
3795 : fflush (stdout);
3796 : */
3797 : #endif
3798 : }
3799 : }
3800 0 : }
3801 :
3802 : /*****************************************************************************
3803 : * GroupPack() --
3804 : *
3805 : * Arthur Taylor / MDL
3806 : *
3807 : * PURPOSE
3808 : * To compute groups for packing the data using complex or second order
3809 : * complex packing.
3810 : *
3811 : * ARGUMENTS
3812 : * Src = The original data. (Input)
3813 : * Dst = The scaled data. (Output)
3814 : * numData = The number of elements in data. (Input)
3815 : * DSF = Decimal Scale Factor for scaling the data. (Input)
3816 : * BSF = Binary Scale Factor for scaling the data. (Input)
3817 : * f_primMiss = Flag saying if we have a primary missing value (In/Out)
3818 : * primMiss = primary missing value. (In/Out)
3819 : * f_secMiss = Flag saying if we have a secondary missing value (In/Out)
3820 : * secMiss = secondary missing value. (In/Out)
3821 : * f_grid = Flag if this is grid data (or vector) (Input)
3822 : * NX = The number of X values. (Input)
3823 : * NY = The number of Y values. (Input)
3824 : * f_sndOrder = Flag if we should do second order diffencing (Output)
3825 : * group = Resulting groups. (Output)
3826 : * numGroup = Number of groups. (Output)
3827 : * Min = Overall minimum. (Output)
3828 : * a1 = if f_sndOrder, the first first order difference (Output)
3829 : * b2 = if f_sndOrder, the first second order difference (Output)
3830 : * groupSize = How many bytes the groups and data will take. (Output)
3831 : * ibit = # of bits for largest minimum value in groups (Output)
3832 : * jbit = # of bits for largest # bits in groups (Output)
3833 : * kbit = # of bits for largest # values in groups (Output)
3834 : *
3835 : * FILES/DATABASES: None
3836 : *
3837 : * RETURNS: int (could use errSprintf())
3838 : * 0 = OK
3839 : * -1 = Primary or Secondary missing value == 0.
3840 : *
3841 : * HISTORY
3842 : * 12/2004 Arthur Taylor (MDL): Updated from "group.c" in "C" tdlpack code.
3843 : * 1/2005 AAT: Cleaned up.
3844 : *
3845 : * NOTES
3846 : *****************************************************************************
3847 : */
3848 0 : static int GroupPack (double *Src, sInt4 **Dst, sInt4 numData,
3849 : int DSF, int BSF, char *f_primMiss, double *primMiss,
3850 : char *f_secMiss, double *secMiss, char f_grid,
3851 : short int NX, short int NY, char *f_sndOrder,
3852 : TDLGroupType ** group, size_t *numGroup,
3853 : sInt4 *Min, sInt4 *a1, sInt4 *b2,
3854 : sInt4 *groupSize, size_t *ibit, size_t *jbit,
3855 : size_t *kbit)
3856 : {
3857 0 : sInt4 *SecDiff = NULL; /* Consists of the 2nd order differences if *
3858 : * requested. */
3859 : sInt4 *Data; /* The scaled data. */
3860 : char f_min; /* Flag saying overallMin is valid. */
3861 : sInt4 overallMin; /* The overall min of the scaled data. */
3862 : sInt4 secMin; /* The overall min of the 2nd order differences */
3863 0 : sInt4 li_primMiss = 0; /* The scaled primary missing value */
3864 0 : sInt4 li_secMiss = 0; /* The scaled secondary missing value */
3865 0 : int minGroup = 20; /* The minimum group size. Equivalent to xFactor?
3866 : * Chose 20 because that was a good estimate of
3867 : * XFactor. */
3868 :
3869 : /* Check consistency of f_primMiss and f_secMiss. */
3870 0 : if (*primMiss == *secMiss) {
3871 0 : *f_secMiss = 0;
3872 : }
3873 0 : if ((*f_secMiss) && (!(*f_primMiss))) {
3874 0 : *f_primMiss = *f_secMiss;
3875 0 : *primMiss = *secMiss;
3876 0 : *f_secMiss = 0;
3877 : }
3878 0 : if (*f_secMiss && (*secMiss == 0)) {
3879 0 : errSprintf ("Error: Secondary missing value not allowed to = 0.\n");
3880 0 : return -1;
3881 : }
3882 0 : if (*f_primMiss && (*primMiss == 0)) {
3883 0 : errSprintf ("Error: Primary missing value not allowed to = 0.\n");
3884 0 : return -1;
3885 : }
3886 :
3887 : /* Check minGroup size. */
3888 0 : if (minGroup > numData) {
3889 0 : minGroup = numData;
3890 : }
3891 :
3892 : /* Scale the data and check if we can change f_prim or f_sec. */
3893 : /* Note: if we use sec_diff, we have a different overall min. */
3894 0 : f_min = 0;
3895 0 : Data = (sInt4 *) malloc (numData * sizeof (sInt4));
3896 : TDL_ScaleData (Src, Data, numData, DSF, BSF, f_primMiss, primMiss,
3897 0 : f_secMiss, secMiss, &f_min, &overallMin);
3898 : /* Note: ScaleData also scales missing values. */
3899 0 : if (*f_primMiss) {
3900 0 : li_primMiss = (sInt4) (*primMiss * SCALE_MISSING + .5);
3901 : }
3902 0 : if (*f_secMiss) {
3903 0 : li_secMiss = (sInt4) (*secMiss * SCALE_MISSING + .5);
3904 : }
3905 :
3906 : /* Reason this is after TDL_ScaleData is we don't want to reoder the
3907 : * caller's copy of the data. */
3908 0 : if (f_grid) {
3909 0 : TDL_ReorderGrid (Data, NX, NY);
3910 : } else {
3911 : /* TDLPack requires the following (see pack2d.f and pack1d.f) */
3912 0 : *f_sndOrder = 0;
3913 : }
3914 : /* TDLPack requires the following (see pack2d.f) */
3915 0 : if (*f_secMiss) {
3916 0 : *f_sndOrder = 0;
3917 : }
3918 : /* If overallMin is "invalid" then they are all prim or sec values. */
3919 : /* See pack2d.f line 336 */
3920 : /* IF ALL VALUES ARE MISSING, JUST PACK; DON'T CONSIDER 2ND ORDER
3921 : * DIFFERENCES. */
3922 0 : if (!f_min) {
3923 0 : *f_sndOrder = 0;
3924 : }
3925 :
3926 : /* This has to be after TDL_ReorderGrid */
3927 0 : if (*f_sndOrder) {
3928 0 : SecDiff = (sInt4 *) malloc (numData * sizeof (sInt4));
3929 0 : if (TDL_GetSecDiff (Data, numData, SecDiff, *f_primMiss, li_primMiss,
3930 : a1, b2, &secMin)) {
3931 : /* Problem finding SecDiff, so we don't bother with it. */
3932 0 : *f_sndOrder = 0;
3933 : } else {
3934 : /* Check if it is worth doing second order differences. */
3935 0 : if (*f_primMiss) {
3936 : *f_sndOrder = TDL_UseSecDiff_Prim (Data, numData, SecDiff,
3937 0 : li_primMiss, minGroup);
3938 : } else {
3939 0 : *f_sndOrder = TDL_UseSecDiff (Data, numData, SecDiff, minGroup);
3940 : }
3941 : }
3942 : }
3943 :
3944 : /* Side affect of GroupIt2: it subtracts OverallMin from Data. */
3945 0 : if (!(*f_sndOrder)) {
3946 : GroupIt (overallMin, Data, numData, group, numGroup, *f_primMiss,
3947 : li_primMiss, *f_secMiss, li_secMiss, groupSize, ibit, jbit,
3948 0 : kbit);
3949 0 : *Min = overallMin;
3950 0 : *a1 = 0;
3951 0 : *b2 = 0;
3952 0 : *Dst = Data;
3953 0 : free (SecDiff);
3954 : } else {
3955 : GroupIt (secMin, SecDiff, numData, group, numGroup, *f_primMiss,
3956 : li_primMiss, *f_secMiss, li_secMiss, groupSize, ibit, jbit,
3957 0 : kbit);
3958 0 : *Min = secMin;
3959 0 : *Dst = SecDiff;
3960 0 : free (Data);
3961 : }
3962 0 : return 0;
3963 : }
3964 :
3965 : /*****************************************************************************
3966 : * WriteTDLPRecord() --
3967 : *
3968 : * Arthur Taylor / MDL
3969 : *
3970 : * PURPOSE
3971 : * Writes a TDLP message to file.
3972 : *
3973 : * ARGUMENTS
3974 : * fp = An opened TDLP file already at the correct location. (Input)
3975 : * Data = The data to write. (Input)
3976 : * DataLen = Length of Data. (Input)
3977 : * DSF = Decimal scale factor to apply to the data (Input)
3978 : * BSF = Binary scale factor to apply to the data (Input)
3979 : * f_primMiss = Flag saying if we have a primary missing value (Input)
3980 : * primMiss = primary missing value. (Input)
3981 : * f_secMiss = Flag saying if we have a secondary missing value (Input)
3982 : * secMiss = secondary missing value. (Input)
3983 : * gds = The grid definition section (Input)
3984 : * comment = Describes the kind of data (max 32 bytes). (Input)
3985 : * refTime = The reference (creation) time of this message. (Input)
3986 : * ID1 = TDLPack ID1 (Input)
3987 : * ID2 = TDLPack ID2 (Input)
3988 : * ID3 = TDLPack ID3 (Input)
3989 : * ID4 = TDLPack ID4 (Input)
3990 : * projSec = The projection in seconds (Input)
3991 : * processNum = The process number that created it (Input)
3992 : * seqNum = The sequence number that created it (Input)
3993 : *
3994 : * FILES/DATABASES:
3995 : * An already opened file pointing to the desired TDLP message.
3996 : *
3997 : * RETURNS: int (could use errSprintf())
3998 : * 0 = OK
3999 : * -1 = comment is too long.
4000 : * -2 = projHr is inconsistent with ID3.
4001 : * -3 = Type of map projection that TDLP can't handle.
4002 : * -4 = Primary or Secondary missing value == 0.
4003 : *
4004 : * HISTORY
4005 : * 12/2004 Arthur Taylor (MDL): Created
4006 : * 1/2005 AAT: Cleaned up.
4007 : *
4008 : * NOTES
4009 : *****************************************************************************
4010 : */
4011 0 : int WriteTDLPRecord (FILE * fp, double *Data, sInt4 DataLen, int DSF,
4012 : int BSF, char f_primMiss, double primMiss,
4013 : char f_secMiss, double secMiss, gdsType *gds,
4014 : char *comment, double refTime, sInt4 ID1,
4015 : sInt4 ID2, sInt4 ID3, sInt4 ID4,
4016 : sInt4 projSec, sInt4 processNum, sInt4 seqNum)
4017 : {
4018 : sInt4 *Scaled; /* The scaled data. */
4019 : TDLGroupType *group; /* The groups used to pack the data. */
4020 : size_t numGroup; /* Number of groups. */
4021 0 : char f_grid = 1; /* Flag if this is gridded data. In theory can handle
4022 : * vector data, but haven't tested it. */
4023 : char f_sndOrder; /* Flag if we should try second order packing. */
4024 : sInt4 overallMin; /* Overall min value of the scaled data. */
4025 : sInt4 a1; /* 2nd order difference: 1st value. */
4026 : sInt4 b2; /* 2nd order difference: 1st 1st order difference */
4027 : sInt4 li_primMiss; /* Scaled primary missing value. */
4028 : sInt4 li_secMiss; /* Scaled secondary missing value. */
4029 : int mbit; /* # of bits for b2. */
4030 : int nbit; /* # of bits for overallMin. */
4031 : size_t ibit; /* # of bits for largest minimum value in groups */
4032 : size_t jbit; /* # of bits for largest # bits in groups */
4033 : size_t kbit; /* # of bits for largest # values in groups */
4034 : int sec1Len; /* Length of section 1. */
4035 : size_t pad; /* Number of bytes to pad the message to get to the
4036 : * correct byte boundary. */
4037 : sInt4 groupSize; /* How many bytes the groups and data will take. */
4038 : sInt4 sec4Len; /* Length of section 4. */
4039 : sInt4 tdlRecSize; /* Size of the TDLP message. */
4040 : sInt4 recSize; /* Actual record size (including FORTRAN bytes). */
4041 : int commentLen; /* Length of comment */
4042 : short int projHr; /* The hours part of the forecast projection. */
4043 : char projMin; /* The minutes part of the forecast projection. */
4044 : sInt4 year; /* The reference year. */
4045 : int month, day; /* The reference month day. */
4046 : int hour, min; /* The reference hour minute. */
4047 : double sec; /* The reference second. */
4048 0 : char f_bitmap = 0; /* Bitmap flag: not implemented in specs. */
4049 0 : char f_simple = 0; /* Simple Pack flag: not implemented in specs. */
4050 : int gridType; /* Which type of grid. (Polar, Mercator, Lambert). */
4051 : int dataCnt; /* Keeps track of which element we are writing. */
4052 : sInt4 max0; /* The max value in a group. Represents primary or *
4053 : * secondary missing value depending on scheme. */
4054 : sInt4 max1; /* The next to max value in a group. Represents *
4055 : * secondary missing value. */
4056 : size_t i, j; /* loop counters */
4057 : sInt4 li_temp; /* Temporary variable (sInt4). */
4058 : short int si_temp; /* Temporary variable (short int). */
4059 : double d_temp; /* Temporary variable (double). */
4060 : char buffer[6]; /* Used to write reserved values */
4061 : uChar pbuf; /* A buffer of bits that weren't written to disk */
4062 : sChar pbufLoc; /* Where in pbuf to add more bits. */
4063 :
4064 0 : commentLen = strlen (comment);
4065 0 : if (commentLen > 32) {
4066 0 : errSprintf ("Error: '%s' is > 32 bytes long\n", comment);
4067 0 : return -1;
4068 : }
4069 0 : projHr = projSec / 3600;
4070 0 : projMin = (projSec % 3600) / 60;
4071 0 : if (projHr != (ID3 - ((ID3 / 1000) * 1000))) {
4072 : errSprintf ("Error: projHr = %d is inconsistent with ID3 = %ld\n",
4073 0 : projHr, ID3);
4074 0 : return -2;
4075 : }
4076 0 : if (f_grid) {
4077 0 : switch (gds->projType) {
4078 : case GS3_POLAR:
4079 0 : gridType = TDLP_POLAR;
4080 0 : break;
4081 : case GS3_LAMBERT:
4082 0 : gridType = TDLP_LAMBERT;
4083 0 : break;
4084 : case GS3_MERCATOR:
4085 0 : gridType = TDLP_MERCATOR;
4086 0 : break;
4087 : default:
4088 : errSprintf ("TDLPack can't handle GRIB projection type %d\n",
4089 0 : gds->projType);
4090 0 : return -3;
4091 : }
4092 : }
4093 :
4094 0 : if (GroupPack (Data, &Scaled, DataLen, DSF, BSF, &f_primMiss, &primMiss,
4095 : &f_secMiss, &secMiss, f_grid, gds->Nx, gds->Ny,
4096 : &f_sndOrder, &group, &numGroup, &overallMin, &a1, &b2,
4097 : &groupSize, &ibit, &jbit, &kbit) != 0) {
4098 0 : return -4;
4099 : }
4100 :
4101 : /* Make sure missing data is properly scaled. */
4102 0 : if (f_primMiss) {
4103 0 : li_primMiss = (sInt4) (primMiss * SCALE_MISSING + .5);
4104 : }
4105 0 : if (f_secMiss) {
4106 0 : li_secMiss = (sInt4) (secMiss * SCALE_MISSING + .5);
4107 : }
4108 :
4109 : /* Compute TDL record size. */
4110 : /* *INDENT-OFF* */
4111 : /* TDL Record size
4112 : * 8 (section 0),
4113 : * 39 + strlen(comment) (section 1),
4114 : * 0 or 28 (depending on if you have a section 2),
4115 : * 0 (section 3),
4116 : * 16 + 5 + 1 + nbit(min val) + 16 + 5 + 5 + 5 + GroupSize() (section 4)
4117 : * 4 (section 5)
4118 : * pad (to even 8 bytes...)
4119 : * Group size uses:
4120 : * ibit * num_groups + jbit * num_groups +
4121 : * kbit * num_groups + group.nbit * group.nvalues */
4122 : /* *INDENT-ON* */
4123 0 : sec1Len = 39 + commentLen;
4124 0 : if (overallMin < 0) {
4125 0 : nbit = power (-1 * overallMin, 0);
4126 : } else {
4127 0 : nbit = power (overallMin, 0);
4128 : }
4129 0 : if (!f_sndOrder) {
4130 0 : if (f_secMiss) {
4131 : sec4Len = 16 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4132 0 : + groupSize) / 8.));
4133 0 : } else if (f_primMiss) {
4134 : sec4Len = 12 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4135 0 : + groupSize) / 8.));
4136 : } else {
4137 : sec4Len = 8 + (sInt4) (ceil ((5 + 1 + nbit + 16 + 5 + 5 + 5
4138 0 : + groupSize) / 8.));
4139 : }
4140 : } else {
4141 0 : if (b2 < 0) {
4142 0 : mbit = power (-1 * b2, 0);
4143 : } else {
4144 0 : mbit = power (b2, 0);
4145 : }
4146 0 : if (f_secMiss) {
4147 : sec4Len =
4148 : 16 +
4149 : (sInt4) (ceil
4150 : ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4151 0 : groupSize) / 8.));
4152 0 : } else if (f_primMiss) {
4153 : sec4Len =
4154 : 12 +
4155 : (sInt4) (ceil
4156 : ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4157 0 : groupSize) / 8.));
4158 : } else {
4159 : sec4Len =
4160 : 8 +
4161 : (sInt4) (ceil
4162 : ((32 + 5 + 1 + mbit + 5 + 1 + nbit + 16 + 5 + 5 + 5 +
4163 0 : groupSize) / 8.));
4164 : }
4165 : }
4166 0 : if (f_grid) {
4167 0 : tdlRecSize = 8 + sec1Len + 28 + 0 + sec4Len + 4;
4168 : } else {
4169 0 : tdlRecSize = 8 + sec1Len + 0 + 0 + sec4Len + 4;
4170 : }
4171 : /* Actual recSize is 8 + round(tdlRecSize) to nearest 8 byte boundary. */
4172 0 : recSize = 8 + (sInt4) (ceil (tdlRecSize / 8.0)) * 8;
4173 0 : pad = (int) ((recSize - 8) - tdlRecSize);
4174 :
4175 : /* --- Start Writing record. --- */
4176 : /* First write FORTRAN record information */
4177 0 : fread(&(recSize), sizeof (sInt4), 1, fp);
4178 0 : li_temp = 0;
4179 0 : FWRITE_BIG (&(li_temp), sizeof (sInt4), 1, fp);
4180 0 : li_temp = recSize - 8; /* FORTRAN rec length. */
4181 0 : FWRITE_BIG (&(li_temp), sizeof (sInt4), 1, fp);
4182 :
4183 : /* Now write section 0. */
4184 0 : fwrite ("TDLP", sizeof (char), 4, fp);
4185 0 : FWRITE_ODDINT_BIG (&(tdlRecSize), 3, fp);
4186 : /* version number */
4187 0 : fputc (0, fp);
4188 :
4189 : /* Now write section 1. */
4190 0 : fputc (sec1Len, fp);
4191 : /* output type specification... */
4192 0 : i = 0;
4193 0 : if (f_grid) {
4194 0 : i |= 1;
4195 : }
4196 0 : if (f_bitmap) {
4197 0 : i |= 2;
4198 : }
4199 0 : fputc (i, fp);
4200 : /* tempTime = gmtime (&(refTime));*/
4201 0 : Clock_PrintDate (refTime, &year, &month, &day, &hour, &min, &sec);
4202 : /* year = tempTime->tm_year + 1900;
4203 : month = tempTime->tm_mon + 1;
4204 : day = tempTime->tm_mday;
4205 : hour = tempTime->tm_hour;
4206 : min = tempTime->tm_min;
4207 : */
4208 0 : si_temp = year;
4209 0 : FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4210 0 : fputc (month, fp);
4211 0 : fputc (day, fp);
4212 0 : fputc (hour, fp);
4213 0 : fputc (min, fp);
4214 0 : li_temp = (year * 1000000L + month * 10000L + day * 100 + hour);
4215 0 : FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4216 0 : FWRITE_BIG (&ID1, sizeof (sInt4), 1, fp);
4217 0 : FWRITE_BIG (&ID2, sizeof (sInt4), 1, fp);
4218 0 : FWRITE_BIG (&ID3, sizeof (sInt4), 1, fp);
4219 0 : FWRITE_BIG (&ID4, sizeof (sInt4), 1, fp);
4220 0 : FWRITE_BIG (&projHr, sizeof (short int), 1, fp);
4221 0 : fputc (projMin, fp);
4222 0 : fputc (processNum, fp);
4223 0 : fputc (seqNum, fp);
4224 0 : i = (DSF < 0) ? 128 - DSF : DSF;
4225 0 : fputc (i, fp);
4226 0 : i = (BSF < 0) ? 128 - BSF : BSF;
4227 0 : fputc (i, fp);
4228 : /* Reserved: 3 bytes of 0. */
4229 0 : li_temp = 0;
4230 0 : fwrite (&li_temp, sizeof (char), 3, fp);
4231 0 : fputc (commentLen, fp);
4232 0 : fwrite (comment, sizeof (char), commentLen, fp);
4233 :
4234 : /* Now write section 2. */
4235 0 : if (f_grid) {
4236 0 : fputc (28, fp);
4237 0 : fputc (gridType, fp);
4238 0 : si_temp = gds->Nx;
4239 0 : FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4240 0 : si_temp = gds->Ny;
4241 0 : FWRITE_BIG (&si_temp, sizeof (short int), 1, fp);
4242 0 : li_temp = (sInt4) (gds->lat1 * 10000. + .5);
4243 0 : if (li_temp < 0) {
4244 0 : pbuf = 128;
4245 0 : li_temp = -1 * li_temp;
4246 : } else {
4247 0 : pbuf = 0;
4248 : }
4249 0 : pbufLoc = 7;
4250 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4251 0 : myAssert (pbufLoc == 8);
4252 0 : myAssert (pbuf == 0);
4253 0 : d_temp = 360 - gds->lon1;
4254 0 : if (d_temp < 0)
4255 0 : d_temp += 360;
4256 0 : if (d_temp > 360)
4257 0 : d_temp -= 360;
4258 0 : li_temp = (sInt4) (d_temp * 10000. + .5);
4259 0 : if (li_temp < 0) {
4260 0 : pbuf = 128;
4261 0 : li_temp = -1 * li_temp;
4262 : }
4263 0 : pbufLoc = 7;
4264 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4265 0 : myAssert (pbufLoc == 8);
4266 0 : myAssert (pbuf == 0);
4267 0 : d_temp = 360 - gds->orientLon;
4268 0 : if (d_temp < 0)
4269 0 : d_temp += 360;
4270 0 : if (d_temp > 360)
4271 0 : d_temp -= 360;
4272 0 : li_temp = (sInt4) (d_temp * 10000. + .5);
4273 0 : if (li_temp < 0) {
4274 0 : pbuf = 128;
4275 0 : li_temp = -1 * li_temp;
4276 : }
4277 0 : pbufLoc = 7;
4278 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4279 0 : myAssert (pbufLoc == 8);
4280 0 : myAssert (pbuf == 0);
4281 0 : li_temp = (sInt4) (gds->Dx * 1000. + .5);
4282 0 : FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4283 0 : li_temp = (sInt4) (gds->meshLat * 10000. + .5);
4284 0 : if (li_temp < 0) {
4285 0 : pbuf = 128;
4286 0 : li_temp = -1 * li_temp;
4287 : }
4288 0 : pbufLoc = 7;
4289 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 23, fp, &pbuf, &pbufLoc);
4290 0 : myAssert (pbufLoc == 8);
4291 0 : myAssert (pbuf == 0);
4292 0 : memset (buffer, 0, 6);
4293 0 : fwrite (buffer, sizeof (char), 6, fp);
4294 : }
4295 :
4296 : /* Now write section 3. */
4297 : /* Bitmap is not supported, skipping. */
4298 0 : myAssert (!f_bitmap);
4299 :
4300 : /* Now write section 4. */
4301 0 : FWRITE_ODDINT_BIG (&(sec4Len), 3, fp);
4302 0 : i = 0;
4303 0 : if (f_secMiss)
4304 0 : i |= 1;
4305 0 : if (f_primMiss)
4306 0 : i |= 2;
4307 0 : if (f_sndOrder)
4308 0 : i |= 4;
4309 0 : if (!f_simple)
4310 0 : i |= 8;
4311 0 : if (!f_grid)
4312 0 : i |= 16;
4313 0 : fputc (i, fp);
4314 0 : li_temp = DataLen;
4315 0 : FWRITE_BIG (&li_temp, sizeof (sInt4), 1, fp);
4316 0 : if (f_primMiss) {
4317 0 : FWRITE_BIG (&(li_primMiss), sizeof (sInt4), 1, fp);
4318 0 : if (f_secMiss) {
4319 0 : FWRITE_BIG (&(li_secMiss), sizeof (sInt4), 1, fp);
4320 : }
4321 : }
4322 0 : if (f_sndOrder) {
4323 0 : if (a1 < 0) {
4324 0 : pbuf = 128;
4325 0 : li_temp = -1 * a1;
4326 : } else {
4327 0 : pbuf = 0;
4328 0 : li_temp = a1;
4329 : }
4330 0 : pbufLoc = 7;
4331 0 : fileBitWrite (&(li_temp), sizeof (li_temp), 31, fp, &pbuf, &pbufLoc);
4332 0 : myAssert (pbufLoc == 8);
4333 0 : myAssert (pbuf == 0);
4334 0 : fileBitWrite (&mbit, sizeof (mbit), 5, fp, &pbuf, &pbufLoc);
4335 0 : if (b2 < 0) {
4336 0 : i = 1;
4337 0 : li_temp = -1 * b2;
4338 : } else {
4339 0 : i = 0;
4340 0 : li_temp = b2;
4341 : }
4342 0 : fileBitWrite (&i, sizeof (i), 1, fp, &pbuf, &pbufLoc);
4343 0 : myAssert (pbufLoc == 2);
4344 : fileBitWrite (&li_temp, sizeof (li_temp), (unsigned short int) mbit,
4345 0 : fp, &pbuf, &pbufLoc);
4346 : }
4347 0 : fileBitWrite (&nbit, sizeof (nbit), 5, fp, &pbuf, &pbufLoc);
4348 0 : if (overallMin < 0) {
4349 0 : i = 1;
4350 0 : li_temp = -1 * overallMin;
4351 : } else {
4352 0 : i = 0;
4353 0 : li_temp = overallMin;
4354 : }
4355 0 : fileBitWrite (&i, sizeof (i), 1, fp, &pbuf, &pbufLoc);
4356 : fileBitWrite (&li_temp, sizeof (li_temp), (unsigned short int) nbit, fp,
4357 0 : &pbuf, &pbufLoc);
4358 0 : fileBitWrite (&numGroup, sizeof (numGroup), 16, fp, &pbuf, &pbufLoc);
4359 0 : fileBitWrite (&ibit, sizeof (ibit), 5, fp, &pbuf, &pbufLoc);
4360 0 : fileBitWrite (&jbit, sizeof (jbit), 5, fp, &pbuf, &pbufLoc);
4361 0 : fileBitWrite (&kbit, sizeof (kbit), 5, fp, &pbuf, &pbufLoc);
4362 0 : for (i = 0; i < numGroup; i++) {
4363 0 : fileBitWrite (&(group[i].min), sizeof (sInt4),
4364 0 : (unsigned short int) ibit, fp, &pbuf, &pbufLoc);
4365 : }
4366 0 : for (i = 0; i < numGroup; i++) {
4367 0 : fileBitWrite (&(group[i].bit), sizeof (char),
4368 0 : (unsigned short int) jbit, fp, &pbuf, &pbufLoc);
4369 : }
4370 : #ifdef DEBUG
4371 0 : li_temp = 0;
4372 : #endif
4373 0 : for (i = 0; i < numGroup; i++) {
4374 0 : fileBitWrite (&(group[i].num), sizeof (sInt4),
4375 0 : (unsigned short int) kbit, fp, &pbuf, &pbufLoc);
4376 : #ifdef DEBUG
4377 0 : li_temp += group[i].num;
4378 : #endif
4379 : }
4380 : #ifdef DEBUG
4381 : /* Sanity check ! */
4382 0 : if (li_temp != DataLen) {
4383 : printf ("Total packed in groups %d != DataLen %d\n", li_temp,
4384 0 : DataLen);
4385 : }
4386 0 : myAssert (li_temp == DataLen);
4387 : #endif
4388 :
4389 0 : dataCnt = 0;
4390 : /* Start going through the data. Grid data has already been reordered.
4391 : * Already taken care of 2nd order differences. Data at this point should
4392 : * contain a stream of either 2nd order differences or 0 order
4393 : * differences. We have also removed overallMin from all non-missing
4394 : * elements in the stream */
4395 0 : if (f_secMiss) {
4396 0 : for (i = 0; i < numGroup; i++) {
4397 0 : max0 = (1 << group[i].bit) - 1;
4398 0 : max1 = (1 << group[i].bit) - 2;
4399 0 : for (j = 0; j < group[i].num; j++) {
4400 0 : if (Scaled[dataCnt] == li_primMiss) {
4401 0 : li_temp = max0;
4402 0 : } else if (Scaled[dataCnt] == li_secMiss) {
4403 0 : li_temp = max1;
4404 : } else {
4405 0 : li_temp = Scaled[dataCnt] - group[i].min;
4406 : }
4407 0 : fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4408 0 : &pbuf, &pbufLoc);
4409 0 : dataCnt++;
4410 : }
4411 : }
4412 0 : } else if (f_primMiss) {
4413 0 : for (i = 0; i < numGroup; i++) {
4414 : /* see what happens when bit == 0. */
4415 0 : max0 = (1 << group[i].bit) - 1;
4416 0 : for (j = 0; j < group[i].num; j++) {
4417 0 : if (group[i].bit != 0) {
4418 0 : if (Scaled[dataCnt] == li_primMiss) {
4419 0 : li_temp = max0;
4420 : } else {
4421 0 : li_temp = Scaled[dataCnt] - group[i].min;
4422 : }
4423 0 : fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4424 0 : &pbuf, &pbufLoc);
4425 : }
4426 0 : dataCnt++;
4427 : }
4428 : }
4429 : } else {
4430 0 : for (i = 0; i < numGroup; i++) {
4431 0 : for (j = 0; j < group[i].num; j++) {
4432 0 : li_temp = Scaled[dataCnt] - group[i].min;
4433 0 : if (group[i].bit != 0) {
4434 0 : fileBitWrite (&(li_temp), sizeof (sInt4), group[i].bit, fp,
4435 0 : &pbuf, &pbufLoc);
4436 : }
4437 0 : dataCnt++;
4438 : }
4439 : }
4440 : }
4441 0 : myAssert (dataCnt == DataLen);
4442 : /* flush the PutBit buffer... */
4443 0 : if (pbufLoc != 8) {
4444 0 : fputc ((int) pbuf, fp);
4445 : }
4446 :
4447 : /* Now write section 5. */
4448 0 : fwrite ("7777", sizeof (char), 4, fp);
4449 :
4450 : /* Deal with padding at end of record... */
4451 0 : for (i = 0; i < pad; i++) {
4452 0 : fputc (0, fp);
4453 : }
4454 : /* Now write FORTRAN record information */
4455 0 : FWRITE_BIG (&(recSize), sizeof (sInt4), 1, fp);
4456 :
4457 0 : free (Scaled);
4458 0 : free (group);
4459 0 : return 0;
4460 : }
|
