Странная ошибка в программе на C ++: удаление программы для прерывания печати - PullRequest
1 голос
/ 10 августа 2010

Это очень странная проблема ... удаление метки в приведенной ниже функции приводит к тому, что она перестает печатать правильные / ожидаемые результаты и выводить значения мусора. (то есть он все еще запускает данные, которые выводит, хотя и неверен). Есть идеи?

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, OR 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, обратите внимание, что они предназначены для сброса данных. Я использовал GDB раньше. Я предполагаю, что если бы я не сказал это, кто-то прокомментирует «Научитесь использовать GDB». В этом случае программа проходит через очень много итераций, GDB не дает мне полезной информации, где распечатки выводятся в выходной файл DO.

Примечание:
На самом деле я обнаружил, что если вы используете урезанную версию анализируемого файла (в разделе данных выше), программа также не выводит правильные данные. Когда я восстановил исходный файл данных, он работал, но файл слишком велик для размещения здесь (я пытался ...), так что это не вариант .... Вместо этого я загрузил полную вставку сюда: http://pastebin.com/JasbSc7B

Ответы [ 2 ]

2 голосов
/ 11 августа 2010

Это странно ... Так что я разобрался с проблемой. Это была довольно дилетантская ошибка. Я забыл вернуть логический успех на data_grab_success.

Однако это не объясняет странную часть - каким-то образом это значение автоматически заполнялось ложным - если я не запускал этот cout. В любом случае, я счастлив, что проблема решена, потому что я сходил с ума, но я озадачен тем, откуда функция возвращает свое значение, если ничего не указано, и как cout может повлиять на это ...

(П.С. В конце концов я понял это, используя GDB ....)

Спасибо Дугану за дельный совет!

1 голос
/ 11 августа 2010

Мой первый пост здесь. Отличный сайт.

Это, конечно, таинственно. Это не ответ, а целый ряд возможностей для изучения, если вы еще этого не сделали:

  1. Что-то необычное во ваших входных данных (NaNs и т. Д.), Которое заставляет математику в операторе cout изменять состояние FPU. Я не могу понять, что это может быть, но звучит так, будто вы находитесь на стадии, где нет камня.

  2. Странная ошибка в используемой вами реализации std :: vector <>, вызванная каким-то образом вашим конкретным шаблоном использования, в результате чего вызовы оператора [] вызывают изменяющиеся состояния побочных эффектов. Расследуйте, оставив iostream вне уравнения и просто вызвав оператор []. Вы можете быть в состоянии сузить это до определенного вызова, который портит вещи. Хотя это невероятно маловероятно.

  3. Охарактеризуйте «некорректность» вывода. Какой-нибудь шаблон в мусоре? Соответствует ли это каким-либо образом правильному выводу?

  4. Обычные подходы к непонятным ошибкам: обрезать код до тех пор, пока у вас не появится простейший случай повторения, который вы можете получить, дать ему инструмент для точного определения момента появления первого плохого результата и т. Д.

  5. Попробуйте воспроизвести ошибку с другой библиотекой, компилятором, ОС или ЦП. Длинный выстрел, но вы никогда не знаете, и если все еще репродукции, по крайней мере, вы заверили себя, что это на самом деле ваш собственный баг, и вы не бьетесь головой о кирпичную стену.

Некоторые из этих советов являются общими, но я надеюсь, что это поможет. Дайте нам знать, когда вы взломаете его!

...