UltraScan III
us_pcsa_process.cpp
Go to the documentation of this file.
1 #include <QApplication>
3 #include <QtCore>
4 #include <float.h>
5 
6 #include "us_pcsa_process.h"
7 #include "us_util.h"
8 #include "us_settings.h"
9 #include "us_astfem_math.h"
10 #include "us_astfem_rsa.h"
11 #include "us_sleep.h"
12 #include "us_math2.h"
13 #include "us_lm.h"
14 #include "us_solve_sim.h"
15 #include "us_constants.h"
16 #include "us_memory.h"
17 
18 // Class to process PCSA simulations
19 US_pcsaProcess::US_pcsaProcess( QList< US_SolveSim::DataSet* >& dsets,
20  QObject* parent /*=0*/ ) : QObject( parent ), dsets( dsets )
21 {
22  // Pointer to experiment data
23  edata = &dsets[ 0 ]->run_data;
24  // Pointer to simulation parameters
25  simparms = &dsets[ 0 ]->simparams;
26  dbg_level = US_Settings::us_debug();
27  maxrss = 0; // max memory
28 
29  mrecs .clear(); // computed model records
30 
31  nscans = edata->scanCount();
32  npoints = edata->pointCount();
33  cresolu = 100;
34  curvtype = CTYPE_NONE;
35  nmtasks = 0;
36  kctask = 0;
37  kstask = 0;
38  varimin = 9.e+9;
39  minvarx = 99999;
40  st_mask = dsets[ 0 ]->solute_type;
41  attr_x = ( st_mask >> 6 ) & 7;
42  attr_y = ( st_mask >> 3 ) & 7;
43  attr_z = st_mask & 7;
44  double zcurr = dsets[ 0 ]->zcoeffs[ 0 ];
45  double vbar20 = dsets[ 0 ]->vbar20;
46  zcurr = ( zcurr == 0.0 ) ? vbar20 : zcurr;
47 
48  parlims[ 0 ] = -1.0;
49  parlims[ 4 ] = st_mask;
50  parlims[ 5 ] = zcurr;
51 }
52 
53 // Get maximum used memory
54 long int US_pcsaProcess::max_rss( void )
55 {
56  return US_Memory::rss_max( maxrss );
57 }
58 
59 // Start a specified PCSA fit run
60 void US_pcsaProcess::start_fit( double xll, double xul, double yll, double yul,
61  int nyp, int res, int typ, int nth,
62  int noi, int lmmxc, int gfits,
63  double gfthr, double alf )
64 {
65 DbgLv(1) << "PC(pcsaProc): start_fit()";
66  abort = false;
67  xlolim = xll;
68  xuplim = xul;
69  ylolim = yll;
70  yuplim = yul;
71  nypts = nyp;
72  cresolu = res;
73  curvtype = typ;
74  nthreads = nth;
75  noisflag = noi;
76  alpha_scn = false;
77  alpha = alf;
78  alpha_fx = 0.0;
79  alpha_lm = 0.0;
80  fi_itermax = gfits;
81  rd_thresh = gfthr;
82  lmmxcall = lmmxc;
83  st_mask = dsets[ 0 ]->solute_type;
84  attr_x = ( st_mask >> 6 ) & 7;
85  attr_y = ( st_mask >> 3 ) & 7;
86  attr_z = st_mask & 7;
87  double vbar20 = dsets[ 0 ]->vbar20;
88  double zcurr = dsets[ 0 ]->zcoeffs[ 0 ];
89  zcurr = ( zcurr == 0.0 ) ? vbar20 : zcurr;
90  parlims[ 0 ] = -1.0;
91  parlims[ 4 ] = st_mask;
92  parlims[ 5 ] = zcurr;
93 
94  if ( alpha < 0.0 )
95  { // Negative alpha acts as flag
96  if ( alpha == (-99.0) )
97  { // This is an alpha scan run
98  alpha_scn = true;
99  alpha = 0.0;
100  }
101 
102  else
103  { // This is a run with regularize in earlier stage(s)
104  alpha = -alpha;
105  alpha_lm = alpha; // In L-M stage for sure
106 
107  if ( alpha > 1.0 )
108  {
109  alpha -= 1.0;
110  alpha_fx = alpha; // Also in Fixed stage
111  alpha_lm = alpha;
112  }
113  }
114  }
115 
116  pfi_rmsd = 1.0e+99;
117  cfi_rmsd = pfi_rmsd;
118  fi_iter = 1;
119 
120 DbgLv(1) << "PC: alf alpha lm fx scn"
121  << alf << alpha << alpha_lm << alpha_fx << alpha_scn;
122  errMsg = tr( "NO ERROR: start" );
123  //maxrss = 0;
124  varimin = 9.e+9;
125  minvarx = 99999;
126 
127  wkstates .resize( nthreads );
128  wthreads .clear();
129  job_queue.clear();
130 
131  rmsds .clear();
132  orig_sols.clear();
133  mrecs .clear();
134 
135 DbgLv(1) << "PC: sll sul" << xlolim << xuplim
136  << " yll yul nyp" << ylolim << yuplim << nypts
137  << " type reso nth noi" << curvtype << cresolu << nthreads << noisflag;
138 int memmb = qRound( (double)US_Memory::rss_now() / 1024.0 );
139 DbgLv(1) << "STR_STR: nowmem" << memmb << "kthr ksol" << nthreads << cresolu;
140 
141  timer.start(); // start a timer to measure run time
142 
143  // experiment data dimensions
144  nscans = edata->scanCount();
145  npoints = edata->pointCount();
146 
147  if ( curvtype == CTYPE_SL || curvtype == CTYPE_HL )
148  { // Determine models for straight-line curves
149  nmtasks = slmodels( curvtype, xlolim, xuplim, ylolim, yuplim, nypts,
150  cresolu );
151  }
152 
153  else if ( curvtype == CTYPE_IS || curvtype == CTYPE_DS )
154  { // Determine models for sigmoid curves
155  nmtasks = sigmodels( curvtype, xlolim, xuplim, ylolim, yuplim, nypts,
156  cresolu );
157  }
158 
159  else if ( curvtype == CTYPE_2O )
160  { // Determine models for 2nd-order power law curves
161  nmtasks = pl2models( xlolim, xuplim, ylolim, yuplim, nypts, cresolu );
162  }
163 
164  else
165  nmtasks = 0;
166 
167  kcsteps = 0;
168  nctotal = nmtasks;
169  emit stage_complete( kcsteps, nctotal );
170 
171  kctask = 0;
172  kstask = 0;
173  nthreads = ( nthreads < nmtasks ) ? nthreads : nmtasks;
174 DbgLv(1) << "PC: nscans npoints" << nscans << npoints;
175 DbgLv(1) << "PC: nctotal nthreads" << nctotal << nthreads;
176  max_rss();
177 DbgLv(1) << "PC: (1)maxrss" << maxrss;
178  int memava;
179  int memtot;
180  int memuse;
181  int mempav = US_Memory::memory_profile( &memava, &memtot, &memuse );
182  int nsoln = orig_sols[ 0 ].size();
183  int dsize = nscans * npoints;
184  int dfact = ( nsoln + 4 ) * sizeof( double );
185  int memneed = qRound( (double)dsize * (double)dfact / ( 1024. * 1024. ) );
186  memneed *= nthreads;
187 DbgLv(1) << "PC:MEM(1): pav" << mempav << "ns dsz dfac mneed"
188  << nsoln << dsize << dfact << memneed;
189  memuse = memuse + memneed;
190  memava = memtot - memuse;
191  mempav = ( memava * 100 ) / memtot;
192 DbgLv(1) << "PC:MEM(1): pav" << mempav << "ava tot use"
193  << memava << memtot << memuse;
194 
195  // Queue all the tasks
196  for ( int ktask = 0; ktask < nmtasks; ktask++ )
197  {
198  WorkPacketPc wtask;
199  int mm = orig_sols[ ktask ].size() - 1;
200  double stry = orig_sols[ ktask ][ 0 ].y;
201  double endy = orig_sols[ ktask ][ mm ].y;
202  wtask.par1 = mrecs[ ktask ].par1;
203  wtask.par2 = mrecs[ ktask ].par2;
204  wtask.par3 = mrecs[ ktask ].par3;
205  wtask.depth = 0;
206 
207  wtask.sim_vals.alpha = alpha_fx;
208 
209  queue_task( wtask, stry, endy, ktask, noisflag, orig_sols[ ktask ] );
210  }
211 
212  // Start the first threads. This will begin the first work units (subgrids).
213  // Thereafter, work units are started in new threads when previous
214  // threads signal that they have completed their work.
215 
216  for ( int ii = 0; ii < nthreads; ii++ )
217  {
218  wthreads << 0;
219 
220  WorkPacketPc wtask = job_queue.takeFirst();
221  submit_job( wtask, ii );
222  }
223 
224  mrecs .clear();
225  max_rss();
226  kstask = nthreads; // count of started tasks is initially thread count
227 DbgLv(1) << "PC: kstask nthreads" << kstask << nthreads << job_queue.size()
228  << "(2)maxrss" << maxrss;
229 
230  emit message_update( pmessage_head() +
231  tr( "Starting computations of %1 models\n using %2 threads ..." )
232  .arg( nmtasks ).arg( nthreads ), false );
233  max_rss();
234 DbgLv(1) << "PC: (2)maxrss" << maxrss;
235  mempav = US_Memory::memory_profile( &memava, &memtot, &memuse );
236  memuse = memuse + memneed;
237  memava = memtot - memuse;
238  mempav = ( memava * 100 ) / memtot;
239 DbgLv(1) << "PC:MEM(2): pav" << mempav << "ava tot use"
240  << memava << memtot << memuse;
241 
242 #if 0
243  if ( mempav < 20 )
244  {
245  stop_fit();
246  emit process_complete( 6 ); // signal memory-related stop
247  abort = true;
248 
249  QMessageBox::warning( 0, tr( "High Memory Usage" ),
250  tr( "The available memory percent of\n"
251  "the total has fallen below 20%.\n"
252  "Total memory is %1 MB;\n"
253  "Available memory is %2 MB.\n\n"
254  "Re-parameterize the fit with adjusted\n"
255  "Curve Resolution Points and/or\n"
256  "Thread Count." )
257  .arg( memtot ).arg( memava ) );
258  }
259 #endif
260 }
261 
262 // Abort a fit run
263 void US_pcsaProcess::stop_fit()
264 {
265  abort = true;
266 
267  for ( int ii = 0; ii < wthreads.size(); ii++ )
268  {
269 DbgLv(1) << "StopFit test Thread" << ii + 1;
270  WorkerThreadPc* wthr = wthreads[ ii ];
271 
272  if ( wthr != 0 && wthr->isRunning() )
273  {
274  wthr->disconnect();
275  wthr->flag_abort();
276 DbgLv(1) << " STOPTHR: thread aborted";
277  }
278 
279  else if ( wthr != 0 && ! wthr->isFinished() )
280  {
281  wthr->disconnect();
282  delete wthr;
283 DbgLv(1) << " STOPTHR: thread deleted";
284  }
285 
286  wthreads[ ii ] = 0;
287  }
288 
289  job_queue.clear();
290 
291  emit message_update( pmessage_head() +
292  tr( "All computations have been aborted." ), false );
293  abort = false;
294 }
295 
296 // Complete a specified PCSA fit run after alpha scan or alpha change
297 void US_pcsaProcess::final_fit( double alf )
298 {
299 DbgLv(1) << "PC(pcsaProc): final_fit()";
300  abort = false;
301  alpha = alf;
302 
303  compute_final();
304 
305  emit process_complete( 9 ); // signal that all processing is complete
306 }
307 
308 // Slot to handle the low-variance record of fixed final calculated solutes
309 void US_pcsaProcess::process_fxfinal( US_ModelRecord& mrec )
310 {
311  if ( abort ) return;
312 
313  if ( ( noisflag & 1 ) != 0 )
314  { // copy TI noise to caller and internal vector
315  ti_noise.minradius = edata->radius( 0 );
316  ti_noise.maxradius = edata->radius( npoints - 1 );
317  ti_noise.minradius = (double)qRound( ti_noise.minradius * 1e+5 ) * 1e-5;
318  ti_noise.maxradius = (double)qRound( ti_noise.maxradius * 1e+5 ) * 1e-5;
319  ti_noise.values.resize( npoints );
320  ti_noise.count = npoints;
321 
322  for ( int rr = 0; rr < npoints; rr++ )
323  {
324  ti_noise.values[ rr ] = mrec.ti_noise[ rr ];
325  }
326  }
327 
328  if ( ( noisflag & 2 ) != 0 )
329  { // copy RI noise to caller and internal vector
330  ri_noise.values.resize( nscans );
331  ri_noise.count = nscans;
332 
333  for ( int ss = 0; ss < nscans; ss++ )
334  {
335  ri_noise.values[ ss ] = mrec.ri_noise[ ss ];
336  }
337  }
338 DbgLv(1) << "FIN_FIN: ti,ri counts" << ti_noise.count << ri_noise.count;
339 
340  //US_SolveSim::DataSet* dset = dsets[ 0 ];
341  int nsolutes = mrec.csolutes.size();
342  //double vbar20 = dset->vbar20;
343  model.components.resize( nsolutes );
344  qSort( mrec.csolutes );
345 
346  // build the final model
347  for ( int cc = 0; cc < nsolutes; cc++ )
348  {
349  // Get standard-space solute values (20,W)
350  US_Model::SimulationComponent mcomp;
351 #if 0
352  mcomp.vbar20 = vbar20;
353  mcomp.s = mrec.csolutes[ cc ].x;
354  mcomp.D = 0.0;
355  mcomp.mw = 0.0;
356  mcomp.f = 0.0;
357  mcomp.f_f0 = mrec.csolutes[ cc ].y;
358 #endif
359  US_ZSolute::set_mcomp_values( mcomp, mrec.csolutes[ cc ],
360  st_mask, true );
361 
362  // Complete other coefficients in standard-space
363  model.calc_coefficients( mcomp );
364 DbgLv(1) << " Bcc comp D" << mcomp.D;
365 
366  model.components[ cc ] = mcomp;
367  }
368 
369 DbgLv(1) << "FIN_FIN: c0 cn" << mrec.csolutes[0].c
370  << mrec.csolutes[qMax(0,nsolutes-1)].c << " nsols" << nsolutes;
371  nscans = edata->scanCount();
372  npoints = edata->pointCount();
373  US_AstfemMath::initSimData( sdata, *edata, 0.0 );
374  US_AstfemMath::initSimData( rdata, *edata, 0.0 );
375  US_DataIO::RawData* simdat = &mrec.sim_data;
376  US_DataIO::RawData* resids = &mrec.residuals;
377 DbgLv(1) << "FIN_FIN: nscans npoints" << nscans << npoints;
378 DbgLv(1) << "FIN_FIN: simdat nsc npt"
379  << simdat->scanCount() << simdat->pointCount();
380 DbgLv(1) << "FIN_FIN: resids nsc npt"
381  << resids->scanCount() << resids->pointCount();
382 DbgLv(1) << "FIN_FIN: rdata nsc npt"
383  << rdata.scanCount() << rdata.pointCount();
384 DbgLv(1) << "FIN_FIN: sdata nsc npt"
385  << sdata.scanCount() << sdata.pointCount();
386 
387  // build residuals data set (experiment minus simulation minus any noise)
388  for ( int ss = 0; ss < nscans; ss++ )
389  {
390  for ( int rr = 0; rr < npoints; rr++ )
391  {
392  sdata.setValue( ss, rr, simdat->value( ss, rr ) );
393  rdata.setValue( ss, rr, resids->value( ss, rr ) );
394  }
395  }
396 
397 int mms=nscans/2;
398 int mmr=npoints/2;
399 DbgLv(1) << "FIN_FIN: edatm sdatm rdatm" << edata->value(mms,mmr)
400  << sdata.value(mms,mmr) << rdata.value(mms,mmr);
401 
402 DbgLv(1) << "FIN_FIN: vari" << mrec.variance << "bmndx" << mrec.taskx;
403 
404  // determine elapsed time
405  time_fg = timer.elapsed();
406  int ktimes = ( time_fg + 500 ) / 1000;
407  int ktimeh = ktimes / 3600;
408  int ktimem = ( ktimes - ktimeh * 3600 ) / 60;
409 double bvol = dsets[0]->simparams.band_volume;
410 bvol=dsets[0]->simparams.band_forming?bvol:0.0;
411 DbgLv(1) << "done: vari bvol" << mrec.variance << bvol
412  << "slun klun" << xlolim << xuplim << ylolim << yuplim << nypts
413  << "ets" << ktimes;
414  ktimes = ktimes - ktimeh * 3600 - ktimem * 60;
415 
416  // compose final status message
417  QString pmsg = tr( "The Solution has converged...\n"
418  "Threads: %1 ; Models: %2 ; Iterations: %3 ."
419  "\nRun time: " )
420  .arg( nthreads ).arg( nmtasks ).arg( fi_iter );
421 
422  if ( ktimeh > 0 )
423  pmsg += tr( "%1 hr., %2 min., %3 sec.;" )
424  .arg( ktimeh ).arg( ktimem ).arg( ktimes );
425 
426  else
427  pmsg += tr( "%1 min., %2 sec.;" )
428  .arg( ktimem ).arg( ktimes );
429 
430  max_rss();
431  int memmb = qRound( (double)maxrss / 1024.0 );
432 int ksol=mrec.isolutes.size();
433 DbgLv(1) << "FIN_FIN: maxmem" << memmb << "kthr ksol" << nthreads << ksol;
434 
435  pmsg += tr( " Maximum memory used: " ) +
436  QString::number( memmb ) + " MB\n\n" +
437  tr( "The best model (RMSD=%1, %2 solutes, index %3) is:\n" )
438  .arg( mrec.rmsd ).arg( nsolutes ).arg( mrec.taskx );
439  if ( curvtype == CTYPE_SL )
440  {
441  double slopel = (double)( mrec.end_y - mrec.str_y )
442  / (double)( xuplim - xlolim );
443  pmsg += tr( " the line from s,f/f0 %1, %2 to %3, %4 (slope %5)." )
444  .arg( xlolim ).arg( mrec.str_y ).arg( xuplim ).arg( mrec.end_y )
445  .arg( slopel );
446  }
447  else if ( curvtype == CTYPE_IS || curvtype == CTYPE_DS )
448  {
449  pmsg += tr( " the curve with par1=%1 and par2=%2." )
450  .arg( mrec.par1 ).arg( mrec.par2 );
451 DbgLv(1) << "FIN_FIN: par1,par2" << mrec.par1 << mrec.par2;
452  }
453 
454  else if ( curvtype == CTYPE_HL )
455  {
456  pmsg += tr( " the line from x,y %1, %2 to %3, %4." )
457  .arg( xlolim ).arg( mrec.str_y ).arg( xuplim ).arg( mrec.end_y );
458  }
459 
460  else if ( curvtype == CTYPE_2O )
461  {
462  pmsg += tr( " the curve with a=%1, b=%2, c=%3\n"
463  " from x,y %4, %5 to %6, %7." )
464  .arg( mrec.par1 ).arg( mrec.par2 ).arg( mrec.par3 )
465  .arg( xlolim ).arg( mrec.str_y ).arg( xuplim ).arg( mrec.end_y );
466  }
467 
468  // Signal final message
469  emit message_update( pmessage_head() + pmsg, false );
470 
471 DbgLv(1) << "FIN_FIN: Run Time: hr min sec" << ktimeh << ktimem << ktimes;
472 DbgLv(1) << "FIN_FIN: maxrss memmb nthr" << maxrss << memmb << nthreads
473  << " nsubg noisf" << nmtasks << noisflag;
474 DbgLv(1) << "FIN_FIN: kcsteps nctotal" << kcsteps << nctotal;
475 DbgLv(1) << "FIN_FIN: alpha_fx" << alpha_fx;
476  pfi_rmsd = cfi_rmsd;
477  cfi_rmsd = mrecs[ 0 ].rmsd;
478  rd_frac = ( fi_iter > 1 ) ?
479  qAbs( ( pfi_rmsd - cfi_rmsd ) / pfi_rmsd ) :
480  rd_thresh * 2;
481  rmsds << cfi_rmsd;
482 DbgLv(1) << "FIN_FIN: rd_frac rd_thresh" << rd_frac << rd_thresh;
483 
484  if ( fi_iter < fi_itermax && rd_frac > rd_thresh )
485  {
486  restart_fit();
487 
488  return;
489  }
490 
491  emit process_complete( 8 ); // signal that L-M is starting
492 
493  // Compute a best model using Levenberg-Marquardt
494  LevMarq_fit();
495 
496  if ( alpha_scn )
497  { // Signal analysis control that an alpha scan may begin
498  emit process_complete( 7 );
499  return;
500  }
501 
502  if ( alpha_lm == 0.0 && alpha != 0.0 )
503  { // L-M nonregularized, but need to do a final regularized computation
504  compute_final();
505  }
506 
507  nsolutes = model.components.size();
508 
509  emit process_complete( 9 ); // signal that all processing is complete
510 }
511 
512 // Public slot to get results upon completion of all refinements
513 bool US_pcsaProcess::get_results( US_DataIO::RawData* da_sim,
514  US_DataIO::RawData* da_res,
515  US_Model* da_mdl,
516  US_Noise* da_tin,
517  US_Noise* da_rin,
518  int& bm_ndx,
519  QStringList& modstats,
520  QVector< US_ModelRecord >& p_mrecs )
521 {
522  bool all_ok = true;
523  model.alphaRP = alpha;
524  mrecs[ 0 ].model = model;
525  abort = false;
526 DbgLv(1) << "PC:GR: IN abort" << abort;
527 
528 // if ( abort ) return false;
529 
530  *da_sim = sdata; // copy simulation data
531  *da_res = rdata; // copy residuals data
532  *da_mdl = model; // copy model
533 
534  if ( ( noisflag & 1 ) != 0 && da_tin != 0 )
535  {
536  *da_tin = ti_noise; // copy any ti noise
537  mrecs[ 0 ].ti_noise = ti_noise.values;
538  }
539 
540  if ( ( noisflag & 2 ) != 0 && da_rin != 0 )
541  {
542  *da_rin = ri_noise; // copy any ri noise
543  mrecs[ 0 ].ri_noise = ri_noise.values;
544  }
545 
546  p_mrecs = mrecs; // copy model records vector
547  bm_ndx = mrecs[ 0 ].taskx;
548 DbgLv(0) << " GET_RES: ti,ri counts" << ti_noise.count << ri_noise.count;
549 DbgLv(0) << " GET_RES: VARI,RMSD" << mrecs[0].variance << mrecs[0].rmsd
550  << "BM_NDX" << bm_ndx << "ALPHA" << alpha;
551 
552  model_statistics( mrecs, modstats );
553 
554 DbgLv(1) << "PC:GR: RTN";
555  return all_ok;
556 }
557 
558 // Public slot to get best model record in preparation for an alpha scan
559 void US_pcsaProcess::get_mrec( US_ModelRecord& p_mrec )
560 {
561  p_mrec = mrecs[ 0 ]; // Copy best model record
562 }
563 
564 // Replace the top-spot model record
565 void US_pcsaProcess::put_mrec( US_ModelRecord& a_mrec )
566 {
567  mrecs[ 0 ] = a_mrec; // Copy best model record
568 
569  model = a_mrec.model;
570  alpha = 0.0;
571  sdata = a_mrec.sim_data;
572  rdata = a_mrec.residuals;
573 }
574 
575 // Replace the model records list and related variables after file load
576 void US_pcsaProcess::put_mrecs( QVector< US_ModelRecord >& a_mrecs )
577 {
578  mrecs = a_mrecs; // Copy model records list
579  int nmrecs = mrecs.size();
580  int ntinoi = mrecs[ 0 ].ti_noise.count();
581  int nrinoi = mrecs[ 0 ].ri_noise.count();
582 DbgLv(1) << "PC:putMRs: nmrecs" << nmrecs << "rmsd" << mrecs[0].rmsd;
583 
584  if ( nmrecs < 1 )
585  return;
586 
587  model = mrecs[ 0 ].model;
588  sdata = mrecs[ 0 ].sim_data;
589  rdata = mrecs[ 0 ].residuals;
590  xlolim = mrecs[ 0 ].xmin;
591  xuplim = mrecs[ 0 ].xmax;
592  ylolim = mrecs[ 0 ].ymin;
593  yuplim = mrecs[ 0 ].ymax;
594  curvtype = mrecs[ 0 ].ctype;
595  int v_ctype = mrecs[ 0 ].v_ctype;
596 DbgLv(1) << "PC:putMRs: ctype v_ctype" << curvtype << v_ctype;
597 
598  if ( ntinoi > 0 )
599  {
600  ti_noise.values = mrecs[ 0 ].ti_noise;
601  ti_noise.minradius = edata->radius( 0 );
602  ti_noise.maxradius = edata->radius( npoints - 1 );
603  ti_noise.count = ntinoi;
604  }
605 
606  if ( nrinoi > 0 )
607  {
608  ri_noise.values = mrecs[ 0 ].ri_noise;
609  ri_noise.count = nrinoi;
610  }
611 
612  nmtasks = ( mrecs[ 0 ].taskx == mrecs[ 1 ].taskx )
613  ? ( nmrecs - 1 ) : nmrecs;
614  nmtasks = ( mrecs[ 1 ].taskx == mrecs[ 2 ].taskx )
615  ? ( nmtasks - 1 ) : nmtasks;
616  nypts = ( v_ctype != CTYPE_ALL ) ? nmtasks : ( nmtasks / 3 );
617  nypts = ( curvtype != CTYPE_HL && curvtype != CTYPE_2O )
618  ? ( qRound( sqrt( (double)nypts ) ) )
619  : nypts;
620  nypts = ( curvtype == CTYPE_2O )
621  ? ( qRound( pow( (double)nypts, 0.3333333 ) ) )
622  : nypts;
623  alpha = 0.0;
624 DbgLv(1) << "PC:putMRs: nypts nmtasks" << nypts << nmtasks;
625 DbgLv(1) << "PC:putMRs: m0 rmsd" << mrecs[0].rmsd;
626 
627  for ( int ii = 0; ii < nmrecs; ii++ )
628  {
629  varimin = qMin( varimin, mrecs[ ii ].variance );
630  }
631 
632 DbgLv(1) << "PC:putMRs: curvtype" << curvtype << "nypts" << nypts;
633 }
634 
635 // Submit a job
636 void US_pcsaProcess::submit_job( WorkPacketPc& wtask, int thrx )
637 {
638  wtask.thrn = thrx + 1;
639 
640  WorkerThreadPc* wthr = new WorkerThreadPc( this );
641  wthreads[ thrx ] = wthr;
642  wkstates[ thrx ] = WORKING;
643  wtask.sim_vals.maxrss = maxrss;
644 
645  wthr->define_work( wtask );
646 
647  connect( wthr, SIGNAL( work_complete( WorkerThreadPc* ) ),
648  this, SLOT( process_job( WorkerThreadPc* ) ) );
649 DbgLv(1) << "SUBMIT_JOB taskx" << wtask.taskx
650  << "sk ek" << wtask.str_y << wtask.end_y << "maxrss" << maxrss;
651 
652  wthr->start();
653 }
654 
655 // Slot to handle the results of a just-completed worker thread.
656 // Accumulate computed solutes.
657 // If there is more work to do, start a new thread for a new work unit.
658 void US_pcsaProcess::process_job( WorkerThreadPc* wthrd )
659 {
660  if ( abort ) return;
661  WorkPacketPc wresult;
662 
663  wthrd->get_result( wresult ); // get results of thread task
664  int thrn = wresult.thrn; // thread number of task
665  int thrx = thrn - 1; // index into thread list
666  int taskx = wresult.taskx; // task index of task
667  int orecx = mrecs.size(); // output record index
668 DbgLv(1) << "PROCESS_JOB thrn" << thrn << "taskx orecx" << taskx << orecx;
669 //DBG-CONC
670 if (dbg_level>0) for (int mm=0; mm<wresult.csolutes.size(); mm++ ) {
671  if ( wresult.csolutes[mm].c > 100.0 ) {
672  DbgLv(1) << "PJ: CONC=" << wresult.csolutes[mm].c
673  << " x,y" << wresult.csolutes[mm].x << wresult.csolutes[mm].y; } }
674 //DBG-CONC
675  double variance = wresult.sim_vals.variance;
676 
677  // Save variance and model record for this task result
678  US_ModelRecord mrec;
679  mrec.taskx = taskx;
680  mrec.str_y = wresult.str_y;
681  mrec.end_y = wresult.end_y;
682  mrec.par1 = wresult.par1;
683  mrec.par2 = wresult.par2;
684  mrec.par3 = wresult.par3;
685  mrec.variance = variance;
686  mrec.rmsd = ( variance > 0.0 ) ? sqrt( variance ) : 99.9;
687  mrec.isolutes = wresult.isolutes;
688  mrec.csolutes = wresult.csolutes;
689  mrec.ctype = curvtype;
690  mrec.xmin = xlolim;
691  mrec.xmax = xuplim;
692  mrec.ymin = ylolim;
693  mrec.ymax = yuplim;
694  maxrss = qMax( maxrss, wresult.sim_vals.maxrss );
695 
696  if ( variance < varimin )
697  { // Handle a new minimum variance record
698  if ( minvarx < orecx )
699  { // Clear vectors from the previous minimum
700  mrecs[ minvarx ].clear_data();
701 DbgLv(1) << "PJ: CLEAR: VARI VMIN ORECX" << variance << varimin << orecx;
702  }
703  // Save information from the minimum-variance model record
704  varimin = variance;
705  minvarx = orecx;
706  mrec.sim_data = wresult.sim_vals.sim_data;
707  mrec.residuals = wresult.sim_vals.residuals;
708  mrec.ti_noise = wresult.ti_noise;
709  mrec.ri_noise = wresult.ri_noise;
710 DbgLv(1) << "PJ: MINVARX=" << minvarx;
711  }
712  else if ( variance > varimin )
713  { // Clear vectors from a model record that is not minimum-variance
714 DbgLv(1) << "PJ: CLEAR: VARI VMIN MVARX" << variance << varimin << minvarx;
715  mrec.clear_data();
716  }
717 
718  mrecs << mrec; // Append to the vector of model records
719 DbgLv(1) << "THR_FIN: taskx minvarx varimin" << taskx << minvarx << varimin;
720 
721  max_rss();
722 DbgLv(1) << "PJ: (4)maxrss" << maxrss;
723  int memava;
724  int memtot;
725  int memuse;
726  int mempav = US_Memory::memory_profile( &memava, &memtot, &memuse );
727  int nsoln = mrec.isolutes.size();
728  int nscann = wresult.sim_vals.residuals.scanCount();
729  int npointn = wresult.sim_vals.residuals.pointCount();
730  int dsize = nscann * npointn;
731  int dfact = ( nsoln + 4 ) * sizeof( double );
732  int memneed = qRound( (double)dsize * (double)dfact / ( 1024. * 1024. ) );
733 DbgLv(1) << "PC:MEM: (4)nsol,nscn,npt,dsize,dfact,memneed"
734  << nsoln << nscann << npointn << dsize << dfact << memneed;
735  memuse = memuse + memneed;
736  memava = memtot - memuse;
737  mempav = ( memava * 100 ) / memtot;
738 DbgLv(1) << "PC:MEM: (4)AvPercent" << mempav;
739 
740  free_worker( thrx ); // Free up this worker thread
741 
742  if ( abort )
743  return;
744 
745  kctask++; // Bump count of completed subgrid tasks
746  emit progress_update( variance ); // Pass progress on to control window
747 DbgLv(1) << "THR_FIN: thrn" << thrn << " taskx orecx" << taskx << orecx
748  << " kct kst" << kctask << kstask;
749 
750  emit message_update( pmessage_head() +
751  tr( "Of %1 models, computations are complete for %2." )
752  .arg( nmtasks ).arg( kctask ), false );
753 
754  if ( kctask >= nmtasks )
755  { // All model tasks are now complete
756  emit stage_complete( kcsteps, nctotal );
757 
758  emit message_update( pmessage_head() +
759  tr( "All models have now been completed. Evaluating..." ), false );
760  }
761 
762  if ( kctask < nmtasks )
763  { // Not the last: add to the queue if necessary
764  // Submit jobs while queue is not empty and a worker thread is ready
765  while ( ! job_queue.isEmpty() &&
766  ( thrx = wkstates.indexOf( READY ) ) >= 0 )
767  {
768  WorkPacketPc wtask = next_job();
769 
770  submit_job( wtask, thrx );
771  kstask++; // Bump count of started worker threads
772 
773  }
774  }
775 
776  else
777  { // We have done the last computation, so determine the low-rmsd result
778 
779  qSort( mrecs ); // Sort model records by variance
780 
781 DbgLv(1) << "THR_FIN: thrn" << thrn << " mrecs size" << mrecs.size()
782  << "mrec0 taskx,vari" << mrecs[0].taskx << mrecs[0].variance;
783  process_fxfinal( mrecs[ 0 ] );
784 //*DBG*
785 if (dbg_level>0) {
786  for (int ii=0;ii<mrecs.size();ii++)
787  {
788  DbgLv(1) << "MREC:rank" << ii << "taskx" << mrecs[ii].taskx
789  << "st_k en_k" << mrecs[ii].str_y << mrecs[ii].end_y
790  << "vari,rmsd" << mrecs[ii].variance << mrecs[ii].rmsd
791  << "ncsols" << mrecs[ii].csolutes.size();
792  }
793 }
794 //*DBG*
795  }
796 }
797 
798 // Build a task and add it to the queue
799 void US_pcsaProcess::queue_task( WorkPacketPc& wtask, double stry, double endy,
800  int taskx, int noisf, QVector< US_ZSolute > isolutes )
801 {
802  wtask.thrn = 0; // thread number (none while queued)
803  wtask.taskx = taskx; // task index
804  wtask.noisf = noisf; // noise flag
805  wtask.state = READY; // initialized state, ready for submit
806  wtask.str_y = stry; // start k value
807  wtask.end_y = endy; // end k value
808  wtask.dsets = dsets; // pointer to experiment data
809  wtask.isolutes = isolutes; // solutes for calc_residuals task
810 
811  wtask.csolutes.clear(); // clear output vectors
812  wtask.ti_noise.clear();
813  wtask.ri_noise.clear();
814 
815  job_queue << wtask; // put the task on the queue
816 }
817 
818 // Free up a worker thread
819 void US_pcsaProcess::free_worker( int tx )
820 {
821  if ( tx >= 0 && tx < wthreads.size() )
822  {
823  if ( wthreads[ tx ] != 0 )
824  delete wthreads[ tx ]; // destroy thread
825 
826  wthreads[ tx ] = 0; // set thread pointer to null
827  wkstates[ tx ] = READY; // mark its slot as available
828  }
829 }
830 
831 QString US_pcsaProcess::pmessage_head()
832 {
833  const char* ctp[] = { "None",
834  "Straight Line",
835  "Increasing Sigmoid",
836  "Unknown-3",
837  "Decreasing Sigmoid",
838  "Unknown-5",
839  "Unknown-6",
840  "All (SL+IS+DS)",
841  "Horizontal Line [ C(s) ]",
842  "Unknown-9",
843  "Unknown-10",
844  "Unknown-11",
845  "Unknown-12",
846  "Unknown-13",
847  "Unknown-14",
848  "Unknown-15",
849  "Second-Order Power Law",
850  "?UNKNOWN?"
851  };
852  QString ctype = QString( ctp[ curvtype ] );
853  QString hmsg = tr( "Analysis of %1 %2 %3-solute models.\n" )
854  .arg( nmtasks ).arg( ctype ).arg( cresolu );
855 
856  for ( int ii = 1; ii < fi_iter; ii++ )
857  hmsg += tr( "Iteration %1 had a best RMSD of %2 ;\n" )
858  .arg( ii ).arg( rmsds[ ii - 1 ] );
859 
860  hmsg += tr( "Grid Fit Iteration %1.\n" ).arg( fi_iter );
861 
862  return hmsg;
863 }
864 
865 // Get next job from queue, insuring we get the lowest depth
866 WorkPacketPc US_pcsaProcess::next_job()
867 {
868  WorkPacketPc wtask;
869  if ( job_queue.size() == 0 ) return wtask;
870 
871  int jobx = 0;
872  wtask = job_queue[ jobx ];
873 
874 if(jobx>0) {
875 DbgLv(1) << "NEXTJ: jobx taskx taskx0"
876  << jobx << wtask.taskx << job_queue[0].taskx; }
877 else {
878 DbgLv(1) << "NEXTJ: jobx taskx" << jobx << wtask.taskx; }
879 DbgLv(1) << "NEXTJ: wtask" << &wtask << &job_queue[jobx];
880 
881  job_queue.removeAt( jobx ); // Remove job from queue
882  return wtask;
883 }
884 
885 // Build all the straight-line models
886 int US_pcsaProcess::slmodels( int ctp, double xlo, double xup,
887  double ylo, double yup, int nyp, int res )
888 {
889 DbgLv(1) << "SLMO: xlo xup ylo yup nyp res" << xlo << xup << ylo << yup
890  << nyp << res;
891  int nmodels = 0;
892  //double vbar20 = dsets[ 0 ]->vbar20;
893  orig_sols.clear();
894 
895  // Compute straight-line model records
896  if ( ctp == CTYPE_SL )
897  nmodels = US_ModelRecord::compute_slines( xlo, xup, ylo, yup, nyp,
898  res, parlims, mrecs );
899  else if ( ctp == CTYPE_HL )
900  nmodels = US_ModelRecord::compute_hlines( xlo, xup, ylo, yup, nyp,
901  res, parlims, mrecs );
902 
903  // Add model records solutes to task solutes list
904  for ( int ii = 0; ii < nmodels; ii++ )
905  {
906  orig_sols << mrecs[ ii ].isolutes;
907  }
908 DbgLv(1) << "SLMO: orig_sols size" << orig_sols.size() << "nmodels" << nmodels;
909 
910  return nmodels;
911 }
912 
913 // Build all the sigmoid models
914 int US_pcsaProcess::sigmodels( int ctp, double xlo, double xup,
915  double ylo, double yup, int nyp, int nlpts )
916 {
917 DbgLv(1) << "SGMO: ctp xlo xup ylo yup nyp nlp" << ctp << xlo << xup
918  << ylo << yup << nyp << nlpts;
919  //double vbar20 = dsets[ 0 ]->vbar20;
920  orig_sols.clear();
921 
922  // Compute sigmoid model records
923  int nmodels = US_ModelRecord::compute_sigmoids( ctp, xlo, xup, ylo, yup,
924  nyp, nlpts, parlims, mrecs );
925 
926  // Add model records solutes to task solutes list
927  for ( int ii = 0; ii < nmodels; ii++ )
928  {
929  orig_sols << mrecs[ ii ].isolutes;
930  }
931 DbgLv(1) << "SGMO: orig_sols size" << orig_sols.size() << "nmodels" << nmodels;
932 
933  return nmodels;
934 }
935 
936 // Build all the 2nd-order power law models
937 int US_pcsaProcess::pl2models( double xlo, double xup, double ylo, double yup,
938  int nyp, int nlpts )
939 {
940 DbgLv(1) << "2OMO: xlo xup ylo yup nyp nlp" << xlo << xup
941  << ylo << yup << nyp << nlpts;
942  //double vbar20 = dsets[ 0 ]->vbar20;
943  orig_sols.clear();
944 
945  // Compute 2nd-order model records
946  int nmodels = US_ModelRecord::compute_2ndorder( xlo, xup, ylo, yup,
947  nyp, nlpts, parlims, mrecs );
948 
949  // Add model records solutes to task solutes list
950  for ( int ii = 0; ii < nmodels; ii++ )
951  {
952  orig_sols << mrecs[ ii ].isolutes;
953  }
954 DbgLv(1) << "2OMO: orig_sols size" << orig_sols.size() << "nmodels" << nmodels;
955 
956  return nmodels;
957 }
958 
959 // Generate the strings of model statistics for a report
960 void US_pcsaProcess::model_statistics( QVector< US_ModelRecord >& mrecs,
961  QStringList& modstats )
962 {
963  const char* ctp[] = { "None",
964  "Straight Line",
965  "Increasing Sigmoid",
966  "Unknown-3",
967  "Decreasing Sigmoid",
968  "Unknown-5",
969  "Unknown-6",
970  "All (SL+IS+DS)",
971  "Horizontal Line [ C(s) ]",
972  "Unknown-9",
973  "Unknown-10",
974  "Unknown-11",
975  "Unknown-12",
976  "Unknown-13",
977  "Unknown-14",
978  "Unknown-15",
979  "Second-Order Power Law",
980  "?UNKNOWN?"
981  };
982  const int lenctp = sizeof( ctp ) / sizeof( ctp[ 0 ] );
983 
984  // Accumulate the statistics
985  int nbmods = nmtasks / 5;
986  nbmods = qMin( nmtasks - 1, qMax( 3, nbmods ) );
987  int nlpts = cresolu;
988  int nblpts = mrecs[ 0 ].isolutes.size();
989  int soltype = dsets[ 0 ]->solute_type;
990  double rmsdmin = 99999.0;
991  double rmsdmax = 0.0;
992  double rmsdavg = 0.0;
993  double brmsmin = 99999.0;
994  double brmsmax = 0.0;
995  double brmsavg = 0.0;
996 DbgLv(1) << "PC:MS: nmtasks mrssiz nbmods" << nmtasks << mrecs.size() << nbmods;
997  double rmsdmed = mrecs[ nmtasks / 2 ].rmsd;
998  double brmsmed = mrecs[ nbmods / 2 ].rmsd;
999  int nsolmin = 999999;
1000  int nsolmax = 0;
1001  int nsolavg = 0;
1002  int nbsomin = 999999;
1003  int nbsomax = 0;
1004  int nbsoavg = 0;
1005 
1006  for ( int ii = 0; ii < nmtasks; ii++ )
1007  {
1008  double rmsd = mrecs[ ii ].rmsd;
1009  int nsols = mrecs[ ii ].csolutes.size();
1010  rmsdmin = qMin( rmsdmin, rmsd );
1011  rmsdmax = qMax( rmsdmax, rmsd );
1012  rmsdavg += rmsd;
1013  nsolmin = qMin( nsolmin, nsols );
1014  nsolmax = qMax( nsolmax, nsols );
1015  nsolavg += nsols;
1016 
1017  if ( ii < nbmods )
1018  {
1019  brmsmin = qMin( brmsmin, rmsd );
1020  brmsmax = qMax( brmsmax, rmsd );
1021  brmsavg += rmsd;
1022  nbsomin = qMin( nbsomin, nsols );
1023  nbsomax = qMax( nbsomax, nsols );
1024  nbsoavg += nsols;
1025  }
1026  }
1027 
1028  rmsdavg /= (double)nmtasks;
1029  nsolavg = ( nsolavg + nmtasks / 2 ) / nmtasks;
1030  brmsavg /= (double)nbmods;
1031  nbsoavg = ( nbsoavg + nbmods / 2 ) / nbmods;
1032 
1033  modstats.clear();
1034 
1035  int curvtx = qMin( curvtype, lenctp - 1 );
1036  modstats << tr( "Curve Type:" )
1037  << QString( ctp[ curvtx ] );
1038  modstats << tr( "Solute Type:" )
1039  << US_ModelRecord::stype_text( soltype );
1040  modstats << tr( "X Range:" )
1041  << QString().sprintf( "%10.4f %10.4f", xlolim, xuplim );
1042  double str_y = mrecs[ 0 ].str_y;
1043  double end_y = mrecs[ 0 ].end_y;
1044 
1045  if ( curvtype == CTYPE_SL || curvtype == CTYPE_HL )
1046  {
1047  double slope = ( end_y - str_y ) / ( xuplim - xlolim );
1048  double yincr = ( yuplim - ylolim ) / (double)( nypts - 1 );
1049  modstats << tr( "Y Range + delta:" )
1050  << QString().sprintf( "%10.4f %10.4f %10.4f",
1051  ylolim, yuplim, yincr );
1052  modstats << tr( "Best curve Y end points + slope:" )
1053  << QString().sprintf( "%10.4f %10.4f %10.4f",
1054  str_y, end_y, slope );
1055 DbgLv(1) << "PC:MS: best str_y,end_y" << str_y << end_y;
1056  }
1057  else if ( curvtype == CTYPE_IS || curvtype == CTYPE_DS )
1058  {
1059  int p1ndx = mrecs[ 0 ].taskx / nypts;
1060  int p2ndx = mrecs[ 0 ].taskx - ( p1ndx * nypts );
1061  double krng = (double)( nypts - 1 );
1062  double p1off = (double)p1ndx / krng;
1063  double par1 = exp( log( 0.001 )
1064  + ( log( 0.5 ) - log( 0.001 ) ) * p1off );
1065  double par2 = (double)p2ndx / krng;
1066  modstats << tr( "Y Range:" )
1067  << QString().sprintf( "%10.4f %10.4f", ylolim, yuplim );
1068  modstats << tr( "Best curve par1 and par2:" )
1069  << QString().sprintf( "%10.4f %10.4f", par1, par2 );
1070  modstats << tr( "Best curve Y end points:" )
1071  << QString().sprintf( "%10.4f %10.4f",
1072  mrecs[ 0 ].str_y, mrecs[ 0 ].end_y );
1073  }
1074  else if ( curvtype == CTYPE_2O )
1075  {
1076  modstats << tr( "Y Range:" )
1077  << QString().sprintf( "%10.4f %10.4f", ylolim, yuplim );
1078  modstats << tr( "Best curve A, B, C:" )
1079  << QString().sprintf( "%10.4f %10.4f %10.4f",
1080  mrecs[ 0 ].par1, mrecs[ 0 ].par2, mrecs[ 0 ].par3 );
1081  modstats << tr( "Best curve Y end points:" )
1082  << QString().sprintf( "%10.4f %10.4f",
1083  mrecs[ 0 ].str_y, mrecs[ 0 ].end_y );
1084  }
1085  else
1086  {
1087  modstats << tr( "Y Range:" )
1088  << QString().sprintf( "%10.4f %10.4f", ylolim, yuplim );
1089  modstats << tr( "Best curve par1 and par2:" )
1090  << QString().sprintf( "%10.4f %10.4f",
1091  mrecs[ 0 ].par1, mrecs[ 0 ].par2 );
1092  modstats << tr( "Best curve Y end points:" )
1093  << QString().sprintf( "%10.4f %10.4f",
1094  mrecs[ 0 ].str_y, mrecs[ 0 ].end_y );
1095  }
1096  modstats << tr( "Number of models:" )
1097  << QString().sprintf( "%5d", nmtasks );
1098  modstats << tr( "Number of curve variations:" )
1099  << QString().sprintf( "%5d", nypts );
1100  modstats << tr( "Solute points per curve:" )
1101  << QString().sprintf( "%5d", nlpts );
1102  if ( nblpts != nlpts )
1103  modstats << tr( "Best-curve points per curve:" )
1104  << QString().sprintf( "%5d", nblpts );
1105  modstats << tr( "Index of best model:" )
1106  << QString().sprintf( "%5d", mrecs[ 0 ].taskx );
1107  modstats << tr( "Best curve calculated solutes:" )
1108  << QString().sprintf( "%5d", mrecs[ 0 ].csolutes.size() );
1109  modstats << tr( "Minimum, Maximum calculated solutes:" )
1110  << QString().sprintf( "%5d %5d", nsolmin, nsolmax );
1111  modstats << tr( "Average calculated solutes:" )
1112  << QString().sprintf( "%5d", nsolavg );
1113  modstats << tr( "Minimum variance:" )
1114  << QString().sprintf( "%12.6e", varimin );
1115  modstats << tr( "Minimum, Maximum rmsd:" )
1116  << QString().sprintf( "%12.8f %12.8f", rmsdmin, rmsdmax );
1117  modstats << tr( "Average, Median rmsd:" )
1118  << QString().sprintf( "%12.8f %12.8f", rmsdavg, rmsdmed );
1119  modstats << tr( "Number of \"better\" models:" )
1120  << QString().sprintf( "%5d", nbmods );
1121  modstats << tr( "%1 Best Min,Max calculated solutes:" ).arg( nbmods )
1122  << QString().sprintf( "%5d %5d", nbsomin, nbsomax );
1123  modstats << tr( "%1 Best Average calculated solutes:" ).arg( nbmods )
1124  << QString().sprintf( "%5d", nbsoavg );
1125  modstats << tr( "%1 Best Minimum, Maximum rmsd:" ).arg( nbmods )
1126  << QString().sprintf( "%12.8f %12.8f", brmsmin, brmsmax );
1127  modstats << tr( "%1 Best Average, Median rmsd:" ).arg( nbmods )
1128  << QString().sprintf( "%12.8f %12.8f", brmsavg, brmsmed );
1129  modstats << tr( "Tikhonov regularization parameter:" )
1130  << QString().sprintf( "%12.3f", alpha );
1131 DbgLv(1) << "PC:MS: RTN";
1132 
1133 }
1134 
1135 // Do curve-fit evaluate function (return RMSD) for a Straight Line model
1136 double US_pcsaProcess::fit_function_SL( double t, double* par )
1137 {
1138  static int ffcall=0; // Fit function call counter
1139  static double epar[ 20 ]; // Static array for holding parameters
1140  static const int nepar = sizeof( epar ) / sizeof( epar[ 0 ] );
1141 
1142  if ( t > 0.0 )
1143  { // If not t[0], return immediately
1144  return 0.0;
1145  }
1146 
1147  if ( t < 0.0 )
1148  { // If t is special flag (-ve.), reset ffcall and return
1149  ffcall = 0;
1150  return 0.0;
1151  }
1152 
1153  QTime ftimer;
1154  ftimer.start();
1155  QList< US_SolveSim::DataSet* > dsets;
1156  void** iparP = (void**)par;
1157  void** parP = (void**)epar;
1158  double par1 = par[ 0 ];
1159  double par2 = par[ 1 ];
1160  int px = ( sizeof( double ) / sizeof( void* ) ) * 2;
1161  if ( ffcall == 0 )
1162  { // On 1st call, copy par array to internal static one
1163  for ( int ii = 0; ii < nepar - 2; ii++ )
1164  epar[ ii ] = par[ ii + 2 ];
1165  parP[ 0 ] = iparP[ px ];
1166  }
1167 
1168  ffcall++; // Bump function call counter
1169  dsets << (US_SolveSim::DataSet*)parP[ 0 ];
1170  int nlpts = (int)epar[ 1 ]; // Get limit parameters
1171  double xmin = epar[ 2 ];
1172  double xmax = epar[ 3 ];
1173  double ylow = epar[ 6 ];
1174  double yhigh = epar[ 7 ];
1175  double p1lo = epar[ 8 ];
1176  double p1hi = epar[ 9 ];
1177  double p2lo = epar[ 10 ];
1178  double p2hi = epar[ 11 ];
1179  int noisfl = (int)epar[ 12 ];
1180  int dbg_lv = (int)epar[ 13 ];
1181  double alpha = epar[ 14 ];
1182  int attr_x = (int)epar[ 16 ];
1183  int attr_y = (int)epar[ 17 ];
1184  //int attr_z = (int)epar[ 18 ];
1185  //int zbase = epar[ 19 ];
1186  double ystart = par1;
1187  double yend = ystart + par2 * ( xmax - xmin );
1188 
1189  // After 1st few calls, test if parameters are within limits
1190  if ( ffcall > 3 )
1191  {
1192  // Leave a little wiggle room on limits
1193  ylow -= 0.1;
1194  yhigh += 0.1;
1195  p1lo -= 0.1;
1196  p1hi += 0.1;
1197  p2lo -= 0.01;
1198  p2hi += 0.02;
1199 
1200  // If this record is beyond any limit, return now with it marked as bad
1201  if ( par1 < p1lo || par2 < p2lo ||
1202  par1 > p1hi || par2 > p2hi ||
1203  ystart < ylow || yend < ylow ||
1204  ystart > yhigh || yend > yhigh )
1205  {
1206 qDebug() << "ffSL: call" << ffcall << "par1 par2" << par1 << par2
1207  << "ys ye" << ystart << yend << "*OUT-OF-LIMITS*";
1208  return 1e+99;
1209  }
1210  }
1211 
1212  double prange = (double)( nlpts - 1 );
1213  double xinc = ( xmax - xmin ) / prange;
1214  double yinc = ( yend - ystart ) / prange;
1215  double vbar20 = dsets[ 0 ]->vbar20;
1216  double xcurr = xmin;
1217  double ycurr = ystart;
1218  double zcurr = dsets[ 0 ]->zcoeffs[ 0 ];
1219  zcurr = ( zcurr == 0.0 ) ? vbar20 : zcurr;
1220 qDebug() << "ffSL: zcurr vb" << zcurr << vbar20;
1221 
1222  if ( attr_x == ATTR_S )
1223  { // X is s: scale by 1e-13
1224  xcurr *= 1e-13;
1225  xinc *= 1e-13;
1226  }
1227  if ( attr_y == ATTR_S )
1228  { // Y is s: scale by 1e-13
1229  ycurr *= 1e-13;
1230  yinc *= 1e-13;
1231  }
1232  US_SolveSim::Simulation sim_vals;
1233  sim_vals.noisflag = noisfl;
1234  sim_vals.dbg_level = dbg_lv;
1235  sim_vals.alpha = alpha;
1236 
1237  for ( int ii = 0; ii < nlpts; ii++ )
1238  { // Fill the input solutes vector
1239  sim_vals.zsolutes << US_ZSolute( xcurr, ycurr, zcurr, 0.0 );
1240  xcurr += xinc;
1241  ycurr += yinc;
1242  }
1243 
1244  // Evaluate the model
1245  double rmsd = evaluate_model( dsets, sim_vals );
1246 
1247  epar[ 15 ] = rmsd;
1248  int ktimms = ftimer.elapsed();
1249 qDebug() << "ffSL: call" << ffcall << "par1 par2" << par1 << par2
1250  << "rmsd" << rmsd << "eval time" << ktimms << "ms.";
1251 if(ffcall<6)
1252 qDebug() << "ffSL: epar0 epar1-9" << parP[0] << epar[1] << epar[2] << epar[3]
1253  << epar[4] << epar[5] << epar[6] << epar[7] << epar[8] << epar[9];
1254 //qDebug() << "ffSL: dsets[0]" << dsets[0] << parP[0]
1255 // << "dsets[0]->vbar20" << dsets[0]->vbar20;
1256 qDebug() << "ffSL: ys ye yl yh" << ystart << yend << ylow << yhigh;
1257 
1258  return rmsd;
1259 }
1260 
1261 // Do curve-fit evaluate function (return RMSD) for a Increasing Sigmoid model
1262 double US_pcsaProcess::fit_function_IS( double t, double* par )
1263 {
1264  static int ffcall=0; // Fit function call counter
1265  static double epar[ 18 ]; // Static array for holding parameters
1266  static const int nepar = sizeof( epar ) / sizeof( epar[ 0 ] );
1267 
1268  if ( t > 0.0 )
1269  { // If not t[0], return immediately
1270  return 0.0;
1271  }
1272 
1273  if ( t < 0.0 )
1274  { // If t is special flag (-ve.), reset ffcall and return
1275  ffcall = 0;
1276  return 0.0;
1277  }
1278 
1279  QTime ftimer;
1280  ftimer.start();
1281  QList< US_SolveSim::DataSet* > dsets;
1282  void** iparP = (void**)par;
1283  void** parP = (void**)epar;
1284  double par1 = par[ 0 ];
1285  double par2 = par[ 1 ];
1286  int px = ( sizeof( double ) / sizeof( void* ) ) * 2;
1287  if ( ffcall == 0 )
1288  { // On 1st call, copy par array to internal static one
1289  for ( int ii = 0; ii < nepar - 2; ii++ )
1290  epar[ ii ] = par[ ii + 2 ];
1291  parP[ 0 ] = iparP[ px ];
1292  }
1293 
1294  ffcall++; // Bump function call counter
1295  dsets << (US_SolveSim::DataSet*)parP[ 0 ];
1296  int nlpts = (int)epar[ 1 ]; // Get limit parameters
1297  double xmin = epar[ 2 ];
1298  double xmax = epar[ 3 ];
1299  double ymin = epar[ 4 ];
1300  double ymax = epar[ 5 ];
1301  double ylow = epar[ 6 ];
1302  double yhigh = epar[ 7 ];
1303  double p1lo = epar[ 8 ];
1304  double p1hi = epar[ 9 ];
1305  double p2lo = epar[ 10 ];
1306  double p2hi = epar[ 11 ];
1307  int noisfl = (int)epar[ 12 ];
1308  int dbg_lv = (int)epar[ 13 ];
1309  double alpha = epar[ 14 ];
1310  int attr_x = (int)epar[ 16 ];
1311  int attr_y = (int)epar[ 17 ];
1312  //int attr_z = (int)epar[ 18 ];
1313  //int zbase = epar[ 19 ];
1314  double ystr = ymin;
1315  double ydif = ymax - ymin;
1316  double xrange = xmax - xmin;
1317  double p1fac = sqrt( 2.0 * qMax( par1, p1lo ) );
1318  double ystart = ystr + ydif * ( 0.5 * erf( ( 0.0 - par2 ) / p1fac ) + 0.5 );
1319  double yend = ystr + ydif * ( 0.5 * erf( ( 1.0 - par2 ) / p1fac ) + 0.5 );
1320 
1321  // After 1st few calls, test if parameters are within limits
1322  if ( ffcall > 3 )
1323  {
1324  // Leave a little wiggle room on limits
1325  ylow -= 0.1;
1326  yhigh += 0.1;
1327  p1lo -= 0.00001;
1328  p1hi += 0.00001;
1329  p2lo -= 0.01;
1330  p2hi += 0.01;
1331 
1332  // If this record is beyond any limit, return now with it marked as bad
1333  if ( par1 < p1lo || par2 < p2lo ||
1334  par1 > p1hi || par2 > p2hi ||
1335  ystart < ylow || yend < ylow ||
1336  ystart > yhigh || yend > yhigh )
1337  {
1338 qDebug() << "ffIS: call" << ffcall << "par1 par2" << par1 << par2
1339  << "ys ye" << ystart << yend << "*OUT-OF-LIMITS*";
1340  return 1e+99;
1341  }
1342  }
1343 
1344  double prange = (double)( nlpts - 1 );
1345  double xinc = 1.0 / prange;
1346  double vbar20 = dsets[ 0 ]->vbar20;
1347  double xcurr = xmin;
1348  double ycurr = ymin;
1349  double zcurr = dsets[ 0 ]->zcoeffs[ 0 ];
1350  zcurr = ( zcurr == 0.0 ) ? vbar20 : zcurr;
1351 qDebug() << "ffIS: zcurr vb" << zcurr << vbar20;
1352  if ( attr_x == ATTR_S )
1353  { // X is s: scale by 1e-13
1354  xmin *= 1e-13;
1355  xrange *= 1e-13;
1356  }
1357  if ( attr_y == ATTR_S )
1358  { // Y is s: scale by 1e-13
1359  ymin *= 1e-13;
1360  ydif *= 1e-13;
1361  }
1362  US_SolveSim::Simulation sim_vals;
1363  sim_vals.noisflag = noisfl;
1364  sim_vals.dbg_level = dbg_lv;
1365  sim_vals.alpha = alpha;
1366 
1367  double xoff = 0.0;
1368 
1369  for ( int ii = 0; ii < nlpts; ii++ )
1370  { // Fill the input solutes vector
1371  double efac = 0.5 * erf( ( xoff - par2 ) / p1fac ) + 0.5;
1372  xcurr = xmin + xoff * xrange;
1373  ycurr = ymin + ydif * efac;
1374  sim_vals.zsolutes << US_ZSolute( xcurr, ycurr, zcurr, 0.0 );
1375  xoff += xinc;
1376  }
1377 
1378  // Evaluate the model
1379  double rmsd = evaluate_model( dsets, sim_vals );
1380 
1381  epar[ 15 ] = rmsd;
1382  int ktimms = ftimer.elapsed();
1383 qDebug() << "ffIS: call" << ffcall << "par1 par2" << par1 << par2
1384  << "rmsd" << rmsd << "eval time" << ktimms << "ms.";
1385 if(ffcall<6)
1386 qDebug() << "ffIS: epar0 epar1-9" << parP[0] << epar[1] << epar[2] << epar[3]
1387  << epar[4] << epar[5] << epar[6] << epar[7] << epar[8] << epar[9];
1388 
1389  return rmsd;
1390 }
1391 
1392 // Do curve-fit evaluate function (return RMSD) for a Decreasing Sigmoid model
1393 double US_pcsaProcess::fit_function_DS( double t, double* par )
1394 {
1395  static int ffcall=0; // Fit function call counter
1396  static double epar[ 18 ]; // Static array for holding parameters
1397  static const int nepar = sizeof( epar ) / sizeof( epar[ 0 ] );
1398 
1399  if ( t > 0.0 )
1400  { // If not t[0], return immediately
1401  return 0.0;
1402  }
1403 
1404  if ( t < 0.0 )
1405  { // If t is special flag (-ve.), reset ffcall and return
1406  ffcall = 0;
1407  return 0.0;
1408  }
1409  QTime ftimer;
1410  ftimer.start();
1411  QList< US_SolveSim::DataSet* > dsets;
1412  void** iparP = (void**)par;
1413  void** parP = (void**)epar;
1414  double par1 = par[ 0 ];
1415  double par2 = par[ 1 ];
1416  int px = ( sizeof( double ) / sizeof( void* ) ) * 2;
1417  if ( ffcall == 0 )
1418  { // On 1st call, copy par array to internal static one
1419  for ( int ii = 0; ii < nepar - 2; ii++ )
1420  epar[ ii ] = par[ ii + 2 ];
1421  parP[ 0 ] = iparP[ px ];
1422  }
1423 
1424  ffcall++; // Bump function call counter
1425  dsets << (US_SolveSim::DataSet*)parP[ 0 ];
1426  int nlpts = (int)epar[ 1 ]; // Get limit parameters
1427  double xmin = epar[ 2 ];
1428  double xmax = epar[ 3 ];
1429  double ymin = epar[ 4 ];
1430  double ymax = epar[ 5 ];
1431  double ylow = epar[ 6 ];
1432  double yhigh = epar[ 7 ];
1433  double p1lo = epar[ 8 ];
1434  double p1hi = epar[ 9 ];
1435  double p2lo = epar[ 10 ];
1436  double p2hi = epar[ 11 ];
1437  int noisfl = (int)epar[ 12 ];
1438  int dbg_lv = (int)epar[ 13 ];
1439  double alpha = epar[ 14 ];
1440  int attr_x = (int)epar[ 16 ];
1441  int attr_y = (int)epar[ 17 ];
1442  //int attr_z = (int)epar[ 18 ];
1443  //int zbase = epar[ 19 ];
1444  double ystr = ymax;
1445  double ydif = ymin - ymax;
1446  double xrange = xmax - xmin;
1447  double p1fac = sqrt( 2.0 * qMax( par1, p1lo ) );
1448  double ystart = ystr + ydif * ( 0.5 * erf( ( 0.0 - par2 ) / p1fac ) + 0.5 );
1449  double yend = ystr + ydif * ( 0.5 * erf( ( 1.0 - par2 ) / p1fac ) + 0.5 );
1450 
1451  // After 1st few calls, test if parameters are within limits
1452  if ( ffcall > 3 )
1453  {
1454  // Leave a little wiggle room on limits
1455  ylow -= 0.1;
1456  yhigh += 0.1;
1457  p1lo -= 0.00001;
1458  p1hi += 0.00001;
1459  p2lo -= 0.01;
1460  p2hi += 0.01;
1461 
1462  // If this record is beyond any limit, return now with it marked as bad
1463  if ( par1 < p1lo || par2 < p2lo ||
1464  par1 > p1hi || par2 > p2hi ||
1465  ystart < ylow || yend < ylow ||
1466  ystart > yhigh || yend > yhigh )
1467  {
1468 qDebug() << "ffDS: call" << ffcall << "par1 par2" << par1 << par2
1469  << "ys ye" << ystart << yend << "*OUT-OF-LIMITS*";
1470  return 1e+99;
1471  }
1472  }
1473 
1474  double prange = (double)( nlpts - 1 );
1475  double xoinc = 1.0 / prange;
1476  double vbar20 = dsets[ 0 ]->vbar20;
1477  double xcurr = xmin;
1478  double ycurr = ymin;
1479  double zcurr = dsets[ 0 ]->zcoeffs[ 0 ];
1480  zcurr = ( zcurr == 0.0 ) ? vbar20 : zcurr;
1481 qDebug() << "ffDS: zcurr vb" << zcurr << vbar20;
1482  if ( attr_x == ATTR_S )
1483  { // X is s: scale by 1e-13
1484  xmin *= 1e-13;
1485  xrange *= 1e-13;
1486  }
1487  if ( attr_y == ATTR_S )
1488  { // Y is s: scale by 1e-13
1489  ymin *= 1e-13;
1490  ydif *= 1e-13;
1491  }
1492  US_SolveSim::Simulation sim_vals;
1493  sim_vals.noisflag = noisfl;
1494  sim_vals.dbg_level = dbg_lv;
1495  sim_vals.alpha = alpha;
1496 
1497  double xoff = 0.0;
1498 
1499  for ( int ii = 0; ii < nlpts; ii++ )
1500  { // Fill the input solutes vector
1501  double efac = 0.5 * erf( ( xoff - par2 ) / p1fac ) + 0.5;
1502  xcurr = xmin + xoff * xrange;
1503  ycurr = ystr + ydif * efac;
1504  sim_vals.zsolutes << US_ZSolute( xcurr * 1e-13, ycurr, vbar20, 0.0 );
1505  xoff += xoinc;
1506  }
1507 
1508  // Evaluate the model
1509  double rmsd = evaluate_model( dsets, sim_vals );
1510 
1511  epar[ 15 ] = rmsd;
1512  int ktimms = ftimer.elapsed();
1513 qDebug() << "ffDS: call" << ffcall << "par1 par2" << par1 << par2
1514  << "rmsd" << rmsd << "eval time" << ktimms << "ms.";
1515 if(ffcall<6)
1516 qDebug() << "ffDS: epar0 epar1-9" << parP[0] << epar[1] << epar[2] << epar[3]
1517  << epar[4] << epar[5] << epar[6] << epar[7] << epar[8] << epar[9];
1518 
1519  return rmsd;
1520 }
1521 
1522 // Do curve-fit evaluate function (return RMSD) for a Horizontal Line model
1523 double US_pcsaProcess::fit_function_HL( double t, double* par )
1524 {
1525  static int ffcall=0; // Fit function call counter
1526  static double epar[ 18 ]; // Static array for holding parameters
1527  static const int nepar = sizeof( epar ) / sizeof( epar[ 0 ] );
1528 
1529  if ( t > 0.0 )
1530  { // If not t[0], return immediately
1531  return 0.0;
1532  }
1533 
1534  if ( t < 0.0 )
1535  { // If t is special flag (-ve.), reset ffcall and return
1536  ffcall = 0;
1537  return 0.0;
1538  }
1539 
1540  QTime ftimer;
1541  ftimer.start();
1542  QList< US_SolveSim::DataSet* > dsets;
1543  void** iparP = (void**)par;
1544  void** parP = (void**)epar;
1545  double par1 = par[ 0 ];
1546  double par2 = par[ 0 ];
1547  int px = ( sizeof( double ) / sizeof( void* ) ) * 2;
1548  if ( ffcall == 0 )
1549  { // On 1st call, copy par array to internal static one
1550  for ( int ii = 0; ii < nepar - 2; ii++ )
1551  epar[ ii ] = par[ ii + 2 ];
1552  parP[ 0 ] = iparP[ px ];
1553  }
1554 
1555  ffcall++; // Bump function call counter
1556  dsets << (US_SolveSim::DataSet*)parP[ 0 ];
1557  int nlpts = (int)epar[ 1 ]; // Get limit parameters
1558  double xmin = epar[ 2 ];
1559  double xmax = epar[ 3 ];
1560  double ylow = epar[ 6 ];
1561  double yhigh = epar[ 7 ];
1562  double p1lo = epar[ 8 ];
1563  double p1hi = epar[ 9 ];
1564  int noisfl = (int)epar[ 12 ];
1565  int dbg_lv = (int)epar[ 13 ];
1566  double alpha = epar[ 14 ];
1567  int attr_x = (int)epar[ 16 ];
1568  int attr_y = (int)epar[ 17 ];
1569  //int attr_z = (int)epar[ 18 ];
1570  //int zbase = epar[ 19 ];
1571  double yval = par1;
1572 
1573  // After 1st few calls, test if parameters are within limits
1574  if ( ffcall > 3 )
1575  {
1576  // Leave a little wiggle room on limits
1577  ylow -= 0.1;
1578  yhigh += 0.1;
1579  p1lo -= 0.1;
1580  p1hi += 0.1;
1581 
1582  // If this record is beyond any limit, return now with it marked as bad
1583  if ( par1 < p1lo || par1 > p1hi ||
1584  yval < ylow || yval > yhigh )
1585  {
1586 qDebug() << "ffHL: call" << ffcall << "par1" << par1 << "yv" << yval
1587  << "*OUT-OF-LIMITS*";
1588  return 1e+99;
1589  }
1590  }
1591 
1592  double prange = (double)( nlpts - 1 );
1593  double xinc = ( xmax - xmin ) / prange;
1594  double vbar20 = dsets[ 0 ]->vbar20;
1595  double xcurr = xmin;
1596  double ycurr = yval;
1597  double zcurr = dsets[ 0 ]->zcoeffs[ 0 ];
1598  zcurr = ( zcurr == 0.0 ) ? vbar20 : zcurr;
1599 qDebug() << "ffHL: zcurr vb" << zcurr << vbar20;
1600  if ( attr_x == ATTR_S )
1601  { // X is s: scale by 1e-13
1602  xcurr *= 1e-13;
1603  xinc *= 1e-13;
1604  }
1605  if ( attr_y == ATTR_S )
1606  { // Y is s: scale by 1e-13
1607  ycurr *= 1e-13;
1608  }
1609  US_SolveSim::Simulation sim_vals;
1610  sim_vals.noisflag = noisfl;
1611  sim_vals.dbg_level = dbg_lv;
1612  sim_vals.alpha = alpha;
1613 
1614  for ( int ii = 0; ii < nlpts; ii++ )
1615  { // Fill the input solutes vector
1616  sim_vals.zsolutes << US_ZSolute( xcurr, ycurr, zcurr, 0.0 );
1617  xcurr += xinc;
1618  }
1619 
1620  // Evaluate the model
1621  double rmsd = evaluate_model( dsets, sim_vals );
1622 
1623  epar[ 15 ] = rmsd;
1624  int ktimms = ftimer.elapsed();
1625 qDebug() << "ffHL: call" << ffcall << "par1 par2" << par1 << par2
1626  << "rmsd" << rmsd << "eval time" << ktimms << "ms.";
1627 if(ffcall<6)
1628 qDebug() << "ffHL: epar0 epar1-9" << parP[0] << epar[1] << epar[2] << epar[3]
1629  << epar[4] << epar[5] << epar[6] << epar[7] << epar[8] << epar[9];
1630 //qDebug() << "ffHL: dsets[0]" << dsets[0] << parP[0]
1631 // << "dsets[0]->vbar20" << dsets[0]->vbar20;
1632 qDebug() << "ffHL: yv yl yh" << yval << ylow << yhigh;
1633 
1634  return rmsd;
1635 }
1636 
1637 // Do Levenberg-Marquardt fit
1638 void US_pcsaProcess::LevMarq_fit( void )
1639 {
1640  const int eslnc = 32; // Estimated straight-line LM eval calls
1641  const int esigc = 44; // Estimated sigmoid LM eval calls
1642  static US_LM::LM_Control control( 1.e-5, 1.e-5, 1.e-5, 1.e-5,
1643  100., 100, 0, 3 );
1644  if ( lmmxcall < 1 )
1645  return;
1646 
1647  static US_LM::LM_Status status;
1648  static int n_par = 2;
1649  static int m_dat = 3;
1650 DbgLv(1) << "LMf: n_par m_dat" << n_par << m_dat;
1651  static double tarray[ 3 ] = { 0.0, 1.0, 2.0 };
1652  static double yarray[ 3 ] = { 0.0, 0.0, 0.0 };
1653  static double parray[ 20 ];
1654  bool LnType = ( curvtype == CTYPE_SL || curvtype == CTYPE_HL );
1655  double minyv = yuplim;
1656  double maxyv = ylolim;
1657 #if 0
1658  double minsl = ( yuplim - ylolim ) / ( xuplim - xlolim );
1659  double maxsl = -minsl;
1660 #endif
1661  double minp1 = LnType ? minyv : 0.5;
1662  double maxp1 = LnType ? maxyv : 0.001;
1663  double minp2 = LnType ? minyv : 1.0;
1664  double maxp2 = LnType ? maxyv : 0.0;
1665  double minp3 = 9.99e+22;
1666  double maxp3 = -9.99e+22;
1667  int npar = ( curvtype != CTYPE_HL ) ? n_par : 1;
1668  control.maxcall = lmmxcall / ( npar + 1 );
1669  lm_done = false;
1670  // Start timer for L-M progress bar, based on estimated duration
1671  kcsteps = 0;
1672  int stepms = 500;
1673  kctask = LnType ?
1674  ( time_fg * nthreads * eslnc + nmtasks / 2 ) / nmtasks :
1675  ( time_fg * nthreads * esigc + nmtasks / 2 ) / nmtasks;
1676  nctotal = ( kctask + stepms / 2 ) / stepms;
1677  kctask = nctotal * stepms;
1678  lmtm_id = startTimer( stepms );
1679  mrecs[ 0 ].variance = 0.1;
1680 
1681  QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );
1682  emit message_update( tr( "\nNow refining the best model with a "
1683  "Levenberg-Marquardt fit ...\n" ), true );
1684 
1685 #if 0
1686  elite_limits( mrecs, minyv, maxyv, minp1, maxp1, minp2, maxp2 );
1687 #endif
1688 #if 1
1689  US_ModelRecord::elite_limits( mrecs, curvtype, minyv, maxyv,
1690  minp1, maxp1, minp2, maxp2, minp3, maxp3 );
1691 #endif
1692 #if 0
1693 const double par1lo=0.001;
1694 const double par1up=0.5;
1695 const double par2lo=0.0;
1696 const double par2up=1.0;
1697 minyv=ylolim;
1698 maxyv=yuplim;
1699 minp1=(minp1+par1lo)*0.5;
1700 maxp1=(maxp1+par1up)*0.5;
1701 minp2=(minp2+par2lo)*0.5;
1702 maxp2=(maxp2+par2up)*0.5;
1703 #endif
1704 
1705  double ibm_rmsd = mrecs[ 0 ].rmsd; // Initial Best Model RMSD
1706  double *t = tarray;
1707  double *y = yarray;
1708  double *par = parray;
1709  void **ppar = (void**)par;
1710  int px = ( sizeof( double ) / sizeof( void* ) ) * 2;
1711  double zbase = dsets[ 0 ]->zcoeffs[ 0 ];
1712  par[ 0 ] = mrecs[ 0 ].par1;
1713  par[ 1 ] = mrecs[ 0 ].par2;
1714 DbgLv(0) << "LMf: par1 par2" << par[0] << par[1];
1715  par[ 2 ] = 0.0;
1716  ppar[ px ] = (void*)dsets[ 0 ];
1717  par[ 3 ] = (double)cresolu;
1718  par[ 4 ] = xlolim;
1719  par[ 5 ] = xuplim;
1720  par[ 6 ] = ylolim;
1721  par[ 7 ] = yuplim;
1722  par[ 8 ] = minyv;
1723  par[ 9 ] = maxyv;
1724  par[ 10 ] = minp1;
1725  par[ 11 ] = maxp1;
1726  par[ 12 ] = minp2;
1727  par[ 13 ] = maxp2;
1728  par[ 14 ] = noisflag;
1729  par[ 15 ] = dbg_level;
1730 // par[ 15 ] = 0;
1731  par[ 16 ] = alpha_lm;
1732  par[ 17 ] = attr_x;
1733  par[ 18 ] = attr_y;
1734  par[ 19 ] = attr_z;
1735  par[ 20 ] = zbase;
1736 DbgLv(1) << "LMf: ppar2" << ppar[4] << dsets[0] << curvtype << ppar[5];
1737 DbgLv(1) << "LMf: alpha" << alpha_lm << par[14];
1738  timer.start(); // start a timer to measure run time
1739 
1740  if ( curvtype == CTYPE_SL )
1741  { // Fit with Levenberg-Marquardt for straight-line curves
1742  control.ftol = 1.e-5;
1743  control.xtol = 1.e-5;
1744  control.gtol = 1.e-5;
1745  control.epsilon = 1.e-5;
1746 DbgLv(0) << "lmcurve_fit (SL) with par1,par2" << par[0] << par[1]
1747  << "ftol,epsl" << control.ftol << control.epsilon;
1748  fit_function_SL( -1.0, par ); // Make sure to reset eval. function
1749 
1750  US_LM::lmcurve_fit_rmsd( npar, par, m_dat, t, y,
1751  &(US_pcsaProcess::fit_function_SL), &control, &status );
1752 
1753 DbgLv(0) << " lmcurve_fit (SL) return: par1,par2" << par[0] << par[1];
1754 DbgLv(0) << " lmcfit status: fnorm nfev info"
1755  << status.fnorm << status.nfev << status.info
1756  << US_LM::lm_statmsg( &status, false );
1757 double ys = par[0];
1758 double ye = ys + par[1] * ( xuplim - xlolim );
1759 //double rmsd = par[10];
1760 double rmsd = sqrt( dsets[0]->model.variance );
1761 int nsol = dsets[0]->model.components.size();
1762 DbgLv(0) << " lmcfit xs,ys xe,ye rmsd" << xlolim << ys << xuplim << ye
1763  << rmsd << "ncsol" << nsol;
1764  }
1765 
1766  else if ( curvtype == CTYPE_IS )
1767  { // Fit with Levenberg-Marquardt for increasing-sigmoid curves
1768  //control.ftol = 30.0 * DBL_EPSILON;
1769  control.ftol = 1.0e-5;
1770  control.xtol = control.ftol;
1771  control.gtol = control.ftol;
1772  control.epsilon = ibm_rmsd * 8.0e-5;
1773  fit_function_IS( -1.0, par ); // Make sure to reset eval. function
1774 DbgLv(0) << "lmcurve_fit (IS) with par1,par2" << par[0] << par[1]
1775  << QString().sprintf( "%14.8e %14.8e", par[0], par[1] )
1776  << "ftol,epsl" << control.ftol << control.epsilon;
1777 
1778  US_LM::lmcurve_fit_rmsd( npar, par, m_dat, t, y,
1779  &(US_pcsaProcess::fit_function_IS), &control, &status );
1780 
1781 DbgLv(0) << " lmcurve_fit (IS) return: par1,par2" << par[0] << par[1]
1782  << QString().sprintf( "%14.8e %14.8e", par[0], par[1] );
1783 DbgLv(0) << " lmcfit status: fnorm nfev info"
1784  << status.fnorm << status.nfev << status.info
1785  << US_LM::lm_statmsg( &status, false );
1786 double rmsd = sqrt( dsets[0]->model.variance );
1787 int nsol = dsets[0]->model.components.size();
1788 DbgLv(0) << " lmcfit rmsd" << rmsd << "#solutes" << nsol;
1789  }
1790 
1791  else if ( curvtype == CTYPE_DS )
1792  { // Fit with Levenberg-Marquardt for decreasing-sigmoid curves
1793  control.ftol = 1.e-16;
1794  control.xtol = 1.e-16;
1795  control.gtol = 1.e-16;
1796  control.epsilon = 1.e-4;
1797  fit_function_DS( -1.0, par ); // Make sure to reset eval. function
1798 DbgLv(0) << "lmcurve_fit (DS) with par1,par2" << par[0] << par[1]
1799  << QString().sprintf( "%14.8e %14.8e", par[0], par[1] )
1800  << "ftol,epsl" << control.ftol << control.epsilon;
1801 
1802  US_LM::lmcurve_fit_rmsd( npar, par, m_dat, t, y,
1803  &(US_pcsaProcess::fit_function_DS), &control, &status );
1804 
1805 DbgLv(0) << " lmcurve_fit (DS) return: par1,par2" << par[0] << par[1]
1806  << QString().sprintf( "%14.8e %14.8e", par[0], par[1] );
1807 DbgLv(0) << " lmcfit status: fnorm nfev info"
1808  << status.fnorm << status.nfev << status.info;
1809 double rmsd = sqrt( dsets[0]->model.variance );
1810 int nsol = dsets[0]->model.components.size();
1811 DbgLv(0) << " lmcfit rmsd" << rmsd << "#solutes" << nsol;
1812  }
1813 
1814  else if ( curvtype == CTYPE_HL )
1815  { // Fit with Levenberg-Marquardt for horizontal-line curves
1816  control.ftol = 1.e-5;
1817  control.xtol = 1.e-5;
1818  control.gtol = 1.e-5;
1819  control.epsilon = 1.e-5;
1820 DbgLv(0) << "lmcurve_fit (HL) with par1" << par[0]
1821  << "ftol,epsl" << control.ftol << control.epsilon;
1822  fit_function_HL( -1.0, par ); // Make sure to reset eval. function
1823 
1824  US_LM::lmcurve_fit_rmsd( npar, par, m_dat, t, y,
1825  &(US_pcsaProcess::fit_function_HL), &control, &status );
1826 
1827 DbgLv(0) << " lmcurve_fit (HL) return: par1" << par[0];
1828 DbgLv(0) << " lmcfit status: fnorm nfev info"
1829  << status.fnorm << status.nfev << status.info
1830  << US_LM::lm_statmsg( &status, false );
1831 double yv = par[0];
1832 double rmsd = sqrt( dsets[0]->model.variance );
1833 int nsol = dsets[0]->model.components.size();
1834 DbgLv(0) << " lmcfit xs,yv xe,yv rmsd" << xlolim << yv << xuplim << yv
1835  << rmsd << "ncsol" << nsol;
1836  }
1837 
1838  else
1839  {
1840  DbgLv( 0 ) << "*ERROR* invalid curvtype" << curvtype;
1841  }
1842 
1843  lm_done = true;
1844  QApplication::restoreOverrideCursor();
1845  US_SolveSim::DataSet* dset = dsets[ 0 ];
1846  double rmsd = sqrt( dset->model.variance );
1847  int nsol = dset->model.components.size();
1848  int nfev = status.nfev;
1849  time_lm = timer.elapsed();
1850  int ktimes = ( time_lm + 500 ) / 1000;
1851  int ktimeh = ktimes / 3600;
1852  int ktimem = ( ktimes - ktimeh * 3600 ) / 60;
1853  ktimes = ktimes - ktimeh * 3600 - ktimem * 60;
1854 DbgLv(0) << " lmcfit time: " << ktimeh << "h" << ktimem
1855  << "m" << ktimes << "s";
1856 DbgLv(0) << " lmcfit LM time(ms): estimated" << kctask
1857  << "actual" << time_lm;
1858  QString fmsg = tr( "The new best model has par1 %1, par2 %2,\n"
1859  " RMSD %3, %4 solutes, %5 LM iters. " )
1860  .arg( par[ 0 ] ).arg( par[ 1 ] ).arg( rmsd ).arg( nsol ).arg( nfev );
1861  if ( ktimeh == 0 )
1862  fmsg = fmsg + tr( "(%1 min., %2 sec.)" )
1863  .arg( ktimem ).arg( ktimes );
1864  else
1865  fmsg = fmsg + tr( "(%1 hr., %2 min., %3 sec.)" )
1866  .arg( ktimeh ).arg( ktimem ).arg( ktimes );
1867  if ( alpha_lm != 0.0 )
1868  fmsg = fmsg + tr( "\nA Tikhonov regularization parameter of %1"
1869  " was used." ).arg( alpha_lm );
1870  emit message_update( fmsg, true );
1871 
1872  // Replace best model in vector and build out model more completely
1873  US_ModelRecord mrec = mrecs[ 0 ];
1874  US_SimulationParameters* spar = &dset->simparams;
1875  int jsp = 0;
1876 
1877  if ( curvtype == CTYPE_SL )
1878  {
1879  mrec.str_y = par[ 0 ];
1880  mrec.end_y = par[ 0 ] + par[ 1 ] * ( xuplim - xlolim );
1881  }
1882 
1883  else if ( curvtype == CTYPE_HL )
1884  {
1885  mrec.str_y = par[ 0 ];
1886  mrec.end_y = par[ 0 ];
1887  par[ 1 ] = 0.0;
1888  }
1889 
1890  mrec.par1 = par[ 0 ];
1891  mrec.par2 = par[ 1 ];
1892  mrec.par3 = 0.0;
1893  mrec.variance = dset->model.variance;
1894  mrec.rmsd = rmsd;
1895  mrec.ctype = curvtype;
1896  mrec.xmin = xlolim;
1897  mrec.xmax = xuplim;
1898  mrec.ymin = ylolim;
1899  mrec.ymax = yuplim;
1900 
1901  for ( int ii = 0; ii < mrec.isolutes.size(); ii++ )
1902  { // Replace x and y in top model input solutes
1903  mrec.isolutes[ ii ].x = spar->mesh_radius[ jsp++ ];
1904  mrec.isolutes[ ii ].y = spar->mesh_radius[ jsp++ ];
1905  }
1906 
1907  mrec.csolutes.clear();
1908  model = dset->model;
1909  double sfactor = 1.0 / dset->s20w_correction;
1910  double dfactor = 1.0 / dset->D20w_correction;
1911 
1912  for ( int ii = 0; ii < nsol; ii++ )
1913  {
1914  // Insert calculated solutes into top model record
1915  US_ZSolute solute;
1916 #if 0
1917  solute.x = model.components[ ii ].s;
1918  solute.y = model.components[ ii ].f_f0;
1919  solute.z = model.components[ ii ].vbar20;
1920  solute.c = model.components[ ii ].signal_concentration;
1921 #endif
1922  US_ZSolute::set_solute_values( model.components[ ii ], solute, st_mask );
1923  mrec.csolutes << solute;
1924 DbgLv(1) << "LMf: ii" << ii << "x y c" << solute.x << solute.y << solute.c;
1925 
1926  // Calculate the remainder of component values
1927  model.components[ ii ].D = 0.0;
1928  model.components[ ii ].mw = 0.0;
1929  model.components[ ii ].f = 0.0;
1930  model.calc_coefficients( model.components[ ii ] );
1931 
1932  // Convert to experiment-space for simulation below
1933  model.components[ ii ].s *= sfactor;
1934  model.components[ ii ].D *= dfactor;
1935 DbgLv(1) << "LMf: s D mw" << model.components[ii].s
1936  << model.components[ii].D << model.components[ii].mw;
1937  }
1938 
1939  // Recalculate final refined simulation and residual
1940  US_DataIO::RawData* simdat = &mrec.sim_data;
1941  US_DataIO::RawData* resids = &mrec.residuals;
1942  US_AstfemMath::initSimData( sdata, *edata, 0.0 );
1943  US_AstfemMath::initSimData( rdata, *edata, 0.0 );
1944 DbgLv(1) << "LMf:simparms: spts meni bott temp bpos"
1945  << dset->simparams.simpoints
1946  << dset->simparams.meniscus
1947  << dset->simparams.bottom
1948  << dset->simparams.temperature
1949  << dset->simparams.bottom_position;
1950 DbgLv(1) << "LMf:model: desc analys vari" << model.description
1951  << model.analysis << model.variance;
1952 
1953  US_Astfem_RSA astfem_rsa( model, dset->simparams );
1954  astfem_rsa.calculate( sdata );
1955 
1956  // Convert model back to standard space and save in model record
1957  for ( int ii = 0; ii < nsol; ii++ )
1958  {
1959  model.components[ ii ].s /= sfactor;
1960  model.components[ ii ].D /= dfactor;
1961 DbgLv(1) << "LMf: s D mw" << model.components[ii].s
1962  << model.components[ii].D << model.components[ii].mw;
1963  }
1964 
1965  mrec.model = model;
1966 
1967  // Fetch any noise saved in dset
1968  bool tino = ( ( noisflag & 1 ) != 0 );
1969  bool rino = ( ( noisflag & 2 ) != 0 );
1970  ti_noise.count = 0;
1971  ri_noise.count = 0;
1972  ti_noise.values.clear();
1973  ri_noise.values.clear();
1974  mrec.ti_noise.clear();
1975  mrec.ri_noise.clear();
1976 
1977  // Compose any noise records
1978  if ( tino )
1979  {
1980  for ( int ii = 0; ii < npoints; ii++ )
1981  ti_noise.values << spar->mesh_radius[ jsp++ ];
1982 
1983  ti_noise.minradius = edata->radius( 0 );
1984  ti_noise.maxradius = edata->radius( npoints - 1 );
1985  ti_noise.count = npoints;
1986  mrec.ti_noise = ti_noise.values;
1987 DbgLv(1) << "LMf: ti count size" << ti_noise.count << ti_noise.values.size();
1988  }
1989 
1990  if ( rino )
1991  {
1992  for ( int ii = 0; ii < nscans; ii++ )
1993  ri_noise.values << spar->mesh_radius[ jsp++ ];
1994 
1995  ri_noise.count = nscans;
1996  mrec.ri_noise = ri_noise.values;
1997 DbgLv(1) << "LMf: ri count size" << ri_noise.count << ri_noise.values.size();
1998  }
1999 
2000  // Insert new refined best model at the top of the list
2001 DbgLv(0) << "LMf:insert-new: old par1 par2" << mrecs[0].par1 << mrecs[0].par2
2002  << "new par1 par2" << mrec.par1 << mrec.par2 << "nmtasks mrec-size"
2003  << nmtasks << mrecs.size();
2004 DbgLv(0) << "LMf: old01 x,y" << mrecs[0].isolutes[0].x << mrecs[0].isolutes[0].y
2005  << mrecs[0].isolutes[1].x << mrecs[0].isolutes[1].y;
2006 DbgLv(0) << "LMf: new01 x,y" << mrec.isolutes[0].x << mrec.isolutes[0].y
2007  << mrec.isolutes[1].x << mrec.isolutes[1].y;
2008  mrecs.insert( 0, mrec );
2009 
2010  // Re-compute simulation and residuals
2011 DbgLv(0) << "LMf: tino rino" << tino << rino;
2012 DbgLv(0) << "LMf: simdat nsc npt"
2013  << simdat->scanCount() << simdat->pointCount();
2014 DbgLv(1) << "LMf: resids nsc npt"
2015  << resids->scanCount() << resids->pointCount();
2016 DbgLv(1) << "LMf: rdata nsc npt"
2017  << rdata.scanCount() << rdata.pointCount();
2018 DbgLv(1) << "LMf: sdata nsc npt"
2019  << sdata.scanCount() << sdata.pointCount();
2020  spar->mesh_radius.clear();
2021  double cvari = 0.0;
2022  double crmsd = 0.0;
2023 
2024  for ( int ss = 0; ss < nscans; ss++ )
2025  {
2026  double rnois = rino ? mrec.ri_noise[ ss ] : 0.0;
2027 
2028  for ( int rr = 0; rr < npoints; rr++ )
2029  {
2030  double tnois = tino ? mrec.ti_noise[ rr ] : 0.0;
2031  double resval = edata->value( ss, rr )
2032  - sdata. value( ss, rr ) - tnois - rnois;
2033  rdata. setValue( ss, rr, resval );
2034 
2035  simdat->setValue( ss, rr, sdata.value( ss, rr ) );
2036  resids->setValue( ss, rr, resval );
2037 if ((ss<3 && rr<3)||((ss+4)>nscans && (rr+4)>npoints)||(rr==(npoints/2)))
2038 DbgLv(1) << "LMf: ss rr" << ss << rr << "edat sdat resv"
2039  << edata->value(ss,rr) << sdata.value(ss,rr) << resval;
2040  cvari += sq( resval );
2041  }
2042  }
2043  cvari /= (double)( nscans * npoints );
2044  crmsd = sqrt( cvari );
2045  emit progress_update( cvari ); // Pass progress on to control window
2046 DbgLv(0) << "LMf: recomputed variance rmsd" << cvari << crmsd;
2047 
2048 }
2049 
2050 void US_pcsaProcess::elite_limits( QVector< US_ModelRecord >& mrecs,
2051  double& minyv, double& maxyv, double& minp1, double& maxp1,
2052  double& minp2, double& maxp2 )
2053 {
2054  const double efrac = 0.1;
2055  // Set up variables that help insure that the par1,par2 extents of elites
2056  // extend at least one step on either side of record 0's par1,par2
2057  double m0p1 = mrecs[ 0 ].par1;
2058  double m0p2 = mrecs[ 0 ].par2;
2059  double m0p1l = m0p1;
2060  double m0p1h = m0p1;
2061  double m0p2l = m0p2;
2062  double m0p2h = m0p2;
2063  bool ln_type = false;
2064 
2065  if ( curvtype == CTYPE_SL || curvtype == CTYPE_HL )
2066  { // Possibly adjust initial par1,par2 limits for lines
2067  ln_type = true;
2068  m0p1 = mrecs[ 0 ].str_y;
2069  m0p2 = mrecs[ 0 ].end_y;
2070  m0p1l = ( m0p1 > ylolim ) ? m0p1 : ( m0p1 * 1.0001 );
2071  m0p1h = ( m0p1 < yuplim ) ? m0p1 : ( m0p1 * 0.9991 );
2072  m0p2l = ( m0p2 > ylolim ) ? m0p2 : ( m0p2 * 1.0001 );
2073  m0p2h = ( m0p2 < yuplim ) ? m0p2 : ( m0p2 * 0.9991 );
2074  }
2075 
2076  if ( curvtype == CTYPE_IS || curvtype == CTYPE_DS )
2077  { // Possibly adjust initial par1,par2 limits for sigmoids
2078  double dif1 = m0p1 - 0.001;
2079  double dif2 = m0p1 - 0.500;
2080  double dif3 = m0p2 - 0.000;
2081  double dif4 = m0p2 - 1.000;
2082  m0p1l = ( dif1 > 1.e-8 ) ? m0p1 : 0.002;
2083  m0p1h = ( dif2 < -1e-8 ) ? m0p1 : 0.499;
2084  m0p2l = ( dif3 > 1.e-8 ) ? m0p2 : 0.001;
2085  m0p2h = ( dif4 < -1e-8 ) ? m0p2 : 0.999;
2086 DbgLv(1) << " ElLim: ADJUST SIGM: m0p1 m0p1h" << m0p1 << m0p1h
2087  << "m0p1<0.500" << (m0p1<0.500) << 0.500 << "(m0p1-0.5)" << (m0p1-0.5);
2088  }
2089 
2090  int nmr = mrecs.size();
2091  int nelite = qRound( efrac * nmr ); // Elite is top 10%
2092  int maxel = nmr / 2;
2093  int minel = qMin( maxel, 4 );
2094  nelite = qMin( nelite, maxel ); // At most half of all
2095  nelite = qMax( nelite, minel ); // At least 4
2096  nelite -= 2; // Less 2 for compare
2097 DbgLv(0) << " ElLim: nmr nelite nmtasks" << nmr << nelite << nmtasks;
2098 DbgLv(1) << " ElLim: in minyv maxyv" << minyv << maxyv;
2099 DbgLv(1) << " ElLim: in min/max p1/p2" << minp1 << maxp1 << minp2 << maxp2;
2100 DbgLv(1) << " ElLim: in m0p1 m0p2" << m0p1 << m0p2;
2101 DbgLv(1) << " ElLim: in m0p1l,m0p1h,m0p2l,m0p2h" << m0p1l << m0p1h
2102  << m0p2l << m0p2h;
2103 
2104  for ( int ii = 0; ii < nmr; ii++ )
2105  {
2106  double str_y = mrecs[ ii ].str_y;
2107  double end_y = mrecs[ ii ].end_y;
2108  double par1 = ln_type ? str_y : mrecs[ ii ].par1;
2109  double par2 = ln_type ? end_y : mrecs[ ii ].par2;
2110 if(ii<3||(ii+4)>nelite)
2111 DbgLv(1) << " ElLim: ii" << ii << "par1 par2" << par1 << par2
2112  << "str_y end_y" << str_y << end_y << "rmsd" << mrecs[ii].rmsd;
2113  minyv = qMin( minyv, str_y );
2114  maxyv = qMax( maxyv, str_y );
2115  minyv = qMin( minyv, end_y );
2116  maxyv = qMax( maxyv, end_y );
2117  minp1 = qMin( minp1, par1 );
2118  maxp1 = qMax( maxp1, par1 );
2119  minp2 = qMin( minp2, par2 );
2120  maxp2 = qMax( maxp2, par2 );
2121 
2122  // We want to break out of the min,max scan loop if the sorted index
2123  // exceeds the elite count. But we continue in the loop if we have not
2124  // yet found min,max par1,par2 values that are at least a step
2125  // on either side of the par1,par2 values for the best model (m0).
2126 if(ii>nelite)
2127 DbgLv(1) << " ElLim: minp1 maxp1 m0p1" << minp1 << maxp1 << m0p1
2128  << "minp2 maxp2 m0p2" << minp2 << maxp2 << m0p2;
2129  if ( ii > nelite &&
2130  minp1 < m0p1l && maxp1 > m0p1h &&
2131  minp2 < m0p2l && maxp2 > m0p2h )
2132  break;
2133  }
2134 
2135 DbgLv(0) << " ElLim: out minyv maxyv" << minyv << maxyv;
2136 DbgLv(0) << " ElLim: out min/max p1/p2" << minp1 << maxp1 << minp2 << maxp2;
2137 }
2138 
2139 // Evaluate a model; return rmsd, model, noises
2140 double US_pcsaProcess::evaluate_model( QList< US_SolveSim::DataSet* >& dsets,
2141  US_SolveSim::Simulation& sim_vals )
2142 {
2143  US_SolveSim::DataSet* dset = dsets[ 0 ];
2144  US_SimulationParameters* spar = &dset->simparams;
2145  spar->mesh_radius.clear();
2146 
2147  // Save the s and k of the input solutes
2148  for ( int ii = 0; ii < sim_vals.zsolutes.size(); ii++ )
2149  {
2150  spar->mesh_radius << sim_vals.zsolutes[ ii ].x;
2151  spar->mesh_radius << sim_vals.zsolutes[ ii ].y;
2152  }
2153 
2154  // Do astfem fit, mostly to get an RMSD
2155  US_SolveSim* solvesim = new US_SolveSim( dsets, 0, false );
2156  solvesim->calc_residuals( 0, 1, sim_vals );
2157 
2158  // Construct a rudimentary model from computed solutes and save it
2159  dset->model = US_Model();
2160  dset->model.variance = sim_vals.variance;
2161  dset->model.components.clear();
2162  int st_mask = dsets[ 0 ]->solute_type;
2163 
2164  for ( int ii = 0; ii < sim_vals.zsolutes.size(); ii++ )
2165  {
2166  US_Model::SimulationComponent mcomp;
2167 #if 0
2168  mcomp.s = sim_vals.zsolutes[ ii ].x;
2169  mcomp.f_f0 = sim_vals.zsolutes[ ii ].y;
2170  mcomp.vbar20 = sim_vals.zsolutes[ ii ].z;
2171  mcomp.signal_concentration = sim_vals.zsolutes[ ii ].c;
2172 #endif
2173  US_ZSolute::set_mcomp_values( mcomp, sim_vals.zsolutes[ ii ],
2174  st_mask, true );
2175 
2176  dset->model.components << mcomp;
2177  }
2178 
2179  // If noise was computed, save it in simparams mesh_radius vector
2180  if ( sim_vals.noisflag != 0 )
2181  {
2182  int ntin = ( ( sim_vals.noisflag & 1 ) == 0 ) ? 0
2183  : sim_vals.ti_noise.size();
2184  int nrin = ( ( sim_vals.noisflag & 2 ) == 0 ) ? 0
2185  : sim_vals.ri_noise.size();
2186 
2187  for ( int ii = 0; ii < ntin; ii++ )
2188  spar->mesh_radius << sim_vals.ti_noise[ ii ];
2189 
2190  for ( int ii = 0; ii < nrin; ii++ )
2191  spar->mesh_radius << sim_vals.ri_noise[ ii ];
2192  }
2193 
2194  // Compute the RMSD and return it
2195  double rmsd = sqrt( sim_vals.variance );
2196  return rmsd;
2197 }
2198 
2199 // Protected slot to filter timer event and handle L-M status timing
2200 void US_pcsaProcess::timerEvent( QTimerEvent *event )
2201 {
2202  int tm_id = event->timerId();
2203 
2204  if ( tm_id != lmtm_id )
2205  { // If other than L-M timing event, pass it on to the normal handler
2206  QObject::timerEvent( event );
2207  return;
2208  }
2209 
2210  // Otherwise, bump progress counter and emit a progress signal
2211  kcsteps++;
2212 
2213  if ( lm_done )
2214  { // If L-M is done, signal that with progress status
2215  nctotal = kcsteps;
2216  emit stage_complete( kcsteps, nctotal );
2217  killTimer( tm_id );
2218  }
2219 
2220  else if ( kcsteps < nctotal )
2221  { // Normal less-than-100% progress
2222  emit stage_complete( kcsteps, nctotal );
2223  }
2224 
2225  else
2226  { // Not done with L-M but task count met, so adjust total
2227  nctotal = ( kcsteps * 12 ) / 10;
2228  emit stage_complete( kcsteps, nctotal );
2229  }
2230 
2231  return;
2232 }
2233 
2234 // Perform one final regularization computation when L-M was non-regularized
2235 void US_pcsaProcess::compute_final()
2236 {
2237 DbgLv(1) << "CFin: alpha" << alpha << "mrecs size" << mrecs.size();
2238  QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );
2239  QTime ftimer;
2240  ftimer.start();
2241  US_ModelRecord mrec = mrecs[ 0 ];
2242  model = mrec.model;
2243  US_SolveSim::Simulation sim_vals;
2244  sim_vals.noisflag = noisflag;
2245  sim_vals.dbg_level = dbg_level;
2246  sim_vals.alpha = alpha;
2247  sim_vals.zsolutes = mrec.isolutes;
2248 DbgLv(1) << "CFin: isolutes size" << sim_vals.zsolutes.size();
2249 DbgLv(1) << "CFin: zsol0 x y z"
2250  << sim_vals.zsolutes[0].x << sim_vals.zsolutes[0].y << sim_vals.zsolutes[0].z;
2251 
2252  // Evaluate the model
2253  US_SolveSim::DataSet* dset = dsets[ 0 ];
2254 
2255  double rmsd = evaluate_model( dsets, sim_vals );
2256 
2257  sdata = sim_vals.sim_data;
2258  rdata = sim_vals.residuals;
2259  int nmsol = dset->model.components.size();
2260  int nisol = mrec.isolutes.size();
2261  int ktimsc = ( ftimer.elapsed() + 500 ) / 1000;
2262 DbgLv(1) << "CFin: nmsol nisol" << nmsol << nisol;
2263 
2264  QString fmsg = tr(
2265  "\nA final best model (RMSD=%1; %2-solute, %3 out; %4 sec.)\n"
2266  " used a Tikhonov regularization parameter of %5 .\n" )
2267  .arg( rmsd ).arg( nisol ).arg( nmsol ).arg( ktimsc ).arg( alpha );
2268  emit message_update( fmsg, true );
2269 
2270  // Replace best model in vector and build out model more completely
2271  mrec.variance = dset->model.variance;
2272  mrec.rmsd = rmsd;
2273 
2274  mrec.csolutes.clear();
2275  model = dset->model;
2276 
2277  for ( int ii = 0; ii < nmsol; ii++ )
2278  {
2279  // Insert calculated solutes into top model record
2280  US_ZSolute solute;
2281 #if 0
2282  solute.x = model.components[ ii ].s;
2283  solute.y = model.components[ ii ].f_f0;
2284  solute.z = model.components[ ii ].vbar20;
2285  solute.c = model.components[ ii ].signal_concentration;
2286 #endif
2287  US_ZSolute::set_solute_values( model.components[ ii ], solute, st_mask );
2288  mrec.csolutes << solute;
2289 DbgLv(1) << "CFin: ii" << ii << "x y c" << solute.x << solute.y << solute.c;
2290 
2291  // Calculate the remainder of component values
2292  model.components[ ii ].D = 0.0;
2293  model.components[ ii ].mw = 0.0;
2294  model.components[ ii ].f = 0.0;
2295  model.calc_coefficients( model.components[ ii ] );
2296 
2297 DbgLv(1) << "CFin: s D mw" << model.components[ii].s
2298  << model.components[ii].D << model.components[ii].mw;
2299  }
2300 
2301 DbgLv(1) << "CFin:model: desc analys vari" << model.description
2302  << model.analysis << model.variance;
2303 
2304  // Replace the top model with the new regularized best model
2305  model.alphaRP = alpha;
2306  mrec.model = model;
2307  mrecs[ 0 ] = mrec;
2308 
2309  // Report new variance
2310  emit progress_update( mrec.variance );
2311  QApplication::restoreOverrideCursor();
2312 DbgLv(0) << "CFin: recomputed variance rmsd" << mrec.variance << rmsd;
2313 }
2314 
2315 // Restart the curve grid iteration sequence
2316 void US_pcsaProcess::restart_fit()
2317 {
2318 #if 0
2319  bool LnType = ( curvtype == CTYPE_SL || curvtype == CTYPE_HL );
2320  bool SgType = ( curvtype == CTYPE_IS || curvtype == CTYPE_DS );
2321 #endif
2322  errMsg = tr( "NO ERROR: start" );
2323  maxrss = 0;
2324  varimin = 9.e+9;
2325  minvarx = 99999;
2326 DbgLv(1) << "RF: nmr" << mrecs.size() << "cfi_rmsd" << cfi_rmsd;
2327 
2328  wkstates .resize( nthreads );
2329  wthreads .clear();
2330  job_queue.clear();
2331 
2332  orig_sols.clear();
2333 
2334 DbgLv(1) << "RF: sll sul" << xlolim << xuplim
2335  << " yll yul nyp" << ylolim << yuplim << nypts
2336  << " type reso nth noi" << curvtype << cresolu << nthreads << noisflag;
2337 
2338  fi_iter++;
2339  double elite_fact = 0.1;
2340 
2341  if ( curvtype == CTYPE_2O )
2342  {
2343  elite_fact = (double)qMax( 1, ( 4 - fi_iter ) ) * 0.1;
2344  }
2345 
2346  mrecs[ 0 ].variance = elite_fact;
2347 
2348  US_ModelRecord::recompute_mrecs( curvtype, xlolim, xuplim, ylolim, yuplim,
2349  nypts, cresolu, parlims, mrecs );
2350 
2351  // Update the solutes list from the newly recomputed model records
2352  for ( int ii = 0; ii < mrecs.size(); ii++ )
2353  {
2354  orig_sols << mrecs[ ii ].isolutes;
2355  }
2356 DbgLv(1) << "SGMO: orig_sols size" << orig_sols.size()
2357  << "nmodels" << mrecs.size();
2358 
2359  if ( curvtype == CTYPE_2O )
2360  {
2361  mrecs[ 0 ].ymin = xlolim;
2362  mrecs[ 0 ].ymax = xuplim;
2363  }
2364 
2365  kcsteps = 0;
2366  nctotal = nmtasks;
2367  emit stage_complete( kcsteps, nctotal );
2368 
2369  kctask = 0;
2370  kstask = 0;
2371  nthreads = ( nthreads < nmtasks ) ? nthreads : nmtasks;
2372 DbgLv(1) << "RF: nscans npoints" << nscans << npoints;
2373 DbgLv(1) << "RF: nctotal nthreads" << nctotal << nthreads;
2374  max_rss();
2375 DbgLv(1) << "RF: (1)maxrss" << maxrss;
2376 
2377  // Queue all the tasks
2378  for ( int ktask = 0; ktask < nmtasks; ktask++ )
2379  {
2380  WorkPacketPc wtask;
2381  int mm = orig_sols[ ktask ].size() - 1;
2382  double stry = orig_sols[ ktask ][ 0 ].y;
2383  double endy = orig_sols[ ktask ][ mm ].y;
2384  wtask.par1 = mrecs[ ktask ].par1;
2385  wtask.par2 = mrecs[ ktask ].par2;
2386  wtask.par3 = mrecs[ ktask ].par3;
2387  wtask.depth = 0;
2388 
2389  wtask.sim_vals.alpha = alpha_fx;
2390 
2391  queue_task( wtask, stry, endy, ktask, noisflag, orig_sols[ ktask ] );
2392  }
2393 
2394  // Start the first threads. This will begin the first work units (subgrids).
2395  // Thereafter, work units are started in new threads when previous
2396  // threads signal that they have completed their work.
2397 
2398  for ( int ii = 0; ii < nthreads; ii++ )
2399  {
2400  wthreads << 0;
2401 
2402  WorkPacketPc wtask = job_queue.takeFirst();
2403  submit_job( wtask, ii );
2404  int mempav = US_Memory::memory_profile();
2405 DbgLv(1) << "PC:MEM: (5)PcAvail" << mempav;
2406  }
2407 
2408  mrecs .clear();
2409  max_rss();
2410  kstask = nthreads; // count of started tasks is initially thread count
2411 DbgLv(1) << "RF: kstask nthreads" << kstask << nthreads << job_queue.size();
2412  int mempav = US_Memory::memory_profile();
2413 DbgLv(1) << "PC:MEM: (6)PcAvail" << mempav;
2414 
2415  emit message_update( pmessage_head() +
2416  tr( "Starting fit iteration %1 computations of %2 models,\n"
2417  " using %3 threads" )
2418  .arg( fi_iter ).arg( nmtasks ).arg( nthreads ), false );
2419 }
2420 
2421 // Clear memory allocation
2422 void US_pcsaProcess::clear_memory()
2423 {
2424 int memmb = qRound( (double)US_Memory::rss_now() / 1024.0 );
2425 DbgLv(1) << "PC:MEM: ClrMem IN rssnow" << memmb;
2426  wkstates .clear();
2427  wthreads .clear();
2428  job_queue.clear();
2429 
2430  orig_sols.clear();
2431  mrecs .clear();
2432 memmb = qRound( (double)US_Memory::rss_now() / 1024.0 );
2433 DbgLv(1) << "PC:MEM: ClrMem OUT rssnow" << memmb;
2434 }
Generated on Sun Jun 21 2015 11:40:42 for UltraScan III by   doxygen 1.8.3.1-20130324