// -------------------------------------------------------------------------------------------------
// Draw IDs with a uniform probability over [1, M]
// -------------------------------------------------------------------------------------------------
	clear all

// -------------------------------------------------------------------------------------------------
// Programs
// -------------------------------------------------------------------------------------------------

// Uniform draw
capture program drop smpl
program define smpl
	syntax newvarname, size(integer)
	gen long `varlist' = 1 + int(uniform() * `size')
	compress `varlist'
end


// -------------------------------------------------------------------------------------------------
// Hard scenario
// -------------------------------------------------------------------------------------------------
	clear
	set trace off
	set more off
	timer clear

	local M1 125000
	local q = 1 // 0.5
	local M2 = `M1' * `q' // We want M1 and M2 to grow at the same rate
	local d 4
	local N = `M2' * `d' // If we require -d- to stay small, then the dataset should be proportional to the Ms

	set obs `N'
	set seed 1234
	gen double x1 = rnormal()
	gen double x2 = uniform()

	smpl id1, size(`M1')
	smpl id2, size(`M2')

	gen double y = 10 + x1 + x2 + sin(id1) + log(10+id2) + 100 * rnormal()
	assert !missing(y)

	su id1 id2
	sort id1 id2

	saveold dataset_hard, version(12) replace
	outsheet using dataset_hard.txt, nolabel noquote replace // tab-separated

// -------------------------------------------------------------------------------------------------

	* Benchmark

	// Load .mlib
	//qui reghdfe y, noabsorb dof(none) nosample noprune keepsing v(-1)
	//

	//reghdfe y x1 x2, a(id1 id2) accel(no) transform(kac) keepsing noprune // cnverge in 2
	//reghdfe y x1 x2, a(id1 id2) accel(sd) transform(kac) keepsing noprune // converge in 2
	//reghdfe y x1 x2, a(id1 id2) accel(cg) transform(sym) keepsing noprune // converge in 2 (more if we set d!=1 in the code)
	set rmsg on
	reghdfe y x1 x2, a(id1 id2) accel(cg) transform(sym) noprune keepsing dof(none) nosample 
	set rmsg off

exit