Changeset 2808


Ignore:
Timestamp:
Aug 13, 2019 8:11:33 PM (3 months ago)
Author:
gegorbet
Message:

density_match initial computation code complete and all plots working

Location:
trunk/programs/us_density_match
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/programs/us_density_match/us_density_match.cpp

    r2807 r2808  
    4848US_Density_Match::US_Density_Match() : US_Widgets()
    4949{
     50   dbg_level = US_Settings::us_debug();
     51
    5052   // Set up the GUI
    51 
    5253   setWindowTitle( tr( "Density Matching" ) );
    5354   setPalette( US_GuiSettings::frameColor() );
    54 
    5555
    5656   // Primary layouts
     
    6565   spec->setContentsMargins( 0, 0, 0, 0 );
    6666
    67    int s_row = 0;
    68    dbg_level = US_Settings::us_debug();
    69 DbgLv(1) << "MD: main: AA";
    70 //   clean_etc_dir();
    71 
    7267   // Top banner
    7368   QLabel* lb_info1      = us_banner( tr( "Model Selection Controls" ) );
    7469
    75 
    7670   us_checkbox( tr( "Save Plot(s)"    ), ck_savepl,  false );
    7771   us_checkbox( tr( "Local Save Only" ), ck_locsave, true  );
    7872
     73   // Distribution information text box
    7974   te_distr_info = us_textedit();
    80    te_distr_info->setText    ( tr( "Run:  runID.triple (method)\n" )
     75   te_distr_info->setText( tr( "Run:  runID.triple (method)\n" )
    8176            + tr( "    analysisID" ) );
    8277   us_setReadOnly( te_distr_info, true );
    8378
     79   // X axis radio buttons and button group
    8480   plot_x      = 0;
    85 
    8681   QLabel* lb_x_axis   = us_label( tr( "Plot X:" ) );
    8782           bg_x_axis   = new QButtonGroup( this );
     
    9893   bg_x_axis->addButton( rb_x_s,    ATTR_S );
    9994   bg_x_axis->addButton( rb_x_d,    ATTR_D );
    100    rb_x_mass->setChecked( true  );
    10195   rb_x_mass->setToolTip( tr( "Set X axis to Molar Mass"                ) );
    10296   rb_x_ff0 ->setToolTip( tr( "Set X axis to Frictional Ratio"          ) );
     
    10599   rb_x_s   ->setToolTip( tr( "Set X axis to Sedimentation Coefficient" ) );
    106100   rb_x_d   ->setToolTip( tr( "Set X axis to Diffusion Coefficient"     ) );
     101   rb_x_s   ->setChecked( true );
    107102   connect( bg_x_axis,  SIGNAL( buttonReleased( int ) ),
    108103            this,       SLOT  ( select_x_axis ( int ) ) );
     
    160155
    161156   // Order plot components on the left side
     157   int s_row = 0;
    162158   spec->addWidget( lb_info1,      s_row++, 0, 1, 8 );
    163159   spec->addLayout( dkdb_cntrls,   s_row++, 0, 1, 8 );
     
    258254}
    259255
     256// Reset information: clear prefilter and distributions
    260257void US_Density_Match::reset( void )
    261258{
     
    275272}
    276273
     274// Save plots and CSV files
    277275void US_Density_Match::save( void )
    278276{
     
    335333}
    336334
    337 // plot the data
     335// Plot the data
    338336void US_Density_Match::plot_data( void )
    339337{
     
    345343   DisSys* tsys   = (DisSys*)&alldis.at( 0 );
    346344   plot_x         = ( plot_x < 0 ) ? plot_x_select() : plot_x;
    347 
    348345
    349346#if 0
     
    378375   double* xx     = v_sedcs[ 0 ].data();
    379376   QString curvtitl( "curve" );
     377   // Plot title
    380378   QString tstr   = QString( tsys->run_name ).section( ".", 0, -2 );
    381379DbgLv(1) << "DaPl: (1)tstr" << tstr;
     
    390388DbgLv(1) << "DaPl: (2)tstr" << tstr;
    391389
     390   // X pointer, curve title, rest of plot title; based on plot type
    392391   if      ( plot_x == ATTR_S )
    393392   {
     
    408407DbgLv(1) << "DaPl: v_vbars" << v_vbars;
    409408      curvtitl       = tr( "Vbar_curve" );
    410       tstr          += tr( "Specific Density" );
     409      tstr          += tr( "Partial Specific Density" );
    411410   }
    412411   else if ( plot_x == ATTR_D )
     
    431430DbgLv(1) << "DaPl: (3)tstr" << tstr;
    432431
     432   // Initial plot settings
    433433   dataPlotClear( data_plot );
    434434   data_plot->replot();
     
    452452   {
    453453      // Plot a single line
     454      data_plot->detachItems( QwtPlotItem::Rtti_PlotLegend, true );
     455      QwtLegend *legend = new QwtLegend;
     456      data_plot->insertLegend( legend, QwtPlot::BottomLegend  );
    454457      QColor colr1( Qt::blue );
    455458      data_curv      = us_curve( data_plot, curvtitl );
     459      data_curv->setItemAttribute( QwtPlotItem::Legend, true );
    456460      data_curv->setPen  ( QPen( QBrush( colr1 ), 3.0, Qt::SolidLine ) );
    457461      data_curv->setStyle( QwtPlotCurve::Lines );
     
    496500}
    497501
     502// Plot data based on current plot type index
    498503void US_Density_Match::plot_data( int cplx )
    499504{
     
    564569}
    565570
     571// Load all the distributions (models)
    566572void US_Density_Match::load_distro()
    567573{
     
    583589   need_save  = false;
    584590
     591   te_distr_info->setText(
     592      QString( models[ 0 ].description ).section( ".", 0, -4 ) );
     593
    585594   for ( int jj = 0; jj < models.count(); jj++ )
    586595   {  // Load each selected distribution model
     
    593602
    594603   // Notify user of need to set D2O-percent, label, density model values
    595    QString qmsg = tr( "%1 models are loaded. In the dialog to follow,\n"
     604   QString qmsg = tr( "%1 models are loaded.\n\nIn the dialog to follow,\n"
    596605                      "you must set D2O Percent values for each,\n"
    597606                      "then review and set Label and Density values for them.\n"
    598                       "Begin with the model(s) with 0% D2O." )
     607                      "Begin with the model with 0% D2O." )
    599608                      .arg( alldis.size() );
    600609   QMessageBox::warning( this, tr( "Model Parameters" ), qmsg );
     
    602611   // Set model distributions parameters
    603612   set_mparms();
     613
     614   // Plot S lines for all models
     615   rb_x_s->setChecked( true );
     616   select_x_axis( ATTR_S );
    604617}
    605618
     
    631644   }
    632645   tsys.editGUID     = model.editGUID;
    633    tsys.plot_x       = plot_x;
    634646DbgLv(1) << "LD: method" << tsys.method << "mdesc" << mdesc;
    635647
     
    703715
    704716   
    705    te_distr_info->setText( tr( "Run:  " ) + tsys.run_name
    706       + " (" + tsys.method + ")\n    " + tsys.analys_name );
     717//   te_distr_info->setText( tr( "Run:  " ) + tsys.run_name
     718//      + " (" + tsys.method + ")\n    " + tsys.analys_name );
    707719   int nsolmc  = model.components.size();
    708720
     
    738750   for ( int jj = 0; jj < nsolnm; jj++ )
    739751   {
    740       sol_nm      = wk_distro[ jj ];
    741       sol_nm.c   /= tot_conc;
    742       tsys.nm_distro << sol_nm;
     752      sol_nm      = wk_distro[ jj ];  // Solute point
     753      sol_nm.c   /= tot_conc;         // Normalized concentration
     754      tsys.nm_distro << sol_nm;       // Saved to "normalized" distro
    743755if ( jj<3   ||  (jj+4)>nsolnm )
    744756 DbgLv(1) << "LD:    jj" << jj << "soln s,d,c"
     
    749761   // Create version of distribution with boundary fraction
    750762
    751    tsys.bo_distro.clear();
     763   tsys.bo_distro.clear();                // Boundary distro
    752764   double sum_co  = 0.0;
    753765
    754766   for ( int jj = 0; jj < nsolnm; jj++ )
    755767   {
    756       sol_nm      = tsys.nm_distro[ jj ];
    757       sum_co     += sol_nm.c;
    758       sol_bf      = sol_nm;
    759       sol_bf.f    = sum_co;
    760       tsys.bo_distro << sol_bf;
     768      sol_nm      = tsys.nm_distro[ jj ]; // Norm'd solute point
     769      sum_co     += sol_nm.c;             // Concentration integration
     770      sol_bf      = sol_nm;               // Boundary fraction solute point
     771      sol_bf.f    = sum_co;               // With boundary fraction for "f"
     772      tsys.bo_distro << sol_bf;           // Save to boundary distro
    761773if ( jj<3   ||  (jj+4)>nsolnm )
    762774 DbgLv(1) << "LD:    jj" << jj << "solb s,d,c,f"
     
    787799      pb_pltall->setText( tr( "Plot All Distros" ) );
    788800#endif
     801   pb_refresh->setEnabled( true );
     802   pb_reset  ->setEnabled( true );
     803
     804   // Update status text box
     805   QString distripl = QString( model.description ).section( ".", -3, -3 );
     806DbgLv(1) << "LD:  model:" << model.description;
     807   QString disdesc  = tsys.label;
     808   QString disinfo  = te_distr_info->toPlainText() + "\n  "
     809                      + distripl + "  " + disdesc + tr( "  loaded" );
     810   te_distr_info->setText( disinfo );
    789811DbgLv(1) << "LD: RETURN";
    790812}
     
    10301052{
    10311053qDebug() << "mdlpar:Set Model Parameters";
     1054   // Open Model Parameters dialog
    10321055   US_ModelParams mpdiag( alldis, this );
    10331056
     1057   // Use parameters returned
    10341058   if ( mpdiag.exec() == QDialog::Accepted )
    10351059   {  // Redo text box summarizing models; calculate vectors
    1036       QString dinfo;
    10371060      QString mdesc     = mdescs[ 0 ].section( mdescs[ 0 ].left( 1 ), 1, 1 );
    1038 //      int kk            = mdesc.indexOf( "-run" );
    1039 //      mdesc             = ( kk > 0 ) ? QString( mdesc ).mid( ( kk + 1 ) )
    1040 //                                     : mdesc;
    10411061      mdesc             = QString( mdesc ).left( 50 );
    1042       dinfo             = tr( "Run:\n  " ) + mdesc + "...\n\n"
     1062      QString dinfo     = tr( "Run:\n  " ) + mdesc + "...\n\n"
    10431063                        + tr( "  D2O Percent  Density  Label  MDescr.\n" );
    10441064
    10451065      for ( int jj = 0; jj < alldis.size(); jj++ )
    1046       {
     1066      {  // Compose and display a distribution line for distro info box
    10471067         double d2opct     = alldis[ jj ].d2opct;
    10481068         double bdens      = alldis[ jj ].bdensity;
     
    10671087{
    10681088DbgLv(1) << "sel_x:  ival" << ival;
     1089   plot_x     = ival;
     1090
    10691091#if 0
    10701092   const QString xlabs[] = {   "mass", "f/f0",   "rh", "vbar",     "s" };
     
    10901112   const double  xincs[] = {     0.01,  0.01, 1000.0,  0.01, 1e-9, 1e-9 };
    10911113#endif
    1092 
    1093    plot_x     = ival;
    10941114
    10951115#if 0
     
    11091129 << xincs[plot_x];
    11101130
    1111    build_bf_dists();
    1112    build_bf_vects();
     1131   build_bf_dists();    // Build the boundary fraction distributions
     1132   build_bf_vects();    // Build the boundary fraction vectors
    11131133
    11141134//   set_limits();
    11151135
    1116    plot_data();
     1136   plot_data();         // Plot data
    11171137}
    11181138
     
    11521172      j1               = j2 - 1;
    11531173      S_Solute sol_bf  = tsys->bo_distro[ j2 ];
    1154       sol_bf.f         = bfrac;
     1174      sol_bf.f         = bfrac;                   // Boundary fraction
    11551175      double x1        = tsys->bo_distro[ j1 ].f;
    11561176      double x2        = tsys->bo_distro[ j2 ].f;
     
    11581178      double y1        = tsys->bo_distro[ j1 ].s;
    11591179      double y2        = tsys->bo_distro[ j2 ].s;
    1160       sol_bf.s         = y1 + ( y2 - y1 ) * xfrac;
     1180      sol_bf.s         = y1 + ( y2 - y1 ) * xfrac; // Interpolated s
    11611181      y1               = tsys->bo_distro[ j1 ].k;
    11621182      y2               = tsys->bo_distro[ j2 ].k;
    1163       sol_bf.k         = y1 + ( y2 - y1 ) * xfrac;
     1183      sol_bf.k         = y1 + ( y2 - y1 ) * xfrac; // Interpolated f/f0
    11641184      y1               = tsys->bo_distro[ j1 ].c;
    11651185      y2               = tsys->bo_distro[ j2 ].c;
    1166       sol_bf.c         = y1 + ( y2 - y1 ) * xfrac;
     1186      sol_bf.c         = y1 + ( y2 - y1 ) * xfrac; // Interpolated concentration
    11671187      y1               = tsys->bo_distro[ j1 ].w;
    11681188      y2               = tsys->bo_distro[ j2 ].w;
    1169       sol_bf.w         = y1 + ( y2 - y1 ) * xfrac;
     1189      sol_bf.w         = y1 + ( y2 - y1 ) * xfrac; // Interpolated wt
    11701190      y1               = tsys->bo_distro[ j1 ].v;
    11711191      y2               = tsys->bo_distro[ j2 ].v;
    1172       sol_bf.v         = y1 + ( y2 - y1 ) * xfrac;
     1192      sol_bf.v         = y1 + ( y2 - y1 ) * xfrac; // Interpolated vbar
    11731193      y1               = tsys->bo_distro[ j1 ].d;
    11741194      y2               = tsys->bo_distro[ j2 ].d;
    1175       sol_bf.d         = y1 + ( y2 - y1 ) * xfrac;
     1195      sol_bf.d         = y1 + ( y2 - y1 ) * xfrac; // Interpolated D
    11761196if (kk<3 || (kk+4)>nsolbf)
    11771197 DbgLv(1) << "BldBf:  kk bfrac" << kk << bfrac << "j1 j2" << j1 << j2
     
    11791199 
    11801200
    1181       tsys->bf_distro << sol_bf;
    1182       bfrac           += bfincr;
     1201      tsys->bf_distro << sol_bf;  // Solute point at boundary fraction
     1202      bfrac           += bfincr;  // Bump to next boundary fraction
    11831203   }
    11841204
     
    11891209void US_Density_Match::build_bf_dists()
    11901210{
    1191    for ( int jj = 0; jj < alldis.size(); jj++ )
    1192    {
     1211   bool diff_avg    = true;
     1212   int ndists       = alldis.size();
     1213   for ( int jj = 0; jj < ndists; jj++ )
     1214   {  // Build bfrac distro for each model
    11931215      build_bf_distro( jj );
     1216   }
     1217
     1218   if ( diff_avg )
     1219   {  // Replace diffusion coefficients with weight averages
     1220      int npoints      = alldis[ 0 ].bf_distro.size();
     1221      for ( int jj = 0; jj < npoints; jj++ )
     1222      {
     1223         double dsum      = 0.0;
     1224         double wsum      = 0.0;
     1225         for ( int ii = 0; ii < ndists; ii++ )
     1226         {  // Accumulate (weighted) sum and sum of weights
     1227//            double dwt       = alldis[ ii ].bf_distro[ jj ].c;
     1228            double dwt       = 1.0;
     1229            dsum            += ( alldis[ ii ].bf_distro[ jj ].d * dwt );
     1230            wsum            += dwt;
     1231         }
     1232         double dval      = dsum / wsum;  // average diffusion coefficient
     1233         for ( int ii = 0; ii < ndists; ii++ )
     1234         {  // Propagate to all models
     1235            alldis[ ii ].bf_distro[ jj ].d = dval;
     1236         }
     1237      }
    11941238   }
    11951239}
     
    12111255   v_bfracs.reserve( npoints );
    12121256   for ( int jj = 0; jj < npoints; jj++ )
    1213    {
     1257   {  // Grab each fraction from 1st model's distro
    12141258      v_bfracs << alldis[ 0 ].bf_distro[ jj ].f;
    12151259   }
     
    12261270      // Save density for each model distribution
    12271271      double density   = alldis[ ii ].bdensity;
    1228       v_dens << density;
     1272      v_dens << density;                // Vector of densities
    12291273      v_sedcs[ ii ].clear();
    12301274      v_difcs[ ii ].clear();
     
    12341278      // Build vectors of s and D for this model
    12351279      for ( int jj = 0; jj < npoints; jj++ )
    1236       {
     1280      {  // Build s,D vectors for this model
    12371281         v_sedcs[ ii ] << alldis[ ii ].bf_distro[ jj ].s;
    12381282         v_difcs[ ii ] << alldis[ ii ].bf_distro[ jj ].d;
     
    12521296      v_seds.clear();
    12531297      for ( int ii = 0; ii < ndists; ii++ )
    1254       {
     1298      {  // Build sed coeffs vector across models
    12551299         v_seds << alldis[ ii ].bf_distro[ jj ].s;
    12561300      }
    1257       double* xx       = v_seds.data();
    1258       double* yy       = v_dens.data();
     1301      double* xx       = v_seds.data();  // X is sed coeffs
     1302      double* yy       = v_dens.data();  // Y is densities
    12591303      double slope, intcept, sigma, corre;
    12601304      US_Math2::linefit( &xx, &yy, &slope, &intcept, &sigma, &corre, ndists );
     
    12771321   for ( int jj = 0; jj < npoints; jj++ )
    12781322   {
    1279    }
     1323      double sedco     = alldis[ 0 ].bf_distro[ jj ].s * 1.0e-13;
     1324      double difco     = alldis[ 0 ].bf_distro[ jj ].d * 1.0e-7;
     1325      double vbari     = v_vbars[ jj ];
     1326      double mmass     = sedco * R_GC * K20 / ( difco * ( 1.0 - vbari * DENS_20W ) );
     1327      mmass            = qAbs( mmass );
     1328      v_mmass << mmass;
     1329   }
     1330DbgLv(1) << "BldVc: mm 0 1 k n" << v_mmass[0] << v_mmass[1]
     1331 << v_mmass[npoints-2] << v_mmass[npoints-1];
    12801332
    12811333   // Compute hydrodynamic radius values and build the vector
     
    12841336   for ( int jj = 0; jj < npoints; jj++ )
    12851337   {
     1338      double difco     = alldis[ 0 ].bf_distro[ jj ].d * 1.0e-7;
     1339      double frico     = R_GC * K20 / ( difco * AVOGADRO );
     1340      double hyrad     = frico / ( 6.0 * M_PI * VISC_20W );
     1341      v_hrads << hyrad;
    12861342   }
    12871343
    12881344   // Compute frictional ratio values and build the vector
     1345   const double a_third = ( 1.0 / 3.0 );
    12891346   v_frats.clear();
    12901347   v_frats.reserve( npoints );
    12911348   for ( int jj = 0; jj < npoints; jj++ )
    12921349   {
    1293    }
    1294 
     1350      double difco     = alldis[ 0 ].bf_distro[ jj ].d * 1.0e-7;
     1351      double vbari     = v_vbars[ jj ];
     1352      double rzero     = pow( ( ( 0.75 / M_PI ) * vbari ), a_third );
     1353      double frico     = R_GC * K20 / ( difco * AVOGADRO );
     1354      double fzero     = 6.0 * M_PI * VISC_20W * rzero;
     1355      double frrat     = frico / fzero;
     1356      v_frats << frrat;
     1357   }
     1358
     1359#if 0
     1360
     1361*** Mi = si*R*T/(Di_avg*(1-vbari*rho))
     1362
     1363i = boundary fraction from 0 to 1
     1364
     1365Mi*vbari/N = Volume of moleculei
     1366
     1367V=4/3 * pi*r_0^3    (3/(4*pi) *v)^1/3 = r_0
     1368
     1369f_0i = 6 * pi * eta * r_0i
     1370
     1371fi = RT/(N*Di)
     1372
     1373*** fi/f_0i
     1374
     1375*** ri = fi/(6 * pi * eta)   <-- hydrodynamic radius
     1376
     1377*** --> distributions to plot, let user choose which
     1378#endif
    12951379#if 0
    12961380   QVector< double >             v_bfracs;
  • trunk/programs/us_density_match/us_density_match.h

    r2807 r2808  
    2525   QList< S_Solute >   bo_distro;      // Boundary distro w/ orig points
    2626   QList< S_Solute >   bf_distro;      // Boundary fractions distro
    27    QString             run_name;
    28    QString             analys_name;
    29    QString             method;
    30    QString             editGUID;
    31    QString             solutionGUID;
    32    QString             label;
    33    int                 distro_type;
    34    int                 plot_x;
    35    int                 solutionID;
    36    double              d2opct;
    37    double              bdensity;
     27   QString             run_name;       // Distro run name
     28   QString             analys_name;    // Distro analysis name
     29   QString             method;         // Model method (e.g., "2DSA")
     30   QString             editGUID;       // Associated edit GUID
     31   QString             solutionGUID;   // Associated solution GUID
     32   QString             label;          // Distro label (channel description)
     33   int                 distro_type;    // Distro type flag
     34   int                 solutionID;     // Associated solution db ID
     35   double              d2opct;         // D2O percent for distro
     36   double              bdensity;       // Distro buffer density
    3837} DisSys;
    3938
     
    6867 
    6968      QwtCounter*   ct_resolu;
    70       QwtCounter*   ct_xreso;
    71       QwtCounter*   ct_yreso;
    72       QwtCounter*   ct_zfloor;
    7369      QwtCounter*   ct_plt_kmin;     
    7470      QwtCounter*   ct_plt_kmax;     
     
    7672      QwtCounter*   ct_plt_smax;     
    7773      QwtCounter*   ct_plt_dlay;     
    78       QwtCounter*   ct_curr_distr;
    79       QwtCounter*   ct_tolerance;
    8074      QwtCounter*   ct_division;
    8175      QwtCounter*   ct_smoothing;
     
    8983      US_Disk_DB_Controls* dkdb_cntrls;
    9084
    91       QPushButton*  pb_pltall;
    92       QPushButton*  pb_stopplt;
    9385      QPushButton*  pb_refresh;
    9486      QPushButton*  pb_reset;
     
    119111      QButtonGroup* bg_x_axis;
    120112
    121       QVector< DisSys >             alldis;
     113      QVector< DisSys >             alldis;    // All distributions
    122114
    123       QVector< double >             v_bfracs;
    124       QVector< double >             v_vbars;
    125       QVector< double >             v_mmass;
    126       QVector< double >             v_hrads;
    127       QVector< double >             v_frats;
    128       QVector< QVector< double > >  v_sedcs;
    129       QVector< QVector< double > >  v_difcs;
     115      QVector< double >             v_bfracs;  // Boundary fraction vector
     116      QVector< double >             v_vbars;   // Vector of vbars per fraction
     117      QVector< double >             v_mmass;   // Vector of molar masses per fraction
     118      QVector< double >             v_hrads;   // Vector of hydro radii per fraction
     119      QVector< double >             v_frats;   // Vector of frict ratios per fraction
     120      QVector< QVector< double > >  v_sedcs;   // Vector of sedi coeff vectors per distro
     121      QVector< QVector< double > >  v_difcs;   // Vector of diff coeff vectors per distro
    130122
    131123      double        resolu;
  • trunk/programs/us_density_match/us_model_params.cpp

    r2806 r2808  
    126126   QFont font( US_GuiSettings::fontFamily(), US_GuiSettings::fontSize() );
    127127   QFontMetrics fm( font );
    128    int fhigh = fm.lineSpacing();
     128//   int fhigh = fm.lineSpacing();
     129//   int lhigh = fhigh * 10 + 12;
    129130   int fwide = fm.width( QChar( '6' ) );
    130    int lhigh = fhigh * 10 + 12;
    131131   int lwide = fwide * ( maxdlen + 2 );
    132132
Note: See TracChangeset for help on using the changeset viewer.