class KoverDataset(object):
	def __init__(self, file):
		self.path = file
		self.dataset_open = partial(_hdf5_open_no_chunk_cache, file)

	@property
	def classification_type(self):
		dataset = self.dataset_open()
		# Backwards compatibilty with pre-2.0.0 Kover datasets
		try:
			classification_type = dataset.attrs["classification_type"]
		except:
			classification_type = "binary"
		return classification_type

	@property
	def compression(self):
		dataset = self.dataset_open()
		return dataset.attrs["compression"]

	@property
	def kmer_filter(self):
		dataset = self.dataset_open()
		return dataset.attrs["filter"]

	@property
	def genome_count(self):
		dataset = self.dataset_open()
		return dataset["genome_identifiers"].shape[0]

	@property
	def genome_identifiers(self):
		dataset = self.dataset_open()
		return dataset["genome_identifiers"]

	@property
	def genome_source(self):
		dataset = self.dataset_open()
		return dataset.attrs["genomic_data"]

	@property
	def genome_source_type(self):
		dataset = self.dataset_open()
		return dataset.attrs["genome_source_type"]

	@property
	def kmer_by_matrix_column(self):
		dataset = self.dataset_open()
		return dataset["kmer_by_matrix_column"]

	@property
	def kmer_count(self):
		dataset = self.dataset_open()
		return dataset["kmer_sequences"].shape[0]

	@property
	def kmer_length(self):
		dataset = self.dataset_open()
		return len(dataset["kmer_sequences"][0])

	@property
	def kmer_matrix(self):
		dataset = self.dataset_open()
		return dataset["kmer_matrix"]

	@property
	def kmer_sequences(self):
		dataset = self.dataset_open()
		return dataset["kmer_sequences"]

	@property
	def phenotype(self):
		dataset = self.dataset_open()
		# Backwards compatibilty with pre-2.0.0 Kover datasets
		try:
			description = dataset.attrs["phenotype_description"]
		except:
			description = dataset.attrs["phenotype_name"]

		try:
			tags = dataset["phenotype_tags"]
		except:
			tags = np.array(['0', '1'])

		return KoverDatasetPhenotype(description=description,
									 tags=tags,
									 metadata=dataset["phenotype"],
									 metadata_source=dataset.attrs["phenotype_metadata_source"])

	@property
	def splits(self):
		dataset = self.dataset_open()
		if "splits" in dataset:
			return [self.get_split(split_name) for split_name in dataset["splits"]]
		else:
			return []

	@property
	def uuid(self):
		dataset = self.dataset_open()
		return dataset.attrs["uuid"]

	def get_split(self, name):
		dataset = self.dataset_open()
		split = dataset["splits"][name]
		return KoverDatasetSplit(name,
								 split.attrs["train_proportion"],
								 # Backwards compatibility with datasets without test_proportion
								 split.attrs["test_proportion"] if "test_proportion" in split.attrs else
																1.0 - split.attrs["train_proportion"],
								 split["train_genome_idx"],
								 split["test_genome_idx"],
								 split["unique_risks"],
								 split["unique_risk_by_kmer"],
								 split["unique_risk_by_anti_kmer"],
								 [] if not "folds" in split
								 else [KoverDatasetFold(fold_name,
														fold["train_genome_idx"],
														fold["test_genome_idx"],
														fold["unique_risks"],
														fold["unique_risk_by_kmer"],
														fold["unique_risk_by_anti_kmer"]) for fold_name, fold in split["folds"].iteritems()],
								 split.attrs["random_seed"]) class DecisionTreeClassifier(object):
	def __init__(self, criterion, max_depth, min_samples_split, class_importance):
		# Validate the node splitting criterion
		supported_criteria = ["gini", "cross-entropy"]
		if criterion not in supported_criteria:
			raise ValueError("The supporting splitting criteria are: %s." % str(supported_criteria))
		self.criterion = criterion

		# Validate the maximum depth of the tree
		if max_depth < 1:
			raise ValueError("The maximum tree depth must be greater than 1.")
		self.max_depth = max_depth

		# Validate the minimum number of examples required to split a node
		if min_samples_split < 2.0:
			raise ValueError( "The minimum number of examples used to split a node must be 2 or greater.")
		self.min_samples_split = int(min_samples_split)

		# A dictionnary that gives the importance of each class for the loss minimization
		self.class_importance = class_importance

	def fit(self, rules, rule_classifications, example_idx, rule_blacklist=None,
			tiebreaker=None, level_callback=None, split_callback=None):

		"""
		Fits the decision tree classifier
		"""

		if level_callback is None:
			level_callback = lambda x: None

		if split_callback is None:
			split_callback = lambda x, y: None

		if tiebreaker is None:
			tiebreaker = lambda x: x

		# Store important information about the training set
		n_total_class_examples = {c: float(len(idx)) for c, idx in example_idx.iteritems()}

		# Compute the class priors based on the importance of making errors on each class.
		# See Section "4.4 Priors and Variable Misclassification Costs" in Breiman et al. 1984.
		priors = {c: 1.0 * n_examples / sum(n_total_class_examples.values()) for c, n_examples in n_total_class_examples.iteritems()}
		denum = sum([self.class_importance[c] * priors[c] for c in priors.iterkeys()])
		altered_priors = {c: 1.0 * self.class_importance[c] * priors[c] / denum for c in priors.iterkeys()}

		# Criteria for node impurity and splitting

		#############################################
		#				GINI IMPURITY				#
		#############################################

		def _gini_impurity(n_examples_by_class, multiply_by_node_proba=False):
			"""
			Computes the gini impurity for a node. Also supports multiplying by
			the node probability if this is required.

			Notes:
				* equation numbers refer to Breiman et al. (1984)
				* j = a class
				* t = current node
				* This code works if n_examples_by_class[c] is an array giving
				  the value for multiple candidate splits, or an integer giving
				  the value for a single split.

			"""
			p_j_t = {c: 1.0 * altered_priors[c] * n_examples_by_class[c] / n_total_class_examples[c] \
						for c in n_examples_by_class.keys()}  # Equation (2.2)

			p_t = sum(p_j_t.values())  # Equation (2.3)

			with np.errstate(divide='ignore', invalid='ignore'):
				p_j_given_t = {c: np.divide(p_j_t[c], p_t) for c in p_j_t.keys()}

			gini_diversity_index = sum([p_j_given_t[i] * p_j_given_t[j] \
									for i in p_j_given_t.keys() for j in p_j_given_t.keys() if i != j])

			return gini_diversity_index * (p_t if multiply_by_node_proba else 1.0)

		def _gini_rule_score(example_idx):
			"""
			example_idx: a dictionnary where the keys are classes and the values
						 give the indices of the examples with the corresponding
						 class that are in the node to split.

			"""
			logging.debug("Scoring rules with the gini impurity strategy")
			# A k-mer rule splits the examples in two groups: those that don't have the k-mer in their
			# genome (left child) and those that have the k-mer in their genome (right child).
			# XXX: We keep only the first half of the rule list and classifications, since sum_rows returns the counts
			#      for presence and absence rules, which we don't need here.
			# TODO: We could have a parameter called return_absence=True, which avoid using twice the RAM.
			last_presence_rule_idx = int(1.0 * len(rules) / 2)  # Just because sum_rows returns pres/abs rules
			# For each class, we compute a vector that gives the number of examples of this class
			# that contain each k-mer. This is the number of examples that will go in the left leaf
			# if a split is made on a given k-mer rule.
			left_n_examples_by_class = \
				{c: np.asarray(rule_classifications.sum_rows(example_idx[c])[: last_presence_rule_idx], dtype=np.float) \
					for c in example_idx.keys()}

			# Similarly, we compute the number of examples that would be sent to the right leaf (don't contain the k-mer)
			right_n_examples_by_class = {c: np.asarray(len(example_idx[c]) - left_n_examples_by_class[c], dtype=np.float) \
											for c in left_n_examples_by_class.keys()}

			# XXX: To maximize the decrease in impurity, we can simply minimize the sum of the gini impurity of the
			#      resulting leaves, so we only compute those quantities.

			# We will compute the gini impurities in blocks
			n_kmers = left_n_examples_by_class[left_n_examples_by_class.keys()[0]].shape[0]
			BLOCK_SIZE = 100000
			gini = np.zeros(n_kmers)
			n_blocks = int(ceil(1.0 * n_kmers / BLOCK_SIZE))

			# Left child:
			for i in range(n_blocks):
				block_examples_by_class = {c: ex[i * BLOCK_SIZE : (i + 1) * BLOCK_SIZE] for c, ex in left_n_examples_by_class.iteritems()}
				gini[i * BLOCK_SIZE : (i + 1) * BLOCK_SIZE] = _gini_impurity(block_examples_by_class, multiply_by_node_proba=True)

			# Right child:
			for i in range(n_blocks):
				block_examples_by_class = {c : ex[i * BLOCK_SIZE : (i + 1) * BLOCK_SIZE] for c, ex in right_n_examples_by_class.iteritems()}
				gini[i * BLOCK_SIZE : (i + 1) * BLOCK_SIZE] += _gini_impurity(block_examples_by_class, multiply_by_node_proba=True)

			# Don't consider rules that lead to empty nodes
			# XXX: we minimize the sum of impurities, so setting to inf will prevent the rules from being selected
			gini[sum(left_n_examples_by_class.values()) == 0] = np.infty
			gini[sum(right_n_examples_by_class.values())  == 0] = np.infty

			return gini

		#############################################
		#				CROSS-ENTROPY				#
		#############################################

		def _cross_entropy(n_class_examples, multiply_by_node_proba=False):
			p_class_node = {c: 1.0*altered_priors[c]*n_class_examples[c]/n_total_class_examples[c] \
											for c in n_class_examples.keys()}

			node_resubstitution_estimate = sum(p_class_node.values())
			with np.errstate(divide='ignore', invalid='ignore'):
				p_class_given_node = {c: np.divide(p_class_node[c], node_resubstitution_estimate) for c in p_class_node.keys()}
				diversity_index = -1.0 * sum([np.nan_to_num(p_class_given_node[c]*np.log(p_class_given_node[c])) \
																						for c in p_class_given_node.keys()])
			return diversity_index * (node_resubstitution_estimate if multiply_by_node_proba else 1.0)

		def _cross_entropy_rule_score(example_idx):
			"""
			example_idx: a dictionnary where the keys are classes and the values
						 give the number of examples of the class in the node to
						 split

			"""
			logging.debug("Scoring rules with the gini impurity strategy")
			# A k-mer rule splits the examples in two groups: those that don't have the k-mer in their
			# genome (left child) and those that have the k-mer in their genome (right child).
			# XXX: We keep only the first half of the rule list and classifications, since sum_rows returns the counts
			# XXX: for presence and absence rules, which we don't need here.
			# TODO: We could have a parameter called return_absence=True, which avoid using twice the RAM.

			last_presence_rule_idx = int(1.0 * len(rules) / 2)  # Just because sum_rows returns pres/abs rules

			left_count = {c:np.asarray(rule_classifications.sum_rows(example_idx[c])[: last_presence_rule_idx],dtype=np.float) \
						for c in example_idx.keys() if example_idx[c].size}
			right_count = {c:np.asarray(example_idx[c].shape[0] - left_count[c], dtype=np.float) for c in left_count.keys()}

			# Left child:
			cross_entropy = _cross_entropy(left_count, multiply_by_node_proba=True)
			# Right child:
			cross_entropy += _cross_entropy(right_count, multiply_by_node_proba=True)

			# Don't consider rules that lead to empty nodes
			cross_entropy[sum(left_count.values()) == 0] = np.infty
			cross_entropy[sum(right_count.values())  == 0] = np.infty

			return cross_entropy

		if self.criterion == "gini":
			get_criterion = _gini_impurity
			score_rules = _gini_rule_score
			choice_func = min
		elif self.criterion == "cross-entropy":
			get_criterion = _cross_entropy
			score_rules = _cross_entropy_rule_score
			choice_func = min
		node_type = ProbabilisticTreeNode

		def _find_best_split(node):
			"""
			Selects the best split according to the splitting criterion

			"""
			example_idx = node.class_examples_idx

			# Score all the rules according to the criterion
			rules_criterion = score_rules(example_idx=example_idx)

			# Remove rules that are blacklisted
			logging.debug("Removing {0:d} blacklisted rules".format(len(rule_blacklist)))
			rules_criterion[rule_blacklist] = (np.infty if choice_func == min else -np.infty)

			# Check if we were able to find a split
			if (choice_func == min and min(rules_criterion) == np.infty) or (choice_func == max and max(rules_criterion) == -np.infty):
				return None, None, None, None

			# Tiebreaker to select a single rule in case of ties
			candidate_rules_idx = np.where(rules_criterion == choice_func(rules_criterion))[0]
			logging.debug("There are %d equivalent rules before the tiebreaker." % len(candidate_rules_idx))
			best_rules_idx = tiebreaker(candidate_rules_idx)
			logging.debug("The are %d equivalent rules after the tiebreaker." % len(best_rules_idx))
			selected_rule_idx = best_rules_idx[0]

			# Dispatch the examples to the children nodes based on the splitting rule's predictions
			rule_preds = rule_classifications.get_columns(selected_rule_idx)

			left_child_example_idx_by_class = {c: example_idx[c][rule_preds[example_idx[c]] == 1] for c in example_idx.keys()}
			right_child_example_idx_by_class = {c: example_idx[c][rule_preds[example_idx[c]] == 0] for c in example_idx.keys()}

			return selected_rule_idx, best_rules_idx, left_child_example_idx_by_class, right_child_example_idx_by_class

		logging.debug("Training start.")

		root = node_type(class_examples_idx=example_idx,
						 depth=0,
						 criterion_value=get_criterion(n_total_class_examples),
						 class_priors=altered_priors,
						 total_n_examples_by_class=n_total_class_examples)

		nodes_to_split = deque([root])
		runtime_infos = {}
		current_depth = -1
		min_samples_split = self.min_samples_split
		if min_samples_split < 2:
			min_samples_split = 2

		while len(nodes_to_split) > 0:
			node = nodes_to_split.popleft()

			# Check if we have reached a new depth
			if node.depth != current_depth:
				current_depth = node.depth

				runtime_infos["depth"] = current_depth

				logging.debug("The tree depth is %d" % current_depth)
				if current_depth > 0:
					# The level callback is called when all the nodes of a level have been created
					level_callback(runtime_infos)
				if current_depth == self.max_depth:
					logging.debug("The maximum tree depth has been reached. No more leaves will be split.")
					break  # We have reached the nodes of the last level, which must remain leaves

			# Check if the node to split is a pure leaf
			if 1.0 in node.class_proportions.values():
				logging.debug("The leaf is pure. It will not be split.")
				continue

			# Check if the HP constraints allows us to split this node
			if node.n_examples < min_samples_split:
				logging.debug("The leaf contains less examples (%d) than the minimum required to split (%d) a node. "
				"It will not be split." % (node.n_examples, min_samples_split))
				continue

			# Find the best rule to split the node
			selected_rule_idx, \
			equivalent_rule_idx, \
			left_child_example_idx_by_class, \
			right_child_example_idx_by_class = _find_best_split(node)

			# If we were incapable of splitting the node into two non-empty leafs
			if selected_rule_idx is None:
				logging.debug("Found no rule to split the node. The node will not be split.")
				continue

			# Perform the split
			node.rule = rules[selected_rule_idx]
			left_child_n_class_example = {c: len(idx) for c, idx in left_child_example_idx_by_class.items()}
			right_child_n_class_example = {c: len(idx) for c, idx in right_child_example_idx_by_class.items()}

			node.left_child = node_type(parent=node,
										class_examples_idx=left_child_example_idx_by_class,
										depth=node.depth + 1,
										criterion_value=get_criterion(left_child_n_class_example),
										class_priors=altered_priors,
										total_n_examples_by_class=n_total_class_examples)

			node.right_child = node_type(parent=node,
										class_examples_idx=right_child_example_idx_by_class,
										depth=node.depth + 1,
										criterion_value=get_criterion(right_child_n_class_example),
										class_priors=altered_priors,
										total_n_examples_by_class=n_total_class_examples)

			# Inject the rule importance in the rule object (unnormalized)
			node.rule.importance = \
		        node.breiman_info.p_t * node.criterion_value - \
				node.left_child.breiman_info.p_t * node.left_child.criterion_value - \
				node.right_child.breiman_info.p_t * node.right_child.criterion_value

			logging.debug("Split with rule %s." % node.rule)
			split_callback(node, equivalent_rule_idx)

			# Add the new child nodes to the splitting queue
			nodes_to_split.append(node.left_child)
			nodes_to_split.append(node.right_child)

			# Update the model in the runtime informations
			runtime_infos["model"] = root

		logging.debug("Done building the tree.")

		# Save the decision tree
		self.decision_tree = root

		logging.debug("Training finished.")

	def predict(self, X):
		if not self._is_fitted():
			raise RuntimeError("The classifier must be fitted before predicting.")
		return self.decision_tree.predict(X)

	def predict_proba(self, X):
		if not self._is_fitted():
			raise RuntimeError("The classifier must be fitted before predicting.")
		return self.decision_tree.predict_proba(X)

	def _is_fitted(self):
		return self.decision_tree is not None def _prune_tree(tree):
	"""
	Prunes a decision tree

	"""
	def _initial_pruning(root):
		# Get all nodes that have two leaves as children
		def _get_leaf_parents(root):
			def __get_leaf_parents(node):
				leaf_parents = []
				if not node.is_leaf:
					if node.left_child.is_leaf and node.right_child.is_leaf:
						leaf_parents.append(node)
					else:
						leaf_parents += __get_leaf_parents(node.left_child)
						leaf_parents += __get_leaf_parents(node.right_child)
				return leaf_parents
			return __get_leaf_parents(root)

		# Perform the initial pruning (Tmax -> T1)
		def __initial_pruning(parents):
			if len(parents) == 0:
				return
			node = parents.pop()
			if np.allclose(node.breiman_info.R_t, node.left_child.breiman_info.R_t + node.right_child.breiman_info.R_t):
				logging.debug("Converting node %s into a leaf" % node)
				del node.rule
				node.rule = None
				del node.left_child
				node.left_child = None
				del node.right_child
				node.right_child = None
				assert node.is_leaf  # Runtime validation
				if not node.is_root and node.parent.left_child.is_leaf and node.parent.right_child.is_leaf:
					logging.debug("Adding the new leaf's parent to the list of leaf parents")
					parents.append(node.parent)
			__initial_pruning(parents)

		# Start the initial pruning
		__initial_pruning(_get_leaf_parents(root))

	def _find_weakest_links(root):
		def __find_weakest_links(node):
			if node.is_leaf:
				return np.inf, [node]
			else:
				RTt = sum(l.breiman_info.R_t for l in node.leaves)
				current_gt = float(node.breiman_info.R_t - RTt) / (len(node.leaves) - 1)
				left_min_gt, left_weakest_links = __find_weakest_links(node.left_child)
				right_min_gt, right_weakest_links = __find_weakest_links(node.right_child)

				if np.allclose(current_gt, min(left_min_gt, right_min_gt)):
					if np.allclose(left_min_gt, right_min_gt):
						return current_gt, [node] + left_weakest_links + right_weakest_links
					else:
						return current_gt, [node] + (left_weakest_links if left_min_gt < right_min_gt else right_weakest_links)
				elif current_gt < min(left_min_gt, right_min_gt):
					return current_gt, [node]
				elif np.allclose(left_min_gt, right_min_gt):
					return left_min_gt, left_weakest_links + right_weakest_links
				elif left_min_gt > right_min_gt:
					return right_min_gt, right_weakest_links
				elif left_min_gt < right_min_gt:
					return left_min_gt, left_weakest_links
				else:
					raise Exception("Unhandled case detected!")

		return __find_weakest_links(root)

	def _sequential_prune(root):
		Tmax = deepcopy(root)
		logging.debug("Initial pruning (Tmax >> T1)")
		_initial_pruning(Tmax)
		T1 = Tmax
		del Tmax

		def __sequential_prune(root):
			root = deepcopy(root)

			# Find the weakest link in the tree
			logging.debug("Computing link strenghts for each node")
			min_gt, weakest_links = _find_weakest_links(root)

			# Prune the weakest link (and save alpha)
			logging.debug("Pruning occurs at alpha %.9f" % min_gt)
			# TODO: faire une fonction qui descend jusqu'aux feuilles et prune en remontant car perte de mémoire ici

			for n in weakest_links:
				del n.rule
				n.rule = None
				del n.left_child
				n.left_child = None
				del n.right_child
				n.right_child = None
				assert n.is_leaf  # Runtime validation

			# Repeat until only the root node remains
			return [(min_gt, root)] + (__sequential_prune(root) if not root.is_leaf else [])

		logging.debug("Pruning sequentially until only the root remains (T1 >> ... >> {root}")
		return [(0, T1)] + (__sequential_prune(T1) if not T1.is_leaf else [])

	logging.debug("Copying model")
	tree = deepcopy(tree)

	logging.debug("Generating the sequence of pruned trees")
	alphas, trees = zip(*_sequential_prune(tree))

	return alphas, trees class CARTModel(BaseModel):
    def __init__(self, class_tags=None):
        super(CARTModel, self).__init__()
        self.decision_tree = None
        self.class_tags = class_tags

    def predict(self, X):
        if self.decision_tree is None:
            raise RuntimeError("A decision tree must be fitted prior to calling predict.")
        predictions = self.decision_tree.predict(X)
        return np.asarray(predictions, dtype=np.uint8)

    def predict_proba(self, X):
        if self.decision_tree is None:
            raise RuntimeError("A decision tree must be fitted prior to calling predict.")
        self.decision_tree.predict_proba(X)

    @property
    def learner(self):
        return cart

    def _to_string(self, node=None,depth=0):
        if node is None:
            if self.decision_tree is None:
                print("No tree has been added to the model")
            node = self.decision_tree

        if self.class_tags is None:
            return str(self.decision_tree)

        tree_str = ""

        # Case : node is a leaf
        if node.is_leaf:
			tree_str += "\n" + ("    " * depth) + str(self.class_tags[node.class_prediction])

        # Case : node has two children
        else:
            # Print right branch
            tree_str += self._to_string(node=node.right_child, depth=depth + 1)

            # Print own value
            tree_str += "\n" + ("    " * depth + "   ") + str("/")
            tree_str += "\n" + ("    " * depth) + str(node.rule)
            tree_str += "\n" + ("    " * depth + "   ") + str("\\")

            # Print left_child branch
            tree_str += self._to_string(node=node.left_child, depth=depth + 1)

        return tree_str

    def __len__(self):
        if self.decision_tree is None:
            return 0
        return len(self.decision_tree)

    @property
    def depth(self):
        if self.decision_tree is None:
            return 0
        return self.decision_tree.tree_depth class LazyKmerRuleList(object):
    """
    By convention, the first half of the list contains presence rules and the second half contains the absence rules in
    the same order.
    """
    def __init__(self, kmer_sequences, kmer_by_rule):
        self.n_rules = kmer_by_rule.shape[0] * 2
        self.kmer_sequences = kmer_sequences
        self.kmer_by_rule = kmer_by_rule

    def __getitem__(self, idx):
        if idx >= self.n_rules:
            raise ValueError("Index %d is out of range for list of size %d" % (idx, self.n_rules))
        if idx >= len(self.kmer_sequences):
            type = "absence"
            kmer_idx = self.kmer_by_rule[idx % len(self.kmer_sequences)]
        else:
            type = "presence"
            kmer_idx = self.kmer_by_rule[idx]
        return KmerRule(idx % len(self.kmer_sequences), self.kmer_sequences[kmer_idx], type)

    def __len__(self):
        return self.n_rules class KmerRuleClassifications(BaseRuleClassifications):
    """
    Methods involving columns account for presence and absence rules
    """
    # TODO: Clean up. Get rid of the code to handle deleted rows. We don't need this.
    def __init__(self, dataset, n_rows, block_size=None):
        self.dataset = dataset
        self.dataset_initial_n_rows = n_rows
        self.dataset_n_rows = n_rows
        self.dataset_removed_rows = []
        self.dataset_removed_rows_mask = np.zeros(self.dataset_initial_n_rows, dtype=np.bool)
        self.block_size = (None, None)

        if block_size is None:
            if self.dataset.chunks is None:
                self.block_size = (1, self.dataset.shape[1])
            else:
                self.block_size = self.dataset.chunks
        else:
            if len(block_size) != 2 or not isinstance(block_size[0], int) or not isinstance(block_size[1], int):
                raise ValueError("The block size must be a tuple of 2 integers.")
            self.block_size = block_size

        # Get the size of the ints used to store the data
        if self.dataset.dtype == np.uint32:
            self.dataset_pack_size = 32
            self.inplace_popcount = inplace_popcount_32
        elif self.dataset.dtype == np.uint64:
            self.dataset_pack_size = 64
            self.inplace_popcount = inplace_popcount_64
        else:
            raise ValueError("Unsupported data type for packed attribute classifications array. The supported data" +
                             " types are np.uint32 and np.uint64.")

        super(BaseRuleClassifications, self).__init__()

    def get_columns(self, columns):
        """
        Columns can be an integer (or any object that implements __index__) or a sorted list/ndarray.
        """
        #TODO: Support slicing, make this more efficient than getting the columns individually.
        columns_is_int = False
        if hasattr(columns, "__index__"):  # All int types implement the __index__ method (PEP 357)
            columns = [columns.__index__()]
            columns_is_int = True
        elif isinstance(columns, np.ndarray):
            columns = columns.tolist()
        elif isinstance(columns, list):
            pass
        else:
            columns = list(columns)

        # Detect where an inversion is needed (columns corresponding to absence rules)
        columns, invert_result = zip(* (((column if column < self.dataset.shape[1] else column % self.dataset.shape[1]),
                                         (True if column > self.dataset.shape[1] else False)) for column in columns))
        columns = list(columns)
        invert_result = np.array(invert_result)

        # Don't return rows that have been deleted
        row_mask = np.ones(self.dataset.shape[0] * self.dataset_pack_size, dtype=np.bool)
        row_mask[self.dataset_initial_n_rows:] = False
        row_mask[self.dataset_removed_rows] = False

        # h5py requires that the column indices are sorted
        unique, inverse = np.unique(columns, return_inverse=True)
        result = _unpack_binary_bytes_from_ints(self.dataset[:, unique.tolist()])[row_mask]
        result = result[:, inverse]
        result[:, invert_result] = 1 - result[:, invert_result]

        if columns_is_int:
            return result.reshape(-1)
        else:
            return result

    def remove_rows(self, rows):
        # Find in which dataset the rows must be removed
        dataset_removed_rows = []
        # TODO: This is inefficient!
        for row_idx in rows:
            current_idx = -1
            n_active_elements_seen = 0
            while n_active_elements_seen <= row_idx:
                current_idx += 1
                if not self.dataset_removed_rows_mask[current_idx]:
                    n_active_elements_seen += 1
            dataset_removed_rows.append(current_idx)
        # Update the dataset removed row lists
        # Update the start and stop indexes
        # Adjust the shape
        # Adjust the number of rows in each dataset
        # Store the sorted relative removed row indexes by dataset
        if len(dataset_removed_rows) > 0:
            self.dataset_removed_rows = sorted(set(self.dataset_removed_rows + dataset_removed_rows))
            self.dataset_removed_rows_mask = np.zeros(self.dataset_initial_n_rows, dtype=np.bool)
            self.dataset_removed_rows_mask[self.dataset_removed_rows] = True
            self.dataset_n_rows = self.dataset_initial_n_rows - len(self.dataset_removed_rows)

    @property
    def shape(self):
        return self.dataset_n_rows, self.dataset.shape[1] * 2

    # TODO: allow summing over multiple lists of rows at a time (saves i/o operations)
    def sum_rows(self, rows):
        """
        Note: Assumes that the rows argument does not contain duplicate elements. Rows will not be considered more than once.
        """
        rows = np.asarray(rows)
        result_dtype = _minimum_uint_size(rows.shape[0])
        result = np.zeros(self.dataset.shape[1] * 2, dtype=result_dtype)

        # Builds a mask to turn off the bits of the rows we do not want to count in the sum.
        def build_row_mask(example_idx, n_examples, mask_n_bits):
            if mask_n_bits not in [8, 16, 32, 64, 128]:
                raise ValueError("Unsupported mask format. Use 8, 16, 32, 64 or 128 bits.")

            n_masks = int(ceil(float(n_examples) / mask_n_bits))
            masks = [0] * n_masks

            for idx in example_idx:
                example_mask = idx / mask_n_bits
                example_mask_idx = mask_n_bits - (idx - mask_n_bits * example_mask) - 1
                masks[example_mask] |= 1 << example_mask_idx

            return np.array(masks, dtype="u" + str(mask_n_bits / 8))

        # Find the rows that occur in each dataset and their relative index
        rows = np.sort(rows)
        dataset_relative_rows = []
        for row_idx in rows:
            # Find which row in the dataset corresponds to the requested row
            # TODO: This is inefficient! Could exploit the fact that rows is sorted to reuse previous iterations.
            current_idx = -1
            n_active_elements_seen = 0
            while n_active_elements_seen <= row_idx:
                current_idx += 1
                if not self.dataset_removed_rows_mask[current_idx]:
                    n_active_elements_seen += 1
            dataset_relative_rows.append(current_idx)

        # Create a row mask for each dataset
        row_mask = build_row_mask(dataset_relative_rows, self.dataset_initial_n_rows, self.dataset_pack_size)
        del dataset_relative_rows

        # For each dataset load the rows for which the mask is not 0. Support column slicing aswell
        n_col_blocks = int(ceil(1.0 * self.dataset.shape[1] / self.block_size[1]))
        rows_to_load = np.where(row_mask != 0)[0]
        n_row_blocks = int(ceil(1.0 * len(rows_to_load) / self.block_size[0]))

        for row_block in xrange(n_row_blocks):
            block_row_mask = row_mask[rows_to_load[row_block * self.block_size[0]:(row_block + 1) * self.block_size[0]]]

            for col_block in xrange(n_col_blocks):

                # Load the appropriate rows/columns based on the block sizes
                block = self.dataset[rows_to_load[row_block * self.block_size[0]:(row_block + 1) * self.block_size[0]],
                                     col_block * self.block_size[1]:(col_block + 1) * self.block_size[1]]

                # Popcount
                if len(block.shape) == 1:
                    block = block.reshape(1, -1)
                self.inplace_popcount(block, block_row_mask)

                # Increment the sum
                result[col_block * self.block_size[1]:min((col_block + 1) * self.block_size[1], self.dataset.shape[1])] += np.sum(block, axis=0)

        # Compute the sum for absence rules
        result[self.dataset.shape[1] : ] = len(rows) - result[: self.dataset.shape[1]]

        return result def _duplicate_last_element(l, length):
    """
    Duplicates the last element of a list until a given length is reached. (In-place)
    """
    l += [l[-1]] * (length - len(l))
    return l def _init_callback_functions(warning_callback = None, error_callback = None, progress_callback = None):
    """
    Initialize the execution callback functions. If a function is provided, it is returned unchanged. Otherwise, a
    default function is returned.

    Returns:
    --------
    warning_callback: function
        The callback function for displaying warnings to the user.
    error_callback: function
        The callback function for displaying errors to the user.
    progress_callback: function
        The callback for logging execution progress.

    """
    # Default callback functions
    if warning_callback is None:
        warning_callback = lambda w: logging.warning(w)
    if error_callback is None:
        def normal_raise(exception):
            raise exception
        error_callback = normal_raise
    if progress_callback is None:
        progress_callback = lambda t, p: None
    return warning_callback, error_callback, progress_callback def _unpack_binary_bytes_from_ints(a):
    """
    Unpacks binary values stored in bytes into ints
    """
    type = a.dtype

    if type == np.uint32:
        pack_size = 32
    elif type == np.uint64:
        pack_size = 64
    else:
        raise ValueError("Supported data types are 32-bit and 64-bit integers.")

    unpacked_n_rows = a.shape[0] * pack_size
    unpacked_n_columns = a.shape[1] if len(a.shape) > 1 else 1
    b = np.zeros((unpacked_n_rows, a.shape[1]) if len(a.shape) > 1 else unpacked_n_rows, dtype=np.uint8)

    packed_rows = 0
    packing_row = 0
    for i in xrange(b.shape[0]):
        if packed_rows == pack_size:
            packed_rows = 0
            packing_row += 1
        tmp = np.left_shift(np.ones(unpacked_n_columns, dtype=type), pack_size - (i - pack_size * packing_row)-1)
        np.bitwise_and(a[packing_row], tmp, tmp)
        b[i] = tmp > 0
        packed_rows += 1

    return b def _parse_kmer_blacklist(blacklist_path, expected_kmer_len):
    data = []
    
    # Fasta file format
    fasta_extensions = [".fasta", ".fa", ".fas", ".fna"]
    if any(blacklist_path.endswith(extension) for extension in fasta_extensions):
        data = _fasta_to_sequences(blacklist_path)
            
    # Other file format (one kmer per line)
    else:
        # Loading data and splitting on every line
        data = [l.rstrip('\n') for l in open(blacklist_path, "r")]
        
        # Filtering for empty strings (text file with one kmer per line)
        data = [x for x in data if x]
        
    def is_valid_kmer(k):
        return len(set(k).difference(["A", "C", "G", "T", "a", "c", "g", "t"])) == 0
    for kmer in data:
        if not is_valid_kmer(kmer):
            raise ValueError("{} is not a valid DNA sequence".format(kmer))
        
    if not(all(len(kmer) == expected_kmer_len for kmer in data)):
        raise ValueError("Extracted k-mers to blacklist do not have all the same length as the dataset k-mers")
    return data def _get_binary_metrics(predictions, answers):
    """
    Compute metrics for binary predictions.
    Parameters:
    -----------
    predictions: numpy_array, shape=(n_predictions, n_examples)
	    The predictions array.
    answers: numpy_array, shape=(n_examples,)
	    The labels vector associated to the examples.
    Returns:
    --------
    metrics: dictionary
	    A dictionary containing for each metric key,
	    a vector of values associated to the predictions.
    """
    if len(predictions.shape) == 1:
	    predictions = predictions.reshape(1, -1)
    metrics = defaultdict(list)
    for i in xrange(predictions.shape[0]):
	p = predictions[i]
	risk = 1.0 * len(p[p != answers]) / len(answers)
	tp = len(np.where(p[answers == 1] == 1)[0])
	fp = len(np.where(p[answers == 0] == 1)[0])
	tn = len(np.where(p[answers == 0] == 0)[0])
	fn = len(np.where(p[answers == 1] == 0)[0])
	precision = 1.0 * tp / (tp + fp) if (tp + fp) != 0 else -np.infty
	sensitivity = recall = 1.0 * tp / (tp + fn) if (tp + fn) != 0 else -np.infty
	specificity = 1.0 * tn / (fp + tn) if (fp + tn) != 0 else -np.infty
	f1_score = 2.0 * precision * recall / (precision + recall) if (precision + recall) > 0.0 else -np.infty
	metrics["risk"].append(risk)
	metrics["tp"].append(tp)
	metrics["fp"].append(fp)
	metrics["tn"].append(tn)
	metrics["fn"].append(fn)
	metrics["precision"].append(precision)
	metrics["sensitivity"].append(sensitivity)
	metrics["recall"].append(recall)
	metrics["specificity"].append(specificity)
	metrics["f1_score"].append(f1_score)
    return metrics def _get_multiclass_metrics(predictions, answers, nb_class):
    """
    Compute metrics for multi-class predictions.
    Parameters:
    -----------
    predictions: numpy_array, shape=(n_predictions, n_examples)
	    The predictions array.
    answers: numpy_array, shape=(n_examples,)
	    The labels vector associated to the examples.
    Returns:
    --------
    metrics: dictionary
	    A dictionary containing for each metric key,
	    a vector of values associated to the predictions.
    """
    if len(predictions.shape) == 1:
	    predictions = predictions.reshape(1, -1)
    metrics = defaultdict(list)
    for i in xrange(predictions.shape[0]):
	    p = predictions[i]
	    risk = 1.0 * len(p[p != answers]) / len(answers)
	    confusion_matrix = []
	    for actual_class in range(nb_class):
		    confusion_matrix.append([len(np.where(p[answers == actual_class] == predicted_class)[0]) \
										for predicted_class in range(nb_class)])
	    metrics["risk"].append(risk)
	    metrics["confusion_matrix"].append(confusion_matrix)
    return metrics
import logging
import numpy as np
from collections import defaultdict
from copy import deepcopy
from functools import partial
from itertools import product
from math import ceil, exp, log as ln, pi, sqrt
from multiprocessing import Pool, cpu_count
from scipy.misc import comb
from ...dataset.ds import KoverDataset
from ..learners.cart import DecisionTreeClassifier, _prune_tree
from ..common.models import CARTModel
from ..common.rules import LazyKmerRuleList, KmerRuleClassifications
from ...utils import _duplicate_last_element, _init_callback_functions, _unpack_binary_bytes_from_ints, _parse_kmer_blacklist
from ..experiments.metrics import _get_binary_metrics, _get_multiclass_metrics
    Copyright (C) 2015  Alexandre Drouin & Gael Letarte St-Pierre

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""




TIEBREAKER_RISK_BLOCK_SIZE = 100000


class BetweenDict(dict):
    """
    The BetweenDict
    By: Joshua Kugler
    Source: http://joshuakugler.com/archives/30-BetweenDict,-a-Python-dict-for-value-ranges.html
    """
    def __init__(self, d = {}):
        for k,v in d.items():
            self[k] = v

    def __getitem__(self, key):
        for k, v in self.items():
            if k[0] <= key < k[1] or (k[0] <= key and k[1] == np.infty) or (k[0] == -np.infty and key < k[1]):
                return v
        raise KeyError("Key '%s' is not between any values in the BetweenDict" % key)

    def __setitem__(self, key, value):
        try:
            if len(key) == 2:
                if key[0] < key[1]:
                    dict.__setitem__(self, (key[0], key[1]), value)
                else:
                    raise RuntimeError('First element of a BetweenDict key '
                                       'must be strictly less than the '
                                       'second element. Got [%.6f, %.6f]' % (key[0], key[1]))
            else:
                raise ValueError('Key of a BetweenDict must be an iterable '
                                 'with length two')
        except TypeError:
            raise TypeError('Key of a BetweenDict must be an iterable '
                             'with length two')

    def __contains__(self, key):
        try:
            return bool(self[key]) or True
        except KeyError:
            return False


def _tiebreaker(best_score_idx, rule_kmer_occurrences):
    """
    The tiebreaker finds the rules in best_score_idx for which the k-mer has the
    largest number of occurrences in the training set. This makes sense because
    it should avoid to select k-mers that isolate a small set of genomes (e.g.,
    a phylogenetic effect) when possible. It should also help find small
    compression sets in the sample compression bound.

    """
    logging.debug("There are %d candidate rules." % best_score_idx.shape[0])
    tie_rule_kmer_occurrences = rule_kmer_occurrences[best_score_idx]
    result = best_score_idx[np.isclose(tie_rule_kmer_occurrences, tie_rule_kmer_occurrences.max())]
    return result


def _split_callback(node, equivalent_rules_idx):
    """
    Called each time a leaf is split into two new leaves. This function is used
    to keep track of the number of equivalent rules at each decision point. It
    simply injects a list of equivalent rule indices into the rule used to split.
    It's kinda hacky, but also really simple.

    """
    logging.debug("Injecting equivalent rules into new tree node")
    node.rule.equivalent_rules_idx = equivalent_rules_idx


def _readdress_tree(tree, rule_new_idx_by_kmer_seq):
    def _readdress(node, kmer_idx):
        if node.rule is not None:
            node.rule.kmer_index = kmer_idx[node.rule.kmer_sequence]
            _readdress(node.left_child, kmer_idx)
            _readdress(node.right_child, kmer_idx)
    new_tree = deepcopy(tree)
    _readdress(new_tree, rule_new_idx_by_kmer_seq)
    return new_tree


def _predictions(decision_tree, kmer_matrix, train_example_idx, test_example_idx, progress_callback=None):
    """
    Makes predictions by loading only the columns of the kmer matrix that are targetted by the model.

    """
    # TODO check if this has the same bug that we fixed for kover SCM
    if progress_callback is None:
        progress_callback = lambda t, p: None
    progress_callback("Testing", 0.0)

    # Case: Standart tree
    if len(decision_tree.rules) > 0:
        model_rules = decision_tree.rules
        kmer_idx_by_rule = np.array([r.kmer_index for r in model_rules])
        kmer_sequence_by_rule = np.array([r.kmer_sequence for r in model_rules])
        sort_by_idx = np.argsort(kmer_idx_by_rule)
        kmer_idx_by_rule = kmer_idx_by_rule[sort_by_idx]
        kmer_sequence_by_rule = kmer_sequence_by_rule[sort_by_idx]
        readdressed_kmer_idx_by_rule = dict((s, i) for i, s in enumerate(kmer_sequence_by_rule))
        readdressed_decision_tree = _readdress_tree(tree=decision_tree, rule_new_idx_by_kmer_seq=readdressed_kmer_idx_by_rule)

        X = np.vstack((_unpack_binary_bytes_from_ints(kmer_matrix[:, idx]) for idx in kmer_idx_by_rule)).T