chapter20.cpp

#include "chapter20.h"

#include "stdafx.h"
#include "Executor.h"
#include "testing.h"
#include "geometry.h"
#include "matlib.h"

using namespace std;

double multiThreadSum(const std::vector<double>& vec, int nThreads) {

	auto executor = Executor::newInstance();

	int batchSize = vec.size() / nThreads;
	vector<double> result(nThreads);
	for (int i = 0; i < nThreads; i++) {
		int startIndex = i*batchSize;
		int endIndex = (i + 1)*batchSize;
		if (i == nThreads - 1) {
			endIndex = vec.size();
		}
		if (endIndex > (int)vec.size()) {
			endIndex = vec.size();
		}
		executor->executeTask([&result,i,&vec,startIndex,endIndex]() {
			double sum = 0.0;
			for (int inner = startIndex; inner < endIndex; inner++) {
				sum += vec[inner];
			}
			result[i] = sum;
		});
	}
	executor->join();

	double sum = 0.0;
	for (int i = 0; i < nThreads; i++) {
		sum += result[i];
	}
	return sum;
}


static void testMultiThreadedSum() {
	vector<double> d(100);
	for (int i = 0; i < 100; i++) {
		d[i] = i + 1;
	}
	ASSERT_APPROX_EQUAL( multiThreadSum(d,7), 5050.0, 0.0001);
}

double randuniform(mt19937& random) {
	return (random() + 0.5) / (random.max() + 1.0);
}

/**
 *   Integrate a 2d function
 */
static double integrate2D(function<double(double, double)> f,
						const Rectangle& rectangle,
						mt19937& random,
						int nPoints) {
	double total = 0.0;
	for (int i = 0; i < nPoints; i++) {
		double u1 = randuniform(random);
		double u2 = randuniform(random);
		double x = rectangle.getLeft() + rectangle.getWidth() * u1;
		double y = rectangle.getBottom() + rectangle.getHeight() * u2;
		total += f(x, y);
	}
	return total * rectangle.getWidth()* rectangle.getHeight() / nPoints;
}

static double integrate2DMultiThreaded(function<double(double, double)> f,
				const Rectangle& rectangle,
				mt19937& random,
				int nPoints,
				int nThreads) {
	auto executor = Executor::newInstance();

	int batchSize = nPoints / nThreads;
	vector<double> result(nThreads);
	for (int i = 0; i < nThreads; i++) {
		int startIndex = i*batchSize;
		int endIndex = (i + 1)*batchSize;
		if (i == nThreads - 1) {
			endIndex = nPoints;
		}
		if (endIndex >nPoints) {
			endIndex = nPoints;
		}
		if (startIndex != endIndex) {
			executor->executeTask([&f, &result, i, &random, rectangle, startIndex, endIndex, nPoints]() {
				mt19937 randomCopy = random;
				randomCopy.discard(2*startIndex);
				result[i] = (integrate2D(f, rectangle, randomCopy, endIndex - startIndex) * (endIndex-startIndex)) / (nPoints);
			});
		}
	}
	executor->join();
	random.discard(nPoints);

	double sum = 0.0;
	for (int i = 0; i < nThreads; i++) {
		sum += result[i];
	}
	return sum;
}

static void testIntegrate2D() {
	Rectangle rect(6.0, 2.0, 1.0, 3.0);
	auto f = [](double x, double y) {
		return sin(x) + cos(y);
	};
	mt19937 random1;
	double integral = integrate2D(f, rect,random1,10000 );
	mt19937 random2;

	double integral2 = integrate2DMultiThreaded(f, rect, random2, 10000, 7);
	ASSERT_APPROX_EQUAL(integral, integral2, 0.00001);
	double expected = 2*((-1 + 4 * sin(1))*sin(2) + sin(6));
	ASSERT_APPROX_EQUAL(integral, expected, 0.02);
}







void testChapter20() {
	TEST(testMultiThreadedSum);
	TEST(testIntegrate2D);

}