#include <RcppArmadillo.h>
#include <vector>
#include <map>
#include "c_optim.h"

namespace rstpm2 {

  using namespace Rcpp;
  using namespace arma;

  // typedefs

  typedef bool constraintfn(int, double *, void *);

  // Hadamard element-wise multiplication for the _columns_ of a matrix with a vector

  mat rmult(mat m, vec v) {
    mat out(m);
    out.each_col() %= v;
    return out;
  }
  mat lmult(vec v, mat m) {
    return rmult(m,v);
  }

  // element-wise multiplication and division for NumericVector and arma::vec

  NumericVector ew_mult(NumericVector nv, vec v) {
    return as<NumericVector>(wrap(as<vec>(wrap(nv)) % v));
  }

  NumericVector ew_div(NumericVector nv, vec v) {
    return as<NumericVector>(wrap(as<vec>(wrap(nv)) / v));
  }

  void Rprint(NumericMatrix m) {
    for (int i=0; i<m.nrow(); ++i) {
      for (int j=0; j<m.ncol(); ++j) 
	Rprintf("%f ", m(i,j));
      Rprintf("\n");
    }
  }

  void Rprint(NumericVector v) {
    for (int i=0; i<v.size(); ++i) 
      Rprintf("%f ", v(i));
    Rprintf("\n");
  }

  void Rprint(vec v) {
    for (size_t i=0; i<v.size(); ++i) 
      Rprintf("%f ", v(i));
    Rprintf("\n");
  }
  void Rprint(rowvec v) {
    for (size_t i=0; i<v.size(); ++i) 
      Rprintf("%f ", v(i));
    Rprintf("\n");
  }

  void Rprint(mat m) {
    for (size_t i=0; i<m.n_rows; ++i) {
      for (size_t j=0; j<m.n_cols; ++j) 
	Rprintf("%f ", m(i,j));
      Rprintf("\n");
    }
  }

  vec pnorm01(vec x) {
    vec out(x.size());
    for (size_t i=0; i<x.size(); ++i)
      out(i) = R::pnorm(x(i),0.0,1.0,1,0);
    return out;
    //return as<vec>(wrap(pnorm(as<NumericVector>(wrap<vec>(x)),0.0,1.0)));
  }

  vec qnorm01(vec x) {
    vec out(x.size());
    for (size_t i=0; i<x.size(); ++i)
      out(i) = R::qnorm(x(i),0.0,1.0,1,0);
    return out;
  }

  vec dnorm01(vec x) {
    vec out(x.size());
    for (size_t i=0; i<x.size(); ++i)
      out(i) = R::dnorm(x(i),0.0,1.0,0);
    return out;
  }

  vec logit(vec p) {
    return log(p/(1-p));
  }

  vec expit(vec x) {
    return 1/(1+exp(-x));
  }

  RcppExport SEXP tapplySum(SEXP s_y, SEXP s_group) {
    NumericVector y(s_y);
    NumericVector group(s_group);
    NumericVector::iterator iy, igroup;
    std::map<double,double> out;
    for (iy = y.begin(), igroup = group.begin(); iy != y.end(); ++iy, ++igroup) { 
	out[*igroup] += *iy;
    }
    return wrap(out);
  }

  // Q: do we need parscale, trace, reltol and criterion in these structs?

  class stpm2 {
  public:
    stpm2(SEXP sexp) {
      List list = as<List>(sexp);
      init = as<NumericVector>(list["init"]);
      X = as<mat>(list["X"]); 
      XD = as<mat>(list["XD"]); 
      bhazard = as<vec>(list["bhazard"]);
      wt = as<vec>(list["wt"]);
      event = as<vec>(list["event"]);
      time = as<vec>(list["time"]);
      delayed = as<int>(list["delayed"]);
      X0.set_size(1,1); X0.fill(0.0);
      wt0.set_size(1); wt0.fill(0.0);
      if (delayed == 1) {
	X0 = as<mat>(list["X0"]); // optional
	XD0 = as<mat>(list["XD0"]); // optional
	wt0 = as<vec>(list["wt0"]); // optional
      }
      parscale = as<vec>(list["parscale"]);
      reltol = as<double>(list["reltol"]);
      kappa = as<double>(list["kappa"]);
      trace = as<int>(list["trace"]);
    }
    // defaults are for the PH model
    virtual vec link(vec S) { return log(-log(S)); }
    virtual vec ilink(vec x) { return exp(-exp(x)); }
    virtual vec h(vec eta, vec etaD) { return etaD % exp(eta); }
    virtual vec H(vec eta, vec etaD) { return exp(eta); }
    virtual mat gradH(vec eta, vec etaD, mat X, mat XD) { return rmult(X,exp(eta)); }
    virtual mat gradh(vec eta, vec etaD, mat X, mat XD) { 
      return rmult(XD, exp(eta)) + rmult(X, etaD % exp(eta));
    }
    virtual double negll(NumericVector beta);
    virtual void grad_negll(NumericVector beta);
    NumericVector init;
    mat X, XD, X0, XD0; 
    vec bhazard,wt,wt0,event,time,parscale;
    double reltol, kappa;
    int delayed, trace;
  };

  double stpm2::negll(NumericVector beta) {
    int n = init.size();
    vec vbeta(&beta[0],n);
    vbeta = vbeta % parscale;
    vec eta = X * vbeta;
    vec etaD = XD * vbeta;
    vec h = this->h(eta,etaD) + bhazard;
    vec H = this->H(eta,etaD) + bhazard;
    double constraint = kappa/2.0 * sum(h % h % (h<0)); // sum(h^2 | h<0)
    vec eps = h*0.0 + 1.0e-16; 
    h = max(h,eps);
    double ll = sum(wt % event % log(h)) - sum(wt % H) - constraint;
    if (delayed == 1) {
      vec eta0 = X0 * vbeta;
      vec etaD0 = XD0 * vbeta;
      mat H0 = this->H(eta0,etaD0);
      ll += sum(wt0 % H0);
    }
    return -ll;  
  }

  void stpm2::grad_negll(NumericVector beta) {
    int n = beta.size();
    vec vbeta = as<vec>(wrap(beta));
    vbeta = vbeta % parscale;
    vec eta = X * vbeta;
    vec etaD = XD * vbeta;
    vec h = this->h(eta,etaD) + bhazard;
    // vec H = exp(eta);
    mat gradH = this->gradH(eta,etaD,X,XD);
    mat gradh = this->gradh(eta,etaD,X,XD);
    mat Xgrad = -gradH + rmult(gradh, event / h);
    mat Xconstraint = - kappa*rmult(gradh,h);
    vec eps = h*0.0 + 1.0e-16; // hack
    //Xgrad.rows(h<=eps) = Xconstraint.rows(h<=eps);
    Xgrad = rmult(Xgrad,wt);
    rowvec vgr = sum(Xgrad,0);
    if (delayed == 1) {
      vec eta0 = X0 * vbeta;
      vec etaD0 = XD0 * vbeta;
      mat gradH0 = this->gradH(eta0, etaD0, X0, XD0); 
      vgr += sum(rmult(gradH0, wt0),0);
    }
    for (int i = 0; i<n; ++i) {
      // gr[i] = -vgr[i]*parscale[i];
    }
  }


  class stpm2PO : public stpm2 {
  public:
    stpm2PO(SEXP sexp) : stpm2(sexp) { }
    vec link(vec S) { return -logit(S); }
    vec ilink(vec x) { return expit(-x); }
    vec h(vec eta, vec etaD) { return etaD % exp(eta) % expit(-eta); }
    vec H(vec eta, vec etaD) { return log(1+exp(eta)); }
    mat gradH(vec eta, vec etaD, mat X, mat XD) { 
      return rmult(X,exp(eta) % expit(-eta));
    }
    mat gradh(vec eta, vec etaD, mat X, mat XD) { 
      return rmult(X, etaD % exp(eta) % expit(-eta)) - 
	rmult(X, exp(2*eta) % etaD % expit(-eta) % expit(-eta)) +
	rmult(XD, exp(eta) % expit(-eta));
    }
  };

  class stpm2Probit : public stpm2 {
  public:
    stpm2Probit(SEXP sexp) : stpm2(sexp) { }
    vec link(vec S) { return -qnorm01(S); }
    vec ilink(vec eta) { return pnorm01(-eta); }
    vec H(vec eta, vec etaD) { return -log(pnorm01(-eta)); }
    vec h(vec eta, vec etaD) { return etaD % dnorm01(eta) / pnorm01(-eta); }
    mat gradH(vec eta, vec etaD, mat X, mat XD) { 
      return rmult(X, dnorm01(eta) / pnorm01(-eta));
    }
    mat gradh(vec eta, vec etaD, mat X, mat XD) { 
      return rmult(X, -eta % dnorm01(eta) % etaD / pnorm01(-eta)) +
	rmult(X, dnorm01(eta) % dnorm01(eta) / pnorm01(-eta) / pnorm01(-eta) % etaD) +
	rmult(XD,dnorm01(eta) / pnorm01(-eta));
    }
  };
  
  class stpm2AH : public stpm2 {
  public:
    stpm2AH(SEXP sexp) : stpm2(sexp) { }
    vec link(vec S) { return -log(S); }
    vec ilink(vec x) { return exp(-x); }
    vec h(vec eta, vec etaD) { return etaD; }
    vec H(vec eta, vec etaD) { return eta; }
    mat gradh(vec eta, vec etaD, mat X, mat XD) { return XD; }
    mat gradH(vec eta, vec etaD, mat X, mat XD) { return X;  }
  };

  struct SmoothLogH {
    int first_para, last_para;
    mat S;
  };

  struct SmoothHaz {
    mat X0, X1, X2, X3;
    vec w;
    double lambda;
  };

  template<class Smooth>
  std::vector<Smooth> read_smoothers(List lsmooth) {}

  template<>
  std::vector<SmoothLogH> read_smoothers(List lsmooth) {
    std::vector<SmoothLogH> smooth;
    for(int i=0; i<lsmooth.size(); ++i) {
      List lsmoothi = as<List>(lsmooth[i]);
      List lsmoothiS = as<List>(lsmoothi["S"]);
      SmoothLogH smoothi =  {
	as<int>(lsmoothi["first.para"]) - 1, 
	as<int>(lsmoothi["last.para"]) - 1, 
	as<mat>(lsmoothiS[0])
      };
      smooth.push_back(smoothi);
    }
    return smooth;
  }

  template<>
  std::vector<SmoothHaz> read_smoothers(List lsmooth) {
    std::vector<SmoothHaz> smooth;
    for(int i=0; i<lsmooth.size(); ++i) {
      List lsmoothi = as<List>(lsmooth[i]);
      SmoothHaz smoothi =  {
	as<mat>(lsmoothi["X0"]), 
	as<mat>(lsmoothi["X1"]), 
	as<mat>(lsmoothi["X2"]), 
	as<mat>(lsmoothi["X3"]), 
	as<vec>(lsmoothi["w"]), 
	as<double>(lsmoothi["lambda"])
      };
      smooth.push_back(smoothi);
    }
    return smooth;
  }

  template<class Smooth, class Stpm2 = stpm2>
  class pstpm2 : public Stpm2 {
  public:
    pstpm2(SEXP sexp) : Stpm2(sexp) {
      List list = as<List>(sexp);
      List lsmooth = as<List>(list["smooth"]);
      sp = as<NumericVector>(list["sp"]);
      reltol_search = as<double>(list["reltol_search"]);
      // reltol = as<double>(list["reltol"]);
      alpha = as<double>(list["alpha"]);
      criterion = as<int>(list["criterion"]);
      smooth = read_smoothers<Smooth>(lsmooth);
    }
    NumericVector sp;
    double alpha, reltol_search; 
    int criterion;
    std::vector<Smooth> smooth;
  };

  /** 
      Extension of stpm2 and pstpm2 to include shared frailties with a cluster variable
  **/
  template<class Stpm2>
  class stpm2Shared : public Stpm2 {
  public:
    stpm2Shared(SEXP sexp) : Stpm2(sexp) {
      List list = as<List>(sexp);
      IntegerVector cluster = as<IntegerVector>(list["cluster"]);
      // wragged array indexed by a map of vectors
      for (int i=0; i<cluster.size(); ++i) {
	clusters[cluster[i]].push_back(i);
      }
    }
    std::map<int,std::vector<int> > clusters;
  };


  /**
     Extension to the BFGS class

     optim() and calc_hessian() assume that data->parscale is an array/vector/NumericVector
     optimWithConstraint() further assumes that data->kappa is a double with an initial value
   */
  template<class Data>
  class BFGS2 : public BFGS {
  public:
    void optim(optimfn fn, optimgr gr, NumericVector init, void * ex,
	       bool apply_parscale = true) {
      Data * data = (Data *) ex;
      n = init.size();
      if (apply_parscale) for (int i = 0; i<n; ++i) init[i] /= data->parscale[i];
      BFGS::optim(fn,gr,init,ex);
      if (apply_parscale) for (int i = 0; i<n; ++i) coef[i] *= data->parscale[i];
      hessian = calc_hessian(gr, ex);
    }
    void optimWithConstraint(optimfn fn, optimgr gr, NumericVector init, void * ex, constraintfn constraint,
			     bool apply_parscale = true) {
      Data * data = (Data *) ex;
      n = init.size();
      if (apply_parscale) for (int i = 0; i<n; ++i) init[i] /= data->parscale[i];
      bool satisfied;
      do {
	BFGS::optim(fn,gr,init,ex);
	satisfied = constraint(n,&coef[0],ex);
	if (!satisfied) data->kappa *= 2.0;   
      } while ((!satisfied)&& data->kappa < 1.0e3);      
      if (apply_parscale) for (int i = 0; i<n; ++i) coef[i] *= data->parscale[i];
      hessian = calc_hessian(gr, ex);
   // Rprintf("kappa=%f\n",data->kappa);
    }
    NumericMatrix calc_hessian(optimgr gr, void * ex) {
      Data * data = (Data *) ex;
      vec parscale(n);
      for (int i=0; i<n; ++i) {
	parscale[i] = data->parscale[i];
	data->parscale[i] = 1.0;
      }
      NumericMatrix hessian = BFGS::calc_hessian(gr,ex);
      for (int i=0; i<n; ++i) data->parscale[i] = parscale[i];
      return hessian;
    }
  };

  template<class Data>
  double fminfn(int n, double * beta, void *ex) {
    Data * data = (Data *) ex;
    vec vbeta(beta,n);
    vbeta = vbeta % data->parscale;
    vec eta = data->X * vbeta;
    vec etaD = data->XD * vbeta;
    vec h = data->h(eta,etaD) + data->bhazard;
    vec H = data->H(eta,etaD);
    double constraint = data->kappa/2.0 * sum(h % h % (h<0)) + 
      data->kappa/2.0 * sum((H/data->time) % (H/data->time) % (H<0));
    vec eps = h*0.0 + 1.0e-16; 
    h = max(h,eps);
    H = max(H,eps);
    double ll = sum(data->wt % data->event % log(h)) - sum(data->wt % H) - constraint;
    if (data->delayed == 1) {
      vec eta0 = data->X0 * vbeta;
      vec etaD0 = data->XD0 * vbeta;
      mat H0 = data->H(eta0,etaD0);
      ll += sum(data->wt0 % H0);
    }
    return -ll;  
  }

  /** 
      Objective function for a Gamma shared frailty
      Assumes that weights == 1.
      The implementation is incomplete, requiring either the use of optimisation without
      gradients or the calculation of the gradients.
  **/
  template<class Data>
  double fminfnShared(int n, double * beta, void *ex) {
    typedef std::vector<int> Index;
    typedef std::map<int,Index> IndexMap;
    Data * data = (Data *) ex;
    vec vbeta(beta,n-1); // theta is the last parameter in beta
    double theta = beta[n-1];
    vbeta = vbeta % data->parscale;
    vec eta = data->X * vbeta;
    vec etaD = data->XD * vbeta;
    vec h = data->h(eta,etaD) + data->bhazard;
    vec H = data->H(eta,etaD) + data->bhazard;
    double constraint = data->kappa/2.0 * sum(h % h % (h<0)); // sum(h^2 | h<0)
    vec eps = h*0.0 + 1.0e-16; 
    h = max(h,eps);
    mat H0;
    //double ll = sum(data->wt % data->event % log(h)) - sum(data->wt % H) - constraint;
    if (data->delayed == 1) {
      vec eta0 = data->X0 * vbeta;
      vec etaD0 = data->XD0 * vbeta;
      H0 = data->H(eta0,etaD0);
      //ll += sum(data->wt0 % H0);
    }
    double ll;
    for (IndexMap::iterator it=data->clusters.begin(); it!=data->clusters.end(); ++it) {
      int mi=0;
      double sumH = 0.0, sumHenter = 0.0;
      for (Index::iterator j=it->second.begin(); j!=it->second.end(); ++j) {
	if (data->event(*j)==1) {
	  mi++;
	  ll += log(h(*j));
	}
	sumH += H(*j);
	if (data->delayed == 1)
	  sumHenter += H0(*j);
      }
      ll -= (1.0/theta+mi)*log(1.0+theta*sumH);
      if (data->delayed == 1)
	ll += 1.0/theta*log(1.0+theta*sumHenter);
      if (mi>0) {
	int k=1;
	for (Index::iterator j=it->second.begin(); j!=it->second.end(); ++j) {
	  if (data->event(*j)==1) {
	    ll += log(1.0+theta*(mi-k));
	    k++;
	  }
	}
      }
    }
    ll -= constraint;
    return -ll;  
  }

  /** 
      Utility function for calculating a gamma frailty log-likelihood within R code 
   **/
  RcppExport SEXP llgammafrailty(SEXP _theta, SEXP _h, SEXP _H, SEXP _d,
				    SEXP _cluster) {
    double theta = as<double>(_theta);
    NumericVector h = as<NumericVector>(_h);
    NumericVector H = as<NumericVector>(_H);
    IntegerVector d = as<IntegerVector>(_d);
    IntegerVector cluster = as<IntegerVector>(_cluster);
    double ll = 0.0;
    // wragged array indexed by a map of vectors
    typedef std::vector<int> Index;
    typedef std::map<int,Index> IndexMap;
    IndexMap clusters;
    for (int i=0; i<cluster.size(); ++i) {
      clusters[cluster[i]].push_back(i);
    }
    for (IndexMap::iterator it=clusters.begin(); it!=clusters.end(); ++it) {
      int mi=0;
      double sumH = 0.0;
      for (Index::iterator j=it->second.begin(); j!=it->second.end(); ++j) {
	if (d[*j]==1) {
	  mi++;
	  ll += log(h[*j]);
	}
	sumH += H[*j];
	// sumHenter += Henter[*j];
      }
      ll -= (1.0/theta+mi)*log(1.0+theta*sumH);
      // ll += 1.0/theta*log(1.0+theta*sumHenter);
      if (mi>0) {
	int k=1;
	for (Index::iterator j=it->second.begin(); j!=it->second.end(); ++j) {
	  if (d[*j]==1) {
	    ll += log(1.0+theta*(mi-k));
	    k++;
	  }
	}
      }
    }
    return wrap(ll);
  }

  /** 
      Utility function for calculating a gamma frailty log-likelihood within R code 
   **/
  RcppExport SEXP llgammafrailtydelayed(SEXP _theta, SEXP _h, SEXP _H, SEXP _H0, SEXP _d,
				    SEXP _cluster) {
    double theta = as<double>(_theta);
    NumericVector h = as<NumericVector>(_h);
    NumericVector H = as<NumericVector>(_H);
    NumericVector H0 = as<NumericVector>(_H0); // assumes that _H0 is zero for data not left truncated
    IntegerVector d = as<IntegerVector>(_d);
    IntegerVector cluster = as<IntegerVector>(_cluster);
    double ll = 0.0;
    // wragged array indexed by a map of vectors
    typedef std::vector<int> Index;
    typedef std::map<int,Index> IndexMap;
    IndexMap clusters;
    for (int i=0; i<cluster.size(); ++i) {
      clusters[cluster[i]].push_back(i);
    }
    for (IndexMap::iterator it=clusters.begin(); it!=clusters.end(); ++it) {
      int mi=0;
      double sumH = 0.0, sumHenter = 0.0;
      for (Index::iterator j=it->second.begin(); j!=it->second.end(); ++j) {
	if (d[*j]==1) {
	  mi++;
	  ll += log(h[*j]);
	}
	sumH += H[*j];
	sumHenter += H0[*j];
      }
      ll -= (1.0/theta+mi)*log(1.0+theta*sumH);
      ll += 1.0/theta*log(1.0+theta*sumHenter);
      if (mi>0) {
	int k=1;
	for (Index::iterator j=it->second.begin(); j!=it->second.end(); ++j) {
	  if (d[*j]==1) {
	    ll += log(1.0+theta*(mi-k));
	    k++;
	  }
	}
      }
    }
    return wrap(ll);
  }


  struct CureModel {
    int n0, n1, n2;
    mat Xshape, Xscale, Xcure;
    vec time, status;
  };

  double fminfn_cureModel(int n, double * beta, void *ex) {
    double ll = 0.0;
    CureModel * data = (CureModel *) ex;
    vec vbeta(beta,n);
    vec shape = exp(data->Xshape * vbeta(span(0,data->n0-1)));
    vec scale = exp(data->Xscale * vbeta(span(data->n0,data->n1-1)));
    vec cure = 1.0/(1.0+exp(-data->Xcure * vbeta(span(data->n1,data->n2-1))));
    for (size_t i=0; i < data->time.size(); ++i) {
      ll += data->status(i)==1.0 ?
	log(1.0-cure(i)) + R::dweibull(data->time(i),shape(i),scale(i),1) :
	log(cure(i)+(1.0-cure(i)) * R::pweibull(data->time(i),shape(i),scale(i),0,0));
    }
    R_CheckUserInterrupt();  /* be polite -- did the user hit ctrl-C? */
    return -ll;
  }

RcppExport SEXP fitCureModel(SEXP stime, SEXP sstatus, SEXP sXshape,
			  SEXP sXscale, SEXP sXcure, SEXP sbeta) {
    mat Xshape = as<mat>(sXshape);
    mat Xscale = as<mat>(sXscale);
    mat Xcure = as<mat>(sXcure);
    vec time = as<vec>(stime);
    vec status = as<vec>(sstatus);
    NumericVector init = as<NumericVector>(sbeta);
    int n0=Xshape.n_cols;
    int n1=n0+Xscale.n_cols;
    int n2=n1+Xcure.n_cols;
    CureModel data = {n0,n1,n2,Xshape,Xscale,Xcure,time,status};
    
    NelderMead nm;
    nm.reltol = 1.0e-8;
    nm.maxit = 1000;
    nm.optim(& fminfn_cureModel, init, (void *) &data);

    for (int i = 0; i<nm.coef.size(); ++i)
      init[i] = nm.coef[i]; // clone
    
    return List::create(_("Fmin")=nm.Fmin, 
			_("coef")=wrap(init),
			_("fail")=nm.fail,
			_("hessian")=wrap(nm.hessian));
  }


  template<class Data>
  bool fminfn_constraint(int n, double * beta, void *ex) {
    Data * data = (Data *) ex;
    vec vbeta(beta,n);
    vbeta = vbeta % data->parscale;
    vec eta = data->X * vbeta;
    vec etaD = data->XD * vbeta;
    vec h = data->h(eta, etaD) + data->bhazard;
    vec H = data->H(eta, etaD);
    return all((h>0) % (H>0));
  }

  template<class Data>
  void fminfn_nlm(int n, double * beta, double * f, void *ex) {
    *f = fminfn<Data>(n, beta, ex);
  };


  // log H penalty - using penalty matrices from mgcv
  template<class Stpm2>
  double pfminfn_SmoothLogH(int n, double * beta, void *ex) {
    typedef SmoothLogH Smooth;
    typedef pstpm2<Smooth,Stpm2> Data;
    double ll = -fminfn<Data>(n,beta,ex);
    Data * data = (Data *) ex;
    vec vbeta(beta,n);
    vbeta = vbeta % data->parscale;
    for (size_t i=0; i < data->smooth.size(); ++i) {
      Smooth smoothi = data->smooth[i];
      ll -= (data->sp)[i]/2 * 
	dot(vbeta.subvec(smoothi.first_para,smoothi.last_para),
	    smoothi.S * vbeta.subvec(smoothi.first_para,smoothi.last_para));
    }
    return -ll;  
  }

  // hazard penalty
  template<class Stpm2>
  double pfminfn_SmoothHaz(int n, double * beta, void *ex) {
    typedef SmoothHaz Smooth;
    typedef pstpm2<Smooth,Stpm2> Data;
    double ll = -fminfn<Data>(n,beta,ex);
    Data * data = (Data *) ex;
    vec vbeta(beta,n);
    vbeta = vbeta % data->parscale;
    for (size_t i=0; i < data->smooth.size(); ++i) {
      Smooth obj = data->smooth[i];
      vec s0 = obj.X0 * vbeta;
      vec s1 = obj.X1 * vbeta;
      vec s2 = obj.X2 * vbeta;
      vec s3 = obj.X3 * vbeta;
      vec h2 = exp(s0) % (s3 + 3*s1%s2 + s1%s1%s1);
      ll -= (data->sp)[i]/2 * obj.lambda * sum(obj.w % h2 % h2);
    }
    return -ll;  
  };

  template<class Data>
  void grfn(int n, double * beta, double * gr, void *ex) {
    Data * data = (Data *) ex;
    vec vbeta(beta,n);
    // if (data->trace) Rprint(vbeta);
    vbeta = vbeta % data->parscale;
    vec eta = data->X * vbeta;
    vec etaD = data->XD * vbeta;
    vec h = data->h(eta,etaD) + data->bhazard;
    vec H = data->H(eta,etaD);
    vec one = ones(h.size());
    vec eps = h*0.0 + 1.0e-16; 
    mat gradH = data->gradH(eta,etaD,data->X,data->XD); 
    mat gradh = data->gradh(eta,etaD,data->X,data->XD); 
    mat Xgrad = -rmult(gradH, one % (H>eps)) + rmult(gradh, data->event / h % (h>eps));
    mat Xconstraint = data->kappa * rmult(gradh,h%(h<eps)) + 
      data->kappa * rmult(gradH,H / data->time / data->time % (H<eps));
    rowvec dconstraint = sum(Xconstraint,0);
    // if (data->trace) Rprint(dconstraint);
    Xgrad = rmult(Xgrad,data->wt);
    rowvec vgr = sum(Xgrad,0);
    if (data->delayed == 1) {
      vec eta0 = data->X0 * vbeta;
      vec etaD0 = data->XD0 * vbeta;
      mat gradH0 = data->gradH(eta0, etaD0, data->X0, data->XD0); 
      vgr += sum(rmult(gradH0, data->wt0),0);
    }
    for (int i = 0; i<n; ++i) {
      gr[i] = -vgr[i]*data->parscale[i] + dconstraint[i];
    }
  }

  template<class Stpm2>
  void pgrfn_SmoothLogH(int n, double * beta, double * gr, void *ex) {
    typedef SmoothLogH Smooth;
    typedef pstpm2<Smooth,Stpm2> Data;
    grfn<Data>(n, beta, gr, ex);
    Data * data = (Data *) ex;
    vec vbeta(beta,n);
    vbeta = vbeta % data->parscale;
    vec eta = data->X * vbeta;
    vec etaD = data->XD * vbeta;
    vec h = data->h(eta,etaD) + data->bhazard;
    vec H = data->H(eta,etaD);
    mat gradh = data->gradh(eta,etaD,data->X,data->XD);
    mat gradH = data->gradH(eta,etaD,data->X,data->XD);
    rowvec vgr(n, fill::zeros);
    for (size_t i=0; i < data->smooth.size(); ++i) {
      SmoothLogH smoothi = data->smooth[i];
      vgr.subvec(smoothi.first_para,smoothi.last_para) += 
	(data->sp)[i] * (smoothi.S * vbeta.subvec(smoothi.first_para,smoothi.last_para)).t();
    }
    for (int i = 0; i<n; ++i) gr[i] += vgr[i]*data->parscale[i];
  }

  template<class Stpm2>
  void pgrfn_SmoothHaz(int n, double * beta, double * gr, void *ex) {
    typedef SmoothHaz Smooth;
    typedef pstpm2<Smooth,Stpm2> Data;
    grfn<Data>(n, beta, gr, ex);
    Data * data = (Data *) ex;
    vec vbeta(beta,n);
    vbeta = vbeta % data->parscale;
    rowvec vgr(n, fill::zeros);
    for (size_t i=0; i < data->smooth.size(); ++i) {
      Smooth obj = data->smooth[i];
      vec s0 = obj.X0 * vbeta;
      vec s1 = obj.X1 * vbeta;
      vec s2 = obj.X2 * vbeta;
      vec s3 = obj.X3 * vbeta;
      vec h2 = exp(s0) % (s3 + 3*s1%s2 + s1%s1%s1);
      mat dh2sq_dbeta = 
	2*lmult(h2 % exp(s0),
		obj.X3+3*(rmult(obj.X1,s2)+rmult(obj.X2,s1))+3*lmult(s1%s1,obj.X1))+
	2*lmult(h2 % h2,obj.X0);
      vgr += (data->sp)[i]*obj.lambda*sum(lmult(obj.w, dh2sq_dbeta),0);
    }
    for (int i = 0; i<n; ++i) gr[i] += vgr[i] * data->parscale[i];
  }

  // oh, for partial template specialisation for functions...

  template<class Smooth,class Stpm2>
  double pfminfn(int n, double * beta, void *ex) { return 0.0; }

  template<>
  double pfminfn<SmoothLogH,stpm2>(int n, double * beta, void *ex) {
    return pfminfn_SmoothLogH<stpm2>(n, beta, ex); }
  template<>
  double pfminfn<SmoothLogH,stpm2PO>(int n, double * beta, void *ex) {
    return pfminfn_SmoothLogH<stpm2PO>(n, beta, ex); }
  template<>
  double pfminfn<SmoothLogH,stpm2Probit>(int n, double * beta, void *ex) {
    return pfminfn_SmoothLogH<stpm2Probit>(n, beta, ex); }
    template<>
  double pfminfn<SmoothLogH,stpm2AH>(int n, double * beta, void *ex) {
    return pfminfn_SmoothLogH<stpm2AH>(n, beta, ex); }

  template<>
  double pfminfn<SmoothHaz,stpm2>(int n, double * beta, void *ex) {
    return pfminfn_SmoothHaz<stpm2>(n, beta, ex); }
  template<>
  double pfminfn<SmoothHaz,stpm2PO>(int n, double * beta, void *ex) {
    return pfminfn_SmoothHaz<stpm2PO>(n, beta, ex); }
  template<>
  double pfminfn<SmoothHaz,stpm2Probit>(int n, double * beta, void *ex) {
    return pfminfn_SmoothHaz<stpm2Probit>(n, beta, ex); }

  template<class Smooth, class Stpm2>
  void pgrfn(int n, double * beta, double * gr, void *ex) { }

  template<>
  void pgrfn<SmoothLogH,stpm2>(int n, double * beta, double * gr, void *ex) {
    pgrfn_SmoothLogH<stpm2>(n, beta, gr, ex); }
  template<>
  void pgrfn<SmoothLogH,stpm2PO>(int n, double * beta, double * gr, void *ex) {
    pgrfn_SmoothLogH<stpm2PO>(n, beta, gr, ex); }
  template<>
  void pgrfn<SmoothLogH,stpm2Probit>(int n, double * beta, double * gr, void *ex) {
    pgrfn_SmoothLogH<stpm2Probit>(n, beta, gr, ex);}
    template<>
  void pgrfn<SmoothLogH,stpm2AH>(int n, double * beta, double * gr, void *ex) {
    pgrfn_SmoothLogH<stpm2AH>(n, beta, gr, ex); }

  template<>
  void pgrfn<SmoothHaz,stpm2>(int n, double * beta, double * gr, void *ex) {
    pgrfn_SmoothHaz<stpm2>(n, beta, gr, ex); }
  template<>
  void pgrfn<SmoothHaz,stpm2PO>(int n, double * beta, double * gr, void *ex) {
    pgrfn_SmoothHaz<stpm2PO>(n, beta, gr, ex); }
  template<>
  void pgrfn<SmoothHaz,stpm2Probit>(int n, double * beta, double * gr, void *ex) {
    pgrfn_SmoothHaz<stpm2Probit>(n, beta, gr, ex); }


  template<class Stpm2>
  SEXP optim_stpm2(SEXP args) {

    Stpm2 data(args);

    BFGS2<Stpm2> bfgs;
    bfgs.reltol = data.reltol;
    bfgs.optimWithConstraint(fminfn<Stpm2>, grfn<Stpm2>, data.init, (void *) &data, fminfn_constraint<Stpm2>);

    return List::create(_("fail")=bfgs.fail,
			_("coef")=wrap(bfgs.coef),
			_("hessian")=wrap(bfgs.hessian));
  }

  RcppExport SEXP optim_stpm2_ph(SEXP args) {
    return optim_stpm2<stpm2>(args); }
  RcppExport SEXP optim_stpm2_po(SEXP args) {
    return optim_stpm2<stpm2PO>(args); }
  RcppExport SEXP optim_stpm2_probit(SEXP args) {
    return optim_stpm2<stpm2Probit>(args); }
  RcppExport SEXP optim_stpm2_ah(SEXP args) {
    return optim_stpm2<stpm2AH>(args); }


  RcppExport SEXP optim_stpm2_nlm(SEXP args) {

    stpm2 data(args);

    data.init = ew_div(data.init,data.parscale);
    int n = data.init.size();
    
    Nlm nlm;
    nlm.gradtl = nlm.steptl = data.reltol;
    //nlm.method=2; nlm.stepmx=0.0;
    bool satisfied;
    do {
      nlm.optim(& fminfn_nlm<stpm2>, & grfn<stpm2>, data.init, (void *) &data);
      satisfied = fminfn_constraint<stpm2>(n,&nlm.coef[0],(void *) &data);
      if (!satisfied) data.kappa *= 2.0;
    } while (!satisfied && data.kappa<1.0e5);

    nlm.coef = ew_mult(nlm.coef, data.parscale);

    return List::create(_("itrmcd")=wrap(nlm.itrmcd),
			_("niter")=wrap(nlm.itncnt),
			_("coef")=wrap(nlm.coef),
			_("hessian")=wrap(nlm.hessian));
  }



  template<class Smooth, class Stpm2>
  double pstpm2_step_first(double logsp, void * ex) {
    typedef pstpm2<Smooth,Stpm2> Data;
    Data * data = (Data *) ex;

    data->sp[0] = exp(logsp);

    BFGS2<Data> bfgs;
    bfgs.reltol = data->reltol;

    bfgs.optimWithConstraint(pfminfn<Smooth,Stpm2>, pgrfn<Smooth,Stpm2>, data->init, ex, fminfn_constraint<Data>, false); // do not apply parscale b4/after

    NumericMatrix hessian0 = bfgs.calc_hessian(grfn<Data>, ex);

    if (data->trace > 0)  {
      Rprintf("Debug on trace calculation. Coef:\n");
      Rprint(bfgs.coef);
      Rprintf("Hessian0:\n");
      Rprint(hessian0);
      Rprintf("Hessian:\n");
      Rprint(bfgs.hessian);
    }

    double edf = trace(solve(as<mat>(bfgs.hessian),as<mat>(hessian0)));
    double negll = bfgs.calc_objective(fminfn<Data>,ex);
    double gcv =  negll + data->alpha*edf;
    double bic =  negll + data->alpha*edf*log(sum(data->event));
    data->init = bfgs.coef;
    if (data->trace > 0)
      Rprintf("sp=%f\tedf=%f\tnegll=%f\tgcv=%f\tbic=%f\talpha=%f\n",data->sp[0],edf,negll,gcv,bic,data->alpha);

    return data->criterion==1 ? gcv : bic;
  }    

  template<class Smooth, class Stpm2>
  double pstpm2_step_multivariate(int n, double * logsp, void * ex) {
    typedef pstpm2<Smooth,Stpm2> Data;
    Data * data = (Data *) ex;

    double lsp = -9.0, usp = 9.0;

    for (int i=0; i < data->sp.size(); ++i)
      //data->sp[i] = exp(logsp[i]);
      data->sp[i] = exp(bound(logsp[i],lsp,usp));

    if (data->trace > 0) {
      for (int i = 0; i < data->sp.size(); ++i)
	Rprintf("sp[%i]=%f\t",i,data->sp[i]);
    }

    BFGS2<Data> bfgs; 
    bfgs.reltol = data->reltol_search;
    bfgs.optimWithConstraint(pfminfn<Smooth,Stpm2>, pgrfn<Smooth,Stpm2>, data->init, ex, fminfn_constraint<Data>, false); // do not apply parscale b4/after

    NumericMatrix hessian0 = bfgs.calc_hessian(grfn<Data>, ex);

    double edf = trace(solve(as<mat>(bfgs.hessian),as<mat>(hessian0)));
    double negll = bfgs.calc_objective(fminfn<Data>, ex);
    double gcv =  negll + data->alpha*edf;
    double bic =  2.0*negll + data->alpha*edf*log(sum(data->event));
    data->init = bfgs.coef;

    // simple boundary constraints
    double constraint = 0.0;
    // double kappa = 1.0;
    for (int i=0; i < data->sp.size(); ++i) {
      if (logsp[i] < lsp)
	constraint +=  pow(logsp[i]-lsp,2.0);
      if (logsp[i] > usp)
	constraint +=  pow(logsp[i]-usp,2.0);
    }
    double objective =  data->criterion==1 ? 
      gcv + data->kappa/2.0*constraint : 
      bic + data->kappa/2.0*constraint;

    if (data->trace > 0) {
      Rprintf("edf=%f\tnegll=%f\tgcv=%f\tbic=%f\tobjective=%f\n",edf,negll,gcv,bic,objective);
    }
    return objective;
  }    


  template<class Smooth, class Stpm2>
  void pstpm2_step_multivariate_nlm(int n, double * logsp, double * f, void *ex) {
    *f = pstpm2_step_multivariate<Smooth,Stpm2>(n, logsp, ex);
  };


  template<class Smooth, class Stpm2>
  SEXP optim_pstpm2_first(SEXP args) { 

    typedef pstpm2<Smooth,Stpm2> Data;
    Data data(args);
    int n = data.init.size();

    for (int i=0; i<n; ++i) data.init[i] /= data.parscale[i];
    double opt_sp = exp(Brent_fmin(log(0.001),log(1000.0),&pstpm2_step_first<Smooth,Stpm2>,&data,1.0e-2));
    data.sp[0] = opt_sp;
    for (int i=0; i<n; ++i) data.init[i] *= data.parscale[i];

    BFGS2<Data> bfgs;
    bfgs.coef = data.init;
    bfgs.reltol = data.reltol;
    bfgs.optimWithConstraint(pfminfn<Smooth,Stpm2>, pgrfn<Smooth,Stpm2>, data.init, (void *) &data, fminfn_constraint<Data>);

    NumericMatrix hessian0 = bfgs.calc_hessian(grfn<Data>, (void *) &data);
    double edf = trace(solve(as<mat>(bfgs.hessian),as<mat>(hessian0)));

    return List::create(_("sp")=wrap(opt_sp),
			_("coef")=wrap(bfgs.coef),
			_("hessian")=wrap(bfgs.hessian),
          _("edf")=wrap(edf));

  }

  RcppExport SEXP optim_pstpm2LogH_first_ph(SEXP args) {
    return optim_pstpm2_first<SmoothLogH,stpm2>(args); }
  RcppExport SEXP optim_pstpm2Haz_first_ph(SEXP args) {
    return optim_pstpm2_first<SmoothHaz,stpm2>(args); }
  RcppExport SEXP optim_pstpm2LogH_first_po(SEXP args) {
    return optim_pstpm2_first<SmoothLogH,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2Haz_first_po(SEXP args) {
    return optim_pstpm2_first<SmoothHaz,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2LogH_first_probit(SEXP args) {
    return optim_pstpm2_first<SmoothLogH,stpm2Probit>(args); }
  RcppExport SEXP optim_pstpm2Haz_first_probit(SEXP args) {
    return optim_pstpm2_first<SmoothHaz,stpm2Probit>(args); }
  RcppExport SEXP optim_pstpm2LogH_first_ah(SEXP args) {
    return optim_pstpm2_first<SmoothLogH,stpm2AH>(args); }
  
  
  template<class Smooth, class Stpm2>
  SEXP optim_pstpm2_multivariate(SEXP args) {

    typedef pstpm2<Smooth,Stpm2> Data;
    Data data(args);
    int n = data.init.size();
    data.kappa = 10.0;

    NelderMead nm;
    nm.reltol = 1.0e-5;
    nm.maxit = 500;
    nm.hessianp = false;
    for (int i=0; i<n; ++i) data.init[i] /= data.parscale[i];

    NumericVector logsp(data.sp.size());
    for (int i=0; i < data.sp.size(); ++i)
      logsp[i] = log(data.sp[i]);

    bool satisfied;
    do {
      nm.optim(pstpm2_step_multivariate<Smooth,Stpm2>, logsp, (void *) &data);
      satisfied = true;
      for (int i=0; i < data.sp.size(); ++i)
	if (logsp[i]< -7.0 || logsp[i] > 7.0) satisfied = false;
      if (!satisfied) data.kappa *= 2.0;
    } while (!satisfied && data.kappa<1.0e5);

    for (int i=0; i<n; ++i) data.init[i] *= data.parscale[i];
    // for (int i=0; i < data.sp.size(); ++i)
    //   data.sp[i] = exp(logsp[i]);
    for (int i=0; i < nm.coef.size(); ++i)
      data.sp[i] = exp(nm.coef[i]);

    BFGS2<Data> bfgs;
    bfgs.coef = data.init;
    bfgs.reltol = data.reltol;
    data.kappa = 1.0;
    bfgs.optimWithConstraint(pfminfn<Smooth,Stpm2>, pgrfn<Smooth,Stpm2>, data.init, (void *) &data, fminfn_constraint<Data>);

    NumericMatrix hessian0 = bfgs.calc_hessian(grfn<Data>, (void *) &data);
    double edf = trace(solve(as<mat>(bfgs.hessian),as<mat>(hessian0)));

    return List::create(_("sp")=wrap(data.sp),
  			_("coef")=wrap(bfgs.coef),
  			_("hessian")=wrap(bfgs.hessian),
        _("edf")=wrap(edf));

  }

  RcppExport SEXP optim_pstpm2LogH_multivariate_ph(SEXP args) { 
    return optim_pstpm2_multivariate<SmoothLogH,stpm2>(args); }
  RcppExport SEXP optim_pstpm2Haz_multivariate_ph(SEXP args) {
    return optim_pstpm2_multivariate<SmoothHaz,stpm2>(args); }
  RcppExport SEXP optim_pstpm2LogH_multivariate_po(SEXP args) {
    return optim_pstpm2_multivariate<SmoothLogH,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2Haz_multivariate_po(SEXP args) {
    return optim_pstpm2_multivariate<SmoothHaz,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2LogH_multivariate_probit(SEXP args) {
    return optim_pstpm2_multivariate<SmoothLogH,stpm2Probit>(args); }
  RcppExport SEXP optim_pstpm2Haz_multivariate_probit(SEXP args) {
    return optim_pstpm2_multivariate<SmoothHaz,stpm2Probit>(args); }
  RcppExport SEXP optim_pstpm2LogH_multivariate_ah(SEXP args) { 
    return optim_pstpm2_multivariate<SmoothLogH,stpm2AH>(args); }

  template<class Smooth, class Stpm2>
  SEXP optim_pstpm2_multivariate_nlm(SEXP args) {

    typedef pstpm2<Smooth,Stpm2> Data;
    Data data(args);
    int n = data.init.size();
    data.kappa = 10.0;

    Nlm nlm;
    nlm.iagflg = 0;
    nlm.gradtl = 1.0e-4;
    nlm.steptl = 1.0e-4;
    nlm.msg = 9 + 4;
    nlm.hessianp = false;
    for (int i=0; i<n; ++i) data.init[i] /= data.parscale[i];

    NumericVector logsp(data.sp.size());
    for (int i=0; i < data.sp.size(); ++i)
      logsp[i] = log(data.sp[i]);

    bool satisfied;
    do {
      nlm.optim(pstpm2_step_multivariate_nlm<Smooth,Stpm2>, (fcn_p) 0, logsp, (void *) &data);
      satisfied = true;
      for (int i=0; i < data.sp.size(); ++i)
	if (logsp[i]< -7.0 || logsp[i] > 7.0) satisfied = false;
      if (!satisfied) data.kappa *= 2.0;
    } while (!satisfied && data.kappa<1.0e5);

    for (int i=0; i<n; ++i) data.init[i] *= data.parscale[i];
    for (int i=0; i < nlm.coef.size(); ++i)
      data.sp[i] = exp(nlm.coef[i]);

    BFGS2<Data> bfgs;
    bfgs.coef = data.init;
    bfgs.reltol = data.reltol;
    data.kappa = 1.0;
    bfgs.optimWithConstraint(pfminfn<Smooth,Stpm2>, pgrfn<Smooth,Stpm2>, data.init, (void *) &data, fminfn_constraint<Data>);

    NumericMatrix hessian0 = bfgs.calc_hessian(grfn<Data>, (void *) &data);
    double edf = trace(solve(as<mat>(bfgs.hessian),as<mat>(hessian0)));

    return List::create(_("sp")=wrap(data.sp),
  			_("coef")=wrap(bfgs.coef),
  			_("hessian")=wrap(bfgs.hessian),
			_("edf")=wrap(edf));
  }

  RcppExport SEXP optim_pstpm2LogH_multivariate_ph_nlm(SEXP args) { 
    return optim_pstpm2_multivariate_nlm<SmoothLogH,stpm2>(args); }
  RcppExport SEXP optim_pstpm2Haz_multivariate_ph_nlm(SEXP args) {
    return optim_pstpm2_multivariate_nlm<SmoothHaz,stpm2>(args); }
  RcppExport SEXP optim_pstpm2LogH_multivariate_po_nlm(SEXP args) {
    return optim_pstpm2_multivariate_nlm<SmoothLogH,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2Haz_multivariate_po_nlm(SEXP args) {
    return optim_pstpm2_multivariate_nlm<SmoothHaz,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2LogH_multivariate_probit_nlm(SEXP args) {
    return optim_pstpm2_multivariate_nlm<SmoothLogH,stpm2Probit>(args); }
  RcppExport SEXP optim_pstpm2Haz_multivariate_probit_nlm(SEXP args) {
    return optim_pstpm2_multivariate_nlm<SmoothHaz,stpm2Probit>(args); }


  template<class Smooth, class Stpm2>
  SEXP optim_pstpm2_fixedsp(SEXP args) { 

    typedef pstpm2<Smooth,Stpm2> Data;
    Data data(args);
    BFGS2<Data> bfgs;
    bfgs.coef = data.init;
    bfgs.reltol = data.reltol;

    bfgs.optimWithConstraint(pfminfn<Smooth,Stpm2>, pgrfn<Smooth,Stpm2>, data.init, (void *) &data, fminfn_constraint<Data>);

    NumericMatrix hessian0 = bfgs.calc_hessian(grfn<Data>, (void *) &data);
    double edf = trace(solve(as<mat>(bfgs.hessian),as<mat>(hessian0)));

    return List::create(_("sp")=wrap(data.sp),
			_("coef")=wrap(bfgs.coef),
			_("hessian")=wrap(bfgs.hessian),
			_("edf")=wrap(edf));
  }

  RcppExport SEXP optim_pstpm2LogH_fixedsp_ph(SEXP args) {
    return optim_pstpm2_fixedsp<SmoothLogH,stpm2>(args); }
  RcppExport SEXP optim_pstpm2Haz_fixedsp_ph(SEXP args) {
    return optim_pstpm2_fixedsp<SmoothHaz,stpm2>(args); }
  RcppExport SEXP optim_pstpm2LogH_fixedsp_po(SEXP args) {
    return optim_pstpm2_fixedsp<SmoothLogH,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2Haz_fixedsp_po(SEXP args) {
    return optim_pstpm2_fixedsp<SmoothHaz,stpm2PO>(args); }
  RcppExport SEXP optim_pstpm2LogH_fixedsp_probit(SEXP args) { 
    return optim_pstpm2_fixedsp<SmoothLogH,stpm2Probit>(args); }
  RcppExport SEXP optim_pstpm2Haz_fixedsp_probit(SEXP args) {
    return optim_pstpm2_fixedsp<SmoothHaz,stpm2Probit>(args); }
  RcppExport SEXP optim_pstpm2LogH_fixedsp_ah(SEXP args) {
    return optim_pstpm2_fixedsp<SmoothLogH,stpm2AH>(args); }

  template<class Smooth, class Stpm2>
  SEXP test_pstpm2(SEXP args) { 

    pstpm2<Smooth,Stpm2> data(args);
    int n = data.init.size();

    NumericVector init = ew_div(data.init,data.parscale);

    double pfmin = pfminfn<Smooth,Stpm2>(n,&init[0],(void *) &data);
    NumericVector gr(n);
    pgrfn<Smooth,Stpm2>(n,&init[0], &gr[0], (void *) &data);

    init = ew_mult(init,data.parscale);

    return List::create(_("sp")=wrap(data.sp),
			_("pfmin")=wrap(pfmin),
			_("gr")=wrap(gr)
			);

  }

  RcppExport SEXP test_pstpm2LogH(SEXP args) {
    return test_pstpm2<SmoothLogH,stpm2>(args);
  }

  RcppExport SEXP test_pstpm2Haz(SEXP args) {
    return test_pstpm2<SmoothHaz,stpm2>(args);
  }

  // R CMD INSTALL ~/src/R/microsimulation
  // R -q -e "require(microsimulation); .Call('test_nmmin',1:2,PACKAGE='microsimulation')"

  // .Call("optim_stpm2",init,X,XD,rep(bhazard,nrow(X)),wt,ifelse(event,1,0),package="rstpm2")

} // anonymous namespace