UltraScan III
us_model.cpp
Go to the documentation of this file.
1 
3 #include "us_model.h"
4 #include "us_constants.h"
5 #include "us_settings.h"
6 #include "us_util.h"
7 #include "us_math2.h"
8 
9 US_Model::SimulationComponent::SimulationComponent()
10 {
11  name = "New Component";
12  analyteGUID.clear();
13  molar_concentration = 0.0;
14  signal_concentration = 1.0;
15  vbar20 = TYPICAL_VBAR;
16  mw = 50000.0;
17  s = 0.0;
18  D = 0.0;
19  f = 1.0;
20  f_f0 = 1.25;
21  extinction = 0.0;
22  sigma = 0.0;
23  delta = 0.0;
24  oligomer = 1;
25  shape = SPHERE;
26  axial_ratio = 10.0;
27  analyte_type = 0; // Protein
28  c0.radius .clear();
29  c0.concentration.clear();
30 }
31 
32 bool US_Model::SimulationComponent::operator==
33  ( const SimulationComponent& sc ) const
34 {
35  if ( analyteGUID != sc.analyteGUID ) return false;
36  if ( name != sc.name ) return false;
37  if ( molar_concentration != sc.molar_concentration ) return false;
38  if ( signal_concentration != sc.signal_concentration ) return false;
39  if ( vbar20 != sc.vbar20 ) return false;
40  if ( mw != sc.mw ) return false;
41  if ( s != sc.s ) return false;
42  if ( D != sc.D ) return false;
43  if ( f != sc.f ) return false;
44  if ( f_f0 != sc.f_f0 ) return false;
45  if ( extinction != sc.extinction ) return false;
46  if ( sigma != sc.sigma ) return false;
47  if ( delta != sc.delta ) return false;
48  if ( oligomer != sc.oligomer ) return false;
49  if ( shape != sc.shape ) return false;
50  if ( axial_ratio != sc.axial_ratio ) return false;
51  if ( analyte_type != sc.analyte_type ) return false;
52 
53  if ( c0.radius != sc.c0.radius ) return false;
54  if ( c0.concentration != sc.c0.concentration ) return false;
55 
56  return true;
57 }
58 
59 US_Model::Association::Association()
60 {
61  k_d = 0.0;
62  k_off = 0.0;
63  rcomps .clear();
64  stoichs.clear();
65 }
66 
67 bool US_Model::Association::operator== ( const Association& a ) const
68 {
69  if ( k_d != a.k_d ) return false;
70  if ( k_off != a.k_off ) return false;
71  if ( rcomps != a.rcomps ) return false;
72  if ( stoichs != a.stoichs ) return false;
73 
74  return true;
75 }
76 
77 // Model constructor
78 US_Model::US_Model()
79 {
80  monteCarlo = false;
81  wavelength = 0.0;
82  variance = 0.0;
83  meniscus = 0.0;
84  alphaRP = 0.0;
85  subGrids = 0;
86  description = "New Model";
87  optics = ABSORBANCE;
88  analysis = MANUAL;
89  global = NONE;
90  nmcixs = 0;
91 
92  coSedSolute = -1;
93  modelGUID .clear();
94  editGUID .clear();
95  requestGUID .clear();
96  components .clear();
97  associations.clear();
98  dataDescrip .clear();
99  mcixmls .clear();
100 }
101 
102 // Equality operator
103 bool US_Model::operator== ( const US_Model& m ) const
104 {
105  if ( monteCarlo != m.monteCarlo ) return false;
106  if ( wavelength != m.wavelength ) return false;
107  if ( variance != m.variance ) return false;
108  if ( meniscus != m.meniscus ) return false;
109  if ( alphaRP != m.alphaRP ) return false;
110  if ( modelGUID != m.modelGUID ) return false;
111  if ( editGUID != m.editGUID ) return false;
112  if ( requestGUID != m.requestGUID ) return false;
113  if ( description != m.description ) return false;
114  if ( optics != m.optics ) return false;
115  if ( analysis != m.analysis ) return false;
116  if ( global != m.global ) return false;
117  if ( coSedSolute != m.coSedSolute ) return false;
118  if ( subGrids != m.subGrids ) return false;
119  if ( dataDescrip != m.dataDescrip ) return false;
120  if ( associations.size() != m.associations.size() ) return false;
121 
122  for ( int i = 0; i < associations.size(); i++ )
123  if ( associations[ i ] != m.associations[ i ] ) return false;
124 
125  if ( components.size() != m.components.size() ) return false;
126 
127  for ( int i = 0; i < components.size(); i++ )
128  if ( components[ i ] != m.components[ i ] ) return false;
129 
130  return true;
131 }
132 
133 // Update any missing coefficient values in the components of the model
134 bool US_Model::update_coefficients()
135 {
136  bool ok = true;
137 
138  // Calculate missing coefficients for each component; note overall success
139  for ( int ii = 0; ii < components.size(); ii++ )
140  ok = ok && calc_coefficients( components[ ii ] );
141 
142  return ok;
143 }
144 
145 // Calculate any missing coefficient values in a model component
146 bool US_Model::calc_coefficients( SimulationComponent& component )
147 {
148  bool ok = true;
149  double vbar; // component vbar
150  double volume; // e.g., vbar * mw / AVOGADRO
151  double vol_fac; // volume factor, e.g., 0.75 / M_PI
152  double radius_sphere; // e.g., pow( volume * vol_fac, 1/3 );
153  double rsph_fac; // radius_sphere factor; e.g., 0.06 * PI * VISC
154  double onethird; // one third ( 1.0 / 3.0 )
155  double c; // concentration
156  double t; // temperature in kelvin
157  double s; // sedimentation coefficient
158  double D; // diffusion coefficient
159  double mw; // molecular weight
160  double f; // frictional coefficient
161  double fv; // frictional coefficient (working value)
162  double f_f0; // frictional ratio
163  double f0; // f-zero
164  double s20w;
165  double buoyancyb;
166 
167  // Ensure that we have a vbar we can use
168  vbar = component.vbar20;
169 
170  if ( vbar <= 0.0 )
171  {
172  vbar = TYPICAL_VBAR;
173  component.vbar20 = vbar;
174  }
175 
176  t = K20; // temperature kelvin of 20 degr. C
177  vol_fac = 0.75 / M_PI; // various factors used in calcs below
178  onethird = 1.0 / 3.0;
179  rsph_fac = 0.06 * M_PI * VISC_20W;
180  buoyancyb = 1.0 - vbar * DENS_20W;
181 
182  s = component.s; // component coefficient values
183  D = component.D;
184  mw = component.mw;
185  f = component.f;
186  f_f0 = component.f_f0;
187  c = component.signal_concentration;
188  fv = f;
189 
190  // Start with already calculated s if possible
191  if ( s != 0.0 )
192  {
193  s20w = s;
194  double ssgn = ( s < 0.0 ) ? -1.0 : 1.0;
195 
196  // First check s and k (f_f0)
198  if ( f_f0 != 0.0 ) // s and f_f0
199  {
200  double numer = 0.02 * s * f_f0 * vbar * VISC_20W;
201  numer *= ssgn;
202  f0 = 0.09 * VISC_20W * M_PI * sqrt( numer / buoyancyb );
203  fv = f_f0 * f0;
204  D = R * t / ( AVOGADRO * fv );
205  mw = s * R * t / ( D * buoyancyb );
206  mw *= ssgn;
207  }
208 
209  // Next check s and D
211  else if ( D != 0.0 ) // s and D
212  {
213  mw = s * R * t / ( D * buoyancyb );
214  mw *= ssgn;
215  volume = vbar * mw / AVOGADRO;
216  radius_sphere = pow( volume * vol_fac, onethird );
217  f0 = radius_sphere * rsph_fac;
218  fv = mw * buoyancyb / ( s20w * AVOGADRO );
219  fv *= ssgn;
220  double ff0sv = f_f0;
221  f_f0 = fv / f0;
222  double ffdif = qAbs( ff0sv - f_f0 );
223  f_f0 = ( ffdif < 1.e-5 ) ? ff0sv : f_f0;
224  }
225 
226  // Then check any other s + combinations
228  else if ( mw != 0.0 ) // s and mw
229  {
230  D = s * R * t / ( buoyancyb * mw );
231  D *= ssgn;
232  fv = mw * buoyancyb / ( s20w * AVOGADRO );
233  fv *= ssgn;
234  volume = vbar * mw / AVOGADRO;
235  radius_sphere = pow( volume * vol_fac, onethird );
236  f0 = radius_sphere * rsph_fac;
237  f_f0 = fv / f0;
238  }
240  else if ( f != 0.0 ) // s and f
241  {
242  D = R * t / ( AVOGADRO * fv );
243  mw = s * R * t / ( D * buoyancyb );
244  mw *= ssgn;
245  volume = vbar * mw / AVOGADRO;
246  radius_sphere = pow( volume * vol_fac, onethird );
247  f0 = radius_sphere * rsph_fac;
248  f_f0 = fv / f0;
249  }
250  //****************************
251  else // do not have 2 valid coeffs
252  ok = false; //****************************
253  }
254 
255  else if ( component.mw != 0.0 )
256  {
257  volume = vbar * mw / AVOGADRO;
258  radius_sphere = pow( volume * vol_fac, onethird );
259  f0 = radius_sphere * rsph_fac;
261  if ( D != 0.0 ) // mw and D
262  {
263  s = D * buoyancyb * mw / ( R * t );
264  fv = mw * buoyancyb / ( s * AVOGADRO );
265  f_f0 = fv / f0;
266  }
268  else if ( f_f0 != 0.0 ) // mw and f_f0
269  {
270  fv = f_f0 * f0;
271  s = mw * buoyancyb / ( AVOGADRO * fv );
272  D = s * R * t / ( buoyancyb * mw );
273  }
275  else if ( f != 0.0 ) // mw and f
276  {
277  f_f0 = fv / f0;
278  s = mw * buoyancyb / ( AVOGADRO * fv );
279  D = s * R * t / ( buoyancyb * mw );
280  }
281  //****************************
282  else // do not have 2 valid coeffs
283  ok = false; //****************************
284  }
285 
286  else if ( component.D != 0.0 )
287  {
288  if ( f_f0 >= 1.0 ) // D and f_f0
289  {
290  fv = R * t / ( AVOGADRO * D );
291  f0 = fv / f_f0;
292  radius_sphere = f0 / ( 0.06 * M_PI * VISC_20W );
293  double volume = ( 4.0 / 3.0 ) * M_PI * pow( radius_sphere, 3.0 );
294  mw = volume * AVOGADRO / vbar;
295  s = mw * buoyancyb / ( AVOGADRO * fv );
296  }
297  //****************************
298  // D and f - not a valid combo
299  //****************************
300  else // do not have 2 valid coeffs
301  ok = false; //****************************
302  }
304  else if ( fv > 0.0 && f_f0 >= 1.0 ) // f and f_f0
305  {
306  f0 = fv / f_f0;
307  D = R * t / ( AVOGADRO * fv );
308  radius_sphere = f0 / ( 0.06 * M_PI * VISC_20W );
309  double volume = ( 4.0 / 3.0 ) * M_PI * pow( radius_sphere, 3.0 );
310  mw = volume * AVOGADRO / vbar;
311  s = mw * buoyancyb / ( AVOGADRO * fv );
312  }
313  //****************************
314  else // do not have 2 valid coeffs
315  ok = false; //****************************
316 
317  double df = qAbs( f - fv ) / fv;
318 
319  if ( df > 0.1 )
320  { // Significant change in f: replace and use old as concentration
321 
322  if ( c == 0.0 )
323  component.signal_concentration = f;
324 
325  f = fv;
326  }
327 
328  component.s = s; // set component properties (some re-calculated)
329  component.D = D;
330  component.mw = mw;
331  component.f_f0 = f_f0;
332  component.f = f;
333 
334  return ok;
335 }
336 
337 // Test the existence of the models directory path and create it if need be
338 bool US_Model::model_path( QString& path, bool is_perm )
339 {
340 #ifdef NO_DB
341  path = is_perm ? "" : "";
342 #else
343  QDir dir;
344  path = is_perm ? US_Settings::dataDir() + "/models"
345  : US_Settings::tmpDir() + "/temp_models";
346 
347  if ( ! dir.exists( path ) )
348  {
349  if ( ! dir.mkpath( path ) )
350  {
351  return false;
352  }
353  }
354 #endif
355 
356  return true;
357 }
358 
359 // Short text string describing the model type
360 QString US_Model::typeText( int subtype )
361 {
362  struct typemap
363  {
364  AnalysisType typeval;
365  QString typedesc;
366  };
367 
368  const typemap tmap[] =
369  {
370  { MANUAL, QObject::tr( "Manual" ) },
371  { TWODSA, QObject::tr( "2DSA" ) },
372  { TWODSA_MW, QObject::tr( "2DSA-MW" ) },
373  { GA, QObject::tr( "GA" ) },
374  { GA_MW, QObject::tr( "GA-MW" ) },
375  { COFS, QObject::tr( "COFS" ) },
376  { FE, QObject::tr( "FE" ) },
377  { PCSA, QObject::tr( "PCSA" ) },
378  { CUSTOMGRID, QObject::tr( "CUSTOMGRID" ) },
379  { DMGA, QObject::tr( "DMGA" ) },
380  { DMGA_CONSTR, QObject::tr( "DMGA_CONSTR" ) }
381  };
382 
383  const int ntmap = sizeof( tmap ) / sizeof( tmap[ 0 ] );
384 
385  QString tdesc = tmap[ 0 ].typedesc;
386 
387  for ( int jj = 0; jj < ntmap; jj++ )
388  { // look for model type match
389 
390  if ( analysis == tmap[ jj ].typeval )
391  { // we have a match: build type description
392  tdesc = tmap[ jj ].typedesc; // set basic model analyis type
393 
394  if ( associations.size() > 0 ) // Reversible Associations subtype
395  tdesc = tdesc + "-RA";
396 
397  if ( alphaRP != 0.0 )
398  tdesc = tdesc + "-TR"; // Tikhonov Regularization subtype
399 
400  if ( global == MENISCUS ) // Fit Meniscus subtype
401  tdesc = tdesc + "-FM";
402 
403  else if ( global == GLOBAL ) // Global subtype
404  tdesc = ( jj > 0 ) ? tdesc + "-GL" : "Global";
405 
406  else if ( global == SUPERGLOBAL ) // SuperGlobal subtype
407  tdesc = ( jj > 0 ) ? tdesc + "-SG" : "SuperGlobal";
408 
409  if ( analysis == PCSA ) // Add sub-type (SL,IS,DS,HL) to PCSA
410  {
411  if ( subtype > 0 )
412  { // Append sub-type based on a given flag
413  if ( subtype == 1 )
414  tdesc = tdesc + "-SL";
415  else if ( subtype == 2 )
416  tdesc = tdesc + "-IS";
417  else if ( subtype == 4 )
418  tdesc = tdesc + "-DS";
419  else if ( subtype == 8 )
420  tdesc = tdesc + "-HL";
421  else
422  tdesc = tdesc + "-IS";
423  }
424 
425  else
426  { // Append sub-type based on already-created description
427  int kk = description.indexOf( "_PCSA" );
428  if ( kk > 0 ) // Append "-SL"|"-IS"|"-DS"|"-HL"
429  tdesc = tdesc + description.mid( kk + 5, 3 );
430  else // By default, assume "-IS"
431  tdesc = tdesc + "-IS";
432  }
433  }
434 
435  if ( monteCarlo ) // Monte Carlo subtype
436  tdesc = tdesc + "-MC";
437 
438  break;
439  }
440  }
441 
442  return tdesc; // return type description text
443 }
444 
445 // Flag whether model component frictional ratios are constant
446 bool US_Model::constant_ff0( void )
447 {
448  double valmin = components[ 0 ].f_f0;
449  double valmax = components[ 0 ].f_f0;
450 
451  for ( int ii = 1; ii < components.size(); ii++ )
452  {
453  valmin = qMin( valmin, components[ ii ].f_f0 );
454  valmax = qMax( valmax, components[ ii ].f_f0 );
455  }
456 
457  return ( ( valmax - valmin ) < 1.0e-3 );
458 }
459 
460 // Flag whether model component vbar values are constant
461 bool US_Model::constant_vbar( void )
462 {
463  double valmin = components[ 0 ].vbar20;
464  double valmax = components[ 0 ].vbar20;
465 
466  for ( int ii = 1; ii < components.size(); ii++ )
467  {
468  valmin = qMin( valmin, components[ ii ].vbar20 );
469  valmax = qMax( valmax, components[ ii ].vbar20 );
470  }
471 
472  return ( ( valmax - valmin ) < 1.0e-4 );
473 }
474 
475 // Flag whether a model component is a reactant
476 bool US_Model::is_reactant( const int compx )
477 {
478  bool is_react = false;
479 
480  for ( int ii = 0; ii < associations.size(); ii++ )
481  {
482  Association* as = &associations[ ii ];
483  int rcx = as->rcomps.indexOf( compx );
484 
485  if ( rcx >= 0 && as->stoichs[ rcx ] >= 0 )
486  {
487  is_react = true;
488  break;
489  }
490  }
491 
492  return is_react;
493 }
494 
495 // Flag whether a model component is a product of a reaction
496 bool US_Model::is_product( const int compx )
497 {
498  bool is_prod = false;
499 
500  for ( int ii = 0; ii < associations.size(); ii++ )
501  {
502  Association* as = &associations[ ii ];
503  int rcx = as->rcomps.indexOf( compx );
504 
505  if ( rcx >= 0 && as->stoichs[ rcx ] < 0 )
506  {
507  is_prod = true;
508  break;
509  }
510  }
511 
512  return is_prod;
513 }
514 
515 // Load a model from DB or local file
516 int US_Model::load( bool db_access, const QString& guid, US_DB2* db )
517 {
518  if ( db_access ) return load_db ( guid, db );
519  else return load_disk( guid );
520 }
521 
522 // Load model from local disk
523 int US_Model::load_disk( const QString& guid )
524 {
525  int error = US_DB2::ERROR; // Error by default
526 
527  QString path;
528 
529  if ( ! model_path( path ) )
530  {
531  message = QObject::tr ( "Could not create model directory" );
532  return error;
533  }
534 
535  QDir f( path );
536  QStringList filter( "M*.xml" );
537  QStringList names = f.entryList( filter, QDir::Files, QDir::Name );
538  QString filename;
539  bool found = false;
540 
541  for ( int i = 0; i < names.size(); i++ )
542  {
543  filename = path + "/" + names[ i ];
544  QFile file( filename );
545 
546  if ( ! file.open( QIODevice::ReadOnly | QIODevice::Text) ) continue;
547 
548  QXmlStreamReader xml( &file );
549 
550  while ( ! xml.atEnd() )
551  {
552  xml.readNext();
553 
554  if ( xml.isStartElement() )
555  {
556  if ( xml.name() == "model" )
557  {
558  QXmlStreamAttributes a = xml.attributes();
559 
560  if ( a.value( "modelGUID" ).toString() == guid )
561  found = true;
562  break;
563  }
564  }
565  }
566 
567  file.close();
568 
569  if ( found ) return load( filename );
570  }
571 
572  message = QObject::tr ( "Could not find analyte guid" );
573  return error;
574 }
575 
576 // Load a model from a local file
577 int US_Model::load( const QString& filename )
578 {
579  QFile file( filename );
580 
581  if ( ! file.open( QIODevice::ReadOnly | QIODevice::Text) )
582  return US_DB2::ERROR;
583 
584  QXmlStreamReader xml( &file );
585 
586  int result = load_stream( xml );
587 
588  if ( result == US_DB2::NO_MODEL && monteCarlo )
589  { // Handle a multi-model stream
590  file.close();
591  file.open( QIODevice::ReadOnly | QIODevice::Text );
592 
593  QTextStream tsi( &file );
594 
595  result = load_multi_model( tsi );
596  }
597 
598  file.close();
599  return result;
600 }
601 
602 // Load a model from an XML string
603 int US_Model::load_string( const QString& mcont )
604 {
605  QXmlStreamReader xml( mcont );
606  int result = load_stream( xml );
607  return result;
608 }
609 
610 // Load a model from an XML stream
611 int US_Model::load_stream( QXmlStreamReader& xml )
612 {
613  QString coSedStr;
614  QString comprStr;
615 
616  components .clear();
617  associations.clear();
618 
619  QXmlStreamAttributes a;
620  bool read_next = true;
621  int nmtag = 0;
622 
623  while ( ! xml.atEnd() )
624  {
625  if ( read_next ) xml.readNext();
626  read_next = true;
627 
628  if ( xml.isStartElement() )
629  {
630  if ( xml.name() == "model" )
631  {
632  nmtag++;
633 
634  if ( nmtag > 1 )
635  { // A second model tag: return to handle multi-model stream
636  monteCarlo = true;
637  return US_DB2::NO_MODEL;
638  }
639 
640  a = xml.attributes();
641 
642  QString mcst = a.value( "monteCarlo" ).toString();
643  monteCarlo = ( ! mcst.isEmpty() && mcst != "0" );
644  QString iter = a.value( "iterations" ).toString();
645  monteCarlo = monteCarlo ? monteCarlo :
646  ( iter.isEmpty() ? false : ( iter.toInt() > 1 ) );
647  wavelength = a.value( "wavelength" ).toString().toDouble();
648  QString vari = a.value( "variance" ).toString();
649  variance = vari.isEmpty() ? 0.0 : vari.toDouble();
650  QString meni = a.value( "meniscus" ).toString();
651  meniscus = meni.isEmpty() ? 0.0 : meni.toDouble();
652  QString alph = a.value( "alphaRP" ).toString();
653  alphaRP = alph.isEmpty() ? 0.0 : alph.toDouble();
654  description = a.value( "description" ).toString();
655  modelGUID = a.value( "modelGUID" ).toString();
656  editGUID = a.value( "editGUID" ).toString();
657  requestGUID = a.value( "requestGUID" ).toString();
658  dataDescrip = a.value( "dataDescrip" ).toString();
659  coSedStr = a.value( "coSedSolute" ).toString();
660  coSedSolute = ( coSedStr.isEmpty() ) ? -1 : coSedStr.toInt();
661  QString subgs = a.value( "subGrids" ).toString();
662  subGrids = subgs.isEmpty() ? subGrids
663  : subgs.toInt();
664  QString anal1 = a.value( "type" ).toString();
665  QString anal2 = a.value( "analysisType" ).toString();
666  analysis = anal1.isEmpty() ? analysis
667  : (AnalysisType)anal1.toInt();
668  analysis = anal2.isEmpty() ? analysis
669  : (AnalysisType)anal2.toInt();
670  global =
671  (GlobalType) a.value( "globalType" ).toString().toInt();
672  optics =
673  (OpticsType) a.value( "opticsType" ).toString().toInt();
674  }
675 
676  else if ( xml.name() == "analyte" )
677  {
678  SimulationComponent sc;
679  a = xml.attributes();
680 
681  sc.analyteGUID = a.value( "analyteGUID" ).toString();
682 
683  sc.name = a.value( "name" ).toString();
684  QString avbar = a.value( "vbar20" ).toString();
685  sc.vbar20 = avbar.isEmpty() ? TYPICAL_VBAR : avbar.toDouble();
686  sc.mw = a.value( "mw" ).toString().toDouble();
687  sc.s = a.value( "s" ).toString().toDouble();
688  sc.D = a.value( "D" ).toString().toDouble();
689  sc.f = a.value( "f" ).toString().toDouble();
690  sc.f_f0 = a.value( "f_f0" ).toString().toDouble();
691  sc.extinction = a.value( "extinction" ).toString().toDouble();
692  QString aaxia = a.value( "axial" ).toString();
693  sc.axial_ratio = aaxia.isEmpty() ? 10.0 : aaxia.toDouble();
694  sc.sigma = a.value( "sigma" ).toString().toDouble();
695  sc.delta = a.value( "delta" ).toString().toDouble();
696 
697  sc.molar_concentration = a.value( "molar" ).toString().toDouble();
698  sc.signal_concentration = a.value( "signal" ).toString().toDouble();
699  sc.oligomer = a.value( "oligomer" ).toString().toInt();
700 
701  if ( sc.oligomer < 1 )
702  {
703  sc.oligomer = a.value( "stoich" ).toString().toInt();
704  sc.oligomer = ( sc.oligomer > 0 ) ? sc.oligomer : 1;
705  }
706 
707  sc.shape =
708  (ShapeType)a.value( "shape" ).toString().toInt();
709  sc.analyte_type = a.value( "type" ).toString().toInt();
710 
711  mfem_scans( xml, sc );
712 
713  read_next = false; // Skip the next read
714 
715  components << sc;
716  }
717 
718  else if ( xml.name() == "association" )
719  {
720  Association as;
721  get_associations( xml, as );
722 
723  associations << as;
724 
725  read_next = false; // Skip the next read
726  }
727  }
728  }
729 
730  if ( US_Settings::us_debug() > 2 ) debug();
731 
732  return US_DB2::OK;
733 }
734 
735 // Load from a multiple-model stream and create an MC composite model
736 int US_Model::load_multi_model( QTextStream& tsi )
737 {
738  int result = US_DB2::OK;
739  QString mline, mdesc, mcont;
740  nmcixs = 0;
741  mcixmls.clear();
742 
743  // Read and save the first three XML lines
744  QString line1 = tsi.readLine() + "\n";
745  QString line2 = tsi.readLine() + "\n";
746  QString line3 = tsi.readLine() + "\n";
747 
748  // Read remaining lines and save description,contents of all models
749  while ( ! tsi.atEnd() )
750  {
751  mline = tsi.readLine();
752  if ( mline.contains( "</ModelData>" ) )
753  break;
754 
755  if ( mline.contains( "<model " ) )
756  { // At model tag, create initial contents
757  mcont = line1 + line2 + line3 + mline + "\n";
758  // Parse description and save it, if first iteration
759  if ( nmcixs == 0 )
760  {
761  int idx = qMax( mline.indexOf( "description=" ), 0 );
762  mdesc = QString( mline ).mid( idx, 99 ).section( "\"", 1, 1 );
763  }
764  }
765 
766  else if ( mline.contains( "</model>" ) )
767  { // At end of model section, save a model content and bump count
768  mcont = mcont + mline + "\n</ModelData>\n";
769  mcixmls << mcont;
770  nmcixs++;
771  }
772 
773  else
774  { // For any other type of line, add it to contents;
775  mcont = mcont + mline + "\n";
776  }
777  }
778 
779  QVector< SimulationComponent > mmcomps;
780  QVector< Association > mmassos;
781  int ncnstv = 0;
782  int ncnstk = 0;
783 
784  // Build composite components and associations
785  for ( int ii = 0; ii < nmcixs; ii++ )
786  {
787  mcont = mcixmls[ ii ];
788  QXmlStreamReader xml( mcont );
789  load_stream( xml );
790 
791  mmcomps << components;
792  mmassos << associations;
793 
794  if ( constant_vbar() ) ncnstv++;
795  if ( constant_ff0() ) ncnstk++;
796  }
797 
798  // Compress and scale components to only unique solute points
799 
800  components .clear();
801  associations.clear();
802  double sclnrm = 1.0 / (double)nmcixs; // Scale for concentrations
803  bool cnst_vb = ( ncnstv >= ncnstk ); // Flag for constant vbar
804  QStringList sklist; // List of all solute points
805  QStringList skvals; // List of unique solute points
806 
807  for ( int ii = 0; ii < mmcomps.size(); ii++ )
808  { // Build list of solute point strings and list of unique ones
809  double sval = mmcomps[ ii ].s * 1.0e+13;
810  double kval = cnst_vb ? mmcomps[ ii ].f_f0 : mmcomps[ ii ].vbar20;
811  QString skval = QString().sprintf( "%10.4f %8.5f", sval, kval );
812  sklist << skval;
813  if ( ! skvals.contains( skval ) )
814  skvals << skval;
815  }
816 
817  int nskl = sklist.size();
818  int nskv = skvals.size();
819  skvals.sort(); // Sort solute points
820  SimulationComponent scomp;
821 
822  for ( int ii = 0; ii < nskv; ii ++ )
823  { // Average concentration at each unique solute point
824  QString skval = skvals[ ii ]; // Identifying solute point string
825  double conc = 0.0;
826 
827  for ( int jj = 0; jj < nskl; jj++ )
828  { // Search all solute points
829  if ( skval == sklist[ jj ] )
830  { // If a match, sum the concentration
831  scomp = mmcomps[ jj ];
832  conc += scomp.signal_concentration;
833  }
834  }
835 
836  scomp.name = QString().sprintf( "SC%04d", ii + 1 );
837  scomp.signal_concentration = conc * sclnrm;
838  components << scomp;
839  }
840 
841  if ( mmassos.size() > 0 )
842  associations << mmassos[ 0 ];
843 
844  QString mdsc1 = QString( mdesc ).section( ".", 0, -3 );
845  QString mdsc2 = QString( mdesc ).section( ".", -2, -2 )
846  .section( "_", 0, -2 );
847  QString mdsc3 = QString( mdesc ).section( ".", -1, -1 );
848  QString miter = QString().sprintf( "_mcN%03i", nmcixs );
849  description = mdsc1 + "." + mdsc2 + miter + "." + mdsc3;
850 //qDebug() << "MDL:LMM: miter" << miter << "desc" << description;
851 
852  return result;
853 }
854 
855 // Write a multiple-model stream
856 void US_Model::write_mm_stream( QTextStream& tso )
857 {
858  if ( ! monteCarlo || nmcixs < 1 ) // Do nothing if no MC iterations
859  return;
860 
861  // Build an output stream from MC iteration input streams
862  for ( int ii = 0; ii < nmcixs; ii++ )
863  {
864  QString mlines = mcixmls[ ii ];
865 
866  // Limit contents to <model>...</model> except for first and last
867  int flx = ( ii == 0 ) ? 0 : 3;
868  int llx = ( ( ii + 1 ) < nmcixs ) ? -3 : -2;
869  QString mcont = mlines.section( "\n", flx, llx ) + "\n";
870  // Concatenate to output stream
871  tso << mcont;
872  }
873 }
874 
875 // Read scan C0 values from an XML stream
876 void US_Model::mfem_scans( QXmlStreamReader& xml, SimulationComponent& sc )
877 {
878  while ( ! xml.atEnd() )
879  {
880  xml.readNext();
881 
882  if ( xml.isStartElement() )
883  {
884  if ( xml.name() == "mfem_scan" )
885  {
886  QXmlStreamAttributes a = xml.attributes();
887 
888  sc.c0.radius << a.value( "radius" ).toString().toDouble();
889  sc.c0.concentration << a.value( "conc" ).toString().toDouble();
890  }
891  else
892  break;
893  }
894  }
895 }
896 
897 // Get associations from an XML stream
898 void US_Model::get_associations( QXmlStreamReader& xml, Association& as )
899 {
900  QXmlStreamAttributes a = xml.attributes();
901  QString skassoc = a.value( "k_assoc" ).toString();
902  QString skeq = a.value( "k_eq" ).toString();
903  QString skdisso = a.value( "K_d" ).toString();
904  QString skoff = a.value( "k_off" ).toString();
905 
906  if ( ! skdisso.isEmpty() )
907  as.k_d = skdisso.toDouble();
908  else
909  {
910  double k_assoc = 0.0;
911  if ( ! skassoc.isEmpty() )
912  k_assoc = skassoc.toDouble();
913  else if ( ! skeq.isEmpty() )
914  k_assoc = skeq.toDouble();
915  as.k_d = ( k_assoc != 0.0 ) ? ( 1.0 / k_assoc ) : 0.0;
916  }
917 
918  if ( ! skoff.isEmpty() )
919  as.k_off = skoff.toDouble();
920 
921  while ( ! xml.atEnd() )
922  {
923  xml.readNext();
924 
925  if ( xml.isStartElement() )
926  {
927  if ( xml.name() == "component" )
928  {
929  a = xml.attributes();
930 
931  as.rcomps << a.value( "index" ).toString().toInt();
932  as.stoichs << a.value( "stoich" ).toString().toInt();
933  }
934  else
935  break;
936  }
937  }
938 }
939 
940 // Load a model from DB (by GUID)
941 int US_Model::load_db( const QString& guid, US_DB2* db )
942 {
943  QStringList q;
944 
945  q << "get_modelID" << guid;
946  db->query( q );
947 
948  if ( db->lastErrno() != US_DB2::OK ) return db->lastErrno();
949 
950  db->next();
951  QString id = db->value( 0 ).toString();
952  return load( id, db );
953 }
954 
955 // Load a model from DB (by DB id)
956 int US_Model::load( const QString& id, US_DB2* db )
957 {
958  QStringList q;
959 
960  q << "get_model_info" << id;
961  db->query( q );
962 
963  if ( db->lastErrno() != US_DB2::OK ) return db->lastErrno();
964 
965  db->next();
966  QByteArray contents = db->value( 2 ).toString().toAscii();
967 
968  // Read the model file into an array in memory
969  QXmlStreamReader xml( contents );
970 
971  int result = load_stream( xml );
972 
973  if ( result == US_DB2::NO_MODEL && monteCarlo )
974  { // Handle a multi-model stream
975  QTextStream tsi( contents );
976 
977  result = load_multi_model( tsi );
978  }
979 
980  return result;
981 }
982 
983 // Write a model to DB or local file
984 int US_Model::write( bool db_access, const QString& filename, US_DB2* db )
985 {
986  if ( db_access ) return write( db );
987  else return write( filename );
988 }
989 
990 // Write a model DB record
991 int US_Model::write( US_DB2* db )
992 {
993  // Create the model xml file in a stream
994  QByteArray temporary;
995  QByteArray contents;
996 
997  if ( ! monteCarlo || nmcixs < 1 )
998  {
999  QXmlStreamWriter xml( &temporary );
1000  write_stream( xml );
1001  }
1002 
1003  else
1004  {
1005  QTextStream tso( &temporary );
1006  write_mm_stream( tso );
1007  }
1008 
1009  db->mysqlEscapeString( contents, temporary, temporary.size() );
1010 qDebug() << "model writedb contsize tempsize" << contents.size() << temporary.size();
1011 
1012  QStringList q;
1013 
1014  // Generate a guid if necessary
1015  // The guid may be valid from a disk read, but is not in the DB
1016  if ( modelGUID.size() != 36 ) modelGUID = US_Util::new_guid();
1017 
1018  q << "get_modelID" << modelGUID;
1019 
1020  db->query( q );
1021 
1022  QString meni = QString::number( meniscus );
1023  QString vari = QString::number( variance );
1024 
1025  if ( db->lastErrno() != US_DB2::OK )
1026  {
1027  q.clear();
1028  q << "new_model" << modelGUID << description << contents
1029  << vari << meni << editGUID
1030  << QString::number( US_Settings::us_inv_ID() );
1031  message = QObject::tr( "created" );
1032  }
1033  else
1034  {
1035  db->next();
1036  QString id = db->value( 0 ).toString();
1037  q.clear();
1038  q << "update_model" << id << description << contents
1039  << vari << meni << editGUID;
1040  message = QObject::tr( "updated" );
1041  }
1042 
1043  int wstat = db->statusQuery( q );
1044 qDebug() << "model writedb message" << message << "wstat" << wstat;
1045  QString path;
1046  model_path( path );
1047  bool newFile;
1048  QString filename = get_filename( path, modelGUID, newFile );
1049  write( filename );
1050 
1051  return wstat;
1052 }
1053 
1054 // Write a model file
1055 int US_Model::write( const QString& filename )
1056 {
1057  QFile file( filename );
1058 
1059  if ( ! file.open( QIODevice::WriteOnly | QIODevice::Text) )
1060  return US_DB2::ERROR;
1061 
1062  if ( ! monteCarlo || nmcixs < 1 )
1063  {
1064  QXmlStreamWriter xml( &file );
1065  write_stream( xml );
1066  }
1067 
1068  else
1069  {
1070  QTextStream tso( &file );
1071  write_mm_stream( tso );
1072  }
1073 
1074  file.close();
1075 
1076  return US_DB2::OK;
1077 }
1078 
1079 // Write to an XML stream
1080 void US_Model::write_stream( QXmlStreamWriter& xml )
1081 {
1082  if ( modelGUID.size() != 36 )
1083  modelGUID = US_Util::new_guid();
1084 
1085  xml.setAutoFormatting( true );
1086 
1087  xml.writeStartDocument();
1088  xml.writeDTD ( "<!DOCTYPE US_Model>" );
1089  xml.writeStartElement( "ModelData" );
1090  xml.writeAttribute ( "version", "1.0" );
1091 
1092  xml.writeStartElement( "model" );
1093  xml.writeAttribute ( "description", description );
1094  xml.writeAttribute ( "modelGUID", modelGUID );
1095  xml.writeAttribute ( "editGUID", editGUID );
1096  xml.writeAttribute ( "wavelength", QString::number( wavelength ) );
1097  if ( variance != 0.0 )
1098  xml.writeAttribute( "variance", QString::number( variance ) );
1099  if ( meniscus != 0.0 )
1100  xml.writeAttribute( "meniscus", QString::number( meniscus ) );
1101  if ( alphaRP != 0.0 )
1102  xml.writeAttribute( "alphaRP", QString::number( alphaRP ) );
1103  xml.writeAttribute ( "coSedSolute", QString::number( coSedSolute ) );
1104  if ( subGrids != 0.0 )
1105  xml.writeAttribute( "subGrids", QString::number( subGrids ) );
1106  xml.writeAttribute ( "opticsType", QString::number( optics ) );
1107  xml.writeAttribute ( "analysisType",QString::number( analysis ) );
1108  xml.writeAttribute ( "globalType", QString::number( global ) );
1109 
1110  if ( requestGUID.length() > 0 )
1111  xml.writeAttribute ( "requestGUID", requestGUID );
1112 
1113  if ( monteCarlo )
1114  xml.writeAttribute ( "monteCarlo", "1" );
1115 
1116  if ( ! dataDescrip.isEmpty() )
1117  {
1118  dataDescrip.replace( '"', "&quot;" ); // Replace '"' character
1119  dataDescrip.replace( '<', "&lt;" ); // Replace '<' character
1120  dataDescrip.replace( '>', "&gt;" ); // Replace '>' character
1121  xml.writeAttribute ( "dataDescrip", dataDescrip );
1122  }
1123 
1124  // Write components
1125  int ncomps = components.size();
1126  bool notmany = ( ncomps < 400 );
1127 
1128  for ( int i = 0; i < ncomps; i++ )
1129  {
1130  SimulationComponent* sc = &components[ i ];
1131  xml.writeStartElement( "analyte" );
1132 
1133  if ( !sc->analyteGUID.isEmpty() )
1134  xml.writeAttribute( "analyteGUID", sc->analyteGUID );
1135  xml.writeAttribute( "name", sc->name );
1136  xml.writeAttribute( "mw", QString::number( sc->mw ) );
1137  xml.writeAttribute( "s", QString::number( sc->s ) );
1138  xml.writeAttribute( "D", QString::number( sc->D ) );
1139  xml.writeAttribute( "f", QString::number( sc->f ) );
1140  xml.writeAttribute( "f_f0", QString::number( sc->f_f0 ) );
1141  QString strVbar = QString::number( sc->vbar20 );
1142  QString strExtinc = QString::number( sc->extinction );
1143  QString strAxial = QString::number( sc->axial_ratio );
1144  QString strSigma = QString::number( sc->sigma );
1145  QString strDelta = QString::number( sc->delta );
1146  QString strOligo = QString::number( sc->oligomer );
1147  QString strShape = QString::number( sc->shape );
1148  QString strType = QString::number( sc->analyte_type );
1149  QString strMolar = QString::number( sc->molar_concentration );
1150  QString strSignal = QString::number( sc->signal_concentration );
1151  if ( notmany )
1152  { // In most cases, don't test values to write
1153  xml.writeAttribute( "vbar20", strVbar );
1154  xml.writeAttribute( "extinction", strExtinc );
1155  xml.writeAttribute( "axial", strAxial );
1156  xml.writeAttribute( "sigma", strSigma );
1157  xml.writeAttribute( "delta", strDelta );
1158  xml.writeAttribute( "oligomer", strOligo );
1159  xml.writeAttribute( "shape", strShape );
1160  xml.writeAttribute( "type", strType );
1161  xml.writeAttribute( "molar", strMolar );
1162  }
1163  else
1164  { // If many components, only write non-default values
1165  if ( sc->vbar20 != TYPICAL_VBAR )
1166  xml.writeAttribute( "vbar20", strVbar );
1167  if ( sc->extinction != 0.0 )
1168  xml.writeAttribute( "extinction", strExtinc );
1169  if ( sc->axial_ratio != 10.0 )
1170  xml.writeAttribute( "axial", strAxial );
1171  if ( sc->sigma != 0.0 )
1172  xml.writeAttribute( "sigma", strSigma );
1173  if ( sc->delta != 0.0 )
1174  xml.writeAttribute( "delta", strDelta );
1175  if ( sc->oligomer != 1 )
1176  xml.writeAttribute( "oligomer", strOligo );
1177  if ( sc->shape != SPHERE )
1178  xml.writeAttribute( "shape", strShape );
1179  if ( sc->analyte_type != 0 )
1180  xml.writeAttribute( "type", strType );
1181  if ( sc->molar_concentration != 0.0 )
1182  xml.writeAttribute( "molar", strMolar );
1183  }
1184  xml.writeAttribute( "signal", strSignal );
1185 
1186  for ( int j = 0; j < sc->c0.radius.size(); j++ )
1187  {
1188  xml.writeStartElement( "mfem_scan" );
1189 
1190  MfemInitial* scan = &sc->c0;
1191  xml.writeAttribute( "radius",
1192  QString::number( scan->radius [ j ] ) );
1193  xml.writeAttribute( "conc",
1194  QString::number( scan->concentration[ j ] ) );
1195  xml.writeEndElement(); // mfem_scan
1196  }
1197 
1198  xml.writeEndElement(); // analyte (SimulationComponent)
1199  }
1200 
1201  // Write associations
1202  for ( int i = 0; i < associations.size(); i++ )
1203  {
1204  Association* as = &associations[ i ];
1205  xml.writeStartElement( "association" );
1206  xml.writeAttribute( "K_d", QString::number( as->k_d ) );
1207  xml.writeAttribute( "k_off", QString::number( as->k_off ) );
1208 
1209  for ( int j = 0; j < as->rcomps.size(); j++ )
1210  {
1211  xml.writeStartElement( "component" );
1212 
1213  QString index = QString::number( as->rcomps [ j ] );
1214  QString stoich = QString::number( as->stoichs[ j ] );
1215 
1216  xml.writeAttribute( "index", index );
1217  xml.writeAttribute( "stoich", stoich );
1218  xml.writeEndElement(); // component
1219  }
1220 
1221  xml.writeEndElement(); // association
1222  }
1223 
1224  xml.writeEndElement(); // model
1225  xml.writeEndElement(); // ModelData
1226  xml.writeEndDocument();
1227 }
1228 
1229 // Get the name of a model file with matching GUID (if any)
1230 QString US_Model::get_filename( const QString& path, const QString& guid,
1231  bool& newFile )
1232 {
1233  QDir f( path );
1234  QStringList filter( "M???????.xml" );
1235  QStringList f_names = f.entryList( filter, QDir::Files, QDir::Name );
1236  QString fnamo;
1237  newFile = true;
1238 
1239  for ( int i = 0; i < f_names.size(); i++ )
1240  {
1241  QFile m_file( path + "/" + f_names[ i ] );
1242 
1243  if ( ! m_file.open( QIODevice::ReadOnly | QIODevice::Text) ) continue;
1244 
1245  QXmlStreamReader xml( &m_file );
1246 
1247  while ( ! xml.atEnd() )
1248  {
1249  xml.readNext();
1250 
1251  if ( xml.isStartElement() )
1252  {
1253  if ( xml.name() == "model" )
1254  {
1255  QXmlStreamAttributes a = xml.attributes();
1256 
1257  if ( a.value( "modelGUID" ).toString() == guid )
1258  { // Break when we have found a file with a matching GUID
1259  fnamo = path + "/" + f_names[ i ];
1260  newFile = false;
1261  break;
1262  }
1263  }
1264  }
1265  }
1266 
1267  m_file.close();
1268  }
1269 
1270  if ( newFile )
1271  { // No matching-GUID file, so look for a break in the number sequence
1272  int number = ( f_names.size() > 0 ) ? // Last used file number
1273  f_names.last().mid( 1, 7 ).toInt() : 0; // and default last sequence
1274 
1275  for ( int ii = 0; ii < number; ii++ )
1276  {
1277  QString fnamck = "M" + QString().sprintf( "%07i", ii + 1 ) + ".xml";
1278 
1279  if ( ! f_names.contains( fnamck ) )
1280  { // There is a hole in the sequence, so re-use this number
1281  number = ii;
1282  break;
1283  }
1284  }
1285 
1286  // File name uses a break in the sequence or one past last used.
1287  fnamo = path + "/M" + QString().sprintf( "%07i", number + 1 ) + ".xml";
1288  }
1289 
1290  return fnamo;
1291 }
1292 
1293 // Create or append to a composite MC model file for a single triple
1294 QString US_Model::composite_mc_file( QStringList& mcfiles, const bool rmvi )
1295 {
1296  const QChar dquo( '"' );
1297  QString empty_str( "" );
1298  QString cmfname = empty_str;
1299  int mc_iters = mcfiles.size();
1300  int mc_comps = 0;
1301 
1302  // Return an empty name if the list is empty
1303  if ( mc_iters < 1 ) return cmfname;
1304 
1305  cmfname = mcfiles[ 0 ];
1306  bool name_desc = cmfname.contains( ".mc" ) || cmfname.contains( "_mc" );
1307 //qDebug() << "MDL:CMF: name_desc" << name_desc;
1308 
1309  // Return the name of an existing composite if it is all of the list
1310  if ( mc_iters == 1 )
1311  {
1312  if ( cmfname.contains( ".mcN" ) || cmfname.contains( "_mcN" ) )
1313  return cmfname;
1314  if ( ! name_desc )
1315  { // If no ".mc" in file name, must check description in contents
1316  QFile filei( cmfname );
1317  if ( ! filei.open( QIODevice::ReadOnly | QIODevice::Text) )
1318  {
1319  qDebug() << "**MC iteration file open error**";
1320  return empty_str;
1321  }
1322 
1323  QTextStream tsi( &filei );
1324  QString mcont = tsi.readAll();
1325  filei.close();
1326  int jj = qMax( 0, mcont.indexOf( "description=" ) );
1327  QString mdesc = QString( mcont ).mid( jj, 99 ).section( dquo, 1, 1 );
1328  if ( mdesc.contains( "_mcN" ) )
1329  return cmfname;
1330  }
1331  }
1332 
1333  // Otherwise, scan the MC iteration file names
1334  for ( int ii = 0; ii < mcfiles.size(); ii++ )
1335  {
1336  if ( name_desc )
1337  { // Handle a file whose name tells that it is a composite
1338  if ( mcfiles[ ii ].contains( ".mcN" ) ||
1339  mcfiles[ ii ].contains( "_mcN" ) )
1340  { // For composite, bump composite count, decrement iters, save name
1341  mc_comps++;
1342  mc_iters--;
1343  cmfname = mcfiles[ ii ];
1344  }
1345  }
1346  else
1347  { // Handle a file whose description tells that it is a composite
1348  QFile filei( mcfiles[ ii ] );
1349  if ( ! filei.open( QIODevice::ReadOnly | QIODevice::Text) )
1350  continue;
1351  QTextStream tsi( &filei );
1352  QString mcont = tsi.readAll();
1353  filei.close();
1354  int jj = qMax( 0, mcont.indexOf( "description=" ) );
1355  QString mdesc = QString( mcont ).mid( jj, 99 ).section( dquo, 1, 1 );
1356  if ( mdesc.contains( "_mcN" ) )
1357  { // For composite, bump composite count, decrement iters, save name
1358  mc_comps++;
1359  mc_iters--;
1360  cmfname = mcfiles[ ii ];
1361 //qDebug() << "MDL:CMF: (A)cmfname" << cmfname;
1362  }
1363  }
1364  }
1365 
1366  if ( mc_comps > 1 || mc_iters < 1 )
1367  { // Return now if we don't have 0 or 1 composite, plus some iter models
1368  qDebug() << "**" << mc_comps << "MC composites, and" << mc_iters
1369  << "MC iterations **";
1370  return empty_str;
1371  }
1372 
1373  QString mditer1 = QString().sprintf( ".mcN%03i", mc_iters );
1374  QString mditer2 = QString( mditer1 ).replace( ".mcN", "_mcN" );
1375  QString mditer = mditer1;
1376 //qDebug() << "MDL:CMF: mditer1 mditer2" << mditer1 << mditer2;
1377 
1378  if ( mc_comps == 0 )
1379  { // No composite exists, so create one from the iteration models
1380  cmfname = mcfiles[ 0 ];
1381  QFile filei( mcfiles[ 0 ] );
1382  if ( ! filei.open( QIODevice::ReadOnly | QIODevice::Text) )
1383  {
1384  qDebug() << "**MC iteration file open error**";
1385  return empty_str;
1386  }
1387  QTextStream tsi( &filei );
1388  QString mcont = tsi.readAll();
1389  filei.close();
1390  int jj = qMax( 0, mcont.indexOf( "modelGUID=" ) );
1391  QString mcguid = QString( mcont ).mid( jj, 99 ).section( dquo, 1, 1 );
1392 
1393  if ( ! name_desc )
1394  { // No description in name ("M00...xml"), so determine composite name
1395  if ( rmvi )
1396  { // Removing iteration files, so reuse first name for composite
1397  filei.remove();
1398  }
1399 
1400  else
1401  { // If not removing iteration files, get a new name for composite
1402  bool newFile = true;
1403  QString path = QString( cmfname ).section( "/", 0, -2 );
1404  mcguid = US_Util::new_guid();
1405  cmfname = get_filename( path, mcguid, newFile );
1406 //qDebug() << "MDL:CMF: (B)cmfname" << cmfname;
1407  }
1408  }
1409 
1410  else
1411  { // Description in name, so create a new name with iters count in it
1412  int moix1 = cmfname.indexOf( ".mc" );
1413  int moix2 = cmfname.indexOf( "_mc" );
1414  int moix = ( moix1 > 0 ) ? moix1 : qMax( 0, moix2 );
1415  QString moiter = QString( cmfname ).mid( moix, 7 );
1416  mditer = ( moix1 > 0 ) ? mditer1 : mditer2;
1417  cmfname = cmfname.replace( moiter, mditer )
1418  .replace( ".mdl.tmp", ".model.xml" );
1419 //qDebug() << "MDL:CMF: moix1 moix2 moiter mditer" << moix1 << moix2 << moiter << mditer;
1420 //qDebug() << "MDL:CMF: (C)cmfname" << cmfname;
1421  }
1422 
1423  QFile fileo( cmfname );
1424  if ( ! fileo.open( QIODevice::WriteOnly | QIODevice::Text ) )
1425  {
1426  qDebug() << "**MC composite file open error**";
1427  return empty_str;
1428  }
1429 
1430  // Output contents of first iteration model except last line
1431  QTextStream tso( &fileo );
1432  jj = qMax( 0, mcont.indexOf( "modelGUID=" ) );
1433  QString miguid = QString( mcont ).mid( jj, 99 ).section( dquo, 1, 1 );
1434  if ( miguid != mcguid )
1435  mcont.replace( miguid, mcguid );
1436  int flx = 0;
1437  int llx = ( mc_iters > 1 ) ? -3 : -2;
1438  tso << mcont.section( "\n", flx, llx ) << "\n";
1439  flx = 3;
1440 
1441  // Output <model>...</model> for 2nd thru last iteration model
1442  for ( int ii = 1; ii < mc_iters; ii++ )
1443  {
1444  QFile filei( mcfiles[ ii ] );
1445  if ( ! filei.open( QIODevice::ReadOnly | QIODevice::Text) )
1446  continue;
1447  QTextStream tsi( &filei );
1448  mcont = tsi.readAll();
1449  filei.close();
1450  jj = qMax( 0, mcont.indexOf( "modelGUID=" ) );
1451  QString miguid = QString( mcont ).mid( jj, 99 ).section( dquo, 1, 1 );
1452  if ( miguid != mcguid )
1453  mcont.replace( miguid, mcguid );
1454  llx = ( ( ii + 1 ) < mc_iters ) ? -3 : -2;
1455  tso << mcont.section( "\n", flx, llx ) << "\n";
1456  }
1457 
1458  fileo.close();
1459  }
1460 
1461  else
1462  { // A composite exists, so append to it from new iteration models
1463  QFile filei( cmfname );
1464  if ( ! filei.open( QIODevice::ReadOnly | QIODevice::Text) )
1465  {
1466  qDebug() << "**MC composite file open error**";
1467  return empty_str;
1468  }
1469  // Read in contents of existing composite model; skip last line
1470  QTextStream tsi( &filei );
1471  QString mcont = tsi.readAll();
1472  filei.close();
1473  mcont = mcont.section( "\n", 0, -3 ) + "\n";
1474  int jj = qMax( 0, mcont.indexOf( "modelGUID=" ) );
1475  QString mcguid = QString( mcont ).mid( jj, 99 ).section( dquo, 1, 1 );
1476 
1477  if ( name_desc )
1478  { // Name has description, so see if it needs to be renamed
1479  int moix1 = cmfname.indexOf( ".mc" );
1480  int moix2 = cmfname.indexOf( "_mc" );
1481  int moix = ( moix1 > 0 ) ? moix1 : qMax( 0, moix2 );
1482  QString moiter = QString( cmfname ).mid( moix, 7 );
1483  int mc_ittot = QString( moiter ).mid( 4 ).toInt() + mc_iters;
1484  mditer = ( ( moix1 > 0 ) ? "." : "_" )
1485  + QString().sprintf( "mcN%03i", mc_ittot );
1486 //qDebug() << "MDL:CMF: (B)moix1 moix2 moiter mditer" << moix1 << moix2 << moiter << mditer;
1487 
1488  if ( moiter != mditer )
1489  { // Iterations changes (almost always), so rename is necessary
1490  cmfname = cmfname.replace( moiter, mditer )
1491  .replace( ".mdl.tmp", ".model.xml" );
1492  filei.rename( cmfname );
1493 //qDebug() << "MDL:CMF: (D)cmfname" << cmfname;
1494  }
1495  else
1496  { // No name change (unlikely), so must delete and re-create
1497  filei.remove();
1498  }
1499  }
1500 
1501  else
1502  { // Name has no description, so no renaming necessary and just recreate
1503  filei.remove();
1504  }
1505 
1506  QFile fileo( cmfname );
1507  if ( ! fileo.open( QIODevice::WriteOnly | QIODevice::Text ) )
1508  {
1509  qDebug() << "**MC composite file open error**";
1510  return empty_str;
1511  }
1512 
1513  // Output contents of previous composite except for last line
1514  QTextStream tso( &fileo );
1515  int flx = 0;
1516  int llx = -3;
1517  tso << mcont.section( "\n", flx, llx ) << "\n";
1518  flx = 3;
1519 
1520  // Append <model>...</model> of all new iteration models
1521  for ( int ii = 0; ii < mc_iters; ii++ )
1522  {
1523  QFile filei( mcfiles[ ii ] );
1524  if ( ! filei.open( QIODevice::ReadOnly | QIODevice::Text) )
1525  continue;
1526  QTextStream tsi( &filei );
1527  mcont = tsi.readAll();
1528  filei.close();
1529  jj = qMax( 0, mcont.indexOf( "modelGUID=" ) );
1530  QString miguid = QString( mcont ).mid( jj, 99 ).section( dquo, 1, 1 );
1531  if ( miguid != mcguid )
1532  mcont.replace( miguid, mcguid );
1533  llx = ( ( ii + 1 ) < mc_iters ) ? -3 : -2;
1534  tso << mcont.section( "\n", flx, llx ) << "\n";
1535  }
1536 
1537  fileo.close();
1538  }
1539 
1540  if ( rmvi )
1541  { // If so requested, remove the MC iteration files
1542  if ( mc_comps == 0 )
1543  { // If list is all iteration files, remove them all
1544  for ( int ii = 0; ii < mcfiles.size(); ii++ )
1545  { // Delete unless file name was reused for composite
1546  if ( mcfiles[ 0 ] != cmfname )
1547  QFile( mcfiles[ ii ] ).remove();
1548  }
1549  }
1550 
1551  else if ( name_desc )
1552  { // If list has composite, name has description; selectively delete
1553  for ( int ii = 0; ii < mcfiles.size(); ii++ )
1554  { // Delete iteration files (not composite)
1555  if ( ! mcfiles[ ii ].contains( "mcN" ) )
1556  QFile( mcfiles[ ii ] ).remove();
1557  }
1558  }
1559 
1560  else
1561  { // If no description in file names, delete based on descr. content
1562  for ( int ii = 0; ii < mcfiles.size(); ii++ )
1563  { // Delete iteration files (not composite)
1564  if ( mcfiles[ ii ] == cmfname )
1565  continue;
1566  QFile filei( mcfiles[ ii ] );
1567  if ( ! filei.open( QIODevice::ReadOnly | QIODevice::Text) )
1568  {
1569  qDebug() << "**MC iteration file open error**";
1570  continue;
1571  }
1572 
1573  QTextStream tsi( &filei );
1574  QString mcont = tsi.readAll();
1575  filei.close();
1576  int jj = qMax( 0, mcont.indexOf( "description=" ) );
1577  QString mdesc = QString(mcont).mid( jj, 99 ).section( dquo, 1, 1 );
1578  if ( ! mdesc.contains( "_mcN" ) )
1579  QFile( mcfiles[ ii ] ).remove();
1580  }
1581  }
1582  }
1583 
1584  return cmfname;
1585 }
1586 
1587 // Get MC iteration xml content strings
1588 int US_Model::mc_iter_xmls( QStringList& mcixs )
1589 {
1590  int kixmls = mcixmls.size();
1591  mcixs = mcixmls;
1592  return kixmls;
1593 }
1594 
1595 // Output debug print model details
1596 void US_Model::debug( void )
1597 {
1598  qDebug() << "model dump";
1599  qDebug() << "desc" << description;
1600  qDebug() << "model guid" << modelGUID;
1601  qDebug() << "variance" << variance;
1602  qDebug() << "meniscus" << meniscus;
1603  qDebug() << "alphaRP" << alphaRP;
1604  qDebug() << "edit guid" << editGUID;
1605  qDebug() << "request guid" << requestGUID;
1606  qDebug() << "waveln" << wavelength;
1607  qDebug() << "monte carlo" << monteCarlo;
1608  qDebug() << "coSed" << coSedSolute;
1609  qDebug() << "subGrids" << subGrids;
1610  qDebug() << "AnalysisType" << (int)analysis;
1611  qDebug() << "GlobalType" << (int)global;
1612  qDebug() << "OpticsType" << (int)optics;
1613  qDebug() << "data description" << dataDescrip;
1614 
1615  for ( int i = 0; i < components.size(); i++ )
1616  {
1617  SimulationComponent* sc = &components[ i ];
1618  qDebug() << " component" << ( i + 1 );
1619  qDebug() << " name" << sc->name;
1620  qDebug() << " vbar20" << sc->vbar20;
1621  qDebug() << " mw" << sc->mw;
1622  qDebug() << " s" << sc->s;
1623  qDebug() << " D" << sc->D;
1624  qDebug() << " f" << sc->f;
1625  qDebug() << " f/f0" << sc->f_f0;
1626  qDebug() << " extinction" << sc->extinction;
1627  qDebug() << " axial" << sc->axial_ratio;
1628  qDebug() << " sigma" << sc->sigma;
1629  qDebug() << " delta" << sc->delta;
1630  qDebug() << " molar C" << sc->molar_concentration;
1631  qDebug() << " signal C" << sc->signal_concentration;
1632  qDebug() << " oligomer" << sc->oligomer;
1633  qDebug() << " shape" << (int)sc->shape;
1634  qDebug() << " type" << (int)sc->analyte_type;
1635 
1636  for ( int j = 0; j < sc->c0.radius.size(); j++ )
1637  {
1638  qDebug() << " c0 r c"
1639  << sc->c0.radius[ j ] << sc->c0.concentration[ j ];
1640  }
1641  }
1642 }
1643 
Generated on Sun Jun 21 2015 11:40:42 for UltraScan III by   doxygen 1.8.3.1-20130324