stim/envi/binary.h

//make sure that this header file is only loaded once
#ifndef RTS_BINARY_H
#define RTS_BINARY_H
#include "../envi/envi_header.h"
#include "../math/vector.h"
#include <fstream>
#include <sys/stat.h>
#include <cstring>
#include <chrono>
#ifdef _WIN32
#include <Windows.h>
#else
#include <unistd.h>
#endif
#ifdef CUDA_FOUND
//CUDA externs
void gpu_permute(char* dest, char* src, size_t sx, size_t sy, size_t sz, size_t d0, size_t d1, size_t d2, size_t typesize);
#include <stim/cuda/cudatools/error.h>
#endif
namespace stim{
/// This class calculates the optimal setting for independent parameter b (batch size) for
///		minimizing the dependent parameter bps (bytes per second)
class stream_optimizer{
protected:
	size_t Bps[2];							//bytes per second for the previous batch
	size_t interval_B;					//number of bytes processed this interval
	size_t interval_ms;					//number of milliseconds spent in the current interval
	size_t n[2];							//current batch size (in bytes)
	size_t h;							//spacing used for finite difference calculations
	size_t dn;							//delta value (in bytes) for setting the batch size (minimum change in batch parameter)
	size_t maxn;						//maximum value for the batch size
	double alpha;						//alpha value controls the factor of the gradient that is used to calculate the next point (speed of convergence)
	bool sample_step;					//calculating the derivative (this class alternates between calculating dBps and B)
	bool forward_diff;					//evaluate the derivative using forward differences
	size_t window_ms;					//size of the interval (in milliseconds) integrated to get a reliable bps value
	// This function rounds x to the nearest value within dB
	size_t round_limit(double n0){
		if(n0 < 0) return dn;					//if n0 is less than zero, return the lowest possible n
		size_t new_n = (size_t)(n0 + 0.5);		//now n0 must be positive, so round it to the nearest integer
		if(new_n > maxn) new_n = maxn;			//limit the returned size of x to within the specified bounds
		size_t lowest = new_n / dn;
		size_t highest = lowest + dn;
		size_t diff[2] = {new_n - lowest, highest - new_n};	//calculate the two differences
		if(diff[0] < diff[1])
			return lowest;
		return highest;
	}
public:
	//constructor initializes a stream optimizer
	stream_optimizer(size_t min_batch_size, size_t max_batch_size, double a = 0.003, size_t probe_step = 5, size_t window = 2000){
		//Bps = 0;						//initialize to zero bytes per second processed
		Bps[0] = Bps[1] = 0;			//initialize the bits per second to 0
		interval_B = 0;					//zero bytes have been processed at initialization
		interval_ms = 0;				//no time has been spent on the batch so far
		dn = min_batch_size;			//set the minimum batch size as the minimum change in batch size
		maxn = max_batch_size;			//set the maximum batch size
		n[0] = max_batch_size;			//set B
		h = (max_batch_size / min_batch_size) / probe_step * dn;
		std::cout<<"h = "<<h<<std::endl;
		if(h < dn) h = dn;
		alpha = a;
		//n[0] = round_limit( (max_batch_size - min_batch_size)/2 );
		window_ms = window;		//minimum integration interval (for getting a reliable bps measure)
		sample_step = true;					//the first step is to calculate the derivative
		forward_diff = true;			//start with the forward difference (since we start at the maximum batch size)
	}
	size_t update(size_t bytes_processed, size_t ms_spent, size_t& data_rate, bool VERBOSE = false){
		interval_B += bytes_processed;		//increment the number of bytes processed
		interval_ms += ms_spent;			//increment the number of milliseconds spent processing
		data_rate = interval_B / interval_ms;
		//if we have sufficient information to evaluate the optimization function at this point
		if(interval_ms < window_ms){					//if insufficient time has passed to get a reliable Bps measurement
			return n[0];
		}
		else{											//if we have collected enough information for a reliable Bps estimate
			
			if(Bps[0] == 0){							//if n[0] hasn't been evaluated yet, this is the first step
				Bps[0] = data_rate;						//set the initial Bps value
				n[1] = n[0] - h;						//set the position of the next sample point
				if(VERBOSE)
					std::cout<<"Bps value at n = "<<n[0]<<" is "<<Bps[0]<<" Bps, probing n = "<<n[1]<<std::endl;
				return n[1];							//return the probe point
			}
			else{
				Bps[1] = data_rate;						//set the Bps for the current point (n[1])
				double Bps_p;							//allocate a variable for the derivative
				//calculate the derivative
				if(n[0] < n[1]){						//if the current point is less than the previous one (probably the most common)
					Bps_p = ((double)Bps[1] - (double)Bps[0]) / (double)h;		//calculate the derivative using the forward finite difference
				}
				else{
					Bps_p = ((double)Bps[0] - (double)Bps[1]) / (double)h;		//calculate the derivative using the backward finite difference
				}
				if(VERBOSE)
					std::cout<<"     probed n = "<<n[1]<<" with "<<Bps[1]<<" Bps, gradient = "<<Bps_p<<" Bps"<<std::endl;
				double new_n_precise = n[0] + alpha * Bps_p;			//calculate the next point (snap to closest integer)
				size_t new_n_nearest = round_limit(new_n_precise);		//calculate the next point (given batch parameters)
				if(new_n_nearest == n[0]){								//if the newest point is the same as the original point
					Bps[0] = Bps[1];									//update the Bps
					//if(n[0] == dn) n[1] = n[0] + h;					//if we're on the left edge, probe forward
					//else n[1] = n[0] - h;								//otherwise probe backwards
					if(VERBOSE)
						std::cout<<"     staying at n = "<<n[0]<<" for now"<<std::endl;
					//return n[1];										//return the probe point
					Bps[0] = 0;											//reset the Bps for the current point
					return n[0];										//return the current point for a re-calculation
				}
				else{													//if the newest point is different from the original point
					n[0] = new_n_nearest;								//move to the new point
					Bps[0] = 0;											//set the Bps to zero (point hasn't been tested)
					if(VERBOSE)
						std::cout<<"     moving to n = "<<n[0]<<std::endl;
					return n[0];										//return the new point
				}
			}
		}
	}
	/*// this function updates the optimizer, given the number of bytes processed in an interval and time spent processing
	size_t update(size_t bytes_processed, size_t ms_spent){
		interval_B += bytes_processed;		//increment the number of bytes processed
		interval_ms += ms_spent;			//increment the number of milliseconds spent processing
		//if we have sufficient information to evaluate the optimization function at this point
		if(interval_ms >= window_ms){					//if sufficient time has passed to get a reliable Bps measurement
			size_t new_Bps = interval_B / interval_ms;	//calculate the current Bps
			if(sample_step){							//if this is a sample step, collect the information for Bps = f(n0)
				Bps = new_Bps;							//set the Bps to the evaluated value
				n[1] = n[0] - dn;								//reduce the batch size by one delta to take a second sample
				if(n[1] == 0){							//if the resulting batch size is zero
					n[1] = 2*dn;						//we're at the left edge: set the new sample point to 2*dn
				}
				interval_B = interval_ms = 0;			//start a new interval at the new sample point
				sample_step = false;					//next step will calculate the new batch size via optimization
				return n[1];								//return the new batch size
			}
			else{								//if we have sufficient information to evaluate the derivative and optimize
				double f = (double)new_Bps;				//we have evaluated the function at this location
				double fprime;
				if(n[1] < n[0] ){									//if the new point is less than the previous point (usually the case)
					fprime = (double)(Bps - new_Bps) / (double)dn;	//calculate the forward difference
				}
				else{												//if the new point is larger (only happens at the minimum limit)
					fprime = (double)(new_Bps - Bps) / (double)dn;	//calculate the backward difference
				}
				size_t bestn = n[1] - (size_t)(f / fprime);			//calculate the best value for B using Newton's method
				n[0] = round_limit( (size_t)bestn );						//set the new dependent point
				sample_step = true;									//the next step will be a sample step
			}
		}
		if(sample_step) return n[0];
		return n[1];										//insufficient information, keep the same batch size
	}*/
	/*size_t update(size_t bytes_processed, size_t ms_spent){
		interval_B += bytes_processed;		//increment the number of bytes processed
		interval_ms += ms_spent;			//increment the number of milliseconds spent processing
		//if( Bps[0] == 0 ){				//if the left boundary hasn't been processed
		//if we have sufficient information to evaluate the optimization function at this point
		if(interval_ms >= window_ms){
			size_t new_Bps = interval_B / interval_ms;	//calculate the current Bps
			if(Bps[0] == 0)							//if the left interval Bps hasn't been calculated
				Bps[0] = interval_B / interval_ms;	//that is the interval being processed
			else
				Bps[1] = interval_B / interval_ms;	//otherwise the right interval is being processed
			if(Bps[0] != 0 && Bps[1] != 0){			//if both intervals have been processed
		}
	}*/
	/*size_t update(size_t bytes_processed, size_t ms_spent, size_t& data_rate, bool VERBOSE){
		size_t time = update(bytes_processed, ms_spent, VERBOSE);
		data_rate = Bps[0];
		return time;
	}*/
};
/** This class manages the streaming of large multidimensional binary files.
 *  Generally these are hyperspectral files with 2 spatial and 1 spectral dimension. However, this class supports
 *  other dimensions via the template parameter D.
 *
 *  @param T is the data type used to store data to disk (generally float or double)
 *  @param D is the dimension of the data (default 3)
 */
template< typename T, unsigned int D = 3 >
class binary{
protected:
	std::fstream file;		//file stream used for reading and writing
	std::string name;		//file name
	unsigned long long R[D];		//resolution
	unsigned long long header;	//header size (in bytes)
	unsigned char* mask;	//pointer to a character array: 0 = background, 1 = foreground (or valid data)
	double progress;		//stores the progress on the current operation (accessible using a thread)
	size_t data_rate;		//data rate (currently in Bps)
	size_t buffer_size;		//available memory for processing large files
	/// Private initialization function used to set default parameters in the data structure.
	void init(){
		std::memset(R, 0, sizeof(unsigned long long) * D);		//initialize the resolution to zero
		header = 0;											//initialize the header size to zero
		mask = NULL;
		progress = 0;										//initialize the progress for any algorithm to zero
		data_rate = 0;										//initialize the data rate to zero
		set_buffer_frac();										//set the maximum buffer size to the default
	}
	/// Private helper function that returns the size of the file on disk using system functions.
	long long int get_file_size(){
#ifdef _WIN32
		struct _stat64 results;
		if(_stat64(name.c_str(), &results) == 0)
			return results.st_size;
#else
		struct stat results;
		if(stat(name.c_str(), &results) == 0)
			return results.st_size;
#endif
		else return 0;
	}
	/// Private helper function that tests to make sure that the calculated data size specified by the structure is the same as the data size on disk.
	bool test_file_size(){
		long long int npts = 1;				//initialize the number of data points to 1
		for(unsigned int i = 0; i<D; i++)	//iterate over each dimension
			npts *= R[i];					//compute the total number of data points in the file
		long long int datasize = npts * sizeof(T);//multiply the sum by the size of the template parameter
		if(datasize + header == get_file_size()) return true;	//if the byte size matches the file size, we're golden
		else return false;					//otherwise return an error
	}
	/// Private helper function that resets the file pointer to the beginning of the data
	void reset(){
		file.seekg(header, std::ios_base::beg);
	}
	/// Private helper file that opens a specified binary file.
	/// @param filename is the name of the binary file to stream
	bool open_file(std::string filename){
		//open the file as binary for reading and writing
		file.open(filename.c_str(), std::ios::in | std::ios::out | std::ios::binary);
		//if the file isn't open, the user may only have read access
		if(!file.is_open()){
			std::cout<<"class STIM::BINARY - failed to open file, trying for read only"<<std::endl;
			file.open(filename.c_str(), std::ios::in | std::ios::binary);
			if(!file.is_open()){
				std::cout<<"               still unable to load the file"<<std::endl;
				return false;
			}
		}
		//if the file is successful
		if(file){
			name = filename;		//set the name
			if(test_file_size())	//test the file size
				return true;
		}
		return false;
	}
public:
	//default constructor
	binary(){
		init();
	}
	double get_progress(){
		return progress;
	}
	void reset_progress(){
		progress = 0;
	}
	size_t get_data_rate(){
		return data_rate;
	}
	//specify the maximum fraction of available memory that this class will use for buffering
	void set_buffer_frac(double mem_frac = 0.5){				//default to 50%
#ifdef _WIN32
		MEMORYSTATUSEX statex;
		statex.dwLength = sizeof (statex);
		GlobalMemoryStatusEx (&statex);
		buffer_size = (size_t)(statex.ullAvailPhys * mem_frac);
#else
		size_t pages = sysconf(_SC_PHYS_PAGES);
		size_t page_size = sysconf(_SC_PAGE_SIZE);
		buffer_size = (size_t)(pages * page_size *  mem_frac);
#endif
	}
	void set_buffer_raw(size_t bytes){
		buffer_size = bytes;
	}
	/// Open a binary file for streaming.
	/// @param filename is the name of the binary file
	/// @param r is a STIM vector specifying the size of the binary file along each dimension
	/// @param h is the length (in bytes) of any header file (default zero)
	bool open(std::string filename, vec<unsigned long long> r, unsigned long long h = 0){
		for(unsigned long long i = 0; i < D; i++)		//set the dimensions of the binary file object
			R[i] = r[i];
		header = h;				//save the header size
		if(!open_file(filename)) return false;	//open the binary file
		//reset();
		return test_file_size();
	}
	/// Creates a new binary file for streaming
	/// @param filename is the name of the binary file to be created
	/// @param r is a STIM vector specifying the size of the file along each dimension
	/// @offset specifies how many bytes to offset the file (used to leave room for a header)
	bool create(std::string filename, vec<unsigned long long> r, unsigned long long offset = 0){
		std::ofstream target(filename.c_str(), std::ios::binary);
		//initialize binary file
		T p = 0;
		for(unsigned long long i =0; i < r[0] * r[1] * r[2]; i++){
			target.write((char*)(&p), sizeof(T));
		}
		for(unsigned long long i = 0; i < D; i++)		//set the dimensions of the binary file object
			R[i] = r[i];
		header = offset;				//save the header size
		if(!open_file(filename)) return false;	//open the binary file
		return test_file_size();
	}
	/// Writes a single page of data to disk. A page consists of a sequence of data of size R[0] * R[1] * ... * R[D-1].
	/// @param p is a pointer to the data to be written
	/// @param page is the page number (index of the highest-numbered dimension)
	bool write_page( T * p, unsigned long long page){
		if(p == NULL){
			std::cout<<"ERROR: unable to write into file, empty pointer"<<std::endl;
			exit(1);
		}
		file.seekg(R[1] * R[0] * page * sizeof(T) + header, std::ios::beg);   //seek to the desired location on disk
		file.write((char *)p, R[0] * R[1] * sizeof(T));				 //write binary data
		return true;
	}
	/// Reads a page from disk. A page consists of a sequence of data of size R[0] * R[1] * ... * R[D-1].
	/// @param p is a pointer to pre-allocated memory equal to the page size
	/// @param page is the index of the page
	bool read_page( T * p, unsigned long long page, bool PROGRESS = false){
		if(PROGRESS) progress = 0;
		if (page >= R[2]){										//make sure the bank number is right
			std::cout<<"ERROR: page out of range"<<std::endl;
			return false;
		}
		file.seekg(R[1] * R[0] * page * sizeof(T) + header, std::ios::beg);   //write into memory from the binary file
		file.read((char *)p, R[0] * R[1] * sizeof(T));
		if(PROGRESS) progress = 100;
		return true;
	}
	///Reads a line Z (slowest dimension) for a given XY value
	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
	/// @param x is the x coordinate
	/// @param y is the y coordinate
	void read_line_2( T* p, size_t n, bool PROGRESS = false){
		unsigned long long i;
		if(PROGRESS) progress = 0;
		if ( n > R[0] * R[1]){							//make sure the sample and line number is right
			std::cout<<"ERROR: sample or line out of range in "<<__FILE__<<" (line "<<__LINE__<<")"<<std::endl;
			exit(1);
		}
		file.seekg(n * sizeof(T), std::ios::beg);           //point to the certain sample and line
		for (i = 0; i < R[2]; i++){							//for each band
			file.read((char *)(p + i), sizeof(T));
			file.seekg((R[1] * R[0] - 1) * sizeof(T), std::ios::cur);    //go to the next band
			if(PROGRESS) progress = (double)i / (double)R[2] * 100;
		}
		if(PROGRESS) progress = 100;
	}
	void read_line_2( T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
		read_line_2(p, y * R[0] + x, PROGRESS);
		/*unsigned long long i;
		if(PROGRESS) progress = 0;
		if ( x >= R[0] || y >= R[1]){							//make sure the sample and line number is right
			std::cout<<"ERROR: sample or line out of range"<<std::endl;
			return false;
		}
		file.seekg((x + y * R[0]) * sizeof(T), std::ios::beg);           //point to the certain sample and line
		for (i = 0; i < R[2]; i++)
		{
			file.read((char *)(p + i), sizeof(T));
			file.seekg((R[1] * R[0] - 1) * sizeof(T), std::ios::cur);    //go to the next band
			if(PROGRESS) progress = (double)i / (double)R[2] * 100;
		}
		if(PROGRESS) progress = 100;
		return true;*/
	}
	///Reads a line X (fastest dimension) for a given YZ value
	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
	/// @param x is the y coordinate
	/// @param y is the z coordinate
	bool read_line_0(T * p, unsigned long long y, unsigned long long z, bool PROGRESS = false){
		//test to make sure the specified value is within range
		if( y >= R[1] || z >= R[2] ){
			std::cout<<"ERROR: sample ("<<y<<", "<<z<<") out of range in "<<__FILE__<<" (line "<<__LINE__<<")"<<std::endl;
			return false;
		}
		file.seekg((z * R[0] * R[1] + y * R[0]) * sizeof(T), std::ios::beg);	//seek to the start of the line
		file.read((char *)p, sizeof(T) * R[0]);									//read the line
		if(PROGRESS) progress = 100;
		return true;
	}
	///Reads a line Y (second fastest dimension) for a given XZ value
	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
	/// @param x is the y coordinate
	/// @param z is the z coordinate
	bool read_line_1(T * p, unsigned long long x, unsigned long long z, bool PROGRESS = false){
		if(PROGRESS) progress = 0;
		//test to make sure the specified value is within range
		if( x >= R[0] || z >= R[2] ){
			std::cout<<"ERROR: sample or line out of range in "<<__FILE__<<" (line "<<__LINE__<<")"<<std::endl;
			return false;
		}
		file.seekg((z * R[0] * R[1] + x) * sizeof(T), std::ios::beg);           //seek to the start of the line
		for (unsigned long long i = 0; i < R[1]; i++){									//for each pixel in the line
			file.read((char *)(p + i), sizeof(T));					//read the pixel
			file.seekg((R[0] - 1) * sizeof(T), std::ios::cur);		//seek to the next pixel in the line
			if(PROGRESS) progress = (double)i / (double)R[1] * 100;
		}
		if(PROGRESS) progress = 100;
		return true;
	}
	/// Reads a plane given a coordinate along the 0-axis (YZ plane)
	/// @param p is a pointer to pre-allocated memory of size R[1] * R[2] * sizeof(T)
	/// @param n is the 0-axis coordinate used to retrieve the plane
	bool read_plane_0(T* p, unsigned long long n, bool PROGRESS = false){
		if(PROGRESS) progress = 0;
		if (n >= R[0]){										//make sure the number is within the possible range
			std::cout<<"ERROR: sample or line out of range in "<<__FILE__<<" (line "<<__LINE__<<")"<<std::endl;
			return false;
		}
		unsigned long long jump = (R[0] - 1) * sizeof(T);		//number of bytes to skip between samples
		//seek to the start of the plane
		file.seekg(n * sizeof(T), std::ios::beg);
		unsigned long long N = R[1] * R[2];
		for(unsigned long long i = 0; i<N; i++){
			file.read((char*)(p+i), sizeof(T));
			file.seekg(jump, std::ios::cur);
			if(PROGRESS) progress = (double)(i+1) / N * 100;
		}
		return true;
	}
	/// Reads a plane given a coordinate along the 1-axis (XZ plane)
	/// @param p is a pointer to pre-allocated memory of size R[0] * R[2] * sizeof(T)
	/// @param n is the 1-axis coordinate used to retrieve the plane
	bool read_plane_1(T* p, unsigned long long n, bool PROGRESS = false){
		if(PROGRESS) progress = 0;
		unsigned long long L = R[0] * sizeof(T);		//caculate the number of bytes in a sample line
		unsigned long long jump = R[0] * (R[1] - 1) * sizeof(T);
		if (n >= R[1]){										//make sure the bank number is right
			std::cout<<"ERROR read_plane_1: page out of range"<<std::endl;
			return false;
		}
		file.seekg(R[0] * n * sizeof(T), std::ios::beg);
		for (unsigned long long i = 0; i < R[2]; i++){
			if(PROGRESS) progress = (double)i / R[2] * 100;
			file.read((char *)(p + i * R[0]), L);
			file.seekg( jump, std::ios::cur);
			std::cout<<i<<"    ";
		}
		if(PROGRESS) progress = 100;
		return true;
	}
	/// Reads a plane given a coordinate along the 2-axis (XY plane)
	/// @param p is a pointer to pre-allocated memory of size R[0] * R[1] * sizeof(T)
	/// @param n is the 2-axis coordinate used to retrieve the plane
	bool read_plane_2(T* p, unsigned long long n, bool PROGRESS = false){
		return read_page(p, n, PROGRESS);
	}
	/// Reads a single pixel, treating the entire data set as a linear array
	/// @param p is a pointer to pre-allocated memory of size sizeof(T)
	/// @param i is the index to the pixel using linear indexing
	bool read_pixel(T* p, unsigned long long i){
		if(i >= R[0] * R[1] * R[2]){
			std::cout<<"ERROR read_pixel: n is out of range"<<std::endl;
			return false;
		}
		file.seekg(i * sizeof(T), std::ios::cur);
		file.read((char*)p, sizeof(T));
	}
	/// Reads a single pixel, given an x, y, z coordinate
	/// @param p is a pointer to pre-allocated memory of size sizeof(T)
	/// @param x is the x (0) axis coordinate
	/// @param y is the y (1) axis coordinate
	/// @param z is the z (2) axis coordinate
	bool read_pixel(T* p, unsigned long long x, unsigned long long y, unsigned long long z){
		if(x < 0 || x >= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){
			std::cout<<"ERROR read_pixel: (x,y,z) is out of range"<<std::endl;
			return false;
		}
		unsigned long long i = z * R[0] * R[1] + y * R[0] + z;
		return read_pixel(p, i);
	}
	/// Reads a block specified by an (x, y, z) position and size using the largest possible contiguous reads
	size_t read(T* dest, size_t x, size_t y, size_t z, size_t sx, size_t sy, size_t sz){
		auto t0 = std::chrono::high_resolution_clock::now();
		size_t size_bytes = sx * sy * sz * sizeof(T);					//size of the block to read in bytes
		size_t start = z * R[0] * R[1] + y * R[0] + x;						//calculate the start postion
		size_t start_bytes = start * sizeof(T);							//start position in bytes
		file.seekg(start * sizeof(T), std::ios::beg);					//seek to the start position
		
		if(sx == R[0] && sy == R[1]){				//if sx and sy result in a contiguous volume along z
			file.read((char*)dest, size_bytes);			//read the block in one pass
		}
		else if(sx == R[0]){												//if sx is contiguous, read each z-axis slice can be read in one pass
			size_t jump_bytes = (R[1] - sy) * R[0] * sizeof(T);		//jump between each slice
			size_t slice_bytes = sx * sy * sizeof(T);				//size of the slice to be read
			for(size_t zi = 0; zi < sz; zi++){						//for each z-axis slice
				file.read((char*)dest, slice_bytes);						//read the slice
				dest += sx * sy;									//move the destination pointer to the next slice
				file.seekg(jump_bytes, std::ios::cur);				//skip to the next slice in the file
			}
		}
		else{
			//in this case, x is not contiguous so the volume must be read line-by-line
			size_t jump_x_bytes = (R[0] - sx) * sizeof(T);				//number of bytes skipped in the x direction
			size_t jump_y_bytes = (R[1] - sy) * R[0] * sizeof(T) + jump_x_bytes;	//number of bytes skipped between slices
			size_t line_bytes = sx * sizeof(T);							//size of the line to be read
			size_t zi, yi;
			for(zi = 0; zi < sz; zi++){									//for each slice
				file.read((char*)dest, line_bytes);							//read the first line
				for(yi = 1; yi < sy; yi++){								//read each additional line
					dest += sx;											//move the pointer in the destination block to the next line
					file.seekg(jump_x_bytes, std::ios::cur);			//skip to the next line in the file
					file.read((char*)dest, line_bytes);						//read the line to the destination block
				}
				file.seekg(jump_y_bytes, std::ios::cur);				//skip to the beginning of the next slice
			}
		}
		auto t1 = std::chrono::high_resolution_clock::now();
		return std::chrono::duration_cast<std::chrono::milliseconds>(t1-t0).count();
	}
	// permutes a block of data from the current interleave to the interleave specified (re-arranged dimensions to the order specified by [d0, d1, d2])
	size_t permute(T* dest, T* src, size_t sx, size_t sy, size_t sz, size_t d0, size_t d1, size_t d2){
		std::chrono::high_resolution_clock::time_point t0, t1;
		t0 = std::chrono::high_resolution_clock::now();
#ifdef CUDA_FOUND
		T* gpu_src;
		HANDLE_ERROR( cudaMalloc(&gpu_src, sx*sy*sz*sizeof(T)) );
		HANDLE_ERROR( cudaMemcpy(gpu_src, src, sx*sy*sz*sizeof(T), cudaMemcpyHostToDevice) );
		T* gpu_dest;
		HANDLE_ERROR( cudaMalloc(&gpu_dest, sx*sy*sz*sizeof(T)) );		
		gpu_permute((char*)gpu_dest, (char*)gpu_src, sx, sy, sz, d0, d1, d2, sizeof(T));
		HANDLE_ERROR( cudaMemcpy(dest, gpu_dest, sx*sy*sz*sizeof(T), cudaMemcpyDeviceToHost) );
		HANDLE_ERROR( cudaFree(gpu_src) );
		HANDLE_ERROR( cudaFree(gpu_dest) );
		t1 = std::chrono::high_resolution_clock::now();
		return std::chrono::duration_cast<std::chrono::milliseconds>(t1-t0).count();
		
#endif
		
		size_t d[3] = {d0, d1, d2};
		size_t s[3] = {sx, sy, sz};
		size_t p[3];
		
		if(d[0] == 0 && d[1] == 1 && d[2] == 2){
			//this isn't actually a permute - just copy the data
			memcpy(dest, src, sizeof(T) * sx * sy * sz);
		}
		else if(d[0] == 0){						//the individual lines are contiguous, so you can memcpy line-by-line
			size_t y, z;
			size_t src_idx, dest_idx;
			size_t x_bytes = sizeof(T) * sx;
			for(z = 0; z < sz; z++){
				p[2] = z;
				for(y = 0; y < sy; y++){
					p[1] = y;
					src_idx = z * sx * sy + y * sx;
					dest_idx = p[d[2]] * s[d[0]] * s[d[1]] + p[d[1]] * s[d[0]];
					memcpy(dest + dest_idx, src + src_idx, x_bytes);
				}
			}
		}
		else{									//loop through every damn point
			size_t x, y, z;
			size_t src_idx, dest_idx;
			size_t src_z, src_y;
			for(z = 0; z < sz; z++){
				p[2] = z;
				src_z = z * sx * sy;
				for(y = 0; y < sy; y++){
					p[1] = y;
					src_y = src_z + y * sx;
					for(x = 0; x < sx; x++){
						p[0] = x;
						src_idx = src_y + x;
						dest_idx = p[d[2]] * s[d[0]] * s[d[1]] + p[d[1]] * s[d[0]] + p[d[0]];
						dest[dest_idx] = src[src_idx];
					}
				}
			}
		}
		t1 = std::chrono::high_resolution_clock::now();
		return std::chrono::duration_cast<std::chrono::milliseconds>(t1-t0).count();
	}
};
}
#endif