UltraScan III
us_pcsa.cpp
Go to the documentation of this file.
1 
3 #include <QApplication>
4 #include <QtSvg>
5 
6 #include "us_pcsa.h"
7 #include "us_resids_bitmap.h"
8 #include "us_plot_control_pc.h"
10 #include "us_license_t.h"
11 #include "us_license.h"
12 #include "us_settings.h"
13 #include "us_gui_settings.h"
14 #include "us_matrix.h"
15 #include "us_constants.h"
16 #include "us_analyte_gui.h"
17 #include "us_passwd.h"
18 #include "us_db2.h"
19 #include "us_data_loader.h"
20 #include "us_util.h"
21 #include "us_investigator.h"
22 #include "us_noise_loader.h"
23 #include "us_loadable_noise.h"
24 
26 // the class US_pcsa.
27 
28 int main( int argc, char* argv[] )
29 {
30  QApplication application( argc, argv );
31 
32  #include "main1.inc"
33 
34  // License is OK. Start up.
35 
36  US_pcsa w;
37  w.show();
38  return application.exec();
39 }
40 
41 // constructor, based on AnalysisBase
43 {
44  setWindowTitle( tr( "Parametrically Constrained Spectrum Analysis" ) );
45  setObjectName( "US_pcsa" );
47  clean_etc_dir();
48  baserss = 0;
49 
50  // Build local and pcsa-specific GUI elements
51  te_results = NULL;
52 
53  QLabel* lb_analysis = us_banner( tr( "Analysis Controls" ) );
54  QLabel* lb_scan = us_banner( tr( "Scan Control" ) );
55 
56  QLabel* lb_status = us_label ( tr( "Status\nInfo:" ) );
58  QLabel* lb_from = us_label ( tr( "Scan focus from:" ) );
59  QLabel* lb_to = us_label ( tr( "to:" ) );
60 
61  pb_exclude = us_pushbutton( tr( "Exclude Scan Range" ) );
62  pb_exclude->setEnabled( false );
63  connect( pb_exclude, SIGNAL( clicked() ), SLOT( exclude() ) );
64 
65  // Effectively disable boundaries to turn off cyan portion of plot2
66  ct_boundaryPercent->disconnect();
67  ct_boundaryPercent->setRange ( 0.0, 300.0, 1.0 );
68  ct_boundaryPercent->setValue ( 300.0 );
69  ct_boundaryPercent->setEnabled( false );
70  ct_boundaryPos ->disconnect();
71  ct_boundaryPos ->setRange ( -50.0, 300.0, 1.0 );
72  ct_boundaryPos ->setValue ( -50.0 );
73  ct_boundaryPos ->setEnabled( false );
74 
75  connect( ct_from, SIGNAL( valueChanged( double ) ),
76  SLOT ( exclude_from( double ) ) );
77  connect( ct_to, SIGNAL( valueChanged( double ) ),
78  SLOT ( exclude_to ( double ) ) );
79  pb_fitcntl = us_pushbutton( tr( "Fit Control" ) );
80  pb_plt3d = us_pushbutton( tr( "3-D Plot" ) );
81  pb_pltres = us_pushbutton( tr( "Residual Plot" ) );
82 
83  connect( pb_fitcntl, SIGNAL( clicked() ), SLOT( open_fitcntl() ) );
84  connect( pb_plt3d, SIGNAL( clicked() ), SLOT( open_3dplot() ) );
85  connect( pb_pltres, SIGNAL( clicked() ), SLOT( open_resplot() ) );
86 
87  // To modify controls layout, first make Base elements invisible
88 
89  QWidget* widg;
90  for ( int ii = 0; ii < controlsLayout->count(); ii++ )
91  if ( ( widg = controlsLayout->itemAt( ii )->widget() ) != 0 )
92  widg->setVisible( false );
93  ct_from ->setVisible( true );
94  ct_to ->setVisible( true );
95  pb_exclude ->setVisible( true );
96  pb_reset_exclude->setVisible( true );
97 
98  // Add variance and rmsd to parameters layout
99  QLabel* lb_vari = us_label ( tr( "Variance:" ) );
100  le_vari = us_lineedit( "0.00000", 0, true );
101  QLabel* lb_rmsd = us_label ( tr( "RMSD:" ) );
102  le_rmsd = us_lineedit( "0.00000", 0, true );
103  int row = parameterLayout->rowCount();
104  parameterLayout->addWidget( lb_vari, row, 0, 1, 1 );
105  parameterLayout->addWidget( le_vari, row, 1, 1, 1 );
106  parameterLayout->addWidget( lb_rmsd, row, 2, 1, 1 );
107  parameterLayout->addWidget( le_rmsd, row++, 3, 1, 1 );
108 
109  // Reconstruct controls layout with some pcsa-specific elements
110  row = 0;
111  controlsLayout->addWidget( lb_scan, row++, 0, 1, 6 );
112  controlsLayout->addWidget( lb_from, row, 0, 1, 2 );
113  controlsLayout->addWidget( ct_from, row++, 2, 1, 4 );
114  controlsLayout->addWidget( lb_to, row, 0, 1, 2 );
115  controlsLayout->addWidget( ct_to, row++, 2, 1, 4 );
116  controlsLayout->addWidget( pb_exclude, row, 0, 1, 2 );
117  controlsLayout->addWidget( pb_reset_exclude, row++, 2, 1, 4 );
118  controlsLayout->addWidget( lb_analysis, row++, 0, 1, 6 );
119  controlsLayout->addWidget( pb_fitcntl, row, 0, 1, 2 );
120  controlsLayout->addWidget( pb_plt3d, row, 2, 1, 2 );
121  controlsLayout->addWidget( pb_pltres, row++, 4, 1, 2 );
122  controlsLayout->addWidget( lb_status, row, 0, 1, 1 );
123  controlsLayout->addWidget( te_status, row, 1, 1, 5 );
124  row += 3;
125 
126  // Set initial status text
127  te_status->setAlignment( Qt::AlignCenter | Qt::AlignVCenter );
128  te_status->setText( tr(
129  "Solution not initiated...\n"
130  "RMSD: 0.000000,\n"
131  "Variance: 0.000000e-05 .\n"
132  "Iterations: 0" ) );
133  us_setReadOnly( te_status, true );
134 
135  connect( pb_help, SIGNAL( clicked() ), SLOT( help() ) );
136  connect( pb_view, SIGNAL( clicked() ), SLOT( view() ) );
137  connect( pb_save, SIGNAL( clicked() ), SLOT( save() ) );
138  connect( pb_close, SIGNAL( clicked() ), SLOT( close() ) );
139 
140  pb_view ->setEnabled( false );
141  pb_save ->setEnabled( false );
142  pb_fitcntl->setEnabled( false );
143  pb_plt3d ->setEnabled( false );
144  pb_pltres ->setEnabled( false );
145  pb_exclude->setEnabled( false );
146  ct_from ->setEnabled( false );
147  ct_to ->setEnabled( false );
148 
149  edata = 0;
150  resplotd = 0;
151  eplotcd = 0;
152  analcd = 0;
153 
154  rbd_pos = this->pos() + QPoint( 100, 100 );
155  epd_pos = this->pos() + QPoint( 400, 200 );
156  acd_pos = this->pos() + QPoint( 500, 50 );
157 
158  dsets.clear();
159  dsets << &dset;
160 }
161 
162 // slot to handle the completion of a PC spectrum analysis stage
163 void US_pcsa::analysis_done( int updflag )
164 {
165  if ( updflag == (-1) )
166  { // fit has been aborted or reset for new fit
167  pb_view ->setEnabled( false );
168  pb_save ->setEnabled( false );
169  pb_plt3d ->setEnabled( false );
170  pb_pltres ->setEnabled( false );
171 
172  qApp->processEvents();
173  return;
174  }
175 
176  if ( updflag == (-2) )
177  { // update RMSD
178  QString mdesc = model.description;
179  QString avari = mdesc.mid( mdesc.indexOf( "VARI=" ) + 5 );
180  double vari = avari.section( " ", 0, 0 ).toDouble();
181  double rmsd = sqrt( vari );
182  le_vari->setText( QString::number( vari ) );
183  le_rmsd->setText( QString::number( rmsd ) );
184 DbgLv(1) << "Analysis Done VARI" << vari;
185 
186  qApp->processEvents();
187  return;
188  }
189 
190  // if here, an analysis is all done
191 
192  bool plotdata = updflag == 1;
193  bool savedata = updflag == 2;
194 
195 DbgLv(1) << "Analysis Done";
196 DbgLv(1) << " model components size" << model.components.size();
197 DbgLv(1) << " ti_noise size" << ti_noise.values.size() << ti_noise.count;
198 DbgLv(1) << " ri_noise size" << ri_noise.values.size() << ri_noise.count;
199 DbgLv(1) << " edat0 sdat0 rdat0"
200  << edata->value(0,0) << sdata.value(0,0) << rdata.value(0,0);
201 
202  pb_plt3d ->setEnabled( true );
203  pb_pltres->setEnabled( true );
204  pb_view ->setEnabled( true );
205  pb_save ->setEnabled( true );
206 
207  data_plot();
208 
209  if ( plotdata )
210  {
211  open_3dplot();
212  open_resplot();
213  }
214 
215  else if ( savedata )
216  {
217  save();
218  }
219 }
220 
221 // load the experiment data, mostly thru AnalysisBase; then disable view,save
222 void US_pcsa::load( void )
223 {
224  US_AnalysisBase2::load(); // load edited experiment data
225 
226  if ( !dataLoaded ) return;
227 
228  pb_view->setEnabled( false ); // disable view,save buttons for now
229  pb_save->setEnabled( false );
230 
231  def_local = ! disk_controls->db();
232  edata = &dataList[ 0 ]; // point to first loaded data
233  baserss = 0;
234 
235  US_Passwd pw;
236  US_DB2* dbP = disk_controls->db()
237  ? new US_DB2( pw.getPasswd() )
238  : NULL;
239 
240  // Get speed steps from disk or DB experiment
241  if ( dbP != NULL )
242  { // Fetch the speed steps for the experiment from the database
243  QStringList query;
244  QString expID;
245  int idExp = 0;
246  query << "get_experiment_info_by_runID"
247  << runID
248  << QString::number( US_Settings::us_inv_ID() );
249  dbP->query( query );
250 
251  if ( dbP->lastErrno() == US_DB2::OK )
252  {
253  dbP->next();
254  idExp = dbP->value( 1 ).toInt();
256 DbgLv(1) << "SS: ss count" << speed_steps.count() << "idExp" << idExp;
257 if (speed_steps.count()>0 )
258 DbgLv(1) << "SS: ss0 w2tfirst w2tlast timefirst timelast"
259  << speed_steps[0].w2t_first << speed_steps[0].w2t_last
260  << speed_steps[0].time_first << speed_steps[0].time_last;
261  }
262  }
263 
264  else
265  { // Read run experiment file and parse out speed steps
266  QString expfpath = directory + "/" + runID + "."
267  + edata->dataType + ".xml";
268 DbgLv(1) << "LD: expf path" << expfpath;
269  QFile xfi( expfpath )
270  ;
271  if ( xfi.open( QIODevice::ReadOnly ) )
272  { // Read and parse "<speedstep>" lines in the XML
273  QXmlStreamReader xmli( &xfi );
274 
275  while ( ! xmli.atEnd() )
276  {
277  xmli.readNext();
278 
279  if ( xmli.isStartElement() && xmli.name() == "speedstep" )
280  {
281  SP_SPEEDPROFILE sp;
283  speed_steps << sp;
284 DbgLv(1) << "LD: sp: rotspeed" << sp.rotorspeed << "t1" << sp.time_first;
285  }
286  }
287 
288  xfi.close();
289  }
290  }
291 
292  exp_steps = ( speed_steps.count() > 0 ); // Flag any multi-step experiment
293 }
294 
295 // plot the data
296 void US_pcsa::data_plot( void )
297 {
298  // Disable base2 cyan boundary portion
299  ct_boundaryPercent->setValue( 300.0 );
300  ct_boundaryPos ->setValue( -50.0 );
301 
302 DbgLv(1) << "Data Plot by Base";
303  US_AnalysisBase2::data_plot(); // plot experiment data
304 DbgLv(1) << "Data Plot from Base";
305 
306  pb_fitcntl->setEnabled( true );
307  ct_from ->setEnabled( true );
308  ct_to ->setEnabled( true );
309 
310  if ( ! dataLoaded ||
311  sdata.scanData.size() != edata->scanData.size() )
312  return;
313 
314  // set up to plot simulation data and residuals
315  int nscans = edata->scanCount();
316  int npoints = edata->pointCount();
317  int count = ( npoints > nscans ) ? npoints : nscans;
318 
319  QVector< double > rvec( count, 0.0 );
320  QVector< double > vvec( count, 0.0 );
321  double* ra = rvec.data();
322  double* va = vvec.data();
323  QString title;
324  QwtPlotCurve* cc;
325  QPen pen_red( Qt::red );
326  QPen pen_cyan( Qt::cyan );
327  QPen pen_plot( Qt::green );
328 
329  bool have_ri = ri_noise.count > 0;
330  bool have_ti = ti_noise.count > 0;
331  double rl = edata->radius( 0 );
332  double vh = edata->value( 0, 0 ) * 0.05;
333  rl -= 0.05;
334  vh += ( vh - edata->value( 0, 0 ) ) * 0.05;
335 
336  // plot the simulation data in red on top of experiment data
337  for ( int ii = 0; ii < nscans; ii++ )
338  {
339  if ( excludedScans.contains( ii ) ) continue;
340 
341  int kk = 0;
342  double rr = 0.0;
343  double vv = 0.0;
344  double rn = have_ri ? ri_noise.values[ ii ] : 0.0;
345 
346  for ( int jj = 0; jj < npoints; jj++ )
347  {
348  double tn = have_ti ? ti_noise.values[ jj ] : 0.0;
349  rr = sdata.radius( jj );
350  vv = sdata.value( ii, jj++ ) + rn + tn;
351  if ( rr > rl )
352  {
353  ra[ kk ] = rr;
354  va[ kk++ ] = vv;
355  }
356  }
357  title = "SimCurve " + QString::number( ii );
358  cc = us_curve( data_plot2, title );
359  cc->setPen( pen_red );
360  cc->setData( ra, va, kk );
361  }
362 
363  data_plot2->replot(); // replot combined exper,simul data
364 
365  data_plot1->detachItems();
366  data_plot1->clear();
367 
368  us_grid( data_plot1 );
369  data_plot1->setAxisTitle( QwtPlot::xBottom, tr( "Radius (cm)" ) );
370  data_plot1->setAxisTitle( QwtPlot::yLeft, tr( "OD Difference" ) );
371  double vari = 0.0;
372 
373  // build vector of radius values
374  for ( int jj = 0; jj < npoints; jj++ )
375  ra[ jj ] = sdata.radius( jj );
376 
377  // build and plot residual points
378  for ( int ii = 0; ii < nscans; ii++ )
379  {
380  double rinoi = have_ri ? ri_noise.values[ ii ] : 0.0;
381 
382  for ( int jj = 0; jj < npoints; jj++ )
383  {
384  double tinoi = have_ti ? ti_noise.values[ jj ] : 0.0;
385  double evalu = edata->value( ii, jj )
386  - sdata .value( ii, jj ) - tinoi - rinoi;
387  va[ jj ] = evalu;
388  vari += sq( evalu );
389  }
390 
391  // plot dots of residuals at current scan
392  title = "resids " + QString::number( ii );
393  cc = us_curve( data_plot1, title );
394  cc->setPen( pen_plot );
395  cc->setStyle( QwtPlotCurve::Dots );
396  cc->setData( ra, va, npoints );
397  }
398 
399  data_plot1->setAxisAutoScale( QwtPlot::xBottom );
400  data_plot1->setAxisAutoScale( QwtPlot::yLeft );
401  data_plot1->setTitle( tr( "Residuals" ) );
402 
403  // plot zero line through residuals plot
404  double xlo = edata->radius( 0 ) - 0.05;
405  double xhi = edata->radius( npoints - 1 ) + 0.05;
406  ra[ 0 ] = xlo;
407  ra[ 1 ] = xhi;
408  va[ 0 ] = 0.0;
409  va[ 1 ] = 0.0;
410  cc = us_curve( data_plot1, "zero-line" );
411  cc->setPen( QPen( QBrush( Qt::red ), 2 ) );
412  cc->setData( ra, va, 2 );
413 
414  // draw the plot
415  data_plot1->setAxisScale( QwtPlot::xBottom, xlo, xhi );
416  data_plot1->replot();
417 
418  // report on variance and rmsd
419  vari /= (double)( nscans * npoints );
420  rmsd = sqrt( vari );
421  le_vari->setText( QString::number( vari ) );
422  le_rmsd->setText( QString::number( rmsd ) );
423 DbgLv(1) << "Data Plot VARI" << vari;
424 }
425 
426 // view data report
427 void US_pcsa::view( void )
428 {
429  // Create the report text
430  QString rtext;
431  QTextStream ts( &rtext );
432  write_report( ts );
433 
434  // Create US_Editor and display report
435  if ( te_results == NULL )
436  {
437  te_results = new US_Editor( US_Editor::DEFAULT, true, QString(), this );
438  te_results->resize( 820, 700 );
439  QPoint p = g.global_position();
440  te_results->move( p.x() + 30, p.y() + 30 );
441  te_results->e->setFont( QFont( US_GuiSettings::fontFamily(),
443  }
444 
445  te_results->e->setHtml( rtext );
446  te_results->show();
447 }
448 
449 // Save data (model,noise), report, and PNG image files
450 void US_pcsa::save( void )
451 {
452 DbgLv(1) << "SV: IN";
453  QString analysisDate = QDateTime::currentDateTime().toUTC()
454  .toString( "yyMMddhhmm" );
455  QString reqGUID = US_Util::new_guid();
456  QString runID = edata->runID;
457  QString editID = edata->editID.startsWith( "20" ) ?
458  edata->editID.mid( 2 ) :
459  edata->editID;
460  QString dates = "e" + editID + "_a" + analysisDate;
461  int mciters = mrecs_mc.size();
462  QString analysisType = "PCSA";
463 DbgLv(1) << "SV: model_stats size" << model_stats.size();
464  QString curvType = model_stats[ 1 ];
465  if ( curvType.contains( "Straight L" ) )
466  analysisType = "PCSA-SL";
467  else if ( curvType.contains( "Increasing Sig" ) )
468  analysisType = "PCSA-IS";
469  else if ( curvType.contains( "Decreasing Sig" ) )
470  analysisType = "PCSA-DS";
471  else if ( curvType.contains( "Horizontal L" ) )
472  analysisType = "PCSA-HL";
473  else if ( curvType.contains( "Second-Order" ) )
474  analysisType = "PCSA-2O";
475 
476  if ( model.alphaRP != 0.0 )
477  analysisType = analysisType + "-TR";
478 
479  if ( mciters > 1 )
480  analysisType = analysisType + "-MC";
481 DbgLv(1) << "SV: analysisType" << analysisType;
482 
483  QString requestID = "local";
484  QString tripleID = edata->cell + edata->channel + edata->wavelength;
485  QString analysisID = dates + "_" + analysisType + "_" + requestID + "_";
486  QString dext = "." + tripleID;
487  QString dext2 = ".e" + editID + "-" + dext.mid( 1 );
488  QString descbase = runID + "." + tripleID + "." + analysisID;
489 
490  QString reppath = US_Settings::reportDir();
491  QString respath = US_Settings::resultDir();
492  QString tmppath = US_Settings::tmpDir() + "/";
493  QString mdlpath;
494  QString noipath;
495  int knois = min( ti_noise.count, 1 )
496  + min( ri_noise.count, 1 ); // noise files per model
497  double meniscus = edata->meniscus;
498  double dwavelen = edata->wavelength.toDouble();
499  QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );
500 
501  // Test existence or create needed subdirectories
502  if ( ! mkdir( reppath, runID ) )
503  {
504  qDebug() << "*** Unable to create or find the report directory ***";
505  return;
506  }
507 
508  if ( ! US_Model::model_path( mdlpath ) )
509  {
510  qDebug() << "*** Unable to create or find the model directory ***";
511  return;
512  }
513 
514  if ( knois > 0 && ! US_Noise::noise_path( noipath ) )
515  {
516  qDebug() << "*** Unable to create or find the noise directory ***";
517  return;
518  }
519 
520  if ( model.components.size() == 0 )
521  {
522  QMessageBox::critical( this, tr( "Zero-Component Model" ),
523  tr( "*ERROR* The model you are attempting to save\n"
524  "has no components. No data was saved." ) );
525  return;
526  }
527 
528  // Save the model and any noise file(s)
529 
530  US_Passwd pw;
531  US_DB2* dbP = def_local ? NULL : new US_DB2( pw.getPasswd() );
532  QDir dirm( mdlpath );
533  QDir dirn( noipath );
534  mdlpath += "/";
535  noipath += "/";
536  QStringList mfilt( "M*.xml" );
537  QStringList nfilt( "N*.xml" );
538  QStringList mdnams = dirm.entryList( mfilt, QDir::Files, QDir::Name );
539  QStringList ndnams = dirn.entryList( nfilt, QDir::Files, QDir::Name );
540  QStringList mnames;
541  QStringList nnames;
542  QStringList tnames;
543  QString mname = "M0000000.xml";
544  QString tname = mname;
545  QString nname = "N0000000.xml";
546  int indx = 1;
547  int kmodels = 0;
548  int knoises = 0;
549  bool have_ti = ( ti_noise.count > 0 );
550  bool have_ri = ( ri_noise.count > 0 );
551 
552  while( indx > 0 )
553  { // build a list of available model file names
554  mname = "M" + QString().sprintf( "%07i", indx++ ) + ".xml";
555  if ( ! mdnams.contains( mname ) )
556  { // Add to the list of new-name models
557  mnames << mname;
558 DbgLv(1) << "SV: kmodels mciters" << kmodels << mciters << "mname" << mname;
559 
560  if ( ++kmodels >= mciters )
561  break;
562  }
563  }
564 DbgLv(1) << "SV: mnames size" << mnames.size();
565 
566  indx = 1;
567 
568  while( indx > 0 )
569  { // build a list of available noise file names
570  nname = "N" + QString().sprintf( "%07i", indx++ ) + ".xml";
571  if ( ! ndnams.contains( nname ) )
572  { // add to the list of new-name noises
573  nnames << nname;
574  if ( ++knoises >= knois )
575  break;
576  }
577  }
578 double tino = ti_noise.count > 0 ? ti_noise.values[0] : 0.0;
579 double tini = ti_noise_in.count > 0 ? ti_noise_in.values[0] : 0.0;
580 double rino = ri_noise.count > 0 ? ri_noise.values[0] : 0.0;
581 double rini = ri_noise_in.count > 0 ? ri_noise_in.values[0] : 0.0;
582 DbgLv(1) << "SV: Pre-sum tno tni" << tino << tini << "rno rni" << rino << rini;
583 
584  // Output model and noises
585 DbgLv(1) << "SV: mrecs size" << mrecs.size();
586  double variance = mrecs[ 0 ].variance;
587  QString iterID = "i01";
588 DbgLv(1) << "SV: variance" << variance;
589 
590  if ( mciters < 2 )
591  { // Output the single non-MC model
592 DbgLv(1) << "SV: non-MC model ncomp" << model.components.size();
593  model.description = descbase + iterID + ".model";
594  mname = mdlpath + mnames[ 0 ];
597  model.requestGUID = reqGUID;
599  model.variance = variance;
600  model.meniscus = meniscus;
601  model.wavelength = dwavelen;
603 
604  for ( int cc = 0; cc < model.components.size(); cc++ )
605  model.components[ cc ].name = QString().sprintf( "SC%04d", cc + 1 );
606 
607  // Output the model
608  if ( dbP != NULL )
609  model.write( dbP );
610  else
611  model.write( mname );
612 
613 DbgLv(1) << "SV: Post-wrdb tno tni"
614  << ti_noise.values[0] << ti_noise_in.values[0];
615  }
616 
617  else
618  { // Output the individual models of a MonteCarlo set
619 DbgLv(1) << "SV: MC models mciters" << mciters;
620  for ( int jmc = 0; jmc < mciters; jmc++ )
621  {
622  iterID = QString().sprintf( "mc%04d", jmc + 1 );
623  model = mrecs_mc[ jmc ].model;
624  model.description = descbase + iterID + ".model";
625  tname = tmppath + descbase + iterID + ".mdl.tmp";
628  model.requestGUID = reqGUID;
630 
631 DbgLv(1) << "SV: MC models jmc" << jmc << "write" << tname;
632  // Output the model
633  model.write( tname );
634 
635  tnames << tname;
636  }
637 
638  // Now output the composite model
639  tname = US_Model::composite_mc_file( tnames, true );
640  mname = mdlpath + mnames[ 0 ];
641 DbgLv(1) << "SV: MC models write complete mname" << mname;
642  QFile tfile( tname );
643  if ( dbP != NULL )
644  { // DB: load, write to DB, delete
645  model.load( tname );
646  model.description = QString( tname ).section( "/", -1, -1 )
647  .replace( ".model.xml", ".model" )
648  .replace( ".mdl.tmp", ".model" );
649  model.write( dbP );
650  tfile.remove();
651 DbgLv(1) << "SV: MC models write DB from tname" << tname;
652  }
653 
654  else
655  { // Local: rename created file to models directory
656  tfile.rename( mname );
657 DbgLv(1) << "SV: MC models write Local by rename tname" << tname;
658  }
659  }
660 
661  int kk = 0;
662  int err = 0;
663 
664  if ( have_ti )
665  { // output a TI noise
666  ti_noise.description = descbase + iterID + ".ti_noise";
670  nname = noipath + nnames[ kk++ ];
671  int nicount = ti_noise_in.count;
672 DbgLv(1) << "SV: TI nicount" << nicount;
673 
674 DbgLv(1) << "SV: Pre-sum tno tni"
675  << ti_noise.values[0] << ti_noise_in.values[0] << "counts sizes"
676  << ti_noise.count << nicount << ti_noise.values.size()
677  << ti_noise_in.values.size();
678  if ( nicount > 0 ) // Sum in any input noise
679  ti_noise.sum_noise( ti_noise_in, true );
680 DbgLv(1) << "SV: Post-sum tno" << ti_noise.values[0];
681 
682  err = ti_noise.write( nname );
683  if ( err != US_DB2::OK )
684  qDebug() << "*ERROR* writing noise" << nname;
685 
686  if ( dbP != NULL )
687  {
688  err = ti_noise.write( dbP );
689  if ( err != US_DB2::OK )
690  qDebug() << "*ERROR* writing noise to DB" << ti_noise.description;
691  }
692 
693  if ( nicount > 0 ) // Remove input noise in case re-plotted
694  {
695  US_Noise noise_rmv = ti_noise_in;
696 
697  for ( int kk = 0; kk < nicount; kk++ )
698  noise_rmv.values[ kk ] *= -1.0;
699 
700  ti_noise.sum_noise( noise_rmv, true );
701  }
702  }
703 
704  if ( have_ri )
705  { // output an RI noise
706  ri_noise.description = descbase + iterID + ".ri_noise";
710  nname = noipath + nnames[ kk++ ];
711  int nicount = ri_noise_in.count;
712 DbgLv(1) << "SV: RI nicount" << nicount;
713 
714  if ( nicount > 0 ) // Sum in any input noise
715  ri_noise.sum_noise( ri_noise_in, true );
716 DbgLv(1) << "SV: Post-sum rno" << ri_noise.values[0];
717 
718  err = ri_noise.write( nname );
719  if ( err != US_DB2::OK )
720  qDebug() << "*ERROR* writing noise" << nname;
721 
722  if ( dbP != NULL )
723  {
724  err = ri_noise.write( dbP );
725  if ( err != US_DB2::OK )
726  qDebug() << "*ERROR* writing noise to DB" << ri_noise.description;
727  }
728 
729  if ( nicount > 0 ) // Remove input noise in case re-plotted
730  {
731  US_Noise noise_rmv = ri_noise_in;
732 
733  for ( int kk = 0; kk < nicount; kk++ )
734  noise_rmv.values[ kk ] *= -1.0;
735 
736  ri_noise.sum_noise( noise_rmv, true );
737  }
738  }
739 tino = ti_noise.count > 0 ? ti_noise.values[0] : 0.0;
740 rino = ri_noise.count > 0 ? ri_noise.values[0] : 0.0;
741 DbgLv(1) << "SV: Post-write tno rno" << tino << rino;
742 
743  if ( dbP != NULL )
744  {
745  delete dbP;
746  dbP = NULL;
747  }
748 
749  reppath = reppath + "/" + runID + "/";
750  respath = respath + "/" + runID + "/";
751  QString filebase = reppath + analysisType + dext + ".";
752  QString htmlFile = filebase + "report.html";
753  QString plot1File = filebase + "velocity.svgz";
754  QString plot2File = filebase + "residuals.png";
755  QString plot3File = filebase + "rbitmap.png";
756  QString plot4File = filebase + "mlines.png";
757  QString ptmp4File = tmppath + "PCSA" + dext + ".mlines."
758  + QString::number( getpid() ) + ".png";
759 DbgLv(1) << "mlines ptmp4File" << ptmp4File;
760 
761  // Write HTML report file
762  QFile rep_f( htmlFile );
763 
764  if ( ! rep_f.open( QIODevice::WriteOnly | QIODevice::Text ) )
765  return;
766 
767  QTextStream ts( &rep_f );
768  write_report( ts );
769  rep_f.close();
770 
771  // Write plots
772  write_plot( plot1File, data_plot2 );
773  write_plot( plot2File, data_plot1 );
774  write_bmap( plot3File );
775 
776  QFile::remove( plot4File );
777  if ( QFile::copy ( ptmp4File, plot4File ) )
778  QFile::remove( ptmp4File );
779 
780  // use a dialog to tell the user what we've output
781  QString wmsg = tr( "Wrote:\n" );
782  wmsg += mname + "\n"; // List the model file
783 
784  if ( knois > 0 )
785  {
786  nname = noipath + nnames[ 0 ]; // List 1st noise file(s)
787  wmsg = wmsg + nname + "\n";
788 
789  if ( knois > 1 )
790  {
791  nname = noipath + nnames[ 1 ];
792  wmsg = wmsg + nname + "\n";
793  }
794  }
795 
796  // list report and plot files
797  wmsg = wmsg + htmlFile + "\n"
798  + plot1File + "\n"
799  + plot2File + "\n"
800  + plot3File + "\n"
801  + plot4File + "\n";
802  QStringList repfiles;
803  update_filelist( repfiles, htmlFile );
804  update_filelist( repfiles, plot1File );
805  update_filelist( repfiles, plot2File );
806  update_filelist( repfiles, plot3File );
807  update_filelist( repfiles, plot4File );
808 
809  if ( disk_controls->db() )
810  { // Write report files to the database
811  reportFilesToDB( repfiles );
812 
813  wmsg += tr( "\nReport files were also saved to the database." );
814  }
815 
816  QApplication::restoreOverrideCursor();
817  QMessageBox::information( this, tr( "Successfully Written" ), wmsg );
818 }
819 
820 // Return pointer to main window edited data
822 {
823  int drow = lw_triples->currentRow();
824  edata = ( drow >= 0 ) ? &dataList[ drow ] : 0;
825 
826  return edata;
827 }
828 
829 // Return pointers to main window data and GUI elements
830 
833 QList< int >* US_pcsa::mw_excllist() { return &excludedScans; }
837 QPointer< QTextEdit > US_pcsa::mw_status_text() { return te_status; }
838 QStringList* US_pcsa::mw_model_stats() { return &model_stats; }
839 QVector< US_ModelRecord >* US_pcsa::mw_mrecs() { return &mrecs; }
840 QVector< US_ModelRecord >* US_pcsa::mw_mrecs_mc() { return &mrecs_mc; }
841 int* US_pcsa::mw_base_rss() { return &baserss; }
842 
843 // Open residuals plot window
845 {
846  if ( resplotd )
847  {
848  rbd_pos = resplotd->pos();
849  resplotd->close();
850  }
851  else
852  rbd_pos = this->pos() + QPoint( 100, 100 );
853 
854  resplotd = new US_ResidPlotPc( this );
855  resplotd->move( rbd_pos );
856  resplotd->setVisible( true );
857  connect( resplotd, SIGNAL( destroyed ( QObject *) ),
858  this, SLOT( child_closed( QObject* ) ) );
859 }
860 
861 // Open 3-D plot control window
863 {
864  if ( eplotcd )
865  {
866  epd_pos = eplotcd->pos();
867  eplotcd->close();
868  }
869  else
870  epd_pos = this->pos() + QPoint( 400, 200 );
871 
872  eplotcd = new US_PlotControlPc( this, &model );
873  eplotcd->move( epd_pos );
874  eplotcd->show();
875  connect( eplotcd, SIGNAL( destroyed ( QObject *) ),
876  this, SLOT( child_closed( QObject* ) ) );
877 }
878 
879 // Open fit analysis control window
881 {
882  int drow = lw_triples->currentRow();
883  if ( drow < 0 ) return;
884 
885  edata = &dataList[ drow ];
886  double avTemp = edata->average_temperature();
887  double vbar20 = US_Math2::calcCommonVbar( solution_rec, 20.0 );
888  double vbartb = US_Math2::calcCommonVbar( solution_rec, avTemp );
889  double buoy = 1.0 - vbar20 * DENS_20W;
890 
891  if ( buoy <= 0.0 )
892  {
893  QMessageBox::critical( this, tr( "Negative Buoyancy Implied" ),
894  tr( "The current vbar20 value (%1) implies a buoyancy\n"
895  "value (%2) that is non-positive.\n\n"
896  "PCSA cannot proceed with this value. Click on the\n"
897  "<Solution> button and change the vbar20 value.\n"
898  "Note that the Solution may be accepted without being saved.\n"
899  "Include negative values in the sedimentation coefficient\n"
900  "range to represent floating data." ).arg( vbar20 ).arg( buoy ) );
901  return;
902  }
903 
905  sd.density = density;
906  sd.viscosity = viscosity;
907  sd.vbar20 = vbar20;
908  sd.vbar = vbartb;
909  sd.manual = manual;
910  US_Math2::data_correction( avTemp, sd );
911 
912  US_Passwd pw;
913  US_DB2* dbP = disk_controls->db() ? new US_DB2( pw.getPasswd() ) : NULL;
914 
915  dset.simparams.initFromData( dbP, dataList[ drow ], !exp_steps );
916 
917  if ( exp_steps )
918  dset.simparams.speed_step = speed_steps;
919 
920  dset.requestID = disk_controls->db() ? "DB" : "Disk";
921  dset.run_data = dataList[ drow ];
922  int atrx = 0;
923  int atry = 1;
924  int atrz = 3;
925  dset.solute_type = ( atrx << 6 ) | ( atry << 3 ) | atrz;
926  dset.viscosity = viscosity;
927  dset.density = density;
928  dset.manual = manual;
929  dset.temperature = avTemp;
930  dset.vbar20 = vbar20;
931  dset.vbartb = vbartb;
932  dset.s20w_correction = sd.s20w_correction;
933  dset.D20w_correction = sd.D20w_correction;
934 DbgLv(1) << "Bottom" << dset.simparams.bottom << "rotorcoeffs"
935  << dset.simparams.rotorcoeffs[0] << dset.simparams.rotorcoeffs[1];
936 
937  if ( dbP != NULL )
938  {
939  delete dbP;
940  dbP = NULL;
941  }
942 
943  if ( analcd != 0 )
944  {
945  acd_pos = analcd->pos();
946  analcd->close();
947  }
948  else
949  acd_pos = this->pos() + QPoint( 500, 50 );
950 
951  analcd = new US_AnalysisControlPc( dsets, this );
952  analcd->move( acd_pos );
953  analcd->show();
954  connect( analcd, SIGNAL( destroyed ( QObject *) ),
955  this, SLOT( child_closed( QObject* ) ) );
956  qApp->processEvents();
957 }
958 
959 // Distribution information HTML string
961 {
962  int ncomp = model.components.size();
963 DbgLv(1) << "distrinfo: ncomp" << ncomp;
964  QString maDesc = model.description;
965  QString runID = edata->runID;
966 
967  if ( ncomp == 0 )
968  return "";
969 
970  if ( maDesc.startsWith( runID ) )
971  { // Saved model: get analysis description from model description
972  maDesc = maDesc.section( ".", -2, -2 ).section( "_", 1, -1 );
973  }
974  else
975  { // No saved model yet: compose analysis description
976  maDesc = "a" + QDateTime::currentDateTime().toUTC()
977  .toString( "yyMMddhhmm" ) + "_PCSA_local_01";
978  }
979 
980  QString mstr = "\n" + indent( 4 )
981  + tr( "<h3>Data Analysis Settings:</h3>\n" )
982  + indent( 4 ) + "<table>\n";
983 
984  mstr += table_row( tr( "Model Analysis:" ),
985  maDesc );
986  mstr += table_row( tr( "Number of Components:" ),
987  QString::number( ncomp ) );
988  mstr += table_row( tr( "Residual RMS Deviation:" ),
989  QString::number( rmsd ) );
990 
991  double sum_mw = 0.0;
992  double sum_s = 0.0;
993  double sum_D = 0.0;
994  double sum_c = 0.0;
995  double mink = 1e+99;
996  double maxk = -1e+99;
997  double minv = 1e+99;
998  double maxv = -1e+99;
999 
1000  for ( int ii = 0; ii < ncomp; ii++ )
1001  {
1002  double conc = model.components[ ii ].signal_concentration;
1003  sum_c += conc;
1004  sum_mw += model.components[ ii ].mw * conc;
1005  sum_s += model.components[ ii ].s * conc;
1006  sum_D += model.components[ ii ].D * conc;
1007  mink = qMin( mink, model.components[ ii ].f_f0 );
1008  maxk = qMax( maxk, model.components[ ii ].f_f0 );
1009  minv = qMin( minv, model.components[ ii ].vbar20 );
1010  maxv = qMax( maxv, model.components[ ii ].vbar20 );
1011  }
1012 
1013  bool cnstvb = ( ( maxk - mink ) / qAbs( maxk )
1014  > ( maxv - minv ) / qAbs( maxv ) );
1015 
1016  mstr += table_row( tr( "Weight Average s20,W:" ),
1017  QString().sprintf( "%6.4e", ( sum_s / sum_c ) ) );
1018  mstr += table_row( tr( "Weight Average D20,W:" ),
1019  QString().sprintf( "%6.4e", ( sum_D / sum_c ) ) );
1020  mstr += table_row( tr( "W.A. Molecular Weight:" ),
1021  QString().sprintf( "%6.4e", ( sum_mw / sum_c ) ) );
1022  mstr += table_row( tr( "Total Concentration:" ),
1023  QString().sprintf( "%6.4e", sum_c ) );
1024 
1025  if ( cnstvb )
1026  mstr += table_row( tr( "Constant Vbar at 20" ) + DEGC + ":",
1027  QString().number( maxv ) );
1028  else
1029  mstr += table_row( tr( "Constant f/f0:" ),
1030  QString().number( maxk ) );
1031 
1032  mstr += indent( 4 ) + "</table>\n\n";
1033 
1034  mstr += indent( 4 ) + tr( "<h3>Distribution Information:</h3>\n" )
1035  + indent( 4 ) + "<table>\n";
1036 
1037  if ( cnstvb )
1038  mstr += table_row( tr( "Molecular Wt." ), tr( "S Apparent" ),
1039  tr( "S 20,W" ), tr( "D Apparent" ),
1040  tr( "D 20,W" ), tr( "f / f0" ),
1041  tr( "Concentration" ) );
1042  else
1043  mstr += table_row( tr( "Molecular Wt." ), tr( "S Apparent" ),
1044  tr( "S 20,W" ), tr( "D Apparent" ),
1045  tr( "D 20,W" ), tr( "Vbar20" ),
1046  tr( "Concentration" ) );
1047 
1048  int drow = lw_triples->currentRow();
1049  edata = &dataList[ drow ];
1050  double avTemp = edata->average_temperature();
1051  double vbar20 = US_Math2::calcCommonVbar( solution_rec, 20.0 );
1052  double vbartb = US_Math2::calcCommonVbar( solution_rec, avTemp );
1054  sd.density = density;
1055  sd.viscosity = viscosity;
1056  sd.manual = manual;
1057  sd.vbar20 = vbar20;
1058  sd.vbar = vbartb;
1059  US_Math2::data_correction( avTemp, sd );
1060 
1061  for ( int ii = 0; ii < ncomp; ii++ )
1062  {
1063  double conc = model.components[ ii ].signal_concentration;
1064  double perc = 100.0 * conc / sum_c;
1065  double s_ap = model.components[ ii ].s;
1066  double D_ap = model.components[ ii ].D;
1067  double f_f0 = model.components[ ii ].f_f0;
1068 
1069  if ( !cnstvb )
1070  {
1071  vbar20 = model.components[ ii ].vbar20;
1072  f_f0 = vbar20;
1073  sd.vbar20 = vbar20;
1074  sd.vbar = US_Math2::adjust_vbar( vbar20, avTemp );
1075  US_Math2::data_correction( avTemp, sd );
1076  }
1077 
1078  s_ap /= sd.s20w_correction;
1079  D_ap /= sd.D20w_correction;
1080 
1081  mstr += table_row(
1082  QString().sprintf( "%10.4e", model.components[ ii ].mw ),
1083  QString().sprintf( "%10.4e", s_ap ),
1084  QString().sprintf( "%10.4e", model.components[ ii ].s ),
1085  QString().sprintf( "%10.4e", D_ap ),
1086  QString().sprintf( "%10.4e", model.components[ ii ].D ),
1087  QString().sprintf( "%10.4e", f_f0 ),
1088  QString().sprintf( "%10.4e (%5.2f %%)", conc, perc ) );
1089  }
1090 
1091  mstr += indent( 4 ) + "</table>\n";
1092 
1093  return mstr;
1094 }
1095 
1096 // Model statistics HTML string
1098 {
1099 DbgLv(1) << "ModStats0" << model_stats[ 0 ];
1100 DbgLv(1) << "ModStats1" << model_stats[ 1 ];
1101  QString mstr = "\n" + indent( 4 )
1102  + tr( "<h3>Model Statistics:</h3>\n" )
1103  + indent( 4 ) + "<table>\n";
1104 
1105  for ( int ii = 0; ii < model_stats.size(); ii += 2 )
1106  {
1107  mstr += table_row( model_stats[ ii ], model_stats[ ii + 1 ] );
1108  }
1109 
1110  mstr += indent( 4 ) + "</table>\n";
1111 
1112  return mstr;
1113 }
1114 
1115 // Write HTML report file
1116 void US_pcsa::write_report( QTextStream& ts )
1117 {
1118  QString curvtype = model_stats[ 1 ];
1119  QString hdr = tr( "Parametrically Constrained Spectrum Analysis" )
1120  + "<br/>( " + curvtype + " )";
1121 
1122  if ( model.alphaRP > 0.0 )
1123  hdr += "<br/>- Tihhonov Regularization";
1124 
1125  ts << html_header( QString( "US_pcsa" ), hdr, edata );
1126  ts << model_statistics();
1127  ts << distrib_info();
1128  ts << indent( 2 ) + "</body>\n</html>\n";
1129 }
1130 
1131 // Write resids bitmap plot file
1132 void US_pcsa::write_bmap( const QString plotFile )
1133 {
1134  // Generate the residuals array
1135  bool have_ri = ri_noise.count > 0;
1136  bool have_ti = ti_noise.count > 0;
1137  int nscans = edata->scanCount();
1138  int npoints = edata->pointCount();
1139  QVector< double > resscn( npoints );
1140  QVector< QVector< double > > resids( nscans );
1141 
1142  for ( int ii = 0; ii < nscans; ii++ )
1143  {
1144  double rnoi = have_ri ? ri_noise.values[ ii ] : 0.0;
1145 
1146  for ( int jj = 0; jj < npoints; jj++ )
1147  {
1148  double tnoi = have_ti ? ti_noise.values[ jj ] : 0.0;
1149  resscn[ jj ] = edata->value( ii, jj )
1150  - sdata.value( ii, jj ) - rnoi - tnoi;
1151  }
1152 
1153  resids[ ii ] = resscn;
1154  }
1155 
1156  // Generate the residuals bitmap plot
1157  US_ResidsBitmap* resbmap = new US_ResidsBitmap( resids );
1158  QPixmap pixmap = QPixmap::grabWidget( resbmap, 0, 0,
1159  resbmap->width(), resbmap->height() );
1160 
1161  // Save the pixmap to the specified file
1162  if ( ! pixmap.save( plotFile ) )
1163  qDebug() << "*ERROR* Unable to write file" << plotFile;
1164 
1165  resbmap->close();
1166  resbmap = 0;
1167 }
1168 
1169 // New triple selected
1170 void US_pcsa::new_triple( int index )
1171 {
1172  edata = &dataList[ index ];
1173 
1174  // Restore pure data type string (other values added in pcsa processing)
1175  edata->dataType = edata->dataType.section( " ", 0, 0 );
1176 
1177  sdata.scanData.clear(); // Clear simulation and upper plot
1178  data_plot1->detachItems();
1179  data_plot1->clear();
1180 
1181  // Temporarily restore loaded noise vectors from triples vectors
1182  ti_noise = tinoises[ index ];
1183  ri_noise = rinoises[ index ];
1184 
1185  US_AnalysisBase2::new_triple( index ); // New triple as in any analysis
1186 
1187  // Move any loaded noise vectors to the "in" versions
1194  tinoises[ index ] = ti_noise_in;
1195  rinoises[ index ] = ri_noise_in;
1196  ti_noise.values.clear();
1197  ri_noise.values.clear();
1198  ti_noise.count = 0;
1199  ri_noise.count = 0;
1200 DbgLv(1) << "NTr: ti noise in n0,type,count,size" << ti_noise_in.values[0]
1202 DbgLv(1) << "NTr: ti noise in minr,maxr"
1204 DbgLv(1) << "NTr: ri noise in n0,type,count,size" << ri_noise_in.values[0]
1206 }
1207 
1208 // Remove any temporary plot file and close all opened windows
1209 void US_pcsa::close( void )
1210 {
1211  QString tripleID = ( edata != 0 )
1213  : "";
1214  QString ptmp4File = US_Settings::tmpDir() + "/PCSA." + tripleID
1215  + ".mlines." + QString::number( getpid() ) + ".png";
1216 
1217  QFile tfile( ptmp4File );
1218  if ( tfile.exists() )
1219  {
1220  tfile.remove();
1221 DbgLv(1) << "pcsa: removed: " << ptmp4File;
1222  }
1223 DbgLv(1) << "pcsa: close d's res epl ana" << resplotd << eplotcd << analcd;
1224 
1225  if ( resplotd != 0 )
1226  resplotd->close();
1227  if ( eplotcd != 0 )
1228  eplotcd->close();
1229  if ( analcd != 0 )
1230  analcd->close();
1231 }
1232 
1233 // Private slot to mark a child widgets as closed, if it has been destroyed
1234 void US_pcsa::child_closed( QObject* o )
1235 {
1236  QString oname = o->objectName();
1237 DbgLv(1) << "pcsa:CC: d's res epl ana" << resplotd << eplotcd << analcd;
1238 
1239  if ( oname.contains( "AnalysisControl" ) )
1240  analcd = 0;
1241  else if ( oname.contains( "ResidPlot" ) )
1242  resplotd = 0;
1243  else if ( oname.contains( "PlotControl" ) )
1244  eplotcd = 0;
1245 DbgLv(1) << "pcsa:CC: return res epl ana" << resplotd << eplotcd << analcd
1246  << "oname" << oname;
1247 }
1248