/*
 * $Id: example0a.cpp,v 1.1 2007/08/15 19:44:55 brook Exp $
 */

/*
 * Copyright (c) 2007 Brook Milligan.
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 */

/*
 * The purpose of this example is to illustrate a traditional solution
 * to a typical likelihood problem, including the pitfalls of using a
 * single type (double) to represent values of very different
 * semantics (probability, likelihood, log(likelihood).  It was
 * created as an example for the library documentation, and included
 * here to ensure that it compiles correctly.  Note that the actual
 * implementations of the classes and functions is meaningless, as
 * they are not relevant to illustrating the structure of the solution
 * encapsulated within the main() function.
 */

#include <cmath>
#include <iostream>
#include <vector>
#include <assert.h>

class GHCN_record { /* ... */ };

typedef std::vector<GHCN_record> GHCN_data_set;

class GlobalClimateModel
{
public:
  virtual ~GlobalClimateModel () { }
  // note:  return probablity (not log(likelihood)) of data
  virtual double probability (const GHCN_record&) const = 0;
};

class GlobalClimateModelFactory
{
public:
  GlobalClimateModel* operator () (const char*) const;
};

// read GHCN data set
std::istream& operator >> (std::istream&, GHCN_data_set&);

int
main ()
{
  GlobalClimateModelFactory climate_model_factory;
  GlobalClimateModel* climate_model (climate_model_factory("GCM"));

  GHCN_data_set data;
  std::cin >> data;
  assert (data.size() == 590543); // XXX - illustrate large size of dataset

				// Note: lnL represents
				// log(likelihood).  The initial
				// condition is required to accumulate
				// the sum of log(likelihood) (rather
				// than the product of likelihood) to
				// avoid underflow with such a large
				// dataset
  double lnL (0);
  for (GHCN_data_set::const_iterator i = data.begin(); i != data.end(); ++i)
    {
      GHCN_record record = *i;
				// Note: p represents the probability
				// (not log(probability)) of the
				// record.
      double p (climate_model->probability(record));
				// Note: remember the following:
				// - likelihood is proportional to probability
				// - product of L <=> sum of log(L)
				// - initial condition for lnL above is
				//   appropriate for adding log(L) here
      lnL += log(p);
    }
				// Note: remember to convert back to
				// the likelihood scale (what happens
				// with underflow? overflow?)
  std::cout << "likelihood:  " << exp(lnL) << std::endl;

  return 0;
}

class GCM : public GlobalClimateModel
{
public:
  GCM () : GlobalClimateModel() { }
  virtual ~GCM () { }
  virtual double probability (const GHCN_record&) const { return 1; }
};

GlobalClimateModel*
GlobalClimateModelFactory::operator () (const char*) const { return new GCM; }

std::istream&
operator >> (std::istream& s, GHCN_data_set& d)
{
  d = GHCN_data_set (590543);
  return s;
}