text:
// Author: Daisuke Kanaizumi
// Affiliation: Department of Applied Mathematics, Waseda University

// verification program for q-Bessel functions
// April 12th, 2018

#ifndef QBESSEL_HPP
#define QBESSEL_HPP

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/constants.hpp>
#include <kv/complex.hpp>
#include <kv/convert.hpp>
#include <kv/defint.hpp>
#include <kv/Heine.hpp>
#include <kv/Pochhammer.hpp>
#include <kv/QHypergeometric.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <limits>
namespace ub = boost::numeric::ublas;
namespace kv{
template <class T> interval<T>Jackson1(const interval<T>& z,const interval<T>& nu,const interval<T>& q){
// verification program for Jackson`s 1st q-Bessel function
interval<T>res,a,b,c,series;
int j,K;
T rad;
if(abs(z)>=2){
throw std::domain_error("Jackson`s 1st q-Bessel function is not defined");
}
if(abs(q)>=1){
throw std::domain_error("absolute value of q must be under 1");
}
 if(z<0){
   throw std::domain_error("z must be positive");
 }
 if (pow(q,nu)<1){
j=1;
K=500;
a=1.;
b=1.;
c=1.;
// Leibniz criterion
 for (int k=1; k<=K; k++){
    j = -1*j;
    b=b*(1-pow(q,k));
    c=c*(1-pow(q,k+nu));
    a = a+j*pow(z/2,2*k)/(b*c);
  }
b=b*(1-pow(q,K+1));
c=c*(1-pow(q,K+nu+1));

rad=abs(pow(z/2,2*K+2)/(b*c)).upper();
series=a+rad*interval<T>(-1.,1.);
res=pow(z/2,nu)*series/Karpelevich(interval<T>(pow(q,nu+1)),interval<T>(q));
}
 else{
res=pow(z/2,nu)*infinite_qPochhammer(interval<T>(pow(q,nu+1)),interval<T>(q))*
  Heine(interval<T>(0),interval<T>(0),interval<T>(pow(q,nu+1)),interval<T>(q),interval<T>(-z*z/4))/Euler(interval<T>(q));
}

 return res;
}
template <class T> complex<interval<T> >Jackson1(const complex<interval<T> >& z,const complex<interval<T> >& nu,const interval<T>& q){
  complex<interval<T> >res,i;

if(abs(q)>=1){
throw std::domain_error("absolute value of q must be under 1");
}
 if(abs(z).upper()>=2){
throw std::domain_error("Jackson`s 1st q-Bessel function is not defined");
}

res=pow(z/2,nu)*infinite_qPochhammer(complex<interval<T> >(pow(q,nu+1)),interval<T>(q))*
  Heine(complex<interval<T> >(0),complex<interval<T> >(0),complex<interval<T> >(pow(q,nu+1)),interval<T>(q),complex<interval<T> >(-z*z/4))/Euler(interval<T>(q));

 return res;
}

template <class T> complex<interval<T> >modified_qBesselI1(const complex<interval<T> >& z,const complex<interval<T> >& nu,const interval<T>& q){
  // verification program for 1st modified q-Bessel I1
  complex<interval<T> >res,i;
  if(abs(q)>=1){
    throw std::domain_error("absolute value of q must be under 1");
  }
  if(abs(z).upper()>=2){
    throw std::domain_error("1st modified q-Bessel function is not defined");
  }
  interval<T>pi;
  i=complex<interval<T> >::i();
  pi=constants<interval<T> >::pi();
  res=exp(-i*nu*pi/2)*Jackson1(complex<interval<T> >(i*z),complex<interval<T> >(nu),interval<T> (q));
  return res;
}

  template <class T> complex<interval<T> >modified_qBesselK1(const complex<interval<T> >& z,const complex<interval<T> >& nu,const interval<T>& q){
    // verification program for 1st modified q-Bessel K1, nu should not be an integer
    complex<interval<T> >res;
    if(abs(q)>=1){
      throw std::domain_error("absolute value of q must be under 1");
    }
    if(abs(z).upper()>=2){
      throw std::domain_error("1st modified q-Bessel function is not defined");
    }
    interval<T>pi;
    pi=constants<interval<T> >::pi();
    res=pi*0.5/sin(pi*nu)*
      (modified_qBesselI1(complex<interval<T> > (z), complex<interval<T> > (-nu), interval<T> (q))-modified_qBesselI1(complex<interval<T> > (z), complex<interval<T> > (nu), interval<T> (q)));
    return res;
  }
  template <class T> interval<T>Jackson2(const interval<T>& z,const interval<T>& nu,const interval<T>& q){
    // verification program for Jackson`s 2nd q-Bessel function
  interval<T>res,a,b,c,series,pq;
int j,K;
T rad;
if(abs(q)>=1){
throw std::domain_error("absolute value of q must be under 1");
}

// if (pow(q,nu)<1 && z*z<4*q){
/*j=1;
K=500;
a=1.;
b=1.;
c=1.;
// Leibniz criterion
 for (int k=1; k<=K; k++){
    j = -1*j;
    b=b*(1-pow(q,k));
    c=c*(1-pow(q,k+nu));
    a = a+j*pow(q,k*(k-1))*pow(pow(q,nu+1)*z*z/4,k)/(b*c);
  }
  b=b*(1-pow(q,K+1));
  c=c*(1-pow(q,K+nu+1));

 rad=abs(pow(q,K*(K+1))*pow(pow(q,nu+1)*z*z/4,K+1)/(b*c)).upper();
 series=a+rad*interval<T>(-1.,1.);
 res=pow(z/2,nu)*series/Karpelevich(interval<T>(pow(q,nu+1)),interval<T>(q));
*/ //}
// else{
   /* implementation by original definition*/
// res=pow(z*0.5,nu)*infinite_qPochhammer(interval<T>(pow(q,nu+1)),interval<T>(q))*
//  _0phi_1(interval<T>(pow(q,nu+1)),interval<T>(q),interval<T>(-z*z*pow(q,nu+1)/4))/Euler(interval<T>(q));

 // alternative implementaion
 // reference
 // H. T Koelink, Hansen-Lommel orthogonality Relations for Jackson`s q-Bessel functions, formula 3.2
 // Journal of Mathematical Analysis and Applications 175, 425-437 (1993)
 pq=pow(q,nu+1);
 res=pow(z/2,nu)*_1phi_1(interval<T> (-z*z/4),interval<T> (0),interval<T>(q), interval<T> (pq))/Euler(interval<T>(q));
 // }

  return res;
}
template <class T> complex<interval<T> >Jackson2(const complex<interval<T> >& z,const complex<interval<T> >& nu,const interval<T>& q){
  complex<interval<T> >res,pq,i;
  i=complex<interval<T> >::i();
  if(abs(q)>=1){
    throw std::domain_error("absolute value of q must be under 1");
  }
  pq=pow(q,nu+1);
  /* implementation by original definition*/

  //  res=pow(z/2,nu)*infinite_qPochhammer(complex<interval<T> >(pq),interval<T>(q))*
  //  _0phi_1(complex<interval<T> >(pq),interval<T>(q),complex<interval<T> >(-z*z*pq/4))/Euler(interval<T>(q));

  // alternative implementaion
  // reference
  // H. T Koelink, Hansen-Lommel Orthogonality Relations for Jackson`s q-Bessel functions, formula 3.2
  // Journal of Mathematical Analysis and Applications 175, 425-437 (1993)
   try{
  res=pow(z/2,nu)*_1phi_1(complex<interval<T> >(-z*z/4),complex<interval<T> >(0),interval<T>(q), complex<interval<T> >(pq))/Euler(interval<T>(q));
  }
  // alternative implementaion
  // reference
  // Y. Chen , M. E. Ismail, K. A. Muttalib. Asymptotics of basic Bessel functions and q-Laguerre polynomials, Lemma 2
  // Journal of Computational and Applied Mathematics, 54(3), 263-272 (1994).
  catch(std::domain_error){
    ub::vector< complex<interval<T> > > a(3);
    ub::vector< complex<interval<T> > > b(2);
    ub::vector< complex<interval<T> > > c(2);
    a(0)=pow(q,(nu+0.5)*0.5);
    a(1)=-pow(q,(nu+0.5)*0.5);
    a(2)=0.;
    b(0)=-sqrt(q);
    b(1)=i*a(0)*z*0.5;
    c(0)=-sqrt(q);
    c(1)=-b(1);

    res=pow(z*0.5,nu)*infinite_qPochhammer(interval<T>(sqrt(q)),interval<T>(q))/Euler(interval<T>(q))*0.5
      *(infinite_qPochhammer(complex<interval<T> >(b(1)),interval<T>(sqrt(q)))
	*QHypergeom(ub::vector<complex<interval<T> > >(a),ub::vector<complex<interval<T> > >(b),interval<T>(sqrt(q)),complex<interval<T> >(sqrt(q),0))+
	infinite_qPochhammer(complex<interval<T> >(c(1)),interval<T>(sqrt(q)))
	*QHypergeom(ub::vector<complex<interval<T> > >(a),ub::vector<complex<interval<T> > >(c),interval<T>(sqrt(q)),complex<interval<T> >(sqrt(q),0)));
    	}

    if((abs(res)).upper()==std::numeric_limits<T>::infinity()){
    ub::vector< complex<interval<T> > > a(3);
    ub::vector< complex<interval<T> > > b(2);
    ub::vector< complex<interval<T> > > c(2);
    a(0)=pow(q,(nu+0.5)*0.5);
    a(1)=-pow(q,(nu+0.5)*0.5);
    a(2)=0.;
    b(0)=-sqrt(q);
    b(1)=i*a(0)*z*0.5;
    c(0)=-sqrt(q);
    c(1)=-b(1);

    res=pow(z*0.5,nu)*infinite_qPochhammer(interval<T>(sqrt(q)),interval<T>(q))/Euler(interval<T>(q))*0.5
      *(infinite_qPochhammer(complex<interval<T> >(b(1)),interval<T>(sqrt(q)))
	*QHypergeom(ub::vector<complex<interval<T> > >(a),ub::vector<complex<interval<T> > >(b),interval<T>(sqrt(q)),complex<interval<T> >(sqrt(q),0))+
	infinite_qPochhammer(complex<interval<T> >(c(1)),interval<T>(sqrt(q)))
	*QHypergeom(ub::vector<complex<interval<T> > >(a),ub::vector<complex<interval<T> > >(c),interval<T>(sqrt(q)),complex<interval<T> >(sqrt(q),0)));
     }
    //std::cout<<qPochhammer(complex<interval<T> >(pow(q,(nu+0.5)*0.5)),interval<T>(q),int(1000))
    //*qPochhammer(complex<interval<T> >(-pow(q,(nu+0.5)*0.5)),interval<T>(q),int(1000))<<std::endl;

    return res;
}
  template <class TT> struct qBesselintegral_nu_int_real {
    TT x, q;
    int nu,n; // Setting parameters
    qBesselintegral_nu_int_real(TT x, TT q, int nu,int n) : x(x),q(q),nu(nu),n(n) {}

    template <class T> T operator() (const T& t) {
      complex<T>  pro;
      T proreal;
      complex<T> i;
      pro=1.;
      i=complex<T>::i();
      for(int k=0;k<=nu-1;k++){
	pro=pro*(1-pow(T(q),k)*exp(2*i*t))*(1-pow(T(q),k)*exp(-2*i*t));
      }
      for(int j=0;j<=n-1;j++){
	pro=pro*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(i*t)/2)*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(-i*t)/2);
      }

      proreal=pro.real();
      return proreal;
    }
  };
  template <class TT> struct qBesselintegral_nu_int_imag {
    TT x, q;
    int nu,n; // Setting parameters
    qBesselintegral_nu_int_imag(TT x, TT q, int nu,int n) : x(x),q(q),nu(nu),n(n) {}

    template <class T> T operator() (const T& t) {
      complex<T>  pro;
      T proimag;
      complex<T> i;
      pro=1.;
      i=complex<T>::i();
      for(int k=0;k<=nu-1;k++){
	pro=pro*(1-pow(T(q),k)*exp(2*i*t))*(1-pow(T(q),k)*exp(-2*i*t));
      }
      for(int j=0;j<=n-1;j++){
	pro=pro*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(i*t)/2)*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(-i*t)/2);
      }

      proimag=pro.imag();
      return proimag;
    }
  };
  template <class T> complex<interval<T> >Jackson2_integral(const interval<T>& z,const int & nu,const interval<T>& q){
      // verification program for Jackson`s 2nd q-Bessel function
      // Integral representation is used

      // References
      // Rahman(1987), An Integral Representation and Some Transformation Properties of q-Bessel Functions, Journal of Mathematical Analysis and Applications 125
      // Zhang(2008), Plancherel-Rotach asymptotics for certain basic hypergeometric series, Advances in Mathematics 217, Lemma 1.1
    int n;
    n=100;
    if(nu<=0){
      throw std::domain_error("nu must be positive");
    }
    if(abs(z)*pow(q,nu/2+n+0.5)/2/(1-q)>=0.5){
      n=n+10;
    }
    complex<interval<T> >res,integral;
    interval<T>num,realint,imagint,pi;
    T numrad;
    numrad=(abs(z)*pow(q,nu/2+n+0.5)*2/(1-q)).upper();
    num=pow((1+numrad*interval<T>(-1.,1.)),2);
    pi=constants<interval<T> >::pi();
    realint=defint(qBesselintegral_nu_int_real<interval<T> >(z,q,nu,n),interval<T>(0.),interval<T>(pi),10,10);
    imagint=defint(qBesselintegral_nu_int_imag<interval<T> >(z,q,nu,n),interval<T>(0.),interval<T>(pi),10,10);
    integral=complex<interval<T> >(realint,imagint);

    res=num*integral*pow(z/2,nu)*infinite_qPochhammer(interval<T>(pow(q,2*nu)),interval<T>(q))/infinite_qPochhammer(interval<T>(pow(q,nu)),interval<T>(q))/(2*pi);

    return res;
  }
    template <class TT> struct qBesselintegral_nu_double_real {
      TT x,q,nu;
      int n; // Setting parameters
      qBesselintegral_nu_double_real(TT x, TT q, TT nu,int n) : x(x),q(q),nu(nu),n(n) {}

      template <class T> T operator() (const T& t) {
      complex<T>  pro;
      T proreal;
      complex<T> i;
      pro=1.;
      i=complex<T>::i();
      for(int k=0;k<=n-1;k++){
	pro=pro*(1-pow(T(q),k)*exp(2*i*t))*(1-pow(T(q),k)*exp(-2*i*t));
      }
      for(int j=0;j<=n-1;j++){
	pro=pro*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(i*t)/2)*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(-i*t)/2);
      }
      for(int l=0;l<=n-1;l++){
	pro=pro/(1-pow(T(q),l+nu)*exp(2*i*t))/(1-pow(T(q),l+nu)*exp(-2*i*t));
      }
      proreal=pro.real();
      return proreal;
      }
  };
  template <class TT> struct qBesselintegral_nu_double_imag {
    TT x,q,nu;
    int n; // Setting parameters
    qBesselintegral_nu_double_imag(TT x, TT q, TT nu,int n) : x(x),q(q),nu(nu),n(n) {}

    template <class T> T operator() (const T& t) {
      complex<T>  pro;
      T proimag;
      complex<T> i;
      pro=1.;
      i=complex<T>::i();
      for(int k=0;k<=n-1;k++){
	pro=pro*(1-pow(T(q),k)*exp(2*i*t))*(1-pow(T(q),k)*exp(-2*i*t));
      }
      for(int j=0;j<=n-1;j++){
	pro=pro*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(i*t)/2)*(1+i*T(x)*pow(T(q),T(nu)/2+0.5+j)*exp(-i*t)/2);
      }
      for(int l=0;l<=n-1;l++){
	pro=pro/(1-pow(T(q),l+nu)*exp(2*i*t))/(1-pow(T(q),l+nu)*exp(-2*i*t));
      }

      proimag=pro.imag();
      return proimag;
    }
  };
  template <class T> complex<interval<T> >Jackson2_integral(const interval<T>& z,const interval<T> & nu,const interval<T>& q){
      // verification program for Jackson`s 2nd q-Bessel function
      // Integral representation is used

      // References
      // Rahman(1987), An Integral Representation and Some Transformation Properties of q-Bessel Functions, Journal of Mathematical Analysis and Applications 125
      // Zhang(2008), Plancherel-Rotach asymptotics for certain basic hypergeometric series, Advances in Mathematics 217, Lemma 1.1
    int n;
    n=100;
    if(nu<=0){
      throw std::domain_error("nu must be positive");
    }
    if(abs(z)*pow(q,nu/2+n+0.5)/2/(1-q)>=0.5){
      n=n+10;
    }
    if(pow(q,nu+n)/2/(1-q)>=0.5){
      n=n+10;
    }

    complex<interval<T> >res,integral;
    interval<T>num1,num2,denom,realint,imagint,pi;
    T numrad1,numrad2,denomrad;
    numrad1=(abs(z)*pow(q,nu/2+n+0.5)*2/(1-q)).upper();
    num1=pow((1+numrad1*interval<T>(-1.,1.)),2);
    numrad2=(pow(q,n)*2/(1-q)).upper();
    num2=pow((1+numrad2*interval<T>(-1.,1.)),2);
    denomrad=(pow(q,nu+n)*2/(1-q)).upper();
    denom=pow((1+denomrad*interval<T>(-1.,1.)),2);
    pi=constants<interval<T> >::pi();
    realint=defint(qBesselintegral_nu_double_real<interval<T> >(z,q,nu,n),interval<T>(0.),interval<T>(pi),10,10);
    imagint=defint(qBesselintegral_nu_double_imag<interval<T> >(z,q,nu,n),interval<T>(0.),interval<T>(pi),10,10);
    integral=complex<interval<T> >(realint,imagint);

    res=num1*num2*denom*integral*pow(z/2,nu)*infinite_qPochhammer(interval<T>(pow(q,2*nu)),interval<T>(q))/infinite_qPochhammer(interval<T>(pow(q,nu)),interval<T>(q))/(2*pi);

    return res;
  }
    template <class TT> struct qBesselintegral_z_complex_real {
      TT q,nu;
      complex<TT> x;
      int n; // Setting parameters
      qBesselintegral_z_complex_real(complex<TT> x, TT q, TT nu,int n) : x(x),q(q),nu(nu),n(n) {}

      template <class T> T operator() (const T& t) {
      complex<T>  pro;
      T proreal;
      complex<T> i;
      pro=1.;
      i=complex<T>::i();
      for(int k=0;k<=n-1;k++){
	pro=pro*(1-pow(T(q),k)*exp(2*i*t))*(1-pow(T(q),k)*exp(-2*i*t));
      }
      for(int j=0;j<=n-1;j++){
	pro=pro*(1+i*complex<T>(x)*pow(T(q),T(nu)/2+0.5+j)*exp(i*t)/2)*(1+i*complex<T>(x)*pow(T(q),T(nu)/2+0.5+j)*exp(-i*t)/2);
      }
      for(int l=0;l<=n-1;l++){
	pro=pro/(1-pow(T(q),l+T(nu))*exp(2*i*t))/(1-pow(T(q),l+T(nu))*exp(-2*i*t));
      }
      proreal=pro.real();
      return proreal;
      }
  };
  template <class TT> struct qBesselintegral_z_complex_imag {
    TT q,nu;
    complex<TT> x;
    int n; // Setting parameters
    qBesselintegral_z_complex_imag(complex<TT> x, TT q, TT nu,int n) : x(x),q(q),nu(nu),n(n) {}

    template <class T> T operator() (const T& t) {
      complex<T>  pro;
      T proimag;
      complex<T> i;
      pro=1.;
      i=complex<T>::i();
      for(int k=0;k<=n-1;k++){
	pro=pro*(1-pow(T(q),k)*exp(2*i*t))*(1-pow(T(q),k)*exp(-2*i*t));
      }
      for(int j=0;j<=n-1;j++){
	pro=pro*(1+i*complex<T>(x)*pow(T(q),T(nu)/2+0.5+j)*exp(i*t)/2)*(1+i*complex<T>(x)*pow(T(q),T(nu)/2+0.5+j)*exp(-i*t)/2);
      }
      for(int l=0;l<=n-1;l++){
	pro=pro/(1-pow(T(q),l+T(nu))*exp(2*i*t))/(1-pow(T(q),l+T(nu))*exp(-2*i*t));
      }

      proimag=pro.imag();
      return proimag;
    }
  };
  template <class T> complex<interval<T> >Jackson2_integral(const complex<interval<T> >& z,const interval<T> & nu,const interval<T>& q){
      // verification program for Jackson`s 2nd q-Bessel function
      // Integral representation is used

      // References
      // Rahman(1987), An Integral Representation and Some Transformation Properties of q-Bessel Functions, Journal of Mathematical Analysis and Applications 125
      // Zhang(2008), Plancherel-Rotach asymptotics for certain basic hypergeometric series, Advances in Mathematics 217, Lemma 1.1
    int n;
    n=100;
    if(nu<=0){
      throw std::domain_error("nu must be positive");
    }
    if(abs(z*pow(q,nu/2+n+0.5))/2/(1-q)>=0.5){
      n=n+10;
    }
    if(abs(pow(q,nu+n))/2/(1-q)>=0.5){
      n=n+10;
    }

    complex<interval<T> >res,integral,num1,denom;

    interval<T> num2;
    interval<T> pi,realint,imagint;
    T numrad1,numrad2,denomrad;
    numrad1=(abs(z*pow(q,nu/2+n+0.5))*2/(1-q)).upper();
    num1=pow(complex_nbd(complex<interval<T> >(1,0),numrad1),2);
    numrad2=(pow(q,n)*2/(1-q)).upper();
    num2=pow((1+numrad2*interval<T>(-1.,1.)),2);
    denomrad=(abs(pow(q,nu+n))*2/(1-q)).upper();
    denom=pow(complex_nbd(complex<interval<T> >(1,0),denomrad),2);
    pi=constants<interval<T> >::pi();
    realint=defint(qBesselintegral_z_complex_real<interval<T> >(z,q,nu,n),interval<T>(0.),interval<T>(pi),10,10);
    imagint=defint(qBesselintegral_z_complex_imag<interval<T> >(z,q,nu,n),interval<T>(0.),interval<T>(pi),10,10);
    integral=complex<interval<T> >(realint,imagint);

    res=num1*num2*denom*integral*pow(z/2,nu)*infinite_qPochhammer(complex<interval<T> >(pow(q,2*nu)),interval<T>(q))/infinite_qPochhammer(complex<interval<T> >(pow(q,nu)),interval<T>(q))/(2*pi);

    return res;
  }
  template <class T> interval<T> J2ratio(const interval<T> & z,const interval<T> & nu,const interval<T>& q){
    interval<T>res;
    res=Jackson2(interval<T>(z),interval<T>(nu),interval<T>(q))/Jackson2(interval<T>(z),interval<T>(nu-1),interval<T>(q));
    return res;
  }
  template <class T> complex<interval<T> >modified_qBesselI2(const complex<interval<T> >& z,const complex<interval<T> >& nu,const interval<T>& q){
    //verification program for 2nd modified q-Bessel function I2
    complex<interval<T> >res,i;
    if(abs(q)>=1){
      throw std::domain_error("absolute value of q must be under 1");
    }
    /*implementation by original definition
      interval<T>pi;
      i=complex<interval<T> >::i();
      pi=constants<interval<T> >::pi();
      res=exp(-i*nu*pi/2)*Jackson2(complex<interval<T> >(i*z),complex<interval<T> >(nu),interval<T> (q));
    */
    // alternative implementation
    // reference
    // Ismail, M. E., & Zhang, R. (2015). $ q $-Bessel Functions and Rogers-Ramanujan Type Identities.
    // arXiv preprint arXiv:1508.06861.
    res=pow(z/2,nu)/Euler(interval<T>(q))
      *_1phi_1(complex<interval<T> >(z*z/4),complex<interval<T> >(0.),interval<T>(q),complex<interval<T> >(pow(q,nu+1)));
    return res;
  }
 template <class T> interval<T> modified_qBesselI2(const interval<T> & z,const interval<T> & nu,const interval<T>& q){
    //verification program for 2nd modified q-Bessel function I2
    interval<T> res;
    if(abs(q)>=1){
      throw std::domain_error("absolute value of q must be under 1");
    }
    /*implementation by original definition
      interval<T>pi;
      i=complex<interval<T> >::i();
      pi=constants<interval<T> >::pi();
      res=exp(-i*nu*pi/2)*Jackson2(complex<interval<T> >(i*z),complex<interval<T> >(nu),interval<T> (q));
    */
    // alternative implementation
    // reference
    // Ismail, M. E., & Zhang, R. (2015). $ q $-Bessel Functions and Rogers-Ramanujan Type Identities.
    // arXiv preprint arXiv:1508.06861.
    res=pow(z/2,nu)/Euler(interval<T>(q))
      *_1phi_1(interval<T> (z*z/4),interval<T> (0.),interval<T>(q),interval<T> (pow(q,nu+1)));
    return res;
  }
  /*template <class T> interval<T> modified_qBesselI2_ae(const interval<T> & z,const interval<T> & nu,const interval<T>& q){
    //verification program for 2nd modified q-Bessel function I2, z>0
    // reference
    // Ismail, M. E., & Zhang, R. (2015). $ q $-Bessel Functions and Rogers-Ramanujan Type Identities.
    // arXiv preprint arXiv:1508.06861.
    interval<T> res;
    if(abs(q)>=1){
      throw std::domain_error("absolute value of q must be under 1");
    }
    res=pow(z/2.,nu)*infinite_qPochhammer(interval<T>(sqrt(q)),interval<T>(q))*0.5/Euler(interval<T>(q))
      *(infinite_qPochhammer(interval<T>(z*0.5*pow(q,(nu+0.5)*0.5)),interval<T>(sqrt(q)))+infinite_qPochhammer(interval<T>(-z*0.5*pow(q,(nu+0.5)*0.5)),interval<T>(sqrt(q))));
    return res;
    }*/
  template <class T> interval<T> Hahn_Exton(const interval<T> & z,const interval<T> & nu,const interval<T>& q){
    // verification program for Hahn-Exton q-Bessel function
    interval<T> res,pq;
    if(abs(q)>=1){
      throw std::domain_error("absolute value of q must be under 1");
    }
    pq=pow(q,nu+1);
    /* implementation by original definition*/
    /*    res=pow(z,nu)*infinite_qPochhammer(interval<T> (pq),interval<T>(q))*
	  _1phi_1(interval<T> (0),interval<T> (pq),interval<T>(q),interval<T> (z*z*q))/Euler(interval<T>(q));*/
       // alternative implementation
       // A. B. Olde Daalhuis, Asymptotic Expansions for q-Gamma, q-Exponential and q-Bessel Functions, formula 4.6
       // Journal of Mathematical Analysis and Applications 186, 896-913 (1994)

    	 res=pow(z,nu)*infinite_qPochhammer(interval<T> (z*z*q),interval<T>(q))*_1phi_1(interval<T> (0),interval<T> (z*z*q),interval<T>(q),interval<T> (pq))/Euler(interval<T> (q));
       return res;
  }
  template <class T> complex<interval<T> >Hahn_Exton(const complex<interval<T> >& z,const complex<interval<T> >& nu,const interval<T>& q){
    // verification program for Hahn-Exton q-Bessel function
    complex<interval<T> >res,pq;
    if(abs(q)>=1){
      throw std::domain_error("absolute value of q must be under 1");
    }
    pq=pow(q,nu+1);
    /* implementation by original definition*/
    /*  res=pow(z,nu)*infinite_qPochhammer(complex<interval<T> >(pq),interval<T>(q))*
       _1phi_1(complex<interval<T> >(0),complex<interval<T> >(pq),interval<T>(q),complex<interval<T> >(z*z*q))/Euler(interval<T>(q));
    */
    // alternative implementation
    // A. B. Olde Daalhuis, Asymptotic Expansions for q-Gamma, q-Exponential and q-Bessel Functions, formula 4.6
    // Journal of Mathematical Analysis and Applications 186, 896-913 (1994)
    res=pow(z,nu)*infinite_qPochhammer(complex<interval<T> >(z*z*q),interval<T>(q))*_1phi_1(complex<interval<T> >(0),complex<interval<T> >(z*z*q),interval<T>(q), complex<interval<T> >(pq))/Euler(interval<T> (q));
     return res;
  }
   template <class T> interval<T> HEratio(const interval<T> & z,const interval<T> & nu,const interval<T>& q){
    interval<T>res;
    res=Hahn_Exton(interval<T>(z),interval<T>(nu),interval<T>(q))/Hahn_Exton(interval<T>(z),interval<T>(nu-1),interval<T>(q));
    return res;
  }

  template <class T> complex<interval<T> >little(const complex<interval<T> >& z,const complex<interval<T> >& nu,const interval<T>& q){
    // verification program for little q-Bessel function
    // reference
    // Koornwinder and Swarttouw, On q-analogues of the Fourier and Hankel transforms, 1992
    // Bouzeffour, New Addition Formula for the Little q-Bessel Functions, arXiv, 2013
    complex<interval<T> >res,pq;
    if(abs(q)>=1){
      throw std::domain_error("absolute value of q must be under 1");
    }
    pq=pow(q,nu+1);
    res=pow(z,nu)*infinite_qPochhammer(complex<interval<T> >(pq),interval<T>(q))*
      _1phi_1(complex<interval<T> >(0),complex<interval<T> >(pq),interval<T>(q),complex<interval<T> >(z))/Euler(interval<T>(q));
    return res;
  }

}

#endif

max_stars_repo_name: Daisuke-Kanaizumi/q-special-functions

max_issues_repo_name: Daisuke-Kanaizumi/q-special-functions

max_forks_repo_name: Daisuke-Kanaizumi/q-special-functions


text:
// Copyright (c) 2015-2021 Daniel Cooke
// Use of this source code is governed by the MIT license that can be found in the LICENSE file.

#ifndef bam_realigner_hpp
#define bam_realigner_hpp

#include <vector>
#include <cstddef>

#include <boost/optional.hpp>

#include "basics/aligned_read.hpp"
#include "core/types/haplotype.hpp"
#include "core/types/genotype.hpp"
#include "core/models/haplotype_likelihood_model.hpp"
#include "containers/mappable_flat_set.hpp"
#include "io/reference/reference_genome.hpp"
#include "io/read/read_reader.hpp"
#include "io/read/read_writer.hpp"
#include "io/variant/vcf_reader.hpp"
#include "utils/memory_footprint.hpp"
#include "utils/thread_pool.hpp"

namespace octopus {

class BAMRealigner
{
public:
    using ReadReader = io::ReadReader;
    using ReadWriter = io::ReadWriter;
    using SampleName = ReadReader::SampleName;
    using SampleList = std::vector<SampleName>;

    struct Config
    {
        HaplotypeLikelihoodModel alignment_model = {};
        bool primary_only = true;
        bool copy_hom_ref_reads = false;
        bool simplify_cigars = false;
        ReadLinkageConfig read_linkage = ReadLinkageConfig {};
        MemoryFootprint max_buffer = *parse_footprint("50M");
        boost::optional<unsigned> max_threads = 1;
    };

    struct Report
    {
        std::size_t n_reads_assigned;
        std::size_t n_reads_unassigned;
    };

    BAMRealigner() = default;
    BAMRealigner(Config config);

    BAMRealigner(const BAMRealigner&)            = delete;
    BAMRealigner& operator=(const BAMRealigner&) = delete;
    BAMRealigner(BAMRealigner&&)                 = delete;
    BAMRealigner& operator=(BAMRealigner&&)      = delete;

    ~BAMRealigner() = default;

    Report realign(ReadReader& src, VcfReader& variants, ReadWriter& dst,
                   const ReferenceGenome& reference, SampleList samples) const;
    Report realign(ReadReader& src, VcfReader& variants, ReadWriter& dst,
                   const ReferenceGenome& reference) const;

private:
    using VcfIterator = VcfReader::RecordIterator;
    using CallBlock   = std::vector<VcfRecord>;
    struct Batch
    {
        MappableFlatSet<Genotype<Haplotype>> genotypes;
        std::vector<AlignedRead> reads;
    };
    using BatchList = std::vector<Batch>;
    using BatchListRegionPair = std::pair<BatchList, boost::optional<GenomicRegion>>;

    Config config_;
    mutable ThreadPool workers_;

    CallBlock read_next_block(VcfIterator& first, const VcfIterator& last, const SampleList& samples) const;
    BatchListRegionPair read_next_batch(VcfIterator& first, const VcfIterator& last, ReadReader& src,
                                        const ReferenceGenome& reference, const SampleList& samples,
                                        const boost::optional<GenomicRegion>& prev_batch_region) const;
    void merge(BatchList& src, BatchList& dst) const;
};

BAMRealigner::Report
realign(io::ReadReader::Path src, VcfReader::Path variants, io::ReadWriter::Path dst,
        const ReferenceGenome& reference);
BAMRealigner::Report
realign(io::ReadReader::Path src, VcfReader::Path variants, io::ReadWriter::Path dst,
        const ReferenceGenome& reference, BAMRealigner::Config config);

} // namespace octopus

#endif

ext: hpp

max_stars_repo_name: Schaudge/octopus

max_issues_repo_name: Schaudge/octopus

max_forks_repo_name: Schaudge/octopus


text:
//          Copyright Gavin Band 2008 - 2012.
// 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)

#include <vector>
#include <unordered_map>
#include <ctime>
#include <boost/function.hpp>
#include <boost/format.hpp>
#include <boost/filesystem.hpp>
#include <boost/timer/timer.hpp>
#include <boost/functional/hash.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/discrete_distribution.hpp>
#include <boost/bind.hpp>

#include "config/package_revision_autogenerated.hpp"

#include "metro/constants.hpp"
#include "metro/ShotgunStochasticSearch.hpp"
#include "metro/regression/Design.hpp"
#include "metro/regression/Logistic.hpp"
#include "metro/IndependentParameterDistribution.hpp"
#include "metro/distributions/LogF.hpp"
#include "metro/distributions/Flat.hpp"
#include "metro/regression/LogUnnormalisedPosterior.hpp"
#include "metro/intersect_ranges.hpp"
#include "metro/fit_model.hpp"
#include "metro/CholeskyStepper.hpp"
#include "metro/SampleRange.hpp"
#include "metro/log_sum_exp.hpp"

#include "appcontext/CmdLineOptionProcessor.hpp"
#include "appcontext/ApplicationContext.hpp"

#include "genfile/MissingValue.hpp"
#include "genfile/bgen/View.hpp"
#include "genfile/bgen/Query.hpp"
#include "genfile/SNPDataSource.hpp"
#include "genfile/db/Error.hpp"
#include "genfile/ToGP.hpp"
#include "genfile/CrossCohortCovariateValueMapping.hpp"
#include "genfile/CommonSNPFilter.hpp"
#include "genfile/VariantIdentifyingDataFilteringSNPDataSource.hpp"

#include "qcdb/FlatFileMultiVariantOutputter.hpp"

namespace bfs = boost::filesystem ;

// #define DEBUG 1
// #define DEBUG_DOSAGESTORE 1

namespace globals {
	std::string const program_name = "multifinemap" ;
	std::string const program_version = package_version ;
	std::string const program_revision =  std::string( package_revision ).substr( 0, 7 ) ;
}

namespace {
	std::string print_state( std::size_t N, std::vector< std::size_t > const& states ) {
		std::string result = "[";
		std::size_t begin = 0 ;
		for( std::size_t i = 0; i <= states.size(); ++i ) {
			std::size_t end = (i==states.size()) ? N : states[i] ;
			for( std::size_t j = begin; j < end; ++j ) {
				result.push_back( '0' ) ;
			}
			if( end != N ) {
				result.push_back( '1' ) ;
			}
			begin = end+1 ;
		}
		result += ']' ;
		return result ;
	}
}

struct MFMOptions: public appcontext::CmdLineOptionProcessor
{
public:
	std::string get_program_name() const { return globals::program_name ; }

	void declare_options( appcontext::OptionProcessor& options ) {
		// Meta-options
		options.set_help_option( "-help" ) ;
		options.set_spec_option( "-spec" ) ;

		// File options
		options.declare_group( "Input file options" ) ;
	    options[ "-g" ]
	        .set_description( 	"Path to BGEN file giving genotypes." )
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 1 )
			.set_maximum_multiplicity( 1 ) ;
	    options[ "-filetype" ]
			.set_description(
				"Specify the filetype of the genotype files specified by -g. "
				"By default, qctool will guess the file type.  Use this option to override that guess. "
				"Possible types are: \"" + genfile::string_utils::join( genfile::SNPDataSource::get_file_types(), "\",\"" ) + "\"." )
			.set_takes_single_value()
			.set_default_value( "guess" ) ;

	    options[ "-assume-chromosome" ]
	        .set_description( "Chromosome to assume if not in genetic data." )
			.set_takes_single_value() ;
	    options[ "-s" ]
	        .set_description( "Path of sample file to input.  If specified, this option must occur as often as the -g option"
							" to specify one sample file per cohort." )
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 1 )
			.set_maximum_multiplicity( 1000 ) ;
		options[ "-outcome" ]
			.set_description( "Name of column containing outcome.  Currently must be a binary outcome." )
			.set_takes_single_value()
			.set_is_required() ;

		options.declare_group( "Input filtering options" ) ;
		options[ "-incl-range" ]
			.set_description( "Specify a range of SNPs (or comma-separated list of ranges of SNPs) to operate on. "
				"Each range should be in the format CC:xxxx-yyyy where CC is the chromosome and xxxx and yyyy are the "
				"start and end coordinates, or just xxxx-yyyy which matches that range from all chromosomes. "
				"You can also omit either of xxxx or yyyy to get all SNPs from the start or to the end of a chromosome." )
			.set_takes_values_until_next_option() ;
		options[ "-excl-range" ]
			.set_description( "Specify a range of SNPs (or comma-separated list of ranges of SNPs) to exclude from operation. "
				"Each range should be in the format CC:xxxx-yyyy where CC is the chromosome and xxxx and yyyy are the "
				"start and end coordinates, or just xxxx-yyyy which matches that range from all chromosomes. "
				"You can also omit either of xxxx or yyyy to get all SNPs from the start or to the end of a chromosome." )
			.set_takes_values_until_next_option() ;
		options[ "-incl-ranges" ]
			.set_description( "Like -incl-range, but read the ranges from the specified file."
				"The file must contain a whitespace-separated list of genomic ranges to include in the analysis." )
			.set_takes_values_until_next_option() ;
		options[ "-excl-ranges" ]
			.set_description( "Like -excl-range, but read the ranges from the specified file."
				"The file must contain a whitespace-separated list of genomic ranges to exclude from the analysis." )
			.set_takes_values_until_next_option() ;

		options.declare_group( "Output file options" ) ;
		options[ "-o" ]
			.set_description( "Path to output file.  Use \"-\" for stdout." )
			.set_takes_single_value()
			.set_default_value( "-" ) ;

		options.declare_group( "Model options" ) ;
		options[ "-incl-samples"]
			.set_description( "Filter out samples whose sample ID does not lie in the given file(s).")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-samples"]
			.set_description( "Filter out samples whose sample ID lies in the given file.")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-incl-samples-where"]
			.set_description( "Include samples by specifying conditions on the values of columns in the sample file.")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-excl-samples-where"]
			.set_description( "Exclude samples by specifying conditions on the values of columns in the sample file.")
			.set_takes_values( 1 )
			.set_minimum_multiplicity( 0 )
			.set_maximum_multiplicity( 100 ) ;
		options[ "-minimum-allele-count" ]
			.set_description( "Skip predictors that have less than this minor allele count" )
			.set_takes_values(1)
			.set_default_value( 10 ) ;
		options[ "-covariates" ]
			.set_description( "Specify the name of one or more covariates to include in the model."
				" These must be columns named in the sample file." )
			.set_takes_values_until_next_option() ;

		options.declare_group( "Model fitting options" ) ;
		options[ "-tolerance" ]
			.set_description( "Tolerance to use when model fitting" )
			.set_takes_values(1)
			.set_default_value( 0.001 ) ;
		options[ "-max-iterations" ]
			.set_description( "Maximum iterations to allow when model fitting" )
			.set_takes_values(1)
			.set_default_value( 100 ) ;
		options[ "-max-causal-variants" ]
			.set_description( "Set the maximum number of variants allowed in a causal configuration")
			.set_takes_values(1)
			.set_default_value(5)
		;
		options[ "-additive" ]
			.set_description( "Specify that only additive effects are allowed for" ) ;
		;
		options.declare_group( "Miscellaneous options" ) ;
		options[ "-analysis-name" ]
			.set_description( "Specify a name to label results from this analysis with." )
			.set_takes_single_value()
			.set_default_value( "qctool analysis" ) ;
		options[ "-analysis-chunk" ]
			.set_description( "Specify a name denoting the current genomic region or chunk on which this is run.  This is intended for use in parallel environments." )
			.set_takes_single_value()
			.set_default_value( genfile::MissingValue() ) ;
		options.declare_group( "Miscellaneous options" ) ;
		options[ "-debug" ]
			.set_description( "Output debugging information." ) ;
		options [ "-log" ]
			.set_description( "Specify that " + globals::program_name + " should write a log file to the given file." )
			.set_takes_single_value() ;
	}
} ;

struct DosageStore {
public:
	typedef boost::function< void( std::size_t const, boost::optional< std::size_t > const ) > ProgressCallback ;
	typedef std::unique_ptr< DosageStore > UniquePtr ;

	static UniquePtr create() {
		return UniquePtr( new DosageStore() ) ;
	}
public:
	DosageStore() {}

	void load( genfile::SNPDataSource::UniquePtr source, ProgressCallback callback ) {
		load_unsafe( *source, callback ) ;
	}

	std::size_t number_of_samples() const { return m_ploidy[0].size() ; }
	std::size_t number_of_variants() const { return m_variants.size() ; }
	genfile::VariantIdentifyingData const& variant( std::size_t i ) const { return m_variants[i] ; }

	std::size_t estimate_memory_usage() const {
		std::size_t result = 0 ;
		for( std::size_t i = 0; i < m_variants.size(); ++i ) {
			result += m_variants[i].estimate_bytes_used() ;
			result += m_ploidy[i].size() ;
			result += m_AB[i].size() ;
			result += m_BB[i].size() ;
		}
		return result ;
	}

	std::string get_summary() const {
		using boost::format ;
		return (
			boost::format( "DosageStore: %d variants, %d samples, approx %.2fMb memory usage.\n" )
				% m_variants.size()
				% (m_variants.size() > 0 ? m_AB[0].size() : 0 )
				% (estimate_memory_usage() / 1000000.0)
		).str() ;
	}

	genfile::VariantIdentifyingData const get_variant( std::size_t i ) const {
		assert( i < m_variants.size() ) ;
		return m_variants[i] ;
	}

	template< typename Callback >
	void get_dosages( std::size_t i, Callback callback ) const {
		assert( i < m_variants.size() ) ;
		double const constant = double(0xFFFFFFFF >> 24) ;
		assert( i < m_variants.size() ) ;
		for( std::size_t j = 0; j < m_AB[i].size(); ++j ) {
			// unpack 8 bit bgen encoding
#if DEBUG_DOSAGESTORE
			if( j < 10 ) {
				std::cerr << "get_dosages(): " << j << ": " << int( m_AB[i][j] ) << ", " << int( m_BB[i][j] ) << ".\n" ;
			}
#endif
			callback( j, double( m_AB[i][j] ) / constant, double( m_BB[i][j] ) / constant ) ;
		}
	}

	std::vector< metro::SampleRange > const nonmissing_samples( std::size_t i ) const {
		assert( i < m_variants.size() ) ;
		return m_nonmissing_samples[i] ;
	}

private:

	/* Store genotype probabilities in 8 bit, BGEN-like encoding */
	struct DosageSetter {
		DosageSetter(
			std::vector< uint8_t >* ploidy,
			std::vector< uint8_t >* AB,
			std::vector< uint8_t >* BB
		):
			m_ploidy( ploidy ),
			m_AB( AB ),
			m_BB( BB ),
			m_bits(8),
			m_constant( double(0xFFFFFFFF >> (32 - m_bits))),
			m_sample_i(0),
			m_order_type( genfile::eUnknownOrderType )
		{}

		void initialise( std::size_t nSamples, std::size_t nAlleles ) {
			assert( nAlleles == 2 ) ;
			m_ploidy->resize( nSamples ) ;
			std::fill( m_ploidy->begin(), m_ploidy->end(), 0x80 ) ;
			m_AB->resize( nSamples ) ;
			std::fill( m_AB->begin(), m_AB->end(), 0 ) ;
			m_BB->resize( nSamples ) ;
			std::fill( m_BB->begin(), m_BB->end(), 0 ) ;
			m_last_nonmissing_sample_i = 0 ;
		}

		bool set_sample( std::size_t n ) {
			m_sample_i = n ;
			m_order_type = genfile::eUnknownOrderType ;
			return true ;
		}

		void set_number_of_entries(
			uint32_t ploidy, std::size_t n,
			genfile::OrderType const order_type,
			genfile::ValueType const value_type
		) {
			assert( ploidy <= 2 ) ;
			assert( (order_type == genfile::ePerOrderedHaplotype || genfile::ePerUnorderedGenotype) && value_type == genfile::eProbability ) ;
			(*m_ploidy)[m_sample_i] = ploidy ;
			m_order_type = order_type ;
		}

		void set_value( std::size_t entry_i, genfile::MissingValue const value ) {
			uint8_t& ploidy = (*m_ploidy)[m_sample_i] ;
			uint8_t& AB = (*m_AB)[m_sample_i] ;
			uint8_t& BB = (*m_BB)[m_sample_i] ;
			ploidy |= 0x80 ;
			AB = BB = 0 ;
			record_missing_sample( m_sample_i ) ;
		}

		void set_value( std::size_t entry_i, double const value ) {
			uint8_t& ploidy = (*m_ploidy)[m_sample_i] ;
			uint8_t& AB = (*m_AB)[m_sample_i] ;
			uint8_t& BB = (*m_BB)[m_sample_i] ;

#if DEBUG_DOSAGESTORE
			if( m_sample_i < 10 ) {
				std::cerr << "DosageSetter::set_value(): " << m_sample_i << ", " << entry_i << ": " << value << ".\n" ;
			}
#endif
			if( !(ploidy & 0x80) ) {
				if( entry_i == 1 ) {
					AB = value * m_constant ;
				} else if (entry_i == 2 ) {
					BB = value * m_constant ;
				}
			}
		}

		void set_value( std::size_t entry_i, genfile::Integer const value ) {
			assert(0) ; // expecting GT field
		}

		void finalise() {
			// add the final sample range if needed.
			record_missing_sample( m_ploidy->size() ) ;
		}

		std::vector< metro::SampleRange > const nonmissing_samples() { return m_nonmissing_samples ; }

		private:
			std::vector< uint8_t >* m_ploidy ;
			std::vector< uint8_t >* m_AB ;
			std::vector< uint8_t >* m_BB ;
			uint32_t m_bits ;
			double m_constant ;
			std::size_t m_sample_i ;
			std::size_t m_last_nonmissing_sample_i ;
			genfile::OrderType m_order_type ;
			std::vector< metro::SampleRange > m_nonmissing_samples ;

		void record_missing_sample( std::size_t sample_i ) {
			if( sample_i >= m_last_nonmissing_sample_i ) {
				m_nonmissing_samples.push_back(
					metro::SampleRange( m_last_nonmissing_sample_i, m_sample_i )
				) ;
			}
			m_last_nonmissing_sample_i = sample_i + 1 ;
		}
	} ;

	void load_unsafe( genfile::SNPDataSource& source, ProgressCallback callback ) {
		if( source.total_number_of_snps() ) {
			std::size_t N = *source.total_number_of_snps() ;
			m_variants.reserve( N ) ;
			m_ploidy.reserve( N ) ;
			m_AB.reserve( N ) ;
			m_BB.reserve( N ) ;
		}

		genfile::VariantIdentifyingData variant ;
		std::vector< uint16_t > data ;
		std::size_t count = 0 ;
		std::vector< uint8_t > ploidy, AB, BB ;

		while( source.get_snp_identifying_data( &variant )) {
			// ok, read dosages.
			std::cerr << "LOADING " << variant << ".\n" ;
			DosageSetter setter( &ploidy, &AB, &BB ) ;
			source.read_variant_data()->get( ":genotypes:", genfile::to_GP_unphased( setter ) ) ;
			m_variants.push_back( variant ) ;
			m_nonmissing_samples.push_back( setter.nonmissing_samples() ) ;
			m_ploidy.push_back( ploidy ) ;
			m_AB.push_back( AB ) ;
			m_BB.push_back( BB ) ;
			//m_data.push_back( data ) ;
			if( callback ) {
				callback( ++count, source.total_number_of_snps() ) ;
			}
		}
	}

	std::vector< genfile::VariantIdentifyingData > m_variants ;
	// To avoid using too much space, we encode dosage data as 16-bit integer
	// and expand to required format when requested.
	std::vector< std::vector< metro::SampleRange > > m_nonmissing_samples ;
	std::vector< std::vector< uint8_t > > m_ploidy ;
	std::vector< std::vector< uint8_t > > m_AB ;
	std::vector< std::vector< uint8_t > > m_BB ;
} ;

struct MFMApplication: public appcontext::ApplicationContext {
	MFMApplication( int argc, char** argv ):
		appcontext::ApplicationContext(
			globals::program_name,
			"version: " + globals::program_version + ", revision: " + globals::program_revision,
			std::auto_ptr< appcontext::OptionProcessor >( new MFMOptions ),
			argc,
			argv,
			"-log"
		)
	{
		process() ;
	}

private:

	void process() {
		try {
			unsafe_process() ;
		} catch( genfile::InputError const& e ) {
			ui().logger() << "!! Error (" << e.what() << "): " << e.format_message() << ".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::FileNotFoundError const& e ) {
			ui().logger() << "\nError: No file matching \"" << e.filespec() << "\" could be found.\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
		catch( genfile::db::Error const& e ) {
			ui().logger() << "!! Error (" << e.what() << ") with the following statement: \""
				<< e.sql()
				<< "\".\n" ;
			throw appcontext::HaltProgramWithReturnCode( -1 ) ;
		}
	}

	void unsafe_process() {
		DosageStore::UniquePtr store = load_genotypes() ;

		genfile::CohortIndividualSource::UniquePtr
			samples = genfile::CohortIndividualSource::create( options().get< std::string >( "-s" ) ) ;

		metro::regression::LogLikelihood::UniquePtr
			ll = create_loglikelihood( samples->size() ) ;

		if( options().check( "-covariates" )) {
			add_covariates( ll->design(), *samples, options().get_values< std::string >( "-covariates" ) ) ;
		}

		set_outcome( *ll, *samples, options().get_value< std::string >( "-outcome" )) ;

		std::cerr << ll->design().get_summary() << "\n" ;

		run_search( *store, *ll ) ;
	}

	DosageStore::UniquePtr load_genotypes() {
		appcontext::UIContext::ProgressContext progress_context = ui().get_progress_context( "Loading data" ) ;
		DosageStore::UniquePtr store = DosageStore::create() ;
		store->load(
			open_snp_data_source( options().get< std::string >( "-g" ) ),
			progress_context
		) ;
		ui().logger() << store->get_summary() ;
		return store ;
	}

	genfile::SNPDataSource::UniquePtr
	open_snp_data_source( std::string const& filename ) const {
		genfile::Chromosome chromosome_hint ;
		std::string chromosome_indicator ;
		if( options().check_if_option_was_supplied( "-assume-chromosome" )) {
			chromosome_indicator = options().get_value< std::string >( "-assume-chromosome" ) ;
		}
		if( chromosome_indicator != "" ) {
			chromosome_hint = chromosome_indicator ;
		}

		genfile::SNPDataSource::UniquePtr source ;

		std::pair< std::string, std::string > uf = genfile::uniformise( filename ) ;

		{
			boost::optional< genfile::vcf::MetadataParser::Metadata > metadata ;
//			if( m_options.check( "-metadata" )) {
//				metadata = genfile::vcf::StrictMetadataParser(
//					m_options.get_value< std::string >( "-metadata" )
//				).get_metadata() ;
//			}
			source = genfile::SNPDataSource::create(
				filename,
				chromosome_hint,
				metadata,
				options().get< std::string >( "-filetype" )
			) ;
		}

		genfile::CommonSNPFilter::UniquePtr snp_filter = construct_snp_filter() ;
		// Filter SNPs if necessary
		if( snp_filter.get() ) {
			genfile::VariantIdentifyingDataFilteringSNPDataSource::UniquePtr snp_filtering_source
				= genfile::VariantIdentifyingDataFilteringSNPDataSource::create(
					source,
					genfile::VariantIdentifyingDataTest::UniquePtr( snp_filter.release() )
				) ;

			source.reset(
				snp_filtering_source.release()
			) ;
		}

		return source ;
	}

	genfile::CommonSNPFilter::UniquePtr construct_snp_filter() const {
		genfile::CommonSNPFilter::UniquePtr snp_filter ;

		if( options().check_if_option_was_supplied_in_group( "Input filtering options" )) {
			snp_filter.reset( new genfile::CommonSNPFilter ) ;
			if( options().check_if_option_was_supplied( "-incl-range" )) {
				std::vector< std::string > specs = options().get_values< std::string >( "-incl-range" ) ;
				for ( std::size_t i = 0; i < specs.size(); ++i ) {
					snp_filter->include_snps_in_range(
						genfile::GenomePositionRange::parse( specs[i] )
					) ;
				}
			}

			if( options().check_if_option_was_supplied( "-excl-range" )) {
				std::vector< std::string > specs = options().get_values< std::string >( "-excl-range" ) ;
				for ( std::size_t i = 0; i < specs.size(); ++i ) {
					snp_filter->exclude_snps_in_range(
						genfile::GenomePositionRange::parse( specs[i] )
					) ;
				}
			}

			if( options().check_if_option_was_supplied( "-incl-ranges" )) {
				std::vector< std::string > files = options().get_values< std::string >( "-incl-ranges" ) ;
				for( std::string const& filename: files ) {
					std::auto_ptr< std::istream > in = genfile::open_text_file_for_input( filename ) ;
					std::string range ;
					while( (*in) >> range ) {
						snp_filter->include_snps_in_range(
							genfile::GenomePositionRange::parse( range )
						) ;
					}
				}
			}

			if( options().check_if_option_was_supplied( "-excl-ranges" )) {
				std::vector< std::string > files = options().get_values< std::string >( "-excl-ranges" ) ;
				for( std::string const& filename: files ) {
					std::auto_ptr< std::istream > in = genfile::open_text_file_for_input( filename ) ;
					std::string range ;
					while( (*in) >> range ) {
						snp_filter->exclude_snps_in_range(
							genfile::GenomePositionRange::parse( range )
						) ;
					}
				}
			}
		}
		return snp_filter ;
	}

	metro::regression::LogLikelihood::UniquePtr create_loglikelihood( std::size_t number_of_samples ) {
		typedef std::vector< std::string > Names ;
		typedef std::vector< metro::SampleRange > SampleRanges ;
		using namespace metro::regression ;

		Eigen::MatrixXd outcome = Eigen::MatrixXd::Zero( number_of_samples, 2 ) ;

		Design::UniquePtr design = Design::create(
			outcome, SampleRanges(), Names({ "outcome=0", "outcome=1"} ),
			Names({
				"add1", "het1",
				"add2", "het2",
				"add3", "het3",
				"add4", "het4",
				"add5", "het5"
			})
		) ;
		LogLikelihood::UniquePtr ll( metro::regression::Logistic::create( design ).release() ) ;
		return ll ;
	}

	void add_covariates(
		metro::regression::Design& design,
		genfile::CohortIndividualSource const& samples,
		std::vector< std::string > const& covariates
	) {
		ui().logger() << "Adding covariates...\n" ;
		{
			std::set< std::string > uniqueCovariates( covariates.begin(), covariates.end() ) ;
			if( uniqueCovariates.size() != covariates.size() ) {
				throw genfile::BadArgumentError(
					"add_covariates()",
					"covariates=\"" + genfile::string_utils::join( covariates, " " ) + "\"",
					"Covariates should not be duplicated."
				) ;
			}
		}
		genfile::CohortIndividualSource::ColumnSpec const& spec = samples.get_column_spec() ;
		for( std::size_t i = 0; i < covariates.size(); ++i ) {
			std::string const& covariateName = covariates[i] ;
			std::size_t const column_index = spec.find_column( covariateName ) ;
			genfile::CohortIndividualSource::SingleColumnSpec columnSpec = spec[ column_index ] ;
			genfile::CrossCohortCovariateValueMapping::UniquePtr mapping
				= genfile::CrossCohortCovariateValueMapping::create( columnSpec, true ) ;
			mapping->add_source( samples ) ;
			if( columnSpec.is_discrete() ) {
				design.add_discrete_covariate(
					covariateName,
					[&samples, &column_index, &mapping] ( std::size_t sample_index ) {
						int value = -1 ;
						genfile::VariantEntry entry = samples.get_entry( sample_index, column_index ) ;
						if( !entry.is_missing() ) {
							value = mapping->get_mapped_value( entry ).as< int >() ;
						}
						return value ;
					},
					[&mapping]( int level ) {
						genfile::VariantEntry entry = mapping->get_unmapped_value( level ) ;
						return entry.as< std::string >() ;
					},
					mapping->get_number_of_distinct_mapped_values()
				) ;
			} else {
				design.add_single_covariate(
					covariateName,
					[&samples, &column_index, &mapping]( std::size_t sample_index ) {
						double value = metro::NA ;
						genfile::VariantEntry entry = samples.get_entry( sample_index, column_index ) ;
						if( !entry.is_missing() ) {
							value = mapping->get_mapped_value( entry ).as< double >() ;
						}
						return value ;
					}
				) ;
			}

			ui().logger()
				<< "++ Added covariate: \"" + covariateName + "\":\n"
				<<  mapping->get_summary( "     " )
				<< "\n\n" ;
		}
	}

	void set_outcome(
		metro::regression::LogLikelihood& ll,
		genfile::CohortIndividualSource const& samples,
		std::string const& column
	) const {
		typedef std::vector< std::string > Names ;
		typedef std::vector< metro::SampleRange > SampleRanges ;
		Eigen::MatrixXd outcome = Eigen::MatrixXd::Zero( samples.size(), 2 ) ;
		SampleRanges nonmissingness ;
		get_binary_outcome( samples, column, &outcome, &nonmissingness ) ;
		ll.design().set_outcome(
			outcome, nonmissingness,
			Names({ column + "=0", column + "=1" })
		) ;
	}

	void get_binary_outcome(
		genfile::CohortIndividualSource const& samples,
		std::string const& column,
		Eigen::MatrixXd* outcome,
		std::vector< metro::SampleRange >* nonmissingness
	) const {
		std::size_t column_i = samples.get_column_spec().find_column( column ) ;
		assert( samples.get_column_spec()[ column_i ].is_discrete() ) ;
		std::size_t last_nonmissing_sample_i = 0 ;
		outcome->setZero( samples.size(), 2 ) ;
		for( std::size_t i = 0; i < samples.size(); ++i ) {
			genfile::VariantEntry const& entry = samples.get_entry( i, column_i ) ;
			if( entry.is_missing() ) {
				nonmissingness->push_back( metro::SampleRange( last_nonmissing_sample_i, i )) ;
				last_nonmissing_sample_i = i+1 ;
			} else {
				int value = entry.as< int >() ;
				assert( value == 0 || value == 1 ) ;
				(*outcome)(i,value) = 1 ;
			}
		}
		if( last_nonmissing_sample_i < samples.size() ) {
			nonmissingness->push_back( metro::SampleRange( last_nonmissing_sample_i, samples.size() )) ;
		}
#if DEBUG
		std::cerr << "get_outcome(): outcome is:\n"
			<< outcome->block(0,0,10,2 ) << "\n"
			<< "nonmissingness: " << *nonmissingness << "\n" ;
#endif
	}

	typedef std::map< std::vector< metro::SampleRange >, std::pair< Eigen::VectorXd, double > > NullLLStore ;

	double test_variant(
		metro::ShotgunStochasticSearch::SelectedStates const& pick,
		DosageStore const& store,
		metro::regression::LogLikelihood& ll,
		NullLLStore& null_ll_store
	) {
#if DEBUG
		std::cerr << "   TESTING: " << print_state( store.number_of_variants(), pick ) << ".\n" ;
#endif
		using namespace metro::regression ;
		using metro::SampleRange ;

		double const minus_infinity = -std::numeric_limits< double >::infinity() ;
		double const log_weight = ( (pick.size() == 0) || (pick.size() > 5) )
			? minus_infinity
			: (pick.size() * std::log( 1.0 / store.number_of_variants() )) ;

		if( log_weight == minus_infinity ) {
			// Disallowed state, don't bother fitting.
			return minus_infinity ;
		}

		boost::format parameter_format( "%s%d/%s" ) ;
		bool const debug = options().check( "-debug" ) ;

		assert( pick.size() < 6 ) ;
		Eigen::MatrixXd predictors = Eigen::MatrixXd::Zero( store.number_of_samples(), 10 ) ;
		std::vector< metro::SampleRange > nonmissing_samples( 1, metro::SampleRange( 0, store.number_of_samples() )) ;

		metro::IndependentParameterDistribution prior( ll.get_parameter_names() ) ;
		metro::IndependentParameterDistribution null_prior( ll.get_parameter_names() ) ;

		double const add_prior_obs = 32.0 ;
		double const het_prior_obs = options().check( "-additive" ) ? 1000.0 : 32.0 ;

		for( std::size_t i = 0; i < 5; ++i ) {
			null_prior.set_prior(
				(parameter_format % "add" % (i+1) % ll.design().get_outcome_name(1) ).str(),
				metro::distributions::LogF::create( 10000.0, 10000.0 )
			) ;
			null_prior.set_prior(
				(parameter_format % "het" % (i+1) % ll.design().get_outcome_name(1) ).str(),
				metro::distributions::LogF::create( 10000.0, 10000.0 )
			) ;
		}

		for( std::size_t i = 0; i < pick.size(); ++i ) {
			store.get_dosages(
				pick[i],
				[&predictors,i] ( int sample, double ab, double bb ){
					predictors( sample, i*2+0) = ab + 2*bb ;
					predictors( sample, i*2+1) = ab ;
				}
			) ;
			nonmissing_samples = metro::impl::intersect_ranges(
				nonmissing_samples,
				store.nonmissing_samples(pick[i])
			) ;
			prior.set_prior(
				(parameter_format % "add" % (i+1) % ll.design().get_outcome_name(1) ).str(),
				metro::distributions::LogF::create( add_prior_obs, add_prior_obs )
			) ;
			prior.set_prior(
				(parameter_format % "het" % (i+1) % ll.design().get_outcome_name(1) ).str(),
				metro::distributions::LogF::create( het_prior_obs, het_prior_obs )
			) ;
		}
		for( std::size_t i = pick.size(); i < 5; ++i ) {
			prior.set_prior(
				(parameter_format % "add" % (i+1) % ll.design().get_outcome_name(1) ).str(),
				metro::distributions::LogF::create( add_prior_obs, add_prior_obs )
			) ;
			prior.set_prior(
				(parameter_format % "het" % (i+1) % ll.design().get_outcome_name(1) ).str(),
				metro::distributions::LogF::create( het_prior_obs, het_prior_obs )
			) ;
		}

		// Set up a trace for model fitting iterations.
		metro::CholeskyStepper::Tracer tracer ;
		if( debug ) {
			tracer = [this] (
				int iteration,
				double ll,
				double target_ll,
				metro::CholeskyStepper::Vector const& point,
				metro::CholeskyStepper::Vector const& first_derivative,
				metro::CholeskyStepper::Vector const& step,
				bool converged
			) {
				this->ui().logger()
					<< " ITERATION: " << iteration << "\n"
					<< "        LL: " << ll << "\n"
					<< "    TARGET: " << target_ll << "\n"
					<< "     POINT: " << point.transpose() << "\n"
					<< "DERIVATIVE: " << first_derivative.transpose() << "\n"
					<< "      STEP: " << step.transpose() << "\n"
					<< " CONVERGED: " << ( converged ? "yes" : "no" ) << "\n" ;
			} ;
		}

		// To compute the BF we need the null model.
		// The null model does not depend on the SNPs, but might depend on the missing samples.
		// We cache the value for speed if the nonmissing samples don't change.
		std::vector< std::string > comments ;
		// null_ll will be log-marginal likelihood of null model
		// i.e. logarithm of prior x likelihood, integrated over parameters.
		double null_ll = 0.0 ;
		Eigen::VectorXd null_parameters ;
		{
			NullLLStore::const_iterator where = null_ll_store.find( nonmissing_samples ) ;
			if( where == null_ll_store.end() ) {
				if( debug ) {
					this->ui().logger() << "---> fitting NULL...\n" ;
				}
				ll.design().set_predictors( Eigen::MatrixXd::Zero( store.number_of_samples(), 10 ), nonmissing_samples ) ;
				LogUnnormalisedPosterior posterior( ll, null_prior ) ;
				metro::CholeskyStepper stopping_condition( 0.01, 100, tracer ) ;
				std::pair< bool, int > fit = metro::fit_model(
					posterior,
					"null",
					Eigen::VectorXd::Zero( ll.identify_parameters().rows() ),
					stopping_condition,
					&comments
				) ;

				null_ll = laplace_approximate( posterior ) ;
				null_parameters = posterior.parameters() ;
				//null_ll = posterior.get_value_of_function() ;
				where = null_ll_store.insert(
					NullLLStore::value_type(
						nonmissing_samples,
						NullLLStore::mapped_type(
							ll.parameters(),
							null_ll
						)
					)
				).first ;
			} else {
				null_parameters = where->second.first ;
				null_ll = where->second.second ;
			}
		}

		ll.design().set_predictors( predictors, nonmissing_samples ) ;
		LogUnnormalisedPosterior posterior( ll, prior ) ;

		double posterior_weight = minus_infinity ;
		double result = minus_infinity ;

		if( log_weight > minus_infinity ) {
			if( debug ) {
				this->ui().logger() << "---> fitting full...\n" ;
			}
			metro::CholeskyStepper stopping_condition( 0.01, 100, tracer ) ;
			std::pair< bool, int > fit = metro::fit_model(
				posterior,
				"full",
				null_parameters,
				stopping_condition,
				&comments
			) ;
			if( fit.first ) {
				result = laplace_approximate( posterior ) - null_ll + log_weight ;
				//result = posterior.get_value_of_function() - null_ll + log_weight ;
			}

#if DEBUG
			std::cerr << print_state( store.number_of_variants(), pick ) << ".\n" ;
			//std::cerr << nonmissing_samples << ".\n" ;
			// std::cerr << "PREDICTORS:\n" << predictors.block( 0, 0, 10, predictors.cols() ) << ".\n" ;
			std::cerr << posterior.get_summary() << "\n" ;
			std::cerr << "  INCLUDED: " << posterior.ll().design().nonmissing_samples() << "\n" ;
			std::cerr << " CONVERGED: " << (fit.first ? "1" : "0") << "\n" ;
			std::cerr << "        LL: " << posterior.get_value_of_function() << "\n" ;
			std::cerr << "      NULL: " << null_ll << "\n" ;
			std::cerr << " pick.size: " << pick.size() << "\n" ;
			std::cerr << "LOG WEIGHT: " << log_weight << "\n" ;
			std::cerr << "    RESULT: " << result << "\n" ;

			std::cerr << "\n+++++++++++++++++++++++++\n" ;
#endif
		}

		return result ;
	}

	// Compute laplace approximation.
	// It is assumed that the function has already been maximised.
	double laplace_approximate( metro::SmoothFunction const& function ) const {
		double const fx = function.get_value_of_function() ;
		metro::SmoothFunction::Matrix const& H = function.get_value_of_second_derivative() ;
		double const log_2pi = std::log( 2.0 * 3.1415926535897932384626 ) ;
		Eigen::ColPivHouseholderQR< metro::SmoothFunction::Matrix > solver( -H ) ;
		return
			fx + 0.5 * ( function.number_of_parameters() * log_2pi - solver.logAbsDeterminant() ) ;
	}

	void run_search( DosageStore& store, metro::regression::LogLikelihood& ll ) {
		// Do a test search for now
		NullLLStore null_ll_store ;

		ui().logger() << "Running shotgun stochastic search...\n" ;

		metro::ShotgunStochasticSearch ss(
			store.number_of_variants(),
			[this, &store, &ll, &null_ll_store]( metro::ShotgunStochasticSearch::SelectedStates const& s ) {
				return this->test_variant( s, store, ll, null_ll_store ) ;
			},
			static_cast<std::uint32_t>(std::time(0))
		) ;

		for( std::size_t i = 0; i < 100; ++i ) {
			boost::timer::cpu_timer timer ;
			metro::ShotgunStochasticSearch::SelectedStates const& pick = ss.update() ;

			ui().logger() << "Iteration " << i << ": "
				<< "took " << boost::timer::format( timer.elapsed(), 2, "%ts" )
				<< ", state is: " << print_state( store.number_of_variants(), pick )
				<< "\n" ;
		}

		ui().logger() << "Outputting search results...\n" ;
		output_search_results( ss, store ) ;

		ui().logger() << "Fine-mapping complete.\n" ;
	}

	void output_search_results(
		metro::ShotgunStochasticSearch const& ss,
		DosageStore const& store
	) {
		qcdb::MultiVariantStorage::UniquePtr outputter
			= qcdb::MultiVariantStorage::create(
				options().get< std::string >( "-o" ),
				5,
				options().get< std::string >( "-analysis-name" ),
				options().get< std::string >( "-analysis-chunk" ),
				get_application_metadata()
		) ;
		outputter->add_variable( "loglikelihood" ) ;
		outputter->add_variable( "logposterior" ) ;
		outputter->add_variable( "posterior" ) ;

		std::vector< genfile::VariantIdentifyingData > variants ;

		auto comparator = []( double a, double b ) { return a > b ; } ;
		typedef std::multimap< double, metro::ShotgunStochasticSearch::SelectedStates, decltype(comparator) > OrderedStates ;
		OrderedStates ordered_states( comparator ) ;
		std::vector< double > posterior_weights ;
		ss.visited_states(
			[&ordered_states, &posterior_weights](
				metro::ShotgunStochasticSearch::SelectedStates const& states,
				double value
			) {
				std::cerr << "GOT: " << states.size() << ", " << value << "\n" ;
				ordered_states.insert( std::make_pair( value, states )) ;
				posterior_weights.push_back( value ) ;
			}
		) ;

		double const total_log_weight = metro::log_sum_exp( posterior_weights ) ;
		ui().logger() << "TOTAL LOG WEIGHT:" << total_log_weight << "\n" ;
		int count = 0 ;
		for( OrderedStates::const_iterator i = ordered_states.begin(); i != ordered_states.end(); ++i, ++count ) {
			metro::ShotgunStochasticSearch::SelectedStates const& states = i->second ;
			variants.resize( states.size() ) ;
			for( std::size_t i = 0; i < states.size(); ++i ) {
				variants[i] = store.variant( states[i] ) ;
			}
			outputter->store_data_for_key(
				variants,
				"rank",
				count
			) ;
			outputter->store_data_for_key(
				variants,
				"loglikelihood",
				i->first
			) ;
			outputter->store_data_for_key(
				variants,
				"logposterior",
				i->first - total_log_weight
			) ;
			outputter->store_data_for_key(
				variants,
				"posterior",
				std::exp( i->first - total_log_weight )
			) ;
		}
		outputter->finalise() ;
	}

#if 0
	void create_query( genfile::bgen::Query* query ) {
		if( options().check( "-incl-range" )) {
			auto const elts = options().get_values< std::string >( "-incl-range" ) ;
			for( std::string const& elt: elts ) {
				query->include_range( parse_range( elt )) ;
			}
		}
		if( options().check( "-excl-range" )) {
			auto const elts = options().get_values< std::string >( "-excl-range" ) ;
			for( std::string const& elt: elts ) {
				query->exclude_range( parse_range( elt )) ;
			}
		}
		if( options().check( "-incl-rsids" )) {
			auto const ids = options().get_values< std::string >( "-incl-rsids" ) ;
			query->include_rsids( ids ) ;
		}

		if( options().check( "-excl-rsids" )) {
			auto const ids = options().get_values< std::string >( "-excl-rsids" ) ;
			query->exclude_rsids( ids ) ;
		}
	}

	genfile::bgen::IndexQuery::GenomicRange parse_range( std::string const& spec ) const {
		std::size_t colon_pos = spec.find( ':' ) ;
		if ( colon_pos == std::string::npos ) {
			throw std::invalid_argument( "spec=\"" + spec + "\"" ) ;
		}

		std::vector< std::string > pieces ;
		pieces.push_back( spec.substr( 0, colon_pos )) ;
		pieces.push_back( spec.substr( colon_pos+1, spec.size() )) ;

		if( pieces.size() != 2 ) {
			throw std::invalid_argument( "spec=\"" + spec + "\"" ) ;
		}

		std::size_t separator_pos = pieces[1].find( '-' ) ;
		if ( separator_pos == std::string::npos ) {
			throw std::invalid_argument( "spec=\"" + spec + "\"" ) ;
		}

		std::string chromosome( pieces[0] ) ;
		int pos1 = (separator_pos == 0) ? 0 : std::stoi( pieces[1].substr( 0, separator_pos ) ) ;
		int pos2 = (separator_pos == (pieces[1].size()-1)) ? std::numeric_limits< int >::max() : std::stoi( pieces[1].substr( separator_pos + 1, pieces[1].size() ) ) ;
		assert( pos1 >= 0 ) ;
		assert( pos2 >= pos1 ) ;

		return genfile::bgen::IndexQuery::GenomicRange( chromosome, pos1, pos2 ) ;
	}
#endif
} ;

int main( int argc, char** argv ) {
	try {
		MFMApplication app( argc, argv ) ;
	}
	catch( appcontext::HaltProgramWithReturnCode const& e ) {
		return e.return_code() ;
	}
	return 0 ;
}

ext: cpp

max_stars_repo_name: gavinband/qctool

max_issues_repo_name: gavinband/qctool

max_forks_repo_name: gavinband/qctool


text:
/* openrtb_bid_request_test.cc
   Jeremy Barnes, 19 February 2013
   Copyright (c) 2013 Datacratic Inc.  All rights reserved.

   Test cases for the OpenRTB bid request parser.
*/

#define BOOST_TEST_MAIN
#define BOOST_TEST_DYN_LINK

#include <boost/test/unit_test.hpp>
#include "rtbkit/plugins/bid_request/openrtb_bid_request_parser.h"
#include "soa/types/json_parsing.h"
#include "rtbkit/openrtb/openrtb_parsing.h"
#include "jml/utils/filter_streams.h"

using namespace std;
using namespace ML;
using namespace Datacratic;
using namespace RTBKIT;

vector<string> samples = {
    "rtbkit/plugins/bid_request/testing/openrtb1_req.json",
    "rtbkit/plugins/bid_request/testing/openrtb2_req.json",
    "rtbkit/plugins/bid_request/testing/openrtb3_req.json",
    "rtbkit/plugins/bid_request/testing/openrtb4_req.json",
    "rtbkit/plugins/bid_request/testing/openrtb_wseat_req.json",
    "rtbkit/plugins/bid_request/testing/openrtb_banner.json",
    "rtbkit/plugins/bid_request/testing/openrtb_expandable_creative.json",
    "rtbkit/plugins/bid_request/testing/openrtb_mobile.json",
    "rtbkit/plugins/bid_request/testing/openrtb_video.json",
    "rtbkit/plugins/bid_request/testing/rubicon_banner1.json",
    "rtbkit/plugins/bid_request/testing/rubicon_banner2.json",
    "rtbkit/plugins/bid_request/testing/rubicon_banner3.json",
    "rtbkit/plugins/bid_request/testing/rubicon_banner4.json",
    "rtbkit/plugins/bid_request/testing/rubicon_desktop.json",
    "rtbkit/plugins/bid_request/testing/rubicon_mobile_app.json",
    "rtbkit/plugins/bid_request/testing/rubicon_mobile_web.json",
    "rtbkit/plugins/bid_request/testing/rubicon_test1.json"
};

vector<string> samples2_2 = {
    "rtbkit/plugins/bid_request/testing/openrtb_2_2_req_imp.json",
    "rtbkit/plugins/bid_request/testing/openrtb_2_2_req_video.json"
};

std::string loadFile(const std::string & filename)
{
    ML::filter_istream stream(filename);

    string result;

    while (stream) {
        string line;
        getline(stream, line);
        result += line + "\n";
    }

    return result;
}

void parseBidRequest(const std::string & filename)
{
    cerr << endl << "loading " << filename << endl;
    StreamingJsonParsingContext context;
    context.init(filename);

    OpenRTB::BidRequest req;
    DefaultDescription<OpenRTB::BidRequest> desc;
    desc.parseJson(&req, context);

    if (!req.unparseable.isNull())
        cerr << "unparseable:" << req.unparseable << endl;

    StreamJsonPrintingContext printContext(cout);
    desc.printJson(&req, printContext);
}

BOOST_AUTO_TEST_CASE( test_parse_openrtb_sample_requests )
{
    for (auto req: samples)
        parseBidRequest(req);
}

void testBidRequest(const std::string & filename, const std::string & version = "2.2")
{
    cerr << endl << "loading " << filename << endl;
    ML::Parse_Context context(filename);
    std::shared_ptr<OpenRTBBidRequestParser> p = OpenRTBBidRequestParser::openRTBBidRequestParserFactory(version);
    auto res = p->parseBidRequest(context, "test", "test");
    cerr << res->toJson() << endl;
}

BOOST_AUTO_TEST_CASE( test_openrtb_sample_requests )
{
    for (auto s: samples)
        testBidRequest(s);

    // Use 2.2 parsing
    for (auto s: samples2_2)
        testBidRequest(s, "2.2");
}

bool jsonDiff(const Json::Value & v1, const Json::Value & v2,
              bool oneOnly = false,
              string path = "")
{
    bool result = true;

    if (v1.type() != v2.type()) {
        cerr << path << ": different types: "
             << v1.toString() << " versus "
             << v2.toString() << endl;
        return false;
    }

    switch (v1.type()) {
    case Json::arrayValue:
        if (v1.size() != v2.size()) {
            cerr << path << ": differing lengths: " << v1.size() << " versus "
                 << v2.size() << endl;
            result = false;
            if (oneOnly) return false;
        }
        for (unsigned i = 0;  i < v1.size();  ++i) {
            if (!jsonDiff(v1[i], v2[i], oneOnly,
                          path + ML::format("[%d]", i))) {
                result = false;
                if (oneOnly) return false;
            }
        }
        return result;

    case Json::objectValue: {
        auto m1 = v1.getMemberNames();
        auto m2 = v2.getMemberNames();

        std::sort(m1.begin(), m1.end());
        std::sort(m2.begin(), m2.end());

        vector<string> intersection, only1, only2;
        std::set_intersection(m1.begin(), m1.end(),
                              m2.begin(), m2.end(),
                              back_inserter(intersection));
        std::set_difference(m1.begin(), m1.end(),
                            m2.begin(), m2.end(),
                            back_inserter(only1));
        std::set_difference(m2.begin(), m2.end(),
                            m1.begin(), m1.end(),
                            back_inserter(only2));

        if (!only1.empty()) {
            cerr << path << ": keys only present in first field: "
                 << only1 << endl;
            result = false;
            if (oneOnly) return false;
        }
        if (!only2.empty()) {
            cerr << path << ": keys only present in second field: "
                 << only1 << endl;
            result = false;
            if (oneOnly) return false;
        }

        for (auto f: intersection) {
            if (!jsonDiff(v1[f], v2[f], oneOnly, path + "." + f)) {
                result = false;
                if (oneOnly) return false;
            }
        }
        return result;
    }

    default:
        if (v1.toString() != v2.toString()) {
            cerr << path << ": difference: " << v1 << " vs " << v2 << endl;
            result = false;
            if (oneOnly) return false;
        }
    }
    return result;
}

void testBidRequestRoundTrip(const std::string & filename,
                             const std::string & reqStr)
{
    // Take an OpenRTB bid request; convert to Datacratic format;
    // serialize that to JSON; reconstitute from JSON back into the
    // original format

    static DefaultDescription<OpenRTB::BidRequest> desc;

    OpenRTB::BidRequest req;

    ML::Parse_Context c(filename);
    {
        StreamingJsonParsingContext context;
        context.init(filename, reqStr.c_str(), reqStr.size());
        desc.parseJson(&req, context);
    }

    string printed;
    {
        std::ostringstream stream;
        StreamJsonPrintingContext printContext(stream);
        desc.printJson(&req, printContext);
        printed = stream.str();
    }

    OpenRTB::BidRequest req2;

    {
        StreamingJsonParsingContext context2;
        context2.init(filename, reqStr.c_str(), reqStr.size());
        desc.parseJson(&req2, context2);
    }

    string printed2;
    {
        std::ostringstream stream;
        StreamJsonPrintingContext printContext(stream);
        desc.printJson(&req2, printContext);
        printed2 = stream.str();
    }

    BOOST_CHECK_EQUAL(printed, printed2);

    std::shared_ptr<OpenRTBBidRequestParser> p = OpenRTBBidRequestParser::openRTBBidRequestParserFactory("2.1");

    // Convert to a standard bid request
    std::unique_ptr<BidRequest> br(p->parseBidRequest(c, "test", "test"));
    // Convert it to JSON
    string s1 = br->toJsonStr();

    std::unique_ptr<BidRequest> br2(BidRequest::parse("rtbkit", s1));

    string s2 = br2->toJsonStr();

    if (s1 != s2) {
        return;
        Json::Value j = br->toJson();
        Json::Value j2 = br2->toJson();
        BOOST_CHECK(jsonDiff(j, j2));
    }
}

void testBidRequestConversion(const std::string &fileName, const std::string &request)
{
    static DefaultDescription<OpenRTB::BidRequest> desc;

    // Parse the BidRequest in OpenRTB format
    OpenRTB::BidRequest originalReq;

    {
        StreamingJsonParsingContext context;
        context.init(fileName, request.c_str(), request.size());
        desc.parseJson(&originalReq, context);
    }

    string jsonBr1;
    {
        std::ostringstream stream;
        StreamJsonPrintingContext printContext(stream);
        desc.printJson(&originalReq, printContext);
        jsonBr1 = stream.str();
    }

    ML::Parse_Context c(fileName);

    std::shared_ptr<OpenRTBBidRequestParser> p = OpenRTBBidRequestParser::openRTBBidRequestParserFactory("2.1");

    // Convert to a standard bid request
    std::unique_ptr<BidRequest> br(p->parseBidRequest(c, "test", "test"));

    // Convert it back to OpenRTB
    OpenRTB::BidRequest br2 = p->toBidRequest(*br);

    string jsonBr2;
    {
        std::ostringstream stream;
        StreamJsonPrintingContext printContext(stream);
        desc.printJson(&br2, printContext);
        jsonBr2 = stream.str();
    }

    auto json = [](const string &jsonStr) {
        Json::Reader reader;
        Json::Value value;
        bool ok;
        ok = reader.parse(jsonStr, value);
        BOOST_CHECK(ok);
        return value;
    };

    BOOST_CHECK(jsonDiff(json(jsonBr1), json(jsonBr2)));
}

BOOST_AUTO_TEST_CASE( test_openrtb_round_trip )
{
    vector<string> reqs;

    for (auto s: samples)
        reqs.push_back(loadFile(s));

    for (unsigned i = 0;  i < reqs.size();  ++i) {
        testBidRequestRoundTrip(samples[i], reqs[i]);
        //testBidRequestConversion(samples[i], reqs[i]);
    }
}

BOOST_AUTO_TEST_CASE( benchmark_openrtb_round_trip )
{
    vector<string> reqs;

    for (auto s: samples)
        reqs.push_back(loadFile(s));

    int done = 0;

    Date before = Date::now();

    for (unsigned i = 0;  i < 1000;  ++i) {

        for (unsigned i = 0;  i < reqs.size();  ++i, ++done)
            testBidRequestRoundTrip(samples[i], reqs[i]);
    }

    double elapsed = Date::now().secondsSince(before);

    cerr << "did " << done << " in " << elapsed << "s at "
         << done / elapsed << "/s" << endl;
}

BOOST_AUTO_TEST_CASE( benchmark_openrtb_parsing )
{
    cerr << "benchmarking OpenRTB parsing" << endl;

    vector<string> reqs;

    for (auto s: samples)
        reqs.push_back(loadFile(s));

    int done = 0;

    Date before = Date::now();

    DefaultDescription<OpenRTB::BidRequest> desc;

    for (unsigned i = 0;  i < 1000;  ++i) {

        for (unsigned i = 0;  i < reqs.size();  ++i, ++done) {

            OpenRTB::BidRequest req;

            {
                StreamingJsonParsingContext context;
                context.init(samples[i], reqs[i].c_str(), reqs[i].size());
                desc.parseJson(&req, context);
            }
        }
    }

    double elapsed = Date::now().secondsSince(before);

    cerr << "did " << done << " in " << elapsed << "s at "
         << done / elapsed << "/s" << endl;
}

BOOST_AUTO_TEST_CASE( benchmark_openrtb_conversion )
{
    cerr << "benchmarking OpenRTB parsing and conversion" << endl;

    vector<string> reqs;

    for (auto s: samples)
        reqs.push_back(loadFile(s));

    int done = 0;

    Date before = Date::now();

    DefaultDescription<OpenRTB::BidRequest> desc;

    std::shared_ptr<OpenRTBBidRequestParser> p = OpenRTBBidRequestParser::openRTBBidRequestParserFactory("2.1");

    for (unsigned i = 0;  i < 1000;  ++i) {

        for (unsigned i = 0;  i < reqs.size();  ++i, ++done) {

            OpenRTB::BidRequest req;

            {
                StreamingJsonParsingContext context;
                context.init(samples[i], reqs[i].c_str(), reqs[i].size());
                std::unique_ptr<BidRequest> br(p->parseBidRequest(*context.context, "openrtb", "openrtb"));

            }
        }
    }

    double elapsed = Date::now().secondsSince(before);

    cerr << "did " << done << " in " << elapsed << "s at "
         << done / elapsed << "/s" << endl;
}

BOOST_AUTO_TEST_CASE( benchmark_canonical_parsing )
{
    cerr << "benchmarking canonical parsing of OpenRTB-derived bid requests" << endl;

    DefaultDescription<OpenRTB::BidRequest> desc;

    vector<string> reqs;

    std::shared_ptr<OpenRTBBidRequestParser> p = OpenRTBBidRequestParser::openRTBBidRequestParserFactory("2.1");

    for (auto s: samples) {
        OpenRTB::BidRequest req;
        {
            StreamingJsonParsingContext context;
            context.init(s);
            std::unique_ptr<BidRequest> br(p->parseBidRequest(*context.context, "openrtb", "openrtb"));
            reqs.push_back(br->toJsonStr());
        }
    }

    int done = 0;

    Date before = Date::now();

    for (unsigned i = 0;  i < 1000;  ++i) {

        for (unsigned i = 0;  i < reqs.size();  ++i, ++done) {
            std::unique_ptr<BidRequest> br2(BidRequest::parse("rtbkit", reqs[i]));
        }
    }

    double elapsed = Date::now().secondsSince(before);

    cerr << "did " << done << " in " << elapsed << "s at "
         << done / elapsed << "/s" << endl;
}

BOOST_AUTO_TEST_CASE( id_provider ) {

    cerr << "id provider test : making sure we parse it correctly and always set it" << endl;

    DefaultDescription<OpenRTB::BidRequest> desc;

    std::shared_ptr<OpenRTBBidRequestParser> p = OpenRTBBidRequestParser::openRTBBidRequestParserFactory("2.1");

    vector<string> reqs;
    bool hasBuyerUid, hasUserId, hasDeviceUA, hasDeviceIP;

    for (auto s: samples) {

        hasBuyerUid = false;
        hasUserId = false;
        hasDeviceUA = false;
        hasDeviceIP = false;

        OpenRTB::BidRequest req;
        std::unique_ptr<BidRequest> datacraticReq;
        {
            StreamingJsonParsingContext context;
            context.init(s);
            datacraticReq.reset(p->parseBidRequest(*context.context, "openrtb", "openrtb"));
        }

        OpenRTB::BidRequest req2;
        {
            StreamingJsonParsingContext context;
            context.init(s);
            desc.parseJson(&req2, context);
        }

        if(req2.user) {
            if(req2.user->buyeruid)
                hasBuyerUid = true;
            // Check if we have userId
            else if(req2.user->id)
                hasUserId = true;
        }

        if(req2.device) {

            // Check if we have the IP
            if(!(req2.device->ip.empty()))
                hasDeviceIP = true;

            // Check if we have the UA
            if(!(req2.device->ua.empty()))
                hasDeviceUA = true;
        }

        // Check if we manage to to get buyeruid
        if(hasBuyerUid) {
            cerr << "BidRequest Buyer Id : " << req2.user->buyeruid << endl;
            cerr << "Datacratic BidRequest Provider Id : " << datacraticReq->userIds.providerId << endl;
            BOOST_CHECK( req2.user->buyeruid == datacraticReq->userIds.providerId );
        }
        // Check what we have for providerId in the bidrequest.
        else if(hasUserId) {
            cerr << "BidRequest User Id : " << req2.user->id << endl;
            cerr << "Datacratic BidRequest Provider Id : " << datacraticReq->userIds.providerId << endl;
            BOOST_CHECK( req2.user->id == datacraticReq->userIds.providerId );
        }
        else if(hasDeviceIP && hasDeviceUA) {
            cerr << "BidRequest Device IP : " << req2.device->ip << endl;
            cerr << "BidRequest Device Agent : " << req2.device->ua << endl;
            cerr << "Datacratic BidRequest Provider Id : " << datacraticReq->userIds.providerId << endl;
            BOOST_CHECK( datacraticReq->userIds.providerId != Id(0) );
        }
        else {
            cerr << "BidRequest User Id : " << req2.user->id << endl;
            cerr << "Datacratic BidRequest Provider Id : " << datacraticReq->userIds.providerId << endl;
            BOOST_CHECK( datacraticReq->userIds.providerId == Id(0) );
        }
    }
}

ext: cc

max_stars_repo_name: etnrlz/rtbkit

max_issues_repo_name: TuanTranEngineer/rtbkit

max_forks_repo_name: TuanTranEngineer/rtbkit


text:
/* ----------------------------------------------------------------------------

 * GTSAM Copyright 2010, Georgia Tech Research Corporation,
 * Atlanta, Georgia 30332-0415
 * All Rights Reserved
 * Authors: Frank Dellaert, et al. (see THANKS for the full author list)

 * See LICENSE for the license information

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

/**
 * @brief wraps PinholeCamera classes to python
 * @author Ellon Paiva Mendes (LAAS-CNRS)
 **/

#include <boost/python.hpp>

#define NO_IMPORT_ARRAY
#include <numpy_eigen/NumpyEigenConverter.hpp>

#include "gtsam/geometry/PinholeCamera.h"
#include "gtsam/geometry/Cal3_S2.h"

using namespace boost::python;
using namespace gtsam;

typedef PinholeBaseK<Cal3_S2> PinholeBaseKCal3_S2;

// Wrapper on PinholeBaseK<Cal3_S2> because it has pure virtual method calibration()
struct PinholeBaseKCal3_S2Callback : PinholeBaseKCal3_S2, wrapper<PinholeBaseKCal3_S2>
{
  const Cal3_S2 & calibration () const {
    return this->get_override("calibration")();
  }
};

BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(project_overloads, PinholeBaseKCal3_S2::project, 1, 4)
BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(range_overloads, PinholeBaseKCal3_S2::range, 1, 3)

// Function pointers to disambiguate project() calls
Point2  (PinholeBaseKCal3_S2::*project1) (const Point3 &pw, OptionalJacobian< 2, 6 > Dpose, OptionalJacobian< 2, 3 > Dpoint, OptionalJacobian< 2, FixedDimension<Cal3_S2>::value > Dcal) const = &PinholeBaseKCal3_S2::project;
Point2  (PinholeBaseKCal3_S2::*project2) (const Unit3 &pw,  OptionalJacobian< 2, 6 > Dpose, OptionalJacobian< 2, 2 > Dpoint, OptionalJacobian< 2, FixedDimension<Cal3_S2>::value > Dcal) const = &PinholeBaseKCal3_S2::project;

// function pointers to disambiguate range() calls
double (PinholeBaseKCal3_S2::*range1) (const Point3 &point, OptionalJacobian< 1, 6 > Dcamera, OptionalJacobian< 1, 3 > Dpoint) const = &PinholeBaseKCal3_S2::range;
double (PinholeBaseKCal3_S2::*range2) (const Pose3 &pose,   OptionalJacobian< 1, 6 > Dcamera, OptionalJacobian< 1, 6 > Dpose)  const = &PinholeBaseKCal3_S2::range;
double (PinholeBaseKCal3_S2::*range3) (const CalibratedCamera &camera, OptionalJacobian< 1, 6 > Dcamera, OptionalJacobian< 1, 6 > Dother) const = &PinholeBaseKCal3_S2::range;

// wrap projectSafe in a function that returns None or a tuple
// TODO(frank): find out how to return an ndarray instead
object project_safe(const PinholeBaseKCal3_S2& camera, const gtsam::Point3& p) {
  auto result = camera.projectSafe(p);
  if (result.second)
    return make_tuple(result.first.x(), result.first.y());
  else
    return object();
}

void exportPinholeBaseK() {
  class_<PinholeBaseKCal3_S2Callback, boost::noncopyable>("PinholeBaseKCal3_S2", no_init)
      .def("calibration", pure_virtual(&PinholeBaseKCal3_S2::calibration),
           return_value_policy<copy_const_reference>())
      .def("project_safe", make_function(project_safe))
      .def("project", project1, project_overloads())
      .def("project", project2, project_overloads())
      .def("backproject", &PinholeBaseKCal3_S2::backproject)
      .def("backproject_point_at_infinity", &PinholeBaseKCal3_S2::backprojectPointAtInfinity)
      .def("range", range1, range_overloads())
      .def("range", range2, range_overloads())
      .def("range", range3, range_overloads());
}

max_stars_repo_name: karamach/gtsam

max_issues_repo_name: yfcube/gtsam

max_forks_repo_name: yfcube/gtsam


text:
/*!
@file
Forward declares `boost::hana::Group`.

@copyright Louis Dionne 2015
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE.md or copy at http://boost.org/LICENSE_1_0.txt)
 */

#ifndef BOOST_HANA_FWD_GROUP_HPP
#define BOOST_HANA_FWD_GROUP_HPP

#include <boost/hana/detail/std/forward.hpp>
#include <boost/hana/fwd/core/datatype.hpp>
#include <boost/hana/fwd/core/operators.hpp>

namespace boost { namespace hana {
    //! @ingroup group-concepts
    //! The `Group` concept represents `Monoid`s where all objects have
    //! an inverse w.r.t. the `Monoid`'s binary operation.
    //!
    //! A [Group][1] is an algebraic structure built on top of a `Monoid`
    //! which adds the ability to invert the action of the `Monoid`'s binary
    //! operation on any element of the set. Specifically, a `Group` is a
    //! `Monoid` `(S, +)` such that every element `s` in `S` has an inverse
    //! (say `s'`) which is such that
    //! @code
    //!     s + s' == s' + s == identity of the Monoid
    //! @endcode
    //!
    //! There are many examples of `Group`s, one of which would be the
    //! additive `Monoid` on integers, where the inverse of any integer
    //! `n` is the integer `-n`. The method names used here refer to
    //! exactly this model.
    //!
    //!
    //! Superclass
    //! ----------
    //! `Monoid`
    //!
    //!
    //! Laws
    //! ----
    //! For all objects `x` of a `Group` `G`, the following laws must be
    //! satisfied:
    //! @code
    //!     plus(x, negate(x)) == zero<G>() // right inverse
    //!     plus(negate(x), x) == zero<G>() // left inverse
    //! @endcode
    //!
    //!
    //! Minimal complete definitions
    //! ----------------------------
    //! 1. `minus`\n
    //! When `minus` is specified, the `negate` method is defaulted by setting
    //! @code
    //!     negate(x) = minus(zero<G>(), x)
    //! @endcode
    //!
    //! 2. `negate`\n
    //! When `negate` is specified, the `minus` method is defaulted by setting
    //! @code
    //!     minus(x, y) = plus(x, negate(y))
    //! @endcode
    //!
    //!
    //! Provided models
    //! ---------------
    //! 1. For non-boolean arithmetic data types\n
    //! A data type `T` is arithmetic if `std::is_arithmetic<T>::%value` is
    //! true. For a non-boolean arithmetic data type `T`, a model of `Group`
    //! is automatically defined by setting
    //! @code
    //!     minus(x, y) = (x - y)
    //!     negate(x) = -x
    //! @endcode
    //!
    //! @note
    //! The rationale for not providing a Group model for `bool` is the same
    //! as for not providing a `Monoid` model.
    //!
    //!
    //! Operators
    //! ---------
    //! For convenience, the following operators are provided as an
    //! equivalent way of calling the corresponding method:
    //! @code
    //!     binary -  ->  minus
    //!     unary -   ->  negate
    //! @endcode
    //!
    //!
    //! Structure-preserving functions
    //! ------------------------------
    //! Let `A` and `B` be two `Group`s. A function `f : A -> B` is said to
    //! be a [Group morphism][2] if it preserves the group structure between
    //! `A` and `B`. Rigorously, for all objects `x, y` of data type `A`,
    //! @code
    //!     f(plus(x, y)) == plus(f(x), f(y))
    //! @endcode
    //! Because of the `Group` structure, it is easy to prove that the
    //! following will then also be satisfied:
    //! @code
    //!     f(negate(x)) == negate(f(x))
    //!     f(zero<A>()) == zero<B>()
    //! @endcode
    //! Functions with these properties interact nicely with `Group`s, which
    //! is why they are given such a special treatment.
    //!
    //!
    //! [1]: http://en.wikipedia.org/wiki/Group_(mathematics)
    //! [2]: http://en.wikipedia.org/wiki/Group_homomorphism
    struct Group { };

    //! Subtract two elements of a group.
    //! @relates Group
    //!
    //! Specifically, this performs the `Monoid` operation on the first
    //! argument and on the inverse of the second argument, thus being
    //! equivalent to:
    //! @code
    //!     minus(x, y) == plus(x, negate(y))
    //! @endcode
    //!
    //!
    //! Cross-type version of the method
    //! --------------------------------
    //! The `minus` method is "overloaded" to handle distinct data types
    //! with certain properties. Specifically, `minus` is defined for
    //! _distinct_ data types `A` and `B` such that
    //! 1. `A` and `B` share a common data type `C`, as determined by the
    //!    `common` metafunction
    //! 2. `A`, `B` and `C` are all `Group`s when taken individually
    //! 3. `to<C> : A -> B` and `to<C> : B -> C` are `Group`-embeddings, as
    //!    determined by the `is_embedding` metafunction.
    //!
    //! The definition of `minus` for data types satisfying the above
    //! properties is obtained by setting
    //! @code
    //!     minus(x, y) = minus(to<C>(x), to<C>(y))
    //! @endcode
    //!
    //!
    //! Example
    //! -------
    //! @snippet example/group.cpp minus
#ifdef BOOST_HANA_DOXYGEN_INVOKED
    constexpr auto minus = [](auto&& x, auto&& y) -> decltype(auto) {
        return tag-dispatched;
    };
#else
    template <typename T, typename U, typename = void>
    struct minus_impl;

    struct _minus {
        template <typename X, typename Y>
        constexpr decltype(auto) operator()(X&& x, Y&& y) const {
            return minus_impl<
                typename datatype<X>::type, typename datatype<Y>::type
            >::apply(
                detail::std::forward<X>(x),
                detail::std::forward<Y>(y)
            );
        }
    };

    constexpr _minus minus{};
#endif

    //! Return the inverse of an element of a group.
    //! @relates Group
    //!
    //!
    //! Example
    //! -------
    //! @snippet example/group.cpp negate
#ifdef BOOST_HANA_DOXYGEN_INVOKED
    constexpr auto negate = [](auto&& x) -> decltype(auto) {
        return tag-dispatched;
    };
#else
    template <typename G, typename = void>
    struct negate_impl;

    struct _negate {
        template <typename X>
        constexpr decltype(auto) operator()(X&& x) const {
            return negate_impl<typename datatype<X>::type>::apply(
                detail::std::forward<X>(x)
            );
        }
    };

    constexpr _negate negate{};
#endif

    template <>
    struct operators::of<Group>
        : decltype(minus), decltype(negate)
    { };
}} // end namespace boost::hana

#endif // !BOOST_HANA_FWD_GROUP_HPP

max_stars_repo_name: josephwinston/hana

max_issues_repo_name: josephwinston/hana

max_forks_repo_name: josephwinston/hana


text: #include <boost/date_time/time_duration.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text: #include <boost/log/keywords/facility.hpp>