twirls/GMM/main.cpp

/*********************************************************************************************
	Description: <EFBFBD>ø<EFBFBD>˹<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD>ͽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ߵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>

	Copyright : All right reserved by ZheYuan.BJ

	Author : Zhang Zhonghai
	Date : 2023/09/12
***********************************************************************************************/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <random>
#include <unordered_map>
#include <omp.h>
#include <time.h>
#include <string>
#include <vector>
#include <queue>
#include <filesystem>
#include <thread>
#ifdef _WIN32
#include <io.h>
#include <process.h>
#define F_OK 0
#else
#include <unistd.h>
#endif
#include <mat.h>
#include "gmm.h"
// #include "CommonLib/thread_pool.h"
#include "CommonLib/matlab_io.h"
#include "CommonLib/kthread.h"
using namespace std;
using std::cout;
using std::vector;
namespace fs = std::filesystem;

/* <20><>x<EFBFBD><78><EFBFBD><EFBFBD><EFBFBD>ŵ<EFBFBD><C5B5><EFBFBD><EFBFBD><EFBFBD>ΪbinWidth<74><68>С<EFBFBD><D0A1>Ͱ<EFBFBD>У<EFBFBD><D0A3><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>matlab<61><62>hist*/
void PutXtoBin(double* x, int xSize, double binWidth, vector<double>& vXBin, vector<double>& vYBin) {
	double maxX = 0.0;
	for (int i = 0; i < xSize; ++i) {
		if (maxX < x[i]) maxX = x[i];
	}
	int binSize = (int)((maxX + binWidth / 2) / binWidth + 1);
	double binMaxVal = (binSize - 1) * binWidth;
	if (binMaxVal > maxX) { // ȷ<><C8B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><D2BB>bin<69><6E><EFBFBD><EFBFBD><EFBFBD><EFBFBD>maxX<78><58><EFBFBD><EFBFBD><EFBFBD>Ҳ<EFBFBD>С<EFBFBD><D0A1>maxX-binWidth
		binSize -= 1;
	}
	vXBin.resize(xSize);
	vYBin.resize(binSize);
	for (int i = 0; i < binSize; ++i) vYBin[i] = 0;
	for (int i = 0; i < xSize; ++i) {
		int binIdx = (int)((x[i] + binWidth / 2) / binWidth);
		if (binIdx >= binSize) binIdx = binSize - 1;
		vYBin[binIdx] += 1;
	}
	// <20><><EFBFBD><EFBFBD>С˳<D0A1><CBB3><EFBFBD><EFBFBD><EFBFBD>޸ĺ<DEB8><C4BA><EFBFBD>x<EFBFBD><78>ֵ<EFBFBD>洢<EFBFBD><E6B4A2>vXBin<69>У<EFBFBD><D0A3><EFBFBD><EFBFBD><EFBFBD>˳<EFBFBD><CBB3><EFBFBD><EFBFBD>ͬ<EFBFBD><CDAC>ѵ<EFBFBD><D1B5><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD>˹<EFBFBD><CBB9><EFBFBD><EFBFBD>ģ<EFBFBD>Ͳ<EFBFBD><CDB2><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һЩ<D2BB><D0A9>ͬ<EFBFBD><CDAC>
	int xIdx = 0;
	for (int i = 0; i < binSize; ++i) {
		for (int j = 0; j < vYBin[i]; ++j) {
			vXBin[xIdx++] = i * binWidth;
		}
	}
}

/* <20><><EFBFBD><EFBFBD>׼<EFBFBD><D7BC>˹ģ<CBB9><C4A3>ѵ<EFBFBD><D1B5><EFBFBD><EFBFBD><EFBFBD>IJ<EFBFBD><C4B2><EFBFBD>ת<EFBFBD><D7AA><EFBFBD><EFBFBD><EFBFBD>Զ<EFBFBD><D4B6><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5>, <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϻ<EFBFBD><CFBA><EFBFBD>Yֵ<59><D6B5><EFBFBD><EFBFBD> */
struct cmpFunc {
	bool operator()(const pair<double, double>& a, const pair<double, double>& b) { return a.first < b.first; }
};
void GMMToFactorEY(GMM& gmm, double binWidth, vector<double> &vYBin, vector<double>& vFactor, vector<double>& vEY) {
	/* <20><>Ҫ<EFBFBD><D2AA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ߵ<EFBFBD>Ȩ<EFBFBD>أ<EFBFBD><D8A3><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϸ<EFBFBD>˹<EFBFBD><CBB9><EFBFBD>ߣ<EFBFBD><DFA3><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ø<EFBFBD><C3B8><EFBFBD><EFBFBD>ܶ<EFBFBD> */
	double zoomFactorSum = 0.0;
	vEY.resize(vYBin.size());
	int topM = (int)(vYBin.size() / 4);
	if (topM < 1) topM = 1;
	/* <20>ö<EFBFBD><C3B6><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽȡǰtopM<70><4D><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ, <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ų<EFBFBD><C5B2><EFBFBD>*/
	priority_queue<pair<double, double>, vector<pair<double, double> >, cmpFunc> pqTopM;

	for (int i = 0; i < vYBin.size(); ++i) {
		double xVal = i * binWidth;
		double probVal = gmm.GetProbability(&xVal);
		vEY[i] = probVal;
		pqTopM.push(make_pair(vYBin[i], probVal));
	}

	for (int i = 0; i < topM; ++i) {
		pair<double, double> topEle = pqTopM.top();
		pqTopM.pop();
		zoomFactorSum += topEle.first / topEle.second;
	}

	double zoomFactor = zoomFactorSum / topM;

	for (int i = 0; i < vEY.size(); ++i) {
		vEY[i] *= zoomFactor;
	}
	vFactor.clear();
	vFactor.push_back(zoomFactor * gmm.Prior(0) / sqrt(2 * M_PI * *gmm.Variance(0)));
	vFactor.push_back(*gmm.Mean(0));
	vFactor.push_back(sqrt(2 * *gmm.Variance(0)));
	vFactor.push_back(zoomFactor * gmm.Prior(1) / sqrt(2 * M_PI * *gmm.Variance(1)));
	vFactor.push_back(*gmm.Mean(1));
	vFactor.push_back(sqrt(2 * *gmm.Variance(1)));
}

/* <20><><EFBFBD><EFBFBD>ƽ<EFBFBD><C6BD><EFBFBD><EFBFBD> */
template <typename T>
T Average(vector<T>& vVal) {
	T sumVal = T(0);
	for (int i = 0; i < vVal.size(); ++i) {
		sumVal += vVal[i];
	}
	return sumVal / vVal.size();
}

/* <20><><EFBFBD><EFBFBD>ƽ<EFBFBD><C6BD><EFBFBD>ľ<EFBFBD>ֵ */
template <typename T>
T SquareAverage(vector<T>& vVal) {
	vector<T> vSquare(vVal.size());
	for (int i = 0; i < vVal.size(); ++i) {
		vSquare[i] = vVal[i] * vVal[i];
	}
	return Average(vSquare);
}

/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>x<EFBFBD><78>y<EFBFBD><79><EFBFBD><EFBFBD><EFBFBD>ؾ<EFBFBD><D8BE><EFBFBD>, <20><><EFBFBD><EFBFBD>ά<EFBFBD>ȱ<EFBFBD><C8B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>*/
double CorrelationDistance(vector<double>& vX, vector<double>& vY) {
	double xMean = Average(vX);
	transform(vX.cbegin(), vX.cend(), vX.begin(), [&](double val) { return val - xMean; });
	double yMean = Average(vY);
	transform(vY.cbegin(), vY.cend(), vY.begin(), [&](double val) { return val - yMean; });

	vector<double> vXY(vX.size());
	for (int i = 0; i < vXY.size(); ++i) {
		vXY[i] = vX[i] * vY[i];
	}
	double uv = Average(vXY);
	double uu = SquareAverage(vX);
	double vv = SquareAverage(vY);

	double dist = 1.0 - uv / sqrt(uu * vv);
	return abs(dist);
}

/* <20><><EFBFBD><EFBFBD>һ<EFBFBD><D2BB>֪ʶ<D6AA><CAB6><EFBFBD><EFBFBD> */
struct ThreadParamKP {
	fs::path matFilePath;
	fs::path outFilePath;
};
//void ThreadProcessData(vector<ThreadParam>& vTP, long idx, int tid) {
void ThreadProcessKP(ThreadParamKP& param) {
	//const ThreadParam& param = vTP[idx];
	const fs::path& matFilePath = param.matFilePath;
	const fs::path& outFilePath = param.outFilePath;
	double* hs = nullptr;
	int rowNum = 0;
	int colNum = 0;

	hs = ReadMtxDouble(matFilePath.string(), "hs", &rowNum, &colNum);

	vector<double>vXBin;
	vector<double>vYBin;
	vector<double>vEY;
	vector<double>vFactor;
	/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݣ<EFBFBD><DDA3><EFBFBD><EFBFBD><EFBFBD>mat<61>ļ<EFBFBD> */
	vector<double>vDist(rowNum);
	vector<double>vFactorAll;
	for (int i = 0; i < rowNum; ++i) {
		PutXtoBin(hs + i * colNum, colNum, 0.2, vXBin, vYBin);

		GMM gmm(1, 2); // 1ά<31><CEAC> 2<><32><EFBFBD><EFBFBD>˹ģ<CBB9><C4A3>
		gmm.Train(vXBin.data(), (int)vXBin.size());

		GMMToFactorEY(gmm, 0.2, vYBin, vFactor, vEY);
		vDist[i] = CorrelationDistance(vYBin, vEY);
		vFactorAll.insert(vFactorAll.end(), vFactor.begin(), vFactor.end());
	}
	/* д<><D0B4>matlab<61>ļ<EFBFBD> */
	MATFile* pMatFile = matOpen(outFilePath.string().c_str(), "w");
	SaveMtxDouble(vFactorAll.data(), pMatFile, "factor", rowNum, 6);
	SaveMtxDouble(vDist.data(), pMatFile, "correlation", rowNum, 1);
	matClose(pMatFile);

	delete[] hs;
}

/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> */
int main(int argc, const char** argv) {

	if (argc < 4) {
		cout << "This program take at least 3 arguments(1.parrent Dir; 2. mat file suffix; 3. out mat filename; [4. thread number])!" << endl;
		return 1;
	}
	string parrentDir(argv[1]); // ֪ʶ<D6AA><CAB6><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD>Ŀ¼<C4BF><C2BC><EFBFBD><EFBFBD>
	string hsMatSuffix(argv[2]); // hs<68><73><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><D3A6>mat<61>ļ<EFBFBD><C4BC>ĺ<EFBFBD>׺<EFBFBD><D7BA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȫ<EFBFBD>ļ<EFBFBD><C4BC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><C4BC><EFBFBD><EFBFBD><EFBFBD>׺<EFBFBD><D7BA><EFBFBD><EFBFBD><EFBFBD>뱣֤Ψһ<CEA8><D2BB>
	fs::path outFileName(argv[3]);
	int numThread = 1;
	if (argc >= 4) numThread = atoi(argv[4]);
	if (numThread < 1) numThread = 1;
	//ThreadPool thPool(numThread);
	clock_t begin, finish;
	begin = clock();
	vector<ThreadParamKP> vTP;
	/* <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>е<EFBFBD>֪ʶ<D6AA><CAB6><EFBFBD><EFBFBD>Ŀ¼<C4BF><C2BC><EFBFBD><EFBFBD>һ<EFBFBD><D2BB><EFBFBD>д<EFBFBD><D0B4><EFBFBD> */
	for (auto& childDir : fs::directory_iterator(parrentDir)) {
		fs::path outFilePath = childDir / outFileName;
		for (auto& file : fs::directory_iterator(childDir)) {
			const string& fileName = file.path().filename().string();
			auto rPos = fileName.rfind(hsMatSuffix);
			if (rPos != string::npos && fileName.size() - rPos == hsMatSuffix.size()) {
				//ThreadParam tParam = { file, outFilePath };
				//thPool.enqueue(ThreadProcessData, tParam);
				vTP.push_back({ file, outFilePath });
			}
		}
	}
	kt_for(numThread, ThreadProcessKP, vTP);
	//thPool.~ThreadPool();
	finish = clock();
	cout << "GMM Total time: " << (double)(finish - begin) / CLOCKS_PER_SEC << " s" << endl;
	return 0;
}