Binning an array of data.

[SpringLC]

Hi,

I have an 1 dimensional array of values. I did plot a histogram of the data in python (using the hist function), and it turned out great. But now I want to plot the histogram with a log-log scale. So I thought I should bin the data myself in C++ and then plot it (the usual way) in python. But I'm a beginner in C++, so I was wondering if someone knows a simple way of explaining to me what to do.

So basicly - I have an array (size 1e8) which I want to divide into equally spaced logarithm (base 10) bins.

I'll post the header file with the code that creates the array itself:
(E_phot[i] is the array of values that should be binned)

  /********************************************************************/
/*						PHOTON PROPERTIES 							*/
/********************************************************************/
/* Initializing photons... */
/* ----- Photon Varaiable Declarations ----- */

/********************************************************************/
/* 1. Photon Energies (BB distribution): */
/********************************************************************/
/* Declaring four random numbers for BB distribution */
double rn1_BB_dist,rn2_BB_dist, rn3_BB_dist, rn4_BB_dist;
/* Decalre Array (Size = N_photons): Photon energies*/
vector<double> E_phot(N_photons) ;
/* Array (size = 50): for jsum in algorithm Appendix A2 */
double jsum[50];
/* Decalaring random number alpha (in algoritm Appendix A2)*/
int alpha;
/* ----- Function: Initializes :
					1. Photon momentum
					2. Photon energies */
void initialize_photons()
{
/********************************************************************/
/** DECLARE A RANDOM NUMBER GENERATOR: MERSENNE TWISTER**/
/********************************************************************/
mt19937 rng;
/*Use hardware parameters to generate a random seed:*/
random_device randomSeed;
/*Set the seed from the random device:*/
rng.seed(randomSeed());

/********************************************************************/
/* 2. Photon Energies (BB distribution): */
/********************************************************************/
/* ---- jsum array ---- */
/* Filling in the first entry of jsum....*/
jsum[0] = 1.0;
/* FOR LOOP: the subsequent entries for jsum: */
for(int i = 1; i < 50; i++)
	{
		jsum[i] = jsum[i-1] + 1.0/((i+1)*(i+1)*(i+1));
	}// end FOR LOOP: the subsequent entries for jsum:
/*Put jsum entries into a vector(double): */
vector<double> jsum_vect(jsum, jsum + 50);
/* Define a iterator for the jsum vector : */
vector<double> :: iterator low ;
/* draw four random numbers for BB photon energy */
uniform_real_distribution<double> BB_dist_rn1(0.0,1.0);
uniform_real_distribution<double> BB_dist_rn2(0.0,1.0);
uniform_real_distribution<double> BB_dist_rn3(0.0,1.0);
uniform_real_distribution<double> BB_dist_rn4(0.0,1.0);
/*FOR LOOP: BB DISTRIBUTION - Energy of photons with Planck distribution */
/** Open file to save photon energies: BB_Photon_energies.txt **/
ofstream myfile;
  myfile.open ("BB_Photon_energies.txt");
  //

for(int i = 0; i <N_photons; i++)
	{
		/* Obtain four random numbers */
		rn1_BB_dist = BB_dist_rn1(rng);
		rn2_BB_dist = BB_dist_rn2(rng);
		rn3_BB_dist = BB_dist_rn3(rng);
		rn4_BB_dist = BB_dist_rn4(rng);
		/* Calculate alpha : Use lower_bound function in <algorithm>*/
		low = lower_bound(jsum_vect.begin(),jsum_vect.end(), 1.202*rn1_BB_dist);
		alpha = low - jsum_vect.begin() + 1;
		/* Asign energy of each BB photon to array E_phot[i] */
		E_phot[i] = - ((kB*T_phot)/alpha)*log(rn2_BB_dist*rn3_BB_dist*rn4_BB_dist);
		/** Save photon energies to file BB_Photon_energies.txt **/
        //cout << E_phot[i] << endl;
        myfile << i << ',' << E_phot[i] << endl;

	}//END FOR LOOP: BB DISTRIBUTION
	/** Close file to save photon energies: BB_Photon_energies.txt **/
	myfile.close();
/********************************************************************/
} // end: Function - void initialize photons()

[lastchance]

Just write log10( E_phot[i] ) to file instead of E_phot[i].

[jonnin]

since you are using a vector you can actually do most of the work here. I dunno what python's function does, but you can use sort on the vector of logs to order them and then iterate over it once to count how many per bucket. You could even write just the count file out (it will be tiny) .. that is the # of items in each bucket .... very easily, and you can plot that from anything (excel, python, maple/matlab, whatever). If python does all the work, you don't need this, but if you want to plot it in something else, it might be useful, and going from 1e8 values to 10 or 20 values or whatever makes dealing with the data a little easier. You can even do a sideways cheeze histogram by writing a file with like one * for every 10 or 100 values in a bucket in a text file, if you just want to LOOK at it once and get a sense of the data.