Weird Error in C++ Program: Removing Printout Breaks Program

Question

This is a very strange problem... removing the cout in the function below causes it to stop printing the correct/expected results and printing garbage values. (i.e. it still RUNS the data it outputs, though, is wrong). Any ideas?

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);
}

Apparently just touching any of vectors with a print statement allows the program to execute correctly (i.e. a printout involving px, py, OR pz)

Here's the full program:

/*********************************
 *
 * 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;
}

And here's a sample set of data:

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

Compile like so:

g++ main.cc -o it_util

Run using the command:

it_util equil2_NVT_PP_318Slabs.log 318 0.0 318.0 > temp.out

FYI before someone comments on my #ifdef "debug" statements, please note that they are for data dumping. I have used GDB before. I'm guessing had I not said this, someone will comment "Learn to use gdb." In this case the program loops through so many iterations, GDB doesn't give me useful info, where printouts dumped to a output file DO.

NOTE:
Actually I discovered that if you use the chopped down version of the file being parsed (in the data section above) the program also doesn't output the correct data. When I restored the original data file it worked, but the file is too large to post here (I tried...) so that is not an option.... Instead I've uploaded a full pastebin to here: http://pastebin.com/JasbSc7B

Answer 1

My first post here. Great site.

That's certainly mysterious. This is not an answer, but a set of harebrained avenues to explore, if you haven't already:

1. Something unusual in your input data (NaNs, etc) that makes the math in the cout statement change the FPU state. I can't fathom what that might be, but it sounds like you're at the no-stone-unturned stage.

2. Bizarre bug in the implementation of std::vector<> you're using, triggered in some way by your specific usage pattern, making the operator[] calls have state-changing side effects. Investigate by leaving iostream out of the equation and just calling operator[]. You may be able to narrow it down to the specific call that messes things up. This is ridiculously unlikely though.

3. Characterize the output's "incorrectness". Any pattern in the garbage? Does it correlate in any way to what the correct output should be?

4. The usual approaches to incomprehensible bugs: trim the code until you have the simplest repro case you can get, instrument it to figure out exactly when the first bad result appears, and so on.

5. Try to repro the bug with a different library, compiler, OS, or CPU. Long shot, but you never know, and if still repros, at least you've reassured yourself that this is in fact your own bug and you aren't banging your head against a brick wall.

Some of this advice is kind of generic, but I hope it helps. Let us know when you crack it!

Answer 2

This is strange... So I figured out the problem. It was a pretty amateurish mistake. I forgot to return the boolean success on data_grab_success.

However that doesn't explain the strange part -- somehow this value was getting filled with false automatically -- unless I ran that cout. Anyways I'm happy the problem is solved because I was going crazy, but I'm perplexed as to where it a function gets its return value if none is specified and how cout could influence that...

(P.S. I finally figured this out by using GDB after all....)

Thanks to Dugan for giving solid advice!