UltraScan III
pcsa_master.cpp
Go to the documentation of this file.
1 #include "us_mpi_analysis.h"
2 #include "us_math2.h"
3 #include "us_util.h"
4 #include "us_tar.h"
5 #include "us_astfem_rsa.h"
6 #include "us_simparms.h"
7 #include "us_constants.h"
8 
9 #define DbgTime() qDebug() << "TM:" << (startTime.msecsTo(QDateTime::currentDateTime())/1000.0)
10 void US_MPI_Analysis::pcsa_master( void )
11 {
12  current_dataset = 0;
13  datasets_to_process = 1; // Process one dataset at a time for now
14 
15 DbgLv(1) << "pcsa_mast: IN";
16  init_pcsa_solutes();
17 DbgLv(1) << "pcsa_mast: init complete";
18  fill_pcsa_queue();
19 DbgLv(1) << "pcsa_mast: fill_pcsa_queue complete";
20 
21  work_rss.resize( gcores_count );
22 
23  alpha = 0.0;
24  mc_iterations = 0;
25  max_iterations = parameters[ "gfit_iterations" ].toInt();
26  int kcurve = 0;
27 
28  while ( true )
29  {
30  int worker;
31  meniscus_value = data_sets[ current_dataset ]->run_data.meniscus;
32 //if ( max_depth > 1 )
33 // DbgLv(1) << " master loop-TOP: jq-empty?" << job_queue.isEmpty() << " areReady?" << worker_status.contains(READY)
34 // << " areWorking?" << worker_status.contains(WORKING);
35 
36  // Give the jobs to the workers
37  while ( ! job_queue.isEmpty() && worker_status.contains( READY ) )
38  {
39  worker = ready_worker();
40 
41  Sa_Job job = job_queue.takeFirst();
42  job.mpi_job.depth = kcurve++;
43 DbgLv(1) << "pcsa_mast: submit_pcsa kc" << kcurve;
44 
45  submit_pcsa( job, worker );
46  }
47 
48  // All done with the pass if no jobs are ready or running
49  if ( job_queue.isEmpty() && ! worker_status.contains( WORKING ) )
50  {
51  kcurve = 0;
52  qSort( mrecs );
53  US_DataIO::EditedData* edata = &data_sets[ current_dataset ]->run_data;
54  QString tripleID = edata->cell + edata->channel + edata->wavelength;
55 
56  simulation_values.variance = mrecs[ 0 ].variance;
57  simulation_values.zsolutes = mrecs[ 0 ].csolutes;
58  simulation_values.ti_noise = mrecs[ 0 ].ti_noise;
59  simulation_values.ri_noise = mrecs[ 0 ].ri_noise;
60 //*DEBUG*
61 wksim_vals = simulation_values;
62 wksim_vals.zsolutes = mrecs[ 0 ].isolutes;
63 calc_residuals( current_dataset, 1, wksim_vals );
64 DbgLv(1) << "final-mr0: variance1 variance2" << simulation_values.variance << wksim_vals.variance;
65 DbgLv(1) << "final-mr0: csolsize1 csolsize2" << simulation_values.zsolutes.size() << wksim_vals.zsolutes.size();
66 DbgLv(1) << "final-mr0: worker_status" << worker_status;
67 int ki=worker_status.count(INIT);
68 int kw=worker_status.count(WORKING);
69 int kr=worker_status.count(READY);
70 int km=mrecs.size();
71 int kg=0;
72 for(int jm=0; jm<km; jm++ ) if(mrecs[jm].rmsd<1.0) kg++;
73 DbgLv(1) << "final-mr0: ki kw kr" << ki << kw << kr << "km kg" << km << kg;
74 //*DEBUG*
75 
76 
77  QString progress =
78  "Iteration: " + QString::number( iterations );
79 
80  if ( datasets_to_process > 1 )
81  progress += "; Datasets: "
82  + QString::number( datasets_to_process );
83  else
84  progress += "; Dataset: "
85  + QString::number( current_dataset + 1 )
86  + " (" + tripleID + ") of "
87  + QString::number( count_datasets );
88 
89  if ( mc_iterations > 1 )
90  progress += "; MonteCarlo: "
91  + QString::number( mc_iteration + 1 );
92 
93  else
94  progress += "; RMSD: "
95  + QString::number( mrecs[ 0 ].rmsd );
96 
97  send_udp( progress );
98 
99  // Iterative refinement
100  if ( max_iterations > 1 )
101  {
102  if ( data_sets.size() > 1 && iterations == 1 )
103  {
104  if ( datasets_to_process == 1 )
105  {
106  qDebug() << " == Grid-Fit Iterations for Dataset"
107  << current_dataset + 1 << "==";
108  }
109  else
110  {
111  qDebug() << " == Grid-Fit Iterations for Datasets 1 to"
112  << datasets_to_process << "==";
113  }
114  }
115  qDebug() << "Iteration:" << iterations
116  << " Variance:" << mrecs[ 0 ].variance
117  << "RMSD:" << mrecs[ 0 ].rmsd;
118 
119  iterate_pcsa();
120  }
121 
122  if ( ! job_queue.isEmpty() ) continue;
123 
124  iterations = 1;
125  max_iterations = 1;
126 DbgLv(1) << " master loop-BOT: dssize" << data_sets.size() << "ds_to_p"
127  << datasets_to_process << "curr_ds" << current_dataset;
128 US_DataIO::EditedData* edat=&data_sets[current_dataset]->run_data;
129 int ks=edat->scanCount() - 10;
130 int kp=edat->pointCount() - 10;
131 int ss=ks/2;
132 int rr=kp/2;
133 DbgLv(1) << " master loop-BOT: ds" << current_dataset+1 << "data l m h"
134  << edat->value(10,10) << edat->value(ss,rr) << edat->value(ks,kp);
135  // Clean up mrecs of any empty-calculated-solutes records
136  clean_mrecs( mrecs );
137 
138  // Manage multiple data sets in global fit
139  if ( is_global_fit && datasets_to_process == 1 )
140  {
141  global_fit();
142  }
143 DbgLv(1) << " master loop-BOT: GF job_queue empty" << job_queue.isEmpty();
144 
145  if ( ! job_queue.isEmpty() ) continue;
146 
147  if ( is_global_fit )
148  write_global();
149 
150  else
151  write_output();
152 
153  if ( ! job_queue.isEmpty() ) continue;
154 
155  // Save information from best model
156  pcsa_best_model();
157 
158  // Tikhonov Regularization
159  tikreg_pcsa();
160 
161  // Monte Carlo
162  montecarlo_pcsa();
163 
164  if ( is_composite_job )
165  { // Composite job: update outputs in TAR and bump dataset count
166  QString tripleID = QString( data_sets[ current_dataset ]->model
167  .description ).section( ".", -3, -3 );
168  current_dataset++;
169 
170  if ( simulation_values.noisflag == 0 )
171  {
172  DbgLv(0) << my_rank << ": Dataset" << current_dataset
173  << "(" << tripleID << ")"
174  << " : model was output.";
175  }
176  else
177  {
178  DbgLv(0) << my_rank << ": Dataset" << current_dataset
179  << "(" << tripleID << ")"
180  << " : model/noise(s) were output.";
181  }
182 
183  if ( current_dataset < count_datasets )
184  {
185  iterations = 1;
186  mc_iteration = 0;
187  alpha = 0.0;
188  mc_iterations = 0;
189  max_iterations = parameters[ "gfit_iterations" ].toInt();
190  kcurve = 0;
191 
192  fill_pcsa_queue();
193 
194  for ( int ii = 1; ii <= my_workers; ii++ )
195  worker_status[ ii ] = READY;
196 
197  continue;
198  }
199  }
200 
201  shutdown_all(); // All done
202  break; // Break out of main loop.
203  }
204 
205  // Wait for worker to send a message
206  MPI_Status status;
207  int sizes[ 4 ];
208 
209  MPI_Recv( sizes,
210  4,
211  MPI_INT,
212  MPI_ANY_SOURCE,
213  MPI_ANY_TAG,
214  my_communicator,
215  &status);
216 
217 //if ( max_depth > 0 )
218 // DbgLv(1) << " master loop-BOTTOM: status TAG" << status.MPI_TAG << MPI_Job::READY << MPI_Job::RESULTS
219 // << " source" << status.MPI_SOURCE;
220  switch( status.MPI_TAG )
221  {
222  case MPI_Job::READY: // Ready for work
223  worker = status.MPI_SOURCE;
224 DbgLv(1) << " master loop-BOTTOM: READY worker" << worker;
225  worker_status[ worker ] = READY;
226  break;
227 
228  case MPI_Job::RESULTS: // Return solute data
229  worker = status.MPI_SOURCE;
230 DbgLv(1) << " master loop-BOTTOM: RESULTS worker" << worker;
231  process_pcsa_results( worker, sizes );
232  work_rss[ worker ] = sizes[ 3 ];
233  break;
234 
235  default: // Should never happen
236  QString msg = "Master PCSA: Received invalid status " +
237  QString::number( status.MPI_TAG );
238  abort( msg );
239  break;
240  }
241 
242  max_rss();
243  }
244 }
245 
246 // Generate the initial set of solutes
247 void US_MPI_Analysis::init_pcsa_solutes( void )
248 {
249 DbgLv(1) << " init_pcsa_sols: IN";
250  calculated_zsolutes.clear();
251  orig_zsolutes .clear();
252  simulation_values.noisflag = parameters[ "tinoise_option" ].toInt() > 0 ?
253  1 : 0;
254  simulation_values.noisflag += parameters[ "rinoise_option" ].toInt() > 0 ?
255  2 : 0;
256  simulation_values.dbg_level = dbg_level;
257  simulation_values.dbg_timing = dbg_timing;
258 DbgLv(0) << "DEBUG_LEVEL" << simulation_values.dbg_level;
259 
260  double x_min = parameters[ "x_min" ].toDouble();
261  double x_max = parameters[ "x_max" ].toDouble();
262  double y_min = parameters[ "y_min" ].toDouble();
263  double y_max = parameters[ "y_max" ].toDouble();
264  double zval = parameters[ "z_value" ].toDouble();
265  QString s_ctyp = parameters[ "curve_type" ];
266  int ctype = US_ModelRecord::ctype_flag( s_ctyp );
267  QString s_styp = parameters[ "solute_type" ];
268  int stype = US_ModelRecord::stype_flag( s_styp );
269  int nkpts = parameters[ "vars_count" ].toInt();
270  int nkpto = sq( nkpts );
271  int nlpts = parameters[ "curves_points" ].toInt();
272 DbgLv(1) << " init_pcsa_sols: nkpts" << nkpts << "nkpto" << nkpto;
273  double vbar20 = data_sets[ current_dataset ]->vbar20;
274  zval = ( zval == 0.0 ) ? vbar20 : zval;
275 DbgLv(1) << " init_pcsa_sols: currds" << current_dataset << "vbar" << vbar20;
276  double *parlims = (double*)pararry;
277  parlims[ 0 ] = -1.0;
278  parlims[ 4 ] = stype;
279  parlims[ 5 ] = zval;
280  mrecs.clear();
281 DbgLv(1) << " init_pcsa_sols: s_ctyp" << s_ctyp << "ctype" << ctype;
282 
283  if ( ctype == CTYPE_SL )
284  {
285  mrecs.reserve( nkpto );
286 DbgLv(1) << " init_pcsa_sols: call compute_slines nkpts" << nkpts;
287  US_ModelRecord::compute_slines( x_min, x_max, y_min, y_max,
288  nkpts, nlpts, parlims, mrecs );
289 DbgLv(1) << " init_pcsa_sols: compute_slines: mrecs sz" << mrecs.size();
290  }
291 
292  else if ( ctype == CTYPE_IS || ctype == CTYPE_DS )
293  {
294  mrecs.reserve( nkpto );
295 DbgLv(1) << " init_pcsa_sols: call compute_sigmoids nkpts" << nkpts;
296  US_ModelRecord::compute_sigmoids( ctype, x_min, x_max, y_min, y_max,
297  nkpts, nlpts, parlims, mrecs );
298 DbgLv(1) << " init_pcsa_sols: compute_sigmoids: mrecs sz" << mrecs.size();
299  }
300 
301  else if ( ctype == CTYPE_HL )
302  {
303  nkpto = nkpts;
304  mrecs.reserve( nkpto );
305  US_ModelRecord::compute_hlines( x_min, x_max, y_min, y_max,
306  nkpts, nlpts, parlims, mrecs );
307  }
308 
309  else if ( ctype == CTYPE_ALL )
310  {
311  nkpto *= 3;
312  int ctype1 = CTYPE_IS;
313  int ctype2 = CTYPE_DS;
314  double *parlims1 = parlims + 12;
315  double *parlims2 = parlims + 24;
316  parlims1[ 0 ] = parlims[ 0 ];
317  parlims1[ 4 ] = parlims[ 4 ];
318  parlims1[ 5 ] = parlims[ 5 ];
319  parlims2[ 0 ] = parlims[ 0 ];
320  parlims2[ 4 ] = parlims[ 4 ];
321  parlims2[ 5 ] = parlims[ 5 ];
322  mrecs.reserve( nkpto );
323 
324  US_ModelRecord::compute_slines( x_min, x_max, y_min, y_max,
325  nkpts, nlpts, parlims, mrecs );
326  US_ModelRecord::compute_sigmoids( ctype1, x_min, x_max, y_min, y_max,
327  nkpts, nlpts, parlims1, mrecs );
328  US_ModelRecord::compute_sigmoids( ctype2, x_min, x_max, y_min, y_max,
329  nkpts, nlpts, parlims2, mrecs );
330  }
331 
332  else if ( ctype == CTYPE_2O )
333  {
334  nkpto *= nkpts;
335 DbgLv(1) << " init_pcsa_sols: nkpts" << nkpts << "nkpto" << nkpto;
336  mrecs.reserve( nkpto );
337  US_ModelRecord::compute_2ndorder( x_min, x_max, y_min, y_max,
338  nkpts, nlpts, parlims, mrecs );
339  }
340 
341  for ( int ii = 0; ii < mrecs.size(); ii++ )
342  {
343  orig_zsolutes << mrecs[ ii ].isolutes;
344  }
345 }
346 
347 // Submit a PCSA job
348 void US_MPI_Analysis::submit_pcsa( Sa_Job& job, int worker )
349 {
350  job.mpi_job.command = MPI_Job::PROCESS;
351  job.mpi_job.length = job.zsolutes.size();
352  job.mpi_job.meniscus_value = 0.0;
353  job.mpi_job.solution = mc_iteration;
354  job.mpi_job.dataset_offset = current_dataset;
355  job.mpi_job.dataset_count = datasets_to_process;
356 int dd=job.mpi_job.depth;
357 if (dd==0) { DbgLv(1) << "Mast: submit: worker" << worker << " sols"
358  << job.mpi_job.length << "mciter cds" << mc_iteration << current_dataset << " depth" << dd; }
359 else { DbgLv(1) << "Mast: submit: worker" << worker << " sols"
360  << job.mpi_job.length << "mciter cds" << mc_iteration << current_dataset << " depth" << dd; }
361 DbgLv(1) << "Mast: submit: len sol offs cnt"
362  << job.mpi_job.length
363  << job.mpi_job.solution
364  << job.mpi_job.dataset_offset
365  << job.mpi_job.dataset_count;
366 
367  // Tell worker that solutes are coming
368 DbgLv(1) << " master loop-SEND PROCESS worker" << worker
369  << "isolsiz" << job.mpi_job.length << "depth" << job.mpi_job.depth;
370  MPI_Send( &job.mpi_job,
371  sizeof( MPI_Job ),
372  MPI_BYTE,
373  worker, // Send to system that needs work
374  MPI_Job::MASTER,
375  my_communicator );
376 DbgLv(1) << "Mast: submit: send #1";
377 
378  // Send solutes
379  MPI_Send( job.zsolutes.data(),
380  job.mpi_job.length * zsolut_doubles,
381  MPI_DOUBLE, // Pass solute vector as hw independent values
382  worker, // to worker
383  MPI_Job::MASTER,
384  my_communicator );
385 DbgLv(1) << "Mast: submit: send #2";
386 
387  worker_depth [ worker ] = job.mpi_job.depth;
388  worker_status[ worker ] = WORKING;
389 }
390 
391 // Process the results from a just-completed worker task
392 void US_MPI_Analysis::process_pcsa_results( const int worker, const int* sizes )
393 {
394  simulation_values.zsolutes.resize( sizes[ 0 ] );
395  simulation_values.variances.resize( datasets_to_process );
396  simulation_values.ti_noise.resize( sizes[ 1 ] );
397  simulation_values.ri_noise.resize( sizes[ 2 ] );
398 
399  MPI_Status status;
400 
401  // Get all simulation_values
402  MPI_Recv( simulation_values.zsolutes.data(),
403  sizes[ 0 ] * zsolut_doubles,
404  MPI_DOUBLE,
405  worker,
406  MPI_Job::TAG0,
407  my_communicator,
408  &status );
409 
410  MPI_Recv( &simulation_values.variance,
411  1,
412  MPI_DOUBLE,
413  worker,
414  MPI_Job::TAG0,
415  my_communicator,
416  &status );
417 
418  MPI_Recv( simulation_values.variances.data(),
419  datasets_to_process,
420  MPI_DOUBLE,
421  worker,
422  MPI_Job::TAG0,
423  my_communicator,
424  &status );
425 
426  MPI_Recv( simulation_values.ti_noise.data(),
427  sizes[ 1 ],
428  MPI_DOUBLE,
429  worker,
430  MPI_Job::TAG0,
431  my_communicator,
432  &status );
433 
434  MPI_Recv( simulation_values.ri_noise.data(),
435  sizes[ 2 ],
436  MPI_DOUBLE,
437  worker,
438  MPI_Job::TAG0,
439  my_communicator,
440  &status );
441 
442 DbgLv(1) << "Mast: process_results: worker" << worker
443  << " solsize" << sizes[0];
444  Result result;
445  result.depth = worker_depth[ worker ];
446  result.worker = worker;
447  result.zsolutes = simulation_values.zsolutes;
448 
449  // Process the result solutes
450  process_pcsa_solutes( result );
451 }
452 
453 // Process the calculated solute vector from a job result
454 void US_MPI_Analysis::process_pcsa_solutes( Result& result )
455 {
456  int jcurve = result.depth;
457  jcurve = qMax( 0, qMin( ( mrecs.size() - 1 ), jcurve ) );
458 DbgLv(1) << "Mast: process_solutes: worker" << result.worker
459  << "jcurve" << jcurve << " solsize" << result.solutes.size();
460 
461  // Update the model record entry with calculated solutes and RMSD
462  US_ModelRecord*
463  cMr = &mrecs[ jcurve ];
464  cMr->csolutes = result.zsolutes;
465  cMr->variance = simulation_values.variance;
466  cMr->rmsd = cMr->variance > 0.0 ? sqrt( cMr->variance ) : 0.0;
467  int noiflg = simulation_values.noisflag;
468 
469  // Add in any noises
470  if ( noiflg > 0 )
471  {
472  if ( ( noiflg & 1 ) > 0 )
473  mrecs[ jcurve ].ti_noise = simulation_values.ti_noise;
474 
475  if ( ( noiflg & 2 ) > 0 )
476  mrecs[ jcurve ].ri_noise = simulation_values.ri_noise;
477  }
478 }
479 
480 // Write the model records file
481 void US_MPI_Analysis::write_mrecs()
482 {
483  US_Model * modelP = &data_sets[ current_dataset ]->model;
484  US_DataIO::EditedData*
485  edata = &data_sets[ current_dataset ]->run_data;
486  mrecs[ 0 ].modelGUID = modelP->modelGUID;
487 DbgLv(1) << "pcsa:wrmr: currds" << current_dataset;
488  mrecs[ 0 ].editGUID = edata->editGUID;
489  mrecs[ 0 ].mrecGUID = US_Util::new_guid();
490 DbgLv(1) << "pcsa:wrmr: editGUID=" << mrecs[0].editGUID
491  << "modelGUID=" << mrecs[0].modelGUID;
492  QString s_desc = QString( modelP->description )
493  .replace( ".model", ".mrecs" );
494  QString s_ctyp = parameters[ "curve_type" ];
495  int ctype = US_ModelRecord::ctype_flag( s_ctyp );
496  QString s_styp = parameters[ "solute_type" ];
497  int stype = US_ModelRecord::stype_flag( s_styp );
498  double x_min = parameters[ "x_min" ].toDouble();
499  double x_max = parameters[ "x_max" ].toDouble();
500  double y_min = parameters[ "y_min" ].toDouble();
501  double y_max = parameters[ "y_max" ].toDouble();
502  if ( ctype != CTYPE_ALL )
503  mrecs[ 0 ].ctype = ctype;
504  mrecs[ 0 ].xmin = x_min;
505  mrecs[ 0 ].xmax = x_max;
506  mrecs[ 0 ].ymin = y_min;
507  mrecs[ 0 ].ymax = y_max;
508  QString fnameo = s_desc + ".xml";
509  QFile fileo( fnameo );
510 
511  if ( fileo.open( QIODevice::WriteOnly | QIODevice::Text ) )
512  {
513  QXmlStreamWriter xmlo( &fileo );
514 
515 DbgLv(1) << "pcsa:wrmr: call wr_mr: s_desc" << s_desc;
516  US_ModelRecord::write_modelrecs( xmlo, mrecs, s_desc, ctype,
517  x_min, x_max, y_min, y_max, stype );
518 DbgLv(1) << "pcsa:wrmr: rtn wr_mr";
519 
520  fileo.close();
521  }
522 
523  // Add the file name of the mrecs file to the output list
524  QFile filea( "analysis_files.txt" );
525 
526  if ( ! filea.open( QIODevice::WriteOnly | QIODevice::Text
527  | QIODevice::Append ) )
528  {
529  abort( "Could not open 'analysis_files.txt' for writing" );
530  return;
531  }
532 
533  QTextStream tsout( &filea );
534  QString tripleID = edata->cell + edata->channel + edata->wavelength;
535  int mc_iter = mc_iterations;
536  int run = ( mc_iter > 1 ) ? mc_iter : 1;
537 
538  if ( mc_iterations > 0 )
539  run = mc_iter;
540 
541  QString runstring = "Run: " + QString::number( run ) + " " + tripleID;
542 
543  tsout << fnameo
544  << ";meniscus_value=" << meniscus_value
545  << ";MC_iteration=" << mc_iter
546  << ";variance=" << simulation_values.variance
547  << ";run=" << runstring
548  << "\n";
549  filea.close();
550 DbgLv(1) << "pcsa:wrmr: wr-to-atxt fnameo=" << fnameo;
551 }
552 
553 // Engineer gfit iterations for PCSA
554 void US_MPI_Analysis::iterate_pcsa()
555 {
556  static double rmsd_last = 0.0;
557  double rmsd_curr = mrecs[ 0 ].rmsd;
558  double thr_dr_rat = parameters[ "thr_deltr_ratio" ].toDouble();
559 
560 DbgLv(1) << "iter_p: iter" << iterations << "of" << max_iterations;
561  // If we have reached the maximum iterations, we are done
562  if ( iterations >= max_iterations ) return;
563 
564  if ( iterations > 1 )
565  { // If after first iteration, check RMSD difference
566  double dr_rat = qAbs( ( rmsd_curr - rmsd_last ) / rmsd_last );
567 DbgLv(1) << "iter_p: rmsd_c rmsd_l" << rmsd_curr << rmsd_last
568  << "dr_rat thr_dr_rat" << dr_rat << thr_dr_rat;
569 
570  // If difference in iteration RMSD less than threshold, we are done
571  if ( dr_rat < thr_dr_rat )
572  {
573  qDebug() << "Virtually identical RMSDs of"
574  << rmsd_last << rmsd_curr
575  << "truncate iterations at" << iterations;
576  max_iterations = iterations;
577  return;
578  }
579  }
580 
581  // Set up to create a narrower set of model record curves
582  double x_min = parameters[ "x_min" ].toDouble();
583  double x_max = parameters[ "x_max" ].toDouble();
584  double y_min = parameters[ "y_min" ].toDouble();
585  double y_max = parameters[ "y_max" ].toDouble();
586  double zval = parameters[ "z_value" ].toDouble();
587  QString s_ctyp = parameters[ "curve_type" ];
588  int ctype = US_ModelRecord::ctype_flag( s_ctyp );
589  QString s_styp = parameters[ "solute_type" ];
590  int stype = US_ModelRecord::stype_flag( s_styp );
591  int nkpts = parameters[ "vars_count" ].toInt();
592  int nlpts = parameters[ "curves_points" ].toInt();
593  double vbar20 = data_sets[ current_dataset ]->vbar20;
594  zval = ( zval == 0.0 ) ? vbar20 : zval;
595  double *parlims = (double*)pararry;
596  parlims[ 0 ] = -1.0;
597  parlims[ 4 ] = stype;
598  parlims[ 5 ] = zval;
599 DbgLv(1) << "iter_p: ctype" << ctype << s_ctyp;
600 
601  // Recompute mrecs to cover the range of the elite (top 10%) records
602 
603  if ( ctype != CTYPE_ALL )
604  { // All non-combination records: re-compute with narrower limits
605 
606  US_ModelRecord::recompute_mrecs( ctype, x_min, x_max, y_min, y_max,
607  nkpts, nlpts, parlims, mrecs );
608 DbgLv(1) << "iter_p: parlims"
609  << parlims[0] << parlims[1] << parlims[2] << parlims[3] << parlims[4]
610  << parlims[5] << parlims[6] << parlims[7] << parlims[8] << parlims[9];
611  } // END: non-combo records
612 
613  else
614  { // Records that are a mix of SL,IS,DS
615  QVector< US_ModelRecord > mrecs_sl;
616  QVector< US_ModelRecord > mrecs_is;
617  QVector< US_ModelRecord > mrecs_ds;
618  double* plims_sl = parlims;
619  double* plims_is = plims_sl + 12;
620  double* plims_ds = plims_is + 24;
621 
622  // Separate entries by curve type
623  filter_mrecs( CTYPE_SL, mrecs, mrecs_sl );
624  filter_mrecs( CTYPE_IS, mrecs, mrecs_is );
625  filter_mrecs( CTYPE_DS, mrecs, mrecs_ds );
626 
627  // Recompute each curve type with narrower limits
628  ctype = CTYPE_SL;
629  US_ModelRecord::recompute_mrecs( ctype, x_min, x_max, y_min, y_max,
630  nkpts, nlpts, plims_sl, mrecs_sl );
631 
632  ctype = CTYPE_IS;
633  US_ModelRecord::recompute_mrecs( ctype, x_min, x_max, y_min, y_max,
634  nkpts, nlpts, plims_is, mrecs_is );
635 
636  ctype = CTYPE_DS;
637  US_ModelRecord::recompute_mrecs( ctype, x_min, x_max, y_min, y_max,
638  nkpts, nlpts, plims_ds, mrecs_ds );
639 
640  // Re-combine into a single vector
641  int kk = 0;
642  int ncurve = sq( nkpts );
643 DbgLv(1) << "iter_p: ncurve" << ncurve;
644 
645  for ( int ii = 0; ii < ncurve; ii++, kk++ )
646  {
647  mrecs[ kk ] = mrecs_sl[ ii ];
648  mrecs[ kk ].taskx = kk;
649  }
650 
651  for ( int ii = 0; ii < ncurve; ii++, kk++ )
652  {
653  mrecs[ kk ] = mrecs_is[ ii ];
654  mrecs[ kk ].taskx = kk;
655  }
656 
657  for ( int ii = 0; ii < ncurve; ii++, kk++ )
658  {
659  mrecs[ kk ] = mrecs_ds[ ii ];
660  mrecs[ kk ].taskx = kk;
661  }
662 
663  ctype = CTYPE_ALL;
664  mrecs[0].v_ctype = ctype; // Vector curve type is "All"
665  plims_is[ 5 ] = parlims[ 5 ]; // Set z value for all types
666  plims_ds[ 5 ] = parlims[ 5 ];
667 DbgLv(1) << "iter_p: plims_sl"
668  << plims_sl[0] << plims_sl[1] << plims_sl[2] << plims_sl[3] << plims_sl[4];
669 DbgLv(1) << "iter_p: plims_is"
670  << plims_is[0] << plims_is[1] << plims_is[2] << plims_is[3] << plims_is[4];
671 DbgLv(1) << "iter_p: plims_ds"
672  << plims_ds[0] << plims_ds[1] << plims_ds[2] << plims_ds[3] << plims_ds[4];
673  } // END: records a mix of SL,IS,DS
674 
675  // Now reset original solutes and fill queue
676  orig_zsolutes.clear();
677  orig_zsolutes.reserve( mrecs.size() );
678 
679  for ( int ii = 0; ii < mrecs.size(); ii++ )
680  {
681  orig_zsolutes << mrecs[ ii ].isolutes;
682  }
683 
684  fill_pcsa_queue();
685 
686  for ( int ii = 1; ii <= my_workers; ii++ )
687  worker_status[ ii ] = READY;
688 
689  rmsd_last = rmsd_curr;
690  iterations++;
691 }
692 
693 // Engineer Tikhonov Regularization for PCSA
694 void US_MPI_Analysis::tikreg_pcsa()
695 {
696  int tikreg = parameters[ "tikreg_option" ].toInt();
697 DbgLv(1) << "tikr: tikreg" << tikreg;
698 
699  if ( tikreg == 0 )
700  return;
701 
702  alpha = ( tikreg == 1 )
703  ? parameters[ "tikreg_alpha" ].toDouble()
704  : alpha_scan();
705 
706  wksim_vals = simulation_values;
707  wksim_vals.zsolutes = mrecs[ 0 ].isolutes;
708  wksim_vals.alpha = alpha;
709  wksim_vals.noisflag = 0;
710  wksim_vals.ri_noise.clear();
711  wksim_vals.ri_noise.clear();
712 DbgLv(1) << "tikr: alpha" << alpha;
713 
714  calc_residuals( current_dataset, 1, wksim_vals );
715 
716  mrec = mrecs[ 1 ];
717  mrec.csolutes = wksim_vals.zsolutes;
718  mrec.variance = wksim_vals.variance;
719  mrec.rmsd = sqrt( mrec.variance );
720  mrec.sim_data = wksim_vals.sim_data;
721  mrec.residuals = wksim_vals.residuals;
722  mrecs[ 1 ] = mrec;
723 
724  write_pcsa_aux_model( 0 );
725 
726  qDebug() << "Tikhonov Regularization RMSD" << mrec.rmsd;
727 DbgLv(1) << "tikr: RMSD" << mrec.rmsd << "csol size" << mrec.csolutes.size();
728 
729  US_DataIO::EditedData* edata = &data_sets[ current_dataset ]->run_data;
730  QString tripleID = edata->cell + edata->channel + edata->wavelength;
731  QString progress = "Regularization: "
732  + QString::number( wksim_vals.alpha );
733 
734  if ( datasets_to_process > 1 )
735  progress += "; Datasets: "
736  + QString::number( datasets_to_process );
737  else
738  progress += "; Dataset: " + QString::number( current_dataset + 1 )
739  + " (" + tripleID + ") of "
740  + QString::number( count_datasets );
741 
742  progress += "; RMSD: " + QString::number( mrec.rmsd );
743 
744  send_udp( progress );
745 }
746 
747 // Engineer Monte Carlo for PCSA
748 void US_MPI_Analysis::montecarlo_pcsa()
749 {
750  int worker;
751  meniscus_value = data_sets[ current_dataset ]->run_data.meniscus;
752  mc_iterations = parameters[ "mc_iterations" ].toInt();
753  mc_iteration = 1;
754  mrec = mrecs[ 0 ];
755  int stype = mrec.stype;
756  US_DataIO::EditedData* edata = &data_sets[ current_dataset ]->run_data;
757  QString tripleID = edata->cell + edata->channel + edata->wavelength;
758  job_queue.clear();
759 
760  if ( mc_iterations < 2 ) return;
761 
762  // Fill the job queue and worker status
763  for ( int ii = 1; ii <= mc_iterations; ii++ )
764  {
765  Sa_Job job;
766  job.zsolutes = mrec.isolutes;
767  job.mpi_job.command = MPI_Job::PROCESS_MC;
768  job.mpi_job.length = job.zsolutes.size();
769  job.mpi_job.meniscus_value = meniscus_value;
770  job.mpi_job.dataset_offset = current_dataset;
771  job.mpi_job.dataset_count = datasets_to_process;
772  job.mpi_job.depth = ii;
773  job.mpi_job.solution = ii;
774 
775  job_queue << job;
776  }
777 
778  for ( int ii = 1; ii <= my_workers; ii++ )
779  {
780  worker_depth [ ii ] = ii;
781  worker_status[ ii ] = READY;
782  }
783 
784  // Build composite simulation,residuals data array
785  int nscan = edata->scanCount();
786  int npoint = edata->pointCount();
787  int mcd_incr = nscan * npoint;
788  int mcd_size = mcd_incr * 2;
789  int kci_send = 0;
790  int kci_recv = 0;
791  mc_data.resize( mcd_size );
792 
793  for ( int ii = 0; ii < nscan; ii++ )
794  {
795  int ks = ii * npoint;
796  int kr = ks + mcd_incr;
797 
798  for ( int jj = 0; jj < npoint; jj++, ks++, kr++ )
799  {
800  mc_data[ ks ] = mrec.sim_data .scanData[ ii ].rvalues[ jj ];
801  mc_data[ kr ] = mrec.residuals.scanData[ ii ].rvalues[ jj ];
802  }
803  }
804 
805  // Send simulation and residuals to the workers
806 
807  for ( worker = 1; worker <= my_workers; worker++ )
808  { // Tell each worker new data coming; worker expects a Send
809  MPI_Job mjob;
810  mjob.command = MPI_Job::NEWDATA;
811  mjob.length = mcd_size;
812  mjob.meniscus_value = meniscus_value;
813  mjob.dataset_offset = current_dataset;
814  mjob.dataset_count = datasets_to_process;
815  mjob.depth = worker;
816  mjob.solution = 10000;
817 
818  MPI_Send( &mjob,
819  sizeof( MPI_Job ),
820  MPI_BYTE,
821  worker,
822  MPI_Job::MASTER,
823  my_communicator );
824  }
825 
826  MPI_Barrier( my_communicator ); // Get everybody synced up
827 
828  MPI_Bcast( mc_data.data(), // Send simulation,residuals data
829  mcd_size,
830  MPI_DOUBLE,
831  MPI_Job::MASTER,
832  my_communicator );
833 
834  MPI_Barrier( my_communicator ); // Get everybody synced up
835 
836  for ( int ii = 1; ii < gcores_count; ii++ )
837  worker_status[ ii ] = INIT;
838 
839  worknext = 1;
840 
841  // Loop to submit and handle MC iteration jobs
842 
843  while ( true )
844  {
845  int mc_iter;
846  int mc_iters;
847  int sizes[ 4 ];
848  QString progress;
849  MPI_Status status;
850 
851  // Give the jobs to the workers
852  while ( ! job_queue.isEmpty() && worker_status.contains( READY ) )
853  {
854  kci_send++;
855  worker = ready_worker();
856 DbgLv(1) << " masterMC loop-ready_worker" << worker;
857  Sa_Job job = job_queue.takeFirst();
858 
859  job.mpi_job.command = MPI_Job::PROCESS_MC;
860  job.mpi_job.depth = kci_send;
861  job.mpi_job.solution = kci_send;
862 
863  // Tell worker that solutes are coming
864 DbgLv(1) << " masterMC loop-SEND PROCESS_MC worker" << worker
865  << "iter" << kci_send;
866  MPI_Send( &job.mpi_job,
867  sizeof( MPI_Job ),
868  MPI_BYTE,
869  worker,
870  MPI_Job::MASTER,
871  my_communicator );
872 
873  // Send solutes
874  MPI_Send( job.zsolutes.data(),
875  job.mpi_job.length * zsolut_doubles,
876  MPI_DOUBLE,
877  worker,
878  MPI_Job::MASTER,
879  my_communicator );
880 
881  worker_depth [ worker ] = job.mpi_job.depth;
882  worker_status[ worker ] = WORKING;
883  }
884 
885  // All done with the pass if no jobs are ready or running
886  if ( job_queue.isEmpty() && ! worker_status.contains( WORKING ) )
887  break; // Break out of MC iteration loop
888 
889  // Wait for worker to send a message
890  MPI_Recv( sizes,
891  4,
892  MPI_INT,
893  MPI_ANY_SOURCE,
894  MPI_ANY_TAG,
895  my_communicator,
896  &status );
897 
898  worker = status.MPI_SOURCE;
899 DbgLv(1) << " masterMC loop-BOTTOM: status TAG" << status.MPI_TAG
900  << " source" << status.MPI_SOURCE
901  << "sizes" << sizes[0] << sizes[1] << sizes[2] << sizes[3];
902  switch( status.MPI_TAG )
903  {
904  case MPI_Job::READY: // Ready for work
905 DbgLv(1) << " masterMC loop-RECV READY worker" << worker;
906  worker_status[ worker ] = READY;
907  break;
908 
909  case MPI_Job::RESULTS_MC: // Return solute data
910  kci_recv++;
911 DbgLv(1) << " masterMC loop-RECV RESULTS_MC worker" << worker << "kci_recv" << kci_recv;
912 
913  // Process PCSA MC results
914  wksim_vals = simulation_values;
915  wksim_vals.zsolutes .resize( sizes[ 0 ] );
916  wksim_vals.variances.resize( datasets_to_process );
917 
918  MPI_Recv( wksim_vals.zsolutes.data(),
919  sizes[ 0 ] * solute_doubles,
920  MPI_DOUBLE,
921  worker,
922  MPI_Job::TAG0,
923  my_communicator,
924  &status );
925 
926  MPI_Recv( &wksim_vals.variance,
927  1,
928  MPI_DOUBLE,
929  worker,
930  MPI_Job::TAG0,
931  my_communicator,
932  &status );
933 
934  MPI_Recv( wksim_vals.variances.data(),
935  datasets_to_process,
936  MPI_DOUBLE,
937  worker,
938  MPI_Job::TAG0,
939  my_communicator,
940  &status );
941 
942  work_rss[ worker ] = sizes[ 3 ];
943 
944  // Output the temporary MC iteration model file
945  mc_iter = worker_depth[ worker ];
946  mc_iters = mc_iterations;
947  mc_iterations++;
948 
949  write_pcsa_aux_model( mc_iter );
950 
951  mc_iterations = mc_iters;
952 DbgLv(1) << " masterMC loop-RECV RESULTS_MC mc_iter" << mc_iter << "of" << mc_iterations;
953 
954  // Send a status message
955  if ( datasets_to_process > 1 )
956  progress = "Datasets: "
957  + QString::number( datasets_to_process );
958  else
959  progress = "Dataset: "
960  + QString::number( current_dataset + 1 )
961  + " (" + tripleID + ") of "
962  + QString::number( count_datasets );
963 
964  progress += "; MonteCarlo: "
965  + QString::number( mc_iter );
966 
967  send_udp( progress );
968  break;
969 
970  case MPI_Job::RESULTS: // Return solute data (unused worker?)
971 DbgLv(1) << " masterMC loop-RECV RESULTS worker" << worker << "wstat" << worker_status[worker];
972  { // Receives just to consume messages
973  QVector< double > dwork;
974  int mxsiz = qMax( sizes[ 0 ], sizes[ 1 ] );
975  mxsiz = qMax( mxsiz, sizes[ 2 ] );
976  mxsiz = qMax( mxsiz, datasets_to_process );
977  dwork.resize( mxsiz );
978  double* wbuf = dwork.data();
979  MPI_Status status;
980 
981  MPI_Recv( wbuf,
982  sizes[ 0 ] * solute_doubles,
983  MPI_DOUBLE,
984  worker,
985  MPI_Job::TAG0,
986  my_communicator,
987  &status );
988  MPI_Recv( wbuf,
989  1,
990  MPI_DOUBLE,
991  worker,
992  MPI_Job::TAG0,
993  my_communicator,
994  &status );
995  MPI_Recv( wbuf,
996  datasets_to_process,
997  MPI_DOUBLE,
998  worker,
999  MPI_Job::TAG0,
1000  my_communicator,
1001  &status );
1002  MPI_Recv( wbuf,
1003  sizes[ 1 ],
1004  MPI_DOUBLE,
1005  worker,
1006  MPI_Job::TAG0,
1007  my_communicator,
1008  &status );
1009  MPI_Recv( wbuf,
1010  sizes[ 2 ],
1011  MPI_DOUBLE,
1012  worker,
1013  MPI_Job::TAG0,
1014  my_communicator,
1015  &status );
1016  }
1017  break;
1018 
1019  default: // Should never happen
1020 DbgLv(1) << " masterMC loop-RECV invalid worker" << worker << "tag" << status.MPI_TAG;
1021  progress = "Master PCSA: Received invalid status "
1022  + QString::number( status.MPI_TAG );
1023  abort( progress );
1024  break;
1025  }
1026 
1027  max_rss();
1028  }
1029 
1030  // Update output with MC model
1031  mc_iteration = mc_iterations;
1032 
1033  update_outputs();
1034 
1035  // Update the MC entry of the mrecs vector
1036  int mrx = ( mrecs[ 2 ].taskx == mrecs[ 0 ].taskx ) ? 2 : 1;
1037  mrec = mrecs[ mrx ];
1038  QString mfilt = "*" + tripleID + "*.mcN*model.xml";
1039  QStringList mfltl( mfilt );
1040  QStringList mlist = QDir( "." ).entryList( mfltl, QDir::Files );
1041  QString mfile;
1042  int nmfile = mlist.size();
1043 
1044  if ( nmfile == 1 )
1045  {
1046  mfile = mlist[ 0 ];
1047  }
1048  else if ( nmfile > 1 )
1049  {
1050  mfile = mlist[ 0 ];
1051  qDebug() << "*WARNING* More than 1" << mfilt << "file exists!";
1052  QDateTime mtime = QFileInfo( mfile ).lastModified();
1053 
1054  for ( int jj = 1; jj < mlist.size(); jj++ )
1055  { // Find the last of multiple files created
1056  QDateTime ftime = QFileInfo( mlist[ 0 ] ).lastModified();
1057  if ( ftime > mtime )
1058  {
1059  mtime = ftime;
1060  mfile = mlist[ 0 ];
1061  }
1062  }
1063  }
1064  else
1065  {
1066  qDebug() << "*ERROR* No" << mfilt << "file exists!";
1067  return;
1068  }
1069 
1070  wmodel.load( mfile );
1071  mrec.csolutes.resize( wmodel.components.size() );
1072 
1073  for ( int jj = 0; jj < wmodel.components.size(); jj++ )
1074  {
1075  US_ZSolute::set_solute_values( wmodel.components[ jj ],
1076  mrec.csolutes[ jj ], stype );
1077  }
1078 
1079  mrecs[ mrx ] = mrec;
1080 
1081  // Write the updated mrecs
1082  write_mrecs();
1083 }
1084 
1085 // Filter model records by a specified curve type
1086 void US_MPI_Analysis::filter_mrecs( const int ctype,
1087  QVector< US_ModelRecord >& mrecs_a, QVector< US_ModelRecord >& mrecs_t )
1088 {
1089  int namrec = mrecs_a.size();
1090  int ntmrec = namrec / 3;
1091  mrecs_t.clear();
1092  mrecs_t.reserve( ntmrec );
1093 
1094  for ( int ii = 0; ii < namrec; ii++ )
1095  if ( mrecs_a[ ii ].ctype == ctype )
1096  mrecs_t << mrecs_a[ ii ];
1097 }
1098 
1099 // Clean up model records to handle records with empty calculated solutes
1100 void US_MPI_Analysis::clean_mrecs( QVector< US_ModelRecord >& mrecs )
1101 {
1102  // First get index to last good
1103  int nmrec = mrecs.size();
1104  int lmrec = nmrec - 1;
1105  int goodx = 0;
1106 
1107  for ( int ii = lmrec; ii > 1; ii-- )
1108  {
1109  if ( mrecs[ ii ].csolutes.size() != 0 )
1110  {
1111  goodx = ii;
1112  break;
1113  }
1114  }
1115 
1116  // If last overall is good, return with no change to records
1117  if ( goodx == lmrec )
1118  return;
1119 
1120 DbgLv(1) << " master:clean_mrecs: goodx lmrec" << goodx << lmrec;
1121  if ( mrecs[ 0 ].v_ctype != CTYPE_ALL )
1122  { // If records of one type, just duplicate solutes from last good
1123  for ( int ii = goodx + 1; ii < nmrec; ii++ )
1124  {
1125  mrecs[ ii ].csolutes = mrecs[ goodx ].csolutes;
1126  mrecs[ ii ].rmsd = mrecs[ goodx ].rmsd;
1127  }
1128  }
1129 
1130  else
1131  { // If mixed records, duplicate solutes from last of each type
1132  int gx_sl = 0;
1133  int gx_is = 0;
1134  int gx_ds = 0;
1135  int ktype = 0;
1136 
1137  // Get indexes of last-good for each type
1138  for ( int ii = goodx; ii > 1; ii-- )
1139  {
1140  if ( mrecs[ ii ].ctype == CTYPE_SL )
1141  {
1142  if ( gx_sl > 0 ) continue;
1143  gx_sl = ii;
1144  ktype++;
1145  }
1146  else if ( mrecs[ ii ].ctype == CTYPE_IS )
1147  {
1148  if ( gx_is > 0 ) continue;
1149  gx_is = ii;
1150  ktype++;
1151  }
1152  else if ( mrecs[ ii ].ctype == CTYPE_DS )
1153  {
1154  if ( gx_ds > 0 ) continue;
1155  gx_ds = ii;
1156  ktype++;
1157  }
1158 
1159  if ( ktype >= 3 ) break;
1160  }
1161 
1162  // Then replace each type with last good of that type
1163  for ( int ii = goodx + 1; ii < nmrec; ii++ )
1164  {
1165  if ( mrecs[ ii ].ctype == CTYPE_SL )
1166  {
1167  mrecs[ ii ].csolutes = mrecs[ gx_sl ].csolutes;
1168  mrecs[ ii ].rmsd = mrecs[ gx_sl ].rmsd;
1169  }
1170  else if ( mrecs[ ii ].ctype == CTYPE_IS )
1171  {
1172  mrecs[ ii ].csolutes = mrecs[ gx_is ].csolutes;
1173  mrecs[ ii ].rmsd = mrecs[ gx_is ].rmsd;
1174  }
1175  else if ( mrecs[ ii ].ctype == CTYPE_DS )
1176  {
1177  mrecs[ ii ].csolutes = mrecs[ gx_ds ].csolutes;
1178  mrecs[ ii ].rmsd = mrecs[ gx_ds ].rmsd;
1179  }
1180  }
1181  }
1182 }
1183 
1184 // Save best pre-regularization/pre-montecarlo model and model record;
1185 // and insert place-holders for any regularization and/or montecarlo
1186 void US_MPI_Analysis::pcsa_best_model()
1187 {
1188  US_Model mdummy;
1189  mrec = mrecs[ 0 ];
1190  simulation_values.zsolutes = mrec.isolutes;
1191  simulation_values.ti_noise.clear();
1192  simulation_values.ri_noise.clear();
1193 
1194  calc_residuals( current_dataset, 1, simulation_values );
1195 
1196  mrec.csolutes = simulation_values.zsolutes;
1197  mrec.variance = simulation_values.variance;
1198  mrec.rmsd = sqrt( mrec.variance );
1199  mrec.model = data_sets[ current_dataset ]->model;
1200  mrec.sim_data = simulation_values.sim_data;
1201  mrec.residuals = simulation_values.residuals;
1202  mrec.ti_noise = simulation_values.ti_noise;
1203  mrec.ri_noise = simulation_values.ri_noise;
1204 DbgLv(1) << "bestm: RMSD" << mrec.rmsd;
1205 
1206  mrecs[ 0 ] = mrec;
1207 
1208  int tikreg = parameters[ "tikreg_option" ].toInt();
1209  mc_iterations = parameters[ "mc_iterations" ].toInt();
1210  int noiflg = simulation_values.noisflag;
1211 DbgLv(1) << "bestm: tikreg mciters" << tikreg << mc_iterations;
1212 
1213  if ( tikreg != 0 )
1214  {
1215  int nadd = ( mc_iterations > 1 ) ? 2 : 1;
1216  mrecs.insert( 1, nadd, mrec );
1217 
1218  mrecs[ 1 ].model = mdummy;
1219  mrecs[ 1 ].csolutes .clear();
1220  mrecs[ 1 ].ti_noise .clear();
1221  mrecs[ 1 ].ri_noise .clear();
1222  mrecs[ 1 ].mrecGUID .clear();
1223  mrecs[ 1 ].mrecGUID .clear();
1224  mrecs[ 1 ].modelGUID.clear();
1225 
1226  if ( nadd > 1 )
1227  {
1228  mrecs[ 2 ].model = mdummy;
1229  mrecs[ 2 ].csolutes .clear();
1230  mrecs[ 2 ].ti_noise .clear();
1231  mrecs[ 2 ].ri_noise .clear();
1232  mrecs[ 2 ].mrecGUID .clear();
1233  mrecs[ 2 ].mrecGUID .clear();
1234  mrecs[ 2 ].modelGUID.clear();
1235  }
1236  }
1237 
1238  else if ( mc_iterations > 1 )
1239  {
1240  mrecs.insert( 1, 1, mrec );
1241  mrecs[ 1 ].model = mdummy;
1242  mrecs[ 1 ].csolutes .clear();
1243  mrecs[ 1 ].ti_noise .clear();
1244  mrecs[ 1 ].ri_noise .clear();
1245  mrecs[ 1 ].mrecGUID .clear();
1246  mrecs[ 1 ].mrecGUID .clear();
1247  mrecs[ 1 ].modelGUID.clear();
1248  }
1249 
1250  if ( ( tikreg != 0 || mc_iterations > 1 ) &&
1251  ( noiflg != 0 ) )
1252  { // Apply computed noise to data
1253  US_DataIO::EditedData* edata = &data_sets[ current_dataset ]->run_data;
1254  double vnoise = 0.0;
1255  int nscan = edata->scanCount();
1256  int npoint = edata->pointCount();
1257  int rcount = mrec.ti_noise.size();
1258  int scount = mrec.ri_noise.size();
1259 DbgLv(1) << "bestm: noise apply ns np sc rc" << nscan << npoint << scount << rcount;
1260  npoint = ( noiflg & 1 ) > 0 ? npoint : rcount;
1261  nscan = ( noiflg & 2 ) > 0 ? nscan : scount;
1262 
1263  if ( rcount != npoint || scount != nscan )
1264  {
1265  qDebug() << "*ERROR* noise count(s) do not match data dimensions!";
1266  return;
1267  }
1268 
1269  if ( ( noiflg & 1 ) > 0 )
1270  { // Apply ti noise
1271 DbgLv(1) << "bestm: noise flag" << noiflg << "TI_NOISE apply";
1272  for ( int jj = 0; jj < rcount; jj++ )
1273  { // Get constant noise value for each reading and apply to scans
1274  vnoise = -1.0 * mrec.ti_noise[ jj ];
1275  for ( int ii = 0; ii < scount; ii++ )
1276  { // Apply to all scans at reading position
1277  edata->scanData[ ii ].rvalues[ jj ] += vnoise;
1278  }
1279  }
1280  }
1281 
1282  if ( ( noiflg & 2 ) > 0 )
1283  { // Apply ri noise
1284 DbgLv(1) << "bestm: noise flag" << noiflg << "RI_NOISE apply";
1285  for ( int ii = 0; ii < scount; ii++ )
1286  { // Get constant noise value for each scan and apply to readings
1287  vnoise = -1.0 * mrec.ri_noise[ ii ];
1288  for ( int jj = 0; jj < rcount; jj++ )
1289  { // Apply to all readings at scan position
1290  edata->scanData[ ii ].rvalues[ jj ] += vnoise;
1291  }
1292  }
1293  }
1294  }
1295 }
1296 
1297 // Write a PCSA auxiliary (TR/MC) model and potentially model records
1298 void US_MPI_Analysis::write_pcsa_aux_model( int iter )
1299 {
1300  wmodel = mrecs[ 0 ].model;
1301  wmodel.modelGUID = US_Util::new_guid();
1302  wmodel.variance = wksim_vals.variance;
1303  QString mdesc = wmodel.description;
1304  QString runID = QString( mdesc ).section( ".", 0, -4 );
1305  QString tripID = QString( mdesc ).section( ".", -3, -3 );
1306  QString asysID = QString( mdesc ).section( ".", -2, -2 );
1307  QString typeExt = QString( ".model" );
1308  QString dates = QString( asysID ).section( "_", 0, 1 );
1309  QString atype = QString( asysID ).section( "_", 2, 2 );
1310  QString reqID = QString( asysID ).section( "_", 3, 3 );
1311  QString iterID = "i01" ;
1312  if ( iter == 0 )
1313  {
1314  atype += "-TR";
1315  wmodel.alphaRP = alpha;
1316  mrecs[ 1 ].modelGUID = wmodel.modelGUID;
1317  }
1318  else
1319  {
1320  atype += "-MC";
1321  wmodel.monteCarlo = true;
1322  iterID = QString().sprintf( "mc%04d", iter );
1323  int jj = ( mrecs[ 2 ].taskx == mrecs[ 0 ].taskx ) ? 2 : 1;
1324  mrecs[ jj ].modelGUID = wmodel.modelGUID;
1325  }
1326  asysID = dates + "_" + atype + "_" + reqID + "_" + iterID;
1327  mdesc = runID + "." + tripID + "." + asysID + typeExt;
1328  wmodel.description = mdesc;
1329  wmodel.components.clear();
1330 DbgLv(1) << "wraux: mdesc" << mdesc;
1331  QString s_styp = parameters[ "solute_type" ];
1332  int stype = US_ModelRecord::stype_flag( s_styp );
1333 
1334  for ( int ii = 0; ii < wksim_vals.zsolutes.size(); ii++ )
1335  {
1336  US_Model::SimulationComponent component;
1337  US_ZSolute::set_mcomp_values( component, wksim_vals.zsolutes[ ii ],
1338  stype, true );
1339  component.name = QString().sprintf( "SC%04d", ii + 1 );
1340 
1341  US_Model::calc_coefficients( component );
1342 
1343  wmodel.components << component;
1344  }
1345 
1346  mrec.model = wmodel;
1347  QString fext = ( iter == 0 ) ? ".model.xml" : ".mdl.tmp";
1348  QString fileid = "." + atype + "." + tripID + "." + iterID + fext;
1349  QString fn = runID + fileid;
1350  int lenfn = fn.length();
1351 
1352  if ( lenfn > 99 )
1353  {
1354  int lenri = runID.length() + 99 - lenfn;
1355  fn = QString( runID ).left( lenri ) + fileid;
1356  }
1357 
1358  // Output the model to a file
1359  wmodel.write( fn );
1360 
1361  // Add the model file name to the output list
1362  QFile fileo( "analysis_files.txt" );
1363 
1364  if ( ! fileo.open( QIODevice::WriteOnly | QIODevice::Text
1365  | QIODevice::Append ) )
1366  {
1367  abort( "Could not open 'analysis_files.txt' for writing" );
1368  return;
1369  }
1370 
1371  QTextStream tsout( &fileo );
1372 
1373  tsout << fn
1374  << ";meniscus_value=" << meniscus_value
1375  << ";MC_iteration=" << iter
1376  << ";variance=" << wmodel.variance
1377  << ";run=Run: 1 " << tripID << "\n";
1378 
1379  fileo.close();
1380 
1381  // If final regularization or final MC iter model, write model records
1382  if ( ( iter == 0 && mc_iterations < 2 ) ||
1383  ( iter > 0 && iter >= mc_iterations ) )
1384  {
1385  write_mrecs();
1386  }
1387 }
1388 
1389 // Scan alphas
1390 double US_MPI_Analysis::alpha_scan()
1391 {
1392  const double salpha = 0.10;
1393  const double ealpha = 0.90;
1394  const double dalpha = 0.01;
1395  const double roundv = dalpha * 0.05;
1396  const double alphadef = salpha + ( ealpha - salpha ) * 0.25;
1397  int nalpha = qRound( ( ealpha - salpha ) / dalpha ) + 1;
1398  QVector< double > alphas;
1399  QVector< double > varias;
1400  QVector< double > xnorms;
1401  QVector< double > sv_nnls_a;
1402  QVector< double > sv_nnls_b;
1403  alphas.clear();
1404  varias.clear();
1405  xnorms.clear();
1406  alphas.reserve( nalpha );
1407  varias.reserve( nalpha );
1408  xnorms.reserve( nalpha );
1409 
1410  double varmx = 0.0;
1411  double xnomx = 0.0;
1412  double v_vari = 0.0;
1413  double v_xnsq = 0.0;
1414  double calpha = salpha;
1415 
1416  for ( int ja = 0; ja < nalpha; ja++ )
1417  { // Loop to populate alpha scan arrays with preliminary values
1418  alphas << calpha;
1419  varias << v_vari;
1420  xnorms << v_xnsq;
1421  calpha += dalpha;
1422  }
1423 
1424  // Do alpha scan
1425  QVector< double > csolutes;
1426  sv_nnls_a.clear();
1427  sv_nnls_b.clear();
1428  mrec = mrecs[ 0 ];
1429  int nscans = data_sets[ current_dataset ]->run_data.scanCount();
1430  int npoints = data_sets[ current_dataset ]->run_data.pointCount();
1431  int nisols = mrec.isolutes.size();
1432 
1433  for ( int ja = 0; ja < nalpha; ja++ )
1434  { // Loop to evaluate each alpha in the specified range
1435  calpha = alphas[ ja ];
1436 
1437  if ( ja == 0 )
1438  { // For the first one do a full compute and save the A,B matrices
1439  US_SolveSim::Simulation sim_vals;
1440  sim_vals.alpha = calpha;
1441  sim_vals.zsolutes = mrec.isolutes;
1442 
1443  US_SolveSim* solvesim = new US_SolveSim( data_sets, 0, false );
1444 
1445  solvesim->calc_residuals( current_dataset, 1, sim_vals, true,
1446  &sv_nnls_a, &sv_nnls_b );
1447 
1448  v_vari = sim_vals.variance;
1449  v_xnsq = sim_vals.xnormsq;
1450  }
1451 
1452  else
1453  { // After 1st, regularize by modifying A matrix diagonal, then NNLS
1454  apply_alpha( calpha, &sv_nnls_a, &sv_nnls_b,
1455  nscans, npoints, nisols, v_vari, v_xnsq );
1456  }
1457 
1458  varias[ ja ] = v_vari;
1459  xnorms[ ja ] = v_xnsq;
1460  varmx = qMax( varmx, v_vari );
1461  xnomx = qMax( xnomx, v_xnsq );
1462 DbgLv(1) << "a v x" << calpha << v_vari << v_xnsq;
1463  }
1464 
1465  sv_nnls_a.clear();
1466  sv_nnls_b.clear();
1467 
1468  int lgv = 0 - (int)qFloor( log10( varmx ) );
1469  int lgx = -1 - (int)qFloor( log10( xnomx ) );
1470  double vscl = qPow( 10.0, lgv );
1471  double xscl = qPow( 10.0, lgx );
1472 DbgLv(1) << "Log-varia Log-xnorm" << lgv << lgx << "vscl xscl" << vscl << xscl;
1473 
1474  for ( int ja = 0; ja < nalpha; ja++ )
1475  { // Scale the variance,normsq points
1476  varias[ ja ] *= vscl;
1477  xnorms[ ja ] *= xscl;
1478  }
1479 
1480  // Compute lines that hint at the elbow point of the curve
1481  double* xx = varias.data();
1482  double* yy = xnorms.data();
1483  double xa[ 5 ];
1484  double ya[ 5 ];
1485  double* xe = (double*)xa;
1486  double* ye = (double*)ya;
1487  double slope; double slop2;
1488  double intcp; double intc2;
1489  double sigma; double corre;
1490  double xlipt; double ylipt; double xcipt; double ycipt;
1491  int nlp = 5;
1492 
1493  // Compute a line fitted to the first few main curve points
1494  while ( nlp < nalpha )
1495  {
1496  double avg = US_Math2::linefit( &xx, &yy, &slope, &intcp,
1497  &sigma, &corre, nlp );
1498 DbgLv(1) << "ASCN:H1: avg" << avg << "nlp" << nlp;
1499 DbgLv(1) << "ASCN:H1: sl" << slope << "in" << intcp << "sg" << sigma
1500  << "co" << corre;
1501  if ( slope < 1e99 && slope > (-1e99) )
1502  break;
1503 
1504  nlp += 2;
1505  }
1506 
1507  // Compute a line fitted to the last few main curve points
1508  int je = nalpha - 1;
1509  for ( int jj = 0; jj < 5; jj++, je-- )
1510  {
1511  xe[ jj ] = xx[ je ];
1512  ye[ jj ] = yy[ je ];
1513  }
1514 
1515  US_Math2::linefit( &xe, &ye, &slop2, &intc2, &sigma, &corre, 5 );
1516 
1517  // Find the intersection point for the 2 fitted lines
1518  US_Math2::intersect( slope, intcp, slop2, intc2, &xlipt, &ylipt );
1519 
1520  // Find the curve point nearest to the intersection point;
1521  // then compute a line from intersection to nearest curve point.
1522  US_Math2::nearest_curve_point( xx, yy, nalpha, true, xlipt, ylipt,
1523  &xcipt, &ycipt, alphas.data(), &calpha );
1524 
1525  // Do a sanity check. If the intersection point is outside the
1526  // rectangle that encloses the curve, we likely have an aberrant curve.
1527  // So, forget elbow fit and default alpha.
1528  double xcvp1 = xx[ 0 ];
1529  double ycvp1 = yy[ 0 ];
1530  double xcvp2 = xx[ je ];
1531  double ycvp2 = yy[ je ];
1532  bool good_fit = ( xlipt >= xcvp1 && ylipt <= ycvp1 &&
1533  xlipt <= xcvp2 && ylipt >= ycvp2 );
1534 DbgLv(1) << "ASCN:T4: cv: x1,y1" << xcvp1 << ycvp1
1535  << "x2,y2" << xcvp2 << ycvp2 << " good_fit" << good_fit;
1536 
1537  calpha = (double)qRound( calpha / roundv ) * roundv;
1538  calpha = good_fit ? calpha : alphadef;
1539 
1540  return calpha;
1541 }
1542 
1543 void US_MPI_Analysis::apply_alpha( const double alpha, QVector< double >* psv_nnls_a,
1544  QVector< double >* psv_nnls_b, const int nscans, const int npoints,
1545  const int nisols, double& variance, double& xnormsq )
1546 {
1547  int ntotal = nscans * npoints;
1548  int narows = ntotal + nisols;
1549  int ncsols = 0;
1550  variance = 0.0;
1551  xnormsq = 0.0;
1552  QVector< double > nnls_a = *psv_nnls_a;
1553  QVector< double > nnls_b = *psv_nnls_b;
1554  QVector< double > nnls_x;
1555  QVector< double > simdat;
1556  nnls_x.fill( 0.0, nisols );
1557  simdat.fill( 0.0, ntotal );
1558 
1559  // Replace alpha in the diagonal of the lower square of A
1560  int dx = ntotal;
1561  int dinc = ntotal + nisols + 1;
1562 
1563  for ( int cc = 0; cc < nisols; cc++ )
1564  {
1565  nnls_a[ dx ] = alpha;
1566  dx += dinc;
1567  }
1568 
1569  // Compute the X vector using NNLS
1570  US_Math2::nnls( nnls_a.data(), narows, narows, nisols,
1571  nnls_b.data(), nnls_x.data() );
1572 
1573  // Construct the output solutes and the implied simulation and xnorm-sq
1574  for ( int cc = 0; cc < nisols; cc++ )
1575  {
1576  double soluval = nnls_x[ cc ]; // Computed concentration, this solute
1577 
1578  if ( soluval > 0.0 )
1579  {
1580  xnormsq += sq( soluval );
1581  ncsols++;
1582  int aa = cc * narows;
1583 
1584  for ( int kk = 0; kk < ntotal; kk++ )
1585  {
1586  simdat[ kk ] += ( soluval * (*psv_nnls_a)[ aa++ ] );
1587  }
1588  }
1589  }
1590 
1591  // Calculate the sum for the variance computation
1592  for ( int kk = 0; kk < ntotal; kk++ )
1593  {
1594  variance += sq( ( (*psv_nnls_b)[ kk ] - simdat[ kk ] ) );
1595  }
1596 
1597  // Return computed variance and xnorm-sq
1598  variance /= (double)ntotal;
1599 }
1600 
1601 // Fill the job queue, using the list of initial solutes
1602 void US_MPI_Analysis::fill_pcsa_queue( void )
1603 {
1604  worker_status.resize( gcores_count );
1605  worker_depth .resize( gcores_count );
1606 
1607  worker_status.fill( INIT );
1608  worker_depth .fill( 0 );
1609  max_depth = 0;
1610  worknext = 1;
1611  max_experiment_size = min_experiment_size;
1612 
1613  // Put all jobs in the queue
1614  job_queue.clear();
1615 
1616  for ( int i = 0; i < orig_zsolutes.size(); i++ )
1617  {
1618  max_experiment_size = qMax( max_experiment_size,
1619  orig_zsolutes[ i ].size() );
1620  Sa_Job job;
1621  job.zsolutes = orig_zsolutes[ i ];
1622  job_queue << job;
1623  }
1624 }
1625 
Generated on Sun Jun 21 2015 11:40:40 for UltraScan III by   doxygen 1.8.3.1-20130324