QMcBeaver: QMCStatistic.cpp Source File

00001 //            QMcBeaver
00002 //
00003 //         Constructed by 
00004 //
00005 //     Michael Todd Feldmann 
00006 //              and 
00007 //   David Randall "Chip" Kent IV
00008 //
00009 // Copyright 2000.  All rights reserved.
00010 //
00011 // drkent@users.sourceforge.net mtfeldmann@users.sourceforge.net
00012 
00013 #include "QMCStatistic.h"
00014 
00015 QMCStatistic::QMCStatistic()
00016 {
00017   zeroOut();
00018 
00019 #ifdef PARALLEL
00020   if( !mpiTypeCreated )
00021     {
00022       mpiTypeCreated = true;
00023       buildMpiType();
00024       buildMpiReduce();
00025     }
00026 #endif
00027 }
00028   
00029 
00030 void QMCStatistic::zeroOut()
00031 {
00032   //sum      = 0.0;
00033   //sum2     = 0.0;
00034   //weights  = 0.0;
00035   nsamples = 0;
00036 }
00037 
00038 unsigned long QMCStatistic::getNumberSamples() const
00039 {
00040   return nsamples;
00041 }
00042 
00043 double QMCStatistic::getAverage() const
00044 {
00045   if( nsamples == 0 || weights == 0 ) return 0;
00046   return (sum/weights);
00047 }
00048 
00049 double QMCStatistic::getVariance() const
00050 {
00051   if( nsamples == 0 || weights == 0 ) return 0;
00052   return (sum2/weights-getAverage()*getAverage());
00053 }
00054 
00055 double QMCStatistic::getSkewness() const
00056 {
00057   if( nsamples < 5 || weights == 0 ) return 0;
00058   double skew = sum3/weights;
00059   skew -= 3.0 * getAverage() * sum2/weights;
00060   skew += 2.0 * getAverage() * getAverage() * getAverage();
00061   double var = getVariance();
00062 
00063   if(fabs(skew) > 1e100 || fabs(var) < 1e-50) return 0.0;
00064   return skew / pow(var,1.5);
00065 }
00066 
00067 double QMCStatistic::getKurtosis() const
00068 {
00069   if( nsamples < 5 || weights == 0 ) return 0;
00070   double kurt = sum4/weights;
00071   kurt -= 4.0 * getAverage() * sum3/weights;
00072   kurt += 6.0 * getAverage() * getAverage() * sum2/weights;
00073   kurt -= 3.0 * getAverage() * getAverage() * getAverage() * getAverage();
00074   double var = getVariance();
00075 
00076   if(fabs(kurt) > 1e100 || fabs(var) < 1e-50) return 0.0;
00077   return kurt / (var * var) - 3.0;
00078 }
00079 
00080 double QMCStatistic::getStandardDeviation() const
00081 {
00082   return sqrt(getVariance());
00083 }
00084 
00085 void QMCStatistic::newSample(long double s, long double weight)
00086 {
00087   if(nsamples == 0)
00088     {
00089       weights =         weight;
00090       sum     =       s*weight;
00091       sum2    =     s*s*weight;
00092       sum3    =   s*s*s*weight;
00093       sum4    = s*s*s*s*weight;
00094     } else {
00095       weights +=         weight;
00096       sum     +=       s*weight;
00097       sum2    +=     s*s*weight;
00098       sum3    +=   s*s*s*weight;
00099       sum4    += s*s*s*s*weight;
00100     }
00101   nsamples++;
00102 }
00103 
00104 void QMCStatistic::quickSample(long double s, long double weight)
00105 {
00106   if(nsamples == 0)
00107     {
00108       weights =       1;
00109       sum     =       s;
00110     } else {
00111       weights +=      1;
00112       sum     +=      s;
00113     }
00114   nsamples++;
00115 }
00116 
00117 void QMCStatistic::operator = ( const QMCStatistic & rhs )
00118 {
00119   weights = rhs.weights;
00120   sum = rhs.sum;
00121   sum2 = rhs.sum2;
00122   sum3 = rhs.sum3;
00123   sum4 = rhs.sum4;
00124   nsamples = rhs.nsamples;
00125 }
00126 
00127 QMCStatistic QMCStatistic::operator + (const QMCStatistic &rhs)
00128 {
00129   QMCStatistic result;
00130   if(nsamples == 0)
00131     {
00132       result = rhs;
00133     } else {
00134       result.weights  = weights+rhs.weights;
00135       result.sum      = sum+rhs.sum;
00136       result.sum2     = sum2+rhs.sum2;
00137       result.sum3     = sum3+rhs.sum3;
00138       result.sum4     = sum4+rhs.sum4;
00139       result.nsamples = nsamples+rhs.nsamples;
00140     }
00141   return result;
00142 }
00143 
00144 void QMCStatistic::reWeight(double w)
00145 {
00146   //weights *= w;
00147   sum     *= w;
00148   sum2    *= w;
00149   sum3    *= w;
00150   sum4    *= w;
00151 }
00152 
00153 void QMCStatistic::toXML(ostream& strm)
00154 {
00155   // Open XML
00156   strm << "<QMCStatistic>" << endl;
00157 
00158   if(nsamples == 0)
00159     {
00160       // sum
00161       strm << "<Sum>\n" << 0.0 << "\n</Sum>" << endl;
00162       // sum sq
00163       strm << "<SumSquared>\n" << 0.0 << "\n</SumSquared>" << endl; 
00164       // weights
00165       strm << "<SumWeights>\n" << 0.0 << "\n</SumWeights>" << endl;
00166     } else {
00167       // sum
00168       strm << "<Sum>\n" << double(sum) << "\n</Sum>" << endl;
00169       // sum sq
00170       strm << "<SumSquared>\n" << double(sum2) << "\n</SumSquared>" << endl; 
00171       // weights
00172       strm << "<SumWeights>\n" << double(weights) << "\n</SumWeights>" << endl;
00173     }
00174   
00175   // nsamples
00176   strm << "<NumberOfSamples>\n" << nsamples << "\n</NumberOfSamples>" << endl;
00177 
00178   // Close XML
00179   strm << "</QMCStatistic>" << endl;
00180 }
00181 
00182 bool QMCStatistic::readXML(istream& strm)
00183 {
00184   string temp;
00185 
00186   // Open XML
00187   strm >> temp;
00188   if (temp != "<QMCStatistic>")
00189     return false;
00190 
00191   // sum
00192   strm >> temp;
00193   if (temp != "<Sum>")
00194     return false;
00195   strm >> temp;
00196   sum = atof(temp.c_str());
00197   strm >> temp;
00198   if (temp != "</Sum>")
00199     return false;
00200   
00201   // sum sq
00202   strm >> temp;
00203   if (temp != "<SumSquared>")
00204     return false;
00205   strm >> temp;
00206   sum2 = atof(temp.c_str());
00207   strm >> temp;
00208   if (temp != "</SumSquared>")
00209     return false;
00210 
00211   // weights
00212   strm >> temp;
00213   if (temp != "<SumWeights>")
00214     return false;
00215   strm >> temp;
00216   weights = atof(temp.c_str());
00217   strm >> temp;
00218   if (temp != "</SumWeights>")
00219     return false;
00220 
00221   // nsamples
00222   strm >> temp;
00223   if (temp != "<NumberOfSamples>")
00224     return false;
00225   strm >> temp;
00226   nsamples = atoi(temp.c_str());
00227   strm >> temp;
00228   if (temp != "</NumberOfSamples>")
00229     return false;
00230 
00231   // Close XML
00232   strm >> temp;
00233   if (temp != "</QMCStatistic>")
00234     return false;
00235 
00236   return true;
00237 }
00238 
00239 ostream& operator << (ostream& strm, const QMCStatistic &rhs)
00240 {
00241   strm.precision(12);
00242   strm.width(20);
00243   strm << scientific << rhs.getAverage() << " +/- ";
00244   if( fabs(rhs.getStandardDeviation()) > 1e-300 )
00245     if( log( fabs( rhs.getStandardDeviation() )) > 10.0)
00246       strm << scientific;
00247   strm.precision(12);
00248   strm.width(20);
00249   strm << scientific << rhs.getStandardDeviation() << " (" << rhs.getNumberSamples() << " samples)" << endl;
00250   
00251   //strm << rhs.getAverage() << " +/- " << rhs.getStandardDeviation();
00252   return strm;
00253 }
00254 
00255 
00256 #ifdef PARALLEL
00257 
00258 bool QMCStatistic::mpiTypeCreated = false;
00259 
00260 void QMCStatistic::buildMpiReduce()
00261 {
00262   MPI_Op_create((MPI_User_function*)Reduce_Function,
00263                 true,&MPI_REDUCE);
00264 }
00265 
00266 void QMCStatistic::buildMpiType()
00267 {
00268   QMCStatistic indata;
00269   int num = 6;
00270   int          block_lengths[num];
00271   MPI_Aint     displacements[num];
00272   MPI_Aint     addresses[num+1];
00273   MPI_Datatype typelist[num];
00274 
00275   typelist[0] = MPI_LONG_DOUBLE; //this is the weights 
00276                                  //(roughly the # of samples)
00277   typelist[1] = MPI_LONG_DOUBLE; //this is the sum
00278   typelist[2] = MPI_LONG_DOUBLE; //this is the sum2
00279   typelist[3] = MPI_LONG_DOUBLE; //this is the sum3
00280   typelist[4] = MPI_LONG_DOUBLE; //this is the sum4
00281   typelist[5] = MPI_LONG;        //this is nsamples
00282   
00283   block_lengths[0] = 1;
00284   block_lengths[1] = 1;
00285   block_lengths[2] = 1;
00286   block_lengths[3] = 1;
00287   block_lengths[4] = 1;
00288   block_lengths[5] = 1;
00289   
00290   MPI_Address(&indata, &addresses[0]);
00291   MPI_Address(&(indata.weights), &addresses[1]);
00292   MPI_Address(&(indata.sum), &addresses[2]);
00293   MPI_Address(&(indata.sum2), &addresses[3]);
00294   MPI_Address(&(indata.sum3), &addresses[4]);
00295   MPI_Address(&(indata.sum4), &addresses[5]);
00296   MPI_Address(&(indata.nsamples), &addresses[6]);
00297 
00298   displacements[0] = addresses[1] - addresses[0];
00299   displacements[1] = addresses[2] - addresses[0];
00300   displacements[2] = addresses[3] - addresses[0];
00301   displacements[3] = addresses[4] - addresses[0];
00302   displacements[4] = addresses[5] - addresses[0];
00303   displacements[5] = addresses[6] - addresses[0];
00304   
00305   MPI_Type_struct(num, block_lengths, displacements, typelist, 
00306                   &MPI_TYPE);
00307   MPI_Type_commit(&MPI_TYPE);
00308 }
00309 
00310 
00311 MPI_Datatype QMCStatistic::MPI_TYPE;
00312 
00313 MPI_Op QMCStatistic::MPI_REDUCE;
00314 
00315 void QMCStatistic::Reduce_Function(QMCStatistic *in, QMCStatistic *inout, 
00316                                    int *len, MPI_Datatype *dptr)
00317 {
00318   *inout = *inout + *in;
00319 }
00320 
00321 
00322 #endif
00323 
00324 
00325 
00326 
00327 
00328 
00329 
00330