<commit_before>#pragma once

#include "Vector3.h"
#include "Database.h"
#include "MagneticField.h"

#include <errno.h>
#include <fcntl.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/mman.h>
#include <sys/types.h>
#include <unistd.h>


namespace quimby {

class invalid_position : std::exception {
public:
    virtual const char* what() const throw() {
        return "Invalid position.";
    }
};

template<size_t N>
class HCube {
	Vector3f elements[N * N * N];

public:
	void init(MagneticField *field, const Vector3f &offsetKpc, float sizeKpc,
			float error, float threshold, size_t maxdepth, size_t depth,
			size_t &idx) {
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;

		const float s = sizeKpc / N;
		const Vector3f o = offsetKpc + Vector3f(s / 2.);
		
        size_t thisidx = idx;
        idx++;

		if (depth == maxdepth) {
			for (size_t iX = 0; iX < N; iX++) {
				for (size_t iY = 0; iY < N; iY++) {
					for (size_t iZ = 0; iZ < N; iZ++) {
						size_t idx_target = iX * N2 + iY * N + iZ;
						field->getField(o + Vector3f(iX *s, iY * s, iZ * s), elements[idx_target]);
					}
				}
			}
    	} else {
			for (size_t n = 0; n < N3; n++) {
				if (depth == 0)
					std::cout << "\n " << n;
				else if(depth == 1)
					std::cout << "." << n;
				std::cout.flush();
				HCube<N> *hc = this + (idx - thisidx);
				size_t i = n / N2;
				size_t j = (n % N2) / N;
				size_t k = n % N;
				size_t tmpidx = idx;
				hc->init(field, offsetKpc + Vector3f(i * s, j * s, k * s), s,
						error, threshold, maxdepth, depth + 1, tmpidx);
				Vector3f mean;
				if (hc->collapse(mean, error, threshold)) {
					setValue(i, j, k, mean);
				} else {
					setCube(i, j, k, idx - thisidx);
					idx = tmpidx;
				}
			}

		}
	}

	
	void init(Database *db, const Vector3f &offsetKpc, float sizeKpc,
			float error, float threshold, size_t maxdepth, size_t depth,
			size_t &idx, size_t target_depth) {
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;

		size_t desired_depth = int(log2(1024) / log2(N));
		size_t remaining_depth = 1 + maxdepth - depth;
		//std::cout << "init: desired_depth=" << desired_depth << ", remaining_depth=" << remaining_depth << ", offset="<< offsetKpc << ", idx=" << idx << std::endl;
		if (remaining_depth <= desired_depth) {
			size_t sample_depth = std::min(remaining_depth, desired_depth);
			size_t starting_depth = std::max(maxdepth - sample_depth, 0ul);

			//std::cout << "Sample. N=" << N << ", depth=" << depth << ", remaining_depth=" << remaining_depth << std::endl;
			// sample data with max depth resolution
			size_t n = N;
			for (size_t i = 1; i < sample_depth; i++)
				n *= N;
			size_t n3 = n*n*n;
			Vector3f *data = new Vector3f[n3];

			// zero
			memset(data, 0, sizeof(Vector3f) * n3);

			// create visitor
			Vector3f lower = offsetKpc;
			//Vector3f upper = lower + Vector3f(sizeKpc, sizeKpc, sizeKpc);

			//std::cout << "Sample: " << lower << " - " << sizeKpc << std::endl;
			SimpleSamplingVisitor v(data, n, lower, sizeKpc);
			//v.showProgress(true);

			db->accept(v);

			// call this->init with sampled data
			init(data, n, sizeKpc, Vector3f(0, 0, 0), sizeKpc, error, threshold, maxdepth, depth, idx, target_depth);

			delete[] data;
		} else {
			const float s = sizeKpc / N;
			size_t thisidx = idx;
			idx++;
			for (size_t n = 0; n < N3; n++) {
				if (depth == 0)
					std::cout << "\n " << n;
				else if(depth == 1)
					std::cout << "." << n;
				std::cout.flush();
				HCube<N> *hc = this + (idx - thisidx);
				size_t i = n / N2;
				size_t j = (n % N2) / N;
				size_t k = n % N;
				size_t tmpidx = idx;
				hc->init(db, offsetKpc + Vector3f(i * s, j * s, k * s), s,
						error, threshold, maxdepth, depth + 1, tmpidx, target_depth);
				Vector3f mean;
				if (hc->collapse(mean, error, threshold) || (depth > target_depth)) {
					setValue(i, j, k, mean);
				} else {
					setCube(i, j, k, idx - thisidx);
					idx = tmpidx;
				}
			}

		}
	}

	void init(const Vector3f *data, size_t dataN, float dataSize,
			const Vector3f &offsetKpc, float sizeKpc, float error,
			float threshold, size_t maxdepth, size_t depth, size_t &idx, size_t target_depth) {
		const float s = sizeKpc / N;
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;
		size_t thisidx = idx;
		idx++;
		if (depth == maxdepth) {
			//std::cout << "init data final" << std::endl;
			float dataStep = dataSize / dataN;
			size_t oX = offsetKpc.x / dataStep;
			size_t oY = offsetKpc.y / dataStep;
			size_t oZ = offsetKpc.z / dataStep;
			const size_t dataN2 = dataN * dataN;
			for (size_t iX = 0; iX < N; iX++) {
				for (size_t iY = 0; iY < N; iY++) {
					for (size_t iZ = 0; iZ < N; iZ++) {
						size_t idx_data = (iX + oX) * dataN2 + (iY + oY) * dataN + (iZ + oZ);
						size_t idx_target = iX * N2 + iY * N + iZ;
						elements[idx_target] = data[idx_data];
					}
				}
			}
		} else {

			for (size_t n = 0; n < N3; n++) {
				HCube<N> *hc = this + (idx - thisidx);
				size_t i = n / N2;
				size_t j = (n % N2) / N;
				size_t k = n % N;
				size_t tmpidx = idx;
				hc->init(data, dataN, dataSize,
						offsetKpc + Vector3f(i * s, j * s, k * s), s, error,
						threshold, maxdepth, depth + 1, tmpidx, target_depth);
				Vector3f mean;
				if (hc->collapse(mean, error, threshold) || (depth > target_depth)) {
					setValue(i, j, k, mean);
				} else {
					setCube(i, j, k, idx - thisidx);
					idx = tmpidx;
				}
			}

		}
	}

	void load(std::istream &in, size_t dataN, float dataSize,
			const Vector3f &offsetKpc, float sizeKpc, float threshold) {
		const float s = sizeKpc / N;
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;

		float dataStep = dataSize / dataN;
		size_t oX = offsetKpc.x / dataStep;
		size_t oY = offsetKpc.y / dataStep;
		size_t oZ = offsetKpc.z / dataStep;
		const size_t dataN2 = dataN * dataN;
		for (size_t iX = 0; iX < N; iX++) {
			for (size_t iY = 0; iY < N; iY++) {
				size_t offsetN = (iX + oX) * dataN2 + (iY + oY) * dataN + oZ;
				in.seekg(offsetN * sizeof(Vector3f), std::ios::beg);
				in.read((char *) (elements + (iX * N2 + iY * N)),
						N * sizeof(Vector3f));
			}
		}
//		const float t2 = threshold * threshold;
//		for (size_t n = 0; n < N3; n++) {<commit_msg>fix target depth<commit_after>#pragma once

#include "Vector3.h"
#include "Database.h"
#include "MagneticField.h"

#include <errno.h>
#include <fcntl.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/mman.h>
#include <sys/types.h>
#include <unistd.h>


namespace quimby {

class invalid_position : std::exception {
public:
    virtual const char* what() const throw() {
        return "Invalid position.";
    }
};

template<size_t N>
class HCube {
	Vector3f elements[N * N * N];

public:
	void init(MagneticField *field, const Vector3f &offsetKpc, float sizeKpc,
			float error, float threshold, size_t maxdepth, size_t depth,
			size_t &idx) {
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;

		const float s = sizeKpc / N;
		const Vector3f o = offsetKpc + Vector3f(s / 2.);
		
        size_t thisidx = idx;
        idx++;

		if (depth == maxdepth) {
			for (size_t iX = 0; iX < N; iX++) {
				for (size_t iY = 0; iY < N; iY++) {
					for (size_t iZ = 0; iZ < N; iZ++) {
						size_t idx_target = iX * N2 + iY * N + iZ;
						field->getField(o + Vector3f(iX *s, iY * s, iZ * s), elements[idx_target]);
					}
				}
			}
    	} else {
			for (size_t n = 0; n < N3; n++) {
				if (depth == 0)
					std::cout << "\n " << n;
				else if(depth == 1)
					std::cout << "." << n;
				std::cout.flush();
				HCube<N> *hc = this + (idx - thisidx);
				size_t i = n / N2;
				size_t j = (n % N2) / N;
				size_t k = n % N;
				size_t tmpidx = idx;
				hc->init(field, offsetKpc + Vector3f(i * s, j * s, k * s), s,
						error, threshold, maxdepth, depth + 1, tmpidx);
				Vector3f mean;
				if (hc->collapse(mean, error, threshold)) {
					setValue(i, j, k, mean);
				} else {
					setCube(i, j, k, idx - thisidx);
					idx = tmpidx;
				}
			}

		}
	}

	
	void init(Database *db, const Vector3f &offsetKpc, float sizeKpc,
			float error, float threshold, size_t maxdepth, size_t depth,
			size_t &idx, size_t target_depth) {
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;

		size_t desired_depth = int(log2(1024) / log2(N));
		size_t remaining_depth = 1 + maxdepth - depth;
		//std::cout << "init: desired_depth=" << desired_depth << ", remaining_depth=" << remaining_depth << ", offset="<< offsetKpc << ", idx=" << idx << std::endl;
		if (remaining_depth <= desired_depth) {
			size_t sample_depth = std::min(remaining_depth, desired_depth);
			size_t starting_depth = std::max(maxdepth - sample_depth, 0ul);

			//std::cout << "Sample. N=" << N << ", depth=" << depth << ", remaining_depth=" << remaining_depth << std::endl;
			// sample data with max depth resolution
			size_t n = N;
			for (size_t i = 1; i < sample_depth; i++)
				n *= N;
			size_t n3 = n*n*n;
			Vector3f *data = new Vector3f[n3];

			// zero
			memset(data, 0, sizeof(Vector3f) * n3);

			// create visitor
			Vector3f lower = offsetKpc;
			//Vector3f upper = lower + Vector3f(sizeKpc, sizeKpc, sizeKpc);

			//std::cout << "Sample: " << lower << " - " << sizeKpc << std::endl;
			SimpleSamplingVisitor v(data, n, lower, sizeKpc);
			//v.showProgress(true);

			db->accept(v);

			// call this->init with sampled data
			init(data, n, sizeKpc, Vector3f(0, 0, 0), sizeKpc, error, threshold, maxdepth, depth, idx, target_depth);

			delete[] data;
		} else {
			const float s = sizeKpc / N;
			size_t thisidx = idx;
			idx++;
			for (size_t n = 0; n < N3; n++) {
				if (depth == 0)
					std::cout << "\n " << n;
				else if(depth == 1)
					std::cout << "." << n;
				std::cout.flush();
				HCube<N> *hc = this + (idx - thisidx);
				size_t i = n / N2;
				size_t j = (n % N2) / N;
				size_t k = n % N;
				size_t tmpidx = idx;
				hc->init(db, offsetKpc + Vector3f(i * s, j * s, k * s), s,
						error, threshold, maxdepth, depth + 1, tmpidx, target_depth);
				Vector3f mean;
				if (hc->collapse(mean, error, threshold) || (depth >= target_depth)) {
					setValue(i, j, k, mean);
				} else {
					setCube(i, j, k, idx - thisidx);
					idx = tmpidx;
				}
			}

		}
	}

	void init(const Vector3f *data, size_t dataN, float dataSize,
			const Vector3f &offsetKpc, float sizeKpc, float error,
			float threshold, size_t maxdepth, size_t depth, size_t &idx, size_t target_depth) {
		const float s = sizeKpc / N;
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;
		size_t thisidx = idx;
		idx++;
		if (depth == maxdepth) {
			//std::cout << "init data final" << std::endl;
			float dataStep = dataSize / dataN;
			size_t oX = offsetKpc.x / dataStep;
			size_t oY = offsetKpc.y / dataStep;
			size_t oZ = offsetKpc.z / dataStep;
			const size_t dataN2 = dataN * dataN;
			for (size_t iX = 0; iX < N; iX++) {
				for (size_t iY = 0; iY < N; iY++) {
					for (size_t iZ = 0; iZ < N; iZ++) {
						size_t idx_data = (iX + oX) * dataN2 + (iY + oY) * dataN + (iZ + oZ);
						size_t idx_target = iX * N2 + iY * N + iZ;
						elements[idx_target] = data[idx_data];
					}
				}
			}
		} else {

			for (size_t n = 0; n < N3; n++) {
				HCube<N> *hc = this + (idx - thisidx);
				size_t i = n / N2;
				size_t j = (n % N2) / N;
				size_t k = n % N;
				size_t tmpidx = idx;
				hc->init(data, dataN, dataSize,
						offsetKpc + Vector3f(i * s, j * s, k * s), s, error,
						threshold, maxdepth, depth + 1, tmpidx, target_depth);
				Vector3f mean;
				if (hc->collapse(mean, error, threshold) || (depth >= target_depth)) {
					setValue(i, j, k, mean);
				} else {
					setCube(i, j, k, idx - thisidx);
					idx = tmpidx;
				}
			}

		}
	}

	void load(std::istream &in, size_t dataN, float dataSize,
			const Vector3f &offsetKpc, float sizeKpc, float threshold) {
		const float s = sizeKpc / N;
		const size_t N3 = N * N * N;
		const size_t N2 = N * N;

		float dataStep = dataSize / dataN;
		size_t oX = offsetKpc.x / dataStep;
		size_t oY = offsetKpc.y / dataStep;
		size_t oZ = offsetKpc.z / dataStep;
		const size_t dataN2 = dataN * dataN;
		for (size_t iX = 0; iX < N; iX++) {
			for (size_t iY = 0; iY < N; iY++) {
				size_t offsetN = (iX + oX) * dataN2 + (iY + oY) * dataN + oZ;
				in.seekg(offsetN * sizeof(Vector3f), std::ios::beg);
				in.read((char *) (elements + (iX * N2 + iY * N)),
						N * sizeof(Vector3f));
			}
		}
//		const float t2 = threshold * threshold;
//		for (size_t n = 0; n < N3; n++) {
{'fraction_non_alphanumeric': 0.16076141725985763, 'fraction_numerical': 0.013916660001599615, 'mean_word_length': 2.6191027496382056, 'pattern_counts': {'":': 0, '<': 181, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5249244', 'n_tokens_mistral': 5851, 'n_tokens_neox': 4984, 'n_words': 1434}