inline void f_harmonic(valarray<double> &r, valarray<double> &k, const int dim)

{
    //constants
    const double m = 1.0;
    const double omega = 1.0;

    //get number of particles
    const int N_p = r.size() / (2 * dim);


    //reset k vector
    for (unsigned int i = 0; i < k.size(); i++)
    {
        k[i] = 0.0;
    }



    for (int d = 0; d < dim; d++)
    {

        for (int i = 0; i < N_p; i++)

        {
            k[(2 * i + 1)*dim + d] = (-1.) * m * omega * omega * r[i * 2 * dim + d];
            k[i * 2 * dim + d] = r[(2 * i + 1) * dim + d] / m;
        }

    }


}

T norm(const valarray<T>& v, const unsigned p) {
   T res = 0, scale = 0., absvi, ssq = 1., tmp;

   if ( v.size() < 1 ) {
      return static_cast<T>(0.);
   }

   else if ( v.size() == 1 ) {
      return abs(v[0]);
   }

   else {
      switch (p) {
         case 0:
            //for (size_t i = 0; i < v.size(); ++i) {
            //   if ( res < abs(v[i]) )
            //      res = abs(v[i]);
            //}
            //return res;
            for (auto it = begin(v); it != end(v); ++it) {
               if ( res < abs(*it) )
                  res = abs(*it);
            }
            return res;
            //return MAX(abs(v.min()), abs(v.max()));

         case 1:
            for (size_t i = 0; i < v.size(); ++i) {
               res += abs(v[i]);
            }
            return res;
            //for ( auto it = begin(v); it != end(v); ++it)
            //   res += abs(*it);
            //return res;

         case 2:
            for (size_t i = 0; i < v.size(); ++i) {
               if ( v[i] != 0. ) {
                  absvi = abs(v[i]);
                  if ( scale < absvi ) {
                     tmp = scale/absvi;
                     ssq = 1. + ssq*tmp*tmp;
                     scale = absvi;
                  }
                  else {
                     tmp = absvi/scale;
                     ssq += tmp*tmp;
                  }
               }
            }
            res = scale*sqrt(ssq);
            return res;

         default:
            cerr << "error: utils:norm(vector): unsupported p-norm: p = " << p 
                 << endl;
            exit(-1);
      }  
   }
}

	void SetUpCL(valarray<cl_float> & data, valarray<cl_float> & data_err, valarray<cl_float> & model, valarray<cl_float> & output)
	{
		unsigned int test_size = data.size();

		assert(data_err.size() == test_size);
		assert(model.size() == test_size);
		assert(output.size() == test_size);

		// Init OpenCL and the routine
		cl = new COpenCL(CL_DEVICE_TYPE_GPU);
		zero = new CRoutine_Zero(cl->GetDevice(), cl->GetContext(), cl->GetQueue());
		zero->SetSourcePath(LIBOI_KERNEL_PATH);
		zero->Init();
		square = new CRoutine_Square(cl->GetDevice(), cl->GetContext(), cl->GetQueue());
		square->SetSourcePath(LIBOI_KERNEL_PATH);
		square->Init();
		r = new CRoutine_Chi(cl->GetDevice(), cl->GetContext(), cl->GetQueue(), zero, square);
		r->SetSourcePath(LIBOI_KERNEL_PATH);
		r->Init(test_size);

		// Make OpenCL buffers for the data, data_err, model, and output.
		data_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);
		data_err_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);
		model_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);
		output_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);

		// Fill the input buffer
		int err = CL_SUCCESS;
		err = clEnqueueWriteBuffer(cl->GetQueue(), data_cl, CL_TRUE, 0, sizeof(cl_float) * test_size, &data[0], 0, NULL, NULL);
		err = clEnqueueWriteBuffer(cl->GetQueue(), data_err_cl, CL_TRUE, 0, sizeof(cl_float) * test_size, &data_err[0], 0, NULL, NULL);
		err = clEnqueueWriteBuffer(cl->GetQueue(), model_cl, CL_TRUE, 0, sizeof(cl_float) * test_size, &model[0], 0, NULL, NULL);
	}

void apply_op(valarray<T>& lhs, valarray<T> const& x, valarray<T> const& y, Op op = Op{}) {
	uint64_t size = std::min(x.size(), y.size());
	size = std::min(size, lhs.size()); // probably not needed
	for (uint64_t k = 0; k < size; k += 1) {
		lhs[k] = op(x[k], y[k]);
	}
}

	void MakeChiZeroBuffers(valarray<cl_float> & data, valarray<cl_float> & data_err, valarray<cl_float> & model, valarray<cl_float> & output, unsigned int test_size)
	{
		unsigned int n = 2*test_size;
		// Create buffers
		data.resize(n);
		data_err.resize(n);
		model.resize(n);
		output.resize(n);

		valarray<cl_float2> temp = CModel::GenerateUVSpiral_CL(test_size);

		// Set data = model to produce a zero chi result.
		for(int i = 0; i < test_size; i++)
		{
			data[i] = temp[i].s0;
			data[test_size + i] = temp[i].s1;

			model[i] = temp[i].s0;
			model[test_size + i] = temp[i].s1;

			// 1% error on amplitudes, 10% error on phases
			data_err[i] = amp_err * data[i];
			data_err[test_size + i] = phi_err * data[i];
		}
	}

static double compareAgainst(
            const valarray<double> &z,     // column vector
            const valarray<double> &h,     // column vector
            const valarray<double> &Q,     // symmetric matrix Q(x)
            const valarray<double> &invQ   // inverse of Q
            )
{
    const int n = z.size();

    assert( Q.size() == n * n );
    assert( h.size() == n );

    //assert( hasPosDet( dmatQ ) );

    valarray<double> Q2PI( n * n );
    
    Q2PI = Q * 2.0 * PI;

    // Compute det(Q2PI) = b * 2^c .
    double b;
    long c;
    det( Q2PI, n, &b, &c );

    // Compute log(det(Q2PI)).
    double dTerm1 = log( b ) + c * log( 2.0 );

    valarray<double>  residual = z - h;

    valarray<double>  term2 = multiply( multiply( transpose( residual, n ), n, invQ, n ), n, residual, 1 );
    double        dTerm2 = term2[0];

    return (dTerm1 + dTerm2) / 2.0;
}

// return dy/dx for each particular equation
valarray<double> dydx( const double x, const valarray<double>& y ) {

    if ( y.size() != 4 ) throw;

    valarray<double> tmpdydx(0.0,y.size());
    try {
        // slices would be easier, but oh well...
        double z = y[0];
        double p = y[1];
        double zbar = y[2];
        double pbar = y[3];
        
        double zdot = ( 1 + z*zbar ) * ( pbar + (pbar*zbar)*z );
        double pdot = 2 * zbar * ( p*pbar + (p*z)*(pbar*zbar) ) +
            ( 1 + z*zbar ) * ( (pbar*zbar)*z );
        double zbardot = ( 1 + z*zbar ) * ( p + (p*z)*zbar );
        double pbardot = 2 * z * ( p*pbar + (p*z)*(pbar*zbar) ) +
            ( 1 + z*zbar ) * ( (p*z)*zbar );

        // recombine into tmpdydx... again, slices would be easier
        tmpdydx[0] = zdot;
        tmpdydx[1] = pdot;
        tmpdydx[2] = zbardot;
        tmpdydx[3] = pbardot;

    }
    catch(out_of_range) {
        cerr << "oops" << endl;
        exit(1);
    }

    return tmpdydx;

}

inline void rk4(valarray<double> &r, double h, double &t, const int dim, void (*f)(valarray<double> &, valarray<double> &, const int))

{
    //define k, r_step, and temp arrays to hold eom evaluations
    valarray<double> k(r.size());
    valarray<double> r_step(r.size());
    valarray<double> temp(r.size());

    const double half_h = h / 2.;

    //1st rk4 step
    f(r, k, dim);
    r_step = h * (1. / 6.) * k;
    temp = r + half_h * k;

    //2nd
    f(temp, k, dim);
    r_step += h * (1. / 3.) * k;
    temp = r + half_h * k;

    //3rd
    f(temp, k, dim);
    r_step += h * (1. / 3.) * k;
    temp = r + h * k;

    //4th
    f(temp, k, dim);

    //advance r in time
    r += r_step + h * (1. / 6.) * k;


    //advance time
    t += h;
}

BroadEstepper1::BroadEstepper1(CanaryOptions& opts,
                               valarray<double> vals_in,
                               CanaryPrior & prior_in,
                               vector<int> & clusters_in)
  : _opts(opts)
  {
  N = vals_in.size();
  G = clusters_in.size();
  _vals.resize(vals_in.size()); _vals = vals_in;
  _prior = prior_in;
  _cvec.resize(G); _cvec = clusters_in;
  _prop.resize(G,1.0/G);
  _mean.resize(G);
  _var.resize(G);

  for (unsigned int k=0; k<G; k++)
    {
    int cpos = _cvec[k];
//  The algorithm converges to a cycle dependent on initial values.  If it
//  converged to a point - as it should - it would make sense to uncomment.
//    _prop[k] = _prior.prop()[cpos];
    _mean[k] = _prior.mean()[cpos];
    _var[k] = _prior.var()[cpos];
    }

  _prob.ReSize(N,G);
  _prob = 0.0;

  update_prob();
  }

double FullCovariance::calcLogdetFromCholesky( const valarray<double>& chol ) const
{
  //------------------------------------------------------------
  // Preliminaries.
  //------------------------------------------------------------

  int nChol = static_cast<int>(sqrt(static_cast<double>(chol.size())));


  //------------------------------------------------------------
  // Validate the current state (debug version only).
  //------------------------------------------------------------

  assert( chol.size() == nChol * nChol );


  //------------------------------------------------------------
  // Compute the logdet.
  //------------------------------------------------------------

  // Sum the logs of the diagonals.
  double logsum = 0.0;
  for ( int i = 0; i < nChol; i++ )
  {
    logsum += log( choleskyCached[i + i * nChol] );
  }

  // Return
  //                      T                       2         ---
  //     logdet( Chol Chol  )  =  log[ det( Chol )  ]  =  2 >    log( chol   )  .
  //                                                        ---           ii
  return 2.0 * logsum;
}

vector<double> convert(const valarray<double>& v1)
{
  vector<double> v2(v1.size());
  for(int i=0;i<v1.size();i++)
    v2[i] = v1[i];
  return v2;
}

double GradientProjection::computeSteepestDescentVector(
        valarray<double> const &b,
        valarray<double> const &x,
        valarray<double> &g) const {
    // find steepest descent direction
    //  g = 2 ( b - A x )
    //    where: A = denseQ + sparseQ
    //  g = 2 ( b - denseQ x) - 2 sparseQ x
    //
    //  except the 2s don't matter because we compute 
    //  the optimal stepsize anyway
    COLA_ASSERT(x.size()==b.size() && b.size()==g.size());
    g = b;
    for (unsigned i=0; i<denseSize; i++) {
        for (unsigned j=0; j<denseSize; j++) {
            g[i] -= (*denseQ)[i*denseSize+j]*x[j];
        }
    }
    // sparse part:
    if(sparseQ) {
        valarray<double> r(x.size());
        sparseQ->rightMultiply(x,r);
        g-=r;
    }
    return computeStepSize(g,g);
}

int drawCairo(const string& fname,
        const valarray<double>& Xin, const valarray<double>& Yin, 
        const Hull& hull) {
#ifdef CAIRO_HAS_SVG_SURFACE
    unsigned n=Xin.size();
    assert(Yin.size()==n);

    // normalise coords to range 0-1
    valarray<double> X=Xin, Y=Yin;
    X-=X.min();
    Y-=Y.min();
    X/=X.max();
    Y/=Y.max();

    Cairo::RefPtr<Cairo::SvgSurface> surface =
        Cairo::SvgSurface::create(fname, width+2*border, height+2*border);

    Cairo::RefPtr<Cairo::Context> cr = Cairo::Context::create(surface);

    cr->save(); // save the state of the context
    cr->set_source_rgba(0.0, 0.0, 0.0, 0.7);
    // draw a circle at each coordinate
    for(unsigned i=0;i<n;i++) {
        dot(cr,xcoord(X[i]),ycoord(Y[i]));
    }

    cr->set_source_rgba(0.0, 0.0, 0.0, 0.3);
    cr->move_to(xcoord(X[hull[0]]),ycoord(Y[hull[0]]));
    for(unsigned i=1;i<hull.size();i++) {
        cr->line_to(xcoord(X[hull[i]]),ycoord(Y[hull[i]]));
    }
    cr->line_to(xcoord(X[hull[0]]),ycoord(Y[hull[0]]));
    cr->stroke();
    cr->set_source_rgba(0.0, 0.0, 0.0, 1.);
    for(vector<unsigned>::const_iterator i=hull.begin();i!=hull.end();++i) {
        unsigned j=*i;
        stringstream ss;
        ss<<j;
        printf("p[%d]=(%f,%f)\n",j,X[j],Y[j]);
        cr->move_to(xcoord(X[j]),ycoord(Y[j]));
        cr->show_text(ss.str());
        cr->stroke();
    }
    cr->restore();

    cr->show_page();

    cout << "Wrote SVG file \"" << fname << "\"" << endl;
    return 0;

#else

    cout << "You must compile cairo with SVG support for this example to work."
        << endl;
    return 1;

#endif

}

  /// FIXME - this could be faster by a factor of 2 - we are sorting TWICE!
  pair<double,double> confidence_interval(const valarray<double>& values,double P) 
  {
    vector<double> values2(values.size());
    for(int i=0;i<values.size();i++)
      values2[i] = values[i];

    return confidence_interval(values2,P);
  }

valarray<int> iterate_zs(valarray<cmplx>& zs, const complex<float>& c,
                         size_t max_iters) {
    valarray<int> ns(zs.size());
#pragma omp parallel for schedule(runtime)
    for (size_t i = 0; i < zs.size(); i++)
        ns[i] = iterate_z(zs[i], c, max_iters);
    return ns;
}

/**
 * generates a random set of n points in X and Y.
 */
void randTest(unsigned n, valarray<double>& X, valarray<double>& Y) {
    X.resize(n);
    Y.resize(n);
    srand(time(NULL));
	for(unsigned i=0;i<n;i++) {
		X[i]=getRand(1.);
		Y[i]=getRand(1.);
	}
}

valarray<double> outer_prod(valarray<double> x, valarray<double> y) {
    valarray<double> result(x.size()*y.size());
    for(unsigned j = 0; j < x.size(); j++) {
        for(unsigned i = 0; i < y.size(); i++) {
            result[j*y.size() + i] = x[j]*y[i];
        }
    }
    return result;
}

valarray<int> iterate_zs(valarray<cmplx>& zs, const complex<double>& c,
                         size_t max_iters) {
    valarray<int> ns(zs.size());
#pragma omp parallel for default(none) shared(zs, c, max_iters, ns) \
                         schedule(runtime)
    for (size_t i = 0; i < zs.size(); i++)
        ns[i] = iterate_z(zs[i], c, max_iters);
    return ns;
}

const valarray<double> multiply( const valarray<double>& X, int nColsX, 
                                 const valarray<double>& Y, int nColsY )
{
  if( X.size() == 0 || Y.size() == 0 )
    return valarray<double>(0);

  //
  // Row Major Matrices
  //
  //    C = alpha * A * B + beta * C  --- (1)
  //
  // A : m by k         lda (stride) = k
  // B : k by n         ldb (stride) = n
  // C : m by n         ldc (stride) = n
  //
  //
  // Column Major Matrices
  //
  //    Z = alpha * X * Y + beta * C  --- (2)
  //    Z = C^t 
  //      = alpha * B^t * A^t + beta * C^t  --- (3)
  //
  // X = B^t : n by k   ldx (stride) = n
  // Y = A^t : k by m   ldy (stride) = k
  // Z = C^t : n by m   ldz (stride) = n
  //
   const int m = nColsY;
   const int k = nColsX;
   const int n = X.size() / k; assert( n*k == X.size() );
 
   valarray<double> Z( n * m );
   const double *pX = &X[0];
   const double *pY = &Y[0];
   double *pZ = &(Z[0]);

   int lda = n;
   int ldb = k;
   int ldc = n;
  
   cblas_dgemm( CblasColMajor, 
		CblasNoTrans,
		CblasNoTrans,
		n,   // #rows of the first matrix
		m,   // #cols of the second matrix
		k,   // #cols of the first which is equal to the second matrix
		ALPHA,
		pX,  // B^t
		lda,   // stride of B^t = #rows in B^t = #cols of B
		pY,  // A^t
		ldb,   // stride of A^t = #rows of A^t = #cols of A
		BETA,
		pZ,  // C^t
		ldc    // stride of C^t = #rows of C^t = #cols of C
              );
   return Z;
}

double CMinimizerThread::ComputeChi2r(valarray<double> & chis, unsigned int n_params)
{
	// Square the numerator, divide by the uncertainties
	chis *= chis;

	double chi2_sum = chis.sum();

	// Now compute the sum dividied by (n_data - n_params - 1)
	return chi2_sum / (chis.size() - n_params - 1);
}