/*
  Copyright (C) 2007 Steven L. Scott

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
*/

#include <LinAlg/Types.hpp>
#include <Models/ModelTypes.hpp>
#include <Models/Sufstat.hpp>
#include <Models/Policies/ParamPolicy_2.hpp>
#include <Models/Policies/SufstatDataPolicy.hpp>
#include <Models/Policies/PriorPolicy.hpp>
#include <Models/Glm/GlmCoefs.hpp>
#include <Models/MvnBase.hpp>

/*
  Base class that can serve as a prior model for GlmCoefs.  Might
  depend on a scalar sigma, or not.  Might depend on XTX, which might
  or might not change from one iteration to the next.  Probably
  depends on a scalar 'prior sample size' kappa as well.
 */

namespace BOOM{

  class GlmMvnSuf;

  class GlmMvnPriorBase
    : public SufstatDataPolicy<GlmCoefs,GlmMvnSuf>
  {

    // model:  beta | b, V, sigsq, gamma ~ N(b_g, sigsq * V_g/kappa)


  public:
    GlmMvnPriorBase(uint dim);
    virtual double pdf(Ptr<Data>, bool logscale)const=0;
    virtual double pdf(const Ptr<GlmCoefs> beta, bool logscale)const=0;

    virtual double sigsq()const;       // default return is 1.0;
    virtual double kappa()const;       // default return is 1.0;

    virtual Vec mu()const=0;      // not conditional on gamma
    virtual Spd siginv()const=0;  // not conditional on gamma
    // conceptually, siginv is the inverse of sigsq * XTX / kappa

  };


  class GlmMvnSuf
    : public SufstatDetails<GlmCoefs>
  {
  public:
    GlmMvnSuf(uint p=0);
    GlmMvnSuf(const std::vector<Ptr<GlmCoefs> > & d);
    GlmMvnSuf * clone()const;

    virtual void clear();
    void Update(const GlmCoefs &beta);

    Spd center_sumsq(const Vec &b)const;
    const Vec & vnobs()const;  // sum of gamma
    const Spd & GTG()const;    // sum of gamma gamma^T
    const Mat & BTG()const;    // sum of beta gamma^T
    uint nobs()const;
  private:
    mutable Spd bbt_;    // sum of beta beta.transpose()
    mutable Spd ggt_;    // sum of gamma * gamma^T
    Mat bgt_;    // sum of beta gamma.transpose()
    Vec vnobs_;  // sum of gamma
    uint nobs_;  // number of observations

    Vec b, gam;
    mutable bool sym_;
    void make_symmetric()const;
  };


}
