/* ----------------------------------------------------------------------------

Copyright (C) 2005, 2006.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

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

This program 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 General Public License for more details.

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

-----------------------------------------------------------------------------*/

#include "libscl.h"
#include "libgsm.h"
#include "gsmusr.h"
using namespace std;  
using namespace scl;
using namespace libgsm;

namespace gsm {

  two_dim_tinker_toy_sci_mod::two_dim_tinker_toy_sci_mod
    (const realmat* dat_ptr, const std::vector<std::string>& pfvec,
     const std::vector<std::string>& alvec, std::ostream& detail)
  : theta(2,1)
  {
    theta[1]=1.0; theta[2]=2.0;
    vector<string>::const_iterator pfv_ptr = alvec.begin();
    len_sim = atoi((++pfv_ptr)->substr(0,12).c_str());
    detail << starbox("/tt2d sci_mod constructor//");
    detail << "\t len_sim = " << len_sim << '\n';
    detail << "\t theta = " << theta << '\n';
  }
  
  
  bool two_dim_tinker_toy_sci_mod::support(const realmat& parm)
  {
    const REAL a1 = 0.1;
    const REAL a2 = 0.1;
    const REAL b1 = 4.0;
    const REAL b2 = 4.0;
    if (parm[1] < a1 || b1 < parm[1]) return false;
    if (parm[2] < a2 || b2 < parm[2]) return false;
    return true;
  }
  
  den_val two_dim_tinker_toy_sci_mod::prior
    (const realmat& parm, const realmat& func)
  {
    if (parm[2] < parm[1] + 1.0) return den_val(false,-REAL_MAX);
    if (func[7] > 3.0) return den_val(false,-REAL_MAX);
    REAL z = -4.0*pow( (func[7] - 3.0), 2);
    return den_val(true,z);
  }
  
  bool two_dim_tinker_toy_sci_mod::gen_sim(realmat& sim, realmat& func)
  {
    if (theta[1] <= 0.0) return false;
    if (theta[2] <= 0.0) return false;
  
    INT_32BIT seed = 770116;
  
    sim.resize(2,len_sim);
    func.resize(7,1);
    realmat mean(2,1,0.0);
    realmat var(2,2,0.0);
    REAL prod = 0.0;
    
    for (INTEGER i=1; i<=len_sim; ++i) {
      sim(1,i) = theta[1] + theta[1]*scl::unsk(seed);
      sim(2,i) = theta[2] + theta[2]*scl::unsk(seed);
      mean[1] += sim(1,i);
      mean[2] += sim(2,i);
      prod += sim(1,i)*sim(2,i);
    }
  
    mean = mean/len_sim;
    prod = prod/len_sim;
  
    realmat e(2,1);
    for (INTEGER i=1; i<=len_sim; ++i) {
      e[1] = sim(1,i)-mean[1];
      e[2] = sim(2,i)-mean[2];
      var += e*T(e);
    }
  
    var = var/(len_sim-1);
  
  
    func[1] = mean[1];
    func[2] = var(1,1);
    func[3] = func[1]/sqrt(func[2]);
    func[4] = mean[2];
    func[5] = var(2,2);
    func[6] = func[4]/sqrt(func[5]);
    func[7] = prod;
  
    return true;
  }
      
  two_dim_tinker_toy_stat_mod::two_dim_tinker_toy_stat_mod
    (const realmat* dat_ptr, const std::vector<std::string>& pfvec,
     const std::vector<std::string>& alvec, std::ostream& detail)
  : eta(5,1)
  {
    eta[1]=1.0; eta[2]=2.0; eta[3]=1.0; eta[4]=0.0; eta[5]=4.0; 
    vector<string>::const_iterator pfv_ptr = alvec.begin();
    len_sim = atoi((++pfv_ptr)->substr(0,12).c_str());
    detail << starbox("/tt2d stat_mod constructor//");
    detail << "\t len_sim = " << len_sim << '\n';
    detail << "\t eta = " << eta << '\n';
  }
  
  bool two_dim_tinker_toy_stat_mod::support(const realmat& parm)
  {
    if (parm[3] <= 0.0) return false;
    if (parm[5] <= 0.0) return false;
    return true;
  }
  
  den_val two_dim_tinker_toy_stat_mod::prior
    (const realmat& parm, const realmat& func)
  {
    return den_val(true,0.0);
  }
    
  bool two_dim_tinker_toy_stat_mod::gen_sim(realmat& sim, realmat& func)
  {
    realmat mu(2,1);
    realmat sig(2,2);
    realmat R(2,2);
  
    mu[1] = eta[1];
    mu[2] = eta[2];
  
    sig(1,1) = eta[3];
    sig(2,1) = eta[4];
    sig(1,2) = eta[4];
    sig(2,2) = eta[5];
  
    if (sig(2,2) <= 0) return false;
    R(2,2) = sqrt(sig(2,2));
    R(1,2) = sig(1,2)/R(2,2);
    REAL v = sig(1,1)-pow(R(1,2),2);
    if (v <= 0.0) return false;
    R(1,1) = sqrt(v);
    R(2,1) = 0.0;
  
    INTEGER rank;
  
    realmat siginv = invpsd(sig,rank);
  
    if (rank != 2) return false;
  
    INT_32BIT seed = 770116;
  
    sim.resize(2,len_sim);
    func.resize(6,1);
    realmat mean(2,1,0.0);
    realmat var(2,2,0.0);
    
    realmat z(2,1);
    for (INTEGER i=1; i<=len_sim; ++i) {
      z[1] = scl::unsk(seed);
      z[2] = scl::unsk(seed);
      realmat y = mu + R*z;
      sim(1,i) = y[1];
      sim(2,i) = y[2];
      mean += sim("",i);
    }
    
    mean = mean/len_sim;
  
    for (INTEGER i=1; i<=len_sim; ++i) {
      var += (sim("",i)-mean)*T(sim("",i)-mean);
    }
  
    var = var/(len_sim-1);
  
    func[1] = mean[1];
    func[2] = var(1,1);
    func[3] = func[1]/sqrt(func[2]);
    func[4] = mean[2];
    func[5] = var(2,2);
    func[6] = func[4]/sqrt(func[5]);
  
    return true;
  }
  
  den_val two_dim_tinker_toy_stat_mod::loglikelihood()
  {
    const REAL half_log_two_pi = 0.91893853320467266954;
  
    INTEGER n = data_ptr->get_cols();
  
    realmat mu(2,1);
    realmat sig(2,2);
    realmat R(2,2);
  
    mu[1] = eta[1];
    mu[2] = eta[2];
  
    sig(1,1) = eta[3];
    sig(2,1) = eta[4];
    sig(1,2) = eta[4];
    sig(2,2) = eta[5];
  
    realmat siginv = invpsd(sig);
  
    realmat sum(1,1,0.0);
    realmat e(2,1);
    for (INTEGER i=1; i<=n; ++i) { 
      e[1] = (*data_ptr)(1,i) - mu[1];
      e[2] = (*data_ptr)(2,i) - mu[2];
      sum += T(e)*(siginv*e);
    }
  
    REAL q = (-0.5)*sum[1];
  
    REAL log_bot = 2.0*half_log_two_pi; 
    log_bot += 0.5*log( sig(1,1)*sig(2,2)-sig(1,2)*sig(2,1) );
    log_bot *= REAL(n);
  
    return den_val(true, q - log_bot);
  }
  
}
