double RandNorm(){
	random_device rnd;
	normal_distribution<double> norm(0.0, 1.0);
	return norm(rnd);
}

double MakeTree2(S s, Node *node, queue<Node *> q1, double theta_true){
	if(node->depth == 0){
		node->theta = RandNorm();
		node->item_no = (int)(rand() * (I + 1.0) / (1.0 + RAND_MAX));
	}else{
		node->theta = Theta_CAT(s, node);
		node->item_no = SelectItem(s, node, node->theta);
		cout << node->hist_item_no[node->depth-1] << "->";
	}

	presented[node->item_no] = 1;
	cout << node->item_no << "(" << node->theta << ")" << " ";

	Node *p_left;
	p_left = (Node *)malloc(sizeof(*node->incorrect));
	p_left = node->incorrect;
	if(node->depth > 0){
		for(int c = 0; c < node->depth; c++){
			p_left->hist_item_no.push_back(node->hist_item_no[c]);
			p_left->hist_response.push_back(node->hist_response[c]);
		}
	}
	p_left->hist_item_no.push_back(node->item_no);
	p_left->hist_response.push_back(0);
	p_left->depth = node->depth + 1;
	q1.push(p_left);

	Node *p_right;
	p_right = (Node *)malloc(sizeof(*node->correct));
	p_right = node->correct;
	if(node->depth > 0){
		for(int c = 0; c < node->depth; c++){
			p_right->hist_item_no.push_back(node->hist_item_no[c]);
			p_right->hist_response.push_back(node->hist_response[c]);
		}
	}
	p_right->hist_item_no.push_back(node->item_no);
	p_right->hist_response.push_back(1);
	p_right->depth = node->depth + 1;
	q1.push(p_right);

	int nodes = 0, rest = 0;
	for(int i = 0; i < I; i++){
		if(presented[i] == 1) nodes++;
	}
	if(nodes == pow(2, node->depth + 1) - 1){
		cout << endl;
		for(int i = 0; i < I; i++){
			if(presented[i] == 0) rest++;
		}
		if(pow(2, node->depth + 1) > rest){
			free(p_left);
			free(p_right);
			return node->depth;
		}
	}

	Node *next;
	next = (Node *)malloc(sizeof(*q1.front()));
	next = q1.front();
	q1.pop();
	double depth_of_tree = MakeTree2(s, next, q1, theta_true);
	free(next);

	if(node->depth == 0){
		double theta_est;
		Node *current;
		current = (Node *)malloc(sizeof(*node));
		current = node;
		for(int d = 0; d <= (int)depth_of_tree; d++){
			if((double) rand() / RAND_MAX < P_X_h(s.a[current->item_no], s.b[current->item_no], theta_true)){
				current = node->correct;
			}else{
				current = node->incorrect;
			}
			theta_est = node->theta;
		}
		free(current);
		return theta_est;
	}

	return depth_of_tree;
}

void BreakTree(Node *node){
	if(node->incorrect != NULL){
		BreakTree(node->incorrect);
	}
	if(node->correct != NULL){
		BreakTree(node->correct);
	}
	free(node);
}

double CAT_sim2(S s){
	double theta_true, theta_est, error;
	theta_true = RandNorm();

	queue<Node *> q1;
	Node root, *p_root = &root;

	p_root->depth = 0;
	p_root->theta = 0.0;
	theta_est = MakeTree2(s, p_root, q1, theta_true);
	BreakTree(p_root);

	error = theta_true - theta_est;
	cout << "theta_true = " << theta_true << ", theta_estimated = " << theta_est << ", error = " << error << endl;
	return error;
}

double CAT_sim_exam(S s){
	int count = 10;
	double error[count] = {};
	double error_sum = 0;

	for(int c = 0; c < count; c++){
		error[c] = CAT_sim2(s);
		error_sum += error[c];
	}
	cout << "MSE: " << error_sum / (double)count << endl;
}