hnsrtreefileobj.cc

R 树

	min.count = -1;
	min.leftVar = 0;
	min.rightVar = 0;

	for(count=minCount; count<=maxCount; count++) {
		double leftVar, rightVar;

		/* calculate the sum of variances on the left side */
		leftVar = 0;

		for(axis=0; axis<dimension; axis++) {
			double sum, sum2, avg, var;

			sum = 0;
			sum2 = 0;

			for(i=0; i<count; i++) {
				double coord = points[array[i]].getCoord(axis);
				sum += coord;
				sum2 += coord * coord;
			}

			avg = sum / count;
			var = sum2 / count - avg * avg;

			leftVar += var;
		}

		/* calculate the sum of variances on the right side */
		rightVar = 0;

		for(axis=0; axis<dimension; axis++) {
			double sum, sum2, avg, var;

			sum = 0;
			sum2 = 0;

			for(i=count; i<numPoints; i++) {
				double coord = points[array[i]].getCoord(axis);
				sum += coord;
				sum2 += coord * coord;
			}
			avg = sum / (numPoints - count);
			var = sum2 / (numPoints - count) - avg * avg;

			rightVar += var;
		}

		if(min.count < 0) {
			min.count = count;
			min.leftVar = leftVar;
			min.rightVar = rightVar;
		}
		else {
			if(leftVar + rightVar < min.leftVar + min.rightVar) {
				min.count = count;
				min.leftVar = leftVar;
				min.rightVar = rightVar;
			}
		}
	}

	/* make split flags */
	i=0;
	while(i<min.count) {
		flags[array[i]] = LEFT;
		i++;
	}
	while(i<numPoints) {
		flags[array[i]] = RIGHT;
		i++;
	}

	if(debug) {
		fprintf(stderr, "HnSRTreeFileObj::splitPoints\n");
		fprintf(stderr, "    LEFT\n");
		for(i=0; i<numPoints; i++) {
			if(flags[i] == LEFT)
				fprintf(stderr,
					"    points[%d].getCoord(%d) = %g\n",
					i, max.axis,
					points[i].getCoord(max.axis));
		}
		fprintf(stderr, "    RIGHT\n");
		for(i=0; i<numPoints; i++) {
			if(flags[i] == RIGHT)
				fprintf(stderr,
					"    points[%d].getCoord(%d) = %g\n",
					i, max.axis,
					points[i].getCoord(max.axis));
		}
	}

	*flags_return = flags;

	delete array;
}

void
HnSRTreeFileObj::splitCores(const HnSRTreeCoreArray &cores,
			    SplitFlag **flags_return)
{
	static SplitFlag *flags = NULL;
	int numCores = cores.length();
	int minCount = numCores * splitFactor / 100;
	int maxCount = numCores - minCount;
	struct {
		int axis;
		double var;
	} max;
	struct {
		int count;
		double leftVar, rightVar;
	} min;
	int axis, count;
	int i, j;
	int *array = new int[numCores];

	/* initialize flags */
	if(flags != NULL)
		delete flags;
	flags = new SplitFlag[numCores];

	for(i=0; i<numCores; i++)
		flags[i] = SPLIT_NONE;

	max.axis = -1;
	max.var = -1;

	/* calculate variance */
	for(axis=0; axis<dimension; axis++) {
		double avg, var;
		double sum = 0;
		double sum2 = 0;
		int total = 0;

		for(i=0; i<numCores; i++) {
			HnPoint center = cores[i].getCenter();
			double coord = center.getCoord(axis);

			/*
			 * NOTE:
			 *	In the calculation of the variance,
			 *	it is supposed that all points are located
			 *	at the center of the sphere.
			 */

			sum += coord * cores[i].getWeight();
			sum2 += coord * coord * cores[i].getWeight();
			total += cores[i].getWeight();
		}

		avg = sum / total;
		var = sum2 / total - avg * avg;

		if(max.axis == -1 || var > max.var) {
			max.axis = axis;
			max.var = var;
		}
	}

	/* sort points */
	for(i=0; i<numCores; i++)
		array[i] = i;

	for(i=0; i<numCores; i++) {
		for(j=i+1; j<numCores; j++) {
			HnPoint center1 = cores[array[i]].getCenter();
			HnPoint center2 = cores[array[j]].getCenter();
			double coord1 =	center1.getCoord(max.axis);
			double coord2 = center2.getCoord(max.axis);

			if(coord1 > coord2) {
				int temp = array[i];
				array[i] = array[j];
				array[j] = temp;
			}
		}
	}

	/* choose split point */
	min.count = -1;
	min.leftVar = 0;
	min.rightVar = 0;

	for(count=minCount; count<=maxCount; count++) {
		double leftVar, rightVar;

		/* calculate the sum of variances on the left side */
		leftVar = 0;

		for(axis=0; axis<dimension; axis++) {
			double sum, sum2, avg, var;
			int total;

			/*
			 * NOTE:
			 *	In the calculation of the variance,
			 *	it is supposed that all points are located
			 *	at the center of the sphere.
			 */

			sum = 0;
			sum2 = 0;
			total = 0;

			for(i=0; i<count; i++) {
				HnPoint center = cores[array[i]].getCenter();
				double coord = center.getCoord(axis);

				sum += coord * cores[array[i]].getWeight();
				sum2 += coord * coord
					* cores[array[i]].getWeight();
				total += cores[array[i]].getWeight();
			}

			avg = sum / total;
			var = sum2 / total - avg * avg;

			leftVar += var;
		}

		/* calculate the sum of variances on the right side */
		rightVar = 0;

		for(axis=0; axis<dimension; axis++) {
			double sum, sum2, avg, var;
			int total;

			/*
			 * NOTE:
			 *	In the calculation of the variance,
			 *	it is supposed that all points are located
			 *	at the center of the sphere.
			 */

			sum = 0;
			sum2 = 0;
			total = 0;

			for(i=count; i<numCores; i++) {
				HnPoint center = cores[array[i]].getCenter();
				double coord = center.getCoord(axis);

				sum += coord * cores[array[i]].getWeight();
				sum2 += coord * coord
					* cores[array[i]].getWeight();
				total += cores[array[i]].getWeight();
			}

			avg = sum / total;
			var = sum2 / total - avg * avg;

			rightVar += var;
		}

		if(min.count < 0) {
			min.count = count;
			min.leftVar = leftVar;
			min.rightVar = rightVar;
		}
		else {
			if(leftVar + rightVar < min.leftVar + min.rightVar) {
				min.count = count;
				min.leftVar = leftVar;
				min.rightVar = rightVar;
			}
		}
	}

	/* make split flags */
	i=0;
	while(i<min.count) {
		flags[array[i]] = LEFT;
		i++;
	}
	while(i<numCores) {
		flags[array[i]] = RIGHT;
		i++;
	}

	if(debug) {
		fprintf(stderr, "HnSRTreeFileObj::selectSpheres:\n");
		fprintf(stderr, "    LEFT\n");
		for(i=0; i<numCores; i++) {
			if(flags[i] == LEFT)
				fprintf(stderr,
					"    cores[%d].getCenter()."
					"getCoord(%d) = %g\n",
					i, max.axis,
					cores[i].getCenter().
					getCoord(max.axis));
		}
		fprintf(stderr, "    RIGHT\n");
		for(i=0; i<numCores; i++) {
			if(flags[i] == RIGHT)
				fprintf(stderr,
					"    cores[%d].getCenter()."
					"getCoord(%d) = %g\n",
					i, max.axis,
					cores[i].getCenter().
					getCoord(max.axis));
		}
	}

	*flags_return = flags;

	delete array;
}

/*
 * Construction
 */

void
HnSRTreeFileObj::getSplitFactors(const HnPointArray &points,
				 const HnDataArray &values,
				 int *numLevels_return, int *numSplits_return,
				 int *numRootSplits_return,
				 int *numNodeSplits_return)
{
	int i;
	int maxLeafCount, maxNodeCount, maxSplitsPerNode;
	int numLeaves, numSplits, numLevels;
	int numRootSplits, numNodeSplits;

	for(i=0; i<values.length(); i++) {
		if(values[i].length() > dataSize) {
			HnAbort("The size of data[%d] (%d) "
				"exceeds the limit (%d).",
				i, values[i].length(), dataSize);
		}
	}

	maxLeafCount =
		HnSRTreeLeaf::getMaxCount(dimension, dataSize, blockSize);
	maxNodeCount = HnSRTreeNode::getMaxCount(dimension, blockSize);

	numLeaves = (points.length() - 1) / maxLeafCount + 1;

	maxSplitsPerNode = (int)(log(maxNodeCount) / log(2));
	numSplits = (int)ceil(log(numLeaves) / log(2));

	numLevels = ((numSplits - 1) / maxSplitsPerNode + 1) + 1;

	if((numRootSplits = numSplits % maxSplitsPerNode) == 0)
		numRootSplits = maxSplitsPerNode;
	numNodeSplits = maxSplitsPerNode;

	*numLevels_return = numLevels;
	*numSplits_return = numSplits;
	*numRootSplits_return = numRootSplits;
	*numNodeSplits_return = numNodeSplits;

	if(debug) {
		fprintf(stderr, "HnSRTreeFileObj::getSplitFactors: "
			"maxLeafCount = %d, maxNodeCount = %d, "
			"maxSplitsPerNode = %d, numLeaves = %d, "
			"numSplits = %d, numLevels = %d, "
			"numRootSplits = %d, numNodeSplits = %d\n",
			maxLeafCount, maxNodeCount, maxSplitsPerNode,
			numLeaves, numSplits, numLevels,
			numRootSplits, numNodeSplits);
	}
}

void
HnSRTreeFileObj::constructTree(HnPointArray &points, HnDataArray &values,
			       int offset, int count,
			       int levelCount, int numLevels,
			       int splitCount, int numSplits,
			       int numRootSplits, int numNodeSplits,
			       HnSRTreeNode &node)
{
	int axis;

	axis = getConstructionAxis(points, offset, count);

	sortPoints(points, values, offset, count, axis);

	if(debug) {
		int i;

		for(i=offset; i<offset + count - 1; i++) {
			if(points[i].getCoord(axis) >
			   points[i+1].getCoord(axis))
				HnAbort("sorting failed.");
		}
	}

	if(splitCount + 1 == numSplits) {
		HnSRTreeLeaf leftLeaf = allocateLeaf();
		HnSRTreeLeaf rightLeaf = allocateLeaf();
		int index;

		index = offset;
		while(index < offset + count / 2) {
			leftLeaf.append(points[index], values[index]);
			index ++;
		}
		while(index < offset + count) {
			rightLeaf.append(points[index], values[index]);
			index ++;
		}

		node.append(leftLeaf.getCore(),
			    leftLeaf.getBoundingRect(),
			    leftLeaf.getOffset());
		node.append(rightLeaf.getCore(),
			    rightLeaf.getBoundingRect(),
			    rightLeaf.getOffset());

		writeLeaf(leftLeaf);
		writeLeaf(rightLeaf);
	}
	else if(splitCount + 1 == numRootSplits + levelCount * numNodeSplits) {
		HnSRTreeNode leftNode = allocateNode();
		HnSRTreeNode rightNode = allocateNode();

		constructTree(points, values,
			      offset, count / 2,
			      levelCount + 1, numLevels,
			      splitCount + 1, numSplits,
			      numRootSplits, numNodeSplits,
			      leftNode);
		constructTree(points, values,
			      offset + count / 2, count - count / 2,
			      levelCount + 1, numLevels, 
			      splitCount + 1, numSplits,
			      numRootSplits, numNodeSplits,
			      rightNode);

		node.append(leftNode.getCore(),
			    leftNode.getBoundingRect(),
			    leftNode.getOffset());
		node.append(rightNode.getCore(),
			    rightNode.getBoundingRect(),
			    rightNode.getOffset());

		writeNode(leftNode);
		writeNode(rightNode);
	}
	else {
		constructTree(points, values,
			      offset, count / 2,
			      levelCount, numLevels,
			      splitCount + 1, numSplits,
			      numRootSplits, numNodeSplits,
			      node);
		constructTree(points, values,
			      offset + count / 2, count - count / 2,
			      levelCount, numLevels,
			      splitCount + 1, numSplits,
			      numRootSplits, numNodeSplits,
			      node);
	}
}

int
HnSRTreeFileObj::getConstructionAxis(HnPointArray &points,
				     int offset, int count)
{
	int axis;
	double *sum = new double[dimension];
	double *sum2 = new double[dimension];
	int i;
	struct {
		int axis;
		double var;
	} max;

	for(axis=0; axis<dimension; axis++) {
		sum[axis] = 0;
		sum2[axis] = 0;
	}

	for(i=offset; i<offset + count; i++) {
		for(axis=0; axis<dimension; axis++) {
			double x = points[i].getCoord(axis);

			sum[axis] += x;
			sum2[axis] += x * x;
		}
	}

	max.axis = -1;
	max.var = 0;

	for(axis=0; axis<dimension; axis++) {
		double mean = sum[axis] / count;
		double var = sum2[axis] / count - mean * mean;

		if(axis == 0) {
			max.axis = 0;
			max.var = var;
		}
		else {
			if(var > max.var) {
				max.axis = axis;
				max.var = var;
			}
		}
	}

	delete sum;
	delete sum2;

	return max.axis;
}

void
HnSRTreeFileObj::sortPoints(HnPointArray &points, HnDataArray &values,
			    int offset, int count, int axis)
{
	if(count < 16) {
		int i, j;

		for(i=offset; i<offset + count; i++) {
			for(j=i+1; j<offset + count; j++) {
				if(points[i].getCoord(axis) >
				   points[j].getCoord(axis)) {
					points.swap(i, j);
					values.swap(i, j);
				}
			}
		}
	}
	else {
		double pivot =
			points[offset + count / 2].getCoord(axis);
		int leftIndex, rightIndex;

		points.swap(offset, offset + count / 2);
		values.swap(offset, offset + count / 2);

		leftIndex = offset;
		rightIndex = offset + count;

		for(;;) {
			do {
				leftIndex ++;
			} while(leftIndex < offset + count &&
				points[leftIndex].getCoord(axis)
				<= pivot);
			do {
				rightIndex --;
			} while(rightIndex >= offset && 
				points[rightIndex].getCoord(axis)
				> pivot);

			if(leftIndex > rightIndex)
				break;

			if(leftIndex == rightIndex)
				HnAbort("impossible situation.");

			points.swap(leftIndex, rightIndex);
			values.swap(leftIndex, rightIndex);
		}

		if(leftIndex != rightIndex + 1)
			HnAbort("impossible situation.");

		if(offset != leftIndex - 1) {
			points.swap(offset, leftIndex - 1);
			values.swap(offset, leftIndex - 1);
		}

		if(debug) {
			int i;

			for(i=offset; i<leftIndex - 1; i++) {
				if(points[i].getCoord(axis) > pivot)
					HnAbort("left pivot is failed.");
			}
			if(points[leftIndex-1].getCoord(axis) != pivot)
				HnAbort("center pivot is failed.");
			for(i=leftIndex; i<offset + count; i++) {
				if(points[i].getCoord(axis) <= pivot)
					HnAbort("right pivot is failed.");
			}
		}

		sortPoints(points, values,
			   offset, leftIndex - 1 - offset, axis);
		sortPoints(points, values,
			   leftIndex, offset + count - leftIndex, axis);
	}
}

/*
 * Remove
 */

void
HnSRTreeFileObj::remove(const HnPoint &point, const HnData &data)
{
	HnSRTreeStack stack;
	HnSRTreeLeaf leaf;