void MolecularDynamics::assign_velocities(Float temperature) {
  ParticleIndexes ips = get_simulation_particle_indexes();
  setup_degrees_of_freedom(ips);
  ParticlesTemp ps = IMP::internal::get_particle(get_model(), ips);

  boost::normal_distribution<Float> mrng(0., 1.);
  boost::variate_generator<RandomNumberGenerator &,
                           boost::normal_distribution<Float> >
      sampler(random_number_generator, mrng);

  for (ParticlesTemp::iterator iter = ps.begin(); iter != ps.end();
       ++iter) {
    Particle *p = *iter;
    LinearVelocity(p).set_velocity(algebra::Vector3D(sampler(), sampler(),
                                                     sampler()));
  }

  Float rescale =
      sqrt(temperature / get_kinetic_temperature(get_kinetic_energy()));

  for (ParticlesTemp::iterator iter = ps.begin(); iter != ps.end();
       ++iter) {
    Particle *p = *iter;
    LinearVelocity v(p);

    algebra::Vector3D velocity = v.get_velocity();
    velocity *= rescale;
    v.set_velocity(velocity);
  }
}

void MolecularDynamics::assign_velocities(Float temperature)
{
  ParticleIndexes ips=get_simulation_particle_indexes();
  setup_degrees_of_freedom(ips);
  ParticlesTemp ps= IMP::internal::get_particle(get_model(), ips);

  boost::normal_distribution<Float> mrng(0., 1.);
  boost::variate_generator<RandomNumberGenerator&,
                           boost::normal_distribution<Float> >
      sampler(random_number_generator, mrng);

  for (ParticlesTemp::iterator iter = ps.begin();
       iter != ps.end(); ++iter) {
    Particle *p = *iter;

    for (int i = 0; i < 3; ++i) {
      p->set_value(vs_[i], sampler());
    }
  }

  Float rescale = sqrt(temperature/
                  get_kinetic_temperature(get_kinetic_energy()));

  for (ParticlesTemp::iterator iter = ps.begin();
       iter != ps.end(); ++iter) {
    Particle *p = *iter;

    for (int i = 0; i < 3; ++i) {
      Float velocity = p->get_value(vs_[i]);
      velocity *= rescale;
      p->set_value(vs_[i], velocity);
    }

  }
}

inline float sample_chisq (rng_t & rng, float nu)
{
    // HACK <float> appears to be broken in libstdc++ 4.6
    //typedef std::chi_squared_distribution<float> chi_squared_distribution_t;
    typedef std::chi_squared_distribution<double> chi_squared_distribution_t;

    chi_squared_distribution_t sampler(nu);
    return sampler(rng);
}

inline float sample_gamma (
        rng_t & rng,
        float alpha,
        float beta = 1.f)
{
    // HACK <float> appears to be broken in libstdc++ 4.6
    //typedef std::gamma_distribution<float> gamma_distribution_t;
    typedef std::gamma_distribution<double> gamma_distribution_t;

    gamma_distribution_t sampler(alpha, beta);
    return sampler(rng);
}

int Cosisim::full_ik( RandomPath begin, RandomPath end,
                      Progress_notifier* progress_notifier ) {
  typedef FullIK_cdf_estimator< geostat_utils::CoKrigingConstraints,
                                geostat_utils::CoKrigingCombiner, 
                                MarkovBayesCovariance > FullIkEstimator;
  FullIkEstimator cdf_estimator( marginal_->p_begin(),marginal_->p_end(),
                                 covariances_.begin(), covariances_.end(),
                                 combiners_.begin(), combiners_.end(),
                                 *kconstraints_ );

  // set up the sampler
  Random_number_generator gen;
  Monte_carlo_sampler_t< Random_number_generator > sampler( gen );

  
  // this vector will contain the indicies of the indicators to be estimated
  // at a given grid node
  std::vector< unsigned int > unestimated_indicators;
  
  for( ; begin != end ; ++begin ) {
    if( !progress_notifier->notify() ) return 1;

    get_current_local_cdf( *begin, unestimated_indicators );
    if( unestimated_indicators.empty() ) {
      sampler( *begin, *ccdf_ ); 
      continue;
    }
   
    geostat_utils::NeighborhoodVector neighborhoods =
      full_ik_find_neighbors( *begin, unestimated_indicators );

    DEBUG_PRINT_LOCATION( "center", begin->location() );
    DEBUG_PRINT_NEIGHBORHOOD( "hard1", &neighborhoods[0] );
    DEBUG_PRINT_NEIGHBORHOOD( "soft1", &neighborhoods[1] );
    DEBUG_PRINT_NEIGHBORHOOD( "hard2", &neighborhoods[2] );
    DEBUG_PRINT_NEIGHBORHOOD( "soft2", &neighborhoods[3] );

    GsTLPoint watch( 52, 95, 13 );
    if( begin->location() == watch ) {
      std::cout << "found" << std::endl;
    }

    clear_ccdf();
    int status = cdf_estimator( *begin,
                                neighborhoods.begin(), neighborhoods.end(),
                                *ccdf_ );

    sampler( *begin, *ccdf_ );
    code_into_indicators( *begin );
  }

  return 0;
}

void pdf::generate_sample(int size, long double theta) {
	static unsigned int calls = 0; // use as seed for generator
	calls++; // ensures unique seed in every call
	sample.clear();
	boost::uniform_real<> unif(theta, theta * theta);
	boost::random::mt19937 gen(calls);
	boost::variate_generator< boost::random::mt19937&, boost::uniform_real<> > sampler(gen, unif);
	for(int i=0;i<size;i++) {
		sample.push_back(sampler());

	}

}

double ProfileNormal::getSample() {
  // Create a Mersenne twister random number generator
  // that is seeded once with #seconds since 1970
  static boost::mt19937 rng(static_cast<unsigned> (std::time(0)));
 
  // select Normal probability distribution
  boost::normal_distribution<double> dist(boost::math::mean(profile), sqrt(boost::math::variance(profile)));
 
  // bind random number generator to distribution, forming a function
  boost::variate_generator<boost::mt19937&, boost::normal_distribution<double> >  sampler(rng, dist);
  
  // sample from the distribution
  return sampler();
}

void MotionSample::recalculate(Pose end_pose, std::array<double, 4> params) {
    // pfcpp::VelocityMotionModelSampler sampler(std::move(params));
    pfcpp::OdometryMotionModelSampler sampler(params);

    m_moves.clear();
    m_moves.push_back(std::make_shared<Pose>(*m_startPose));
    m_moves.push_back(std::make_shared<Pose>(end_pose));

    m_samples.clear();
    for (int i = 0; i < 400; ++i) {
        auto p = sampler(*m_startPose, end_pose);
        m_samples.push_back(std::make_shared<Pose>(p));
    }
}

/**
 * Computes the curvature of the image and the force of the ACWE
 * @param {vector<cl::Event>} The events that this kernel depends on 
 * @param {bool} Indicates if we need to use all the bands of the image (by computing the avg)
 */
cl::Event ActiveContours::compCurvAndF(vector<cl::Event> vecEvPrev, bool useAllBands) {

    cl::Event evCurvF;
    if (WRITE) {
        cout << endl << " ----------------- Computing Curvature and F ---------" << endl;
    }
    try {
        cl::Context* context = clMan.getContext();
        cl::CommandQueue* queue = clMan.getQueue();
        cl::Program* program = clMan.getProgram();

        cl::Sampler sampler(*context, CL_FALSE,
                CL_ADDRESS_REPEAT, CL_FILTER_NEAREST, &err);

        cl::Kernel kernelCurvAndF(*program, (char*) "CurvatureAndF");
        kernelCurvAndF.setArg(0, img_phi);
        kernelCurvAndF.setArg(1, img_in);
        kernelCurvAndF.setArg(2, img_curv_F);
        kernelCurvAndF.setArg(3, sampler);
        kernelCurvAndF.setArg(4, buf_avg_in_out); 
        kernelCurvAndF.setArg(5, (int)useAllBands); 

        queue->enqueueNDRangeKernel(
                kernelCurvAndF,
                cl::NullRange,
                cl::NDRange((size_t) width, (size_t) height),
                cl::NDRange((size_t) grp_size_x, (size_t) grp_size_y),
                &vecEvPrev,
                &evCurvF);

    } catch (cl::Error ex) {
        clMan.printError(ex);
    }
    return evCurvF;
}

TEST(McmcNutsBaseNuts, transition) {

  rng_t base_rng(0);

  int model_size = 1;
  double init_momentum = 1.5;

  stan::mcmc::ps_point z_init(model_size);
  z_init.q(0) = 0;
  z_init.p(0) = init_momentum;

  stan::mcmc::mock_model model(model_size);
  stan::mcmc::mock_nuts sampler(model, base_rng);

  sampler.set_nominal_stepsize(1);
  sampler.set_stepsize_jitter(0);
  sampler.sample_stepsize();
  sampler.z() = z_init;

  std::stringstream output_stream;
  stan::interface_callbacks::writer::stream_writer writer(output_stream);
  std::stringstream error_stream;
  stan::interface_callbacks::writer::stream_writer error_writer(error_stream);

  stan::mcmc::sample init_sample(z_init.q, 0, 0);

  stan::mcmc::sample s = sampler.transition(init_sample, writer, error_writer);

  EXPECT_EQ(31.5, s.cont_params()(0));
  EXPECT_EQ(0, s.log_prob());
  EXPECT_EQ(1, s.accept_stat());
  EXPECT_EQ("", output_stream.str());
  EXPECT_EQ("", error_stream.str());
}

void mFakeMonkey::spike(){
	boost::mutex::scoped_lock lock(monkey_lock);
	boost::exponential_distribution<double> dist = boost::exponential_distribution<double>(spiking_rate->getValue().getFloat() / 1000000); 	
	variate_generator<boost::mt19937&, boost::exponential_distribution<double> > sampler = variate_generator<boost::mt19937&, boost::exponential_distribution<double> >(rng,dist);
		
	if(spike_node != NULL){
		spike_node->cancel();
	}
	
	if(spike_var != 0) {
		*spike_var = 1;
	}
	
	float delay = sampler();
	spike_node = scheduler->scheduleUS(FILELINE,
									   delay,
									   0,
									   1, 
									   &fake_monkey_spike, 
									   (void *)this, 
									   M_DEFAULT_IODEVICE_PRIORITY,
									   M_DEFAULT_IODEVICE_WARN_SLOP_US,
									   M_DEFAULT_IODEVICE_FAIL_SLOP_US,
									   M_MISSED_EXECUTION_CATCH_UP);
}

sk_sp<GrTextureProxy> GrTextureProducer::refTextureProxy(GrMipMapped willNeedMips,
                                                         SkColorSpace* dstColorSpace,
                                                         sk_sp<SkColorSpace>* proxyColorSpace) {
    GrSamplerState::Filter filter =
            GrMipMapped::kNo == willNeedMips ? GrSamplerState::Filter::kNearest
                                             : GrSamplerState::Filter::kMipMap;
    GrSamplerState sampler(GrSamplerState::WrapMode::kClamp, filter);

    int mipCount = SkMipMap::ComputeLevelCount(this->width(), this->height());
    bool willBeMipped = GrSamplerState::Filter::kMipMap == sampler.filter() && mipCount &&
                        fContext->contextPriv().caps()->mipMapSupport();

    auto result =
            this->onRefTextureProxyForParams(sampler, dstColorSpace, proxyColorSpace,
                                             willBeMipped, nullptr);

    // Check to make sure that if we say the texture willBeMipped that the returned texture has mip
    // maps, unless the config is not copyable.
    SkASSERT(!result || !willBeMipped || result->mipMapped() == GrMipMapped::kYes ||
             !fContext->contextPriv().caps()->isConfigCopyable(result->config()));

    // Check that no scaling occured and we returned a proxy of the same size as the producer.
    SkASSERT(!result || (result->width() == this->width() && result->height() == this->height()));
    return result;
}

void guess_EBV_profile(TMCMCOptions &options, TLOSMCMCParams &params, unsigned int N_regions) {
	TNullLogger logger;
	
	unsigned int N_steps = options.steps / 8;
	if(N_steps < 40) { N_steps = 40; }
	unsigned int N_samplers = options.samplers;
	unsigned int N_threads = options.N_threads;
	unsigned int ndim = N_regions + 1;
	
	TAffineSampler<TLOSMCMCParams, TNullLogger>::pdf_t f_pdf = &lnp_los_extinction;
	TAffineSampler<TLOSMCMCParams, TNullLogger>::rand_state_t f_rand_state = &gen_rand_los_extinction;
	
	std::cout << "Generating Guess ..." << std::endl;
	
	TParallelAffineSampler<TLOSMCMCParams, TNullLogger> sampler(f_pdf, f_rand_state, ndim, N_samplers*ndim, params, logger, N_threads);
	sampler.set_scale(1.05);
	sampler.set_replacement_bandwidth(0.75);
	
	sampler.step(int(N_steps*30./100.), true, 0., 0.5, 0.);
	sampler.step(int(N_steps*20./100), true, 0., 1., 0., true);
	sampler.step(int(N_steps*30./100.), true, 0., 0.5, 0.);
	sampler.step(int(N_steps*20./100), true, 0., 1., 0., true);
	
	sampler.print_stats();
	std::cout << std::endl << std::endl;
	
	sampler.get_chain().get_best(params.EBV_prof_guess);
	for(size_t i=0; i<ndim; i++) {
		//params.EBV_prof_guess[i] = log(params.EBV_prof_guess[i]);
		std::cout << "\t" << params.EBV_prof_guess[i] << std::endl;
	}
	std::cout << std::endl;
}

    TexOGL * TexOGL::CreateTexWidthImage( Image * _pImage,bool _mipmap )
    {
        GLint error;
        TexOGL * pTex = new TexOGL;
        pTex->m_desc.ePixelFormat = PIXEL_FORMAT_RGBA8888;
        pTex->m_desc.eTexClass = TEX_CLASS_STATIC_RAW;
		pTex->m_desc.size.width = _pImage->nWidth;
		pTex->m_desc.size.height = _pImage->nHeight;
        glGenTextures(1,&pTex->m_texture);
        error = glGetError();
        glBindTexture(GL_TEXTURE_2D,pTex->m_texture);        
        glTexImage2D(GL_TEXTURE_2D, 0,GL_RGBA, _pImage->nWidth, _pImage->nHeight, 0, GL_RGBA, GL_UNSIGNED_BYTE, _pImage->pData);
        if(_mipmap)
        {
            glGenerateMipmap(GL_TEXTURE_2D);            
        }
        // 新版api的创建方式
        //glTexStorage2D(GL_TEXTURE_2D,4,GL_RGBA8,64,64);
        //glTexSubImage2D(GL_TEXTURE_2D,0,0,0,64,64,GL_RGBA,GL_UNSIGNED_BYTE,bitmap);
        error = glGetError();
        SamplerOGL sampler(_mipmap);
        pTex->ApplySamplerState(&sampler);
        assert(GL_NO_ERROR == error);
        return pTex;
    }

TEST(McmcStaticBaseStaticHMC, set_nominal_stepsize_and_L) {

  rng_t base_rng(0);

  std::vector<double> q(5, 1.0);
  std::vector<int> r(2, 2);

  stan::mcmc::mock_model model(q.size());

  stan::mcmc::mock_static_hmc sampler(model, base_rng);

  double old_epsilon = 1.0;
  int old_L = 10;

  sampler.set_nominal_stepsize_and_L(old_epsilon, old_L);
  EXPECT_EQ(old_epsilon, sampler.get_nominal_stepsize());
  EXPECT_EQ(old_L, sampler.get_L());
  EXPECT_EQ(true, sampler.get_T() > 0);

  sampler.set_nominal_stepsize_and_L(-0.1, 5);
  EXPECT_EQ(old_epsilon, sampler.get_nominal_stepsize());
  EXPECT_EQ(old_L, sampler.get_L());

  sampler.set_nominal_stepsize_and_T(5.0, -1);
  EXPECT_EQ(old_epsilon, sampler.get_nominal_stepsize());
  EXPECT_EQ(old_L, sampler.get_L());
}

	color getPhong(
		const normal& i_N, const vector& i_V, const float cosinePower, 
		const eiBool i_keyLightsOnly, const eiBool unshadowed)
	{
		color C = 0.0f;
		vector R = reflect( normalize(i_V), normalize(i_N) );
		LightSampler	sampler(this, P, i_N, PI/2.0f );


			float isKeyLight = 1;
			//if( i_keyLightsOnly != 0 )
			//{
			//	lightsource( "iskeylight", isKeyLight );
			//}
			if( isKeyLight != 0 )
			{
				const float nonspecular = 0.0f;
				//lightsource( "__nonspecular", nonspecular );
				if( nonspecular < 1 )
				{
					//SAMPLE_LIGHT_2(color, C, 0.0f,
					//	C += Cl()*pow(max<float>(0.0f,R%Ln),cosinePower)*(1.0f-nonspecular);
					//);
					while (sampler.sample())
					{
						vector Ln = normalize(L);
						C += Cl*pow(max<float>(0.0f,R%Ln),cosinePower)*(1.0f-nonspecular);
					}
				}
			}

		return C;
	}

int main(){
	int i_dim = 2;
	double max[2] = {1.0,1.0};
	double min[2] = {-1.0,-1.0};

	//We create an MCSampler class 
	MCSample sampler(i_dim,min,max,&pdf);

	for(int i = 0; i < 10; i++){

		double * sample = sampler.genSample();
		double d_radius2 = 0.0;

		for(int j = 0; j < i_dim; j++){
			d_radius2 += sample[j]*sample[j];
		}

		std::cout<<"d_radius: "<<std::sqrt(d_radius2)<<", "<<sample[0]<<", "<<sample[1]<<std::endl;

		delete [] sample;

	}

	return 0;
}

TEST(McmcBaseStaticHMC, set_nominal_stepsize) {
  
  rng_t base_rng(0);
  
  std::vector<double> q(5, 1.0);
  std::vector<int> r(2, 2);
  
  stan::mcmc::mock_model model(q.size());

  std::stringstream output, error;
  
  stan::mcmc::mock_static_hmc sampler(model, base_rng, &output, &error);
  
  double old_epsilon = 1.0;
  
  sampler.set_nominal_stepsize(old_epsilon);
  EXPECT_EQ(old_epsilon, sampler.get_nominal_stepsize());
  EXPECT_EQ(true, sampler.get_L() > 0);
  
  sampler.set_nominal_stepsize(-0.1);
  EXPECT_EQ(old_epsilon, sampler.get_nominal_stepsize());
  
  EXPECT_EQ("", output.str());
  EXPECT_EQ("", error.str());
}

int main(int argc, char *argv[])
{
    if (argc != 3) {
        printf("Usage: %s <filename> <duration>\n", argv[0]);
        exit(1);
    }

    try {
        CDFSampler sampler(argv[1], atof(argv[2]));
        int i = 0;
        while (1) {
            double sample = sampler.sample();
            printf("Sample %d: %lf\n", ++i, sample);
            if (sample < 0.0) {
                /* for our purposes, the CDF has no negative x values. */
                throw CDFErr("negative sample!  BUG IN CDF_SAMPLER CODE.");
            }
            sleep(1);
        }
    } catch (CDFErr &e) {
        printf("Error: %s\n", e.str.c_str());
        exit(1);
    }
    
    return 0;
}

void HmmWithVonMisesObservations::doGenerateObservationsSymbol(const unsigned int state, const size_t n, dmatrix_type &observations) const
{
	if (!targetDist_) targetDist_.reset(new VonMisesTargetDistribution());
	targetDist_->setParameters(mus_[state], kappas_[state]);

#if 0
	// when using univariate normal proposal distribution.
	{
		// FIXME [modify] >> these parameters are incorrect.
		const double sigma = 1.55;
		const double k = 1.472;
		proposalDist_->setParameters(mus_[state], sigma, k);
	}
#else
	// when using univariate uniform proposal distribution.
	{
		const double lower = 0.0;
		const double upper = MathConstant::_2_PI;
		const UnivariateUniformProposalDistribution::vector_type mean_dir(1, mus_[state]);
		const double k = targetDist_->evaluate(mean_dir) * (upper - lower) * 1.05;
		proposalDist_->setParameters(lower, upper, k);
	}
#endif

	swl::RejectionSampling sampler(*targetDist_, *proposalDist_);

	swl::RejectionSampling::vector_type x(D_, 0.0);
	const std::size_t maxIteration = 1000;

	// the range of each observation, [0, 2 * pi)
	const bool retval = sampler.sample(x, maxIteration);
	assert(retval);
	observations(n, 0) = x[0];
}