UltraScan III
us_solutedata.cpp
Go to the documentation of this file.
1 
3 #include "us_solutedata.h"
4 #include "us_ga_init.h"
5 #include "us_defines.h"
6 #include "us_math2.h"
7 #include "us_settings.h"
8 
9 // bucket vertex LessThan routine
10 bool buck_vx_lessthan( const bucket &buck1, const bucket &buck2 )
11 { // TRUE iff (sx1<sx2) || (sx1==sx2 && sy1<sy2)
12  return ( buck1.x_min < buck2.x_min ) ||
13  ( ( buck1.x_min == buck2.x_min ) && ( buck1.y_min < buck2.y_min ) );
14 }
15 
16 // Holds S_Solute data for US_GA_Initialize
17 US_SoluteData::US_SoluteData( void ) : QObject()
18 {
19  bndx = -1;
21 
22  allbucks.clear();
23 }
24 
26 {
27  int rc = 0;
28  bndx = -1;
29  allbucks.clear();
30  return rc;
31 }
32 
33 int US_SoluteData::sortBuckets( QList< bucket >* buks )
34 {
35  int rc = 0;
36  qSort( buks->begin(), buks->end() );
37  return rc;
38 }
39 
41 {
42  int rc = 0;
43  qSort( allbucks.begin(), allbucks.end() );
44  return rc;
45 }
46 
47 
49 {
50  int rc = 0;
51 
52  limitBucket( buk );
53  bndx = allbucks.size();
54  allbucks.append( buk );
55  return rc;
56 }
57 
58 int US_SoluteData::appendBucket( QRectF brect, QPointF bpnt,
59  qreal bconc, int bstat )
60 {
61  bucket buk = createBucket( brect, bpnt, bconc, bstat );
62 
63  return appendBucket( buk );
64 }
65 
66 int US_SoluteData::appendBucket( QRectF brect, QPointF bpnt, qreal bconc )
67 {
68  return appendBucket( brect, bpnt, bconc, 2 );
69 }
70 
71 int US_SoluteData::setBucket( int ix, QRectF brect, QPointF bpnt,
72  qreal bconc, int bstat )
73 {
74  bucket buk = createBucket( brect, bpnt, bconc, bstat );
75  int rc = 0;
76  int bsiz = allbucks.size();
77  int lx = bsiz - 1;
78  int jx = bndx + 1;
79  jx = ( jx < lx ) ? jx : lx;
80 
81  if ( ix > (-1) && ix < bsiz )
82  {
83  bndx = ix;
84  }
85  else if ( ix == (-1) )
86  {
87  bndx = lx;
88  }
89  else if ( ix == (-2) )
90  {
91  bndx = jx;
92  }
93  else if ( ix > lx )
94  {
95  bndx = bsiz;
96  allbucks.append( buk );
97  rc = 1;
98  }
99 
100  allbucks.replace( bndx, buk );
101  return rc;
102 }
103 
105 {
106  return allbucks.size();
107 }
108 
109 bucket US_SoluteData::createBucket( QRectF brect, QPointF bpnt,
110  qreal bconc, int bstat )
111 {
112  bucket buk;
113 
114  buk.x = bpnt.x();
115  buk.y = bpnt.y();
116  buk.x_min = brect.left();
117  buk.x_max = brect.right();
118  buk.y_min = brect.bottom();
119  buk.y_max = brect.top();
120  buk.conc = bconc;
121  buk.status = bstat;
122 
123  limitBucket( buk );
124 
125  return buk;
126 }
128 {
129  bucket buk;
130  return buk;
131 }
133 {
134  bucket buk;
135  return buk;
136 }
137 
139 {
140  bucket buk;
141  int bsiz = allbucks.size();
142  int lx = bsiz - 1;
143  int jx = bndx + 1;
144  jx = ( jx < lx ) ? jx : lx;
145 
146  if ( ix > (-1) && ix < bsiz )
147  { // get bucket at specified index and save the index
148  buk = allbucks.at( ix );
149  bndx = ix;
150  }
151  else if ( ix == (-1) )
152  { // get last bucket and set its index
153  buk = allbucks.at( lx );
154  bndx = lx;
155  }
156  else if ( ix == (-2) )
157  { // get next bucket and set its index
158  buk = allbucks.at( jx );
159  bndx = jx;
160  }
161  else if ( ix == (-3) )
162  { // get current bucket
163  buk = allbucks.at( bndx );
164  }
165 
166  return buk;
167 }
168 
170 {
171  return bucketAt( 0 );
172 }
173 
175 {
176  return bucketAt( -1 );
177 }
178 
180 {
181  return bucketAt( -2 );
182 }
183 
185 {
186  return bucketAt( -3 );
187 }
188 
190 {
191  bucket buk = bucketAt( ix );
192  qreal x1 = buk.x_min;
193  qreal x2 = buk.x_max;
194  qreal y1 = buk.y_min;
195  qreal y2 = buk.y_max;
196 
197  if ( x1 > x2 )
198  {
199  x1 = x2;
200  x2 = buk.x_min;
201  }
202 
203  if ( y1 > y2 )
204  {
205  y1 = y2;
206  y2 = buk.y_min;
207  }
208 
209  return QRectF( QPointF( x1, y1 ), QPointF( x2, y2 ) );
210 }
211 
212 QPointF US_SoluteData::bucketPoint( int ix, bool nearest )
213 {
214  bucket buk = bucketAt( ix ); // get specified bucket
215  qreal x1 = buk.x_min; // get vertices
216  qreal x2 = buk.x_max;
217  qreal y1 = buk.y_min;
218  qreal y2 = buk.y_max;
219 
220  QPointF mpt( ( x1 + x2 ) / 2.0, ( y1 + y2 ) / 2.0 );
221  QPointF spt = mpt; // get mid-point
222  QPointF& pt = spt;
223 
224  if ( nearest ) // get nearby solute point
225  findNearestPoint( pt );
226 
227  return pt;
228 }
229 
230 QPointF US_SoluteData::bucketPoint( int ix )
231 {
232  return bucketPoint( ix, true ); // solute point near bucket
233 }
234 
236 {
237  bucket buk = bucketAt( ix );
238  qreal xext = buk.x_max - buk.x_min;
239  xext = ( xext > 0.0 ) ? xext : -xext;
240  qreal yext = buk.y_max - buk.y_min;
241  yext = ( yext > 0.0 ) ? yext : -yext;
242  QSizeF siz( xext, yext );
243  return siz;
244 }
245 
246 QString US_SoluteData::bucketLine( int ix )
247 {
248  const QString cts[] = { "s", "ff0", "MW", "vbar", "D", "f" };
249 
250  bucket buk = bucketAt( ix );
251  int kx = ( ix < 0 ) ? bucketsCount() : ( ix + 1 );
252  QString slo = cts[ attr_x ] + "_min";
253  QString shi = cts[ attr_x ] + "_max";
254  QString klo = cts[ attr_y ] + "_min";
255  QString khi = cts[ attr_y ] + "_max";
256 
257  limitBucket( buk );
258 
259  if ( attr_x < attr_y )
260  return QString(
261  "Solute Bin %1: %2=%3, %4=%5, %6=%7, %8=%9" ).arg( kx )
262  .arg( slo ).arg( buk.x_min ).arg( shi ).arg( buk.x_max )
263  .arg( klo ).arg( buk.y_min ).arg( khi ).arg( buk.y_max );
264  else
265  return QString(
266  "Solute Bin %1: %2=%3, %4=%5, %6=%7, %8=%9" ).arg( kx )
267  .arg( klo ).arg( buk.y_min ).arg( khi ).arg( buk.y_max )
268  .arg( slo ).arg( buk.x_min ).arg( shi ).arg( buk.x_max );
269 }
270 
271 void US_SoluteData::setDistro( QList< S_Solute >* a_distro,
272  int a_ax, int a_ay, int a_az )
273 {
274  distro = a_distro;
275  attr_x = a_ax;
276  attr_y = a_ay;
277  attr_z = a_az;
278 }
279 
280 int US_SoluteData::findNearestPoint( QPointF& midpt )
281 {
282  qreal xm = midpt.x();
283  qreal ym = midpt.y();
284  qreal xn = xm;
285  qreal yn = ym;
286  qreal dl = xm * xm + ym * ym;
287  qreal dd = dl * 10.0;
288  int kx = -1;
289 
290  for ( int ii = 0; ii < distro->size(); ii++ )
291  {
292  qreal xs = distro->at( ii ).s;
293  qreal ys = distro->at( ii ).k;
294  qreal xd = xs - xm;
295  qreal yd = ys - ym;
296  dd = xd * xd + yd * yd;
297 
298  if ( dd < dl )
299  {
300  kx = ii;
301  dl = dd;
302  xn = xs;
303  yn = ys;
304  }
305  }
306 
307  if ( kx >= 0 )
308  {
309  midpt.setX( xn );
310  midpt.setY( yn );
311  }
312 
313  return kx;
314 }
315 
317 {
318  allbucks.removeAt( ix );
319  return allbucks.size();
320 }
321 
322 // do automatic calculation of bins
323 int US_SoluteData::autoCalcBins( int mxsols, qreal wsbuck, qreal hfbuck )
324 {
325 #define _MIN_VHR_ 0.01
326  QList< bucket > tbuk1;
327  QList< bucket > tbuk2;
328  QList< bucket >* buks1 = &tbuk1;
329  QList< bucket >* buks2 = &tbuk2;
330  QList< qreal > cvals;
331  bucket buk;
332  int bukflip = 0;
333  int nisols = distro->size();
334  int nssols = nisols;
335  int ntsols = mxsols;
336  qreal cval;
337  qreal cutlo;
338  qreal wbuckh = wsbuck / 2.0;
339  qreal hbuckh = hfbuck / 2.0;
340 
341  if ( ntsols < 2 || ntsols > nisols )
342  {
343  ntsols = nisols;
344  }
345 
346 
347  // initialize work list from solute distribution and get prelim stats
348  for ( int jj = 0; jj < nisols; jj++ )
349  {
350  qreal sval = distro->at( jj ).s;
351  qreal kval = distro->at( jj ).k;
352  cval = distro->at( jj ).c;
353  buk.x = sval;
354  buk.x_min = sval - wbuckh;
355  buk.x_max = sval + wbuckh;
356  buk.y = kval;
357  buk.y_min = kval - hbuckh;
358  buk.y_max = kval + hbuckh;
359  buk.conc = cval;
360  buk.status = 2;
361  limitBucket( buk );
362  buks1->append( buk );
363  cvals.append( cval );
364  }
365 
366  // sort the concentrate values to find cut-off values
367  qSort( cvals.begin(), cvals.end() );
368 
369  // (possibly) pare down the list based on cut-off C values
370  cutlo = cvals.at( nisols - ntsols ); // low val in sorted concen vals
371 
372  for ( int jj = 0; jj < nssols; jj++ )
373  {
374  buk = buks1->at( jj );
375  cval = buk.conc;
376  if ( cval >= cutlo )
377  {
378  buks2->append( buk );
379  }
380  }
381 
382  nssols = buks2->size();
383  buks1->clear();
384  buks1 = bukflip ? &tbuk1 : &tbuk2;
385  buks2 = bukflip ? &tbuk2 : &tbuk1;
386  bukflip = 1 - bukflip;
387 
388  // now perform as many passes as needed to handle all overlaps
389  bool overs = true;
390  int npass = 0;
391  int tstat = 1; // set up to skip reduced bins during 1st few passes
392 
393  while ( overs )
394  {
395  int novls = 0;
396  npass++;
397  buks2->append( buks1->at( 0 ) );
398 
399  for ( int jj = 1; jj < nssols; jj++ )
400  { // examine each of the current buckets
401  qreal horzr;
402  qreal vertr;
403  buk = buks1->at( jj );
404 
405  // don't break up bins that are already reduced for 1st few passes
406  if ( buk.status == tstat )
407  {
408  buks2->append( buk );
409  continue;
410  }
411 
412  qreal cx1 = buk.x_min;
413  qreal cx2 = buk.x_max;
414  qreal cy1 = buk.y_min;
415  qreal cy2 = buk.y_max;
416 
417  for ( int kk = 0; kk < jj; kk++ )
418  { // compare current to each previous bucket
419  bucket buk2 = buks1->at( kk );
420  qreal px1 = buk2.x_min;
421  qreal px2 = buk2.x_max;
422  qreal py1 = buk2.y_min;
423  qreal py2 = buk2.y_max;
424 
425  if ( cx1 < px2 )
426  { // possible horizontal overlap
427 
428  if ( cy1 < py2 && cy2 > py2 )
429  { // current overlaps in its lower left
430  novls++;
431  buk.status = 1;
432  buk2 = buk;
433  horzr = ( px2 - cx1 ) / wsbuck;
434  vertr = ( py2 - cy1 ) / hfbuck;
435  if ( vertr < horzr )
436  { // split vertically
437  buk2.x_min = qMin( cx1, px2 );
438  buk2.x_max = qMax( px2, cx1 );
439  buk2.y_min = qMin( py2, cy2 );
440  buk2.y_max = qMax( cy2, py2 );
441  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
442  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
443  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
444  buks2->append( buk2 );
445 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
446 
447  buk.x_min = qMin( px2, cx2 );
448  buk.x_max = qMax( cx2, px2 );
449  buk.y_min = cy1;
450  buk.y_max = cy2;
451  }
452 
453  else
454  { // split horizontally
455  buk2.x_min = cx1;
456  buk2.x_max = cx2;
457  buk2.y_min = qMin( py2, cy2 );
458  buk2.y_max = qMax( cy2, py2 );
459  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
460  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
461  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
462  buks2->append( buk2 );
463 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
464 
465  buk.x_min = qMin( px2, cx2 );
466  buk.x_max = qMax( cx2, px2 );
467  buk.y_min = qMin( cy1, py2 );
468  buk.y_max = qMax( py2, cy1 );
469  }
470 DbgLv(2) << " LL OVL: novls " << novls;
471  break;
472  }
473 
474  else if ( cy2 > py1 && cy1 < py1 )
475  { // current overlaps in its upper left
476  novls++;
477  buk.status = 1;
478  buk2 = buk;
479  horzr = ( px2 - cx1 ) / wsbuck;
480  vertr = ( cy2 - py1 ) / hfbuck;
481  if ( vertr < horzr )
482  { // split vertically
483  buk2.x_min = qMin( cx1, px2 );
484  buk2.x_max = qMax( px2, cx1 );
485  buk2.y_min = qMin( cy1, py1 );
486  buk2.y_max = qMax( py1, cy1 );
487  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
488  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
489  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
490  buks2->append( buk2 );
491 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
492 
493  buk.x_min = qMin( px2, cx2 );
494  buk.x_max = qMax( cx2, px2 );
495  buk.y_min = cy1;
496  buk.y_max = cy2;
497  }
498 
499  else
500  { // split horizontally
501  buk2.x_min = qMin( px2, cx2 );
502  buk2.x_max = qMax( cx2, px2 );
503  buk2.y_min = qMin( py1, cy2 );
504  buk2.y_max = qMax( cy2, py1 );
505  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
506  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
507  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
508  buks2->append( buk2 );
509 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
510 
511  buk.x_min = qMin( cx1, cx2 );
512  buk.x_max = qMax( cx2, cx1 );
513  buk.y_min = qMin( cy1, py1 );
514  buk.y_max = qMax( py1, cy1 );
515  }
516 DbgLv(2) << " UL OVL: novls " << novls;
517  break;
518  }
519 
520  else if ( cy1 >= py1 && cy2 <= py2 )
521  { // current overlaps in its middle left
522  novls++;
523  buk.status = 1;
524  buk.x_min = qMin( px2, cx2 );
525  buk.x_max = qMax( cx2, px2 );
526 DbgLv(2) << " UL OVL: novls " << novls;
527  break;
528  }
529 
530  else if ( cy2 > py2 && cy1 < py1 )
531  { // current overlaps in both upper and lower left
532  novls++;
533  buk.status = 1;
534  buk2 = buk;
535  horzr = ( px2 - cx1 ) / wsbuck;
536  vertr = ( py2 - py1 ) / hfbuck;
537  if ( vertr < horzr )
538  { // split vertically (into 3 total parts)
539  buk2.x_min = qMin( cx1, px2 ); // top
540  buk2.x_max = qMax( px2, cx1 );
541  buk2.y_min = qMin( py1, cy2 );
542  buk2.y_max = qMax( cy2, py1 );
543  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
544  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
545  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
546  buks2->append( buk2 );
547 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
548 
549  buk2.y_min = qMin( cy1, py1 ); // bottom
550  buk2.y_max = qMax( py1, cy1 );
551  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
552  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
553  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
554  buks2->append( buk2 );
555 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
556 
557  buk.x_min = qMin( px2, cx2 ); // right
558  buk.x_max = qMax( cx2, px2 );
559  buk.y_min = cy1;
560  buk.y_max = cy2;
561  }
562 
563  else
564  { // split horizontally (into 3 total parts)
565  buk2.x_min = qMin( cx1, cx2 ); // top
566  buk2.x_max = qMax( cx2, cx1 );
567  buk2.y_min = qMin( py2, cy2 );
568  buk2.y_max = qMax( cy2, py2 );
569  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
570  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
571  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
572  buks2->append( buk2 );
573 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
574 
575  buk2.y_min = qMin( cy1, py1 ); // bottom
576  buk2.y_max = qMax( py1, cy1 );
577  horzr = ( buk2.x_max - buk2.x_min ) / wsbuck;
578  vertr = ( buk2.y_max - buk2.y_min ) / hfbuck;
579  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
580  buks2->append( buk2 );
581 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
582 
583  buk.x_min = qMin( px2, cx2 ); // middle
584  buk.x_max = qMax( cx2, px2 );
585  buk.y_min = qMin( py1, py2 );
586  buk.y_max = qMax( py2, py1 );
587  }
588 DbgLv(2) << " UL OVL: novls " << novls;
589  break;
590  }
591 
592  }
593 
594  else if ( cx1 == px1 )
595  { // possible pure vertical overlap
596 
597  if ( cy1 < py1 && cy2 > py1 )
598  { // current overlaps in its upper part
599  novls++;
600  buk.status = 1;
601  buk.y_min = qMin( cy1, py1 );
602  buk.y_max = qMax( py1, cy1 );
603 DbgLv(2) << " UV OVL: novls " << novls;
604  break;
605  }
606  }
607  }
608  horzr = ( buk.x_max - buk.x_min ) / wsbuck;
609  vertr = ( buk.y_max - buk.y_min ) / hfbuck;
610  if ( horzr > _MIN_VHR_ && vertr > _MIN_VHR_ )
611  buks2->append( buk );
612 else DbgLv(2) << "BUCKET TOO THIN H,V " << horzr << "," << vertr;
613  }
614  // get new bucket count; flip-flop input,output lists
615  nssols = buks2->size();
616  buks1->clear();
617  buks1 = bukflip ? &tbuk1 : &tbuk2;
618  buks2 = bukflip ? &tbuk2 : &tbuk1;
619  bukflip = 1 - bukflip;
620 
621  if ( novls > 0 )
622  { // if there were overlaps, re-sort buckets by vertex for next pass
623  qSort( buks1->begin(), buks1->end(), buck_vx_lessthan );
624  }
625 
626  else if ( tstat == 1 )
627  { // start new set of passes where reduced buckets can be reduced more
628  tstat = 4;
629  }
630 
631  else
632  { // otherwise, we must be done!
633  overs = false;
634  }
635  }
636 
637  // do a re-sort based on center point
638  qSort( buks1->begin(), buks1->end() );
639 
640  // copy the final bucket list to the master
641 
642  allbucks.clear();
643 
644  for ( int jj = 0; jj < nssols; jj++ )
645  {
646  allbucks.append( buks1->at( jj ) );
647  }
648 
649  tbuk1.clear();
650  tbuk2.clear();
651  cvals.clear();
652  return nssols;
653 }
654 
655 // save bucket information to file for use by GA
656 int US_SoluteData::saveGAdata( QString& fname, int xtype, int ytype,
657  int ztype, double fixval )
658 {
659  int rc = 0;
660  int nbuk = allbucks.size();
661  QString line;
662  bucket buk;
663  const QString cts[] = { "s", "ff0", "MW", "vbar", "D", "f" };
664  const char hfmt[] = "%d %d %d %d %.3e"; // Header format
665  const char lfmt[] = "%.3e, %.3e, %.3e, %.3e"; // Lines format
666  const char* cffmt = "%.3e"; // Comment fixed format
667 
668  QFile fileo( fname );
669 
670  if ( fileo.open( QIODevice::WriteOnly | QIODevice::Text ) )
671  {
672  QTextStream ts( &fileo );
673 
674  // Line 1 with count,xtype,ytype,ztype,fixed
675  line.sprintf( hfmt, nbuk, xtype, ytype, ztype, fixval );
676  line += " # buckets=" + QString::number( nbuk )
677  + " x=" + cts[ xtype ] + " y=" + cts[ ytype ]
678  + " fixed=" + cts[ ztype ] + "="
679  + ( ( ztype == US_GA_Initialize::ATTR_V )
680  ? "0.0=(data_set_vbar)"
681  : QString().sprintf( cffmt, fixval ) );
682  ts << line << endl;
683 
684  for ( int jj = 0; jj < nbuk; jj++ )
685  {
686  buk = allbucks.at( jj );
687 
688  limitBucket( buk );
689 
690  line.sprintf( lfmt, buk.x_min, buk.x_max, buk.y_min, buk.y_max );
691  ts << line << endl;
692  }
693  fileo.close();
694  }
695  else
696  rc = 1;
697 
698  return rc;
699 }
700 
701 
702 // build the data lists for Monte Carlo analysis
704 {
705  int rc = 0;
706  int nsol = distro->size();
707  int nbuk = allbucks.size();
708  bucket buk; // bucket record
709  S_Solute d_sol; // solute record
710  SimComp simc; // simulation component record
711  QList< SimComp > bcomp; // sim component list
712  qreal bsmin; // bucket vertices
713  qreal bsmax;
714  qreal bkmin;
715  qreal bkmax;
716  qreal ssval; // component s,f_f0 values
717  qreal sfval;
718 
719  // build component list from solute lists
720  component.clear();
721 
722  for ( int jj = 0; jj < nsol; jj++ )
723  {
724  d_sol = distro->at( jj ); // get solute record ("x" is "s")
725 
726  simc.s = d_sol.si * 1e-13; // compose simulation component
727  simc.w = d_sol.w;
728  simc.f = d_sol.ki;
729  simc.c = d_sol.c;
730  simc.d = d_sol.d;
731  simc.v = d_sol.v;
732 
733  component.append( simc ); // and add it to list
734  }
735 
736  // build list of component solute data for each bucket
737  MC_solute.clear();
738 
739  for ( int jj = 0; jj < nbuk; jj++ )
740  {
741  // get bucket dimensions
742  buk = allbucks.at( jj ); // get bucket and its vertices
743  bsmin = buk.x_min;
744  bsmax = buk.x_max;
745  bkmin = buk.y_min;
746  bkmax = buk.y_max;
747  bcomp.clear();
748 
749  for ( int kk = 0; kk < nsol; kk++ )
750  { // add solute points that fall within bin dimensions
751  ssval = distro->at( kk ).s; // "x,y" of solute point
752  sfval = distro->at( kk ).k;
753  if ( ssval >= bsmin && ssval <= bsmax &&
754  sfval >= bkmin && sfval <= bkmax )
755  { // solute is in this bin: save component for bin
756  bcomp.append( component.at( kk ) );
757  }
758  }
759 
760  MC_solute.append( bcomp );
761  }
762 
763  return rc;
764 }
765 
766 // complete analysis and report Monte Carlo statistics
767 int US_SoluteData::reportDataMC( QString& fname, int mc_iters )
768 {
769  int rc = 0;
770  int nbuk = MC_solute.size();
771  bucket buk; // bucket record
772  QList< SimComp > bcomp; // sim component list
773  QList< double > vals;
774 
775  QFile fileo( fname );
776 
777  if ( fileo.open( QIODevice::WriteOnly | QIODevice::Text ) )
778  { // output Monte Carlo statistics to a report file
779  QList< qreal > valus;
780  QList< qreal > concs;
781  QList< qreal > csums;
782  QTextStream ts( &fileo );
783  QString str1;
784 
785  ts << "*****************************************"
786  "**********************************\n\n";
787  if ( mc_iters > 1 )
788  ts << "Monte Carlo Analysis Statistical Results "
789  "(from Genetic Algorithm Analysis):\n\n";
790  else
791  ts << "Global Analysis Statistical Results "
792  "(from Genetic Algorithm Analysis):\n\n";
793  ts << "*****************************************"
794  "**********************************\n\n";
795  ts << "Summary:\n";
796 
797  if ( nbuk < 1 )
798  {
799  return 2;
800  }
801  qreal vsum = 0.0;
802  qreal vsiz = 0.0;
803 
804  for ( int kk = 0; kk < nbuk; kk++ )
805  { // accumulate statistics for each bin
806  bcomp = MC_solute.at( kk );
807  ts << "\nSolute " << ( kk + 1 ) << ":\n";
808  int ksol = bcomp.size();
809  vsiz = (double)ksol;
810 
811  if ( ksol < 3 )
812  { // Summary prints for a bin that has only a point or two
813  ts << tr( "This solute bin does not have sufficient points to"
814  "\ncalculate meaningful statistics.\n" );
815 
816  double tconc = 0.0;
817  for ( int jj = 0; jj < ksol; jj++ )
818  tconc += bcomp.at( jj ).c;
819 
820  ts << tr( "Molecular weight: " );
821  vsum = 0.0;
822  for ( int jj = 0; jj < ksol; jj++ )
823  vsum += ( bcomp.at( jj ).w * bcomp.at( jj ).c );
824  ts << ( vsum / tconc ) << endl;
825 
826  ts << tr( "Sedimentation coefficient: " );
827  vsum = 0.0;
828  for ( int jj = 0; jj < ksol; jj++ )
829  vsum += ( bcomp.at( jj ).s * bcomp.at( jj ).c );
830  ts << ( vsum / tconc ) << endl;
831 
832  ts << tr( "Diffusion coefficient: " );
833  vsum = 0.0;
834  for ( int jj = 0; jj < ksol; jj++ )
835  vsum += ( bcomp.at( jj ).d * bcomp.at( jj ).c );
836  ts << ( vsum / tconc ) << endl;
837 
838  ts << tr( "Frictional ratio: " );
839  vsum = 0.0;
840  for ( int jj = 0; jj < ksol; jj++ )
841  vsum += ( bcomp.at( jj ).f * bcomp.at( jj ).c );
842  ts << ( vsum / tconc ) << endl;
843 
844  ts << tr( "Partial specific volume: " );
845  vsum = 0.0;
846  for ( int jj = 0; jj < ksol; jj++ )
847  vsum += ( bcomp.at( jj ).v * bcomp.at( jj ).c );
848  ts << ( vsum / tconc ) << endl;
849 
850  ts << tr( "Partial concentration: " );
851  vsum = tconc;
852  ts << ( vsum / vsiz ) << endl;
853  csums.append( vsum );
854  }
855 
856  else
857  { // Summary prints for the typical bin with many points
858  qreal vtotal = 0.0;
859 
860  concs.clear();
861  for ( int jj = 0; jj < ksol; jj++ )
862  concs.append( bcomp.at( jj ).c );
863 
864  valus.clear();
865  for ( int jj = 0; jj < ksol; jj++ )
866  valus.append( bcomp.at( jj ).w );
867  outputStats( ts, valus, concs, false,
868  tr( "Molecular weight: " ) );
869 
870  valus.clear();
871  for ( int jj = 0; jj < ksol; jj++ )
872  valus.append( bcomp.at( jj ).s );
873  outputStats( ts, valus, concs, false,
874  tr( "Sedimentation coefficient: " ) );
875 
876  valus.clear();
877  for ( int jj = 0; jj < ksol; jj++ )
878  valus.append( bcomp.at( jj ).d );
879  outputStats( ts, valus, concs, false,
880  tr( "Diffusion coefficient: " ) );
881 
882  valus.clear();
883  for ( int jj = 0; jj < ksol; jj++ )
884  valus.append( bcomp.at( jj ).f );
885  outputStats( ts, valus, concs, false,
886  tr( "Frictional ratio: " ) );
887 
888  valus.clear();
889  for ( int jj = 0; jj < ksol; jj++ )
890  valus.append( bcomp.at( jj ).v );
891  outputStats( ts, valus, concs, false,
892  tr( "Partial specific volume: " ) );
893 
894  for ( int jj = 0; jj < ksol; jj++ )
895  vtotal += bcomp.at( jj ).c;
896  csums.append( vtotal );
897  ts << tr( "Partial concentration: " ) <<
898  str1.sprintf( " %6.4e\n", vtotal );
899 
900  }
901  }
902 
903  qreal concsum = 0.0;
904 
905  for ( int jj = 0; jj < distro->size(); jj++ )
906  concsum += distro->at( jj ).c;
907 
908  ts << tr( "\nMonte Carlo iterations: " ) << mc_iters;
909  if ( mc_iters < 2 )
910  ts << tr( " (Global or High-Solutes-per-Bucket)" );
911  ts << tr( "\n\nRelative Concentrations:\n\n" );
912  ts << tr( "Total concentration: " ) << concsum << " OD\n";
913 
914  for ( int jj = 0; jj < csums.size(); jj++ )
915  {
916  qreal pcconc = 100.0 * csums.at( jj ) / concsum;
917  ts << tr( "Relative percentage of Solute " ) << ( jj + 1 )
918  << ": " << QString().sprintf( "%7.3f", pcconc ) << " %\n";
919  }
920 
921  ts << tr( "\n\nDetailed Results:\n" );
922 
923  for ( int kk = 0; kk < nbuk; kk++ )
924  { // output detailed statistics for all the bins
925  bcomp = MC_solute.at( kk );
926  ts << "\n*****************************************"
927  "**********************************";
928  ts << "\nSolute " << ( kk + 1 ) << ":";
929  int ksol = bcomp.size();
930  vsiz = (double)ksol;
931 
932  if ( ksol < 3 )
933  { // just print the values for a sparse bin
934  ts << tr( "\nThis solute bin does not have sufficient points to"
935  "\ncalculate a meaningful distribution.\n\n" );
936 
937  ts << tr( "Molecular weight: " );
938  for ( int jj = 0; jj < ksol; jj++ )
939  ts << bcomp.at( jj ).w << " ";
940  ts << endl;
941 
942  ts << tr( "Sedimentation coefficient: " );
943  for ( int jj = 0; jj < ksol; jj++ )
944  ts << bcomp.at( jj ).s << " ";
945  ts << endl;
946 
947  ts << tr( "Diffusion coefficient: " );
948  for ( int jj = 0; jj < ksol; jj++ )
949  ts << bcomp.at( jj ).d << " ";
950  ts << endl;
951 
952  ts << tr( "Frictional ratio: " );
953  for ( int jj = 0; jj < ksol; jj++ )
954  ts << bcomp.at( jj ).f << " ";
955  ts << endl;
956 
957  ts << tr( "Partial specific volume: " );
958  for ( int jj = 0; jj < ksol; jj++ )
959  ts << bcomp.at( jj ).v << " ";
960  ts << endl;
961 
962  ts << tr( "Concentration: " );
963  for ( int jj = 0; jj < ksol; jj++ )
964  ts << bcomp.at( jj ).c << " ";
965  ts << endl;
966  }
967 
968  else
969  { // calculate and output detailed statistics for the bin
970  concs.clear();
971  for ( int jj = 0; jj < ksol; jj++ )
972  concs.append( bcomp.at( jj ).c );
973 
974  valus.clear();
975  for ( int jj = 0; jj < ksol; jj++ )
976  valus.append( bcomp.at( jj ).w );
977  outputStats( ts, valus, concs, true,
978  tr( "Molecular Weight" ) );
979 
980  valus.clear();
981  for ( int jj = 0; jj < ksol; jj++ )
982  valus.append( bcomp.at( jj ).s );
983  outputStats( ts, valus, concs, true,
984  tr( "Sedimentation Coefficient" ) );
985 
986  valus.clear();
987  for ( int jj = 0; jj < ksol; jj++ )
988  valus.append( bcomp.at( jj ).d );
989  outputStats( ts, valus, concs, true,
990  tr( "Diffusion Coefficient" ) );
991 
992  valus.clear();
993  for ( int jj = 0; jj < ksol; jj++ )
994  valus.append( bcomp.at( jj ).f );
995  outputStats( ts, valus, concs, true,
996  tr( "Frictional Ratio" ) );
997 
998  valus.clear();
999  for ( int jj = 0; jj < ksol; jj++ )
1000  valus.append( bcomp.at( jj ).v );
1001  outputStats( ts, valus, concs, true,
1002  tr( "Partial Specific Volume" ) );
1003  }
1004  }
1005 
1006  fileo.close();
1007  }
1008  else
1009  rc = 1;
1010  return rc;
1011 }
1012 
1013 // output statistics for an array of values of a given type
1014 void US_SoluteData::outputStats( QTextStream& ts, QList< qreal >& vals,
1015  QList< qreal >& concs, bool details, QString title )
1016 {
1017  QString str1;
1018  QString str2;
1019  int nvals = vals.size();
1020  int nbins = 50;
1021  qreal vsiz = (qreal)nvals;
1022  qreal binsz = 50.0;
1023  qreal vlo = 9.9e30;
1024  qreal vhi = -9.9e30;
1025 
1026  QVector< qreal > xpvec( qMax( nvals, nbins ) );
1027  QVector< qreal > ypvec( qMax( nvals, nbins ) );
1028 
1029  qreal *xplot = xpvec.data();
1030  qreal *yplot = ypvec.data();
1031  qreal vsum = 0.0;
1032  qreal vm2 = 0.0;
1033  qreal vm3 = 0.0;
1034  qreal vm4 = 0.0;
1035  qreal vmean;
1036  qreal mode_cen;
1037  qreal mode_lo;
1038  qreal mode_hi;
1039  qreal conf99lo;
1040  qreal conf99hi;
1041  qreal conf95lo;
1042  qreal conf95hi;
1043  qreal skew;
1044  qreal kurto;
1045  qreal vmedi;
1046  qreal slope;
1047  qreal vicep;
1048  qreal sigma;
1049  qreal corr;
1050  qreal sdevi;
1051  qreal sderr;
1052  qreal vari;
1053  qreal area;
1054  qreal bininc;
1055  qreal val;
1056  qreal conc;
1057  qreal vctot = 0.0;
1058 
1059  // get basic min,max,mean information
1060 
1061  for ( int jj = 0; jj < nvals; jj++ )
1062  {
1063  val = vals.at( jj );
1064  conc = concs.at( jj );
1065  vsum += ( val * conc );
1066  vctot += conc;
1067  vlo = qMin( vlo, val );
1068  vhi = qMax( vhi, val );
1069  xplot[jj] = (qreal)jj;
1070  yplot[jj] = val;
1071  }
1072 
1073  vmean = vsum / vctot;
1074 
1075  // get difference information
1076 
1077  for ( int jj = 0; jj < nvals; jj++ )
1078  {
1079  val = vals.at( jj );
1080  qreal dif = val - vmean;
1081  qreal dsq = dif * dif;
1082  vm2 += dsq; // diff squared
1083  vm3 += ( dsq * dif ); // cubed
1084  vm4 += ( dsq * dsq ); // to the 4th
1085  }
1086 
1087  vm2 /= vsiz;
1088  vm3 /= vsiz;
1089  vm4 /= vsiz;
1090  skew = vm3 / pow( vm2, 1.5 );
1091  kurto = vm4 / pow( vm2, 2.0 ) - 3.0;
1092  vmedi = ( vlo + vhi ) / 2.0;
1093 
1094  // do line fit (mainly for corr value)
1095 
1096  US_Math2::linefit( &xplot, &yplot, &slope, &vicep, &sigma, &corr, nvals );
1097 
1098  // standard deviation and error
1099 
1100  sdevi = pow( vm2, 0.5 );
1101  sderr = sdevi / pow( vsiz, 0.5 );
1102  vari = vm2;
1103  area = 0.0;
1104 
1105  bininc = ( vhi - vlo ) / binsz;
1106 
1107  // mode and confidence
1108 
1109  for ( int ii = 0; ii < nbins; ii++ )
1110  {
1111  xplot[ii] = vlo + bininc * (qreal)ii;
1112  yplot[ii] = 0.0;
1113 
1114  for ( int jj = 0; jj < nvals; jj++ )
1115  {
1116  val = vals.at( jj );
1117 
1118  if ( val >= xplot[ ii ] && val < ( xplot[ ii ] + bininc ) )
1119  {
1120  yplot[ii] += ( concs.at( jj ) );
1121  }
1122  }
1123 
1124  area += yplot[ ii ] * bininc;
1125  }
1126 
1127  double fvdif = qAbs( ( vhi - vlo ) / vlo );
1128  bool is_constant = ( fvdif < 0.0001 );
1129  val = -1.0;
1130  int thisb = 0;
1131 
1132  for ( int ii = 0; ii < nbins; ii++ )
1133  {
1134  if ( yplot[ii] > val )
1135  {
1136  val = yplot[ ii ];
1137  thisb = ii;
1138  }
1139  }
1140 
1141  mode_lo = xplot[ thisb ];
1142  mode_hi = mode_lo + bininc;
1143  mode_cen = ( mode_lo + mode_hi ) / 2.0;
1144  conf99lo = vmean - 2.576 * sdevi;
1145  conf99hi = vmean + 2.576 * sdevi;
1146  conf95lo = vmean - 1.960 * sdevi;
1147  conf95hi = vmean + 1.960 * sdevi;
1148 
1149  if ( details && is_constant )
1150  {
1151  ts << "\n\n" << tr( "Results for the " ) << title << ":\n\n";
1152 
1153  if ( vhi == vlo )
1154  {
1155  ts << tr( "Constant Value: " )
1156  << str1.sprintf( "%6.4e\n", vhi );
1157  }
1158 
1159  else
1160  {
1161  ts << tr( "Maximum Value: " )
1162  << str1.sprintf( "%6.4e\n", vhi );
1163  ts << tr( "Minimum Value: " )
1164  << str1.sprintf( "%6.4e\n", vlo );
1165  ts << tr( "(Nearly) Constant Value: " )
1166  << str1.sprintf( "%6.4e\n", vmean );
1167  }
1168  }
1169 
1170  else if ( details )
1171  { // Details
1172  ts << "\n\n" << tr( "Results for the " ) << title << ":\n\n";
1173  ts << tr( "Maximum Value: " )
1174  << str1.sprintf( "%6.4e\n", vhi );
1175  ts << tr( "Minimum Value: " )
1176  << str1.sprintf( "%6.4e\n", vlo );
1177  ts << tr( "Mean Value: " )
1178  << str1.sprintf( "%6.4e\n", vmean );
1179  ts << tr( "Median Value: " )
1180  << str1.sprintf( "%6.4e\n", vmedi );
1181  ts << tr( "Skew Value: " )
1182  << str1.sprintf( "%6.4e\n", skew );
1183  ts << tr( "Kurtosis Value: " )
1184  << str1.sprintf( "%6.4e\n", kurto );
1185  ts << tr( "Lower Mode Value: " )
1186  << str1.sprintf( "%6.4e\n", mode_lo );
1187  ts << tr( "Upper Mode Value: " )
1188  << str1.sprintf( "%6.4e\n", mode_hi );
1189  ts << tr( "Mode Center: " )
1190  << str1.sprintf( "%6.4e\n", mode_cen );
1191  ts << tr( "95% Confidence Limits: " )
1192  << str1.sprintf( "%6.4e, -%6.4e\n",
1193  ( conf95hi - mode_cen ), ( mode_cen - conf95lo ) );
1194  ts << tr( "99% Confidence Limits: " )
1195  << str1.sprintf( "%6.4e, -%6.4e\n",
1196  ( conf99hi - mode_cen ), ( mode_cen - conf99lo ) );
1197  ts << tr( "Standard Deviation: " )
1198  << str1.sprintf( "%6.4e\n", sdevi );
1199  ts << tr( "Standard Error: " )
1200  << str1.sprintf( "%6.4e\n", sderr );
1201  ts << tr( "Variance: " )
1202  << str1.sprintf( "%6.4e\n", vari );
1203  ts << tr( "Correlation Coefficent: " )
1204  << str1.sprintf( "%6.4e\n", corr );
1205  ts << tr( "Number of Bins: " )
1206  << qRound( binsz ) << "\n";
1207  ts << tr( "Distribution Area: " )
1208  << str1.sprintf( "%6.4e\n", area );
1209 
1210  str1.sprintf( "%e", conf95lo ).append( tr( " (low), " ) );
1211  str2.sprintf( "%e", conf95hi ).append( tr( " (high)\n" ) );
1212  ts << tr( "95% Confidence Interval: " ) << str1 << str2;
1213 
1214  str1.sprintf( "%e", conf99lo ).append( tr( " (low), " ) );
1215  str2.sprintf( "%e", conf99hi ).append( tr( " (high)\n" ) );
1216  ts << tr( "99% Confidence Interval: " ) << str1 << str2;
1217  }
1218 
1219  else if ( is_constant )
1220  { // Summary (where value is constant)
1221  ts << title << str1.sprintf( " %6.4e (**constant**)\n", vmean );
1222  }
1223 
1224  else
1225  { // Summary
1226  ts << title << str1.sprintf( " %6.4e (%6.4e, %6.4e)\n",
1227  vmean, conf95lo, conf95hi );
1228  }
1229 }
1230 
1231 // Insure vertexes of a bucket do not exceed physically possible limits
1233 {
1234  if ( attr_x != US_GA_Initialize::ATTR_D &&
1236  {
1237  if ( buk.x_min > 0.0 )
1238  { // All-positive s's start at 0.1 at least
1239  buk.x_min = qMax( 0.1, buk.x_min );
1240  buk.x_max = qMax( ( buk.x_min + 0.0001 ), buk.x_max );
1241  }
1242 
1243  else if ( buk.x_max <= 0.0 )
1244  { // All-negative s's end at -0.1 at most
1245  buk.x_max = qMin( -0.1, buk.x_max );
1246  buk.x_min = qMin( ( buk.x_max - 0.0001 ), buk.x_min );
1247  }
1248 
1249  else if ( ( buk.x_min + buk.x_max ) >= 0.0 )
1250  { // Mostly positive clipped to all positive starting at 0.1
1251  buk.x_min = 0.1;
1252  buk.x_max = qMax( 0.2, buk.x_max );
1253  }
1254 
1255  else
1256  { // Mostly negative clipped to all negative ending at -0.1
1257  buk.x_min = qMin( -0.2, buk.x_min );
1258  buk.x_max = -0.1;
1259  }
1260  }
1261 
1262  if ( attr_y != US_GA_Initialize::ATTR_D &&
1264  {
1266  buk.y_min = qMax( 1.0, buk.y_min );
1267  else
1268  buk.y_min = qMax( 0.1, buk.y_min );
1269  buk.y_max = qMax( ( buk.y_min + 0.0001 ), buk.y_max );
1270  }
1271 }
1272 
1273 // Count the number of overlaps in the current list of buckets
1275 {
1276  int nbuks = allbucks.size();
1277  int novlps = 0;
1278  QList< QRectF > bucket_rects;
1279 DbgLv(1) << "countOv nbuks" << nbuks;
1280 DbgLv(1) << "CountO dbg_level>=1";
1281 
1282  // Create the list of bucket rectangles
1283  for ( int ii = 0; ii < nbuks; ii++ )
1284  bucket_rects << bucketRect( ii );
1285 
1286  // Count the number of overlaps
1287  for ( int ii = 0; ii < nbuks; ii++ )
1288  {
1289  QRectF bukrect = bucket_rects[ ii ];
1290 
1291  // Compare this bucket to all buckets that follow it
1292  for ( int jj = ii + 1; jj < nbuks; jj++ )
1293  {
1294  if ( bukrect.intersects( bucket_rects[ jj ] ) )
1295  { // Buckets overlap: test if overlap is virtually a line
1296  QRectF buki = bukrect.intersected( bucket_rects[ jj ] );
1297  double bwid = qAbs( buki.left() - buki.right() );
1298  double bhgt = qAbs( buki.bottom() - buki.top() );
1299 DbgLv(1) << "OVERLAP? ii jj" << ii << jj;
1300 DbgLv(1) << " buck i" << bukrect.left() << bukrect.right()
1301  << bukrect.top() << bukrect.bottom();
1302 QRectF brj=bucket_rects[jj];
1303 DbgLv(1) << " buck j" << brj.left() << brj.right()
1304  << brj.top() << brj.bottom();
1305 DbgLv(1) << " buki" << buki.left() << buki.right()
1306  << buki.top() << buki.bottom();
1307 DbgLv(1) << " buki left-right" << bwid << " bottom-top" << bhgt;
1308  if ( bwid < 1e-10 || bhgt < 1e-10 )
1309  { // Buckets only share a line: move a bit
1310  bucketSeparate( ii, jj, bucket_rects );
1311  }
1312 
1313  else
1314  { // Overlap is significant: mark it;
1315  novlps++;
1316  }
1317  }
1318  }
1319  }
1320 
1321 DbgLv(1) << " final NOVLPS" << novlps;
1322  return novlps;
1323 }
1324 
1325 // Modify buckets overlapping by virtually a line
1326 void US_SoluteData::bucketSeparate( int ii, int jj,
1327  QList< QRectF >& bucket_rects )
1328 {
1329  QRectF bukreci = bucket_rects[ ii ];
1330  QRectF bukrecj = bucket_rects[ jj ];
1331  double xdiff = qAbs( bukreci.right() - bukrecj.left() );
1332  double ydiff = qAbs( bukreci.bottom() - bukrecj.top() );
1333 
1334  if ( xdiff < 1e-10 )
1335  { // Overlap by just vertical line
1336  double xleft = bukrecj.left() + 1e-6;
1337  bukrecj.setLeft( xleft );
1338  }
1339 
1340  if ( ydiff < 1e-10 )
1341  { // Overlap by just horizontal line
1342  double ybottom = bukrecj.bottom() - 1e-6;
1343  bukrecj.setBottom( ybottom );
1344  }
1345 
1346  bucket bukj = bucketAt( jj );
1347  bukj.x_min = bukrecj.left();
1348  bukj.y_max = bukrecj.bottom();
1349  allbucks .replace( jj, bukj );
1350  bucket_rects.replace( jj, bukrecj );
1351 }
1352 
1353 // Return solute count of fullest bucket
1355 {
1356  int mxsol = 0;
1357  int nbuks = MC_solute.size();
1358 
1359  if ( nbuks == 0 )
1360  {
1361  buildDataMC();
1362  nbuks = MC_solute.size();
1363  }
1364 DbgLv(1) << "countFB-absize" << nbuks;
1365 
1366  for ( int ii = 0; ii < nbuks; ii++ )
1367  {
1368 DbgLv(1) << "countFB ii" << ii;
1369  mxsol = qMax( mxsol, MC_solute.at( ii ).size() );
1370 DbgLv(1) << "countFB mxsol" << mxsol;
1371  }
1372 
1373  return mxsol;
1374 }
1375