UltraScan III
us_rpscan.cpp
Go to the documentation of this file.
1 
3 #include "us_rpscan.h"
4 #include "us_pcsa.h"
5 #include "us_settings.h"
6 #include "us_gui_settings.h"
7 #include "us_math2.h"
8 #include "us_sleep.h"
9 
10 #include <qwt_double_interval.h>
11 #include <qwt_scale_widget.h>
12 #include <qwt_scale_draw.h>
13 
14 // Constructor: Regularization Parameter Scan widget
15 US_RpScan::US_RpScan( QList< US_SolveSim::DataSet* >&dsets,
16  US_ModelRecord& mr, int& thr, double& alf, QWidget* p )
17  : US_WidgetsDialog( p, 0 ), dsets( dsets ), mrec( mr ), nthr( thr ),
18  alpha( alf )
19 {
20  alpha = 0.0;
21 
22 DbgLv(1) << "TRP: RpScan IN";
23  // lay out the GUI
24  setObjectName( "US_RpScan" );
25  setAttribute( Qt::WA_DeleteOnClose, false );
26  setWindowTitle( tr( "Regularization Parameter (Alpha) Scan" ) );
27  setPalette( US_GuiSettings::frameColor() );
28 
29  QSize p1size( 560, 480 );
30 
31  dbg_level = US_Settings::us_debug();
32  v_line = NULL;
33 
34  mainLayout = new QHBoxLayout( this );
35  leftLayout = new QVBoxLayout();
36  rightLayout = new QVBoxLayout();
37  pltctrlsLayout = new QGridLayout();
38  buttonsLayout = new QHBoxLayout();
39 
40  mainLayout->setSpacing ( 2 );
41  mainLayout->setContentsMargins( 2, 2, 2, 2 );
42 
43  QLabel* lb_datctrls = us_banner( tr( "PCSA Model Parameters" ) );
44  QLabel* lb_mtype = us_label( tr( "Model Type:" ) );
45  QLabel* lb_npoints = us_label( tr( "Points per Line:" ) );
46  QLabel* lb_mdlpar1 = us_label( tr( "Best Model Par 1:" ) );
47  QLabel* lb_mdlpar2 = us_label( tr( "Best Model Par 2:" ) );
48  QLabel* lb_scnctrls = us_banner( tr( "Alpha Scan Parameters" ) );
49  QLabel* lb_stralpha = us_label( tr( "Starting Alpha:" ) );
50  QLabel* lb_endalpha = us_label( tr( "Ending Alpha:" ) );
51  QLabel* lb_incalpha = us_label( tr( "Alpha Increment:" ) );
52  QLabel* lb_selalpha = us_label( tr( "Selected Alpha:" ) );
53  QLabel* lb_status = us_label( tr( "Status:" ) );
54 
55  QPushButton* pb_scan = us_pushbutton( tr( "Start Scan" ) );
56  QPushButton* pb_help = us_pushbutton( tr( "Help" ) );
57  QPushButton* pb_accept = us_pushbutton( tr( "Accept" ) );
58  QPushButton* pb_cancel = us_pushbutton( tr( "Cancel" ) );
59 
60  int nppln = mrec.isolutes.size();
61 DbgLv(1) << "TRP: nppln" << nppln;
62  // lay out the GUI
63  double par1 = mrec.par1;
64  double par2 = mrec.par2;
65  le_mtype = us_lineedit( tr( "Straight Line" ), -1, true );
66  le_npoints = us_lineedit( QString::number( nppln ), -1, true );
67  le_mdlpar1 = us_lineedit( QString::number( par1 ), -1, true );
68  le_mdlpar2 = us_lineedit( QString::number( par2 ), -1, true );
69  ct_stralpha = us_counter( 3, 0.000, 1000.0, 0.020 );
70  ct_endalpha = us_counter( 3, 0.000, 1000.0, 0.800 );
71  ct_incalpha = us_counter( 3, 0.000, 1000.0, 0.020 );
72  le_selalpha = us_lineedit( QString::number( alpha ), -1, false );
73  ct_stralpha->setStep( 0.001 );
74  ct_endalpha->setStep( 0.001 );
75  ct_incalpha->setStep( 0.001 );
76 
77  b_progress = us_progressBar( 0, 100, 0 );
78  le_stattext = us_lineedit( "", -1, true );
79 
80  // Adjust the size of line counters and rmsd text
81  lb_mtype ->adjustSize();
82  QFont font( US_GuiSettings::fontFamily(),
83  US_GuiSettings::fontSize() );
84  QFontMetrics fmet( font );
85  int fwidth = fmet.maxWidth();
86  int rheight = lb_mtype->height();
87  int cminw = fwidth * 9;
88  int csizw = cminw + fwidth;
89  ct_stralpha->resize( csizw, rheight );
90  ct_endalpha->resize( csizw, rheight );
91  ct_incalpha->resize( csizw, rheight );
92  ct_stralpha->setMinimumWidth( cminw );
93  ct_endalpha->setMinimumWidth( cminw );
94  ct_incalpha->setMinimumWidth( cminw );
95 DbgLv(1) << "TRP: csizw cminw" << csizw << cminw;
96 
97  // Add elements to the controls layout
98  int row = 0;
99  pltctrlsLayout->addWidget( lb_datctrls, row++, 0, 1, 6 );
100  pltctrlsLayout->addWidget( lb_mtype, row, 0, 1, 3 );
101  pltctrlsLayout->addWidget( le_mtype, row++, 3, 1, 3 );
102  pltctrlsLayout->addWidget( lb_npoints, row, 0, 1, 3 );
103  pltctrlsLayout->addWidget( le_npoints, row++, 3, 1, 3 );
104  pltctrlsLayout->addWidget( lb_mdlpar1, row, 0, 1, 3 );
105  pltctrlsLayout->addWidget( le_mdlpar1, row++, 3, 1, 3 );
106  pltctrlsLayout->addWidget( lb_mdlpar2, row, 0, 1, 3 );
107  pltctrlsLayout->addWidget( le_mdlpar2, row++, 3, 1, 3 );
108  pltctrlsLayout->addWidget( lb_scnctrls, row++, 0, 1, 6 );
109  pltctrlsLayout->addWidget( lb_stralpha, row, 0, 1, 3 );
110  pltctrlsLayout->addWidget( ct_stralpha, row++, 3, 1, 3 );
111  pltctrlsLayout->addWidget( lb_endalpha, row, 0, 1, 3 );
112  pltctrlsLayout->addWidget( ct_endalpha, row++, 3, 1, 3 );
113  pltctrlsLayout->addWidget( lb_incalpha, row, 0, 1, 3 );
114  pltctrlsLayout->addWidget( ct_incalpha, row++, 3, 1, 3 );
115  pltctrlsLayout->addWidget( lb_selalpha, row, 0, 1, 3 );
116  pltctrlsLayout->addWidget( le_selalpha, row, 3, 1, 2 );
117  pltctrlsLayout->addWidget( pb_scan, row++, 5, 1, 1 );
118  pltctrlsLayout->addWidget( lb_status, row, 0, 1, 1 );
119  pltctrlsLayout->addWidget( b_progress, row++, 1, 1, 5 );
120  pltctrlsLayout->addWidget( le_stattext, row++, 0, 1, 6 );
121 // row += 7;
122  QString ctype = tr( "Increasing Sigmoid" );
123  if ( mrec.str_y == mrec.par1 )
124  {
125  ctype = mrec.str_y != mrec.end_y
126  ? tr( "Straight Line" )
127  : tr( "Horizontal Line" );
128  }
129  else if ( mrec.str_y > mrec.end_y )
130  ctype = tr( "Decreasing Sigmoid" );
131  le_mtype->setText( ctype );
132 
133  buttonsLayout ->addWidget( pb_help );
134  buttonsLayout ->addWidget( pb_cancel );
135  buttonsLayout ->addWidget( pb_accept );
136 
137  // Complete layouts and set up signals/slots
138  plotLayout1 = new US_Plot( data_plot1,
139  tr( "Alpha Scan Points" ),
140  tr( "Variance (x 1e-5)" ),
141  tr( "Norm of X (solute concentrations)" ) );
142 
143  data_plot1->setCanvasBackground( Qt::black );
144  data_plot1->setMinimumSize( p1size );
145  data_plot1->enableAxis( QwtPlot::xBottom, true );
146  data_plot1->enableAxis( QwtPlot::yLeft, true );
147 
148  pick = new US_PlotPicker( data_plot1 );
149  pick->setRubberBand ( QwtPicker::CrossRubberBand );
150  pick->setMousePattern( QwtEventPattern::MouseSelect1,
151  Qt::LeftButton );
152 
153  rightLayout->addLayout( plotLayout1 );
154 
155  leftLayout ->addLayout( pltctrlsLayout );
156  leftLayout ->addStretch();
157  leftLayout ->addLayout( buttonsLayout );
158 
159  mainLayout->addLayout( leftLayout );
160  mainLayout->addLayout( rightLayout );
161  mainLayout->setStretchFactor( leftLayout, 3 );
162  mainLayout->setStretchFactor( rightLayout, 5 );
163 
164  connect( pb_scan, SIGNAL( clicked() ),
165  this, SLOT ( scan() ) );
166  connect( pb_help, SIGNAL( clicked() ),
167  this, SLOT ( help() ) );
168  connect( pb_cancel, SIGNAL( clicked() ),
169  this, SLOT ( reject_it() ) );
170  connect( pb_accept, SIGNAL( clicked() ),
171  this, SLOT ( accept_it() ) );
172 
173 DbgLv(1) << "TRP: p1size" << p1size;
174  data_plot1 ->resize( p1size );
175  ct_stralpha->resize( csizw, rheight );
176  ct_endalpha->resize( csizw, rheight );
177  ct_incalpha->resize( csizw, rheight );
178  adjustSize();
179 DbgLv(1) << "TRP: actual csizw" << ct_incalpha->width();
180 DbgLv(1) << "TRP: mrec.taskx" << mrec.taskx;
181 DbgLv(1) << "TRP: mrec.str_y" << mrec.str_y;
182 DbgLv(1) << "TRP: mrec.end_y" << mrec.end_y;
183 DbgLv(1) << "TRP: mrec.par1 " << mrec.par1;
184 DbgLv(1) << "TRP: mrec.par2 " << mrec.par2;
185 DbgLv(1) << "TRP: mrec.vari " << mrec.variance;
186 DbgLv(1) << "TRP: mrec.rmsd " << mrec.rmsd;
187 DbgLv(1) << "TRP: mrec.isolutes.size" << mrec.isolutes.size();
188 DbgLv(1) << "TRP: mrec.csolutes.size" << mrec.csolutes.size();
189  qApp->processEvents();
190 }
191 
192 // Cancel button clicked
193 void US_RpScan::reject_it()
194 {
195  reject();
196  close();
197 }
198 
199 // Accept button clicked
200 void US_RpScan::accept_it()
201 {
202  accept();
203  alpha = le_selalpha->text().toDouble();
204  close();
205 }
206 
207 // Scan alphas
208 void US_RpScan::scan()
209 {
210  QTime timer;
211  int time_sc;
212  double calpha = ct_stralpha->value();
213  double ealpha = ct_endalpha->value();
214  double dalpha = ct_incalpha->value();
215  nalpha = qRound( ( ealpha - calpha ) / dalpha ) + 1;
216  alphas.clear();
217  varias.clear();
218  xnorms.clear();
219  b_progress->reset();
220  b_progress->setMaximum( nalpha );
221  double varmx = 0.0;
222  double xnomx = 0.0;
223  double v_vari = 0.0;
224  double v_xnsq = 0.0;
225  QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );
226  le_stattext->setText( tr( "Beginning Alpha Scan ..." ) );
227  timer.start(); // Start a timer to measure run time
228  qApp->processEvents();
229 
230  for ( int ja = 0; ja < nalpha; ja++ )
231  { // Loop to populate alpha scan arrays with preliminary values
232  alphas << calpha;
233  varias << v_vari;
234  xnorms << v_xnsq;
235  calpha += dalpha;
236  }
237 
238 int kalpha=nalpha;
239  nalpha = alphas.size();
240 DbgLv(1) << "ASC: nalpha" << nalpha << kalpha << "nthr" << nthr;
241 
242  if ( nthr == 1 )
243  { // Non-multi-threaded scan
244  QVector< double > csolutes;
245  sv_nnls_a.clear();
246  sv_nnls_b.clear();
247  int nscans = dsets[ 0 ]->run_data.scanCount();
248  int npoints = dsets[ 0 ]->run_data.pointCount();
249  int nisols = mrec.isolutes.size();
250  le_stattext->setText( tr( "Beginning Single-Thread Alpha Scan ..." ) );
251 
252  for ( int ja = 0; ja < nalpha; ja++ )
253  { // Loop to evaluate each alpha in the specified range
254  calpha = alphas[ ja ];
255 
256  if ( ja == 0 )
257  { // For the first one do a full compute and save the A,B matrices
258  US_SolveSim::Simulation sim_vals;
259  sim_vals.alpha = calpha;
260  sim_vals.noisflag = 0;
261  sim_vals.dbg_level = 0;
262  sim_vals.zsolutes = mrec.isolutes;
263 
264  US_SolveSim* solvesim = new US_SolveSim( dsets, 0, false );
265 
266  solvesim->calc_residuals( 0, 1, sim_vals, true,
267  &sv_nnls_a, &sv_nnls_b );
268 
269  v_vari = sim_vals.variance;
270  v_xnsq = sim_vals.xnormsq;
271  }
272 
273  else
274  { // After 1st, regularize by modifying A matrix diagonal, then NNLS
275  apply_alpha( calpha, &sv_nnls_a, &sv_nnls_b,
276  nscans, npoints, nisols, v_vari, v_xnsq );
277  }
278 
279  varias[ ja ] = v_vari;
280  xnorms[ ja ] = v_xnsq;
281  int kalf = ja + 1;
282 
283  b_progress->setValue( kalf );
284 
285  varmx = qMax( varmx, v_vari );
286  xnomx = qMax( xnomx, v_xnsq );
287  QString astat = tr( "Of %1 models, computed %2 (Alpha %3)" )
288  .arg( nalpha ).arg( kalf ).arg( calpha );
289  le_stattext->setText( astat );
290 DbgLv(1) << "a v x" << calpha << v_vari << v_xnsq;
291  qApp->processEvents();
292  }
293 
294  sv_nnls_a.clear();
295  sv_nnls_b.clear();
296  }
297 
298  else
299  { // Multi-threaded alpha scan
300  le_stattext->setText( tr( "Beginning %1-Thread Alpha Scan ..." )
301  .arg( nthr ) );
302  qApp->processEvents();
303  wthreads .clear();
304  sv_nnls_a.clear();
305  sv_nnls_b.clear();
306 
307  // Run a complete model computation to get the A,B matrices
308  US_SolveSim::Simulation sim_vals;
309  calpha = alphas[ 0 ];
310  sim_vals.alpha = calpha;
311  sim_vals.noisflag = 0;
312  sim_vals.dbg_level = 0;
313  sim_vals.zsolutes = mrec.isolutes;
314 
315  US_SolveSim* solvesim = new US_SolveSim( dsets, 0, false );
316 
317  solvesim->calc_residuals( 0, 1, sim_vals, true, &sv_nnls_a, &sv_nnls_b );
318 
319  v_vari = sim_vals.variance;
320  v_xnsq = sim_vals.xnormsq;
321 qDebug() << "A0: alpha" << calpha << "vari xnsq" << v_vari << v_xnsq
322  << "ncsols" << sim_vals.zsolutes.size();;
323  varias[ 0 ] = v_vari;
324  xnorms[ 0 ] = v_xnsq;
325  nasubm = 1;
326  nacomp = 1;
327 
328  QString astat =
329  tr( "Of %1 models, %2 done (Alpha %3, thread %4)" )
330  .arg( nalpha ).arg( nacomp ).arg( alphas[ 0 ] ).arg( 0 );
331  le_stattext->setText( astat );
332  b_progress->setValue( 1 );
333 
334  for ( int jt = 0; jt < nthr; jt++ )
335  { // Submit the first tasks using available threads
336  WorkPacketPc wtask;
337  calpha = alphas[ jt + 1 ];
338  US_SolveSim::Simulation sim_vals;
339  sim_vals.alpha = calpha;
340  sim_vals.noisflag = 0;
341  sim_vals.dbg_level = qMax( 0, ( dbg_level - 1 ) );
342  sim_vals.zsolutes = mrec.isolutes;
343  wtask.thrn = jt + 1;
344  wtask.taskx = nasubm;
345  wtask.str_y = mrec.str_y;
346  wtask.end_y = mrec.end_y;
347  wtask.par1 = mrec.par1;
348  wtask.par2 = mrec.par2;
349  wtask.sim_vals = sim_vals;
350  wtask.isolutes = sim_vals.zsolutes;
351  wtask.dsets = dsets;
352  wtask.depth = 1;
353  wtask.state = 0;
354  wtask.noisf = sim_vals.noisflag;
355  wtask.psv_nnls_a = &sv_nnls_a;
356  wtask.psv_nnls_b = &sv_nnls_b;
357  wtask.csolutes.clear();
358 DbgLv(1) << "ASC: jt" << jt << "alpha" << calpha;
359 
360  WorkerThreadPc* wthr = new WorkerThreadPc( this );
361 
362  wthr->define_work( wtask );
363  wthreads << wthr;
364  connect( wthr, SIGNAL( work_complete( WorkerThreadPc* ) ),
365  this, SLOT( process_job( WorkerThreadPc* ) ) );
366  wthr->start();
367  nasubm++;
368 DbgLv(1) << "ASC: defined: nasubm" << nasubm;
369  qApp->processEvents();
370  }
371 
372  // Keep testing to see if all tasks are complete
373  while( nacomp < nalpha )
374  {
375  US_Sleep::msleep( 500 );
376  qApp->processEvents();
377  }
378 
379  // Once all threads are done, determine max variance,norm
380  for ( int ja = 0; ja < nalpha; ja++ )
381  {
382  varmx = qMax( varmx, varias[ ja ] );
383  xnomx = qMax( xnomx, xnorms[ ja ] );
384  }
385 
386  // Delete worker threads and free saved A,B matrices
387  for ( int jt = 0; jt < nthr; jt++ )
388  delete wthreads[ jt ];
389 
390  sv_nnls_a.clear();
391  sv_nnls_b.clear();
392  }
393 
394  lgv = 0 -(int)qFloor( log10( varmx ) );
395  lgx = -1 -(int)qFloor( log10( xnomx ) );
396  vscl = qPow( 10.0, lgv );
397  xscl = qPow( 10.0, lgx );
398 DbgLv(1) << "Log-varia Log-xnorm" << lgv << lgx << "vscl xscl" << vscl << xscl;
399 
400  for ( int ja = 0; ja < nalpha; ja++ )
401  { // Scale the variance,normsq points
402  varias[ ja ] *= vscl;
403  xnorms[ ja ] *= xscl;
404  }
405 
406  b_progress->setValue( nalpha );
407  QApplication::restoreOverrideCursor();
408  qApp->processEvents();
409  // Determine elapsed time
410  time_sc = timer.elapsed();
411  int ktimes = ( time_sc + 500 ) / 1000;
412  int ktimeh = ktimes / 3600;
413  int ktimem = ( ktimes - ktimeh * 3600 ) / 60;
414  ktimes = ktimes - ktimeh * 3600 - ktimem * 60;
415 
416  // Compose and display completed-scan status message
417  QString astat = tr( "%1 alphas evaluated in" ).arg( nalpha );
418 
419  if ( ktimeh > 0 )
420  astat += tr( " %1 hr., %2 min., %3 sec." )
421  .arg( ktimeh ).arg( ktimem ).arg( ktimes );
422 
423  else
424  astat += tr( " %1 min., %2 sec." ).arg( ktimem ).arg( ktimes );
425 
426  astat += tr( ", in %1 threads" ).arg( nthr );
427  le_stattext->setText( astat );
428  qApp->processEvents();
429 
430  plot_data();
431 }
432 
433 // Plot the data
434 void US_RpScan::plot_data()
435 {
436  data_plot1->detachItems();
437 
438  if ( nalpha > 1 )
439  {
440  int je = nalpha - 1;
441  double ainc = alphas[ 1 ] - alphas[ 0 ];
442  data_plot1->setTitle(
443  tr( "Alpha Scan Points\nAlphas from %1 to %2 by %3" )
444  .arg( alphas[ 0 ] ).arg( alphas[ je ] ).arg( ainc ) );
445  }
446  else
447  {
448  data_plot1->setTitle( tr( "Alpha Scan Points" ) );
449  }
450 
451  data_plot1->setAxisTitle( QwtPlot::xBottom,
452  tr( "Variance (x 1e%1)" ).arg( lgv ) ),
453  data_plot1->setAxisTitle( QwtPlot::yLeft,
454  ( lgx == 0 ) ?
455  tr( "Norm of X (solute concentrations)" ) :
456  tr( "Norm of X (solute concentrations x 1e%1)" ).arg( lgx ) );
457 
458  grid = us_grid( data_plot1 );
459 
460  QwtPlotCurve* curvln;
461  QwtPlotCurve* curvpt;
462 
463  double* xx = varias.data();
464  double* yy = xnorms.data();
465  QPen pen_plot( US_GuiSettings::plotCurve(), 1 );
466 
467  // Draw the variance,xnorm line
468  curvln = us_curve( data_plot1, tr( "Curve V-X" ) );
469  curvln->setPen( pen_plot );
470  curvln->setData( xx, yy, nalpha );
471 
472  // Show the variance,xnorm points
473  curvpt = us_curve( data_plot1, tr( "Alpha Points" ) );
474  QwtSymbol sym;
475  sym.setStyle( QwtSymbol::Ellipse );
476  sym.setPen ( QPen( Qt::blue ) );
477  sym.setBrush( QBrush( Qt::white ) );
478  sym.setSize ( 8 );
479  curvpt->setStyle( QwtPlotCurve::NoCurve );
480  curvpt->setSymbol( sym );
481  curvpt->setData( xx, yy, nalpha );
482 
483  // Compute and show lines that hint at the elbow point of the curve
484  if ( nalpha > 6 )
485  {
486  double xh[ 2 ];
487  double yh[ 2 ];
488  double xa[ 5 ];
489  double ya[ 5 ];
490  double* xe = (double*)xa;
491  double* ye = (double*)ya;
492  QwtPlotCurve* curvh1;
493  QwtPlotCurve* curvh2;
494  QwtPlotCurve* curvh3;
495  double slope; double slop2;
496  double intcp; double intc2;
497  double sigma; double corre;
498  double xl1p1; double xl1p2; double yl1p1; double yl1p2;
499  double xl2p1; double xl2p2; double yl2p1; double yl2p2;
500  double xl3p1; double xl3p2; double yl3p1; double yl3p2;
501  QPen pen_red ( Qt::red, 0, Qt::DashLine );
502  QPen pen_cyan( Qt::cyan, 0, Qt::DashLine );
503 
504  // Compute a line fitted to the first few main curve points
505  int nlp = 5;
506 
507  while ( nlp < nalpha )
508  {
509  double avg = US_Math2::linefit( &xx, &yy, &slope, &intcp,
510  &sigma, &corre, nlp );
511 DbgLv(1) << "TRP:H1: avg" << avg << "nlp" << nlp;
512 DbgLv(1) << "TRP:H1: sl" << slope << "in" << intcp << "sg" << sigma
513  << "co" << corre;
514  if ( slope < 1e99 && slope > (-1e99) )
515  break;
516 
517  nlp += 2;
518  }
519 
520 DbgLv(1) << "TRP:H1: intcp slope" << intcp << slope;
521  yl1p1 = yy[ 0 ];
522  xl1p1 = ( slope == 0.0 ) ? xx[ 0 ] : ( yl1p1 - intcp ) / slope;
523  yl1p2 = yy[ nalpha - 5 ];
524  xl1p2 = ( slope == 0.0 ) ? xl1p1 : ( yl1p2 - intcp ) / slope;
525 DbgLv(1) << "TRP:H1: l1: x1,y1,x2,y2" << xl1p1 << yl1p1 << xl1p2 << yl1p2;
526 
527  // Compute a line fitted to the last few main curve points
528  int je = nalpha - 1;
529  for ( int jj = 0; jj < 5; jj++, je-- )
530  {
531  xe[ jj ] = xx[ je ];
532  ye[ jj ] = yy[ je ];
533 DbgLv(1) << "TRP:H2: jj x y" << jj << xe[jj] << ye[jj] << "je" << je;
534  }
535  US_Math2::linefit( &xe, &ye, &slop2, &intc2, &sigma, &corre, 5 );
536 DbgLv(1) << "TRP:H2: intcp slope" << intc2 << slop2;
537  xl2p1 = xe[ 0 ];
538  yl2p1 = slop2 * xl2p1 + intc2;
539  xl2p2 = xx[ 4 ];
540  yl2p2 = slop2 * xl2p2 + intc2;
541 DbgLv(1) << "TRP:H2: l2: x1,y1,x2,y2" << xl2p1 << yl2p1 << xl2p2 << yl2p2;
542 
543  // Find the intersection point for the 2 fitted lines
544  US_Math2::intersect( slope, intcp, slop2, intc2, &xl3p1, &yl3p1 );
545 DbgLv(1) << "TRP:H3: l3: x1,y1" << xl3p1 << yl3p1;
546 
547  // Find the curve point nearest to the intersection point;
548  // then compute a line from intersection to nearest curve point.
549  US_Math2::nearest_curve_point( xx, yy, nalpha, true, xl3p1, yl3p1,
550  &xl3p2, &yl3p2, alphas.data(), &alpha );
551 
552  // Do a sanity check. If the intersection point is outside the
553  // rectangle that encloses the curve, we likely have an aberrant curve.
554  // So, skip plotting hint lines that are probably just confusing.
555  double xcvp1 = xx[ 0 ];
556  double ycvp1 = yy[ 0 ];
557  double xcvp2 = xx[ je ];
558  double ycvp2 = yy[ je ];
559  bool do_plot = ( xl3p1 >= xcvp1 && yl3p1 <= ycvp1 &&
560  xl3p1 <= xcvp2 && yl3p1 >= ycvp2 );
561 DbgLv(1) << "TRP:T4: cv: x1,y1" << xcvp1 << ycvp1
562  << "x2,y2" << xcvp2 << ycvp2 << " do_plot" << do_plot;
563 
564  if ( do_plot )
565  {
566  // Plot the alpha hint lines
567  curvh1 = us_curve( data_plot1, tr( "Curve Hint 1" ) );
568  curvh1->setPen( pen_red );
569  xh[ 0 ] = xl1p1; // Line fitted to first few points
570  yh[ 0 ] = yl1p1;
571  xh[ 1 ] = xl3p1;
572  yh[ 1 ] = yl3p1;
573  curvh1->setData( xh, yh, 2 );
574 
575  curvh2 = us_curve( data_plot1, tr( "Curve Hint 2" ) );
576  curvh2->setPen( pen_red );
577  xh[ 0 ] = xl2p1; // Line fitted to last few points
578  yh[ 0 ] = yl2p1;
579  xh[ 1 ] = xl3p1;
580  yh[ 1 ] = yl3p1;
581  curvh2->setData( xh, yh, 2 );
582 
583  curvh3 = us_curve( data_plot1, tr( "Curve Hint 3" ) );
584  curvh3->setPen( pen_cyan );
585  xh[ 0 ] = xl3p1; // Line between hint intersection and main curve
586  yh[ 0 ] = yl3p1;
587  xh[ 1 ] = xl3p2;
588  yh[ 1 ] = yl3p2;
589 DbgLv(1) << "TRP:H3: x1,y1,x2,y2" << xl3p1 << yl3p1 << yl3p2 << yl3p2;
590  curvh3->setData( xh, yh, 2 );
591  }
592 
593  else // If no good intersection point, default the half-way alpha
594  alpha = alphas[ nalpha / 2 ];
595 
596  // Use nearest curve point for default alpha
597  le_selalpha->setText( QString().sprintf( "%.3f", alpha ) );
598  }
599 
600  connect( pick, SIGNAL( cMouseUp( const QwtDoublePoint& ) ),
601  this, SLOT ( mouse ( const QwtDoublePoint& ) ) );
602 
603  data_plot1->replot();
604 }
605 
606 // Handle a mouse click near to a curve point location
607 void US_RpScan::mouse( const QwtDoublePoint& p )
608 {
609  double xloc = p.x();
610  double yloc = p.y();
611  double xcurv;
612  double ycurv;
613  double xlabl;
614  double ylabl;
615 
616  US_Math2::nearest_curve_point( varias.data(), xnorms.data(), nalpha, true,
617  xloc, yloc, &xcurv, &ycurv, alphas.data(), &alpha );
618 
619  alpha = (double)qRound( alpha * 1000.0 ) / 1000.0;
620  QString salpha = QString().sprintf( "%.3f", alpha );
621  le_selalpha->setText( salpha );
622 
623  // Mark selected curve point and give it a label
624  xlabl = xcurv + qAbs( varias[ 0 ] - varias[ nalpha - 1 ] ) / 10.0;
625  ylabl = ycurv + qAbs( xnorms[ 0 ] - xnorms[ nalpha - 1 ] ) / 10.0;
626 
627  data_plot1->detachItems( QwtPlotItem::Rtti_PlotMarker );
628  QwtPlotMarker* msymbo = new QwtPlotMarker;
629  QBrush sbrush( Qt::cyan );
630  QPen spen ( sbrush, 2.0 );
631  msymbo->setValue( xcurv, ycurv );
632  msymbo->setSymbol(
633  QwtSymbol( QwtSymbol::Cross, sbrush, spen, QSize( 8, 8 ) ) );
634  msymbo->attach ( data_plot1 );
635  QwtPlotMarker* marker = new QwtPlotMarker;
636  QwtText label;
637  label.setText ( salpha );
638  label.setColor( Qt::cyan );
639  marker->setValue( xlabl, ylabl );
640  marker->setLabel( label );
641  marker->setLabelAlignment( Qt::AlignHCenter | Qt::AlignVCenter );
642  marker->attach ( data_plot1 );
643 
644  data_plot1->replot();
645 }
646 
647 void US_RpScan::process_job( WorkerThreadPc* wthrd )
648 {
649 DbgLv(1) << "SCPJ: IN";
650  WorkPacketPc wresult;
651  nacomp++; // Bump alphas-complete count
652  b_progress->setValue( nacomp );
653  wthrd->get_result( wresult ); // Get results of thread task
654  int thrn = wresult.thrn; // Thread number of task
655  int taskx = wresult.taskx; // Task index of task
656 DbgLv(1) << "SCPJ: thrn" << thrn << "taskx" << taskx;
657  double variance = wresult.sim_vals.variance;
658  double xnormsq = wresult.sim_vals.xnormsq;
659 DbgLv(1) << "SCPJ: a v x" << alphas[taskx] << variance << xnormsq;
660  varias[ taskx ] = variance;
661  xnorms[ taskx ] = xnormsq;
662 // delete wthrd;
663  QString astat =
664  tr( "Of %1 models, %2 done (Alpha %3, thread %4)" )
665  .arg( nalpha ).arg( nacomp ).arg( alphas[ taskx ] ).arg( thrn );
666  le_stattext->setText( astat );
667 DbgLv(1) << "SCPJ: astat" << astat;
668 
669  if ( nacomp == nalpha )
670  { // We are done with all tasks
671  qApp->processEvents();
672 DbgLv(1) << "SCPJ: ALL-DONE";
673  return;
674  }
675 
676  if ( nasubm < nalpha )
677  { // We need to set up to submit another task
678  WorkPacketPc wtask = wresult;
679  wtask.taskx = nasubm;
680  wtask.sim_vals.alpha = alphas[ nasubm ];
681  wtask.dsets = dsets;
682  wtask.depth = 1;
683  wtask.state = 0;
684  wtask.psv_nnls_a = &sv_nnls_a;
685  wtask.psv_nnls_b = &sv_nnls_b;
686 DbgLv(1) << "SCPJ: new subm nnls_a nnls_b" << wtask.psv_nnls_a
687  << wtask.psv_nnls_b << "b size" << sv_nnls_b.size();
688  wtask.csolutes.clear();
689 DbgLv(1) << "SCPJ: new subm taskx alpha" << nasubm << wtask.sim_vals.alpha;
690 
691 // WorkerThreadPc* wthr = new WorkerThreadPc( this );
692  wthrd->disconnect();
693  WorkerThreadPc* wthr = wthrd;
694 
695  wthr->define_work( wtask );
696 DbgLv(1) << "SCPJ: new subm work defined sv_nnls_a" << wtask.psv_nnls_b;
697 // wthr->disconnect();
698  connect( wthr, SIGNAL( work_complete( WorkerThreadPc* ) ),
699  this, SLOT( process_job( WorkerThreadPc* ) ) );
700  wthr->start();
701  nasubm++;
702 DbgLv(1) << "SCPJ: new subm tsk started" << nasubm;
703  }
704  qApp->processEvents();
705 }
706 
707 void US_RpScan::apply_alpha( const double alpha, QVector< double >* psv_nnls_a,
708  QVector< double >* psv_nnls_b, const int nscans, const int npoints,
709  const int nisols, double& variance, double& xnormsq )
710 {
711  int ntotal = nscans * npoints;
712  int narows = ntotal + nisols;
713  int ncsols = 0;
714  variance = 0.0;
715  xnormsq = 0.0;
716  QVector< double > nnls_a = *psv_nnls_a;
717  QVector< double > nnls_b = *psv_nnls_b;
718  QVector< double > nnls_x;
719  QVector< double > simdat;
720  nnls_x.fill( 0.0, nisols );
721  simdat.fill( 0.0, ntotal );
722 qDebug() << "AA: ns np ni na" << nscans << npoints << nisols << narows;
723 
724  // Replace alpha in the diagonal of the lower square of A
725  int dx = ntotal;
726  int dinc = ntotal + nisols + 1;
727 //qDebug() << "AA: alf alfd" << alpha << alphad << "dx dinc" << dx << dinc;
728 qDebug() << "AA: alf" << alpha << "dx dinc" << dx << dinc;
729 
730  for ( int cc = 0; cc < nisols; cc++ )
731  {
732  nnls_a[ dx ] = alpha;
733  dx += dinc;
734  }
735 
736  // Compute the X vector using NNLS
737  US_Math2::nnls( nnls_a.data(), narows, narows, nisols,
738  nnls_b.data(), nnls_x.data() );
739 
740  // Construct the output solutes and the implied simulation and xnorm-sq
741  for ( int cc = 0; cc < nisols; cc++ )
742  {
743  double soluval = nnls_x[ cc ]; // Computed concentration, this solute
744 
745  if ( soluval > 0.0 )
746  {
747  xnormsq += sq( soluval );
748  ncsols++;
749  int aa = cc * narows;
750 
751  for ( int kk = 0; kk < ntotal; kk++ )
752  {
753  simdat[ kk ] += ( soluval * (*psv_nnls_a)[ aa++ ] );
754  }
755  }
756  }
757 
758  // Calculate the sum for the variance computation
759  for ( int kk = 0; kk < ntotal; kk++ )
760  {
761  variance += sq( ( (*psv_nnls_b)[ kk ] - simdat[ kk ] ) );
762  }
763 qDebug() << "AA: ntotal" << ntotal << "varisum" << variance;
764 
765  // Return computed variance and xnorm-sq
766  variance /= (double)ntotal;
767 qDebug() << "AA: alpha" << alpha << "vari xnsq" << variance << xnormsq
768  << "ncsols" << ncsols;
769 int mm = npoints / 2;
770 qDebug() << "AA: mm=" << mm << "a[m] b[m] s[m]" << (*psv_nnls_a)[mm]
771  << (*psv_nnls_b)[mm] << simdat[mm];
772 }
773 
Generated on Sun Jun 21 2015 11:40:43 for UltraScan III by   doxygen 1.8.3.1-20130324