UltraScan III
dmga_master.cpp
Go to the documentation of this file.
1 #include "us_mpi_analysis.h"
2 #include "us_util.h"
3 #include "us_math2.h"
4 #include "us_settings.h"
5 
6 void US_MPI_Analysis::dmga_master( void )
7 {
8  current_dataset = 0;
9  datasets_to_process = data_sets.size();
10  max_depth = 0;
11  calculated_solutes.clear();
12 
13  // Set noise and debug flags
14  simulation_values.noisflag = 0;
15  simulation_values.dbg_level = dbg_level;
16  simulation_values.dbg_timing = dbg_timing;
17 DbgLv(0) << "DEBUG_LEVEL" << simulation_values.dbg_level;
18 
19  if ( data_sets.size() == 1 )
20  {
21  US_AstfemMath::initSimData( simulation_values.sim_data,
22  data_sets[ 0 ]->run_data, 0.0 );
23  US_AstfemMath::initSimData( simulation_values.residuals,
24  data_sets[ 0 ]->run_data, 0.0 );
25  }
26  else
27  {
28  int ntscan = data_sets[ 0 ]->run_data.scanCount();
29  for ( int ii = 1; ii < data_sets.size(); ii++ )
30  ntscan += data_sets[ ii ]->run_data.scanCount();
31  simulation_values.sim_data .scanData.resize( ntscan );
32  simulation_values.residuals.scanData.resize( ntscan );
33  }
34 
35  // Initialize best fitness
36  best_dgenes .reserve( gcores_count );
37  best_fitness.reserve( gcores_count );
38 
39  // Read in the constraints model and build constraints
40  QString cmfname = "../" + parameters[ "DC_model" ];
41  wmodel.load( cmfname ); // Load the constraints model
42  constraints.load_constraints( &wmodel ); // Build the constraints object
43  constraints.get_work_model ( &wmodel ); // Get the base work model
44 DbgLv(1) << "dmga_master: cmfname" << cmfname;
45 DbgLv(1) << "dmga_master: wmodel #comps" << wmodel.components.size()
46  << "#assoc" << wmodel.associations.size();
47 
48  // Report on the constraints attribute values and ranges
49 
50  DbgLv(0) << parameters[ "DC_model" ] << "Constraints --";
51 
52  QString attrtype = tr( "UNKNOWN" );
53  US_dmGA_Constraints::AttribType atype;
54  QVector< US_dmGA_Constraints::Constraint > cnsv;
55 
56  for ( int mcompx = 0; mcompx < wmodel.components.size(); mcompx++ )
57  {
58  int compno = mcompx + 1;
59  int ncompc = constraints.comp_constraints( mcompx, &cnsv, NULL );
60 
61  DbgLv(0) << " Component" << compno << ":";
62 //DbgLv(0) << "mcompx ncompc c0atype" << mcompx << ncompc << cnsv[0].atype
63 // << "_FF0" << US_dmGA_Constraints::ATYPE_FF0;
64 
65  for ( int cx = 0; cx < ncompc; cx++ )
66  {
67  atype = cnsv[ cx ].atype;
68 //DbgLv(0) << " cx" << cx << "atype" << atype;
69  bool floats = cnsv[ cx ].floats;
70  bool logscl = cnsv[ cx ].logscl;
71  double vmin = cnsv[ cx ].low;
72  double vmax = cnsv[ cx ].high;
73  attrtype = tr( "UNKNOWN" );
74 //DbgLv(0) << " cx" << cx << "fl ls mn mx" << floats << logscl << vmin << vmax;
75 
76  if ( atype == US_dmGA_Constraints::ATYPE_S )
77  attrtype = tr( "Segmentation Coefficient" );
78  else if ( atype == US_dmGA_Constraints::ATYPE_FF0 )
79  attrtype = tr( "Frictional Ratio" );
80  else if ( atype == US_dmGA_Constraints::ATYPE_MW )
81  attrtype = tr( "Molecular Weight" );
82  else if ( atype == US_dmGA_Constraints::ATYPE_D )
83  attrtype = tr( "Diffusion Coefficient" );
84  else if ( atype == US_dmGA_Constraints::ATYPE_F )
85  attrtype = tr( "Frictional Coefficient" );
86  else if ( atype == US_dmGA_Constraints::ATYPE_VBAR )
87  attrtype = tr( "Vbar (Specific Density)" );
88  else if ( atype == US_dmGA_Constraints::ATYPE_CONC )
89  attrtype = tr( "Signal Concentration" );
90  else if ( atype == US_dmGA_Constraints::ATYPE_EXT )
91  attrtype = tr( "Extinction" );
92 
93 //DbgLv(0) << " cx" << cx << "attrtype" << attrtype;
94  if ( floats )
95  {
96 //DbgLv(0) << " FLOATS";
97  if ( logscl )
98  {
99 //DbgLv(0) << " LOGSCL";
100  DbgLv(0) << " " << attrtype << "floats from "
101  << vmin << "to" << vmax << "(log scale)";
102  }
103  else
104  {
105 //DbgLv(0) << " !LOGSCL";
106  DbgLv(0) << " " << attrtype << "floats from "
107  << vmin << "to" << vmax;
108  }
109  }
110  else
111  {
112 //DbgLv(0) << " !FLOATS";
113  DbgLv(0) << " " << attrtype << "is fixed at "
114  << vmin;
115  }
116  }
117  }
118 
119  for ( int assocx = 0; assocx < wmodel.associations.size(); assocx++ )
120  {
121  int compno = assocx + 1;
122  int nassoc = constraints.assoc_constraints( assocx, &cnsv, NULL );
123  int ncomp = wmodel.associations[ assocx ].rcomps.size();
124 
125  DbgLv(0) << " Reaction" << compno << ":";
126 
127  if ( ncomp == 2 )
128  {
129  DbgLv(0) << " Reactant is component"
130  << wmodel.associations[ assocx ].rcomps[ 0 ] + 1
131  << ", Product is component"
132  << wmodel.associations[ assocx ].rcomps[ 1 ] + 1;
133  }
134 
135  else
136  {
137  DbgLv(0) << " Reactants are components"
138  << wmodel.associations[ assocx ].rcomps[ 0 ] + 1
139  << " and" << wmodel.associations[ assocx ].rcomps[ 1 ] + 1
140  << ", Product is component"
141  << wmodel.associations[ assocx ].rcomps[ 2 ] + 1;
142  }
143 
144  for ( int rx = 0; rx < nassoc; rx++ )
145  {
146  atype = cnsv[ rx ].atype;
147  bool floats = cnsv[ rx ].floats;
148  bool logscl = cnsv[ rx ].logscl;
149  double vmin = cnsv[ rx ].low;
150  double vmax = cnsv[ rx ].high;
151  attrtype = tr( "UNKNOWN" );
152 
153  if ( atype == US_dmGA_Constraints::ATYPE_KD )
154  attrtype = tr( "K_Dissociation" );
155  else if ( atype == US_dmGA_Constraints::ATYPE_KOFF )
156  attrtype = tr( "k_Off Rate" );
157 
158  if ( floats )
159  {
160  if ( logscl )
161  {
162  DbgLv(0) << " " << attrtype << "floats from "
163  << vmin << "to" << vmax << "(log scale)";
164  }
165  else
166  {
167  DbgLv(0) << " " << attrtype << "floats from "
168  << vmin << "to" << vmax;
169  }
170  }
171  else
172  {
173  DbgLv(0) << " " << attrtype << "is fixed at "
174  << vmin;
175  }
176  }
177  }
178 DbgLv(0) << " wmodel as1 K_d k_off" << wmodel.associations[0].k_d
179  << wmodel.associations[0].k_off;
180 
181  Fitness empty_fitness;
182  empty_fitness.fitness = LARGE;
183  // Get base Gene and set up mutation controls
184  dgene = wmodel;
185  nfloatc = constraints.float_constraints( &cns_flt );
186  nfvari = ( 1 << nfloatc ) - 1;
187  dgmarker.resize( nfloatc );
188  do_astfem = ( wmodel.components[ 0 ].sigma == 0.0 &&
189  wmodel.components[ 0 ].delta == 0.0 &&
190  wmodel.coSedSolute < 0 &&
191  data_sets[ 0 ]->compress == 0.0 );
192  lfvari.fill( 0, nfvari );
193 //DbgLv(0) << " wmodel DUMP";
194 //wmodel.debug();
195 
196  // Initialize arrays
197  for ( int ii = 0; ii < gcores_count; ii++ )
198  {
199  best_dgenes << dgene;
200 
201  empty_fitness.index = ii;
202  best_fitness << empty_fitness;
203  }
204 
205  QDateTime time = QDateTime::currentDateTime();
206 
207  // Handle Monte Carlo iterations. There will always be at least 1.
208  while ( true )
209  {
210  // Compute all generations of an MC iteration
211 
212  dmga_master_loop();
213 
214  // Get the best-fit gene
215 
216  qSort( best_fitness );
217 DbgLv(1) << "GaMast: EOML: fitness sorted bfx" << best_fitness[0].index
218  << "bestdg size" << best_dgenes.size();
219 
220  Fitness* bfit = &best_fitness[ 0 ];
221  DbgLv(0) << "Last generation best RMSD" << sqrt( bfit->fitness );
222 
223  if ( minimize_opt == 2 )
224  { // If gradient search method for all terminations, minimize now
225  in_gsm = TRUE;
226  double fitness = bfit->fitness;
227  dgene = best_dgenes[ bfit->index ];
228 
229  bfit->fitness = minimize_dmga( dgene, fitness );
230 
231  DbgLv(0) << "Post-minimization RMSD" << sqrt( bfit->fitness );
232  best_dgenes[ bfit->index ] = dgene;
233  in_gsm = FALSE;
234  }
235 
236  // Compute the variance (fitness) for the final best-fit model
237 
238 DbgLv(1) << "GaMast: EOML: call calc_residuals";
239  calc_residuals_dmga( 0, data_sets.size(), simulation_values, dgene );
240 DbgLv(1) << "GaMast: EOML: variance" << simulation_values.variance;
241 
242  // Output the model
243 
244  if ( data_sets.size() == 1 )
245  { // Output the single-data model
246 DbgLv(1) << "GaMast: EOML: CALL write_output";
247  write_output();
248  }
249  else
250  { // Output the global model
251  write_global();
252  }
253 
254  // Handle any MonteCarlo iteration logic
255 
256 DbgLv(1) << "GaMast: mc_iter iters" << mc_iteration << mc_iterations;
257  mc_iteration++;
258  if ( mc_iterations > 1 )
259  {
260  qDebug() << "Fit RMSD" << sqrt( simulation_values.variance )
261  << " of MC_Iteration" << mc_iteration;
262  if ( mc_iteration < mc_iterations )
263  {
264  // Set scaled_data the first time
265  if ( mc_iteration <= mgroup_count )
266  {
267  scaled_data = simulation_values.sim_data;
268  }
269 
270  time_mc_iterations();
271 
272 DbgLv(1) << "GaMast: set_gaMC call";
273  set_dmga_MonteCarlo();
274 DbgLv(1) << "GaMast: set_gaMC return";
275  }
276  else
277  break;
278  }
279  else
280  {
281 DbgLv(1) << "GaMast: FFrmsd: variance" << simulation_values.variance;
282  qDebug() << "Final Fit RMSD" << sqrt( simulation_values.variance );
283  break;
284  }
285  }
286 
287  DbgLv(0) << my_rank << ": Master signalling FINISHED to all Demes";
288 
289  // Report on the attribute values in the final model
290 
291  US_Model* fmodel = &data_sets[ 0 ]->model;
292 //*DEBUG
293 double vsum=0.0;
294 for (int ii=0; ii<simulation_values.sim_data.scanCount(); ii++)
295  for (int jj=0; jj<simulation_values.sim_data.pointCount(); jj++ )
296  vsum += sq(simulation_values.sim_data.value(ii,jj));
297 DbgLv(0) << "VSUM=" << vsum;
298 int hh=simulation_values.sim_data.pointCount()/2;
299 int ss=simulation_values.sim_data.scanCount();
300 DbgLv(0) << "SDAT 0-3"
301  << simulation_values.sim_data.value(0,hh)
302  << simulation_values.sim_data.value(1,hh)
303  << simulation_values.sim_data.value(2,hh)
304  << simulation_values.sim_data.value(3,hh);
305 DbgLv(0) << "SDAT *-n"
306  << simulation_values.sim_data.value(ss-4,hh)
307  << simulation_values.sim_data.value(ss-3,hh)
308  << simulation_values.sim_data.value(ss-2,hh)
309  << simulation_values.sim_data.value(ss-1,hh);
310 DbgLv(0) << "MDL comp1 s,k,w,D,f"
311  << fmodel->components[0].s / data_sets[0]->s20w_correction
312  << fmodel->components[0].f_f0
313  << fmodel->components[0].mw
314  << fmodel->components[0].D / data_sets[0]->D20w_correction
315  << fmodel->components[0].f;
316 DbgLv(0) << "MDL comp2 s,k,w,D,f"
317  << fmodel->components[1].s / data_sets[0]->s20w_correction
318  << fmodel->components[1].f_f0
319  << fmodel->components[1].mw
320  << fmodel->components[1].D / data_sets[0]->D20w_correction
321  << fmodel->components[1].f;
322 DbgLv(0) << "MDL asoc1 Kd,koff"
323  << fmodel->associations[0].k_d
324  << fmodel->associations[0].k_off;
325 DbgLv(0) << "DS dens visc manual temp"
326  << data_sets[0]->density
327  << data_sets[0]->viscosity
328  << data_sets[0]->manual
329  << data_sets[0]->temperature;
330 US_SimulationParameters* simpar = &data_sets[0]->simparams;
331 DbgLv(0) << "SIMP simpt mType gType rreso menis"
332  << simpar->simpoints
333  << simpar->meshType
334  << simpar->gridType
335  << simpar->radial_resolution
336  << simpar->meniscus;
337 DbgLv(0) << "SIMP bott temp rnoise tinoise rinoise"
338  << simpar->bottom
339  << simpar->temperature
340  << simpar->rnoise
341  << simpar->tinoise
342  << simpar->rinoise;
343 DbgLv(0) << "SIMP bform bvol bottpos rcoeffs"
344  << simpar->band_forming
345  << simpar->band_volume
346  << simpar->bottom_position
347  << simpar->rotorcoeffs[0]
348  << simpar->rotorcoeffs[1];
349 US_SimulationParameters::SpeedProfile* spstep = &simpar->speed_step[0];
350 DbgLv(0) << "STEP0 durmin dlymin w2tf w2tl timf timl"
351  << spstep->duration_minutes
352  << spstep->delay_minutes
353  << spstep->w2t_first
354  << spstep->w2t_last
355  << spstep->time_first
356  << spstep->time_last;
357 DbgLv(0) << "STEP0 speed accel accelf"
358  << spstep->rotorspeed
359  << spstep->acceleration
360  << spstep->acceleration_flag;
361 //*DEBUG
362  DbgLv(0) << fmodel->description;
363  DbgLv(0) << " Final Model Attribute Values --";
364 
365  for ( int mcompx = 0; mcompx < fmodel->components.size(); mcompx++ )
366  {
367  int compno = mcompx + 1;
368  int ncompc = constraints.comp_constraints( mcompx, &cnsv, NULL );
369 
370  DbgLv(0) << " Component" << compno << ":";
371 //DbgLv(0) << " mcompx" << mcompx << "ncompc" << ncompc;
372 //DbgLv(0) << " at[n-2] at[n-1]" << cnsv[ncompc-2].atype << cnsv[ncompc-1].atype;
373 
374  for ( int cx = 0; cx < ncompc; cx++ )
375  {
376  atype = cnsv[ cx ].atype;
377 
378  if ( atype == US_dmGA_Constraints::ATYPE_S )
379  attrtype = tr( "Segmentation Coefficient" );
380  else if ( atype == US_dmGA_Constraints::ATYPE_FF0 )
381  attrtype = tr( "Frictional Ratio" );
382  else if ( atype == US_dmGA_Constraints::ATYPE_MW )
383  attrtype = tr( "Molecular Weight" );
384  else if ( atype == US_dmGA_Constraints::ATYPE_D )
385  attrtype = tr( "Diffusion Coefficient" );
386  else if ( atype == US_dmGA_Constraints::ATYPE_F )
387  attrtype = tr( "Frictional Coefficient" );
388  else if ( atype == US_dmGA_Constraints::ATYPE_VBAR )
389  attrtype = tr( "Vbar (Specific Density)" );
390  else if ( atype == US_dmGA_Constraints::ATYPE_CONC )
391  attrtype = tr( "Signal Concentration" );
392  else if ( atype == US_dmGA_Constraints::ATYPE_EXT )
393  attrtype = tr( "Extinction" );
394  else
395  attrtype = tr( "Unknown" );
396 
397  double aval = constraints.fetch_attrib( fmodel->components[ mcompx ],
398  atype );
399  DbgLv(0) << " " << attrtype << "has a value of"
400  << aval;
401  }
402  }
403 
404  for ( int assocx = 0; assocx < fmodel->associations.size(); assocx++ )
405  {
406  int compno = assocx + 1;
407  int nassoc = constraints.assoc_constraints( assocx, &cnsv, NULL );
408  int ncomp = fmodel->associations[ assocx ].rcomps.size();
409 
410  DbgLv(0) << " Reaction" << compno << ":";
411 
412  if ( ncomp == 2 )
413  {
414  DbgLv(0) << " Reactant is component"
415  << fmodel->associations[ assocx ].rcomps[ 0 ] + 1
416  << ", Product is component"
417  << fmodel->associations[ assocx ].rcomps[ 1 ] + 1;
418  }
419 
420  else
421  {
422  DbgLv(0) << " Reactants are components"
423  << fmodel->associations[ assocx ].rcomps[ 0 ] + 1
424  << " and" << fmodel->associations[ assocx ].rcomps[ 1 ] + 1
425  << ", Product is component"
426  << fmodel->associations[ assocx ].rcomps[ 2 ] + 1;
427  }
428 
429  for ( int rx = 0; rx < nassoc; rx++ )
430  {
431  atype = cnsv[ rx ].atype;
432 
433  if ( atype == US_dmGA_Constraints::ATYPE_KD )
434  attrtype = tr( "K_Dissociation" );
435  else if ( atype == US_dmGA_Constraints::ATYPE_KOFF )
436  attrtype = tr( "k_Off Rate" );
437 
438  double aval = constraints.fetch_attrib( fmodel->associations[ assocx ],
439  atype );
440  DbgLv(0) << " " << attrtype << "has a value of"
441  << aval;
442  }
443  }
444 
445  MPI_Job job;
446 
447  // Send finish to workers ( in the tag )
448  for ( int worker = 1; worker <= my_workers; worker++ )
449  {
450  MPI_Send( &job, // MPI #0
451  sizeof( job ),
452  MPI_BYTE,
453  worker,
454  FINISHED,
455  my_communicator );
456  }
457 }
458 
459 void US_MPI_Analysis::dmga_master_loop( void )
460 {
461  static const double DIGIT_FIT = 1.0e+4;
462  static const int min_generation = 10;
463  static const int _KS_BASE_ = 6;
464  static const int _KS_STEP_ = 3;
465  QString dbgtext = parameters[ "debug_text" ];
466  QString s_ksbase = par_key_value( dbgtext, "ksame_base" );
467  QString s_ksstep = par_key_value( dbgtext, "ksame_step" );
468  int ks_base = s_ksbase.isEmpty() ? _KS_BASE_ : s_ksbase.toInt();
469  int ks_step = s_ksstep.isEmpty() ? _KS_STEP_ : s_ksstep.toInt();
470  ks_base = qMax( 2, ks_base );
471  ks_step = qMax( 1, ks_step );
472 // static const int max_same_count = my_workers * 5;
473  int max_same_count = ( ks_base + ( nfloatc / ks_step ) * ks_step )
474  * my_workers;
475  int avg_generation = 0;
476  bool early_termination = false;
477  int fitness_same_count = 0;
478  double best_overall_fitness = LARGE;
479  int tag;
480  int workers = my_workers;
481 DbgLv(1) << "dmga_master start loop: gcores_count fitsize" << gcores_count
482  << best_fitness.size() << "best_overall" << best_overall_fitness;
483 DbgLv(0) << "dmga_master start loop: nfloatc max_same_count"
484  << nfloatc << max_same_count << "ks_base ks_step" << ks_base << ks_step;
485 
486  // Reset best fitness for each worker
487  for ( int i = 0; i < gcores_count; i++ )
488  {
489  best_fitness[ i ].fitness = LARGE;
490  best_fitness[ i ].index = i;
491  }
492 
493  QList < DGene > emigres; // Holds genes passed as emmigrants
494  QVector< int > v_generations( gcores_count, 0 );
495  int sum = 0;
496  int avg = 0;
497  long rsstotal = 0L;
498  double fit_power = 5;
499  double fit_digit = 1.0e4;
500  double fitness_round = 1.0e5;
501 
502 
503  while ( workers > 0 )
504  {
505  MPI_GA_MSG msg;
506  MPI_Status status;
507  int worker;
508 
509  MPI_Recv( &msg, // Get a message MPI #1
510  sizeof( msg ),
511  MPI_BYTE,
512  MPI_ANY_SOURCE,
513  MPI_ANY_TAG,
514  my_communicator,
515  &status );
516 
517  worker = status.MPI_SOURCE;
518 
519  max_rss();
520 
521  switch ( status.MPI_TAG )
522  {
523  case GENERATION:
524  v_generations[ worker ] = msg.generation;
525 
526  sum = 0;
527  for ( int i = 1; i <= my_workers; i++ )
528  sum += v_generations[ i ];
529 
530  avg = qRound( sum / my_workers ) + 1;
531 
532  if ( avg > avg_generation )
533  {
534  avg_generation = avg;
535  int mc_iter = mgroup_count < 2 ? ( mc_iteration + 1 )
536  : mc_iteration;
537 
538  QString progress =
539  "Avg. Generation: " + QString::number( avg_generation );
540 
541  if ( data_sets.size() > 1 )
542  {
543  if ( datasets_to_process == 1 )
544  progress += "; Dataset: "
545  + QString::number( current_dataset + 1 )
546  + " of " + QString::number( count_datasets );
547  else
548  progress += "; Datasets: "
549  + QString::number( datasets_to_process );
550  }
551 
552  progress += "; MonteCarlo: " + QString::number( mc_iter );
553  if ( best_overall_fitness != LARGE && best_overall_fitness > 0.0 )
554  progress += "; RMSD: "
555  + QString::number( sqrt( best_overall_fitness ) );
556 
557  send_udp( progress );
558  }
559 
560  // Get the best gene for the current generation from the worker
561 DbgLv(1) << " MAST: work" << worker << "Recv#2 nfloatc" << nfloatc;
562  MPI_Recv( dgmarker.data(), // MPI #2
563  nfloatc,
564  MPI_DOUBLE,
565  worker,
566  GENE,
567  my_communicator,
568  MPI_STATUS_IGNORE );
569 
570 DbgLv(1) << " MAST: work" << worker << " dgm size" << dgmarker.size()
571  << "bestdg size" << best_dgenes.size();
572  dgene_from_marker( dgmarker, dgene );
573  best_dgenes[ worker ] = dgene;
574  max_rss();
575 
576  // Compute a current-deme best fitness value that is rounded
577  // to 4 significant digits
578  fit_power = (double)qRound( log10( msg.fitness ) );
579  fit_digit = pow( 10.0, -fit_power ) * DIGIT_FIT;
580  fitness_round = (double)qRound64( msg.fitness * fit_digit )
581  / fit_digit;
582 
583 DbgLv(1) << " MAST: work" << worker << "fit msg,round,bestw,besto"
584  << msg.fitness << fitness_round << best_fitness[worker].fitness
585  << best_overall_fitness;
586  // Set deme's best fitness
587  if ( fitness_round < best_fitness[ worker ].fitness )
588  best_fitness[ worker ].fitness = fitness_round;
589 DbgLv(1) << "master: worker/fitness/best gene" << worker << msg.fitness << dgene_key( best_dgenes[worker] );
590 
591  if ( ! early_termination )
592  { // Handle normal pre-early-termination updates
593  if ( avg_generation == 1 && mc_iterations == 1 &&
594  best_overall_fitness == LARGE )
595  { // Report first best-fit RMSD
596  DbgLv(0) << "First Best Fit RMSD" << sqrt( fitness_round );
597  }
598 DbgLv(1) << " MAST: work" << worker << "fit besto,round" << best_overall_fitness << fitness_round
599  << "fit_power fit_digit msgfit" << fit_power << fit_digit << msg.fitness;
600 
601  if ( fitness_round < best_overall_fitness )
602  { // Update over-all best fitness value (rounded)
603  best_overall_fitness = fitness_round;
604  fitness_same_count = 0;
605  }
606  else
607  { // Bump the count of consecutive same best overall fitness
608  fitness_same_count++;
609  }
610 
611 
612  if ( fitness_same_count > max_same_count &&
613  avg_generation > min_generation )
614  { // Mark early termination at threshold same-fitness count
615  DbgLv(0) << "Fitness has not improved in the last"
616  << fitness_same_count
617  << "deme results - Early Termination.";
618  early_termination = true;
619  }
620 
621  }
622 //DbgLv(1) << " best_overall_fitness" << best_overall_fitness
623 // << "fitness_same_count" << fitness_same_count
624 // << " early_term?" << early_termination;
625 if((worker%10)==1)
626 DbgLv(1) << worker << ": best_overall_fitness" << best_overall_fitness
627  << "fitness_same_count" << fitness_same_count
628  << " early_term?" << early_termination;
629 
630  // Tell the worker to either stop or continue
631  tag = early_termination ? FINISHED : GENERATION;
632 DbgLv(1) << "dgmast: Send#3 tag" << tag << "( FINISHED,GENERATION" << FINISHED << GENERATION << " )";
633 
634  MPI_Send( &msg, // MPI #3
635  0, // Only interested in the tag
636  MPI_BYTE,
637  worker,
638  tag,
639  my_communicator );
640  break;
641 
642  case FINISHED:
643 DbgLv(1) << "dgmast: Recv FINISH msg.size rsstotal" << msg.size << rsstotal;
644  rsstotal += (long)msg.size;
645  workers--;
646  break;
647 
648  case EMMIGRATE:
649  {
650  // First get a set of genes as a concatenated marker vector.
651  int gene_count = msg.size;
652  int doubles_count = gene_count * nfloatc;
653  QVector< double > emmigrants( doubles_count ) ;
654 DbgLv(1) << "dgmast: Recv#4 EMMIG: gene_count doubles_count" << gene_count << doubles_count;
655 
656  MPI_Recv( emmigrants.data(), // MPI #4
657  doubles_count,
658  MPI_DOUBLE,
659  worker,
660  EMMIGRATE,
661  my_communicator,
662  MPI_STATUS_IGNORE );
663 
664  // Add the genes to the emmigres list
665  int emgx = emigres.size();
666 DbgLv(1) << "dgmast: Recv#4 EMMIG: emgx" << emgx;
667  marker_to_dgenes( emmigrants, emigres, emgx, gene_count );
668 DbgLv(1) << "dgmast: Recv#4 EMMIG: emigres size" << emigres.size();
669 
670 //*DEBUG*
671 //if(emigres[0][0].s<0.0)
672 // DbgLv(0) << "MAST: **GENE s" << emigres[0][0].s << " Emigrant";
673 //*DEBUG*
674 
675  max_rss();
676 
677  // Don't send any back if the pool is too small
678  if ( emigres.size() < gene_count * 5 ) doubles_count = 0;
679 
680  // Get immigrants from emmigres
681  QVector< double > immigrants;
682  dgmarker.resize( nfloatc );
683 
684  if ( doubles_count > 0 )
685  {
686  // Prepare a marker vector from the emmigrant list
687  for ( int ii = 0; ii < gene_count; ii++ )
688  {
689  dgene = emigres.takeAt( u_random( emigres.size() ) );
690 
691  marker_from_dgene( dgmarker, dgene );
692 
693  immigrants += dgmarker;
694  }
695  }
696  else
697  {
698 DbgLv(1) << "dgmast: Send#5 IMMIG: count==0 migr data" << immigrants.data();
699  immigrants.resize( 1 );
700  }
701 DbgLv(1) << "dgmast: Send#5 IMMIG: immigrants size" << immigrants.size() << "doubles_count" << doubles_count
702  << "migr data" << immigrants.data();
703 
704  MPI_Send( immigrants.data(), // MPI #5
705  doubles_count,
706  MPI_DOUBLE,
707  worker,
708  IMMIGRATE,
709  my_communicator );
710 DbgLv(1) << "dgmast: Send#5 IMMIG: complete";
711 //*DEBUG*
712 //if(immigrants[0].s<0.0)
713 // DbgLv(0) << "MAST: **GENE s" << immigrants[0].s << " Immigrant-to-send";
714 //*DEBUG*
715  break;
716  }
717  }
718 
719  max_rss();
720  }
721 
722 DbgLv(1) << "Master maxrss" << maxrss << " worker total rss" << rsstotal
723  << "rank" << my_rank;
724  maxrss += rsstotal;
725 
726  if ( early_termination )
727  { // Report when we have reached early termination of generations
728  int mc_iter = mgroup_count < 2 ? ( mc_iteration + 1 ) : mc_iteration;
729  DbgLv(0) << "Early termination at average generation" << avg
730  << ", MC" << mc_iter;
731  }
732 }
733 
734 void US_MPI_Analysis::dmga_global_fit( void )
735 {
736  // This is almost the same as 2dsa global_fit.
737  double concentration = 0.0;
738 
739  // The first dataset is done automatically.
740  for ( int solute = 0; solute < simulation_values.solutes.size(); solute++ )
741  {
742  concentration += simulation_values.solutes[ solute ].c;
743  }
744 
745  // Point to current dataset
746  US_DataIO::EditedData* data = &data_sets[ current_dataset ]->run_data;
747 
748  int scan_count = data->scanCount();
749  int radius_points = data->pointCount();
750  int index = 0;
751 
752  QVector< double > scaled_data( scan_count * radius_points );
753 
754  // Scale the data
755  for ( int s = 0; s < scan_count; s++ )
756  {
757  for ( int r = 0; r < radius_points; r++ )
758  {
759  scaled_data[ index++ ] = data->value( s, r ) / concentration;
760  }
761  }
762 
763  // Send the scaled data to the workers
764  MPI_Job job;
765  job.length = scaled_data.size();
766  job.dataset_offset = current_dataset;
767  job.dataset_count = 1;
768 
769  // Tell each worker that new data coming
770  // Can't use a broadcast because the worker is expecting a Send
771 
772  for ( int worker = 1; worker <= my_workers; worker++ )
773  {
774  MPI_Send( &job, // MPI #7
775  sizeof( MPI_Job ),
776  MPI_BYTE,
777  worker,
778  UPDATE,
779  my_communicator );
780  }
781 
782  // Get everybody synced up
783  MPI_Barrier( my_communicator );
784 
785  MPI_Bcast( scaled_data.data(), // MPI #8
786  scaled_data.size(),
787  MPI_DOUBLE,
788  MPI_Job::MASTER,
789  my_communicator );
790 
791  // Go to the next dataset
792  current_dataset++;
793 
794  // If all datasets have been scaled, do all datasets from now on
795  if ( current_dataset >= data_sets.size() )
796  {
797  datasets_to_process = data_sets.size();
798  current_dataset = 0;
799  }
800 }
801 
802 // Use DMGA residuals as the standard deviation for varying data in MonteCarlo
803 void US_MPI_Analysis::set_dmga_gaussians( void )
804 {
805  sigmas.clear();
806  res_data = &simulation_values.residuals;
807  int ks = 0;
808 
809  for ( int ee = 0; ee < data_sets.size(); ee++ )
810  {
811  int nscan = data_sets[ ee ]->run_data.scanCount();
812  int npoint = data_sets[ ee ]->run_data.pointCount();
813 
814  // Place the residuals magnitudes into a single vector for convenience
815  for ( int ss = 0; ss < nscan; ss++, ks++ )
816  {
817  for ( int rr = 0; rr < npoint; rr++ )
818  {
819  sigmas << qAbs( res_data->value( ks, rr ) );
820  }
821  }
822  }
823 }
824 
825 // Set up next MonteCarlo iteration for dmGA
826 void US_MPI_Analysis::set_dmga_MonteCarlo( void )
827 {
828 DbgLv(1) << "sdMC: mciter" << mc_iteration;
829  // This is almost the same as ga set_monteCarlo
830  if ( mc_iteration <= mgroup_count )
831  {
832  max_depth = 0; // Make the datasets compatible
833  calculated_solutes.clear();
834 
835  set_dmga_gaussians();
836  }
837 
838  mc_data.resize( total_points );
839  int index = 0;
840 
841  // Get a randomized variation of the concentrations
842  // Use a gaussian distribution with the residual as the standard deviation
843  for ( int ee = 0; ee < data_sets.size(); ee++ )
844  {
845  US_DataIO::EditedData* edata = &data_sets[ ee ]->run_data;
846 
847  int scan_count = edata->scanCount();
848  int radius_points = edata->pointCount();
849 
850  for ( int ss = 0; ss < scan_count; ss++ )
851  {
852  for ( int rr = 0; rr < radius_points; rr++ )
853  {
854  double variation = US_Math2::box_muller( 0.0, sigmas[ index ] );
855  mc_data[ index ] = scaled_data.value( ss, rr ) + variation;
856  index++;
857  }
858  }
859  }
860 DbgLv(1) << "sdMC: mc_data set index" << index
861  << "total_points" << total_points;
862 
863  // Broadcast Monte Carlo data to all workers
864  MPI_Job job;
865  job.command = MPI_Job::NEWDATA;
866  job.length = total_points;
867  job.dataset_offset = 0;
868  job.dataset_count = data_sets.size();
869 
870  // Tell each worker that new data coming
871  // Can't use a broadcast because the worker is expecting a Send
872  for ( int worker = 1; worker <= my_workers; worker++ )
873  {
874  MPI_Send( &job, // MPI #9
875  sizeof( job ),
876  MPI_BYTE,
877  worker,
878  UPDATE,
879  my_communicator );
880  }
881 
882  // Get everybody synced up
883  MPI_Barrier( my_communicator );
884 
885  MPI_Bcast( mc_data.data(), // MPI #10
886  total_points,
887  MPI_DOUBLE,
888  MPI_Job::MASTER,
889  my_communicator );
890 }
891 
892 // Get a marker vector of doubles from a dmga gene (model)
893 void US_MPI_Analysis::marker_from_dgene( QVector< double >& dgm, DGene& dg )
894 {
895  // Fetch the current value for each floating attribute and store in marker
896  for ( int ii = 0; ii < nfloatc; ii++ )
897  {
898  fetch_attr_value( dgm[ ii ], dg, cns_flt[ ii ].atype,
899  cns_flt[ ii ].mcompx );
900  }
901 }
902 
903 // Get a dmga gene (model) from a marker vector of doubles
904 void US_MPI_Analysis::dgene_from_marker( QVector< double >& dgm, DGene& dg )
905 {
906  // Initial gene is the base work model
907  dg = wmodel;
908 
909  // Store the current value for each floating attribute into the gene
910  for ( int ii = 0; ii < nfloatc; ii++ )
911  {
912  store_attr_value( dgm[ ii ], dg, cns_flt[ ii ].atype,
913  cns_flt[ ii ].mcompx );
914  }
915 }
916 
917 // Get a marker vector of doubles from multiple dmga genes (models)
918 void US_MPI_Analysis::dgenes_to_marker( QVector< double >& dgm,
919  QList< DGene >& dgenes, const int stgx, const int ngenes )
920 {
921  int mx = 0; // Initial marker index
922  int gx = stgx; // Initial genes index
923  DGene wkgene = wmodel; // Base gene
924 
925  for ( int kg = 0; kg < ngenes; kg++ )
926  { // Fetch each gene and store its float attribute values
927  wkgene = dgenes[ gx++ ];
928 
929  for ( int ii = 0; ii < nfloatc; ii++ )
930  { // Fetch each gene float attribute value into the marker vector
931  fetch_attr_value( dgm[ mx++ ], wkgene, cns_flt[ ii ].atype,
932  cns_flt[ ii ].mcompx );
933  }
934  }
935 }
936 
937 // Get multiple dmga genes (models) from a marker vector of doubles
938 void US_MPI_Analysis::marker_to_dgenes( QVector< double >& dgm,
939  QList< DGene >& dgenes, const int stgx, const int ngenes )
940 {
941  int mx = 0; // Initial marker index
942  int gx = stgx; // Initial genes index
943  int igsize = dgenes.size(); // Initial genes list size
944  DGene wkgene = wmodel; // Base gene
945 
946  for ( int kg = 0; kg < ngenes; kg++ )
947  { // Build and store each gene
948 
949  for ( int ii = 0; ii < nfloatc; ii++, gx++ )
950  { // Store each float attribute value into the gene
951  store_attr_value( dgm[ mx++ ], wkgene, cns_flt[ ii ].atype,
952  cns_flt[ ii ].mcompx );
953  }
954 
955  if ( gx >= igsize )
956  {
957  dgenes << wkgene; // Append the constructed gene
958  igsize++;
959  }
960  else
961  dgenes[ gx ] = wkgene; // Store the constructed gene
962  }
963 }
964 
965 // Store a component/association attribute value
966 bool US_MPI_Analysis::store_attr_value( double& aval, US_Model& model,
967  US_dmGA_Constraints::AttribType& atype, int& mcompx )
968 {
969  bool is_ok = true;
970  US_Model::SimulationComponent* sc = NULL;
971  US_Model::Association* as = NULL;
972 
973  if ( atype < US_dmGA_Constraints::ATYPE_KD )
974  sc = &model.components [ mcompx ];
975  else
976  as = &model.associations[ mcompx ];
977 
978  switch( atype )
979  {
980  case US_dmGA_Constraints::ATYPE_S:
981  sc->s = aval;
982  break;
983  case US_dmGA_Constraints::ATYPE_FF0:
984  sc->f_f0 = aval;
985  break;
986  case US_dmGA_Constraints::ATYPE_MW:
987  sc->mw = aval;
988  break;
989  case US_dmGA_Constraints::ATYPE_D:
990  sc->D = aval;
991  break;
992  case US_dmGA_Constraints::ATYPE_F:
993  sc->f = aval;
994  break;
995  case US_dmGA_Constraints::ATYPE_VBAR:
996  sc->vbar20 = aval;
997  break;
998  case US_dmGA_Constraints::ATYPE_CONC:
999  sc->signal_concentration = aval;
1000  if ( sc->extinction != 0.0 )
1001  sc->molar_concentration = aval / sc->extinction;
1002  break;
1003  case US_dmGA_Constraints::ATYPE_EXT:
1004  sc->extinction = aval;
1005  if ( aval != 0.0 )
1006  sc->molar_concentration = sc->signal_concentration / aval;
1007  break;
1008  case US_dmGA_Constraints::ATYPE_KD:
1009  as->k_d = aval;
1010  break;
1011  case US_dmGA_Constraints::ATYPE_KOFF:
1012  as->k_off = aval;
1013  break;
1014  default:
1015  is_ok = false;
1016  break;
1017  }
1018 
1019  return is_ok;
1020 }
1021 
1022 // Fetch a component/association attribute value
1023 bool US_MPI_Analysis::fetch_attr_value( double& aval, US_Model& model,
1024  US_dmGA_Constraints::AttribType& atype, int& mcompx )
1025 {
1026  bool is_ok = true;
1027  US_Model::SimulationComponent* sc = NULL;
1028  US_Model::Association* as = NULL;
1029 
1030  if ( atype < US_dmGA_Constraints::ATYPE_KD )
1031  sc = &model.components [ mcompx ];
1032  else
1033  as = &model.associations[ mcompx ];
1034 
1035  switch( atype )
1036  {
1037  case US_dmGA_Constraints::ATYPE_S:
1038  aval = sc->s;
1039  break;
1040  case US_dmGA_Constraints::ATYPE_FF0:
1041  aval = sc->f_f0;
1042  break;
1043  case US_dmGA_Constraints::ATYPE_MW:
1044  aval = sc->mw;
1045  break;
1046  case US_dmGA_Constraints::ATYPE_D:
1047  aval = sc->D;
1048  break;
1049  case US_dmGA_Constraints::ATYPE_F:
1050  aval = sc->f;
1051  break;
1052  case US_dmGA_Constraints::ATYPE_VBAR:
1053  aval = sc->vbar20;
1054  break;
1055  case US_dmGA_Constraints::ATYPE_CONC:
1056  aval = sc->signal_concentration;
1057  break;
1058  case US_dmGA_Constraints::ATYPE_EXT:
1059  aval = sc->extinction;
1060  break;
1061  case US_dmGA_Constraints::ATYPE_KD:
1062  aval = as->k_d;
1063  break;
1064  case US_dmGA_Constraints::ATYPE_KOFF:
1065  aval = as->k_off;
1066  break;
1067  default:
1068  is_ok = false;
1069  break;
1070  }
1071 
1072  return is_ok;
1073 }
1074 
1075 // Get an apply-ready model from a dmga gene
1076 void US_MPI_Analysis::model_from_dgene( US_Model& model, DGene& dgene )
1077 {
1078  model = dgene; // Initialize the model
1079 
1080  model.update_coefficients(); // Compute missing coefficients
1081 //*DEBUG*
1082 if(my_rank<1) {
1083 DbgLv(2) << my_rank << "MFG:g-comp1 s,k,w,d,f"
1084  << dgene.components[0].s
1085  << dgene.components[0].f_f0
1086  << dgene.components[0].mw
1087  << dgene.components[0].D
1088  << dgene.components[0].f;
1089 DbgLv(2) << my_rank << "MFG:m-comp1 s,k,w,d,f"
1090  << model.components[0].s
1091  << model.components[0].f_f0
1092  << model.components[0].mw
1093  << model.components[0].D
1094  << model.components[0].f;
1095 }
1096 //*DEBUG*
1097 
1098  for ( int ii = 0; ii < model.associations.size(); ii++ )
1099  { // Modify some coefficients based on associations
1100  US_Model::Association* as = &model.associations[ ii ];
1101 
1102  int nrco = as->rcomps.size(); // Components involved
1103  int rpx = nrco - 1; // Last comp index
1104  int rc1 = as->rcomps[ 0 ]; // Comp indexes
1105  int rc2 = nrco == 2 ? rc1 : as->rcomps[ 1 ];
1106  int rcp = as->rcomps[ rpx ];
1107  int st1 = as->stoichs[ 0 ]; // Stoichiometries
1108  int st2 = qAbs( as->stoichs[ 1 ] );
1109  int stp = qAbs( as->stoichs[ rpx ] );
1110  double ff0p = qMax( 1.0, model.components[ rcp ].f_f0 );
1111  model.components[ rcp ] = dgene.components[ rcp ]; // Raw product comp
1112 if(my_rank<2)
1113 DbgLv(1) << my_rank << "MFG: ii" << ii << "nrco,rc1,rc2,rcp"
1114  << nrco << rc1 << rc2 << rcp << "st1,st2,stp" << st1 << st2 << stp;
1115 
1116  // Reset concentration for reactant(s) and product
1117  double cval = model.components[ rc1 ].signal_concentration;
1118  model.components[ rc2 ].signal_concentration = cval;
1119  //model.components[ rcp ].signal_concentration = cval;
1120  model.components[ rcp ].signal_concentration = 0.0;
1121 if(my_rank<2)
1122 DbgLv(1) << my_rank << "MFG: orig cval wval vval" << cval
1123  << model.components[rcp].mw << model.components[rcp].vbar20;
1124 
1125  // Reset molecular weight and vbar for product
1126  double wval = model.components[ rc1 ].mw * (double)st1;
1127  double vsum = model.components[ rc1 ].vbar20 * wval;
1128  double esum = model.components[ rc1 ].extinction * (double)st1;
1129  double wsum = wval;
1130 
1131  if ( nrco > 2 )
1132  {
1133  wval = model.components[ rc2 ].mw * (double)st2;
1134  vsum += model.components[ rc2 ].vbar20 * wval;
1135  esum += model.components[ rc2 ].extinction * (double)st2;
1136  wsum += wval;
1137  }
1138 
1139  model.components[ rcp ].vbar20 = vsum / wsum;
1140  model.components[ rcp ].mw = wsum;
1141  model.components[ rcp ].f_f0 = ff0p;
1142  model.components[ rcp ].s = 0.0;
1143  model.components[ rcp ].D = 0.0;
1144  model.components[ rcp ].f = 0.0;
1145  model.components[ rcp ].extinction = esum;
1146 
1147  // Recompute coefficients with specified vbar,mw,f/f0
1148  model.calc_coefficients( model.components[ rcp ] );
1149 if(my_rank<2)
1150 DbgLv(1) << my_rank << "MFG: rcp" << rcp << "c.s c.k c.mw c.vb"
1151  << model.components[rcp].s << model.components[rcp].f_f0
1152  << model.components[rcp].mw << model.components[rcp].vbar20;
1153  }
1154 }
1155 
Generated on Sun Jun 21 2015 11:40:40 for UltraScan III by   doxygen 1.8.3.1-20130324