/** negative log likelihood of negative binomial with mean and tau 
\brief Negative binomial with mean=mu and variance = mu*tau
\author Mollie Brooks
\param x observed counts
\param mu is the predicted mean
\param tau is the overdispersion parameter like in the quasi-poisson. should be >1
\return negative log likelihood \f$ -( \ln(\Gamma(x+k))-\ln(\Gamma(k))-\ln(x!)+k\ln(k)+x\ln(\mu)-(k+x)\ln(k+\mu) )\f$
where \f$ k=\mu/(10^{-120}+\tau-1.0) \f$
\ingroup STATLIB
**/
dvariable dnbinom_tau(const dvector& x, const dvar_vector& mu, const dvar_vector& tau)
{
	//the observed counts are in x
	//mu is the predicted mean
	//tau is the overdispersion parameter
	RETURN_ARRAYS_INCREMENT();
	int i,imin,imax;
	imin=x.indexmin();
	imax=x.indexmax();
	dvariable loglike;
	loglike=0.;

	for(i = imin; i<=imax; i++)
	{
		if (value(tau(i))<1.0)
		{
			cerr<<"tau("<<i<<") is <=1.0 in dnbinom_tau()";
			return(0.0);
		}
		
		loglike += log_negbinomial_density(x(i), mu(i), tau(i));
	}
	RETURN_ARRAYS_DECREMENT();
	return(-loglike);
}

/**
 * Description not yet available.
 * \param
 */
void lvector::fill_multinomial(const int& seed, const dvector& p)
  // Fils a dvector with random numbers drawn from a multinomial distribution
  {
    double sum=mean(p)*p.size();
    int pmin=p.indexmin();
    int pmax=p.indexmax();
    dvector tmp(pmin,pmax);
    dvector tmp1(pmin,pmax);
    dvector choose(indexmin(),indexmax());
    choose.fill_randu(seed);
    tmp=p/sum;
    tmp1(pmin)=tmp(pmin);
    for (int j=pmin+1;j<=pmax-1;j++)
    {
      tmp1(j)=tmp1(j-1)+tmp(j);
    }
    tmp1(pmax)=1.0;

    for (int i=indexmin();i<=indexmax();i++)
    {
      int j=pmin;
      while (choose(i)>tmp1(j))
      {
        j++;
      }
      (*this)(i)=j;
    }
  }

/**
 * Description not yet available.
 * \param
 */
dvector solve_trans(const lower_triangular_dmatrix& M,const dvector& y)
{
  int mmin=M.indexmin();
  int mmax=M.indexmax();

  if (y.indexmin() !=mmin || y.indexmax() !=mmax)
  {
    cerr << "incompatible size in solve_trans" << endl;
    ad_exit(1);
  }
  dvector x(mmin,mmax);
  int i,j;

  for (i=mmax;i>=mmin;i--)
  {
    double sum=0.0;
    for (j=i+1;j<=mmax;j++)
    {
      sum+=M(j,i)*x(j);
    }
    x(i)=(y(i)-sum)/M(i,i);
  }

  return x;
}

/**
Determine if the lower and upper bounds of two evctors match in a specified
function.
\param v1 a data vector
\param v2 a data vector
\param function_nam a pointer to the name of the function in question.
*/
void shape_check(const dvector& v1, const dvector& v2,
  const char *function_name)
{
  if (v1.indexmin() != v2.indexmin() || v1.indexmax() != v2.indexmax())
  {
    cerr << " Vector sizes do no match in" << function_name << "\n";
    ad_exit(1);
  }
}

/**
Return the total sum of the elements in values.
\param values dvector
*/
double sum(const dvector& values)
{
  double total = 0.0;
  for (int i = values.indexmin(); i <= values.indexmax(); ++i)
  {
    total += values.elem(i);
  }
  return total;
}

/**
Return dvector results of squaring elements in a values;
constant vector object.

\ingroup misc
\param values of constant object to be squared.
\return vector of the same length as #values containing \f$values_i^2\f$
*/
dvector square(const dvector& values)
{
  dvector results;
  results.allocate(values);
  for (int i = values.indexmin(); i <= values.indexmax(); ++i)
  {
    results(i) = square(values(i));
  }
  return results;
}

matrix::matrix(const dvector &vec, bool isRow)
	: data(vector<dvector>(isRow? 1: vec.size())), colNum(isRow? vec.size(): 1) {
		if(isRow) {
			data.front() = vec;
		} else {
			for(long i = 0; i < (signed) vec.size(); i++) {
				data.at(i).resize(1);
				data.at(i).front() = vec.at(i);
			}
		}
}

dvector cubic_spline_function::operator () (const dvector& u)
{
  int mmin=u.indexmin();
  int mmax=u.indexmax();
  dvector z(mmin,mmax);
  for (int i=mmin;i<=mmax;i++)
  {
    z(i)=splint(x,y,y2,u(i));
  }
  return z;
}

/**
Return ivector filled with flags for 1 positive and -1 negative
for values in dvector v.
*/
ivector sgn(const dvector& v)
{
  int mmin = v.indexmin();
  int mmax = v.indexmax();
  ivector ret(mmin, mmax);
  for (int i = mmin; i <= mmax; ++i)
  {
    ret(i)= v(i) > 0 ? 1 : -1;
  }
  return ret;
}

dvariable robust_regression(dvector& obs, dvar_vector& pred, 
  const double& cutoff) 
{
  if (obs.indexmin() != pred.indexmin() || obs.indexmax() != pred.indexmax() )
  {
    cerr << "Index limits on observed vector are not equal to the Index\n\
 limits on the predicted vector in robust_reg_likelihood function\n";
  }
  RETURN_ARRAYS_INCREMENT(); //Need this statement because the function
			     //returns a variable type
  int min=obs.indexmin();
  int max=obs.indexmax();
  dvariable sigma_hat;
  dvariable sigma_tilde; 
  int nobs = max-min+1;
  double width=3.0;
  double pcon=0.05;
  double width2=width*width;
  dvariable zpen;
  zpen=0.0;
  double a,a2;
  a=cutoff; 
     // This bounds a between 0.05 and 1.75
  a2=a*a;
  dvariable tmp,tmp2,tmp4,sum_square,v_hat;
  dvar_vector diff_vec = obs-pred;
  tmp = norm(diff_vec);
  sum_square = tmp * tmp;
  v_hat = 1.e-80 + sum_square/nobs;
  sigma_hat=pow(v_hat,.5);
  sigma_tilde=a*sigma_hat;
  double b=2.*pcon/(width*sqrt(PI));
  dvariable log_likelihood;
  dvariable div1;
  dvariable div2,div4;
  div1 = 2*(a2*v_hat);
  div2 = width2*(a2*v_hat);
  div4 = div2*div2;
  log_likelihood = 0;
  for (int i=min; i<=max; i++)
  {
    tmp=diff_vec[i];
    tmp2=tmp*tmp;
    tmp4=tmp2*tmp2;
    log_likelihood -= log((1-pcon)*exp(-tmp2/div1)+b/(1.+tmp4/div4) );
  }
  log_likelihood += nobs*log(a2*v_hat)/2.;
  log_likelihood += zpen;
  RETURN_ARRAYS_DECREMENT(); // Need this to decrement the stack increment
			     // caused by RETURN_ARRAYS_INCREMENT();
  return(log_likelihood);  
}

/**
 * Description not yet available.
 * \param
 */
dmatrix outer_prod(const dvector& v1, const dvector& v2)
{
    dmatrix tmp(v1.indexmin(),v1.indexmax(), v2.indexmin(), v2.indexmax() );

    for (int i=v1.indexmin(); i<=v1.indexmax(); i++)
    {
        for (int j=v2.indexmin(); j<=v2.indexmax(); j++)
        {
            tmp.elem(i,j)=v1.elem(i)*v2.elem(j);
        }
    }
    return(tmp);
}

/**
Construct ivector with same dimensions as u.
*/
ivector::ivector(const dvector& u)
{
  allocate(u);
  for (int i=indexmin();i<=indexmax();i++)
  {
#ifdef OPT_LIB
    elem(i) = static_cast<int>(u.elem(i));
#else
    double ui = u.elem(i);
    assert(ui <= INT_MAX);
    elem(i) = static_cast<int>(ui);
#endif
  }
}

  bool brightRGB::getMax(const image& img,dvector& dest) const{

    // image empty?
    if (img.empty()) {
      setStatusString("image empty");
      dest.resize(0);
      return false;
    }

    const rgbPixel transColor = getParameters().transColor;
    ivector maxV(3,-1);
    image::const_iterator it = img.begin();
    if(getParameters().transparent) {
      while(it != img.end()) {
	if(*it != transColor) {
	  if((*it).getRed() > maxV.at(0))
	    maxV.at(0) = (*it).getRed();
	  if((*it).getGreen() > maxV.at(1))
	    maxV.at(1) = (*it).getGreen();
	  if((*it).getBlue() > maxV.at(2))
	    maxV.at(2) = (*it).getBlue();
	}
	it++;
      }
      // only transparent pixels?
      if (maxV.at(0)==-1) {
        setStatusString("only transparent pixels");
        dest.resize(0);
        return false;
      }
    } else { // no transparent color
      while(it != img.end()) {
	if((*it).getRed() > maxV.at(0))
	  maxV.at(0) = (*it).getRed();
	if((*it).getGreen() > maxV.at(1))
	  maxV.at(1) = (*it).getGreen();
	if((*it).getBlue() > maxV.at(2))
	  maxV.at(2) = (*it).getBlue();
	it++;
      }
    }
    if(maxV.at(0) == -1)
      return false;
    dest.castFrom(maxV);
    // normalize to 0..1
    dest.divide(255);
    return true;
  };

/**
 * Description not yet available.
 * \param
 */
double ghk(const dvector& lower,const dvector& upper,const dmatrix& Sigma,
  const dmatrix& eps,int _i)
{
  int n=lower.indexmax();
  dmatrix ch=choleski_decomp(Sigma);
  dvector l(1,n);
  dvector u(1,n);

  ghk_test(eps,_i); // test for valid i range
  double weight=1.0;
  int k=_i;
  {
    l=lower;
    u=upper;
    for (int j=1;j<=n;j++)
    {
      l(j)/=ch(j,j);
      u(j)/=ch(j,j);
      double Phiu=cumd_norm(u(j));
      double Phil=cumd_norm(l(j));
      weight*=Phiu-Phil;
      double eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil);
      for (int i=j+1;i<=n;i++)
      {
        double tmp=ch(i,j)*eta;
        l(i)-=tmp;
        u(i)-=tmp;
      }
    }
  }
  return weight;
}

  bool MLP::calcGradient(const dmatrix& inputs,
                         const ivector& ids,
                         dvector& grad) {

    if (inputs.rows() != ids.size()) {
      setStatusString("Number of vectors not consistent with number of ids");
      return false;
    }

    dvector tmp;
    int i;
    double tmpError;

    totalError = 0;
    calcGradient(inputs.getRow(0),ids.at(0),grad);
    computeActualError(ids.at(0),totalError);

    for (i=1;i<inputs.rows();++i) {
      calcGradient(inputs.getRow(i),ids.at(i),tmp);
      computeActualError(ids.at(i),tmpError);
      grad.add(tmp);
      totalError+=tmpError;
    }

    return true;
  }

void param_init_bounded_matrix_vector::set_scalefactor(const dvector& s)
{
  int mmin=indexmin();
  int mmax=indexmax();
  if (s.indexmin()!=mmin || s.indexmax() != mmax)
  {
    cerr << "non matching vector bounds in"
     " init_bounded_matrix_vector::set_scalefactor" << endl;
    ad_exit(1);
  }

  for (int i=mmin;i<=mmax;i++)
  {
    (*this)(i).set_scalefactor(s(i));
  }
}

    inline void copy_from_host(dvector<bool> &out, const std::vector<bool> &in)
    {
        std::vector<unsigned char> temp;
        temp.assign(in.begin(), in.end());

        out.copy_from_host((const bool *)&temp[0], temp.size());
    }

Particle::Particle(const dvector& pos, const prob_t fitness)
{
    mPos = pos;
    mVel.resize(pos.size());
    mPBestPos = pos;
    mFitness  = fitness;
    mPBestFitness = fitness;
}

/**
 * Description not yet available.
 * \param
 */
 dvar_vector::dvar_vector(const dvector& t)
 {
   if (!t)
   {
     allocate();
   }
   else
   {
     va=NULL;
     allocate(t.indexmin(),t.indexmax());
     initialize();
     for (int i = indexmin(); i <= indexmax(); i++)
     {
       va[i].x=(t.v)[i];
     }
   }
 }

   /*
    * compute the error using the last propagated input and the given
    * pattern
    */
  bool MLP::computePatternError(const int id,
                                const dvector& outUnits,
                                double& error) const {

    const int lastIdx = outUnits.size();
    int j;

    double tmp;
    error = 0.0;
    for (j=0;j<lastIdx;++j) {
      tmp = (outUnits.at(j)-((j==id)?on:off));
      error += tmp*tmp;
    }

    error *= 0.5;

    return true;
  }