UltraScan III
us_pcsa_modelrec.cpp
Go to the documentation of this file.
1 
3 #include "us_pcsa_modelrec.h"
4 #include "us_solute.h"
5 #include "us_settings.h"
6 
7 // Construct model record
8 US_ModelRecord::US_ModelRecord( void )
9 {
10  taskx = -1;
11  stype = 11;
12  ctype = CTYPE_NONE;
13  v_ctype = CTYPE_NONE;
14  mrecID = 0;
15  editID = 0;
16  modelID = 0;
17  str_y = 0.0;
18  end_y = 0.0;
19  par1 = 0.0;
20  par2 = 0.0;
21  variance = 9999.9;
22  rmsd = 9999.9;
23  xmin = 0.0;
24  xmax = 0.0;
25  ymin = 0.0;
26  ymax = 0.0;
27  mrecGUID = "";
28  editGUID = "";
29  modelGUID = "";
30  isolutes.clear();
31  csolutes.clear();
32  clear_data();
33 }
34 
35 // Model record destructor
36 US_ModelRecord::~US_ModelRecord()
37 {
38  isolutes.clear();
39  csolutes.clear();
40 
41  clear_data();
42 }
43 
44 // Public slot to clear data vectors
45 void US_ModelRecord::clear_data( void )
46 {
47  sim_data .scanData.clear();
48  sim_data .xvalues .clear();
49  residuals.scanData.clear();
50  residuals.xvalues .clear();
51  ti_noise .clear();
52  ri_noise .clear();
53 }
54 
55 // Static public function to compute straight lines
56 int US_ModelRecord::compute_slines( double& xmin, double& xmax,
57  double& ymin, double& ymax, int& nypts, int& nlpts,
58  double* parlims, QVector< US_ModelRecord >& mrecs )
59 {
60  //int dbg_level = US_Settings::us_debug();
61  US_ModelRecord mrec;
62  int stype = parlims[ 4 ];
63  int attr_x = ( stype >> 6 ) & 7;
64  int attr_y = ( stype >> 3 ) & 7;
65  mrec.ctype = CTYPE_SL;
66  mrec.stype = stype;
67 
68  double prng = (double)( nlpts - 1 );
69  double yrng = (double)( nypts - 1 );
70  double xrng = xmax - xmin;
71  double xinc = xrng / prng;
72  if ( parlims[ 0 ] < 0.0 )
73  {
74  parlims[ 0 ] = ymin;
75  parlims[ 1 ] = ymax;
76  parlims[ 2 ] = ymin;
77  parlims[ 3 ] = ymax;
78  }
79 
80  double yslo = parlims[ 0 ];
81  double yshi = parlims[ 1 ];
82  double yelo = parlims[ 2 ];
83  double yehi = parlims[ 3 ];
84  double zval = parlims[ 5 ];
85  double ysinc = ( yshi - yslo ) / yrng;
86  double yeinc = ( yehi - yelo ) / yrng;
87  double ystr = yslo;
88  double xscl = ( attr_x == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
89  double yscl = ( attr_y == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
90  int mndx = mrecs.size();
91  int nmodl = nypts * nypts;
92 
93  // Generate straight lines
94  for ( int ii = 0; ii < nypts; ii++ )
95  { // Loop over k start values
96  double yend = yelo;
97 
98  for ( int jj = 0; jj < nypts; jj++ )
99  { // Loop over y end values
100  double xval = xmin;
101  double yval = ystr;
102  double yinc = ( yend - ystr ) / prng;
103 
104  mrec.isolutes.clear();
105  US_ZSolute isol( 0.0, 0.0, zval, 0.0 );
106 
107  for ( int kk = 0; kk < nlpts; kk++ )
108  { // Loop over points on a line
109  isol.x = xval * xscl;
110  isol.y = yval * yscl;
111  mrec.isolutes << isol;
112  xval += xinc;
113  yval += yinc;
114  } // END: points-per-line loop
115 
116  mrec.stype = stype;
117  mrec.taskx = mndx;
118  mrec.str_y = ystr;
119  mrec.end_y = yend;
120  mrec.par1 = ystr;
121  mrec.par2 = ( yend - ystr ) / xrng;
122 //DbgLv(1) << "MR: ii jj" << ii << jj << "ys ye p1 p2" << ystr << yend
123 // << mrec.par1 << mrec.par2;
124  mrecs << mrec;
125 
126  yend += yeinc;
127  mndx++;
128  } // END: k-end loop
129 
130  ystr += ysinc;
131  } // END: k-start loop
132 
133  return nmodl;
134 }
135 
136 // Static public function to compute sigmoid curves
137 int US_ModelRecord::compute_sigmoids( int& ctype, double& xmin, double& xmax,
138  double& ymin, double& ymax, int& nypts, int& nlpts,
139  double* parlims, QVector< US_ModelRecord >& mrecs )
140 {
141  if ( parlims[ 0 ] < 0.0 )
142  {
143  parlims[ 0 ] = 0.001;
144  parlims[ 1 ] = 0.5;
145  parlims[ 2 ] = 0.0;
146  parlims[ 3 ] = 1.0;
147  }
148  double p1lo = parlims[ 0 ];
149  double p1up = parlims[ 1 ];
150  double p2lo = parlims[ 2 ];
151  double p2up = parlims[ 3 ];
152  int stype = parlims[ 4 ];
153  double zval = parlims[ 5 ];
154 
155  US_ModelRecord mrec;
156  int attr_x = ( stype >> 6 ) & 7;
157  int attr_y = ( stype >> 3 ) & 7;
158  mrec.ctype = ctype;
159  mrec.stype = stype;
160 
161  double xrng = xmax - xmin;
162  double p1llg = log( p1lo );
163  double p1ulg = log( p1up );
164  double prng = (double)( nlpts - 1 );
165  double yrng = (double)( nypts - 1 );
166  double p1inc = ( p1ulg - p1llg ) / yrng;
167  double p2inc = ( p2up - p2lo ) / yrng;
168  double xinc = 1.0 / prng;
169  double xscl = ( attr_x == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
170  double yscl = ( attr_y == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
171  double ystr = ymin; // Start,Diff of 'IS'
172  double ydif = ymax - ymin;
173  if ( ctype == CTYPE_DS )
174  {
175  ystr = ymax; // Start,Diff of 'DS'
176  ydif = -ydif;
177  }
178  double p1vlg = p1llg;
179  int mndx = mrecs.size();
180  int nmodl = nypts * nypts;
181 
182  for ( int ii = 0; ii < nypts; ii++ )
183  { // Loop over par1 values (logarithmic progression)
184  double p1val = exp( p1vlg );
185  double p1fac = sqrt( 2.0 * p1val );
186  double p2val = p2lo;
187 
188  for ( int jj = 0; jj < nypts; jj++ )
189  { // Loop over par2 value (linear progression)
190  double xoff = 0.0;
191  double yval = ystr;
192  double yval0 = yval;
193  double xval = xmin;
194 
195  mrec.isolutes.clear();
196  US_ZSolute isol( 0.0, 0.0, zval, 0.0 );
197 
198  for ( int kk = 0; kk < nlpts; kk++ )
199  { // Loop over points on a curve
200  double efac = 0.5 * erf( ( xoff - p2val ) / p1fac ) + 0.5;
201  yval = ystr + ydif * efac;
202  xval = xmin + xoff * xrng;
203  xoff += xinc;
204  if ( kk == 0 )
205  yval0 = yval;
206 
207  isol.x = xval * xscl;
208  isol.y = yval * yscl;
209  mrec.isolutes << isol;
210  } // END: points-on-curve loop
211 
212  mrec.stype = stype;
213  mrec.taskx = mndx;
214  mrec.str_y = yval0;
215  mrec.end_y = yval;
216  mrec.par1 = p1val;
217  mrec.par2 = p2val;
218  mrecs << mrec;
219 
220  mndx++;
221  p2val += p2inc;
222  } // END: par2 values loop
223 
224  p1vlg += p1inc;
225  } // END: par1 values loop
226 
227  return nmodl;
228 }
229 
230 // Static public function to compute horizontal lines [ C(s) ]
231 int US_ModelRecord::compute_hlines( double& xmin, double& xmax,
232  double& ymin, double& ymax, int& nypts, int& nlpts,
233  double* parlims, QVector< US_ModelRecord >& mrecs )
234 {
235  US_ModelRecord mrec;
236  mrec.ctype = CTYPE_HL;
237 
238  double yrng = (double)( nypts - 1 );
239  double prng = (double)( nlpts - 1 );
240  double xinc = ( xmax - xmin ) / prng;
241  int stype = parlims[ 4 ];
242  mrec.stype = stype;
243 
244  if ( parlims[ 0 ] < 0.0 )
245  {
246  parlims[ 0 ] = ymin;
247  parlims[ 1 ] = ymax;
248  parlims[ 2 ] = ymin;
249  parlims[ 3 ] = ymax;
250  }
251 
252  int attr_x = ( stype >> 6 ) & 7;
253  int attr_y = ( stype >> 3 ) & 7;
254  double yval = parlims[ 0 ];
255  double xscl = ( attr_x == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
256  double yscl = ( attr_y == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
257  double yinc = ( parlims[ 1 ] - yval ) / yrng;
258  double zval = parlims[ 5 ];
259  int mndx = mrecs.size();
260 
261  // Generate horizontal lines
262  for ( int ii = 0; ii < nypts; ii++ )
263  { // Loop over k start and end values
264  double xval = xmin;
265 
266  mrec.isolutes.clear();
267  US_ZSolute isol( 0.0, 0.0, zval, 0.0 );
268 
269  for ( int kk = 0; kk < nlpts; kk++ )
270  { // Loop over points on a line
271  isol.x = xval * xscl;
272  isol.y = yval * yscl;
273  mrec.isolutes << isol;
274  xval += xinc;
275  } // END: points-per-line loop
276 
277  mrec.stype = stype;
278  mrec.taskx = mndx;
279  mrec.str_y = yval;
280  mrec.end_y = yval;
281  mrec.par1 = yval;
282  mrec.par2 = 0.0;
283  yval += yinc;
284  mndx++;
285 
286  mrecs << mrec;
287  } // END: k-value loop
288 
289  return nypts;
290 }
291 
292 // Static public function to compute Second-order Power Law curves
293 int US_ModelRecord::compute_2ndorder( double& xmin, double& xmax,
294  double& ymin, double& ymax, int& nypts, int& nlpts,
295  double* parlims, QVector< US_ModelRecord >& mrecs )
296 {
297  int dbg_level = US_Settings::us_debug();
298  //double p1lo = -1.0;
299  double p1lo = -2.0;
300  double p1up = -0.05;
301  int stype = parlims[ 4 ];
302  double zval = parlims[ 5 ];
303  int lyptx = nypts - 1;
304 
305  US_ModelRecord mrec;
306  int attr_x = ( stype >> 6 ) & 7;
307  int attr_y = ( stype >> 3 ) & 7;
308  mrec.stype = stype;
309  mrec.ctype = CTYPE_2O;
310 
311  double xrng = xmax - xmin;
312  double xterm = log10( 5.0 ) / log10( qMax( 1.00001, xrng * 0.5 ) );
313  double lprng = (double)( nlpts - 1 );
314  double yprng = (double)( lyptx );
315  double xinc = xrng / lprng;
316  double xscl = ( attr_x == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
317  double yscl = ( attr_y == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
318  double yrng = ymax - ymin;
319  double yinc = yrng / yprng;
320 
321  double bmin = p1lo * xterm;
322  double bmax = p1up * xterm;
323  double amax = pow( xmin, bmax ) * ( yrng - yinc );
324  double amin = amax / yprng;
325  double cmin = ymin;
326  double cmax = ymin + yinc;
327  double ainc = ( amax - amin ) / yprng;
328  double binc = ( bmax - bmin ) / yprng;
329  double cinc = ( cmax - cmin ) / yprng;
330 //amin = xmin^-max(bvec)*(((ymax - ymin) - (ymax - ymin) / var) / var); % lower limit of variable a
331 //amax = xmin^-max(bvec)*(ymax - ymin - (ymax - ymin) / var); % upper limit of variable a
332 DbgLv(1) << "MR: p1l p1u" << p1lo << p1up;
333 
334  if ( parlims[ 0 ] > 0.0 )
335  { // Recomputing: Adjust A,B,C ranges so previous best in on the grid
336  amin = parlims[ 0 ];
337  amax = parlims[ 1 ];
338  bmin = parlims[ 2 ];
339  bmax = parlims[ 3 ];
340  cmin = parlims[ 6 ];
341  cmax = parlims[ 7 ];
342  double abest = parlims[ 8 ];
343  double bbest = parlims[ 9 ];
344  double cbest = parlims[ 10 ];
345  amin = qMin( abest, amin );
346  cmin = qMin( cbest, cmin );
347  ainc = ( amax - amin ) / yprng;
348  binc = ( bmax - bmin ) / yprng;
349  cinc = ( cmax - cmin ) / yprng;
350 
351  double vbx = ( abest - amin ) / ainc; // Adjust so A-best on grid
352  double vbi = qFloor( vbx );
353  double vbf = vbx - vbi;
354 DbgLv(1) << "MR: NEW abest vbx vbi vbf" << abest << vbx << vbi << vbf;
355  if ( vbi == 0.0 )
356  amin = abest - ainc;
357  else if ( vbf < 0.5 )
358  amin = abest - ainc * vbi;
359  else
360  amin = amin - ainc * vbf;
361 
362  vbx = ( bbest - bmin ) / binc; // Adjust so B-best on grid
363  vbi = qFloor( vbx );
364  vbf = vbx - vbi;
365 DbgLv(1) << "MR: NEW bbest vbx vbi vbf" << bbest << vbx << vbi << vbf;
366 DbgLv(1) << "MR: OLD bmin bmax binc" << bmin << bmax << binc;
367  if ( vbi == 0.0 )
368  bmin = bbest - binc;
369  else if ( vbf < 0.5 )
370  bmin = bbest - binc * vbi;
371  else
372  bmin = bmin - binc * vbf;
373 
374  vbx = ( cbest - cmin ) / cinc; // Adjust so C-best on grid
375  vbi = qFloor( vbx );
376  vbf = vbx - vbi;
377 DbgLv(1) << "MR: NEW cbest vbx vbi vbf" << cbest << vbx << vbi << vbf;
378  if ( vbi == 0.0 )
379  cmin = cbest - cinc;
380  else if ( vbf < 0.5 )
381  cmin = cbest - cinc * vbi;
382  else
383  cmin = cmin - cinc * vbf;
384 
385  amax = amin + ainc * yprng;
386  bmax = bmin + binc * yprng;
387  cmax = cmin + cinc * yprng;
388 DbgLv(1) << "MR: NEW amin amax ainc" << amin << amax << ainc;
389 DbgLv(1) << "MR: NEW bmin bmax binc" << bmin << bmax << binc;
390 DbgLv(1) << "MR: NEW cmin cmax cinc" << cmin << cmax << cinc;
391  }
392 
393  QVector< double > avec;
394  QVector< double > bvec;
395  QVector< double > cvec;
396  avec.resize( nypts );
397  bvec.resize( nypts );
398  cvec.resize( nypts );
399 
400  for ( int jj = 0; jj < nypts; jj++ )
401  {
402  double poff = (double)jj;
403  avec[ jj ] = amin + ainc * poff;
404  bvec[ jj ] = bmin + binc * poff;
405  cvec[ jj ] = cmin + cinc * poff;
406  }
407 
408 DbgLv(1) << "MR: NEWER bmin bmax" << bmin << bmax;
409  double* avals = avec.data();
410  double* bvals = bvec.data();
411  double* cvals = cvec.data();
412  parlims[ 0 ] = avals[ 0 ];
413  parlims[ 1 ] = avals[ lyptx ];
414  parlims[ 2 ] = bvals[ 0 ];
415  parlims[ 3 ] = bvals[ lyptx ];
416  parlims[ 6 ] = cvals[ 0 ];
417  parlims[ 7 ] = cvals[ lyptx ];
418 /***
419  * MathLab version of Second-order Power Law variations
420  *
421 xmin = 1; % lower limit of x
422 xmax = 300; % upper limit of x
423 resolution = 1000; % resolution of x
424 
425 x = linspace(xmin,xmax,resolution); % vector of x
426 
427 ymin = 0.1; % lower limit of y
428 ymax = 0.5; % upper limit of y
429 
430 var = 5; % variations count of the three variables
431 
432 bmin = -1; % lower reference limit of variable b
433 bmax = -0.05; % upper reference limit of variable b
434 bvec = linspace(bmin,bmax,var); % reference vector of variable b
435 
436 for m = 1 : var % transfer of reference vector b to any x-range
437  bvec(m) = bvec(m) * log10(5)/log10((xmax-xmin)/2); % (xmax - xmin)/2 has to be larger than 1!
438 end
439 
440 %% remark: bmax is required for the calculation of avec
441 
442 amin = xmin^-max(bvec)*(((ymax - ymin) - (ymax - ymin) / var) / var); % lower limit of variable a
443 amax = xmin^-max(bvec)*(ymax - ymin - (ymax - ymin) / var); % upper limit of variable a
444 avec = linspace(amin,amax,var); % vector of variable a
445 
446 cmin = ymin; % lower limit of variable c
447 cmax = ymin+((ymax - ymin) / var); % upper limit of variable c
448 cvec = linspace(cmin,cmax,var); % vector of variable c
449 
450 %% create figure
451 figure(1)
452 clf;
453 xlabel('x')
454 ylabel('y');
455 
456 %% perform parametrization
457 for m = 1 : var % amin to amax:
458  for n = 1 : var % bmin to bmax:
459  for k = 1 : var % cmin to cmax:
460  for p = 1 : size(x) % xmin to xmax with resolution:
461  y = avec(m)*(x.^bvec(n))+cvec(k);
462  %% plot data
463  axis([0 max(x) 0 1])
464  hold on
465  plot(x,y)
466  hold off
467  end
468  end
469  end
470 end
471 
472 % reference plot
473 y = 0.2994 * x.^-0.4685 + 0.1765;
474 hold on
475 plot(x,y,'r--o')
476 hold off
477  ***/
478  int mndx = mrecs.size();
479  int nmodl = nypts * nypts * nypts;
480  double yvmin = 1e+99;
481  double yvmax = -1e+99;
482 
483  for ( int ii = 0; ii < nypts; ii++ )
484  { // Loop over a values
485  double aval = avals[ ii ];
486 
487  for ( int jj = 0; jj < nypts; jj++ )
488  { // Loop over b values
489  double bval = bvals[ jj ];
490 
491  for ( int kk = 0; kk < nypts; kk++ )
492  { // Loop over c values
493  double cval = cvals[ kk ];
494  double xval = xmin;
495  double yval = ymin;
496  mrec.str_y = aval * pow( xval, bval ) + cval;
497 
498  mrec.isolutes.clear();
499  US_ZSolute isol( 0.0, 0.0, zval, 0.0 );
500 DbgLv(1) << "MR: a b c" << aval << bval << cval
501  << " ii jj kk" << ii << jj << kk;
502 
503  for ( int ll = 0; ll < nlpts; ll++ )
504  { // Loop over points on a curve
505  yval = aval * pow( xval, bval ) + cval;
506 if(ll<4 || (ll+4)>nlpts)
507 DbgLv(1) << "MR: ll" << ll << "x y" << xval << yval;
508 
509  isol.x = xval * xscl;
510  isol.y = yval * yscl;
511  mrec.isolutes << isol;
512 
513  xval += xinc;
514  } // END: points-on-curve loop
515 
516  mrec.stype = stype;
517  mrec.taskx = mndx;
518  mrec.end_y = yval;
519  mrec.par1 = aval;
520  mrec.par2 = bval;
521  mrec.par3 = cval;
522  mrecs << mrec;
523 
524  yvmin = qMin( yvmin, qMin( mrec.str_y, mrec.end_y ) );
525  yvmax = qMax( yvmax, qMax( mrec.str_y, mrec.end_y ) );
526 
527  mndx++;
528  } // END: C values loop
529  } // END: B values loop
530  } // END: A values loop
531 DbgLv(1) << "MR: yvmin yvmax" << yvmin << yvmax;
532 
533  return nmodl;
534 }
535 
536 // Static public function to load model records from an XML stream
537 int US_ModelRecord::load_modelrecs( QXmlStreamReader& xml,
538  QVector< US_ModelRecord >& mrecs, QString& descr, int& ctype,
539  double& xmin, double& xmax, double& ymin, double& ymax, int& stype )
540 {
541 // int dbg_level = US_Settings::us_debug();
542  int nmrecs = 0;
543  int nisols = 0;
544  int ncsols = 0;
545  int attr_x = ( stype >> 6 ) & 7;
546  int attr_y = ( stype >> 3 ) & 7;
547  double xscl = ( attr_x == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
548  double yscl = ( attr_y == US_ZSolute::ATTR_S ) ? 1.0e-13 : 1.0;
549  bool old_vers = false;
550  const double ov = 1.0;
551  US_ModelRecord mrec;
552  QStringList sctypes;
553  sctypes << "SL" << "IS" << "DS" << "All" << "HL" << "2O";
554  mrecs.clear();
555 
556  while ( ! xml.atEnd() )
557  {
558  xml.readNext();
559  QString xmlname = xml.name().toString();
560 
561  if ( xml.isStartElement() )
562  {
563  QXmlStreamAttributes attrs = xml.attributes();
564 
565  if ( xmlname == "modelrecords" )
566  {
567  double vers = attrs.value( "version" ).toString().toDouble();
568  old_vers = ( vers <= ov );
569  QString s_type1 = attrs.value( "type" ).toString();
570  QString s_type = attrs.value( "curve_type" ).toString();
571  s_type = s_type.isEmpty() ? s_type1 : s_type;
572  ctype = sctypes.contains( s_type ) ?
573  ctype_flag( s_type ) : s_type.toInt();
574  if ( old_vers )
575  {
576  ctype = ( ctype == 2 ) ? CTYPE_DS : ctype;
577  ctype = ( ctype == 1 ) ? CTYPE_IS : ctype;
578  ctype = ( ctype == 0 ) ? CTYPE_SL : ctype;
579  }
580 //DbgLv(1) << "STYPE" << s_type << "CTYPE" << ctype;
581  QString xlo = attrs.value( "xmin" ).toString();
582  QString xhi = attrs.value( "xmax" ).toString();
583  QString ylo = attrs.value( "ymin" ).toString();
584  QString yhi = attrs.value( "ymax" ).toString();
585  QString slo = attrs.value( "smin" ).toString();
586  QString shi = attrs.value( "smax" ).toString();
587  QString klo = attrs.value( "kmin" ).toString();
588  QString khi = attrs.value( "kmax" ).toString();
589  QString s_stype = attrs.value( "solute_type" ).toString();
590  stype = stype_flag( s_stype );
591  xlo = xlo.isEmpty() ? slo : xlo;
592  xhi = xhi.isEmpty() ? shi : xhi;
593  ylo = ylo.isEmpty() ? klo : ylo;
594  yhi = yhi.isEmpty() ? khi : yhi;
595  xmin = xlo.toDouble();
596  xmax = xhi.toDouble();
597  ymin = ylo.toDouble();
598  ymax = yhi.toDouble();
599  nisols = attrs.value( "curve_points" ).toString().toInt();
600  mrec.v_ctype = ctype;
601  mrec.ctype = ctype;
602  mrec.xmin = xmin;
603  mrec.xmax = xmax;
604  mrec.ymin = ymin;
605  mrec.ymax = ymax;
606  mrec.stype = stype;
607  descr = attrs.value( "description" ).toString();
608  QString s_eID = attrs.value( "editID" ).toString();
609  QString s_mID = attrs.value( "modelID" ).toString();
610  mrec.mrecGUID = attrs.value( "mrecGUID" ).toString();
611  mrec.editGUID = attrs.value( "editGUID" ).toString();
612  mrec.modelGUID = attrs.value( "modelGUID" ).toString();
613  mrec.editID = s_eID.isEmpty() ? 0 : s_eID.toInt();
614  mrec.modelID = s_mID.isEmpty() ? 0 : s_mID.toInt();
615  }
616 
617  else if ( xmlname == "modelrecord" )
618  {
619  QString s_ctype = attrs.value( "type" ).toString();
620  mrec.ctype = sctypes.contains( s_ctype ) ?
621  ctype_flag( s_ctype ) : s_ctype.toInt();
622  if ( old_vers )
623  {
624  mrec.ctype = ( mrec.ctype == 2 ) ? CTYPE_DS : mrec.ctype;
625  mrec.ctype = ( mrec.ctype == 1 ) ? CTYPE_IS : mrec.ctype;
626  mrec.ctype = ( mrec.ctype == 0 ) ? CTYPE_SL : mrec.ctype;
627  }
628  mrec.taskx = attrs.value( "taskx" ).toString().toInt();
629  QString stry = attrs.value( "start_y" ).toString();
630  QString endy = attrs.value( "end_y" ).toString();
631  QString strk = attrs.value( "start_k" ).toString();
632  QString endk = attrs.value( "end_k" ).toString();
633  mrec.str_y = stry.isEmpty() ? strk.toDouble()
634  : stry.toDouble();
635  mrec.end_y = endy.isEmpty() ? endk.toDouble()
636  : endy.toDouble();
637  mrec.par1 = attrs.value( "par1" ).toString().toDouble();
638  mrec.par2 = attrs.value( "par2" ).toString().toDouble();
639  mrec.par3 = attrs.value( "par3" ).toString().toDouble();
640  mrec.rmsd = attrs.value( "rmsd" ).toString().toDouble();
641  QString s_cvpt = attrs.value( "curve_points" ).toString();
642  int kisols = s_cvpt.isEmpty() ? nisols : s_cvpt.toInt();
643 
644  if ( nmrecs > 0 )
645  {
646  mrec.mrecID = 0;
647  mrec.mrecGUID = "";
648  mrec.editID = 0;
649  mrec.editGUID = "";
650  QString s_mID = attrs.value( "modelID" ).toString();
651  mrec.modelID = s_mID.isEmpty() ? 0 : s_mID.toInt();
652  mrec.modelGUID = attrs.value( "modelGUID" ).toString();
653  }
654 
655  ncsols = 0;
656  mrec.isolutes.resize( kisols );
657  mrec.csolutes.clear();
658  ymin = qMin( ymin, mrec.str_y );
659  ymax = qMax( ymax, mrec.end_y );
660  }
661 
662  else if ( xmlname == "c_solute" )
663  {
664  US_ZSolute csolute;
665  QString xval = attrs.value( "x" ).toString();
666  QString sval = attrs.value( "s" ).toString();
667  QString yval = attrs.value( "y" ).toString();
668  QString kval = attrs.value( "k" ).toString();
669  csolute.x = sval.isEmpty() ? xval.toDouble()
670  : sval.toDouble();
671  csolute.y = kval.isEmpty() ? yval.toDouble()
672  : kval.toDouble();
673  csolute.z = attrs.value( "z" ).toString().toDouble();
674  csolute.c = attrs.value( "c" ).toString().toDouble();
675  xmin = qMin( xmin, csolute.x );
676  xmax = qMax( xmax, csolute.x );
677  ymin = qMin( ymin, csolute.y );
678  ymax = qMax( ymax, csolute.y );
679  csolute.x *= xscl;
680  csolute.y *= yscl;
681 
682  mrec.csolutes << csolute;
683  ncsols++;
684  }
685  }
686 
687  else if ( xml.isEndElement() && xmlname == "modelrecord" )
688  {
689  nmrecs++;
690  mrecs << mrec;
691  mrec.csolutes.clear();
692  }
693  }
694  return nmrecs;
695 }
696 
697 // Static public function to write model records to an XML stream
698 int US_ModelRecord::write_modelrecs( QXmlStreamWriter& xml,
699  QVector< US_ModelRecord >& mrecs, QString& descr, int& ctype,
700  double& xmin, double& xmax, double& ymin, double& ymax, int& stype )
701 {
702 // int dbg_level = US_Settings::us_debug();
703  int nmrecs = mrecs.size();
704  US_ModelRecord mrec;
705  mrec = mrecs[ 0 ];
706  ctype = ( ctype == CTYPE_NONE ) ? mrec.ctype : ctype;
707  int v_ctype = ctype;
708  xmin = mrec.xmin;
709  xmax = mrec.xmax;
710  ymin = mrec.ymin;
711  ymax = mrec.ymax;
712  int nisols = mrec.isolutes.size();
713  int attr_x = ( stype >> 6 ) & 7;
714  int attr_y = ( stype >> 3 ) & 7;
715  double xscl = ( attr_x == US_ZSolute::ATTR_S ) ? 1.0e+13 : 1.0;
716  double yscl = ( attr_y == US_ZSolute::ATTR_S ) ? 1.0e+13 : 1.0;
717  xml.setAutoFormatting( true );
718  xml.writeStartDocument( "1.0" );
719  xml.writeComment( "DOCTYPE PcsaModelRecords" );
720  xml.writeCharacters( "\n" );
721  xml.writeStartElement( "modelrecords" );
722  xml.writeAttribute( "version", "1.2" );
723  xml.writeAttribute( "description", descr );
724  xml.writeAttribute( "curve_type", ctype_text( v_ctype ) );
725  xml.writeAttribute( "xmin", QString::number( xmin ) );
726  xml.writeAttribute( "xmax", QString::number( xmax ) );
727  xml.writeAttribute( "ymin", QString::number( ymin ) );
728  xml.writeAttribute( "ymax", QString::number( ymax ) );
729  xml.writeAttribute( "curve_points", QString::number( nisols ) );
730  xml.writeAttribute( "solute_type", stype_text( stype ) );
731  if ( mrec.editID > 0 )
732  xml.writeAttribute( "editID", QString::number( mrec.editID ) );
733  if ( mrec.modelID > 0 )
734  xml.writeAttribute( "modelID", QString::number( mrec.modelID ) );
735  if ( ! mrec.mrecGUID.isEmpty() )
736  xml.writeAttribute( "mrecGUID", mrec.mrecGUID );
737  if ( ! mrec.editGUID.isEmpty() )
738  xml.writeAttribute( "editGUID", mrec.editGUID );
739  if ( ! mrec.modelGUID.isEmpty() )
740  xml.writeAttribute( "modelGUID", mrec.modelGUID );
741 
742  for ( int mr = 0; mr < nmrecs; mr++ )
743  {
744  mrec = mrecs[ mr ];
745  int kisols = mrec.isolutes.size();
746  int ncsols = mrec.csolutes.size();
747  xml.writeStartElement( "modelrecord" );
748  xml.writeAttribute( "taskx", QString::number( mrec.taskx ) );
749  xml.writeAttribute( "type", ctype_text( mrec.ctype ) );
750  xml.writeAttribute( "start_y", QString::number( mrec.str_y ) );
751  xml.writeAttribute( "end_y", QString::number( mrec.end_y ) );
752  xml.writeAttribute( "par1", QString::number( mrec.par1 ) );
753  xml.writeAttribute( "par2", QString::number( mrec.par2 ) );
754  xml.writeAttribute( "par3", QString::number( mrec.par3 ) );
755  xml.writeAttribute( "rmsd", QString::number( mrec.rmsd ) );
756  if ( kisols != nisols )
757  xml.writeAttribute( "curve_points", QString::number( kisols ) );
758 
759  if ( mr == 0 )
760  {
761  if ( mrec.editID > 0 )
762  xml.writeAttribute( "editID", QString::number( mrec.editID ) );
763  if ( ! mrec.mrecGUID.isEmpty() )
764  xml.writeAttribute( "mrecGUID", mrec.mrecGUID );
765  if ( ! mrec.editGUID.isEmpty() )
766  xml.writeAttribute( "editGUID", mrec.editGUID );
767  }
768 
769  if ( mrec.modelID > 0 )
770  xml.writeAttribute( "modelID", QString::number( mrec.modelID ) );
771  if ( ! mrec.modelGUID.isEmpty() )
772  xml.writeAttribute( "modelGUID", mrec.modelGUID );
773 
774  for ( int cc = 0; cc < ncsols; cc++ )
775  {
776  xml.writeStartElement( "c_solute" );
777  double xval = mrec.csolutes[ cc ].x * xscl;
778  double yval = mrec.csolutes[ cc ].y * yscl;
779  xml.writeAttribute( "x", QString::number( xval ) );
780  xml.writeAttribute( "y", QString::number( yval ) );
781  xml.writeAttribute( "z", QString::number( mrec.csolutes[ cc ].z ) );
782  xml.writeAttribute( "c", QString::number( mrec.csolutes[ cc ].c ) );
783  xml.writeEndElement();
784  }
785 
786  xml.writeEndElement();
787  }
788 
789  xml.writeEndElement();
790  xml.writeEndDocument();
791 
792  return nmrecs;
793 }
794 
795 
796 // Static public function to determine model records elite (top 10%) limits
797 void US_ModelRecord::elite_limits( QVector< US_ModelRecord >& mrecs,
798  int& ctype, double& minyv, double& maxyv,
799  double& minp1, double& maxp1, double& minp2, double& maxp2,
800  double& minp3, double& maxp3 )
801 {
802  int dbg_level = US_Settings::us_debug();
803  double efrac = 0.1;
804  // Set up variables that help insure that the par1,par2 extents of elites
805  // extend at least one step on either side of record 0's par1,par2
806  double m0p1 = mrecs[ 0 ].par1;
807  double m0p2 = mrecs[ 0 ].par2;
808  double m0p1l = m0p1;
809  double m0p1h = m0p1;
810  double m0p2l = m0p2;
811  double m0p2h = m0p2;
812  bool ln_type = false;
813 
814  if ( ctype == CTYPE_SL || ctype == CTYPE_HL )
815  { // Possibly adjust initial par1,par2 limits for lines
816  ln_type = true;
817  m0p1 = mrecs[ 0 ].str_y;
818  m0p2 = mrecs[ 0 ].end_y;
819  m0p1l = ( m0p1 > maxyv ) ? m0p1 : ( m0p1 * 1.0001 );
820  m0p1h = ( m0p1 < minyv ) ? m0p1 : ( m0p1 * 0.9991 );
821  m0p2l = ( m0p2 > maxyv ) ? m0p2 : ( m0p2 * 1.0001 );
822  m0p2h = ( m0p2 < minyv ) ? m0p2 : ( m0p2 * 0.9991 );
823  }
824 
825  if ( ctype == CTYPE_IS || ctype == CTYPE_DS )
826  { // Possibly adjust initial par1,par2 limits for sigmoids
827  double dif1 = m0p1 - 0.001;
828  double dif2 = m0p1 - 0.500;
829  double dif3 = m0p2 - 0.000;
830  double dif4 = m0p2 - 1.000;
831  m0p1l = ( dif1 > 1.e-8 ) ? m0p1 : 0.002;
832  m0p1h = ( dif2 < -1e-8 ) ? m0p1 : 0.499;
833  m0p2l = ( dif3 > 1.e-8 ) ? m0p2 : 0.001;
834  m0p2h = ( dif4 < -1e-8 ) ? m0p2 : 0.999;
835 //DbgLv(1) << " ElLim: ADJUST SIGM: m0p1 m0p1h" << m0p1 << m0p1h
836 // << "m0p1<0.500" << (m0p1<0.500) << 0.500 << "(m0p1-0.5)" << (m0p1-0.5);
837  }
838 
839  if ( ctype == CTYPE_2O )
840  { // Set 2nd-order par1,par2 limits so scan stops at elite cutoff point
841  m0p1l = 1e+99;
842  m0p1h = -1e+99;
843  m0p2l = 1e+99;
844  m0p2h = -1e+99;
845  efrac = mrecs[ 0 ].variance;
846  }
847 
848  int nmrec = mrecs.size();
849  int nelite = qRound( efrac * nmrec ); // Elite is top 10%
850  int maxel = nmrec / 2;
851  int minel = qMin( maxel, 4 );
852  nelite = qMin( nelite, maxel ); // At most half of all
853  nelite = qMax( nelite, minel ); // At least 4
854  nelite -= 2; // Less 2 for compare
855 //DbgLv(1) << " ElLim: nmrec nelite" << nmrec << nelite;
856 //DbgLv(1) << " ElLim: in minyv maxyv" << minyv << maxyv;
857 //DbgLv(1) << " ElLim: in min/max p1/p2" << minp1 << maxp1 << minp2 << maxp2;
858 //DbgLv(1) << " ElLim: in m0p1 m0p2" << m0p1 << m0p2;
859 //DbgLv(1) << " ElLim: in m0p1l,m0p1h,m0p2l,m0p2h" << m0p1l << m0p1h
860 // << m0p2l << m0p2h;
861 
862  for ( int ii = 0; ii < nmrec; ii++ )
863  {
864  double str_y = mrecs[ ii ].str_y;
865  double end_y = mrecs[ ii ].end_y;
866  double par1 = ln_type ? str_y : mrecs[ ii ].par1;
867  double par2 = ln_type ? end_y : mrecs[ ii ].par2;
868  double par3 = ln_type ? 0.0 : mrecs[ ii ].par3;
869 if(ii<3||(ii+4)>nelite)
870 DbgLv(1) << " ElLim: ii" << ii << "par1 par2 par3" << par1 << par2 << par3
871  << "str_y end_y" << str_y << end_y << "rmsd" << mrecs[ii].rmsd;
872  minyv = qMin( minyv, qMin( str_y, end_y ) );
873  maxyv = qMax( maxyv, qMax( str_y, end_y ) );
874  minp1 = qMin( minp1, par1 );
875  maxp1 = qMax( maxp1, par1 );
876  minp2 = qMin( minp2, par2 );
877  maxp2 = qMax( maxp2, par2 );
878  minp3 = qMin( minp3, par3 );
879  maxp3 = qMax( maxp3, par3 );
880 
881  // We want to break out of the min,max scan loop if the sorted index
882  // exceeds the elite count. But we continue in the loop if we have not
883  // yet found min,max par1,par2 values that are at least a step
884  // on either side of the par1,par2 values for the best model (m0).
885 if(ii>nelite)
886 DbgLv(1) << " ElLim: minp1 maxp1 m0p1l m0p1h" << minp1 << maxp1 << m0p1l
887  << m0p1h << "minp2 maxp2 m0p2l m0p2h" << minp2 << maxp2 << m0p2l << m0p2h
888  << "minyv maxyv" << minyv << maxyv;
889  if ( ii > nelite &&
890  minp1 < m0p1l && maxp1 > m0p1h &&
891  minp2 < m0p2l && maxp2 > m0p2h )
892  break;
893  }
894 
895 //DbgLv(1) << " ElLim: out minyv maxyv" << minyv << maxyv;
896 //DbgLv(1) << " ElLim: out min/max p1/p2" << minp1 << maxp1 << minp2 << maxp2;
897 }
898 
899 // Static public function to recompute model records vector for new iteration
900 int US_ModelRecord::recompute_mrecs( int& ctype, double& xmin, double& xmax,
901  double& ymin, double& ymax, int& nypts, int& nlpts, double* parlims,
902  QVector< US_ModelRecord >&mrecs )
903 {
904  int dbg_level = US_Settings::us_debug();
905  int nmrec = mrecs.size();
906  bool LnType = ( ctype == CTYPE_SL || ctype == CTYPE_HL );
907  bool SgType = ( ctype == CTYPE_IS || ctype == CTYPE_DS );
908 
909  double minyv = ymax;
910  double maxyv = ymin;
911  double minp1 = LnType ? minyv : ( SgType ? 0.500 : 1e+99 );
912  double maxp1 = LnType ? maxyv : ( SgType ? 0.001 : -1e+99 );
913  double minp2 = LnType ? minyv : ( SgType ? 1.0 : 1e+99 );
914  double maxp2 = LnType ? maxyv : ( SgType ? 0.0 : -1e+99 );
915  double minp3 = 1e+99;
916  double maxp3 = -1e+99;
917 //DbgLv(1) << "RF: 2)nmrec" << mrecs.size();
918 
919  elite_limits( mrecs, ctype, minyv, maxyv,
920  minp1, maxp1, minp2, maxp2, minp3, maxp3 );
921 
922  // Recompute the new min,max so that the old best is a point
923  // on the new grid to be tested
924  double p1best = LnType ? mrecs[ 0 ].str_y : mrecs[ 0 ].par1;
925  double p2best = LnType ? mrecs[ 0 ].end_y : mrecs[ 0 ].par2;
926  double p3best = mrecs[ 0 ].par3;
927  double yprng = (double)( nypts - 1 );
928  double p1rng = maxp1 - minp1;
929  double p2rng = maxp2 - minp2;
930  double p1inc = p1rng / yprng;
931  double p2inc = p2rng / yprng;
932  double p1dif = qRound( ( p1best - minp1 ) / p1inc ) * p1inc;
933  double p2dif = qRound( ( p2best - minp2 ) / p2inc ) * p2inc;
934 //DbgLv(1) << "RF: rcomp: p12 rng" << p1rng << p2rng << "p12 inc"
935 // << p1inc << p2inc << "p12 dif" << p1dif << p2dif;
936 //DbgLv(1) << "RF: rcomp: mmp1 mmp2" << minp1 << maxp1 << minp2 << maxp2;
937  if ( ctype != CTYPE_2O )
938  { // Adjust par1,par2 limits
939  minp1 = p1best - p1dif;
940  minp2 = p2best - p2dif;
941  maxp1 = minp1 + p1inc * yprng;
942  maxp2 = minp2 + p2inc * yprng;
943  }
944 //DbgLv(1) << "RF: rcomp: nypts" << nypts;
945 
946  mrecs .clear();
947 
948  if ( LnType )
949  { // Determine models for straight-line or horizontal-line curves
950  parlims[ 0 ] = qMax( ymin, minp1 );
951  parlims[ 1 ] = qMin( ymax, maxp1 );
952  parlims[ 2 ] = qMax( ymin, minp2 );
953  parlims[ 3 ] = qMin( ymax, maxp2 );
954  parlims[ 6 ] = 0.0;
955  parlims[ 7 ] = 0.0;
956  parlims[ 8 ] = p1best;
957  parlims[ 9 ] = p2best;
958  parlims[ 10 ] = 0.0;
959 //DbgLv(1) << "RF: slin: mmp1 mmp2" << minp1 << maxp1 << minp2 << maxp2;
960 //DbgLv(1) << "RF: slin: p1,p2 best" << p1best << p2best;
961 //DbgLv(1) << "RF: slin: mnmx p1 p2" << parlims[0] << parlims[1]
962 // << parlims[2] << parlims[3] << "kmax" << kmax;
963 
964  if ( ctype == CTYPE_SL )
965  {
966  nmrec = compute_slines( xmin, xmax, ymin, ymax,
967  nypts, nlpts, parlims, mrecs );
968  }
969 
970  else if ( ctype == CTYPE_HL )
971  {
972  nmrec = compute_hlines( xmin, xmax, ymin, ymax,
973  nypts, nlpts, parlims, mrecs );
974  }
975  }
976 
977  else if ( SgType )
978  { // Determine models for sigmoid curves
979 //DbgLv(1) << "RF: sigm: mnmx p1 p2" << minp1 << maxp1 << minp2 << maxp2;
980  parlims[ 0 ] = qMax( 0.001, minp1 );
981  parlims[ 1 ] = qMin( 0.500, maxp1 );
982  parlims[ 2 ] = qMax( 0.000, minp2 );
983  parlims[ 3 ] = qMin( 1.000, maxp2 );
984  parlims[ 6 ] = 0.0;
985  parlims[ 7 ] = 0.0;
986  parlims[ 8 ] = p1best;
987  parlims[ 9 ] = p2best;
988  parlims[ 10 ] = 0.0;
989 //DbgLv(1) << "RF: sigm: p1,p2 best" << p1best << p2best;
990 
991  nmrec = compute_sigmoids( ctype, xmin, xmax, ymin, ymax,
992  nypts, nlpts, parlims, mrecs );
993 //DbgLv(1) << "RF: sigm: nmrec" << nmrec;
994  }
995 
996  else if ( ctype == CTYPE_2O )
997  { // Determine models for 2nd-order Power Law curves
998  if ( minp1 == maxp1 )
999  {
1000  minp1 *= 0.99;
1001  maxp1 *= 1.01;
1002  }
1003  if ( minp2 == maxp2 )
1004  {
1005  minp2 *= 1.01;
1006  maxp2 *= 0.99;
1007  }
1008  if ( minp3 == maxp3 )
1009  {
1010  minp3 *= 0.99;
1011  maxp3 *= 1.01;
1012  }
1013  parlims[ 0 ] = minp1;
1014  parlims[ 1 ] = maxp1;
1015  parlims[ 2 ] = minp2;
1016  parlims[ 3 ] = maxp2;
1017  parlims[ 6 ] = minp3;
1018  parlims[ 7 ] = maxp3;
1019  parlims[ 8 ] = p1best;
1020  parlims[ 9 ] = p2best;
1021  parlims[ 10 ] = p3best;
1022  ymin = minyv;
1023  ymax = maxyv;
1024 DbgLv(1) << "RF: 2ord: nmrec" << nmrec << "ymin ymax" << ymin << ymax
1025  << "minp1 maxp1" << minp1 << maxp1 << "minp2 maxp2" << minp2 << maxp2
1026  << "parlims0" << parlims[0];
1027 
1028  nmrec = compute_2ndorder( xmin, xmax, ymin, ymax,
1029  nypts, nlpts, parlims, mrecs );
1030 DbgLv(1) << "RF: 2ord: parlims1-4" << parlims[0] << parlims[1]
1031  << parlims[2] << parlims[3];
1032  }
1033 
1034  else if ( ctype == CTYPE_ALL )
1035  {
1036  }
1037 
1038  else
1039  nmrec = 0;
1040 
1041  return nmrec;
1042 }
1043 
1044 // Static public function to return integer curve-type flag for given string
1045 int US_ModelRecord::ctype_flag( const QString s_ctype )
1046 {
1047  int i_ctype = CTYPE_NONE;
1048  i_ctype = ( s_ctype == "SL" ) ? CTYPE_SL : i_ctype;
1049  i_ctype = ( s_ctype == "IS" ) ? CTYPE_IS : i_ctype;
1050  i_ctype = ( s_ctype == "DS" ) ? CTYPE_DS : i_ctype;
1051  i_ctype = ( s_ctype == "HL" ) ? CTYPE_HL : i_ctype;
1052  i_ctype = ( s_ctype == "2O" ) ? CTYPE_2O : i_ctype;
1053  i_ctype = ( s_ctype == "All" ) ? CTYPE_ALL : i_ctype;
1054 
1055  return i_ctype;
1056 }
1057 
1058 // Static public function to return curve-type text for given integer flag
1059 QString US_ModelRecord::ctype_text( const int i_ctype )
1060 {
1061  QString s_ctype = "None";
1062  s_ctype = ( i_ctype == CTYPE_SL ) ? "SL" : s_ctype;
1063  s_ctype = ( i_ctype == CTYPE_IS ) ? "IS" : s_ctype;
1064  s_ctype = ( i_ctype == CTYPE_DS ) ? "DS" : s_ctype;
1065  s_ctype = ( i_ctype == CTYPE_HL ) ? "HL" : s_ctype;
1066  s_ctype = ( i_ctype == CTYPE_2O ) ? "2O" : s_ctype;
1067  s_ctype = ( i_ctype == CTYPE_ALL ) ? "All" : s_ctype;
1068 
1069  return s_ctype;
1070 }
1071 
1072 // Static public function to return integer solute-type flag for given string
1073 int US_ModelRecord::stype_flag( const QString s_stype )
1074 {
1075  QString snum = QString( s_stype ).section( ".", 0, 0 );
1076  int i_stype = ( snum.mid( 0, 1 ).toInt() << 6 )
1077  + ( snum.mid( 1, 1 ).toInt() << 3 )
1078  + ( snum.mid( 2, 1 ).toInt() );
1079 
1080  return i_stype;
1081 }
1082 
1083 // Static public function to return solute-type text for given integer flag
1084 QString US_ModelRecord::stype_text( const int i_stype )
1085 {
1086  const char atyp[] = { 's', 'k', 'w', 'v', 'd' };
1087  int ixv = ( i_stype >> 6 ) & 7;
1088  int iyv = ( i_stype >> 3 ) & 7;
1089  int izv = i_stype & 7;
1090  QString s_stype = QString().sprintf( "%03o", i_stype ) + "."
1091  + QString( atyp[ ixv ] )
1092  + QString( atyp[ iyv ] )
1093  + QString( atyp[ izv ] );
1094 
1095  return s_stype;
1096 }
1097 
Generated on Sun Jun 21 2015 11:40:42 for UltraScan III by   doxygen 1.8.3.1-20130324