twirls/MexFunc/AllClusterRandSim.cpp

#include <mex.h>
#include <mat.h>
#include <iostream>
#include <algorithm>
#include <string>
#include <unordered_set>
#include <ctime>
#include <vector>
#include <queue>
#include <memory>
#include <thread>
#include <mutex>
#include <condition_variable>
#include <future>
#include <functional>
#include <stdexcept>
#include <unordered_map>
#include <set>
#include <fstream>
#include <random>
#include <cmath>
#include <stdlib.h>
#include <limits.h>
#include <atomic>

using std::cout;
using std::endl;
using namespace std;

#define STRING_BUF_SIZE 204800

class ThreadPool {
public:
	ThreadPool(size_t);
	template<class F, class... Args>
	auto enqueue(F&& f, Args&&... args)
		->std::future<typename std::result_of<F(Args...)>::type>;
	~ThreadPool();
private:
	// need to keep track of threads so we can join them
	std::vector< std::thread > workers;
	// the task queue
	std::queue< std::function<void()> > tasks;

	// synchronization
	std::mutex queue_mutex;
	std::condition_variable condition;
	bool stop;
};

// the constructor just launches some amount of workers
inline ThreadPool::ThreadPool(size_t threads)
	: stop(false)
{
	for (size_t i = 0;i < threads;++i)
		workers.emplace_back(
			[this]
			{
				for (;;)
				{
					std::function<void()> task;

					{
						std::unique_lock<std::mutex> lock(this->queue_mutex);
						this->condition.wait(lock,
							[this] { return this->stop || !this->tasks.empty(); });
						if (this->stop && this->tasks.empty())
							return;
						task = std::move(this->tasks.front());
						this->tasks.pop();
					}

					task();
				}
			}
			);
}

// add new work item to the pool
template<class F, class... Args>
auto ThreadPool::enqueue(F&& f, Args&&... args)
-> std::future<typename std::result_of<F(Args...)>::type>
{
	using return_type = typename std::result_of<F(Args...)>::type;

	auto task = std::make_shared< std::packaged_task<return_type()> >(
		std::bind(std::forward<F>(f), std::forward<Args>(args)...)
		);

	std::future<return_type> res = task->get_future();
	{
		std::unique_lock<std::mutex> lock(queue_mutex);

		// don't allow enqueueing after stopping the pool
		if (stop)
			throw std::runtime_error("enqueue on stopped ThreadPool");

		tasks.emplace([task]() { (*task)(); });
	}
	condition.notify_one();
	return res;
}

// the destructor joins all threads
inline ThreadPool::~ThreadPool()
{
	{
		std::unique_lock<std::mutex> lock(queue_mutex);
		stop = true;
	}
	condition.notify_all();
	for (std::thread& worker : workers)
		worker.join();
}

/* <20><>ȡһάdouble<6C><65><EFBFBD><EFBFBD> */
void Read1DDouble(const  mxArray* pMxArray, vector<double>& vDat) {
	int rowNum, colNum;
	double* matData;

	rowNum = (int)mxGetM(pMxArray);
	colNum = (int)mxGetN(pMxArray);
	// cout << rowNum << " " << colNum << endl;
	matData = (double*)mxGetData(pMxArray); //<2F><>ȡָ<C8A1><D6B8>
	vDat.resize(rowNum * colNum);
	for (int i = 0; i < vDat.size(); ++i) vDat[i] = matData[i];
}

/* <20><>ȡ<EFBFBD><C8A1>άdouble<6C><65><EFBFBD><EFBFBD> */
void Read2DDouble(const  mxArray* pMxArray, vector<vector<double>>& vvDat) {
	int rowNum, colNum;
	double* matData;

	rowNum = (int)mxGetM(pMxArray);
	colNum = (int)mxGetN(pMxArray);
	vvDat.resize(rowNum);
	matData = (double*)mxGetData(pMxArray); //<2F><>ȡָ<C8A1><D6B8>
	for (int i = 0; i < rowNum; ++i) {
		vvDat[i].resize(colNum);
		for (int j = 0; j < colNum; ++j) {
			vvDat[i][j] = matData[j * rowNum + i];
		}
	}
}

// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>д<EFBFBD><D0B4>mxArray, <20><>Ϊ<EFBFBD><CEAA><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD><C4B7><EFBFBD>ֵ
mxArray* writeToMatDouble(const double* data, int rowNum, int colNum) {
	mxArray* pWriteArray = NULL;//matlab<61><62>ʽ<EFBFBD><CABD><EFBFBD><EFBFBD>
	int len = rowNum * colNum;
	//<2F><><EFBFBD><EFBFBD>һ<EFBFBD><D2BB>rowNum*colNum<75>ľ<EFBFBD><C4BE><EFBFBD>  
	pWriteArray = mxCreateDoubleMatrix(rowNum, colNum, mxREAL);
	//<2F><>data<74><61>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD>pWriteArrayָ<79><D6B8>
	memcpy((void*)(mxGetPr(pWriteArray)), (void*)data, sizeof(double) * len);
	return pWriteArray; // <20><>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
}

#define START_IDX(tid, threadNum, arrLen) (arrLen) * (tid) / (threadNum)
#define END_IDX(tid, threadNum, arrLen) (arrLen) * (tid + 1) / (threadNum)

// <20>̲߳<DFB3><CCB2><EFBFBD>
struct TPRandSim {
	vector<vector<double>>* pvvSum;
	vector<vector<double>>* pvvSqSum;
	vector<vector<int>>* pvvRandPos;
	vector<double>* pvIx;
	vector<int>* pvCateNum;
	vector<vector<double>>* pvvRound;
	vector<vector<double>>* pvvX;
	int loopNum;
	int maxCategoryNum;
	int tid;
	int numThread;
};

// <20><><EFBFBD>߳<EFBFBD><DFB3><EFBFBD><EFBFBD>ں<EFBFBD><DABA><EFBFBD>
void ThreadRandSim(TPRandSim param) {
	vector<vector<double>>& vvSum = *param.pvvSum;
	vector<vector<double>>& vvSqSum = *param.pvvSqSum;
	vector<vector<int>>& vvRandPos = *param.pvvRandPos;
	vector<double>& vIx = *param.pvIx;
	vector<int>& vCateNum = *param.pvCateNum;
	vector<vector<double>>& vvRound = *param.pvvRound;
	vector<vector<double>>& vvX = *param.pvvX;
	int tid = param.tid;
	int numThread = param.numThread;
	int loopNum = param.loopNum;
	int maxCategoryNum = param.maxCategoryNum;
	int rowNum = vvX.size();
	int colNum = vvX[0].size();

	clock_t begin = clock(), mid = clock(), finish;
	/* <20><><EFBFBD><EFBFBD>ģ<EFBFBD><C4A3> */
	int startIdx = START_IDX(tid, numThread, loopNum);
	int endIdx = END_IDX(tid, numThread, loopNum);

	// if (tid == 0) cout << startIdx << '\t' << endIdx << endl;

	for (int idx = startIdx; idx < endIdx; ++idx) { // ģ<><C4A3><EFBFBD>ִ<EFBFBD>
		// if (tid == 0 && idx % 100 == 0) {
		// 	finish = clock();
		// 	cout << idx << ": time: " << (double)(finish - mid) / CLOCKS_PER_SEC << " s" << endl << flush;
		// 	mid = finish;
		// }
		auto& vRandPos = vvRandPos[idx];
		for (int i = 0; i < rowNum; ++i) {
			const int hRowIdx = vRandPos[i]; // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>֮<EFBFBD><D6AE><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
			const int cateIdx = vIx[hRowIdx] - 1; // <20><><EFBFBD><EFBFBD><EFBFBD>ı<EFBFBD><C4B1><EFBFBD>
			//if (cateIdx == 0) {
				auto& vRound = vvRound[cateIdx];
				for (int j = 0; j < colNum; ++j) {
					vRound[j] += vvX[i][j];
				}
			//}
		}
		for (int c = 0; c < maxCategoryNum; ++c) {
			auto& vRound = vvRound[c];
			const int numPositive = vCateNum[c];
			for (int j = 0; j < colNum; ++j) {
				const double val = vRound[j] / numPositive;
				vvSum[c][j] += val;
				vvSqSum[c][j] += val * val;
			}
		}
		for (auto& vRound : vvRound) for (auto &val : vRound) val = 0;
	}
	finish = clock();
	// cout << tid << ": Random simulation time: " << (double)(finish - begin) / CLOCKS_PER_SEC << " s" << endl << flush;
}

/* <20><><EFBFBD>߳̽<DFB3><CCBD><EFBFBD>shuffle<6C><65><EFBFBD><EFBFBD> */
struct TPShuffle {
	vector<vector<int>>* pvvRandPos;
	int loopNum;
	int tid;
	int numThread;
};
void ThreadShuffle(TPShuffle param) {
	clock_t begin = clock(), finish;
	vector<vector<int>>& vvRandPos = *param.pvvRandPos;
	int tid = param.tid;
	int numThread = param.numThread;
	int loopNum = param.loopNum;
	std::random_device rd;
	int startIdx = START_IDX(tid, numThread, loopNum);
	int endIdx = END_IDX(tid, numThread, loopNum);
	for (int roundIdx = startIdx; roundIdx < endIdx; ++roundIdx) {
		vector<int>& vRandPos = vvRandPos[roundIdx];
		for (int i = 0; i < vRandPos.size(); ++i) vRandPos[i] = i;
		std::shuffle(vRandPos.begin(), vRandPos.end(), std::default_random_engine(rd()));
	}
	finish = clock();
	// cout << tid << ": thread shuffle time: " << (double)(finish - begin) / CLOCKS_PER_SEC << " s" << endl << flush;
}

/* <20><><EFBFBD>ں<EFBFBD><DABA><EFBFBD> */
/*
<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
<EFBFBD><EFBFBD><EFBFBD>룺
1. x <EFBFBD><EFBFBD>ά<EFBFBD><EFBFBD><EFBFBD>ݣ<EFBFBD>double<EFBFBD><EFBFBD><EFBFBD>ͣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD>ֵ䳤<EFBFBD>ȣ<EFBFBD>ÿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>г<EFBFBD><EFBFBD>ֵĴ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>5<EFBFBD><EFBFBD>
2. h <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD>׸<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵΪ1<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڸ<EFBFBD>֪ʶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD>ǣ<EFBFBD><EFBFBD><EFBFBD>Ϊ0<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
3. numThread
4. loopNum
<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
vs z score,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>һά
ps <EFBFBD><EFBFBD>vs<EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD>
*/
void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
	if (nrhs < 2) {
		cout << "At least 2 arguments should be given for this function!" << endl;
		return;
	}
	clock_t begin = clock(), mid, finish;
	vector<double> vIx; // <20><><EFBFBD><EFBFBD>
	vector<vector<double>> vvX;


	Read2DDouble(prhs[0], vvX);
	Read1DDouble(prhs[1], vIx);
	int rowNum = vvX.size();
	int colNum = vvX[0].size();

	// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
	int maxCategoryNum = 1;
	for (auto category : vIx) if (maxCategoryNum < category) maxCategoryNum = category;

	int numThread = 1;
	int loopNum = 1000;

	if (nrhs > 2) {
		double* pNumThread = (double*)mxGetData(prhs[2]);
		numThread = (int)pNumThread[0];
		if (numThread < 1) numThread = 1;
	}
	if (nrhs > 3) {
		double* pLoopNum = (double*)mxGetData(prhs[3]);
		loopNum = (int)pLoopNum[0];
		if (loopNum < 1000) loopNum = 1000;
	}
	// <20>߳<EFBFBD><DFB3><EFBFBD><EFBFBD>õ<EFBFBD><C3B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
	vector<int> vCateNum(maxCategoryNum); // ÿ<><C3BF><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><C4B8><EFBFBD>
	vector<vector<double>> vvTs(maxCategoryNum, vector<double>(colNum)); //<2F><>¼<EFBFBD><C2BC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>е<EFBFBD><D0B5>о<EFBFBD>ֵ
	vector<vector<vector<double>>> vvvSum(numThread, vector<vector<double>>(maxCategoryNum, vector<double>(colNum)));
	vector<vector<vector<double>>> vvvSqSum(numThread, vector<vector<double>>(maxCategoryNum, vector<double>(colNum)));
	vector<vector<vector<double>>> vvvRound(numThread, vector<vector<double>>(maxCategoryNum, vector<double>(colNum)));

	finish = clock();
	cout << "Load Data time: " << (double)(finish - begin) / CLOCKS_PER_SEC << " s" << endl << flush;
	// cout << numThread << '\t' << loopNum << endl;

	// maxCategoryNum = 1;

	/* <20><><EFBFBD><EFBFBD>ʵ<EFBFBD>ʷ<EFBFBD><CAB7><EFBFBD><EFBFBD>ľ<EFBFBD>ֵ */
	mid = clock();
	for (int i = 0; i < rowNum; ++i) {
		const int cateIdx = vIx[i] - 1;
		vCateNum[cateIdx] ++;
		for (int j = 0; j < colNum; ++j) vvTs[cateIdx][j] += vvX[i][j];
	}
	for (int i = 0; i < maxCategoryNum; ++i) {
		for (int j = 0; j < colNum; ++j) {
			vvTs[i][j] /= vCateNum[i];
		}
	}

	// for (auto c : vCateNum) cout << c << endl;

	// <20>Ȱ<EFBFBD>loopNum<75><6D><EFBFBD><EFBFBD><EFBFBD><EFBFBD>shuffle<6C><65><EFBFBD><EFBFBD>
	mid = clock();
	vector<vector<int>> vvRandPos(loopNum, vector<int>(rowNum));
	vector<std::thread> vTHShuffle;
	for (int i = 0; i < numThread; ++i) {
		TPShuffle param = { &vvRandPos, loopNum, i, numThread };
		vTHShuffle.push_back(std::thread(ThreadShuffle, param));
	}
	for (auto& t : vTHShuffle) t.join();
	finish = clock();
	cout << "Shuffle time: " << (double)(finish - begin) / CLOCKS_PER_SEC << " s" << endl << flush;

	vector<std::thread> vTHRandSim;
	for (int i = 0; i < numThread; ++i) {
		TPRandSim tParam = { &vvvSum[i], &vvvSqSum[i], &vvRandPos, &vIx, &vCateNum, &vvvRound[i], & vvX, loopNum, maxCategoryNum, i, numThread};
		vTHRandSim.push_back(std::thread(ThreadRandSim, tParam));
	}
	for (auto& t : vTHRandSim) t.join();

	// <20>ϲ<EFBFBD><CFB2><EFBFBD><EFBFBD><EFBFBD><EFBFBD>̵߳<DFB3><CCB5><EFBFBD><EFBFBD><EFBFBD>
	auto& vvSum = vvvSum[0];
	auto& vvSqSum = vvvSqSum[0];
	for (int t = 1; t < numThread; ++t) {
		for (int i = 0; i < maxCategoryNum; ++i) {
			for (int j = 0; j < colNum; ++j) {
				vvSum[i][j] += vvvSum[t][i][j];
				vvSqSum[i][j] += vvvSqSum[t][i][j];
			}
		}
	}

	finish = clock();
	cout << "Random simulation time: " << (double)(finish - mid) / CLOCKS_PER_SEC << " s" << endl << flush;

	/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> */
	vector<vector<double>> vvVs(maxCategoryNum, vector<double>(colNum));

	// <20><><EFBFBD>м<EFBFBD><D0BC><EFBFBD>ƽ<EFBFBD><C6BD>ֵ
	vector<vector<double>> vvMean(maxCategoryNum, vector<double>(colNum));
	vector<vector<double>> vvStd(maxCategoryNum, vector<double>(colNum));

	for (int c = 0; c < maxCategoryNum; ++c) {
		auto& vMean = vvMean[c];
		auto& vStd = vvStd[c];
		auto& vVs = vvVs[c];
		auto& vSum = vvSum[c];
		auto& vSqSum = vvSqSum[c];
		auto& vTs = vvTs[c];

		for (int i = 0; i < colNum; ++i) vMean[i] = vSum[i] / loopNum;
		for (int i = 0; i < colNum; ++i) {
			const double meanVal = vSum[i] / loopNum;
			vMean[i] = meanVal; // <20><>ֵ
			const double sqDiff = vSqSum[i] + loopNum * meanVal * meanVal - 2 * meanVal * vSum[i];
			vStd[i] = sqrt(sqDiff / (loopNum - 1)); // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
		}

		// <20><><EFBFBD><EFBFBD>vs
		for (int i = 0; i < vVs.size(); ++i) {
			vVs[i] = (vTs[i] - vMean[i]) / vStd[i];
		}
	}

	// auto& vVs = vvVs[0];
	// ofstream ofs("d:\\result_new.txt");
	// for (int i = 0; i < colNum; ++i) {
	// 	ofs << vVs[i] << endl;
	// }
	// ofs.close();

	/* д<><D0B4><EFBFBD><EFBFBD><EFBFBD><EFBFBD> */
	if (nlhs > 0) { // vs
		vector<double> vVsData(maxCategoryNum* colNum);
		for (int i = 0; i < maxCategoryNum; ++i) {
			for (int j = 0; j < colNum; ++j) {
				vVsData[j * maxCategoryNum + i] = vvVs[i][j];
			}
		}
		plhs[0] = writeToMatDouble(vVsData.data(), maxCategoryNum, colNum);
	}

	finish = clock();
	cout << "All Cluster Random simulation Total time: " << (double)(finish - begin) / CLOCKS_PER_SEC << " s" << endl << flush;
}

// <20><>c++<2B><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
void mexFunctionWrap(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
	return mexFunction(nlhs, plhs, nrhs, prhs);
}