1

これは非常に奇妙な問題です...以下の関数で cout を削除すると、正しい/期待される結果の出力とガベージ値の出力が停止します。(つまり、出力するデータをまだ実行していますが、間違っています)。何か案は?

bool extract_tension(std::vector<double> &interfacial_tension_trap,
       std::vector<double> &interfacial_tension_simp,
       const std::string data,
       const unsigned int num_slabs,
       const double z_min, const double z_max)
{

   //start @ first number
   unsigned int start = 17;
   unsigned int end = 17;

   std::vector<double> px;
   std::vector<double> py;
   std::vector<double> pz;

   std::vector<double> pn_minus_pt;

   double result_simp=0.0;
   double result_trap=0.0;

   //skip timestep entry
   end=get_next_space(start, data);

   for(unsigned int counter=0; counter<num_slabs;counter++)
   {
     start = end+2;

     end=get_next_space(start, data);
     px.push_back(atof(data.substr(start,(end-start+1)).c_str()));
     //skip the space
     start = end+2;
     end=get_next_space(start, data);
     py.push_back(atof(data.substr(start,(end-start+1)).c_str()));
     //skip the space
     start = end+2;
     end=get_next_space(start, data);
     pz.push_back(atof(data.substr(start,(end-start+1)).c_str()));

     //calculate pressure difference
     // WARNING : Unit conversion ahead
     // NAMD outputs pressure in bars and distance in Angstroms
     // we want an integrated result of mN/m, instead.
     // 1 Angstrom = 1e-10 m
     // 1 bar = 1e8 mN/m^2
     // net conversion -- 1e-2 
     pn_minus_pt.push_back((pz[counter]-0.5*(px[counter]+py[counter]))*0.01);
       std::cout << "Current del_P : " 
   << (pz[counter]-0.5*(px[counter]+py[counter]))*0.01
   << std::endl;
   }
   calculate_trapezoid(pn_minus_pt, num_slabs, z_min, z_max, result_trap);
   interfacial_tension_trap.push_back(result_trap);
   calculate_simpson(pn_minus_pt, num_slabs, z_min, z_max, result_simp);
   interfacial_tension_simp.push_back(result_simp);
}

どうやら、印刷ステートメントでベクターのいずれかに触れるだけで、プログラムが正しく実行されるようになります (つまり、px、py、または pz を含む印刷)。

完全なプログラムは次のとおりです。

/*********************************
 *
 * NAME: Interfacial Tension Calculator
 * VERSION: 0.1
 * AUTHOR: Jason R. Mick
 * GROUP: Wayne State University, Potoff Group
 * COPYRIGHT: (c) Jason R. Mick 2010
 * DATE: August 9, 2010
 *
 * CHANGE LOG
 * VERSION    DATE         COMMENTS
 *----------------------------------------------------------------------
 *  0.1   Aug. 9, 2010  Finished basic code, sans debugging  
 *  0.5   Aug  10, 2010 Compiled and tested code fixed error in Simpson's
 *                      method where results were being divided rather
 *                      than multiplied.                       
 *
 *
 * FULL NOTES:
 *----------------------------------------------------------------------
 * You can compile this program by typing:
 * g++ main.cc -o it_util
 *
 * You can run this program by typing:
 * it_util <filename>.log <# slabs> <z-min> <z-max>
 *
 * where z-min and z-max represent the z-axis boundaries of the system,
 * e.g.--
 * it_util my_file.log 140 0.0 80.0
 *  
 * This program only works with NAMD *.log file output
 * The pressure profile MUST be turned on in your *.conf file
 * for the pressure profile info to dump to the *.log file.  This
 * program requires that info.
 *
 * This program can handle 1,000+ slabs, but it has a limit to the
 * character buffer and thus VERY large slab counts may cause it to fail.
 *
 * A standard Composite Simpson numerical integration method is used,
 * which assumes a non-smooth data set.
 *
 * The interfacial tension is integrated at each step and then averaged
 * so pertinent statistics can be gathered.
 *
 * You can redirect the output to store the interfacial tension
 * statistics as follows:
 * it_util <filename>.log <# slabs> <z-min> <z-max> > <my_file>.out
 * 
 *******************************************/

#include <stdio.h>
#include <math.h>

#include <iostream>
#include <vector>
#include <fstream>

#include <sys/stat.h> 

//Turn on to enable all interfacial 
//tension results to be printed, pre-averaging
//#define DEBUG true

void start_integrations(const std::string filename, 
   const unsigned int num_slabs,
   const double z_min, const double z_max);


int main ( int argc, char *argv[] )
{
  struct stat file_info; 
  std::string filename = argv[1];
  int slab_count;
  double z_min;
  double z_max;

  if ( argc != 5 ) /* argc should be 3 for correct execution */
  {
    /*Print out proper args syntax */
    std::cout << "ERROR: Missing arguments!" << std::endl 
       << "Proper syntax:" << std::endl
       << "it_util <my_file>.log <# of slabs> <z-coord start>"
       << "<z-coord end>"
       << std::endl;
  }
  if(stat(argv[1],&file_info)==0)
  {
    try
    {
      slab_count = atoi(argv[2]);
      if (slab_count > 2)
      {
 try
 {
   z_min = atof(argv[3]);
   try 
   {
     z_max = atof(argv[4]);
     start_integrations(filename, 
          static_cast<unsigned int>(slab_count),
          z_min,
          z_max);
   }
   catch( char * str ) 
   {
     /*invalid integer third input*/
     std::cout << "Invalid input -- fourth argument was invalid "
        << "decimal number, should be standard " << std::endl
        << "decimal type entry..." << std::endl
        << "I.E." << std::endl
        << "it_util my_file.log 140 0.0 80.0" << std::endl;

   }
 }
 catch( char * str ) 
 {
   /*invalid integer third input*/
   std::cout << "Invalid input -- third argument was invalid "
      << "decimal number, should be standard " << std::endl
      << "decimal type entry..." << std::endl
      << "I.E." << std::endl
      << "it_util my_file.log 140 0.0 80.0" << std::endl;

 }
      }
      else
      { 
 /*invalid integer secondary input*/
 std::cout << "Invalid input -- second argument was invalid integer, "
    << "should be unsigned integer 2 or greater..." << std::endl
    << "I.E." << std::endl
    << "it_util my_file.log 140 0.0 80.0" << std::endl;
      }
    }
    catch( char * str ) 
    {
      /*non integer secondary input*/
      std::cout << "Invalid input -- second argument was non-integer, "
  << "should be unsigned integer 2 or greater..." << std::endl
  << "I.E." << std::endl
  << "it_util my_file.log 140 0.0 80.0" << std::endl;
    }
  }
  else
  {
    /*invalid filename case...*/
    std::cout << "File " << filename << "does not exist!" << std::endl
       << "Please choose valid file!" << std::endl;
  }

  return 1;
}

bool calculate_simpson(const std::vector<double> my_values, 
        const unsigned int num_points, 
        const double x_min, const double x_max, 
        double &results)
{
   bool ret_val = false;
   bool is_even = true;
   double h;

   if (my_values.size() >= 2)
   {
      h = (x_max-x_min)/num_points;
      results+=my_values.front();
      for (unsigned int counter=1; counter<num_points-1;counter++)
      {
         if (is_even)
         {
            results+=4*my_values[counter];
         }
         else
         {
            results+=2*my_values[counter];
         }
         is_even = !is_even;
      }
      results+=my_values.back();
      results*=(h/3);
      ret_val=true;
   }
   return ret_val;
}

bool calculate_trapezoid(const std::vector<double> my_values, 
    const unsigned int num_points, 
                         const double x_min, const double x_max, 
    double &results)
{
   bool ret_val = false;

   double x_incr = (x_max-x_min)/(num_points-1);

   if (my_values.size() >= 2)
   {      
      for (unsigned int counter=1; counter<num_points-1; counter++)
      {
         results+=(x_incr/2)*(my_values[counter]+my_values[counter-1]);
      }
   }
   return ret_val;
}

unsigned int get_next_space(const unsigned int start,
       const std::string data)
{
   unsigned int counter=start;

   while (data.length() > counter &&
   data.substr(counter,1).compare(" ") != 0)
   {    
     counter++;
   }

   //if end of string, add one
   if ( data.length() == counter)
     counter++;
   return (counter-1);
}

bool extract_tension(std::vector<double> &interfacial_tension_trap,
       std::vector<double> &interfacial_tension_simp,
       const std::string data,
       const unsigned int num_slabs,
       const double z_min, const double z_max)
{

   //start @ first number
   unsigned int start = 17;
   unsigned int end = 17;

   std::vector<double> px;
   std::vector<double> py;
   std::vector<double> pz;

   std::vector<double> pn_minus_pt;

   double result_simp=0.0;
   double result_trap=0.0;

   //skip timestep entry
   end=get_next_space(start, data);

   for(unsigned int counter=0; counter<num_slabs;counter++)
   {
     start = end+2;

     end=get_next_space(start, data);
     px.push_back(atof(data.substr(start,(end-start+1)).c_str()));
     //skip the space
     start = end+2;
     end=get_next_space(start, data);
     py.push_back(atof(data.substr(start,(end-start+1)).c_str()));
     //skip the space
     start = end+2;
     end=get_next_space(start, data);
     pz.push_back(atof(data.substr(start,(end-start+1)).c_str()));

     //calculate pressure difference
     // WARNING : Unit conversion ahead
     // NAMD outputs pressure in bars and distance in Angstroms
     // we want an integrated result of mN/m, instead.
     // 1 Angstrom = 1e-10 m
     // 1 bar = 1e8 mN/m^2
     // net conversion -- 1e-2 
     pn_minus_pt.push_back((pz[counter]-0.5*(px[counter]+py[counter]))*0.01);
       std::cout << "Current del_P : " 
   << (pz[counter]-0.5*(px[counter]+py[counter]))*0.01
   << std::endl;
   }
   calculate_trapezoid(pn_minus_pt, num_slabs, z_min, z_max, result_trap);
   interfacial_tension_trap.push_back(result_trap);
   calculate_simpson(pn_minus_pt, num_slabs, z_min, z_max, result_simp);
   interfacial_tension_simp.push_back(result_simp);
}

double average_vector(std::vector<double> my_vector)
{
   double average_val=0.0;

   for(unsigned int counter=0; counter< my_vector.size(); counter++)
   {
     average_val+=my_vector[counter]/my_vector.size();
   }

   return average_val;
}

double std_dev_vector(std::vector<double> my_vector)
{
   double std_deviation=0.0;
   double average_val = average_vector(my_vector);

   for(unsigned int counter=0; counter< my_vector.size(); counter++)
   {
     std_deviation+=(my_vector[counter]-average_val)*
       (my_vector[counter]-average_val);
   }
   std_deviation=sqrt(std_deviation);

   return std_deviation;
}

void start_integrations(const std::string filename, 
   const unsigned int num_slabs,
   const double z_min, const double z_max)
{
   std::ifstream in_file;
   std::vector<double> interfacial_tension_trap;
   std::vector<double> interfacial_tension_simp;
   std::string current_line;
   char * cstr_line;
   bool data_grab_success = true;

   in_file.open(filename.c_str(), std::ifstream::in);
   while (!in_file.eof() && data_grab_success)
   {
     cstr_line=(char *) malloc(sizeof(char)*65536);
     //get new line
     in_file.getline(cstr_line,65536);
     current_line = cstr_line;
     free(cstr_line);
     if (current_line.substr(0,15).compare("PRESSUREPROFILE")==0)
     {
       //pressure profile found!

       //process line to get the interfacial tension, check that it succeeded
       data_grab_success = extract_tension(interfacial_tension_trap,
       interfacial_tension_simp,
        current_line,
        num_slabs,
        z_min,
        z_max);
     }
   }
   in_file.close();

   //print stats
   std::cout << "Interfacial Tension (Trapezoid Method): " 
      << average_vector(interfacial_tension_trap) << std::endl
      << "Standard Deviation (Trapezoid Method): " 
      << std_dev_vector(interfacial_tension_trap) << std::endl
      << "Interfacial Tension (Composite Simpson's Method): " 
      << average_vector(interfacial_tension_simp) << std::endl
      << "Standard Deviation (Composite Simpson's Method): " 
      << std_dev_vector(interfacial_tension_simp) << std::endl;
}

データのサンプル セットは次のとおりです。

Removed... see explanation at end of post for link to data to use.

次のようにコンパイルします。

g++ main.cc -o it_util

次のコマンドを使用して実行します。

it_util equil2_NVT_PP_318Slabs.log 318 0.0 318.0 > temp.out

#ifdef "debug" ス​​テートメントについて誰かがコメントする前に参考までに、それらはデータ ダンプ用であることに注意してください。以前にGDBを使用したことがあります。私がこれを言わなかったら、誰かが「gdb の使い方を学んでください」とコメントするでしょう。この場合、プログラムは非常に多くの反復をループするため、GDB は有用な情報を提供しません。印刷出力は出力ファイル DO にダンプされます。

注:
実際、(上記のデータ セクションで) 解析されているファイルの切り詰められたバージョンを使用すると、プログラムも正しいデータを出力しないことがわかりました。元のデータ ファイルを復元すると機能しましたが、ファイルが大きすぎてここに投稿できませんでした (試してみました...)。そのため、それはオプションではありません....代わりに、完全なペーストビンをここにアップロードしました: http:/ /pastebin.com/JasbSc7B

4

2 に答える 2

2

これは奇妙です...だから私は問題を理解しました。かなりの素人ミスでした。data_grab_success でブール値の成功を返すのを忘れていました。

ただし、それは奇妙な部分を説明していません-どういうわけか、この値は自動的にfalseで満たされていました-そのcoutを実行しない限り。とにかく、私は夢中になっていたので問題が解決したことをうれしく思いますが、何も指定されていない場合に関数が戻り値を取得する場所と、cout がそれにどのように影響するかについて当惑しています...

(PS結局、GDBを使用して最終的にこれを理解しました....)

的確なアドバイスをくれた Dugan に感謝します。

于 2010-08-10T22:42:14.340 に答える
1

ここに私の最初の投稿。素晴らしいサイト。

確かに謎です。これは答えではありませんが、まだ行っていない場合は、探索するための頭の痛い一連の方法です。

  1. 入力データ (NaN など) に異常があり、cout ステートメントの計算によって FPU の状態が変更されます。それが何であるかはわかりませんが、あなたは石のない段階にいるようです.

  2. 使用している std::vector<> の実装に奇妙なバグがあり、特定の使用パターンによって何らかの方法でトリガーされ、operator[] 呼び出しに状態変更の副作用が生じます。式から iostream を除外し、単に operator[] を呼び出して調査します。物事を台無しにする特定の呼び出しに絞り込むことができる場合があります。ただし、これはばかげてありそうにありません。

  3. 出力の「不正確さ」を特徴付けます。ゴミにパターンはありますか?それは正しい出力がどうあるべきかと何らかの形で相関していますか?

  4. 理解できないバグへの通常のアプローチ: 取得できる最も単純な再現ケースが得られるまでコードをトリムし、最初の悪い結果がいつ表示されるかを正確に把握するようにインストルメント化します。

  5. 別のライブラリ、コンパイラ、OS、または CPU でバグを再現してみてください。ロングショットですが、あなたは決して知りません.それでも再現する場合、少なくともこれは実際にはあなた自身のバグであり、レンガの壁に頭をぶつけていないことを安心させました.

このアドバイスの一部は一般的なものですが、お役に立てば幸いです。クラックしたらお知らせください!

于 2010-08-10T21:46:57.750 に答える