UltraScan III
us_solve_sim.cpp
Go to the documentation of this file.
1 #include "us_solve_sim.h"
3 #include "us_util.h"
4 #include "us_settings.h"
5 #include "us_math2.h"
6 #include "us_constants.h"
7 #include "us_memory.h"
8 
9 // Define level-conditioned debug print that includes thread/processor
10 #ifdef DbgLv
11 #undef DbgLv
12 #define DbgLv(a) if(dbg_level>=a)qDebug()<<"SS-w:"<<thrnrank<<":"
14 #endif
15 
16 double zerothr = 0.020;
17 double linethr = 0.050;
18 double maxod = 1.50;
19 double mfactor = 3.00;
20 double mfactex = 1.00;
21 double minnzsc = 0.005;
22 
23 // Create a Solve-Simulation object
24 US_SolveSim::US_SolveSim( QList< DataSet* >& data_sets, int thrnrank,
25  bool signal_wanted ) : QObject(), data_sets( data_sets ),
26  thrnrank( thrnrank ), signal_wanted( signal_wanted )
27 {
28  abort = false; // Default: no abort
29  dbg_level = 0; // Default: no debug prints
30  dbg_timing = false; // Default: no debug timing prints
31  banddthr = false; // Default: no bandform data_threshold
32 
33  // If band-forming, possibly read in threshold control values
34  if ( data_sets[ 0 ]->simparams.band_forming )
35  {
36  QString bfcname = US_Settings::appBaseDir() + "/etc/bandform.config";
37  QFile bfcfile( bfcname );
38  if ( bfcfile.open( QIODevice::ReadOnly ) )
39  {
40  QTextStream tsi( &bfcfile );
41  // Read in, in order: zero-thresh, linear-thresh, max-od, min-nz-scale
42  while ( ! tsi.atEnd() )
43  {
44  QString fline = tsi.readLine();
45  int cmx = fline.indexOf( "#" ) - 1;
46  double val = cmx > 0
47  ? fline.left( cmx ).simplified().toDouble() : 0.0;
48  if ( cmx <= 0 )
49  continue;
50  else if ( fline.contains( "zero threshold" ) )
51  zerothr = val;
52  else if ( fline.contains( "linear threshold" ) )
53  linethr = val;
54  else if ( fline.contains( "maximum OD" ) )
55  maxod = val;
56  else if ( fline.contains( "minimum non-zero" ) )
57  minnzsc = val;
58  else if ( fline.contains( "sim multiply factor" ) )
59  mfactor = val;
60  else if ( fline.contains( "exp multiply factor" ) )
61  mfactex = val;
62  }
63  bfcfile.close();
64 
65  banddthr = true;
66  }
67 if(thrnrank==1) DbgLv(1) << "CR:zthr lthr mxod mnzc mfac mfex bthr"
68  << zerothr << linethr << maxod << minnzsc << mfactor << mfactex << banddthr;
69  }
70 }
71 
72 // Create a Simulation object
73 US_SolveSim::Simulation::Simulation()
74 {
75  variance = 0.0;
76  xnormsq = 0.0;
77  alpha = 0.0;
78  variances.clear();
79  ti_noise .clear();
80  ri_noise .clear();
81  solutes .clear();
82  zsolutes .clear();
83  dbg_level = 0;
84  dbg_timing = false;
85  noisflag = 0;
86 }
87 
88 // Static function to check the grid size implied by data and model
89 bool US_SolveSim::checkGridSize( QList< DataSet* >& data_sets,
90  double s_max, QString& smsg )
91 {
92  const long tstep_max = 10000L;
93  bool too_large = false;
94  smsg = QString( "" );
95  US_SimulationParameters* sparams = &data_sets[ 0 ]->simparams;
96  double meniscus = sparams->meniscus;
97  double bottom = sparams->bottom;
98  int rpm_max = sparams->speed_step[ 0 ].rotorspeed;
99  int time_max = 0;
100  int rsimpts = sparams->simpoints - 1;
101 
102  for ( int dd = 0; dd < data_sets.size(); dd++ )
103  {
104  sparams = &data_sets[ dd ]->simparams;
105  QVector< US_SimulationParameters::SpeedProfile > speed_step
106  = sparams->speed_step;
107 
108  for ( int ss = 0; ss < speed_step.size(); ss++ )
109  {
110  int duration = qRound( speed_step[ ss ].duration_hours * 3600.0
111  + speed_step[ ss ].duration_minutes * 60.0 );
112  time_max = qMax( time_max, duration );
113  rpm_max = qMax( rpm_max, speed_step[ ss ].rotorspeed );
114  }
115  }
116 
117  double s_show = s_max * 1.0e13;
118  double omega_s = sq( rpm_max * M_PI / 30.0 );
119  double lgbmrat = log( bottom / meniscus );
120  double somgfac = s_max * omega_s;
121  double dt = lgbmrat / ( somgfac * rsimpts );
122  long tsteps = (long)( time_max / dt ) + 1L;
123 //if(thrnrank==1) DbgLv(1) << "CK: tsteps" << tsteps << "dt omg rpm" << dt << omega_s << rpm_max
124 // << "bot men rat sfac spts" << bottom << meniscus << lgbmrat << somgfac << rsimpts;
125 
126  if ( tsteps > tstep_max )
127  {
128  too_large = true;
129  smsg = tr( "The combination of rotor speed, the ratio of sedimentation\n"
130  "over diffusion coefficient, and length of experiment\n"
131  "caused the program to exceed available memory.\n"
132  "Please sufficiently reduce any of these to bring\n"
133  "the program into a feasible range.\n\n"
134  "Related data and simulation parameters include:\n"
135  " Maximum speed = %1 RPM ;\n"
136  " Maximum S value = %2 x 1e-13 ;\n"
137  " Computed grid time steps and radius points = %3 ;\n"
138  " Program-imposed grid time steps upper limit = %4 ." )
139  .arg( rpm_max ).arg( s_show ).arg( tsteps ).arg( tstep_max );
140  }
141 
142  return too_large;
143 }
144 
145 // Check the grid size implied by data and model (class method)
146 bool US_SolveSim::check_grid_size( double s_max, QString& smsg )
147 {
148  return checkGridSize( data_sets, s_max, smsg );
149 }
150 
151 // Do the real work of a 2dsa/ga thread/processor: simulation from solutes set
152 void US_SolveSim::calc_residuals( int offset, int dataset_count,
153  Simulation& sim_vals, bool padAB,
154  QVector< double >* ASave, QVector< double >* BSave )
155 {
156  QVector< double > sv_nnls_a;
157  QVector< double > sv_nnls_b;
158  dbg_level = sim_vals.dbg_level; // Debug level
159  dbg_timing = sim_vals.dbg_timing; // Debug-timing flag
160  noisflag = sim_vals.noisflag; // Noise-calculation flag
161  int nsolutes = sim_vals.solutes.size(); // Input solutes count
162  int nzsol = sim_vals.zsolutes.size(); // Input zsolutes count
163  calc_ti = ( ( noisflag & 1 ) != 0 ); // Calculate-TI flag
164  calc_ri = ( ( noisflag & 2 ) != 0 ); // Calculate-RI flag
165  startCalc = QDateTime::currentDateTime(); // Start calc time for timings
166  int ntotal = 0; // Total points count
167  int ntinois = 0; // TI noise value count
168  int nrinois = 0; // RI noise value count
169  double alphad = sim_vals.alpha; // Alpha for diagonal
170  int lim_offs = offset + dataset_count; // Offset limit
171  d_offs = offset; // Initial data offset
172  bool use_zsol = ( nzsol > 0 ); // Flag use ZSolutes
173  nsolutes = use_zsol ? nzsol : nsolutes; // Count of used solutes
174 #if 0
175 #ifdef NO_DB
176  US_Settings::set_us_debug( dbg_level );
177 #endif
178 #endif
179 
180  if ( banddthr )
181  {
182  mfactor = data_sets[ 0 ]->simparams.cp_width;
183  if ( mfactor < 0.0 )
184  {
185  mfactor = 1.0;
186  zerothr = 0.0;
187  linethr = 0.0;
188  maxod = 1.0e+20;
189  }
190 if(thrnrank==1) DbgLv(1) << "CR:zthr lthr mxod mfac"
191  << zerothr << linethr << maxod << mfactor;
192  }
193 
194  US_DataIO::EditedData wdata;
195 
196  for ( int ee = offset; ee < lim_offs; ee++ )
197  { // Count scan,point totals for all data sets
198  DataSet* dset = data_sets[ ee ];
199  int npoints = dset->run_data.pointCount();
200  int nscans = dset->run_data.scanCount();
201  ntotal += ( nscans * npoints );
202  ntinois += npoints;
203  nrinois += nscans;
204  }
205 
206  bool tikreg = ( sim_vals.alpha != 0.0 );
207  int narows = tikreg ? ( ntotal + nsolutes ) : ntotal;
208 
209  // Set up and clear work vectors
210  int navals = narows * nsolutes; // Size of "A" matrix
211 //DbgLv(1) << " CR:na nb nx" << navals << ntotal << nsolutes;
212 
213 #if 1
214  QVector< double > nnls_a( navals, 0.0 );
215 #endif
216 #if 0
217  QVector< double > nnls_a;
218  int iasize = navals;
219  if ( navals > 20000000 )
220  {
221  iasize = narows * ( nsolutes / 2 );
222  }
223  nnls_a.reserve( iasize );
224  nnls_a.fill( 0.0, iasize );
225 #endif
226  QVector< double > nnls_b( narows, 0.0 );
227  QVector< double > nnls_x( nsolutes, 0.0 );
228  QVector< double > tinvec( ntinois, 0.0 );
229  QVector< double > rinvec( nrinois, 0.0 );
230 
231  // Set up for a single-component model
232  US_Model model;
233  model.components.resize( 1 );
234 
235  US_Model::SimulationComponent zcomponent; // Zeroed component to init models
236  zcomponent.s = 0.0;
237  zcomponent.D = 0.0;
238  zcomponent.mw = 0.0;
239  zcomponent.f = 0.0;
240  zcomponent.f_f0 = 0.0;
241 
242  if ( banddthr )
243  { // If band forming, hold data within thresholds
244 //mfactex=mfactor;
245  wdata = data_sets[ d_offs ]->run_data;
246  data_threshold( &wdata, zerothr, linethr, maxod, mfactex );
247  }
248 
249 DebugTime("BEG:calcres");
250  // Populate b array with experiment data concentrations
251  int kk = 0;
252  int kodl = 0;
253 #if 0
254  double s0max = 0.0;
255 #endif
256 
257  for ( int ee = offset; ee < lim_offs; ee++ )
258  {
259  US_DataIO::EditedData* edata = &data_sets[ ee ]->run_data;
260  edata = banddthr ? &wdata : edata;
261  int npoints = edata->pointCount();
262  int nscans = edata->scanCount();
263  double odlim = edata->ODlimit;
264 
265  for ( int ss = 0; ss < nscans; ss++ )
266  {
267  for ( int rr = 0; rr < npoints; rr++ )
268  { // Fill the B matrix with experiment data (or zero)
269  double evalue = edata->value( ss, rr );
270 
271  if ( evalue >= odlim )
272  { // Where ODlimit is exceeded, substitute 0.0 and bump count
273  evalue = 0.0;
274  kodl++;
275  }
276 
277 #if 0
278  else if ( ss == 0 )
279  { // Find max of scan 0 values
280  s0max = qMax( s0max, evalue );
281  }
282 #endif
283 
284  nnls_b[ kk++ ] = evalue;
285  }
286  }
287  }
288 DbgLv(1) << " CR:B fill kodl" << kodl;
289 
290 #if 0
291  // If needed, scale the alpha used in A-matrix appendix diagonals
292  if ( tikreg )
293  {
294  alphad = ( s0max == 0.0 ) ? sim_vals.alpha
295  : ( sim_vals.alpha * sqrt( s0max ) );
296  }
297 #endif
298 
299  if ( abort ) return;
300 
301  QList< US_DataIO::RawData > simulations;
302  simulations.reserve( nsolutes * dataset_count );
303 
304  // Simulate data using models, each with a single s,f/f0 component
305  int increp = nsolutes / 10; // Progress report increment
306  increp = ( increp < 10 ) ? 10 : increp;
307  int kstep = 0; // Progress step count
308  kk = 0; // nnls_a output index
309  int ksols = 0;
310  int stype = data_sets[ offset ]->solute_type;
311 DbgLv(1) << " CR:BF STYPE" << stype;
312 
313  if ( use_zsol )
314  qSort( sim_vals.zsolutes );
315  else
316  qSort( sim_vals.solutes );
317 
318  if ( stype == 0 )
319  { // Normal case of varying f/f0 with constant vbar
320  int attr_x = 0; // Default X is s
321  int attr_y = 1; // Default Y is f/f0
322  int attr_z = 3; // Default Z is vbar
323  int smask = ( attr_x << 6 ) | ( attr_y << 3 ) | attr_z;
324 DataSet* dset=data_sets[0];
325 DbgLv(2) << " CR:BF s20wcorr D20wcorr" << dset->s20w_correction
326  << dset->D20w_correction << "manual" << dset->solution_rec.buffer.manual
327  << "vbar20" << dset->vbar20;
328 
329  for ( int cc = 0; cc < nsolutes; cc++ )
330  { // Solve for each solute
331  if ( abort ) return;
332  int bx = 0;
333 
334  for ( int ee = offset; ee < lim_offs; ee++ )
335  { // Solve for each data set
336  DataSet* dset = data_sets[ ee ];
337  US_DataIO::EditedData* edata = banddthr ? &wdata : &dset->run_data;
338  US_DataIO::RawData simdat;
339  int npoints = edata->pointCount();
340  int nscans = edata->scanCount();
341  zcomponent.vbar20 = dset->vbar20;
342 
343  // Set model with standard space s and k
344  if ( use_zsol )
345  {
346  US_ZSolute::set_mcomp_values( model.components[ 0 ],
347  sim_vals.zsolutes[ cc ], smask );
348  }
349  else
350  {
351  model.components[ 0 ] = zcomponent;
352  set_comp_attr( model.components[ 0 ], sim_vals.solutes[ cc ],
353  attr_x );
354  set_comp_attr( model.components[ 0 ], sim_vals.solutes[ cc ],
355  attr_y );
356  }
357 
358  // Fill in the missing component values
359  model.update_coefficients();
360 //DbgLv(1) << "CR: cc" << cc << "model s,k,w,v,d,c"
361 // << model.components[0].s << model.components[0].f_f0
362 // << model.components[0].mw << model.components[0].vbar20
363 // << model.components[0].D << model.components[0].signal_concentration;
364 
365  // Convert to experimental space
366  model.components[ 0 ].s /= dset->s20w_correction;
367  model.components[ 0 ].D /= dset->D20w_correction;
368 //DbgLv(1) << "CR: exp.space model s,k,w,v,d,c"
369 // << model.components[0].s << model.components[0].f_f0
370 // << model.components[0].mw << model.components[0].vbar20
371 // << model.components[0].D << model.components[0].signal_concentration;
372 
373  // Initialize simulation data with the experiment's grid
374 DbgLv(2) << " CR:111 rss now" << US_Memory::rss_now() << "cc" << cc;
375  US_AstfemMath::initSimData( simdat, *edata, 0.0 );
376 DbgLv(2) << " CR:112 rss now" << US_Memory::rss_now() << "cc" << cc;
377 if (dbg_level>1 && thrnrank<2 && cc==0) {
378  model.debug(); dset->simparams.debug(); }
379 
380  // Calculate Astfem_RSA solution (Lamm equations)
381  US_Astfem_RSA astfem_rsa( model, dset->simparams );
382 DbgLv(2) << " CR:113 rss now" << US_Memory::rss_now() << "cc" << cc;
383 
384  astfem_rsa.set_debug_flag( dbg_level );
385 
386  astfem_rsa.calculate( simdat );
387 
388 #if 0
389 if (dbg_level>0 && thrnrank==1 && cc==0) {
390  model.debug(); dset->simparams.debug(); }
391 DbgLv(1) << "CR: simdat nsc npt" << simdat.scanCount() << simdat.pointCount();
392 int nsc=simdat.scanCount();
393 int npt=simdat.pointCount();
394 int ms=nsc/2;
395 int ls=nsc-1;
396 int mp=npt/2;
397 int lp=npt-1;
398 DbgLv(1) << "CR: edat:"
399  << edata->value( 0, 0) << edata->value( 0,mp) << edata->value( 0,lp)
400  << edata->value(ms, 0) << edata->value(ms,mp) << edata->value(ms,lp)
401  << edata->value(ls, 0) << edata->value(ls,mp) << edata->value(ls,lp);
402 DbgLv(1) << "CR: sdat:"
403  << simdat.value( 0, 0) << simdat.value( 0,mp) << simdat.value( 0,lp)
404  << simdat.value(ms, 0) << simdat.value(ms,mp) << simdat.value(ms,lp)
405  << simdat.value(ls, 0) << simdat.value(ls,mp) << simdat.value(ls,lp);
406 #endif
407 DbgLv(2) << " CR:114 rss now" << US_Memory::rss_now() << "cc" << cc;
408  if ( abort ) return;
409 
410  if ( banddthr )
411  { // If band forming, hold data within thresholds; skip if all-zero
412  if ( data_threshold( &simdat, zerothr, linethr, maxod, mfactor) )
413  continue;
414 
415  ksols++;
416  }
417 
418  simulations << simdat; // Save simulation (each dataset,solute)
419 DbgLv(2) << " CR:115 rss now" << US_Memory::rss_now() << "cc" << cc;
420 
421  // Populate the A matrix for the NNLS routine with simulation
422 DbgLv(1) << " CR: A fill kodl" << kodl << "bndthr ksol" << banddthr << ksols;
423 int ks=kk;
424  if ( kodl == 0 )
425  { // Normal case of no ODlimit substitutions
426  for ( int ss = 0; ss < nscans; ss++ )
427  for ( int rr = 0; rr < npoints; rr++ )
428  nnls_a[ kk++ ] = simdat.value( ss, rr );
429 DbgLv(2) << " CR:116 rss now" << US_Memory::rss_now() << "cc" << cc;
430 //if(lim_offs>1&&(thrnrank==1||thrnrank==11))
431 DbgLv(1) << "CR: ks kk" << ks << kk
432  << "nnA s...k" << nnls_a[ks] << nnls_a[ks+1] << nnls_a[kk-2] << nnls_a[kk-1]
433  << "cc ee" << cc << ee;
434  }
435  else
436  { // Special case where ODlimit substitutions are in B matrix
437  for ( int ss = 0; ss < nscans; ss++ )
438  {
439  for ( int rr = 0; rr < npoints; rr++ )
440  { // Fill A with simulations (or zero where B has zero)
441  if ( nnls_b[ bx++ ] != 0.0 )
442  nnls_a[ kk++ ] = simdat.value( ss, rr );
443  else
444  nnls_a[ kk++ ] = 0.0;
445  }
446  }
447  }
448 
449  if ( tikreg )
450  { // For Tikhonov Regularization append to each column
451  for ( int aa = 0; aa < nsolutes; aa++ )
452  {
453  nnls_a[ kk++ ] = ( aa == cc ) ? alphad : 0.0;
454  }
455  }
456 
457 DbgLv(1) << "CR: NNLS A filled ee" << ee << lim_offs;
458 DbgLv(1) << "CR: NNLS &simdat" << &simdat;
459 DbgLv(1) << "CR: NNLS &model " << &model;
460 DbgLv(1) << "CR: NNLS &nnls_a" << &nnls_a;
461 DbgLv(1) << "CR: NNLS &simulations" << &simulations;
462 DbgLv(1) << "CR: NNLS astfem_rsa" << &astfem_rsa;
463  } // Each data set
464 DbgLv(1) << "CR: NNLS A filled lo" << lim_offs;
465 
466  if ( signal_wanted && ++kstep == increp )
467  { // If asked for and step at increment, signal progress
468  emit work_progress( increp );
469  kstep = 0; // Reset step count
470  }
471 DbgLv(1) << "CR: NNLS A filled EMIT complete";
472 
473 #if 0
474  if ( kk >= iasize && iasize < navals )
475  {
476 int npav = US_Memory::memory_profile();
477 DbgLv(0) << " CR:na iasize navals" << iasize << navals << "kk narows npav" << kk << narows << npav;
478  iasize = qMin( navals, ( iasize + narows * 4 ) );
479  nnls_a.resize( iasize );
480 DbgLv(0) << " CR:na (3)size" << nnls_a.size() << "npav" << US_Memory::memory_profile();
481  }
482 #endif
483 DbgLv(1) << "CR: NNLS A filled cc" << cc << nsolutes;
484  } // Each solute
485 DbgLv(1) << "CR: NNLS A filled";
486  } // Constant vbar
487 
488  else if ( stype == 1 || stype > 9 )
489  { // Special case of varying vbar with constant f/f0 (or other)
490  int attr_x = 0; // Default X is s
491  int attr_y = 3; // Default Y is vbar
492  int attr_z = 1; // Default fixed is f/f0
493  int smask = ( attr_x << 6 ) | ( attr_y << 3 ) | attr_z;
494 
495  if ( stype > 9 )
496  { // Explicitly given attribute types
497  attr_x = ( stype >> 6 ) & 7;
498  attr_y = ( stype >> 3 ) & 7;
499  attr_z = stype & 7;
500  smask = stype;
501  }
502 DbgLv(1) << "CR: attr_ x,y,z" << attr_x << attr_y << attr_z << stype << smask
503  << "use_zsol" << use_zsol;
504 
505  if ( ! use_zsol )
506  {
507  zcomponent.vbar20 = data_sets[ 0 ]->vbar20;
508  set_comp_attr( zcomponent, sim_vals.solutes[ 0 ], attr_z );
509  }
510 
511  for ( int cc = 0; cc < nsolutes; cc++ )
512  { // Solve for each solute
513  if ( abort ) return;
514  int bx = 0;
515 
516  for ( int ee = offset; ee < lim_offs; ee++ )
517  { // Solve for each data set
518  US_Math2::SolutionData sd;
519  DataSet* dset = data_sets[ ee ];
520  US_DataIO::EditedData* edata = banddthr ? &wdata : &dset->run_data;
521  US_DataIO::RawData simdat;
522  int npoints = edata->pointCount();
523  int nscans = edata->scanCount();
524  model.components[ 0 ] = zcomponent;
525 
526  // Set model with standard space s,k,v (or other 3 attributes)
527  if ( use_zsol )
528  {
529  US_ZSolute::set_mcomp_values( model.components[ 0 ],
530  sim_vals.zsolutes[ cc ], smask );
531 //DbgLv(1) << "CR: cc" << cc << "model s,k,w,v,d,c"
532 // << model.components[0].s << model.components[0].f_f0
533 // << model.components[0].mw << model.components[0].vbar20
534 // << model.components[0].D << model.components[0].signal_concentration;
535  }
536  else
537  {
538  set_comp_attr( model.components[ 0 ], sim_vals.solutes[ cc ],
539  attr_x );
540  set_comp_attr( model.components[ 0 ], sim_vals.solutes[ cc ],
541  attr_y );
542  }
543 
544  // Fill in the missing component values
545  model.update_coefficients();
546 
547  // Convert to experimental space
548  double avtemp = dset->temperature;
549  sd.viscosity = dset->viscosity;
550  sd.density = dset->density;
551  sd.manual = dset->manual;
552  sd.vbar20 = model.components[ 0 ].vbar20;
553  sd.vbar = US_Math2::adjust_vbar20( sd.vbar20, avtemp );
554  US_Math2::data_correction( avtemp, sd );
555 
556  model.components[ 0 ].s /= sd.s20w_correction;
557  model.components[ 0 ].D /= sd.D20w_correction;
558 //DbgLv(1) << "CR: exp.space model s,k,w,v,d,c"
559 // << model.components[0].s << model.components[0].f_f0
560 // << model.components[0].mw << model.components[0].vbar20
561 // << model.components[0].D << model.components[0].signal_concentration;
562 
563  // Initialize simulation data with the experiment's grid
564  US_AstfemMath::initSimData( simdat, *edata, 0.0 );
565 if (dbg_level>1 && thrnrank==1 && cc==0) {
566  model.debug(); dset->simparams.debug(); }
567 DbgLv(1) << "CR: simdat nsc npt" << simdat.scanCount() << simdat.pointCount();
568 
569  // Calculate Astfem_RSA solution (Lamm equations)
570  US_Astfem_RSA astfem_rsa( model, dset->simparams );
571 
572  astfem_rsa.set_debug_flag( dbg_level );
573 
574  astfem_rsa.calculate( simdat );
575 #if 0
576 int nsc=simdat.scanCount();
577 int npt=simdat.pointCount();
578 int ms=nsc/2;
579 int ls=nsc-1;
580 int mp=npt/2;
581 int lp=npt-1;
582 DbgLv(1) << "CR: edat:"
583  << edata->value( 0, 0) << edata->value( 0,mp) << edata->value( 0,lp)
584  << edata->value(ms, 0) << edata->value(ms,mp) << edata->value(ms,lp)
585  << edata->value(ls, 0) << edata->value(ls,mp) << edata->value(ls,lp);
586 DbgLv(1) << "CR: sdat:"
587  << simdat.value( 0, 0) << simdat.value( 0,mp) << simdat.value( 0,lp)
588  << simdat.value(ms, 0) << simdat.value(ms,mp) << simdat.value(ms,lp)
589  << simdat.value(ls, 0) << simdat.value(ls,mp) << simdat.value(ls,lp);
590 #endif
591  if ( abort ) return;
592 
593  if ( banddthr )
594  { // If band forming, hold data within thresholds; skip if all-zero
595  if ( data_threshold( &simdat, zerothr, linethr, maxod, mfactor ) )
596  continue;
597 
598  ksols++;
599  }
600 
601  simulations << simdat; // Save simulation (each datset,solute)
602 
603  // Populate the A matrix for the NNLS routine with simulation
604 DbgLv(1) << " CR: A fill kodl" << kodl << "bndthr ksol" << banddthr << ksols;
605 //int ks=kk;
606  if ( kodl == 0 )
607  { // Normal case of no ODlimit substitutions
608  for ( int ss = 0; ss < nscans; ss++ )
609  for ( int rr = 0; rr < npoints; rr++ )
610  nnls_a[ kk++ ] = simdat.value( ss, rr );
611  }
612  else
613  { // Special case where ODlimit substitutions are in B matrix
614  for ( int ss = 0; ss < nscans; ss++ )
615  {
616  for ( int rr = 0; rr < npoints; rr++ )
617  { // Fill A with simulations (or zero where B has zero)
618  if ( nnls_b[ bx++ ] != 0.0 )
619  nnls_a[ kk++ ] = simdat.value( ss, rr );
620  else
621  nnls_a[ kk++ ] = 0.0;
622  }
623  }
624  }
625 //DbgLv(1) << "CR: ks kk" << ks << kk
626 // << "nnA s...k" << nnls_a[ks] << nnls_a[ks+1] << nnls_a[kk-2] << nnls_a[kk-1]
627 // << "cc ee" << cc << ee << "kodl" << kodl;
628 
629  if ( tikreg )
630  { // For Tikhonov Regularization append to each column
631  for ( int aa = 0; aa < nsolutes; aa++ )
632  {
633  nnls_a[ kk++ ] = ( aa == cc ) ? alphad : 0.0;
634  }
635  }
636 
637  } // Each data set
638 
639  if ( signal_wanted && ++kstep == increp )
640  { // If asked for and step at increment, signal progress
641  emit work_progress( increp );
642  kstep = 0; // Reset step count
643  }
644  } // Each solute
645  } // Constant f/f0
646 
647  else
648  { // Special case of custom grid
649  int attr_x = 0; // Set X is s
650  int attr_y = 4; // Set Y is D
651  int attr_z = 3; // Set Z is vbar
652  int smask = ( attr_x << 6 ) | ( attr_y << 3 ) | attr_z;
653 
654  for ( int cc = 0; cc < nsolutes; cc++ )
655  { // Solve for each solute
656  if ( abort ) return;
657  int bx = 0;
658 
659  for ( int ee = offset; ee < lim_offs; ee++ )
660  { // Solve for each data set
661  US_Math2::SolutionData sd;
662  DataSet* dset = data_sets[ ee ];
663  US_DataIO::EditedData* edata = banddthr ? &wdata : &dset->run_data;
664  US_DataIO::RawData simdat;
665  int npoints = edata->pointCount();
666  int nscans = edata->scanCount();
667  double avtemp = dset->temperature;
668  model.components[ 0 ] = zcomponent;
669 
670  if ( use_zsol )
671  {
672  US_ZSolute::set_mcomp_values( model.components[ 0 ],
673  sim_vals.zsolutes[ cc ], smask );
674  }
675  else
676  {
677  set_comp_attr( model.components[ 0 ], sim_vals.solutes[ cc ],
678  attr_x );
679  set_comp_attr( model.components[ 0 ], sim_vals.solutes[ cc ],
680  attr_y );
681  set_comp_attr( model.components[ 0 ], sim_vals.solutes[ cc ],
682  attr_z );
683  }
684 
685  // Fill in the missing component values
686  model.update_coefficients();
687 
688  // Convert to experimental space
689  sd.viscosity = dset->viscosity;
690  sd.density = dset->density;
691  sd.manual = dset->manual;
692  sd.vbar20 = model.components[ 0 ].vbar20;
693  sd.vbar = US_Math2::adjust_vbar20( sd.vbar20, avtemp );
694 
695  US_Math2::data_correction( avtemp, sd );
696 
697  model.components[ 0 ].s /= sd.s20w_correction;
698  model.components[ 0 ].D /= sd.D20w_correction;
699 
700  // Initialize simulation data with the experiment's grid
701  US_AstfemMath::initSimData( simdat, *edata, 0.0 );
702 if (dbg_level>1 && thrnrank==1 && cc==0) {
703  model.debug(); dset->simparams.debug(); }
704 
705  // Calculate Astfem_RSA solution (Lamm equations)
706  US_Astfem_RSA astfem_rsa( model, dset->simparams );
707 
708  astfem_rsa.set_debug_flag( dbg_level );
709 
710  astfem_rsa.calculate( simdat );
711  if ( abort ) return;
712 
713  if ( banddthr )
714  { // If band forming, hold data within thresholds; skip if all-zero
715  if ( data_threshold( &simdat, zerothr, linethr, maxod, mfactor ) )
716  continue;
717 
718  ksols++;
719  }
720 
721  simulations << simdat; // Save simulation (each datset,solute)
722 
723  // Populate the A matrix for the NNLS routine with simulation
724  if ( kodl == 0 )
725  { // Normal case of no ODlimit substitutions
726  for ( int ss = 0; ss < nscans; ss++ )
727  for ( int rr = 0; rr < npoints; rr++ )
728  nnls_a[ kk++ ] = simdat.value( ss, rr );
729  }
730  else
731  { // Special case where ODlimit substitutions are in B matrix
732  for ( int ss = 0; ss < nscans; ss++ )
733  {
734  for ( int rr = 0; rr < npoints; rr++ )
735  { // Fill A with simulations (or zero where B has zero)
736  if ( nnls_b[ bx++ ] != 0.0 )
737  nnls_a[ kk++ ] = simdat.value( ss, rr );
738  else
739  nnls_a[ kk++ ] = 0.0;
740  }
741  }
742  }
743 
744  if ( tikreg )
745  { // For Tikhonov Regularization append to each column
746  for ( int aa = 0; aa < nsolutes; aa++ )
747  {
748  nnls_a[ kk++ ] = ( aa == cc ) ? alphad : 0.0;
749  }
750  }
751  } // Each data set
752 
753  if ( signal_wanted && ++kstep == increp )
754  { // If asked for and step at increment, signal progress
755  emit work_progress( increp );
756  kstep = 0; // Reset step count
757  }
758  } // Each solute
759  }
760 
761  nsolutes = banddthr ? ksols : nsolutes;
762 
763  if ( signal_wanted && kstep > 0 ) // If signals and steps done, report
764  emit work_progress( kstep );
765 
766  if ( abort ) return;
767 
768  int kstodo = nsolutes / 50; // Set steps count for NNLS
769  kstodo = max( kstodo, 2 );
770 DbgLv(2) << " CR:200 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
771 
772 DebugTime("BEG:clcr-nn");
773  if ( calc_ti )
774  { // Compute TI Noise (and, optionally, RI Noise)
775  if ( abort ) return;
776  // Compute a_tilde, the average experiment signal at each time
777  QVector< double > a_tilde( nrinois, 0.0 );
778 
779  if ( calc_ri )
780  compute_a_tilde( a_tilde, nnls_b );
781 
782  // Compute a_bar, the average experiment signal at each radius
783  QVector< double > a_bar( ntinois, 0.0 );
784  compute_a_bar( a_bar, a_tilde, nnls_b );
785 
786  // Compute L_tildes, the average signal at each radius (if RI noise)
787  QVector< double > L_tildes( nrinois * nsolutes, 0.0 );
788 
789  if ( calc_ri )
790  compute_L_tildes( nrinois, nsolutes, L_tildes, nnls_a );
791 
792  // Compute L_bars
793  QVector< double > L_bars( ntinois * nsolutes, 0.0 );
794  compute_L_bars( nsolutes, nrinois, ntinois, ntotal,
795  L_bars, nnls_a, L_tildes );
796 
797  // Set up small_a, small_b for alternate nnls
798 DbgLv(1) << " set SMALL_A+B";
799  QVector< double > small_a( sq( nsolutes ), 0.0 );
800  QVector< double > small_b( nsolutes, 0.0 );
801 
802  ti_small_a_and_b( nsolutes, ntotal, ntinois,
803  small_a, small_b, a_bar, L_bars, nnls_a, nnls_b );
804  if ( abort ) return;
805 
806  // Do NNLS to compute concentrations (nnls_x)
807 DbgLv(1) << " noise small NNLS";
808  US_Math2::nnls( small_a.data(), nsolutes, nsolutes, nsolutes,
809  small_b.data(), nnls_x.data() );
810  if ( abort ) return;
811 
812  // This is Sum( concentration * Lamm ) for the models after NNLS
813  QVector< double > L( ntotal, 0.0 );
814  compute_L( ntotal, nsolutes, L, nnls_a, nnls_x );
815 
816  // Now L contains the best fit sum of L equations
817  // Compute L_tilde, the average model signal at each radius
818  QVector< double > L_tilde( nrinois, 0.0 );
819 
820  if ( calc_ri )
821  compute_L_tilde( L_tilde, L );
822 
823  // Compute L_bar, the average model signal at each radius
824  QVector< double > L_bar( ntinois, 0.0 );
825  compute_L_bar( L_bar, L, L_tilde );
826 
827  // Compute ti noise
828  for ( int ii = 0; ii < ntinois; ii++ )
829  tinvec[ ii ] = a_bar[ ii ] - L_bar[ ii ];
830 
831  if ( calc_ri )
832  { // Compute RI_noise
833  for ( int ii = 0; ii < nrinois; ii++ )
834  rinvec[ ii ] = a_tilde[ ii ] - L_tilde[ ii ];
835  }
836 
837  if ( signal_wanted )
838  emit work_progress( kstodo ); // Report noise NNLS steps done
839  } // End tinoise and optional rinoise calculation
840 
841  else if ( calc_ri )
842  { // Compute RI noise (when RI only)
843  if ( abort ) return;
844  // Compute a_tilde, the average experiment signal at each time
845  QVector< double > a_tilde( nrinois, 0.0 );
846  compute_a_tilde( a_tilde, nnls_b );
847 
848  // Compute L_tildes, the average signal at each radius
849  QVector< double > L_tildes( nrinois * nsolutes, 0.0 );
850  compute_L_tildes( nrinois, nsolutes, L_tildes, nnls_a );
851 
852  // Set up small_a, small_b for the nnls
853  QVector< double > small_a( sq( nsolutes ), 0.0 );
854  QVector< double > small_b( nsolutes, 0.0 );
855  if ( abort ) return;
856 
857  ri_small_a_and_b( nsolutes, ntotal, nrinois,
858  small_a, small_b, a_tilde, L_tildes, nnls_a, nnls_b );
859  if ( abort ) return;
860 
861  US_Math2::nnls( small_a.data(), nsolutes, nsolutes, nsolutes,
862  small_b.data(), nnls_x.data() );
863  if ( abort ) return;
864 
865  // This is sum( concentration * Lamm ) for the models after NNLS
866  QVector< double > L( ntotal, 0.0 );
867  compute_L( ntotal, nsolutes, L, nnls_a, nnls_x );
868 
869  // Now L contains the best fit sum of L equations
870  // Compute L_tilde, the average model signal at each radius
871  QVector< double > L_tilde( nrinois, 0.0 );
872  compute_L_tilde( L_tilde, L );
873 
874  // Compute ri_noise (Is this correct????)
875  for ( int ii = 0; ii < nrinois; ii++ )
876  rinvec[ ii ] = a_tilde[ ii ] - L_tilde[ ii ];
877 
878  if ( signal_wanted )
879  emit work_progress( kstodo ); // Report noise NNLS steps done
880  } // End rinoise alone calculation
881 
882  else
883  { // No TI or RI noise
884  if ( abort ) return;
885 DbgLv(2) << " CR:210 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
886 
887  if ( ASave != NULL && BSave != NULL )
888  {
889  sv_nnls_a = nnls_a;
890  sv_nnls_b = nnls_b;
891 DbgLv(1) << "CR: sv_nnls_a size" << sv_nnls_a.size() << nnls_a.size();
892  }
893 
894  US_Math2::nnls( nnls_a.data(), narows, narows, nsolutes,
895  nnls_b.data(), nnls_x.data() );
896 DbgLv(2) << " CR:211 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
897 if(lim_offs>1&&(thrnrank==1||thrnrank==11)) DbgLv(1) << "CR: narows nsolutes" << narows << nsolutes;
898 if(lim_offs>1&&(thrnrank==1||thrnrank==11)) DbgLv(1) << "CR: a0 a1 b0 b1"
899  << nnls_a[0] << nnls_a[1] << nnls_b[0] << nnls_b[1];
900 
901  if ( abort ) return;
902 
903  if ( signal_wanted )
904  emit work_progress( kstodo ); // Report no-noise NNLS steps done
905  // Note: ti_noise and ri_noise are already zero
906 
907  } // End of core calculations
908 
909  sim_vals.maxrss = US_Memory::rss_max( sim_vals.maxrss );
910 //DbgLv(1) << " CR:na rss now,max" << US_Memory::rss_now() << sim_vals.maxrss;
911 DbgLv(1) << " CR:na rss now,max" << US_Memory::rss_now() << sim_vals.maxrss
912  << &sim_vals;
913 
914  nnls_a.clear();
915  nnls_b.clear();
916 DebugTime("END:clcr-nn");
917 
918 DbgLv(1) << "CR: kstodo" << kstodo;
919  if ( abort ) return;
920 
921  // Size simulation data and initialize concentrations to zero
922  int kscans = 0;
923  int jscan = 0;
924 
925  for ( int ee = offset; ee < lim_offs; ee++ )
926  {
927  US_DataIO::RawData tdata; // Init temp sim data with edata's grid
928  US_AstfemMath::initSimData( tdata, data_sets[ ee ]->run_data, 0.0 );
929 
930  int nscans = tdata.scanData.size();
931  kscans += nscans;
932  sim_vals.sim_data.scanData.resize( kscans );
933 
934  if ( ee == offset )
935  {
936  sim_vals.sim_data = tdata; // Init (first) dataset sim data
937  jscan += nscans;
938  }
939  else
940  {
941  for ( int ss = 0; ss < nscans; ss++ )
942  { // Append zeroed-scans sim_data for multiple data sets
943 DbgLv(1) << "CR: ss jscan" << ss << jscan << "ee kscans nscans" << ee << kscans << nscans;
944  sim_vals.sim_data.scanData[ jscan++ ] = tdata.scanData[ ss ];
945  }
946  }
947 DbgLv(2) << " CR:221 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
948  }
949 if(lim_offs>1&&(thrnrank==1||thrnrank==11))
950 DbgLv(1) << "CR: jscan kscans" << jscan << kscans;
951 
952  // This is a little tricky. The simulations structure was created above
953  // in a loop that alternates experiments with each solute. In the loop
954  // below, the simulation structure alters that order to group all solutes
955  // for each experiment together
956 
957 if(thrnrank==1) DbgLv(1) << "CR: nsolutes" << nsolutes;
958 if(lim_offs>1&&(thrnrank==1||thrnrank==11)) DbgLv(1) << "CR: nsolutes" << nsolutes;
959  for ( int cc = 0; cc < nsolutes; cc++ )
960  {
961  double soluval = nnls_x[ cc ]; // Computed concentration, this solute
962 if(thrnrank==1) DbgLv(1) << "CR: cc soluval" << cc << soluval;
963 if(lim_offs>1&&(thrnrank==1||thrnrank==11)) DbgLv(1) << "CR: cc soluval" << cc << soluval;
964 
965  if ( soluval > 0.0 )
966  { // If concentration non-zero, need to sum in simulation data
967 if(lim_offs>1&&(thrnrank==1||thrnrank==11))
968 DbgLv(1) << "CR: cc soluval" << cc << soluval;
969  int scnx = 0;
970 
971  for ( int ee = offset; ee < lim_offs; ee++ )
972  {
973  // Input sims (ea.dset, ea.solute); out sims (sum.solute, ea.dset)
974  int sim_ix = cc * dataset_count + ee - offset;
975  US_DataIO::RawData* idata = &simulations[ sim_ix ];
976  US_DataIO::RawData* sdata = &sim_vals.sim_data;
977  int npoints = idata->pointCount();
978  int nscans = idata->scanCount();
979 if(lim_offs>1&&(thrnrank==1||thrnrank==11))
980 DbgLv(1) << "CR: ee sim_ix np ns scnx" << ee << sim_ix << npoints << nscans << scnx;
981 
982  for ( int ss = 0; ss < nscans; ss++, scnx++ )
983  {
984  for ( int rr = 0; rr < npoints; rr++ )
985  {
986  sdata->scanData[ scnx ].rvalues[ rr ] +=
987  ( soluval * idata->value( ss, rr ) );
988  }
989  }
990 DbgLv(2) << " CR:231 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
991 //*DEBUG*
992 int ss=nscans/2;
993 int rr=npoints/2;
994 if( thrnrank==1 ) {
995  DbgLv(1) << "CR: scnx ss rr" << scnx << ss << rr;
996  if(use_zsol) {
997  DbgLv(1) << "CR: x y z" << sim_vals.zsolutes[cc].x*1.0e+13
998  << sim_vals.zsolutes[cc].y << sim_vals.zsolutes[cc].z << "sval" << soluval
999  << "idat sdat" << idata->value(ss,rr) << sdata->value(ss,rr);
1000  if (soluval>100.0) {
1001  double drval=0.0; double dmax=0.0; double dsum=0.0;
1002  for ( int ss=0;ss<nscans;ss++ ) { for ( int rr=0; rr<npoints; rr++ ) {
1003  drval=idata->value(ss,rr); dmax=qMax(dmax,drval); dsum+=drval; }}
1004  DbgLv(1) << "CR:B x y" << sim_vals.zsolutes[cc].x*1.0e+13
1005  << sim_vals.zsolutes[cc].y << "sval" << soluval << "amax asum" << dmax << dsum; }
1006  } else {
1007  DbgLv(1) << "CR: s k v" << sim_vals.solutes[cc].s*1.0e+13
1008  << sim_vals.solutes[cc].k << sim_vals.solutes[cc].v << "sval" << soluval
1009  << "idat sdat" << idata->value(ss,rr) << sdata->value(ss,rr);
1010  if (soluval>100.0) {
1011  double drval=0.0; double dmax=0.0; double dsum=0.0;
1012  for ( int ss=0;ss<nscans;ss++ ) { for ( int rr=0; rr<npoints; rr++ ) {
1013  drval=idata->value(ss,rr); dmax=qMax(dmax,drval); dsum+=drval; }}
1014  DbgLv(1) << "CR:B s k" << sim_vals.solutes[cc].s*1.0e+13
1015  << sim_vals.solutes[cc].k << "sval" << soluval << "amax asum" << dmax << dsum; }
1016  }
1017 }
1018 //*DEBUG*
1019  }
1020  }
1021  }
1022 
1023  double rmsds[ dataset_count ];
1024  int kntva[ dataset_count ];
1025  double variance = 0.0;
1026  int tinoffs = 0;
1027  int rinoffs = 0;
1028  int ktotal = 0;
1029  int scnx = 0;
1030  int kdsx = 0;
1031  US_DataIO::RawData* sdata = &sim_vals.sim_data;
1032  US_DataIO::RawData* resid = &sim_vals.residuals;
1033  sim_vals.residuals = *sdata;
1034 DbgLv(2) << " CR:301 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
1035 
1036  // Calculate residuals and rmsd values
1037  for ( int ee = offset; ee < lim_offs; ee++, kdsx++ )
1038  {
1039  DataSet* dset = data_sets[ ee ];
1040  US_DataIO::EditedData* edata = banddthr ? &wdata : &dset->run_data;
1041 // US_AstfemMath::initSimData( sim_vals.residuals, *edata, 0.0 );
1042  int npoints = edata->pointCount();
1043  int nscans = edata->scanCount();
1044  int kntcs = 0;
1045  double varidset = 0.0;
1046 
1047  for ( int ss = 0; ss < nscans; ss++, scnx++ )
1048  { // Create residuals dset: exp - sim - noise(s)
1049 
1050  for ( int rr = 0; rr < npoints; rr++ )
1051  {
1052 
1053  double resval = edata->value( ss, rr )
1054  - sdata->value( scnx, rr )
1055  - tinvec[ rr + tinoffs ]
1056  - rinvec[ ss + rinoffs ];
1057  resid->setValue( scnx, rr, resval );
1058 
1059  double r2 = sq( resval );
1060  variance += r2;
1061  varidset += r2;
1062  kntcs++;
1063  }
1064  }
1065 
1066  rmsds[ kdsx ] = varidset; // Variance for a data set
1067  kntva[ kdsx ] = kntcs;
1068  ktotal += kntcs;
1069  tinoffs += npoints;
1070  rinoffs += nscans;
1071 DbgLv(1) << "CR: kdsx variance" << kdsx << varidset << variance;
1072  }
1073 
1074  sim_vals.variances.resize( dataset_count );
1075  variance /= (double)ktotal; // Total sets variance
1076  sim_vals.variance = variance;
1077  kk = 0;
1078 
1079  for ( int ee = offset; ee < lim_offs; ee++ )
1080  { // Scale variances for each data set
1081  sim_vals.variances[ kk ] = rmsds[ kk ] / (double)( kntva[ kk ] );
1082 DbgLv(1) << "CR: kk variance" << sim_vals.variances[kk];
1083  kk++;
1084  }
1085 
1086  // Store solutes for return
1087  kk = 0;
1088 
1089  if ( use_zsol )
1090  { // Use xyZ type solutes
1091  for ( int cc = 0; cc < nsolutes; cc++ )
1092  {
1093  if ( nnls_x[ cc ] > 0.0 )
1094  { // Store solutes with non-zero concentrations
1095  sim_vals.zsolutes[ cc ].c = nnls_x[ cc ];
1096  sim_vals.zsolutes[ kk++ ] = sim_vals.zsolutes[ cc ];
1097  }
1098  }
1099  }
1100  else
1101  { // Use old type solutes
1102  for ( int cc = 0; cc < nsolutes; cc++ )
1103  {
1104  if ( nnls_x[ cc ] > 0.0 )
1105  { // Store solutes with non-zero concentrations
1106  sim_vals.solutes[ cc ].c = nnls_x[ cc ];
1107  sim_vals.solutes[ kk++ ] = sim_vals.solutes[ cc ];
1108  }
1109  }
1110  }
1111 DbgLv(2) << " CR:310 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
1112 
1113  // Truncate solutes at non-zero count
1114 DbgLv(1) << "CR: out solutes size" << kk;
1115  if ( use_zsol )
1116  {
1117  sim_vals.zsolutes.resize( qMax( kk, 1 ) );
1118 DbgLv(1) << "CR: jj solute-c" << 0 << sim_vals.zsolutes[0].c;
1119 DbgLv(1) << "CR: jj solute-c" << 1 << (kk>1?sim_vals.zsolutes[1].c:0.0);
1120 DbgLv(1) << "CR: jj solute-c" << kk-2 << (kk>1?sim_vals.zsolutes[kk-2].c:0.0);
1121 DbgLv(1) << "CR: jj solute-c" << kk-1 << (kk>0?sim_vals.zsolutes[kk-1].c:0.0);
1122  }
1123  else
1124  {
1125  sim_vals.solutes.resize( qMax( kk, 1 ) );
1126 DbgLv(1) << "CR: jj solute-c" << 0 << sim_vals.solutes[0].c;
1127 DbgLv(1) << "CR: jj solute-c" << 1 << (kk>1?sim_vals.solutes[1].c:0.0);
1128 DbgLv(1) << "CR: jj solute-c" << kk-2 << (kk>1?sim_vals.solutes[kk-2].c:0.0);
1129 DbgLv(1) << "CR: jj solute-c" << kk-1 << (kk>0?sim_vals.solutes[kk-1].c:0.0);
1130  }
1131  if ( abort ) return;
1132 
1133  // Fill noise objects with any calculated vectors
1134  if ( calc_ti )
1135  sim_vals.ti_noise << tinvec;
1136 
1137  if ( calc_ri )
1138  sim_vals.ri_noise << rinvec;
1139 
1140  // Compute and return the xnorm-squared value of concentrations
1141  nsolutes = ( use_zsol ) ? sim_vals.zsolutes.size()
1142  : sim_vals.solutes.size();
1143 DbgLv(1) << "CR: nsolutes" << nsolutes;
1144  double xnorm = 0.0;
1145  for ( int jj = 0; jj < nsolutes; jj++ )
1146  {
1147  double cval = use_zsol ? sim_vals.zsolutes[ jj ].c
1148  : sim_vals.solutes [ jj ].c;
1149  xnorm += sq( cval );
1150  }
1151 DbgLv(2) << " CR:320 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
1152 
1153  sim_vals.xnormsq = xnorm;
1154 DbgLv(1) << "CR: xnormsq" << xnorm;
1155 
1156  // If specified, return the A and B matrices (with or without padding)
1157  if ( ASave != NULL && BSave != NULL )
1158  {
1159  if ( padAB )
1160  { // Return A and B with padding for regularization
1161  if ( tikreg )
1162  { // If regularization was used, return A and B with padding intact
1163  (*ASave) = sv_nnls_a;
1164  (*BSave) = sv_nnls_b;
1165  }
1166 
1167  else
1168  { // If no regularization, return A and B with padding added
1169  ASave->clear();
1170  BSave->clear();
1171  int kk = 0;
1172 
1173  for ( int cc = 0; cc < nsolutes; cc++ )
1174  { // Save and pad A matrix, a column at a time
1175  for ( int jj = 0; jj < ntotal; jj++ )
1176  { // Save an A column
1177  (*ASave) << sv_nnls_a[ kk++ ];
1178  }
1179 
1180  for ( int jj = 0; jj < nsolutes; jj++ )
1181  {
1182  (*ASave) << 0.0; // Pad an A column (zeroes for now)
1183  }
1184  }
1185 
1186  for ( int jj = 0; jj < ntotal; jj++ )
1187  { // Save the B matrix (vector)
1188  (*BSave) << sv_nnls_b[ jj ];
1189  }
1190 
1191  for ( int cc = 0; cc < nsolutes; cc++ )
1192  {
1193  (*BSave) << 0.0; // Pad the B vector with zeroes
1194  }
1195  }
1196 
1197  // If we are doing padding, it is in preparation for regularization.
1198  // So, if any noise was calculated, it should be added to the
1199  // B vector.
1200  if ( calc_ri || calc_ti )
1201  {
1202  int nscans = data_sets[ 0 ]->run_data.scanCount();
1203  int npoints = data_sets[ 0 ]->run_data.pointCount();
1204  int kk = 0;
1205 
1206  for ( int ss = 0; ss < nscans; ss++ )
1207  {
1208  double rinoi = calc_ri ? rinvec[ ss ] : 0.0;
1209 
1210  for ( int rr = 0; rr < npoints; rr++ )
1211  {
1212  double tinoi = calc_ti ? tinvec[ rr ] : 0.0;
1213  (*BSave)[ kk++ ] += ( tinoi + rinoi );
1214  }
1215  }
1216  }
1217 
1218  }
1219 
1220  else
1221  { // Return A and B without any regularization padding
1222  if ( tikreg )
1223  { // If regularization was used, remove padding
1224  ASave->clear();
1225  BSave->clear();
1226  int kk = 0;
1227 
1228  for ( int cc = 0; cc < nsolutes; cc++ )
1229  {
1230  for ( int jj = 0; jj < ntotal; jj++ )
1231  { // Save an A column
1232  (*ASave) << sv_nnls_a[ kk++ ];
1233  }
1234 
1235  kk += nsolutes; // Bump past regularization padding
1236  }
1237 
1238  for ( int jj = 0; jj < ntotal; jj++ )
1239  { // Save the B matrix (vector)
1240  (*BSave) << sv_nnls_b[ jj ];
1241  }
1242  }
1243 
1244  else
1245  { // If no regularization, return A and B as they are
1246  (*ASave) = sv_nnls_a;
1247  (*BSave) = sv_nnls_b;
1248 DbgLv(1) << "CR: ASv: sv_nnls_a size" << sv_nnls_a.size() << ASave->size();
1249  }
1250  }
1251  }
1252 DbgLv(2) << " CR:777 rss now" << US_Memory::rss_now() << "thrn" << thrnrank;
1253 
1254 DebugTime("END:calcres");
1255 }
1256 
1257 // Set abort flag
1258 void US_SolveSim::abort_work()
1259 {
1260  abort = true;
1261 }
1262 
1263 // Compute a_tilde, the average experiment signal at each time
1264 void US_SolveSim::compute_a_tilde( QVector< double >& a_tilde,
1265  const QVector< double >& nnls_b )
1266 {
1267  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1268  int npoints = edata->pointCount();
1269  int nscans = edata->scanCount();
1270  int jb = 0;
1271  double avgscale = 1.0 / (double)npoints;
1272 
1273  for ( int ss = 0; ss < nscans; ss++ )
1274  {
1275  for ( int rr = 0; rr < npoints; rr++ )
1276  a_tilde[ ss ] += nnls_b[ jb++ ];
1277 
1278  a_tilde[ ss ] *= avgscale;
1279  }
1280 }
1281 
1282 // Compute L_tildes, the average signal at each radius
1283 void US_SolveSim::compute_L_tildes( int nrinois,
1284  int nsolutes,
1285  QVector< double >& L_tildes,
1286  const QVector< double >& nnls_a )
1287 {
1288  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1289  int npoints = edata->pointCount();
1290  int nscans = edata->scanCount();
1291  double avgscale = 1.0 / (double)npoints;
1292  int a_index = 0;
1293 
1294  for ( int cc = 0; cc < nsolutes; cc++ )
1295  {
1296  int t_index = cc * nrinois;
1297 
1298  for ( int ss = 0; ss < nscans; ss++ )
1299  {
1300  for ( int rr = 0; rr < npoints; rr++ )
1301  L_tildes[ t_index ] += nnls_a[ a_index++ ];
1302 
1303  L_tildes[ t_index++ ] *= avgscale;
1304  }
1305  }
1306 }
1307 
1308 // Compute L_tilde, the average model signal at each radius
1309 void US_SolveSim::compute_L_tilde( QVector< double >& L_tilde,
1310  const QVector< double >& L )
1311 {
1312  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1313  int npoints = edata->pointCount();
1314  int nscans = edata->scanCount();
1315  double avgscale = 1.0 / (double)npoints;
1316  int index = 0;
1317 
1318  for ( int ss = 0; ss < nscans; ss++ )
1319  {
1320  for ( int rr = 0; rr < npoints; rr++ )
1321  L_tilde[ ss ] += L[ index++ ];
1322 
1323  L_tilde[ ss ] *= avgscale;
1324  }
1325 }
1326 
1327 void US_SolveSim::compute_L( int ntotal,
1328  int nsolutes,
1329  QVector< double >& L,
1330  const QVector< double >& nnls_a,
1331  const QVector< double >& nnls_x )
1332 {
1333  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1334  int npoints = edata->pointCount();
1335  int nscans = edata->scanCount();
1336 
1337  for ( int cc = 0; cc < nsolutes; cc++ )
1338  {
1339  double concentration = nnls_x[ cc ];
1340 
1341  if ( concentration > 0.0 )
1342  {
1343  int r_index = cc * ntotal;
1344  int count = 0;
1345 
1346  for ( int ss = 0; ss < nscans; ss++ )
1347  {
1348  for ( int rr = 0; rr < npoints; rr++ )
1349  {
1350  L[ count++ ] += ( concentration * nnls_a[ r_index++ ] );
1351  }
1352  }
1353  }
1354  }
1355 }
1356 
1357 void US_SolveSim::ri_small_a_and_b( int nsolutes,
1358  int ntotal,
1359  int nrinois,
1360  QVector< double >& small_a,
1361  QVector< double >& small_b,
1362  const QVector< double >& a_tilde,
1363  const QVector< double >& L_tildes,
1364  const QVector< double >& nnls_a,
1365  const QVector< double >& nnls_b )
1366 {
1367 DebugTime("BEG:ri_smab");
1368  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1369  int npoints = edata->pointCount();
1370  int nscans = edata->scanCount();
1371  int kstodo = sq( nsolutes ) / 10; // progress steps to report
1372  int incprg = nsolutes / 20; // increment between reports
1373  incprg = max( incprg, 1 );
1374  incprg = min( incprg, 10 );
1375  int jstprg = ( kstodo * incprg ) / nsolutes; // steps for each report
1376  int kstep = 0; // progress counter
1377 
1378  for ( int cc = 0; cc < nsolutes; cc++ )
1379  {
1380  int jsa2 = cc * ntotal;
1381  int jst2 = cc * nrinois;
1382  int jjna = jsa2;
1383  int jjnb = 0;
1384  int jjlt = jst2;
1385  double sum_b = small_b[ cc ];
1386 
1387  for ( int ss = 0; ss < nscans; ss++ )
1388  {
1389  // small_b[ cc ] +=
1390  // ( edata->value( ss, rr ) - a_tilde[ ss ] )
1391  // *
1392  // ( nnls_a[ cc * ntotal + ss * npoints + rr ]
1393  // -
1394  // L_tildes[ cc * nrinois + ss ] );
1395  double atil = a_tilde [ ss ];
1396  double Ltil = L_tildes[ jjlt++ ];
1397 
1398  for ( int rr = 0; rr < npoints; rr++ )
1399  {
1400  sum_b += ( ( nnls_b[ jjnb++ ] - atil )
1401  * ( nnls_a[ jjna++ ] - Ltil ) );
1402  }
1403 
1404  }
1405 
1406  small_b[ cc ] = sum_b;
1407 
1408  for ( int kk = 0; kk < nsolutes; kk++ )
1409  {
1410  //small_a[ kk * nsolutes + cc ] +=
1411  // ( nnls_a[ kk * ntotal + ss * npoints + rr ]
1412  // -
1413  // L_tildes[ kk * nrinois + ss ]
1414  // )
1415  // *
1416  // ( nnls_a[ cc * ntotal + ss * npoints + rr ]
1417  // -
1418  // L_tildes[ cc * nrinois + ss ] );
1419  int jjma = kk * nsolutes + cc;
1420  int jja1 = kk * ntotal;
1421  int jja2 = jsa2;
1422  int jjt1 = kk * nrinois;
1423  int jjt2 = jst2;
1424  double sum_a = small_a[ jjma ];
1425 
1426  for ( int ss = 0; ss < nscans; ss++ )
1427  {
1428  double Ltil1 = L_tildes[ jjt1++ ];
1429  double Ltil2 = L_tildes[ jjt2++ ];
1430 
1431  for ( int rr = 0; rr < npoints; rr++ )
1432  {
1433  sum_a += ( ( nnls_a[ jja1++ ] - Ltil1 )
1434  * ( nnls_a[ jja2++ ] - Ltil2 ) );
1435  }
1436  }
1437 
1438  small_a[ jjma ] = sum_a;
1439  }
1440 
1441  if ( signal_wanted && ++kstep == incprg )
1442  {
1443  emit work_progress( jstprg );
1444  kstodo -= jstprg;
1445  kstep = 0;
1446  }
1447 
1448  if ( abort ) return;
1449  }
1450 
1451  if ( signal_wanted && kstodo > 0 )
1452  emit work_progress( kstodo );
1453 DebugTime("END:ri_smab");
1454 }
1455 
1456 void US_SolveSim::ti_small_a_and_b( int nsolutes,
1457  int ntotal,
1458  int ntinois,
1459  QVector< double >& small_a,
1460  QVector< double >& small_b,
1461  const QVector< double >& a_bar,
1462  const QVector< double >& L_bars,
1463  const QVector< double >& nnls_a,
1464  const QVector< double >& nnls_b )
1465 {
1466 DebugTime("BEG:ti-smab");
1467  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1468  int npoints = edata->pointCount();
1469  int nscans = edata->scanCount();
1470  int kstodo = sq( nsolutes ) / 10; // progress steps to report
1471  int incprg = nsolutes / 20; // increment between reports
1472  incprg = max( incprg, 1 );
1473  incprg = min( incprg, 10 );
1474  int jstprg = ( kstodo * incprg ) / nsolutes; // steps for each report
1475  int kstep = 0; // progress counter
1476 
1477  for ( int cc = 0; cc < nsolutes; cc++ )
1478  {
1479  int jjsa = cc;
1480  int jssa = cc * ntotal;
1481  int jssb = cc * ntinois;
1482  int jjna = jssa;
1483  int jjnb = 0;
1484 
1485  //small_b[ cc ] +=
1486  // ( edata->value( ss, rr ) - a_bar[ rr ] )
1487  // *
1488  // ( nnls_a[ cc * ntotal + ss * npoints + rr ]
1489  // -
1490  // L_bars[ cc * ntinois + rr ] );
1491  double sum_b = small_b[ cc ];
1492 
1493  for ( int ss = 0; ss < nscans; ss++ )
1494  {
1495  int jjlb = jssb;
1496 
1497  for ( int rr = 0; rr < npoints; rr++ )
1498  {
1499  sum_b += ( ( nnls_b[ jjnb++ ] - a_bar [ rr ] )
1500  * ( nnls_a[ jjna++ ] - L_bars[ jjlb++ ] ) );
1501  }
1502  }
1503 
1504  small_b[ cc ] = sum_b;
1505 
1506  //small_a[ kk * nsolutes + cc ] +=
1507  // ( nnls_a[ kk * ntotal + ss * npoints + rr ]
1508  // -
1509  // L_bars[ kk * ntinois + rr ] )
1510  // *
1511  // ( nnls_a[ cc * ntotal + ss * npoints + rr ]
1512  // -
1513  // L_bars[ cc * ntinois + rr ] );
1514  for ( int kk = 0; kk < nsolutes; kk++ )
1515  {
1516  int jjna1 = kk * ntotal;
1517  int jjna2 = jssa;
1518  int jslb1 = kk * ntinois;
1519  int jslb2 = jssb;
1520  double sum_a = small_a[ jjsa ];
1521 
1522  for ( int ss = 0; ss < nscans; ss++ )
1523  {
1524  int jjlb1 = jslb1;
1525  int jjlb2 = jslb2;
1526 
1527  for ( int rr = 0; rr < npoints; rr++ )
1528  {
1529  sum_a += ( ( nnls_a[ jjna1++ ] - L_bars[ jjlb1++ ] )
1530  * ( nnls_a[ jjna2++ ] - L_bars[ jjlb2++ ] ) );
1531  }
1532  }
1533 
1534  small_a[ jjsa ] = sum_a;
1535  jjsa += nsolutes;
1536  }
1537 
1538  if ( signal_wanted && ++kstep == incprg )
1539  {
1540  emit work_progress( jstprg );
1541  kstodo -= jstprg;
1542  kstep = 0;
1543  }
1544 
1545  if ( abort ) return;
1546  }
1547 
1548  if ( signal_wanted && kstodo > 0 )
1549  emit work_progress( kstodo );
1550 DebugTime("END:ti-smab");
1551 }
1552 
1553 void US_SolveSim::compute_L_bar( QVector< double >& L_bar,
1554  const QVector< double >& L,
1555  const QVector< double >& L_tilde )
1556 {
1557  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1558  int npoints = edata->pointCount();
1559  int nscans = edata->scanCount();
1560  double avgscale = 1.0 / (double)nscans;
1561 
1562  for ( int rr = 0; rr < npoints; rr++)
1563  {
1564  // Note L_tilde is always zero when rinoise has not been requested
1565  for ( int ss = 0; ss < nscans; ss++ )
1566  L_bar[ rr ] += ( L[ ss * npoints + rr ] - L_tilde[ ss ] );
1567 
1568  L_bar[ rr ] *= avgscale;
1569  }
1570 }
1571 
1572 // Calculate the average measured concentration at each radius point
1573 void US_SolveSim::compute_a_bar( QVector< double >& a_bar,
1574  const QVector< double >& a_tilde,
1575  const QVector< double >& nnls_b )
1576 {
1577  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1578  int npoints = edata->pointCount();
1579  int nscans = edata->scanCount();
1580  double avgscale = 1.0 / (double)nscans;
1581 
1582  for ( int rr = 0; rr < npoints; rr++ )
1583  {
1584  int jb = rr;
1585 
1586  // Note: a_tilde is always zero when rinoise has not been requested
1587  for ( int ss = 0; ss < nscans; ss++ )
1588  {
1589  a_bar[ rr ] += ( nnls_b[ jb ] - a_tilde[ ss ] );
1590  jb += npoints;
1591  }
1592 
1593  a_bar[ rr ] *= avgscale;
1594  }
1595 }
1596 
1597 // Calculate the average simulated concentration at each radius point
1598 void US_SolveSim::compute_L_bars( int nsolutes,
1599  int nrinois,
1600  int ntinois,
1601  int ntotal,
1602  QVector< double >& L_bars,
1603  const QVector< double >& nnls_a,
1604  const QVector< double >& L_tildes )
1605 {
1606  US_DataIO::EditedData* edata = &data_sets[ d_offs ]->run_data;
1607  int npoints = edata->pointCount();
1608  int nscans = edata->scanCount();
1609  double avgscale = 1.0 / (double)nscans;
1610 
1611  for ( int cc = 0; cc < nsolutes; cc++ )
1612  {
1613  int solute_offset = cc * ntotal;
1614 
1615  for ( int rr = 0; rr < npoints; rr++ )
1616  {
1617  int r_index = cc * ntinois + rr;
1618 
1619  for ( int ss = 0; ss < nscans; ss++ )
1620  {
1621  // Note: L_tildes is always zero when rinoise has not been
1622  // requested
1623 
1624  int n_index = solute_offset + ss * npoints + rr;
1625  int s_index = cc * nrinois + ss;
1626 
1627  L_bars[ r_index ] += ( nnls_a[ n_index ] - L_tildes[ s_index ] );
1628  }
1629 
1630  L_bars[ r_index ] *= avgscale;
1631  }
1632  }
1633 }
1634 
1635 // Debug message with thread/processor number and elapsed time value
1636 void US_SolveSim::DebugTime( QString mtext )
1637 {
1638  if ( dbg_timing )
1639  {
1640  qDebug() << "w" << thrnrank << "TM:" << mtext
1641  << startCalc.msecsTo( QDateTime::currentDateTime() ) / 1000.0;
1642  }
1643 }
1644 
1645 // Modify amplitude of data by thresholds and return flag if all-zero result
1646 bool US_SolveSim::data_threshold( US_DataIO::RawData* sdata,
1647  double zerothr, double linethr, double maxod, double mfactor )
1648 {
1649  int nnzro = 0;
1650  int nzset = 0;
1651  int nntrp = 0;
1652  int nclip = 0;
1653  int npoints = sdata->pointCount();
1654  int nscans = sdata->scanCount();
1655  double clipout = mfactor * maxod;
1656  double thrfact = mfactor / (double)( linethr - zerothr );
1657 double maxs=0.0;
1658 double maxsi=0.0;
1659 
1660  for ( int ss = 0; ss < nscans; ss++ )
1661  {
1662  for ( int rr = 0; rr < npoints; rr++ )
1663  {
1664  double avalue = sdata->value( ss, rr );
1665 maxsi=qMax(maxsi,avalue);
1666 
1667  if ( avalue < zerothr )
1668  { // Less than zero threshold: set to zero
1669  avalue = 0.0;
1670  nzset++;
1671  }
1672 
1673  else if ( avalue < linethr )
1674  { // Between zero and linear threshold: set to interpolated value
1675  avalue *= ( ( avalue - zerothr ) * thrfact );
1676  nntrp++;
1677  }
1678 
1679  else if ( avalue < maxod )
1680  { // Under maximum OD: set to factor times input
1681  avalue *= mfactor;
1682  }
1683 
1684  else
1685  { // Over maximum OD; set to factor times maximum
1686  avalue = clipout;
1687  nclip++;
1688  }
1689 
1690  if ( avalue != 0.0 )
1691  nnzro++;
1692 maxs=qMax(maxs,avalue);
1693 
1694  sdata->setValue( ss, rr, avalue );
1695  }
1696  }
1697 
1698  int lownnz = qRound( minnzsc * (double)( nscans * npoints ) );
1699  nnzro = ( nnzro < lownnz ) ? 0 : nnzro;
1700 DbgLv(1) << " CR:THR: nnzro zs nt cl" << nnzro << nzset << nntrp << nclip;
1701 //if(nnzro>0) {DbgLv(1) << "CR:THR: maxs" << maxs << maxsi << "mfact" << mfactor;}
1702 //else {DbgLv(1) << "CR:THz: maxs" << maxs << nnzro << "mfact" << mfactor;}
1703 
1704  return ( nnzro == 0 );
1705 }
1706 
1707 // Modify amplitude by thresholds and flag if all-zero (for experiment data)
1708 bool US_SolveSim::data_threshold( US_DataIO::EditedData* edata,
1709  double zerothr, double linethr, double maxod, double mfactor )
1710 {
1711  int nnzro = 0;
1712  int npoints = edata->pointCount();
1713  int nscans = edata->scanCount();
1714  double clipout = mfactor * maxod;
1715  double thrfact = mfactor / (double)( linethr - zerothr );
1716 
1717  for ( int ss = 0; ss < nscans; ss++ )
1718  {
1719  for ( int rr = 0; rr < npoints; rr++ )
1720  {
1721  double avalue = edata->value( ss, rr );
1722 
1723  if ( avalue < zerothr )
1724  { // Less than zero threshold: set to zero
1725  avalue = 0.0;
1726  }
1727 
1728  else if ( avalue < linethr )
1729  { // Between zero and linear threshold: set to interpolated value
1730  avalue *= ( ( avalue - zerothr ) * thrfact );
1731  }
1732 
1733  else if ( avalue < maxod )
1734  { // Under maximum OD: set to factor times input
1735  avalue *= mfactor;
1736  }
1737 
1738  else
1739  { // Over maximum OD; set to factor times maximum
1740  avalue = clipout;
1741  }
1742 
1743  if ( avalue != 0.0 )
1744  nnzro++;
1745 
1746  edata->setValue( ss, rr, avalue );
1747  }
1748  }
1749 
1750  return ( nnzro == 0 );
1751 }
1752 
1753 // Set a model component attribute value
1754 void US_SolveSim::set_comp_attr( US_Model::SimulationComponent& component,
1755  US_Solute& solute, int attr_type )
1756 {
1757  switch ( attr_type )
1758  {
1759  default:
1760  case ATTR_S: // Sedimentation Coefficient
1761  component.s = solute.s;
1762  break;
1763  case ATTR_K: // Frictional Ratio
1764  component.f_f0 = solute.k;
1765  break;
1766  case ATTR_W: // Molecular Weight
1767  component.mw = solute.d;
1768  break;
1769  case ATTR_V: // Partial Specific Volume (vbar)
1770  component.vbar20 = solute.v;
1771  break;
1772  case ATTR_D: // Diffusion Coefficient
1773  component.D = solute.d;
1774  break;
1775  case ATTR_F: // Frictional Coefficient
1776  component.f = solute.d;
1777  break;
1778  }
1779 }
1780 
Generated on Sun Jun 21 2015 11:40:43 for UltraScan III by   doxygen 1.8.3.1-20130324