Domanda

Questo è un problema molto strano ... rimuovendo il cout nella funzione di seguito le cause per interrompere la stampa della corretta / risultati attesi e la stampa di valori della spazzatura. (Cioè esso gestisce ancora dati emette, però, è sbagliato). Tutte le idee?

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

A quanto pare basta toccare qualsiasi di vettori con un comunicato stampa permette al programma di eseguire correttamente (vale a dire una stampa che coinvolge px, py, o PZ)

Ecco il programma completo:

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

Ed ecco un set di dati di esempio:

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

Compila questo modo:

g++ main.cc -o it_util

Esegui usando il comando:

it_util equil2_NVT_PP_318Slabs.log 318 0.0 318.0 > temp.out

FYI prima commenti qualcuno sul mio #ifdef dichiarazioni "di debug", si ricorda che essi sono per lo scarico dei dati. Ho usato GDB prima. Sto indovinando se non avessi detto questo, qualcuno commenterà "imparare ad usare gdb". In questo caso i loop di programma attraverso tante iterazioni, GDB non mi danno informazioni utili, in cui le stampe scaricati in un file di output DO.

NOTA:
In realtà ho scoperto che se si utilizza la versione abbattuto del file da analizzato (nella sezione dati sopra) il programma non lo fa anche in uscita i dati corretti. Quando ho ripristinato il file di dati originale ha funzionato, ma il file è troppo grande per postare qui (ho provato ...) in modo che non è un'opzione .... Invece ho caricato un pastebin completo qui: http://pastebin.com/JasbSc7B

È stato utile?

Soluzione

Questo è strano ... Così ho capito il problema. E 'stato un errore da dilettanti abbastanza. Ho dimenticato di restituire il successo booleano su data_grab_success.

Tuttavia questo non spiega la strana parte - in qualche modo questo valore è stato sempre pieno di falsa automaticamente - a meno che non mi sono imbattuto che cout. In ogni modo sono felice il problema è risolto perché stavo impazzendo, ma sto perplesso da dove una funzione prende il suo valore di ritorno, se non viene specificato e come potrebbe influenzare cout che ...

(P.S. Alla fine ho capito questo utilizzando GDB dopo tutto ....)

Grazie a Dugan per dare consigli solido!

Altri suggerimenti

Il mio primo post qui. Grande sito.

Questo è certamente misterioso. Questa non è una risposta, ma un insieme di viali harebrained da esplorare, se non l'hai già:

  1. Qualcosa di insolito nei dati di input (NaNs, ecc) che rende la matematica nel cambiamento cout stato FPU. Io non riesco a scandagliare ciò che potrebbe essere, ma suona come se fossi in fase di non-pietra-intentato.

  2. Bizarre bug nell'implementazione di std :: vector <> si sta utilizzando, innescato in qualche modo dal vostro modello di utilizzo specifico, rendendo l'operatore [] chiamate hanno effetti collaterali che cambia di stato. Indagare lasciando iostream fuori dall'equazione e basta chiamare operatore []. Si può essere in grado di restringere verso il basso per la specifica chiamata che scombina le cose. Questo è ridicolmente anche se improbabile.

  3. caratterizzano l'uscita di "scorrettezza". Ogni modello nella spazzatura? Ha correlato in alcun modo a ciò che l'uscita corretta dovrebbe essere?

  4. I normali approcci bug incomprensibili: ritagliare il codice fino ad avere il caso Repro più semplice si può ottenere, strumento per capire esattamente quando le prime appare risultato cattivi, e così via

  5. Prova a Repro il bug con un diverso biblioteca, compilatore, sistema operativo, o CPU. Campo lungo, ma non si sa mai, e se Repros ancora, almeno ti hanno rassicurato lei che questo è in realtà il proprio bug e non si sbattere la testa contro un muro di mattoni.

Alcuni di questi consigli è una specie di generica, ma spero che aiuta. Fateci sapere quando si crepa esso!

Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top