text:
/* Copyright (c) 2015 USC, IRIS, Computer vision Lab */
#include "FaceServices2.h"
#include <fstream>
//#include "opencv2/contrib/contrib.hpp"
#include <opencv2/calib3d.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <Eigen/SparseLU>
//#include <Eigen/SPQRSupport>
#include <omp.h>

using namespace std;
using namespace cv;

FaceServices2::FaceServices2(void)
{
	//omp_set_num_threads(8);
	prevEF = 1000;
	mstep = 0.0001;
	countFail = 0;
	maxVal = 4;
	mlambda = 0.005;
	PREV_USE_THRESH = 3.141592/9;
}

void FaceServices2::setUp(int w, int h, float f){
	memset(_k,0,9*sizeof(float));
	_k[8] = 1;
	_k[0] = -f;
	_k[4] = f;
	_k[2] = w/2.0f;
	_k[5] = h/2.0f;
	faces = festimator.getFaces() - 1;
	cv::Mat shape = festimator.getShape(cv::Mat(99,1,CV_32F));
	tex = shape*0 + 128;
	im_render = new FImRenderer(cv::Mat::zeros(h,w,CV_8UC3));
	im_render->loadMesh(shape,shape*0,faces);
}

// Yuval
void FaceServices2::init(int w, int h, float f)
{
    // Initialize camera matrix
    memset(_k, 0, 9 * sizeof(float));
    _k[8] = 1;
    _k[0] = -f;
    _k[4] = f;
    _k[2] = w / 2.0f;
    _k[5] = h / 2.0f;

    // Initialize shape and texture
    if (faces.empty())
        faces = festimator.getFaces() - 1;
    cv::Mat shape = festimator.getShape(cv::Mat(99, 1, CV_32F));
    tex = shape * 0 + 128;

    // Initialize image renderer
    if (im_render == nullptr)
    {
        im_render = new FImRenderer(cv::Mat::zeros(h, w, CV_8UC3));
        im_render->loadMesh(shape, shape * 0, faces);
    }
    else im_render->init(cv::Mat::zeros(h, w, CV_8UC3));
}

bool FaceServices2::projectCheckVis(FImRenderer* im_render, cv::Mat shape, float* r, float *t, cv::Mat refDepth, bool* &visible){
	float zNear_ = im_render->zNear;
	float zFar_ = im_render->zFar;

	cv::Mat k_m( 3, 3, CV_32F, _k );

	int nV = shape.rows;
	if (visible == 0) visible = new bool[shape.rows];
	cv::Mat rVec( 3, 1, CV_32F, r );
	cv::Mat tVec( 3, 1, CV_32F, t );
	cv::Mat rMat;
	cv::Rodrigues(rVec, rMat);
	cv::Mat new3D = rMat * shape.t() + cv::repeat(tVec,1,nV);

	for (int i=0;i<nV;i++){
		visible[i] = false;
		float Z = new3D.at<float>(2,i);
		float x = -new3D.at<float>(0,i)/Z*_k[4] + _k[2];
		float y = new3D.at<float>(1,i)/Z*_k[4] + _k[5];
		if (x > 0 && y > 0 & x < refDepth.cols-1 && y <refDepth.rows-1) {
			for (int dx =-1;dx<2;dx++){
				for (int dy =-1;dy<2;dy++){
					float dd = refDepth.at<float>(y+dy,x+dx);
					dd = - zNear_*zFar_   / ( zFar_ - dd * ( zFar_ - zNear_ ));
					if (fabs(Z - dd) < 5){
						visible[i] = true;
					}
				}
			}
		}
	}

	return true;
}

std::vector<cv::Point2f> FaceServices2::projectCheckVis2(FImRenderer* im_render, cv::Mat shape, float* r, float *t, cv::Mat refDepth, bool* &visible){
	float zNear_ = im_render->zNear;
	float zFar_ = im_render->zFar;

	cv::Mat k_m( 3, 3, CV_32F, _k );
	std::vector<cv::Point2f> out;
	int nV = shape.rows;
	if (visible == 0) visible = new bool[shape.rows];
	cv::Mat rVec( 3, 1, CV_32F, r );
	cv::Mat tVec( 3, 1, CV_32F, t );
	cv::Mat rMat;
	cv::Rodrigues(rVec, rMat);
	cv::Mat new3D = rMat * shape.t() + cv::repeat(tVec,1,nV);

	for (int i=0;i<nV;i++){
		visible[i] = false;
		float Z = new3D.at<float>(2,i);
		float x = -new3D.at<float>(0,i)/Z*_k[4] + _k[2];
		float y = new3D.at<float>(1,i)/Z*_k[4] + _k[5];
		out.push_back(cv::Point2f(x,y));
		if (x > 0 && y > 0 & x < refDepth.cols-1 && y <refDepth.rows-1) {
			for (int dx =-1;dx<2;dx++){
				for (int dy =-1;dy<2;dy++){
					float dd = refDepth.at<float>(y+dy,x+dx);
					dd = - zNear_*zFar_   / ( zFar_ - dd * ( zFar_ - zNear_ ));
					if (fabs(Z - dd) < 2.5){
						visible[i] = true;
					}
				}
			}
		}
	}

	return out;
}

float FaceServices2::updateHessianMatrix(bool part, cv::Mat alpha, float* renderParams, cv::Mat faces, cv::Mat colorIm,std::vector<int> lmInds, cv::Mat landIm, BFMParams &params, cv::Mat &prevR, cv::Mat &prevT, cv::Mat exprW ){
	int M = alpha.rows;
	int EM = exprW.rows;
	float step;
	cv::Mat k_m( 3, 3, CV_32F, _k );
	cv::Mat distCoef = cv::Mat::zeros( 1, 4, CV_32F );
	params.hessDiag.release();
	params.hessDiag = cv::Mat::zeros(2*M+EM+RENDER_PARAMS_COUNT,1,CV_32F);
	//printf("sno_step --------\n");
	//std::cout << "alpha " << alpha.t() <<std::endl;
	//std::cout << "beta " << beta.t() <<std::endl;
	//printf("renderParams ");
	//for (int i=0;i<RENDER_PARAMS_COUNT;i++)
	//	printf("(%d) %f, ",i,renderParams[i]);
	//printf("\n");
	cv::Mat alpha2, expr2;
	float renderParams2[RENDER_PARAMS_COUNT];

	cv::Mat rVec(3,1,CV_32F,renderParams+RENDER_PARAMS_R);
	cv::Mat tVec(3,1,CV_32F,renderParams+RENDER_PARAMS_T);
	cv::Mat rVec2(3,1,CV_32F,renderParams2+RENDER_PARAMS_R);
	cv::Mat tVec2(3,1,CV_32F,renderParams2+RENDER_PARAMS_T);

	float currEF = eF(part, alpha, lmInds, landIm, renderParams, exprW);
	//printf("currF: %f\n",currEF);
	cEF = currEF;

	// expr
	step = mstep*5;
	if (params.optimizeExpr) {
		for (int i=0;i<EM; i++){
			expr2.release(); expr2 = exprW.clone();
			expr2.at<float>(i,0) += step;
			float tmpEF1 = eF(part, alpha, lmInds, landIm, renderParams,expr2);
			expr2.at<float>(i,0) -= 2*step;
			float tmpEF2 = eF(part, alpha, lmInds, landIm, renderParams,expr2);
			params.hessDiag.at<float>(2*M+i,0) = params.sF[FEATURES_LANDMARK] * (tmpEF1 - 2*currEF + tmpEF2)/(step*step)
				+ params.sExpr * 2/(0.25f*29) ;
		}
	}
	// r
	//step = 0.05;
	step = mstep*2;
	//step = 0.02;
	//step = 0.01;
	if (params.doOptimize[RENDER_PARAMS_R]) {
		for (int i=0;i<3; i++){
			memcpy(renderParams2,renderParams,RENDER_PARAMS_COUNT*sizeof(float));
			renderParams2[RENDER_PARAMS_R+i] += step;
			float tmpEF1 = eF(part, alpha, lmInds, landIm, renderParams2,exprW);

			renderParams2[RENDER_PARAMS_R+i] -= 2*step;
			float tmpEF2 = eF(part, alpha, lmInds, landIm, renderParams2,exprW);
			params.hessDiag.at<float>(2*M+EM+i,0) = params.sF[FEATURES_LANDMARK] * (tmpEF1 - 2*currEF + tmpEF2)/(step*step) + 2.0f/params.sR[RENDER_PARAMS_R+i];

		}
	}
	// t
	step = mstep*10;
	//step = 0.05;
	//step = 0.1;
	if (params.doOptimize[RENDER_PARAMS_T]) {
		for (int i=0;i<3; i++){
			memcpy(renderParams2,renderParams,RENDER_PARAMS_COUNT*sizeof(float));
			renderParams2[RENDER_PARAMS_T+i] += step;
			float tmpEF1 = eF(part, alpha, lmInds, landIm, renderParams2,exprW);
			renderParams2[RENDER_PARAMS_T+i] -= 2*step;
			float tmpEF2 = eF(part, alpha, lmInds, landIm, renderParams2,exprW);
			params.hessDiag.at<float>(2*M+EM+RENDER_PARAMS_T+i,0) = params.sF[FEATURES_LANDMARK] * (tmpEF1 - 2*currEF + tmpEF2)/(step*step)
				+ 2.0f/params.sR[RENDER_PARAMS_T+i];
		}
	}
	return 0;
}

FaceServices2::~FaceServices2(void)
{
}

cv::Mat FaceServices2::computeGradient(bool part, cv::Mat alpha, float* renderParams, cv::Mat faces,cv::Mat colorIm, std::vector<int> lmInds, cv::Mat landIm, BFMParams &params, std::vector<int> &inds, cv::Mat exprW, cv::Mat &prevR, cv::Mat &prevT){
	int M = alpha.rows;
	int EM = exprW.rows;
	int nTri = 40;
	float step;
	cv::Mat k_m( 3, 3, CV_32F, _k );
	cv::Mat distCoef = cv::Mat::zeros( 1, 4, CV_32F );
	cv::Mat out(2*M+EM+RENDER_PARAMS_COUNT,1,CV_32F);

	cv::Mat alpha2, expr2;
	float renderParams2[RENDER_PARAMS_COUNT];
	cv::Mat rVec(3,1,CV_32F,renderParams+RENDER_PARAMS_R);
	cv::Mat tVec(3,1,CV_32F,renderParams+RENDER_PARAMS_T);
	cv::Mat rVec2(3,1,CV_32F,renderParams2+RENDER_PARAMS_R);
	cv::Mat tVec2(3,1,CV_32F,renderParams2+RENDER_PARAMS_T);

	float currEF = eF(part, alpha, lmInds, landIm, renderParams,exprW);
	cEF = currEF;

	#pragma omp parallel for
	for (int target=0;target<EM+6; target++){
	  if (target < EM) {
		// expr
		float step = mstep*5;
		if (params.optimizeExpr) {
				int i = target;
				std::vector<cv::Point2f> pPoints;
				cv::Mat expr2 = exprW.clone();
				expr2.at<float>(i,0) += step;
				float tmpEF = eF(part, alpha, lmInds, landIm, renderParams,expr2);
				out.at<float>(2*M+i,0) = params.sF[FEATURES_LANDMARK] * (tmpEF - currEF)/step
					+ params.sExpr * 2*exprW.at<float>(i,0)/(0.25f*29);
				//if (cEF > prevEF) printf("%f, ",tmpEF);
		}
	   }
	   else if (target < EM+3) {
		// r
		//step = 0.05;
		//step = 0.01;
		float step = mstep*2;
		//step = 0.01;
		if (params.doOptimize[RENDER_PARAMS_R]) {
			int i = target-EM;
				float renderParams2[RENDER_PARAMS_COUNT];
				memcpy(renderParams2,renderParams,RENDER_PARAMS_COUNT*sizeof(float));
				renderParams2[RENDER_PARAMS_R+i] += step;
				float tmpEF = eF(part, alpha, lmInds, landIm, renderParams2,exprW);
				out.at<float>(2*M+EM+i,0) = params.sF[FEATURES_LANDMARK] * (tmpEF - currEF)/step;
				if (prevR.rows == 0)
					out.at<float>(2*M+EM+i,0) += 2*(renderParams[RENDER_PARAMS_R+i] - params.initR[RENDER_PARAMS_R+i])/params.sR[RENDER_PARAMS_R+i];
				else {
					float val = renderParams[RENDER_PARAMS_R+i] - prevR.at<float>(i,0);
					if (val > PREV_USE_THRESH)
						out.at<float>(2*M+EM+i,0) += 2*(renderParams[RENDER_PARAMS_R+i] - params.initR[RENDER_PARAMS_R+i])/params.sR[RENDER_PARAMS_R+i];
					else {
						//printf("usePrev %f (%d %d)\n",prevR.at<float>(i,0),REG_FROM_CURR,REG_FROM_PREV);
						out.at<float>(2*M+EM+i,0) += 2*REG_FROM_CURR*(renderParams[RENDER_PARAMS_R+i] - params.initR[RENDER_PARAMS_R+i])/params.sR[RENDER_PARAMS_R+i] + 2*REG_FROM_PREV*val/params.sR[RENDER_PARAMS_R+i];
					}
				}
				//if (cEF > prevEF) printf("%f, ",tmpEF);
		}
	  }
	  else {
		// t
		float step = mstep*10;
		//step = 1;
		if (params.doOptimize[RENDER_PARAMS_T]) {
			int i = target-EM-3;
				float renderParams2[RENDER_PARAMS_COUNT];
				memcpy(renderParams2,renderParams,RENDER_PARAMS_COUNT*sizeof(float));
				renderParams2[RENDER_PARAMS_T+i] += step;
				float tmpEF = eF(part, alpha, lmInds, landIm, renderParams2,exprW);
				out.at<float>(2*M+EM+RENDER_PARAMS_T+i,0) = params.sF[FEATURES_LANDMARK] * (tmpEF - currEF)/step
					+ 2*(renderParams[RENDER_PARAMS_T+i] - params.initR[RENDER_PARAMS_T+i])/params.sR[RENDER_PARAMS_T+i];
				//if (cEF > prevEF) printf("%f, ",tmpEF);
		}
	  }
	}
	return out;
}

float FaceServices2::eF(bool part, cv::Mat alpha, std::vector<int> inds, cv::Mat landIm, float* renderParams, cv::Mat exprW){
	Mat k_m(3,3,CV_32F,_k);
	//printf("%f\n",renderParams[RENDER_PARAMS_R+1]);
	cv::Mat mLM;
	if (!part)
		mLM = festimator.getLMByAlpha(alpha,-renderParams[RENDER_PARAMS_R+1], inds, exprW);
	else
		mLM = festimator.getLMByAlphaParts(alpha,-renderParams[RENDER_PARAMS_R+1], inds, exprW);
	//write_plyFloat("vismLM.ply",mLM.t());
	//printf("inds\n");
	//for (int i=0;i<inds.size(); i++)
//		printf("%d\n", inds[i]);
	//printf("renderParams %f %f %f   %f %f %f\n", renderParams[0],renderParams[1],renderParams[2],renderParams[3],renderParams[4],renderParams[5]);
	//std::cout << "landIm " << landIm << std::endl;
	//cv::Mat ads = cv::Mat::zeros(720,540,CV_8UC3);

	cv::Mat rVec(3,1,CV_32F, renderParams + RENDER_PARAMS_R);
	cv::Mat tVec(3,1,CV_32F, renderParams + RENDER_PARAMS_T);
	std::vector<cv::Point2f> allImgPts;
	cv::Mat distCoef = cv::Mat::zeros( 1, 4, CV_32F );

	cv::projectPoints( mLM, rVec, tVec, k_m, distCoef, allImgPts );
	float err = 0;
	for (int i=0;i<mLM.rows;i++){
		float val = landIm.at<float>(i,0) - allImgPts[i].x;
		err += val*val;
		val = landIm.at<float>(i,1) - allImgPts[i].y;
		err += val*val;
		//cv::circle(ads,cv::Point(landIm.at<float>(i,0),landIm.at<float>(i,1)),2,cv::Scalar(255,0,0));
		//cv::circle(ads,cv::Point(allImgPts[i].x,allImgPts[i].y),2,cv::Scalar(0,0,255));
	}
	//imwrite("sas.png",ads);
	//getchar();
	return sqrt(err/mLM.rows);
}

void FaceServices2::sno_step2(bool part, cv::Mat &alpha, float* renderParams, cv::Mat faces,cv::Mat colorIm, std::vector<int> lmInds, cv::Mat landIm, BFMParams &params, cv::Mat &exprW, cv::Mat &prevR, cv::Mat &prevT){
	//printf("sno_step --------\n");
	//std::cout << "alpha " << alpha.t() <<std::endl;
	//std::cout << "beta " << beta.t() <<std::endl;
	//printf("renderParams ");
	//for (int i=0;i<RENDER_PARAMS_COUNT;i++)
	//	printf("%f, ",renderParams[i]);
	//printf("\n");
	float lambda = 0.05;
	std::vector<int> inds;
	cv::Mat dE = computeGradient(part, alpha, renderParams, faces, colorIm, lmInds, landIm, params,inds, exprW, prevR, prevT);
	params.gradVec.release(); params.gradVec = dE.clone();
	cv::Mat dirMove = dE*0;

	int M = alpha.rows;
	int EM = exprW.rows;
	if (params.optimizeExpr){
		for (int i=0;i<EM;i++)
			if (abs(params.hessDiag.at<float>(2*M+i,0)) > 0.0000001) {
				dirMove.at<float>(2*M+i,0) = - lambda*dE.at<float>(2*M+i,0)/abs(params.hessDiag.at<float>(2*M+i,0));
			}
	}

	for (int i=0;i<RENDER_PARAMS_COUNT;i++) {
		if (params.doOptimize[i]){
			if (abs(params.hessDiag.at<float>(2*M+EM+i,0)) > 0.0000001) {
				dirMove.at<float>(2*M+EM+i,0) = - lambda*dE.at<float>(2*M+EM+i,0)/abs(params.hessDiag.at<float>(2*M+EM+i,0));
			}
		}
	}
	//float pc = 1;
	float pc = line_search(part, alpha, renderParams, dirMove,inds, faces, colorIm, lmInds, landIm, params, exprW, prevR, prevT, 10);
	//printf("pc = %f\n",pc);
	if (pc == 0) countFail++;
	else {
		if (params.optimizeExpr){
			for (int i=0;i<EM;i++) {
				exprW.at<float>(i,0) += pc*dirMove.at<float>(i+2*M,0);
				if (exprW.at<float>(i,0) > 3) exprW.at<float>(i,0) = 3;
				else if (exprW.at<float>(i,0) < -3) exprW.at<float>(i,0) = -3;
			}
		}

		for (int i=0;i<RENDER_PARAMS_COUNT;i++) {
			if (params.doOptimize[i]){
				renderParams[i] += pc*dirMove.at<float>(2*M+EM+i,0);
				if (i == RENDER_PARAMS_CONTRAST || (i>=RENDER_PARAMS_AMBIENT && i<RENDER_PARAMS_DIFFUSE+3) ) {
					if (renderParams[i] > 1.0) renderParams[i] = 1.0;
					if (renderParams[i] < 0) renderParams[i] = 0;

				}
				else if (i >= RENDER_PARAMS_GAIN && i<=RENDER_PARAMS_GAIN+3) {
					if (renderParams[i] > 3.0) renderParams[i]  = 3;
					if (renderParams[i] < 0.3) renderParams[i] = 0.3;
				}
			}
		}
	}
	prevEF = cEF;
}

float FaceServices2::line_search(bool part, cv::Mat &alpha, float* renderParams, cv::Mat &dirMove,std::vector<int> inds, cv::Mat faces,cv::Mat colorIm,std::vector<int> lmInds, cv::Mat landIm, BFMParams &params, cv::Mat &exprW, cv::Mat &prevR, cv::Mat &prevT, int maxIters){
	float step = 1.0f;
	float sstep = 2.0f;
	float minStep = 0.0001f;
	//float cCost = computeCost(cEF, cEI,cETE, cECE, alpha, beta, renderParams,params);
	cv::Mat alpha2, exprW2;
	float renderParams2[RENDER_PARAMS_COUNT];
	alpha2 = alpha.clone();
	exprW2 = exprW.clone();
	memcpy(renderParams2,renderParams,sizeof(float)*RENDER_PARAMS_COUNT);

	cv::Mat k_m( 3, 3, CV_32F, _k );
	cv::Mat distCoef = cv::Mat::zeros( 1, 4, CV_32F );
	std::vector<cv::Point2f> pPoints;
	cv::Mat rVec2(3,1,CV_32F,renderParams2+RENDER_PARAMS_R);
	cv::Mat tVec2(3,1,CV_32F,renderParams2+RENDER_PARAMS_T);

	int M = alpha.rows;
	int EM = exprW.rows;
	float ssize = 0;
	for (int i=0;i<dirMove.rows;i++) ssize += dirMove.at<float>(i,0)*dirMove.at<float>(i,0);
	ssize = sqrt(ssize);
	//printf("ssize: %f\n",ssize);
	if (ssize > (2*M+EM+RENDER_PARAMS_COUNT)/5.0f) {
		step = (2*M+EM+RENDER_PARAMS_COUNT)/(5.0f * ssize);
		ssize = (2*M+EM+RENDER_PARAMS_COUNT)/5.0f;
	}
	if (ssize < minStep){
		return 0;
	}
	int tstep = floor(log(ssize/minStep));
	if (tstep < maxIters) maxIters = tstep;

	float curCost = computeCost(cEF, alpha, renderParams, params, exprW, prevR, prevT );

	bool hasNoBound = false;
	int iter = 0;
	for (; iter<maxIters; iter++){
		if (params.optimizeExpr){
			for (int i=0;i<EM;i++) {
				float tmp = exprW.at<float>(i,0) + step*dirMove.at<float>(2*M+i,0);
				if (tmp >= 3) exprW2.at<float>(i,0) = 3;
				else if (tmp <= -3) exprW2.at<float>(i,0) = -3;
				else {
					exprW2.at<float>(i,0) = tmp;
					hasNoBound = true;
				}
			}
		}

		for (int i=0;i<RENDER_PARAMS_COUNT;i++) {
			if (params.doOptimize[i]){
				float tmp = renderParams[i] + step*dirMove.at<float>(2*M+EM+i,0);
				if (i == RENDER_PARAMS_CONTRAST || (i>=RENDER_PARAMS_AMBIENT && i<RENDER_PARAMS_DIFFUSE+3) ) {
					if (tmp > 1.0) renderParams2[i] = 1.0f;
					else if (tmp < -1.0) renderParams2[i] = -1.0f;
					else {
						renderParams2[i] = tmp;
						hasNoBound = true;
					}
				}
				else if (i >= RENDER_PARAMS_GAIN && i<=RENDER_PARAMS_GAIN+3) {
					if (tmp >= 3.0) renderParams2[i] = 3.0f;
					else if (tmp <= -3.0) renderParams2[i] = -3.0f;
					else {
						renderParams2[i] = tmp;
						hasNoBound = true;
					}
				}
				else renderParams2[i] = tmp;
			}
		}
		if (!hasNoBound) {
			iter = maxIters; break;
		}
		float tmpEF = cEF;
		if (params.sF[FEATURES_LANDMARK] > 0) tmpEF = eF(part,alpha2, lmInds,landIm,renderParams2, exprW2);
		float tmpCost = computeCost(tmpEF, alpha2, renderParams2, params,exprW2, prevR, prevT);
		//printf("tmpCost %f\n",tmpCost);
		if (tmpCost < curCost) {
			break;
		}
		else {
			step = step/sstep;
			//printf("step %f\n",step);
		}
	}
	//getchar();
	if (iter >= maxIters) return 0;
	else return step;
}

float FaceServices2::computeCost(float vEF, cv::Mat &alpha, float* renderParams, BFMParams &params, cv::Mat &exprW, cv::Mat &prevR, cv::Mat &prevT ){
	float val = params.sF[FEATURES_LANDMARK]*vEF;
	int M = alpha.rows;
	if (params.optimizeExpr){
		for (int i=0;i<exprW.rows;i++)
			val += params.sExpr * exprW.at<float>(i,0)*exprW.at<float>(i,0)/(0.5f*29);
	}

	for (int i=0;i<3;i++) {
		if (params.doOptimize[i]){
			if (prevR.rows == 0) val += (renderParams[i] - params.initR[i])*(renderParams[i] - params.initR[i])/params.sR[i];
			else {
				float vval = renderParams[i] - prevR.at<float>(i,0);
				if (vval > PREV_USE_THRESH)
					val += (renderParams[i] - params.initR[i])*(renderParams[i] - params.initR[i])/params.sR[i];
				else
					val += REG_FROM_CURR*(renderParams[i] - params.initR[i])*(renderParams[i] - params.initR[i])/params.sR[i] + REG_FROM_PREV*vval*vval/params.sR[i];
			}
		}
	}
	for (int i=3;i<RENDER_PARAMS_COUNT;i++) {
		if (params.doOptimize[i]){
			val += (renderParams[i] - params.initR[i])*(renderParams[i] - params.initR[i])/params.sR[i];
		}
	}
	return val;
}

bool FaceServices2::loadReference(string refDir, string model_file, cv::Mat &alpha, cv::Mat &beta, float* renderParams, int &M, cv::Mat &exprW, int &EM){
	string fname(model_file);
	size_t sep = model_file.find_last_of("\\/");
	if (sep != std::string::npos)
		fname = model_file.substr(sep + 1, model_file.size() - sep - 1);
	fname = refDir + fname;

	alpha = cv::Mat::zeros(M,1,CV_32F);
	beta = cv::Mat::zeros(M,1,CV_32F);
	char fpath[250];
	char text[250];
	sprintf(fpath,"%s.alpha",fname.c_str());
	//printf("fpath %s\n",fpath);
	FILE* f = fopen(fpath,"r");
	if (f == 0) return false;

	text[0] = '\0';
	for (int i=0;i<M;i++){
		fgets(text,250,f);
		int l = strlen(text);
		if (text[l-1] <'0' || text[l-1] > '9') text[l-1] = '\0';
		alpha.at<float>(i,0) = atof(text);
	}
	fclose(f);

	sprintf(fpath,"%s.beta",fname.c_str());
	f = fopen(fpath,"r");
	text[0] = '\0';
	for (int i=0;i<M;i++){
		fgets(text,250,f);
		int l = strlen(text);
		if (text[l-1] <'0' || text[l-1] > '9') text[l-1] = '\0';
		beta.at<float>(i,0) = atof(text);
	}
	fclose(f);

	sprintf(fpath,"%s.rend",fname.c_str());
	f = fopen(fpath,"r");
	text[0] = '\0';
	for (int i=0;i<RENDER_PARAMS_COUNT;i++){
		fgets(text,250,f);
		int l = strlen(text);
		if (text[l-1] <'0' || text[l-1] > '9') text[l-1] = '\0';
		renderParams[i] = atof(text);
	}
	fclose(f);

	sprintf(fpath,"%s.expr",fname.c_str());
	//printf("fpath %s\n",fpath);
	f = fopen(fpath,"r");
	if (f == 0) return false;

	text[0] = '\0';
	for (int i=0;i<EM;i++){
		fgets(text,250,f);
		int l = strlen(text);
		if (text[l-1] <'0' || text[l-1] > '9') text[l-1] = '\0';
		exprW.at<float>(i,0) = atof(text);
	}
	fclose(f);
	return true;
}

bool FaceServices2::loadReference2(string refDir, string model_file, cv::Mat &alpha, cv::Mat &beta, int &M){
	string fname(model_file);
	size_t sep = model_file.find_last_of("\\/");
	if (sep != std::string::npos)
		fname = model_file.substr(sep + 1, model_file.size() - sep - 1);
	fname = refDir + fname;

	alpha = cv::Mat::zeros(M,1,CV_32F);
	beta = cv::Mat::zeros(M,1,CV_32F);
	char fpath[250];
	char text[250];
	sprintf(fpath,"%s.alpha",fname.c_str());
	//printf("fpath %s\n",fpath);
	FILE* f = fopen(fpath,"r");
	if (f == 0) return false;

	text[0] = '\0';
	for (int i=0;i<M;i++){
		fgets(text,250,f);
		int l = strlen(text);
		if (text[l-1] <'0' || text[l-1] > '9') text[l-1] = '\0';
		alpha.at<float>(i,0) = atof(text);
	}
	fclose(f);

	sprintf(fpath,"%s.beta",fname.c_str());
	f = fopen(fpath,"r");
	text[0] = '\0';
	for (int i=0;i<M;i++){
		fgets(text,250,f);
		int l = strlen(text);
		if (text[l-1] <'0' || text[l-1] > '9') text[l-1] = '\0';
		beta.at<float>(i,0) = atof(text);
	}
	fclose(f);
	return true;
}
void FaceServices2::initRenderer(cv::Mat &colorIm){
    cv::Mat faces = festimator.getFaces() - 1;
    cv::Mat faces_fill = festimator.getFaces_fill() - 1;
    cv::Mat colors;
    cv::Mat shape = festimator.getShape(cv::Mat::zeros(1,1,CV_32F));
    cv::Mat tex = festimator.getTexture(cv::Mat::zeros(1,1,CV_32F));

    im_render = new FImRenderer(cv::Mat::zeros(colorIm.rows,colorIm.cols,CV_8UC3));
    im_render->loadMesh(shape,tex,faces_fill);
}

void FaceServices2::mergeIm(cv::Mat* output,cv::Mat bg,cv::Mat depth){
	for (int i=0;i<bg.rows;i++){
		for (int j=0;j<bg.cols;j++){
			if (depth.at<float>(i,j) >= 0.9999)
				output->at<Vec3b>(i, j) = bg.at<Vec3b>(i,j);
		}
	}
}

bool FaceServices2::estimatePoseExpr(cv::Mat colorIm, cv::Mat lms, cv::Mat alpha, cv::Mat &vecR, cv::Mat &vecT, cv::Mat& K, cv::Mat &exprW, const char* outputDir, bool with_expr){
	char text[200];
	float renderParams[RENDER_PARAMS_COUNT];
	float renderParams2[RENDER_PARAMS_COUNT];
	Mat k_m(3,3,CV_32F,_k);
    K = k_m.clone();    // Yuval
	BFMParams params;
	params.init();
	exprW = cv::Mat::zeros(29,1,CV_32F);
	cv::Mat prevR;
	cv::Mat prevT;

	int M = 99;
	cv::Mat shape = festimator.getShape(alpha);

	Mat landModel0 = festimator.getLM(shape,0);
	int nLM = landModel0.rows;
	Mat landIm = cv::Mat( 60,2,CV_32F);
	Mat landModel = cv::Mat( 60,3,CV_32F);
	for (int i=0;i<60;i++){
		landModel.at<float>(i,0) = landModel0.at<float>(i,0);
		landModel.at<float>(i,1) = landModel0.at<float>(i,1);
		landModel.at<float>(i,2) = landModel0.at<float>(i,2);
		landIm.at<float>(i,0) = lms.at<float>(i,0);
		landIm.at<float>(i,1) = lms.at<float>(i,1);
	}
	festimator.estimatePose3D0(landModel,landIm,k_m,vecR,vecT);
	float yaw = -vecR.at<float>(1,0);
	landModel0 = festimator.getLM(shape,yaw);
	std::vector<int> lmVisInd;
	for (int i=0;i<60;i++){
		if ((yaw > M_PI/10 && i > 7) || (yaw < -M_PI/10 && i < 9) || i > 16 || abs(yaw) <= M_PI/10)
			lmVisInd.push_back(i);
	}
	cv::Mat tmpIm = colorIm/5;
	landModel = cv::Mat( lmVisInd.size(),3,CV_32F);
	landIm = cv::Mat::zeros( lmVisInd.size(),2,CV_32F);
	for (int i=0;i<lmVisInd.size();i++){
		int ind = lmVisInd[i];
		landModel.at<float>(i,0) = landModel0.at<float>(ind,0);
		landModel.at<float>(i,1) = landModel0.at<float>(ind,1);
		landModel.at<float>(i,2) = landModel0.at<float>(ind,2);
		landIm.at<float>(i,0) = lms.at<float>(ind,0);
		landIm.at<float>(i,1) = lms.at<float>(ind,1);
		//cv::circle(tmpIm,Point(landIm.at<float>(i,0),landIm.at<float>(i,1)),2,Scalar(255,0,0),2);
	}
	//sprintf(text,"%s/withLM.png",outputDir);
	//imwrite(text,tmpIm);
	//sprintf(text,"%s/%s_withLM.ply",outputDir,filename.c_str());
	//write_plyFloat(text,landModel.t());
	//getchar();

	festimator.estimatePose3D0(landModel,landIm,k_m,vecR,vecT);

    // Yuval
    if (!with_expr) return true;

	for (int i=0;i<3;i++)
		params.initR[RENDER_PARAMS_R+i] = vecR.at<float>(i,0);
	for (int i=0;i<3;i++)
		params.initR[RENDER_PARAMS_T+i] = vecT.at<float>(i,0);
	memcpy(renderParams,params.initR,sizeof(float)*RENDER_PARAMS_COUNT);

	for (int i=60;i<68;i++) lmVisInd.push_back(i);
	landIm = cv::Mat::zeros( lmVisInd.size(),2,CV_32F);
	for (int i=0;i<lmVisInd.size();i++){
		int ind = lmVisInd[i];
		landIm.at<float>(i,0) = lms.at<float>(ind,0);
		landIm.at<float>(i,1) = lms.at<float>(ind,1);
	}

	float bCost, cCost, fCost;
	int bestIter = 0;
	bCost = 10000.0f;

	memset(params.sF,0,sizeof(float)*NUM_EXTRA_FEATURES);

	params.sI = 0.0;
	params.sF[FEATURES_LANDMARK] = 8.0f;
	Mat alpha0;
	int iter=0;
	int badCount = 0;
	memset(params.doOptimize,true,sizeof(bool)*6);

	int EM = 29;
	float renderParams_tmp[RENDER_PARAMS_COUNT];

	params.optimizeAB[0] = params.optimizeAB[1] = false;
	for (;iter<60;iter++) {
			if (iter%20 == 0) {
				cCost = updateHessianMatrix(false, alpha,renderParams,faces,colorIm,lmVisInd,landIm,params, prevR, prevT, exprW);
				if (countFail > 10) {
					countFail = 0;
					break;
				}
				prevEF = cEF;
			}
			sno_step2(false, alpha, renderParams, faces,colorIm,lmVisInd,landIm,params,exprW,prevR, prevT);
		}
	iter = 60;
	memset(params.doOptimize,false,sizeof(bool)*6);countFail = 0;
    //int total_iter = 200;
    int total_iter = 400;   // Yuval
	for (;iter<total_iter;iter++) {
			if (iter%60 == 0) {
				cCost = updateHessianMatrix(false, alpha,renderParams,faces,colorIm,lmVisInd,landIm,params, prevR, prevT, exprW);
				if (countFail > 10) {
					countFail = 0;
					break;
				}
				prevEF = cEF;
			}
			sno_step2(false, alpha, renderParams, faces,colorIm,lmVisInd,landIm,params,exprW,prevR, prevT);
		}

	///for (int i=0;i<3; i++) vecR.at<float>(i,0) = renderParams[i];
	///for (int i=0;i<3; i++) vecT.at<float>(i,0) = renderParams[i+3];
	//std::cout << "vecR " << vecR << std::endl;
	//std::cout << "vecT " << vecT << std::endl;
	//std::cout << "exprW " << exprW << std::endl;

    return true;
}

bool FaceServices2::updatePoseExpr(cv::Mat colorIm, cv::Mat lms, cv::Mat alpha, cv::Mat &vecR, cv::Mat &vecT, cv::Mat &exprW, char* outputDir, cv::Mat &prevR, cv::Mat &prevT )
{
	char text[200];
	float renderParams[RENDER_PARAMS_COUNT];
	float renderParams2[RENDER_PARAMS_COUNT];
	Mat k_m(3,3,CV_32F,_k);
	BFMParams params;
	params.init();

	int M = 99;
	float yaw = -vecR.at<float>(1,0);
	Mat landModel0 = festimator.getLM(shape,yaw);
	std::vector<int> lmVisInd;
	for (int i=0;i<60;i++){
		if (i > 16 || std::abs(yaw) <= M_PI/10 || (yaw > M_PI/10 && i > 7) || (yaw < -M_PI/10 && i < 9))
			lmVisInd.push_back(i);
	}
	cv::Mat landModel = cv::Mat( lmVisInd.size(),3,CV_32F);
	cv::Mat landIm = cv::Mat::zeros( lmVisInd.size(),2,CV_32F);
	for (int i=0;i<lmVisInd.size();i++){
		int ind = lmVisInd[i];
		landModel.at<float>(i,0) = landModel0.at<float>(ind,0);
		landModel.at<float>(i,1) = landModel0.at<float>(ind,1);
		landModel.at<float>(i,2) = landModel0.at<float>(ind,2);
		landIm.at<float>(i,0) = lms.at<float>(ind,0);
		landIm.at<float>(i,1) = lms.at<float>(ind,1);
		//cv::circle(tmpIm,Point(landIm.at<float>(i,0),landIm.at<float>(i,1)),2,Scalar(255,0,0),2);
	}
	//sprintf(text,"%s/withLM.png",outputDir);
	//imwrite(text,tmpIm);
	//sprintf(text,"%s/%s_withLM.ply",outputDir,filename.c_str());
	//write_plyFloat(text,landModel.t());
	//getchar();

	festimator.estimatePose3D0(landModel,landIm,k_m,vecR,vecT);
	for (int i=60;i<68;i++) lmVisInd.push_back(i);

	for (int i=0;i<3;i++)
		params.initR[RENDER_PARAMS_R+i] = vecR.at<float>(i,0);
	for (int i=0;i<3;i++)
		params.initR[RENDER_PARAMS_T+i] = vecT.at<float>(i,0);
	memcpy(renderParams,params.initR,sizeof(float)*RENDER_PARAMS_COUNT);

	landIm = cv::Mat::zeros( lmVisInd.size(),2,CV_32F);
	for (int i=0;i<lmVisInd.size();i++){
		int ind = lmVisInd[i];
		landIm.at<float>(i,0) = lms.at<float>(ind,0);
		landIm.at<float>(i,1) = lms.at<float>(ind,1);
	}

	float bCost, cCost, fCost;
	int bestIter = 0;
	bCost = 10000.0f;

	memset(params.sF,0,sizeof(float)*NUM_EXTRA_FEATURES);

	params.sI = 0.0;
	params.sF[FEATURES_LANDMARK] = 8.0f;
	int iter=0;
	int badCount = 0;
	memset(params.doOptimize,true,sizeof(bool)*6);

	int EM = 29;
	float renderParams_tmp[RENDER_PARAMS_COUNT];

	params.optimizeAB[0] = params.optimizeAB[1] = false;

    for (; iter < 30; iter++)
    {
        if (iter % 10 == 0) {
            cCost = updateHessianMatrix(false, alpha, renderParams, faces, colorIm, lmVisInd, landIm, params, prevR, prevT, exprW);
            if (countFail > 10) {
                countFail = 0;

                break;
            }
            prevEF = cEF;
        }
        sno_step2(false, alpha, renderParams, faces, colorIm, lmVisInd, landIm, params, exprW, prevR, prevT);

    }

    return true;
}
cv::Mat FaceServices2::testRender(cv::Mat colorIm, cv::Mat alpha,cv::Mat vecR,cv::Mat vecT,cv::Mat exprW, char* outDir, bool updateColor){
	char text[200];
	int M = 99;
	float renderParams[RENDER_PARAMS_COUNT];
	for (int i =0;i<3;i++)
		renderParams[i] = vecR.at<float>(i,0);
	for (int i =0;i<3;i++)
		renderParams[i+3] = vecT.at<float>(i,0);

	////cv::Mat shape = festimator.getShape(alpha,exprW);
	//cv::Mat faces_fill = festimator.getFaces_fill() - 1;
	// Ambient
	//renderParams[6] = 0.29225;
	//renderParams[7] = 0.29225;
	//renderParams[8] = 0.29225;
	renderParams[6] = 0.49225;
	renderParams[7] = 0.49225;
	renderParams[8] = 0.49225;
	// Diffuse
	renderParams[9] = 0.50754;
	renderParams[10] = 0.50754;
	renderParams[11] = 0.50754;
	// LIGHT
	renderParams[12] = 3.1415/4;
	renderParams[13] = 3.1415/4;
	// OTHERS
	renderParams[14] = 1;
	renderParams[15] = renderParams[16] = renderParams[17] = RENDER_PARAMS_GAIN_DEFAULT;
	renderParams[18] = renderParams[19] = renderParams[20] = RENDER_PARAMS_OFFSET_DEFAULT;

	Mat k_m(3,3,CV_32F,_k);
	//RenderServices rs;
	float* r = renderParams + RENDER_PARAMS_R;
	float* t = renderParams + RENDER_PARAMS_T;
	//printf("t %f %f %f\n",t[0],t[1],t[2]);

		shape = festimator.getShape2(alpha,exprW);
	//if (updateColor) {
		cv::Mat colors;
		tex = shape*0 + 128;

		im_render->copyShape(shape);
		bool* visible = new bool[im_render->face_->mesh_.nVertices_];
		bool* noShadow = new bool[im_render->face_->mesh_.nVertices_];
		memset(noShadow,true,im_render->face_->mesh_.nVertices_*sizeof(bool));
		memset(visible,true,im_render->face_->mesh_.nVertices_*sizeof(bool));

		rs.estimateColor(shape,tex,faces,visible,noShadow,renderParams,colors);
		im_render->copyColors(colors);
		delete visible; delete noShadow;
		im_render->loadModel();
	//}
	cv::Mat refRGB = cv::Mat::zeros(colorIm.rows,colorIm.cols,CV_8UC3);
	cv::Mat refDepth = cv::Mat::zeros(colorIm.rows,colorIm.cols,CV_32F);
	im_render->render(r,t,_k[4],refRGB,refDepth);
	return refRGB;
}

void FaceServices2::testRenderWDepth(cv::Mat colorIm, cv::Mat alpha,cv::Mat vecR,cv::Mat vecT,cv::Mat exprW, cv::Mat &refRGB, cv::Mat &refDepth){
	char text[200];
	int M = 99;
	float renderParams[RENDER_PARAMS_COUNT];
	for (int i =0;i<3;i++)
		renderParams[i] = vecR.at<float>(i,0);
	for (int i =0;i<3;i++)
		renderParams[i+3] = vecT.at<float>(i,0);

	////cv::Mat shape = festimator.getShape(alpha,exprW);
	//cv::Mat faces_fill = festimator.getFaces_fill() - 1;
	// Ambient
	//renderParams[6] = 0.29225;
	//renderParams[7] = 0.29225;
	//renderParams[8] = 0.29225;
	renderParams[6] = 0.49225;
	renderParams[7] = 0.49225;
	renderParams[8] = 0.49225;
	// Diffuse
	renderParams[9] = 0.50754;
	renderParams[10] = 0.50754;
	renderParams[11] = 0.50754;
	// LIGHT
	renderParams[12] = 3.1415/4;
	renderParams[13] = 3.1415/4;
	// OTHERS
	renderParams[14] = 1;
	renderParams[15] = renderParams[16] = renderParams[17] = RENDER_PARAMS_GAIN_DEFAULT;
	renderParams[18] = renderParams[19] = renderParams[20] = RENDER_PARAMS_OFFSET_DEFAULT;

	Mat k_m(3,3,CV_32F,_k);
	//RenderServices rs;
	float* r = renderParams + RENDER_PARAMS_R;
	float* t = renderParams + RENDER_PARAMS_T;
	//printf("t %f %f %f\n",t[0],t[1],t[2]);

		shape = festimator.getShape2(alpha,exprW);
	//if (updateColor) {
		cv::Mat colors;
		tex = shape*0 + 128;

		im_render->copyShape(shape);
		bool* visible = new bool[im_render->face_->mesh_.nVertices_];
		bool* noShadow = new bool[im_render->face_->mesh_.nVertices_];
		memset(noShadow,true,im_render->face_->mesh_.nVertices_*sizeof(bool));
		memset(visible,true,im_render->face_->mesh_.nVertices_*sizeof(bool));

		rs.estimateColor(shape,tex,faces,visible,noShadow,renderParams,colors);
		im_render->copyColors(colors);
		delete visible; delete noShadow;
		im_render->loadModel();
	//}
	refRGB = cv::Mat::zeros(colorIm.rows,colorIm.cols,CV_8UC3);
	refDepth = cv::Mat::zeros(colorIm.rows,colorIm.cols,CV_32F);
	im_render->render(r,t,_k[4],refRGB,refDepth);
}

void FaceServices2::nextMotion(int &currFrame, cv::Mat &vecR, cv::Mat &vecT, cv::Mat &exprWeights){
    float stepYaw = 1;
    float PPI = 3.141592;
    int maxYaw = 70/stepYaw;
    float stepPitch = 1;
    int maxPitch = 45/stepPitch;

    int totalFrames = 4 * (maxYaw + maxPitch);
    currFrame = (currFrame + 1) % totalFrames;
    vecR = vecR*0 + 0.00001;
    // Rotate left
    if (currFrame < maxYaw) vecR.at<float>(1,0) = -currFrame*stepYaw * PPI/180;
    else if (currFrame < 2*maxYaw) vecR.at<float>(1,0) = -stepYaw * (2*maxYaw - currFrame) * PPI/180;
    // Rotate right
    else if (currFrame < 3*maxYaw) vecR.at<float>(1,0) = (currFrame-2*maxYaw) * stepYaw * PPI/180;
    else if (currFrame < 4*maxYaw) vecR.at<float>(1,0) = stepYaw * (4*maxYaw - currFrame) * PPI/180;

    // Rotate up
    else if (currFrame < 4*maxYaw + maxPitch) vecR.at<float>(0,0) = -(currFrame-4*maxYaw)*stepPitch * PPI/180;
    else if (currFrame < 4*maxYaw + 2*maxPitch) vecR.at<float>(0,0) = -stepPitch * (4*maxYaw+2*maxPitch - currFrame) * PPI/180;
    // Rotate right
    else if (currFrame < 4*maxYaw + 3*maxPitch) vecR.at<float>(0,0) = (currFrame-4*maxYaw-2*maxPitch) * stepPitch * PPI/180;
    else if (currFrame < 4*maxYaw + 4*maxPitch) vecR.at<float>(0,0) = stepPitch * (4*maxPitch+4*maxYaw - currFrame) * PPI/180;
    //std::cout << "vecR " << vecR.t() << std::endl;
}


text:
#define BOOST_FILESYSTEM_NO_DEPRECATED

#include "SrtmTile.h"

#include <boost/filesystem/operations.hpp>
#include <boost/interprocess/file_mapping.hpp>
#include <boost/interprocess/mapped_region.hpp>
#include <boost/move/unique_ptr.hpp>
#include <boost/thread.hpp>
#include <fmt/format.h>
#include <macgyver/Exception.h>

using FileMapping = boost::interprocess::file_mapping;
using MappedRegion = boost::interprocess::mapped_region;

using MutexType = boost::shared_mutex;
using ReadLock = boost::shared_lock<MutexType>;
using WriteLock = boost::unique_lock<MutexType>;
using UpgradeReadLock = boost::upgrade_lock<MutexType>;
using UpgradeWriteLock = boost::upgrade_to_unique_lock<MutexType>;

namespace Fmi
{
// ----------------------------------------------------------------------
/*!
 * \brief Validate a .hgt filename
 */
// ----------------------------------------------------------------------

bool SrtmTile::valid_path(const std::string &path)
{
  try
  {
    boost::filesystem::path p(path);

    // Avoid locale locs by not using regex here
    auto name = p.filename().string();
    return (name.size() == 1 + 2 + 1 + 3 + 4 && (name[0] == 'N' || name[0] == 'S') &&
            std::isdigit(name[1]) && std::isdigit(name[2]) && (name[3] == 'E' || name[3] == 'W') &&
            std::isdigit(name[4]) && std::isdigit(name[5]) && std::isdigit(name[6]) &&
            name[7] == '.' && name[8] == 'h' && name[9] == 'g' && name[10] == 't');
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

// ----------------------------------------------------------------------
/*!
 * \brief Validate .hgt file size
 *
 * The file size should be of the form N*N*2 for a square tile
 * covering a 1x1 degree tile.
 */
// ----------------------------------------------------------------------

bool SrtmTile::valid_size(const std::string &path)
{
  try
  {
    auto sz = boost::filesystem::file_size(path);
    auto n = static_cast<std::size_t>(sqrt(sz / 2.0));
    return (sz == 2 * n * n);
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

// ----------------------------------------------------------------------
/*!
 * \brief Implementation details
 */
// ----------------------------------------------------------------------

class SrtmTile::Impl
{
 public:
  explicit Impl(const std::string &path);
  const std::string &path() const { return itsPath; }
  std::size_t size() const { return itsSize; }
  int longitude() const { return itsLon; }
  int latitude() const { return itsLat; }
  int value(std::size_t i, std::size_t j);

 private:
  std::string itsPath;
  std::size_t itsSize;
  int itsLon;
  int itsLat;

  MutexType itsMutex;
  boost::movelib::unique_ptr<FileMapping> itsFileMapping;
  boost::movelib::unique_ptr<MappedRegion> itsMappedRegion;
};

// ----------------------------------------------------------------------
/*!
 * \brief Map the tile to memory
 */
// ----------------------------------------------------------------------

SrtmTile::Impl::Impl(const std::string &path) : itsPath(path)
{
  try
  {
    if (!valid_path(path))
      throw Fmi::Exception::Trace(BCP, "Not a valid filename for a .hgt file: '" + path + "'");
    if (!valid_size(path))
      throw Fmi::Exception::Trace(BCP,
                                  "Not a valid size of form 2*N*N for a .hgt file: '" + path + "'");

    auto sz = boost::filesystem::file_size(path);
    itsSize = static_cast<std::size_t>(sqrt(sz / 2.0));

    // Sample filename : S89E172.hgt

    boost::filesystem::path p(path);
    auto name = p.filename().string();

    int sign = (name[0] == 'S' ? -1 : 1);
    itsLat = sign * std::stoi(name.substr(1, 2));

    sign = (name[3] == 'W' ? -1 : 1);
    itsLon = sign * std::stoi(name.substr(4, 3));
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

// ----------------------------------------------------------------------
/*!
 * \brief Return the value at the given tile coordinate
 */
// ----------------------------------------------------------------------

int SrtmTile::Impl::value(std::size_t i, std::size_t j)
{
  try
  {
    // Sanity checks
    if (i >= itsSize || j >= itsSize)
      throw Fmi::Exception::Trace(
          BCP,
          fmt::format("SrtmFile indexes {},{} is out of range, size of tile is {}", i, j, itsSize));

    // We use lazy initialization to reduce core sizes in case
    // private mapped files are not disabled in core dumps.
    // We use private mapping intentionally, this class is expected
    // to be used in servers where no other process maps the same
    // files, and we wish to reserve shared mapping for core dumps.

    // We do not intend to write to the DEM files, but we use read_write
    // mode to get private mappings instead of shared ones.

    UpgradeReadLock readlock(itsMutex);

    if (!itsFileMapping)
    {
      UpgradeWriteLock writelock(readlock);
      if (!itsFileMapping)
      {
        // std::cout << "Mapping " << itsPath << std::endl;
        itsFileMapping = boost::movelib::make_unique<FileMapping>(itsPath.c_str(),
                                                                  boost::interprocess::read_only);

        itsMappedRegion = boost::movelib::make_unique<MappedRegion>(
            *itsFileMapping, boost::interprocess::read_only, 0, 2 * itsSize * itsSize);

        // We do not expect any normal access patterns, so disable prefetching
        itsMappedRegion->advise(boost::interprocess::mapped_region::advice_random);
      }
    }

    // Now the data is definitely available. Note: data runs from
    // north down, but we index if from bottom up

    auto *ptr = reinterpret_cast<std::int16_t *>(itsMappedRegion->get_address());
    std::int16_t big_endian = ptr[i + (itsSize - j - 1) * itsSize];
    std::int16_t little_endian = ((big_endian >> 8) & 0xff) + ((big_endian & 0xff) << 8);
    return little_endian;
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

// ----------------------------------------------------------------------
/*!
 * \brief The destructor is needed for the impl-idiom
 */
// ----------------------------------------------------------------------

SrtmTile::~SrtmTile() = default;

// ----------------------------------------------------------------------
/*!
 * \brief Read a SRTM tile
 */
// ----------------------------------------------------------------------

SrtmTile::SrtmTile(const std::string &path) : impl(new SrtmTile::Impl(path)) {}
// ----------------------------------------------------------------------
/*!
 * \brief Return the original path of the tile
 */
// ----------------------------------------------------------------------

const std::string &SrtmTile::path() const
{
  try
  {
    return impl->path();
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}
// ----------------------------------------------------------------------
/*!
 * \brief Return the size of the tile
 */
// ----------------------------------------------------------------------

std::size_t SrtmTile::size() const
{
  try
  {
    return impl->size();
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}
// ----------------------------------------------------------------------
/*!
 * \brief Return the lower left corner longitude of the tile
 */
// ----------------------------------------------------------------------

int SrtmTile::longitude() const
{
  try
  {
    return impl->longitude();
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}
// ----------------------------------------------------------------------
/*!
 * \brief Return the lower left corner latitude of the tile
 */
// ----------------------------------------------------------------------

int SrtmTile::latitude() const
{
  try
  {
    return impl->latitude();
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}
// ----------------------------------------------------------------------
/*!
 * \brief Return the value at the given tile coordinate
 *
 * Lack of interpolation support is intentional, it is provided
 * at the SrtmMatrix level where adjacent tiles are available
 */
// ----------------------------------------------------------------------

int SrtmTile::value(std::size_t i, std::size_t j) const
{
  try
  {
    return impl->value(i, j);
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}
}  // namespace Fmi

ext: cpp

max_stars_repo_name: fmidev/smartmet-library-gis

max_issues_repo_name: fmidev/smartmet-library-gis

max_forks_repo_name: fmidev/smartmet-library-gis


text:
#include <boost/python/class.hpp>
#include <boost/python/args.hpp>
#include <boost/python/return_value_policy.hpp>
#include <boost/python/copy_const_reference.hpp>
#include <boost/python/return_internal_reference.hpp>
#include <boost/python/with_custodian_and_ward.hpp>
#include <cctbx/sgtbx/wyckoff.h>
#include <scitbx/boost_python/is_polymorphic_workaround.h>

SCITBX_BOOST_IS_POLYMORPHIC_WORKAROUND(cctbx::sgtbx::wyckoff::mapping)

namespace cctbx { namespace sgtbx { namespace boost_python {

namespace {

  struct wyckoff_position_wrappers
  {
    typedef wyckoff::position w_t;

    static void
    wrap()
    {
      using namespace boost::python;
      typedef return_value_policy<copy_const_reference> ccr;
      class_<w_t>("wyckoff_position", no_init)
        .def("multiplicity", &w_t::multiplicity)
        .def("letter", &w_t::letter)
        .def("special_op", &w_t::special_op, ccr())
        .def("point_group_type", &w_t::point_group_type)
        .def("unique_ops", &w_t::unique_ops, (arg("space_group")))
      ;
    }
  };

  struct wyckoff_mapping_wrappers
  {
    typedef wyckoff::mapping w_t;

    static void
    wrap()
    {
      using namespace boost::python;
      typedef return_value_policy<copy_const_reference> ccr;
      typedef return_internal_reference<> rir;
      class_<w_t>("wyckoff_mapping", no_init)
        .def("unit_cell", &w_t::unit_cell, rir())
        .def("original_site", &w_t::original_site)
        .def("position", &w_t::position, rir())
        .def("sym_op", &w_t::sym_op, ccr())
        .def("representative_site", &w_t::representative_site)
        .def("exact_site", &w_t::exact_site)
        .def("distance_moved", &w_t::distance_moved)
        .def("special_op", &w_t::special_op)
      ;
    }
  };

  struct wyckoff_table_wrappers
  {
    typedef wyckoff::table w_t;

    static void
    wrap()
    {
      using namespace boost::python;
      typedef return_value_policy<copy_const_reference> ccr;
      typedef return_internal_reference<> rir;
      class_<w_t>("wyckoff_table", no_init)
        .def(init<space_group_type const&>((arg("space_group_type"))))
        .def("space_group_type", &w_t::space_group_type, rir())
        .def("size", &w_t::size)
        .def("position",
           (wyckoff::position const&(w_t::*)(std::size_t) const)
           &w_t::position, rir())
        .def("position",
           (wyckoff::position const&(w_t::*)(char) const)
           &w_t::position, (arg("letter")), rir())
        .def("lookup_index", &w_t::lookup_index, (arg("letter")))
        .def("mapping",
           (wyckoff::mapping(w_t::*)(site_symmetry const&) const)
           &w_t::mapping, (arg("site_symmetry")),
           with_custodian_and_ward_postcall<0,1>())
        .def("mapping",
          (wyckoff::mapping(w_t::*)(
            uctbx::unit_cell const&,
            fractional<> const&,
            double) const) &w_t::mapping, (
              arg("unit_cell"),
              arg("original_site"),
              arg("special_position_radius")=0.5),
                with_custodian_and_ward_postcall<0,1>())
      ;
    }
  };

} // namespace <anoymous>

  void wrap_wyckoff()
  {
    wyckoff_position_wrappers::wrap();
    wyckoff_mapping_wrappers::wrap();
    wyckoff_table_wrappers::wrap();
  }

}}} // namespace cctbx::sgtbx::boost_python


text:
// Serializable Map - SerializableMap.hpp (Wrapper)
#ifndef SERIALIZABLE_MAP_HPP
#define SERIALIZABLE_MAP_HPP

#include <map>
#include <utility>
#include <string>
#include <fstream>
#include <iostream>
#ifndef __CINT__
#include <boost/noncopyable.hpp>
#include <boost/optional.hpp>
#else
namespace boost;
class boost::noncopyable;
#endif /* __CINT __ */

#include "../HeaderFiles/DBHandler.h"

namespace DB
{

namespace DATA
{

template <typename KEY_TYPE, typename VALUE_TYPE >
class SerializableMap : public boost::noncopyable
{
protected:
  DB::SERVICES::DBHandler &m_dbHandlerInstance;
  std::map<KEY_TYPE, VALUE_TYPE > m_map;
  std::string m_fileNameToSaveMap;

public:
  SerializableMap(std::string p_fileNameToSaveMap, SERVICES::DBHandler &l_dBHandler = SERVICES::DBHandler::getInstance());
  virtual ~SerializableMap(void);

  virtual void clearMap(void);
  virtual void eraseElementFromMap(KEY_TYPE p_key);
  virtual void showMap(void) const = 0;
  virtual boost::optional<VALUE_TYPE > valueFromMap(KEY_TYPE p_key) = 0;
  virtual size_t sizeOfMap(void) const;
  virtual void saveMapToFile(std::string p_fileName) = 0;
  virtual bool removeFileWithMap(void);
  virtual boost::optional<std::pair<KEY_TYPE, VALUE_TYPE > > firstElement(void) const;
  virtual bool empty(void) const;
};

template <typename KEY_TYPE, typename VALUE_TYPE >
SerializableMap<KEY_TYPE, VALUE_TYPE >::SerializableMap(std::string p_fileNameToSaveMap,
                                                        SERVICES::DBHandler& l_dBHandler) :
                                                                                m_dbHandlerInstance(l_dBHandler),
                                                                                m_fileNameToSaveMap(p_fileNameToSaveMap)
{
}

template <typename KEY_TYPE, typename VALUE_TYPE >
SerializableMap<KEY_TYPE, VALUE_TYPE >::~SerializableMap()
{
}

template <typename KEY_TYPE, typename VALUE_TYPE >
void SerializableMap<KEY_TYPE, VALUE_TYPE >::clearMap()
{
  m_map.clear();
}

template <typename KEY_TYPE, typename VALUE_TYPE >
void SerializableMap<KEY_TYPE, VALUE_TYPE >::eraseElementFromMap(KEY_TYPE p_key)
{
  typedef typename std::map<KEY_TYPE, VALUE_TYPE >::iterator IteratorType;
  IteratorType l_iter = m_map.find(p_key);

  if(l_iter != m_map.end())
  {
    m_map.erase(l_iter);
  }
}

template <typename KEY_TYPE, typename VALUE_TYPE >
size_t SerializableMap<KEY_TYPE, VALUE_TYPE >::sizeOfMap() const
{
  return m_map.size();
}

template <typename KEY_TYPE, typename VALUE_TYPE >
bool SerializableMap<KEY_TYPE, VALUE_TYPE >::removeFileWithMap()
{
  if(remove(m_fileNameToSaveMap.c_str()) != 0)
  {
    std::string l_errorMessage = "Error while deleting file : " + m_fileNameToSaveMap;

    perror(l_errorMessage.c_str());
  }
  return 0;
}

template <typename KEY_TYPE, typename VALUE_TYPE >
boost::optional<std::pair<KEY_TYPE, VALUE_TYPE > > SerializableMap<KEY_TYPE, VALUE_TYPE >::firstElement() const
{
  if(!m_map.empty())
  {
    return std::make_pair(m_map.begin()->first, m_map.begin()->second);
  }
  return boost::none;
}

template <typename KEY_TYPE, typename VALUE_TYPE >
bool SerializableMap<KEY_TYPE, VALUE_TYPE >::empty() const
{
  return m_map.empty();
}

} // namespace DATA

} // namespace DB

#endif // SERIALIZABLE_MAP_HPP

ext: hpp

max_stars_repo_name: kamilrakoczy/j-pet-framework

max_issues_repo_name: kamilrakoczy/j-pet-framework

max_forks_repo_name: kamilrakoczy/j-pet-framework


text:
//
// lager - library for functional interactive c++ programs
// Copyright (C) 2017 Juan Pedro Bolivar Puente
//
// This file is part of lager.
//
// lager is free software: you can redistribute it and/or modify
// it under the terms of the MIT License, as detailed in the LICENSE
// file located at the root of this source code distribution,
// or here: <https://github.com/arximboldi/lager/blob/master/LICENSE>
//

#pragma once

#include <boost/preprocessor/punctuation/remove_parens.hpp>
#include <boost/preprocessor/seq/for_each.hpp>

#include <boost/functional/hash.hpp>

#include <functional>

#define LAGER_DERIVE_IMPL_HASH_ITER__(r__, data__, elem__)                     \
    ::boost::hash_combine(seed, x.elem__);

#define LAGER_DERIVE_IMPL_HASH(r__, ns__, name__, members__)                   \
    namespace std {                                                            \
    template <>                                                                \
    struct hash<ns__::name__>                                                  \
    {                                                                          \
        std::size_t operator()(const ns__::name__& x) const                    \
        {                                                                      \
            auto seed = std::size_t{};                                         \
            BOOST_PP_SEQ_FOR_EACH_R(                                           \
                r__, LAGER_DERIVE_IMPL_HASH_ITER__, _, members__)              \
            return seed;                                                       \
        }                                                                      \
    };                                                                         \
    }                                                                          \
    //

#define LAGER_DERIVE_TEMPLATE_IMPL_HASH(r__, ns__, tpl__, name__, members__)   \
    namespace std {                                                            \
    template <BOOST_PP_REMOVE_PARENS(tpl__)>                                   \
    struct hash<ns__::BOOST_PP_REMOVE_PARENS(name__)>                          \
    {                                                                          \
        std::size_t operator()(const ns__::BOOST_PP_REMOVE_PARENS(name__) &    \
                               x) const                                        \
        {                                                                      \
            auto seed = std::size_t{};                                         \
            BOOST_PP_SEQ_FOR_EACH_R(                                           \
                r__, LAGER_DERIVE_IMPL_HASH_ITER__, _, members__)              \
            return seed;                                                       \
        }                                                                      \
    };                                                                         \
    }                                                                          \
    //

#define LAGER_DERIVE_NESTED_IMPL_HASH(r__, name__, members__)                  \
    static_assert(false, "can't implement HASH for LAGER_DERIVE_NESTED, sorry!")

max_stars_repo_name: gomez-addams/lager

max_issues_repo_name: gomez-addams/lager

max_forks_repo_name: gomez-addams/lager


text:
//==============================================================================
//         Copyright 2003 - 2011 LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2011 LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_SIMD_TOOLBOX_ARITHMETIC_FUNCTIONS_SIMD_COMMON_FMA_HPP_INCLUDED
#define BOOST_SIMD_TOOLBOX_ARITHMETIC_FUNCTIONS_SIMD_COMMON_FMA_HPP_INCLUDED
#include <boost/simd/toolbox/arithmetic/functions/fma.hpp>
#include <boost/simd/include/functions/simd/multiplies.hpp>
#include <boost/simd/include/functions/simd/plus.hpp>
// TO DO define a proper fma policy

namespace boost { namespace simd { namespace ext
{
  BOOST_SIMD_FUNCTOR_IMPLEMENTATION(boost::simd::tag::fma_, tag::cpu_,
                      (A0)(X),
                      ((simd_<arithmetic_<A0>,X>))
                      ((simd_<arithmetic_<A0>,X>))
                      ((simd_<arithmetic_<A0>,X>))
                     )
  {
    typedef A0 result_type;
    BOOST_SIMD_FUNCTOR_CALL_REPEAT(3)
    {
      return a0*a1+a2;
    }
  };
} } }
#endif


text:
/* vim:set ts=3 sw=3 sts=3 et: */

/***********************************************************************

Copyright (c) 2012 Marcus Holland-Moritz

Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:

The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

***********************************************************************/

#include <iostream>
#include <stdexcept>
#include <syslog.h>

#include <boost/date_time/posix_time/posix_time.hpp>

#include "log.h"

namespace aird {

log_level::type log_level::str2level(const std::string& str)
{
   if (str == "trace") return TRACE;
   if (str == "debug") return DEBUG;
   if (str == "info") return INFO;
   if (str == "warn") return WARN;
   if (str == "error") return ERROR;
   if (str == "fatal") return FATAL;
   throw std::runtime_error("invalid logging level: " + str);
}

const char *log_level::level2str(type level)
{
   switch (level)
   {
      case TRACE: return "trace";
      case DEBUG: return "debug";
      case INFO:  return "info";
      case WARN:  return "warn";
      case ERROR: return "error";
      case FATAL: return "fatal";
      default: return "***";
   }
}

appender::~appender()
{
}

root_logger::root_logger(log_level::type level)
   : m_level(level)
{
}

void root_logger::log(log_level::type level, const std::string& msg) const
{
   if (level >= m_level)
   {
      for (std::vector< boost::shared_ptr<appender> >::const_iterator it = m_app.begin(); it != m_app.end(); ++it)
      {
         (*it)->log(level, msg);
      }
   }
}

void root_logger::add_appender(boost::shared_ptr<appender> app)
{
   m_app.push_back(app);
}

console_appender::console_appender(log_level::type thresh)
   : m_thresh(thresh)
{
}

void console_appender::log(log_level::type level, const std::string& msg) const
{
   if (level >= m_thresh)
   {
      boost::posix_time::ptime now = boost::posix_time::microsec_clock::local_time();
      std::cerr << to_simple_string(now) << " <" << log_level::level2str(level) << "> " << msg << std::endl;
   }
}

syslog_appender::syslog_appender(const std::string& ident, log_level::type thresh)
   : m_thresh(thresh)
{
   ::openlog(ident.c_str(), 0, LOG_DAEMON);
}

syslog_appender::~syslog_appender()
{
   ::closelog();
}

void syslog_appender::log(log_level::type level, const std::string& msg) const
{
   if (level >= m_thresh)
   {
      int priority;

      switch (level)
      {
         case log_level::INFO:  priority = LOG_INFO; break;
         case log_level::WARN:  priority = LOG_WARNING; break;
         case log_level::ERROR: priority = LOG_ERR; break;
         case log_level::FATAL: priority = LOG_CRIT; break;
         default:               priority = LOG_DEBUG; break;
      }

      ::syslog(priority, "%s", msg.c_str());
   }
}

}

max_stars_repo_name: mhx/aird

max_issues_repo_name: mhx/aird

max_forks_repo_name: mhx/aird


text:
#pragma once
#include "process.hpp"
#include "source_diff.hpp"
#include "source_spellcheck.hpp"
#include "tooltips.hpp"
#include <boost/filesystem.hpp>
#include <boost/property_tree/xml_parser.hpp>
#include <limits>
#include <set>
#include <string>
#include <tuple>
#include <unordered_map>
#include <vector>

namespace Source {
  /// Workaround for buggy Gsv::LanguageManager::get_default()
  class LanguageManager {
  public:
    static Glib::RefPtr<Gsv::LanguageManager> get_default();
  };
  /// Workaround for buggy Gsv::StyleSchemeManager::get_default()
  class StyleSchemeManager {
  public:
    static Glib::RefPtr<Gsv::StyleSchemeManager> get_default();
  };

  Glib::RefPtr<Gsv::Language> guess_language(const boost::filesystem::path &file_path);

  class Offset {
  public:
    Offset() = default;
    Offset(unsigned line, unsigned index, boost::filesystem::path file_path_ = {}) : line(line), index(index), file_path(std::move(file_path_)) {}
    operator bool() const { return !file_path.empty(); }
    bool operator<(const Offset &rhs) const {
      return file_path < rhs.file_path ||
             (file_path == rhs.file_path && (line < rhs.line ||
                                             (line == rhs.line && index < rhs.index)));
    }
    bool operator==(const Offset &rhs) const { return (file_path == rhs.file_path && line == rhs.line && index == rhs.index); }

    unsigned line = 0;
    unsigned index = 0;
    boost::filesystem::path file_path;
  };

  class FixIt {
  public:
    enum class Type {
      insert,
      replace,
      erase
    };

    FixIt(std::string source_, std::string path_, std::pair<Offset, Offset> offsets_);

    std::string string(BaseView &view);

    Type type;
    std::string source;
    std::string path;
    std::pair<Offset, Offset> offsets;

    bool operator<(const FixIt &rhs) const {

      return type < rhs.type ||
             (type == rhs.type && (source < rhs.source ||
                                   (source == rhs.source && (path < rhs.path ||
                                                             (path == rhs.path && offsets < rhs.offsets)))));
    }
    bool operator==(const FixIt &rhs) const {
      return type == rhs.type && source == rhs.source && path == rhs.path && offsets == rhs.offsets;
    }
  };

  class View : public SpellCheckView, public DiffView {
  public:
    static std::set<View *> non_deleted_views;
    static std::set<View *> views;

    View(const boost::filesystem::path &file_path, const Glib::RefPtr<Gsv::Language> &language, bool is_generic_view = false);
    ~View() override;

    bool save() override;
    void configure() override;

    std::function<void()> non_interactive_completion;
    std::function<void(bool, bool)> format_style;
    std::function<Offset()> get_declaration_location;
    std::function<Offset()> get_type_declaration_location;
    std::function<std::vector<Offset>()> get_implementation_locations;
    std::function<std::vector<Offset>()> get_declaration_or_implementation_locations;
    std::function<std::vector<std::pair<Offset, std::string>>()> get_usages;
    std::function<std::string()> get_method;
    std::function<std::vector<std::pair<Offset, std::string>>()> get_methods;
    std::function<std::vector<std::string>()> get_token_data;
    std::function<std::string()> get_token_spelling;
    std::function<void(const std::string &text)> rename_similar_tokens;
    std::function<void()> goto_next_diagnostic;
    std::function<std::vector<FixIt>()> get_fix_its;
    std::function<void()> toggle_comments;
    std::function<std::tuple<Source::Offset, std::string, size_t>()> get_documentation_template;
    std::function<void(int)> toggle_breakpoint;

    void hide_tooltips() override;
    void hide_dialogs() override;

    void scroll_to_cursor_delayed(bool center, bool show_tooltips) override;

    void extend_selection();
    void shrink_selection();

    void show_or_hide(); /// Show or hide text selection

    bool soft_reparse_needed = false;
    bool full_reparse_needed = false;
    virtual void soft_reparse(bool delayed = false) { soft_reparse_needed = false; }
    virtual void full_reparse() { full_reparse_needed = false; }

    static std::unique_ptr<TinyProcessLib::Process> prettier_background_process;

  protected:
    std::atomic<bool> parsed = {true};
    Tooltips diagnostic_tooltips;
    Tooltips type_tooltips;
    sigc::connection delayed_tooltips_connection;

    Glib::RefPtr<Gtk::TextTag> similar_symbol_tag;
    sigc::connection delayed_tag_similar_symbols_connection;
    virtual void apply_similar_symbol_tag() {}
    bool similar_symbol_tag_applied = false;
    Glib::RefPtr<Gtk::TextTag> clickable_tag;
    sigc::connection delayed_tag_clickable_connection;
    virtual void apply_clickable_tag(const Gtk::TextIter &iter) {}
    bool clickable_tag_applied = false;

    Glib::RefPtr<Gtk::TextTag> hide_tag;

    virtual void show_diagnostic_tooltips(const Gdk::Rectangle &rectangle) { diagnostic_tooltips.show(rectangle); }
    void add_diagnostic_tooltip(const Gtk::TextIter &start, const Gtk::TextIter &end, bool error, std::function<void(Tooltip &)> &&set_buffer);
    void clear_diagnostic_tooltips();
    std::set<int> diagnostic_offsets;
    void place_cursor_at_next_diagnostic();
    virtual void show_type_tooltips(const Gdk::Rectangle &rectangle) {}
    gdouble on_motion_last_x = 0.0;
    gdouble on_motion_last_y = 0.0;

    std::vector<FixIt> fix_its;

    /// Returns true if code iter (not part of comment or string, and not whitespace) is found.
    /// Iter will not be moved if iter is already a code iter.
    bool backward_to_code(Gtk::TextIter &iter);
    /// Returns true if code iter (not part of comment or string, and not whitespace) is found.
    /// Iter will not be moved if iter is already a code iter.
    bool forward_to_code(Gtk::TextIter &iter);
    /// Iter will not be moved if iter is already a code iter (not part of comment or string, and not whitespace) or at line start
    bool backward_to_code_or_line_start(Gtk::TextIter &iter);
    /// Iter will not be moved if iter is already a code iter (not part of comment or string, and not whitespace) or at line end
    bool forward_to_code_or_line_end(Gtk::TextIter &iter);
    /// If closing bracket is found, continues until the open bracket.
    /// Returns if open bracket is found that has no corresponding closing bracket.
    /// Else, return at start of line.
    Gtk::TextIter get_start_of_expression(Gtk::TextIter iter);
    /// Iter will not be moved if iter is already at close symbol.
    bool find_close_symbol_forward(Gtk::TextIter iter, Gtk::TextIter &found_iter, unsigned int positive_char, unsigned int negative_char);
    /// Iter will not be moved if iter is already at open symbol.
    bool find_open_symbol_backward(Gtk::TextIter iter, Gtk::TextIter &found_iter, unsigned int positive_char, unsigned int negative_char);
    long symbol_count(Gtk::TextIter iter, unsigned int positive_char, unsigned int negative_char = std::numeric_limits<unsigned int>::max());
    bool is_templated_function(Gtk::TextIter iter, Gtk::TextIter &parenthesis_end_iter);
    /// If insert is at an possible argument. Also based on last key press.
    bool is_possible_argument();

    bool on_key_press_event(GdkEventKey *event) override;
    bool on_key_press_event_basic(GdkEventKey *event);
    bool on_key_press_event_bracket_language(GdkEventKey *event);
    bool on_key_press_event_smart_brackets(GdkEventKey *event);
    bool on_key_press_event_smart_inserts(GdkEventKey *event);
    bool on_button_press_event(GdkEventButton *event) override;

    bool interactive_completion = true;

  private:
    void setup_signals();
    void setup_format_style(bool is_generic_view);

    Gsv::DrawSpacesFlags parse_show_whitespace_characters(const std::string &text);

    Gsv::GutterRendererText *line_renderer;

    bool use_fixed_continuation_indenting = true;
    guint previous_non_modifier_keyval = 0;

    bool keep_previous_extended_selections = false;
    std::vector<std::pair<Gtk::TextIter, Gtk::TextIter>> previous_extended_selections;
  };
} // namespace Source

max_stars_repo_name: jlangvand/jucipp

max_issues_repo_name: jlangvand/jucipp

max_forks_repo_name: jlangvand/jucipp


text:
#include <bio/fasta.h>
USING_BIO_NS

#include <boost/test/unit_test.hpp>
#include <boost/progress.hpp>
using namespace boost;
using boost::unit_test::test_suite;

#include <fstream>
using namespace std;

#define FASTA_FILE "C:\\data\\ensembl\\Homo_sapiens.NCBI35.may.dna.chromosome.21.fa"

void check_parse_fasta()
{
	cout << "******* check_parse_fasta(): Parsing fasta file " FASTA_FILE << endl;

	ifstream file(FASTA_FILE);
	BOOST_CHECK(file.good());

	stringstream sequence;
	progress_timer timer;
	parse_fasta(file, sequence);
}

void check_parse_fasta_2()
{
	const std::string file = "c:\\data\\Helene\\MyoD_DRR.fa";
	cout << "******* check_parse_fasta_2(): Parsing fasta file " << file << "\n";

	ifstream stream(file.c_str());
	BOOST_CHECK(stream.good());

	fasta_file_map_t sequences;
	parse_fasta_2(stream, sequences);
}

void
register_fasta_tests(boost::unit_test::test_suite * test)
{
	//test->add( BOOST_TEST_CASE( &check_parse_fasta ), 0);
	test->add( BOOST_TEST_CASE( &check_parse_fasta_2 ), 0);
}

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
// Distributed under the MIT License.
// See LICENSE.txt for details.

#pragma once

#include <array>
#include <boost/functional/hash.hpp>  // IWYU pragma: keep
#include <cstddef>
#include <functional>
#include <limits>
#include <unordered_map>
#include <utility>
#include <vector>

#include "DataStructures/ApplyMatrices.hpp"
#include "DataStructures/DataVector.hpp"
#include "DataStructures/Index.hpp"
#include "DataStructures/Matrix.hpp"
#include "DataStructures/Tags.hpp"       // IWYU pragma: keep
#include "DataStructures/Variables.hpp"  // IWYU pragma: keep
#include "Domain/Structure/Direction.hpp"
#include "Domain/Structure/DirectionMap.hpp"
#include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
#include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
#include "Evolution/DiscontinuousGalerkin/Limiters/WenoOscillationIndicator.hpp"
#include "NumericalAlgorithms/LinearOperators/MeanValue.hpp"
#include "Utilities/Algorithm.hpp"
#include "Utilities/EqualWithinRoundoff.hpp"
#include "Utilities/ErrorHandling/Assert.hpp"
#include "Utilities/Gsl.hpp"

/// \cond
template <size_t VolumeDim>
class Element;
template <size_t VolumeDim>
class ElementId;
template <size_t>
class Mesh;
/// \endcond

namespace Limiters::Weno_detail {

// Caching class that holds various precomputed terms used in the constrained-
// fit algebra on each element.
//
// The terms to precompute and cache depend on the configuration of neighbors to
// the element (how many neighbors, in which directions, of what resolution).
// Each instance of the caching class represents one particular neighbor
// configuration, and typical simulations will create many instances of
// the caching class: one instance for each neighbor-element configuration
// present in the computational domain.
//
// Because the current implementation of the HWENO fitting makes the simplifying
// restrictions of no h/p-refinement, the structure of the caching is also
// simplified. In particular, with no h/p-refinement, the allowable neighbor
// configurations satisfy,
// 1. the element has at most one neighbor per dimension, AND
// 2. the mesh on every neighbor is the same as on the element.
// With these restrictions, the number of independent configurations to be
// cached is greatly reduced. Because an element has 2*VolumeDim boundaries,
// it has 2^(2*VolumeDim) possible configurations of internal/external
// boundaries, and therefore there are 2^(2*VolumeDim) configurations to cache.
// Different elements with the same configuration of internal/external
// boundaries vs. direction can share the same caching-class instance.
//
// Each instance of the caching class holds several terms, some of which also
// depend on the neighbor configuration. The restriction of no h/p-refinement
// therefore simplifies each cache instance, as well as reducing the necessary
// number of instances.
//
// The most complicated term to cache is the A^{-1} matrix. The complexity
// arises because the matrix can take many values depending on (runtime) choices
// made for each individual HWNEO fit: which element is the primary neighbor,
// and which element(s) are the excluded neighbors.
// Fits with one excluded neighbor are overwhelmingly the most likely, and the
// cache is designed for this particular case. Fits with no excluded neighbors
// can arise in elements with only one neighbor (always the primary neighbor);
// these cases are also handled. However, the very rare case in which more than
// one neighbor is excluded is not handled by the cache; the caller must compute
// A^{-1} from scratch if this scenario arises.
template <size_t VolumeDim>
class ConstrainedFitCache {
 public:
  ConstrainedFitCache(const Element<VolumeDim>& element,
                      const Mesh<VolumeDim>& mesh) noexcept;

  // Valid calls must satisfy these constraints on directions_to_exclude:
  // - the vector is empty, OR
  // - the vector contains one element, and this element is a direction
  //   different from the direction to the primary neighbor.
  // The very rare case where more than one neighbor is excluded from the fit is
  // not cached. In this case, A^{-1} must be computed from scratch.
  const Matrix& retrieve_inverse_a_matrix(
      const Direction<VolumeDim>& primary_direction,
      const std::vector<Direction<VolumeDim>>& directions_to_exclude) const
      noexcept;

  DataVector quadrature_weights;
  DirectionMap<VolumeDim, Matrix> interpolation_matrices;
  DirectionMap<VolumeDim, DataVector>
      quadrature_weights_dot_interpolation_matrices;
  // The many possible values of A^{-1} are stored in a map of maps. The outer
  // map indexes over the primary neighbor, and the inner map indexes over the
  // excluded neighbor.
  // This data structure is not perfect: configurations with only one neighbor
  // (always the primary neighbor) lead to no excluded neighbors, and there is
  // no natural place to hold A^{-1} in the maps. For this case, we simply store
  // the data in the normally-nonsensical slot where
  // excluded_neighbor == primary_neighbor.
  DirectionMap<VolumeDim, DirectionMap<VolumeDim, Matrix>> inverse_a_matrices;
};

// Return the appropriate cache for the given element and mesh.
template <size_t VolumeDim>
const ConstrainedFitCache<VolumeDim>& constrained_fit_cache(
    const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh) noexcept;

// Compute the inverse of the matrix A_st for the constrained fit.
// See the documentation of `hweno_modified_neighbor_solution`.
template <size_t VolumeDim>
Matrix inverse_a_matrix(
    const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
    const DataVector& quadrature_weights,
    const DirectionMap<VolumeDim, Matrix>& interpolation_matrices,
    const DirectionMap<VolumeDim, DataVector>&
        quadrature_weights_dot_interpolation_matrices,
    const Direction<VolumeDim>& primary_direction,
    const std::vector<Direction<VolumeDim>>& directions_to_exclude) noexcept;

// Not all secondary neighbors (i.e., neighbors that are not the primary) are
// included in the HWENO constrained fit. In particular, for the tensor
// component specified by `Tag` and `tensor_index`, the secondary neighbor whose
// mean is most different from the troubled cell's mean is excluded.
//
// Zhu2016 indicate that when multiple secondary neighbors share the property of
// being "most different" from the troubled cell, then all of these secondary
// neighbors should be excluded from the minimization. The probability of this
// occuring is exceedingly low, but we handle it anyway. We return the excluded
// secondary neighbors in a vector.
//
// Note that if there is only one neighbor, it is the primary neighbor, and so
// there are no secondary neighbors to exclude. We return an empty vector. This
// scenario can arise in various test cases, but is unlikely to arise in science
// cases.
template <typename Tag, size_t VolumeDim, typename Package>
std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>
secondary_neighbors_to_exclude_from_fit(
    const double local_mean, const size_t tensor_index,
    const std::unordered_map<
        std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
        boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
        neighbor_data,
    const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
        primary_neighbor) noexcept {
  // Rare case: with only one neighbor, there is no secondary to exclude
  if (UNLIKELY(neighbor_data.size() == 1)) {
    return std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>{};
  }

  // Identify element with maximum mean difference
  const auto mean_difference =
      [&tensor_index, &local_mean](
          const std::pair<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>,
                          Package>& neighbor_and_data) noexcept {
        return fabs(get<::Tags::Mean<Tag>>(
                        neighbor_and_data.second.means)[tensor_index] -
                    local_mean);
      };

  double max_difference = std::numeric_limits<double>::lowest();
  std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>
      neighbor_max_difference{};
  for (const auto& neighbor_and_data : neighbor_data) {
    const auto& neighbor = neighbor_and_data.first;
    if (neighbor == primary_neighbor) {
      continue;
    }
    const double difference = mean_difference(neighbor_and_data);
    if (difference > max_difference) {
      max_difference = difference;
      neighbor_max_difference = neighbor;
    }
  }

  std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>
      neighbors_to_exclude{{neighbor_max_difference}};

  // See if other elements share this maximum mean difference. This loop should
  // only rarely find other neighbors with the same maximal mean difference to
  // add to the vector, so it will usually not change the vector.
  for (const auto& neighbor_and_data : neighbor_data) {
    const auto& neighbor = neighbor_and_data.first;
    if (neighbor == primary_neighbor or neighbor == neighbor_max_difference) {
      continue;
    }
    const double difference = mean_difference(neighbor_and_data);
    if (UNLIKELY(equal_within_roundoff(difference, max_difference))) {
      neighbors_to_exclude.push_back(neighbor);
    }
  }

  ASSERT(not alg::found(neighbors_to_exclude, primary_neighbor),
         "Logical inconsistency: trying to exclude the primary neighbor.");

  return neighbors_to_exclude;
}

// Compute the vector b_s for the constrained fit. For details, see the
// documentation of `hweno_modified_neighbor_solution` below.
template <typename Tag, size_t VolumeDim, typename Package>
DataVector b_vector(
    const size_t tensor_index, const Mesh<VolumeDim>& mesh,
    const DataVector& quadrature_weights,
    const DirectionMap<VolumeDim, Matrix>& interpolation_matrices,
    const DirectionMap<VolumeDim, DataVector>&
        quadrature_weights_dot_interpolation_matrices,
    const std::unordered_map<
        std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
        boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
        neighbor_data,
    const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
        primary_neighbor,
    const std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>&
        neighbors_to_exclude) noexcept {
  const size_t number_of_grid_points = mesh.number_of_grid_points();
  DataVector b(number_of_grid_points, 0.);

  for (const auto& neighbor_and_data : neighbor_data) {
    // Generally, neighbors_to_exclude contains just one neighbor to exclude,
    // and the loop over neighbor_data will hit this (in 3D) roughly 1/6 times.
    // So the condition is somewhat, but not overwhelmingly, unlikely:
    if (UNLIKELY(alg::found(neighbors_to_exclude, neighbor_and_data.first))) {
      continue;
    }

    const auto& direction = neighbor_and_data.first.first;
    ASSERT(interpolation_matrices.contains(direction),
           "interpolation_matrices does not contain key: " << direction);
    ASSERT(
        quadrature_weights_dot_interpolation_matrices.contains(direction),
        "quadrature_weights_dot_interpolation_matrices does not contain key: "
            << direction);

    const auto& neighbor_mesh = mesh;
    const auto& neighbor_quadrature_weights = quadrature_weights;
    const auto& interpolation_matrix = interpolation_matrices.at(direction);
    const auto& quadrature_weights_dot_interpolation_matrix =
        quadrature_weights_dot_interpolation_matrices.at(direction);

    const auto& neighbor_tensor_component =
        get<Tag>(neighbor_and_data.second.volume_data)[tensor_index];

    // Add terms from the primary neighbor
    if (neighbor_and_data.first == primary_neighbor) {
      for (size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
        for (size_t s = 0; s < number_of_grid_points; ++s) {
          b[s] += neighbor_tensor_component[r] *
                  neighbor_quadrature_weights[r] * interpolation_matrix(r, s);
        }
      }
    }
    // Add terms from the secondary neighbors
    else {
      const double quadrature_weights_dot_u = [&]() noexcept {
        double result = 0.;
        for (size_t r = 0; r < neighbor_mesh.number_of_grid_points(); ++r) {
          result +=
              neighbor_tensor_component[r] * neighbor_quadrature_weights[r];
        }
        return result;
      }();
      b += quadrature_weights_dot_u *
           quadrature_weights_dot_interpolation_matrix;
    }
  }
  return b;
}

// Solve the constrained fit problem that gives the HWENO modified solution,
// for one particular tensor component. For details, see documentation of
// `hweno_modified_neighbor_solution` below.
template <typename Tag, size_t VolumeDim, typename Package>
void solve_constrained_fit(
    const gsl::not_null<DataVector*> constrained_fit_result,
    const DataVector& u, const size_t tensor_index, const Mesh<VolumeDim>& mesh,
    const Element<VolumeDim>& element,
    const std::unordered_map<
        std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, Package,
        boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
        neighbor_data,
    const std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>&
        primary_neighbor,
    const std::vector<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>&
        neighbors_to_exclude) noexcept {
  ASSERT(not alg::found(neighbors_to_exclude, primary_neighbor),
         "Logical inconsistency: trying to exclude the primary neighbor.");
  ASSERT(not neighbors_to_exclude.empty() or neighbor_data.size() == 1,
         "The HWENO constrained fit algorithm expects at least one neighbor \n"
         "to exclude from the fit (unless if the element has a single \n"
         "neighbor, which would automatically be the primary neighbor).");

  // Get the cache of linear algebra quantities for this element
  const ConstrainedFitCache<VolumeDim>& cache =
      constrained_fit_cache(element, mesh);

  // Because we don't support h-refinement, the direction is the only piece
  // of the neighbor information that we actually need.
  const Direction<VolumeDim> primary_direction = primary_neighbor.first;
  const std::vector<Direction<VolumeDim>> directions_to_exclude =
      [&neighbors_to_exclude]() noexcept {
        std::vector<Direction<VolumeDim>> result(neighbors_to_exclude.size());
        for (size_t i = 0; i < result.size(); ++i) {
          result[i] = neighbors_to_exclude[i].first;
        }
        return result;
      }();

  const DataVector& w = cache.quadrature_weights;
  const DirectionMap<VolumeDim, Matrix>& interp_matrices =
      cache.interpolation_matrices;
  const DirectionMap<VolumeDim, DataVector>& w_dot_interp_matrices =
      cache.quadrature_weights_dot_interpolation_matrices;

  // Use cache if possible, or compute matrix if we are in edge case
  const Matrix& inverse_a =
      LIKELY(directions_to_exclude.size() < 2)
          ? cache.retrieve_inverse_a_matrix(primary_direction,
                                            directions_to_exclude)
          : inverse_a_matrix(mesh, element, w, interp_matrices,
                             w_dot_interp_matrices, primary_direction,
                             directions_to_exclude);

  const DataVector b = b_vector<Tag>(tensor_index, mesh, w, interp_matrices,
                                     w_dot_interp_matrices, neighbor_data,
                                     primary_neighbor, neighbors_to_exclude);

  const size_t number_of_points = b.size();
  const DataVector inverse_a_times_b = apply_matrices(
      std::array<std::reference_wrapper<const Matrix>, 1>{{inverse_a}}, b,
      Index<1>(number_of_points));
  const DataVector inverse_a_times_w = apply_matrices(
      std::array<std::reference_wrapper<const Matrix>, 1>{{inverse_a}}, w,
      Index<1>(number_of_points));

  // Compute Lagrange multiplier:
  // Note: we take w as an argument (instead of as a lambda capture), because
  //       some versions of Clang incorrectly warn about capturing w.
  const double lagrange_multiplier = [&number_of_points, &inverse_a_times_b,
                                      &inverse_a_times_w,
                                      &u](const DataVector& local_w) noexcept {
    double numerator = 0.;
    double denominator = 0.;
    for (size_t s = 0; s < number_of_points; ++s) {
      numerator += local_w[s] * (inverse_a_times_b[s] - u[s]);
      denominator += local_w[s] * inverse_a_times_w[s];
    }
    return -numerator / denominator;
  }(w);

  // Compute solution:
  *constrained_fit_result =
      inverse_a_times_b + lagrange_multiplier * inverse_a_times_w;
}

/*!
 * \ingroup LimitersGroup
 * \brief Compute the HWENO solution for one tensor
 *
 * The HWENO limiter reconstructs a new solution from the linear combination of
 * the local DG solution and a "modified" solution from each neighbor element.
 * This function computes the modified solution for a particular tensor and
 * neighbor, following Section 3 of \cite Zhu2016.
 *
 * The modified solution associated with a particular neighbor (the "primary"
 * neighbor) is obtained by solving a constrained fit over the local element,
 * the primary neighbor, and the other ("secondary") neighbors of the local
 * element. This fit seeks to minimize in a least-squared sense:
 * 1. The distance between the modified solution and the original solution on
 *    the primary neighbor.
 * 2. The distance between the cell average of the modified solution and the
 *    cell average of the original solution on each secondary neighbor. Note
 *    however that one secondary neighbor (or, rarely, several) is excluded from
 *    this minimization: the secondary neighbor(s) where the original solution
 *    has the most different cell average from the local element. This helps to
 *    prevent an outlier (e.g., near a shock) from biasing the fit.
 *
 * The constraint on the minimization is the following: the cell average of the
 * modified solution on the local element must equal the cell average of the
 * local element's original solution.
 *
 * Below we give the mathematical form of the constraints described above and
 * show how these are translated into a numerical algorithm.
 *
 * Consider an element \f$I_0\f$ with neighbors \f$I_1, I_2, ...\f$. For a
 * given tensor component \f$u\f$, the values on each of these elements are
 * \f$u^{(0)}, u^{(1)}, u^{(2)}, ...\f$. Taking for the sake of example the
 * primary neighbor to be \f$I_1\f$, the modified solution \f$\phi\f$ must
 * minimize
 *
 * \f[
 * \chi^2 = \int_{I_1} (\phi - u^{(1)})^2 dV
 *          + \sum_{\ell \in L} \left(
 *                              \int_{I_{\ell}} ( \phi - u^{(\ell)} ) dV
 *                              \right)^2,
 * \f]
 *
 * subject to the constaint
 *
 * \f[
 * C = \int_{I_0} ( \phi - u^{(0)} ) dV = 0.
 * \f]
 *
 * where \f$\ell\f$ ranges over a subset \f$L\f$ of the secondary neighbors.
 * \f$L\f$ excludes the one (or more) secondary neighbor(s) where the mean of
 * \f$u\f$ is the most different from the mean of \f$u^{(0)}\f$. Typically, only
 * one secondary neighbor is excluded, so \f$L\f$ contains two fewer neighbors
 * than the total number of neighbors to the element. Note that in 1D, this
 * implies that \f$L\f$ is the empty set; for each modified solution, one
 * neighbor is the primary neighbor and the other is the excluded neighbor.
 *
 * The integrals are evaluated by quadrature. We denote the quadrature weights
 * by \f$w_s\f$ and the values of some data \f$X\f$ at the quadrature
 * nodes by \f$X_s\f$. We use subscripts \f$r,s\f$ to denote quadrature nodes
 * on the neighbor and local elements, respectively. The minimization becomes
 *
 * \f[
 * \chi^2 = \sum_r w^{(1)}_r ( \phi_r - u^{(1)}_r )^2
 *          + \sum_{\ell \in L} \left(
 *                              \sum_r w^{(\ell)}_r ( \phi_r - u^{(\ell)}_r )
 *                              \right)^2,
 * \f]
 *
 * subject to the constraint
 *
 * \f[
 * C = \sum_s w^{(0)}_s ( \phi_s - u^{(0)}_s ) = 0.
 * \f]
 *
 * Note that \f$\phi\f$ is a function defined on the local element \f$I_0\f$,
 * and so is fully represented by its values \f$\phi_s\f$ at the quadrature
 * points on this element. When evaluating \f$\phi\f$ on element \f$I_{\ell}\f$,
 * we obtain the function values \f$\phi_r\f$ by polynomial extrapolation,
 * \f$\phi_r = \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s\f$, where
 * \f$\mathcal{I}^{(\ell)}_{rs}\f$ is the interpolation/extrapolation matrix
 * that interpolates data defined at grid points \f$x_s\f$ and evaluates it at
 * grid points \f$x_r\f$ of \f$I_{\ell}\f$. Thus,
 *
 * \f[
 * \chi^2 = \sum_r w^{(1)}_r
 *                 \left(
 *                 \sum_s \mathcal{I}^{(1)}_{rs} \phi_s - u^{(1)}_r
 *                 \right)^2
 *          + \sum_{\ell \in L} \left(
 *                              \sum_r w^{(\ell)}_r
 *                                     \left(
 *                                     \sum_s \mathcal{I}^{(\ell)}_{rs} \phi_s
 *                                             - u^{(\ell)}_r
 *                                     \right)
 *                      \right)^2.
 * \f]
 *
 * The solution to this optimization problem is found in the standard way,
 * using a Lagrange multiplier \f$\lambda\f$ to impose the constraint:
 *
 * \f[
 * 0 = \frac{d}{d \phi_s} \left( \chi^2 + \lambda C \right).
 * \f]
 *
 * Working out the differentiation with respect to \f$\phi_s\f$ leads to the
 * linear problem that must be inverted to obtain the solution,
 *
 * \f[
 * 0 = A_{st} \phi_t - b_s - \lambda w^{(0)}_s,
 * \f]
 *
 * where
 *
 * \f{align*}{
 * A_{st} &= \sum_r \left( w^{(1)}_r
 *                  \mathcal{I}^{(1)}_{rs} \mathcal{I}^{(1)}_{rt}
 *                  \right)
 *          + \sum_{\ell \in L} \left(
 *                              \sum_r \left(
 *                                     w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rt}
 *                                     \right)
 *                              \cdot
 *                              \sum_r \left(
 *                                     w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
 *                                     \right)
 *                              \right)
 * \\
 * b_s &= \sum_r \left( w^{(1)}_r u^{(1)}_r \mathcal{I}^{(1)}_{rs} \right)
 *         + \sum_{\ell \in L} \left(
 *                             \sum_r \left( w^{(\ell)}_r u^{(\ell)}_r \right)
 *                             \cdot
 *                             \sum_r \left(
 *                                    w^{(\ell)}_r \mathcal{I}^{(\ell)}_{rs}
 *                                    \right)
 *                             \right).
 * \f}
 *
 * Finally, the solution to the constrained fit is
 *
 * \f{align*}{
 * \lambda &= - \frac{ \sum_s w^{(0)}_s
 *                            \left( (A^{-1})_{st} b_t - u^{(0)}_s \right)
 *              }{ \sum_s w^{(0)}_s (A^{-1})_{st} w^{(0)}_t }
 * \\
 * \phi_s &= (A^{-1})_{st} ( b_t + \lambda w^{(0)}_t ).
 * \f}
 *
 * Note that the matrix \f$A\f$ does not depend on the values of the tensor
 * \f$u\f$, so its inverse \f$A^{-1}\f$ can be precomputed and stored.
 *
 * \warning
 * Note also that the implementation currently does not support h- or
 * p-refinement; this is checked by some assertions. The implementation is
 * untested for grids where elements are curved, and it should not be expected
 * to work in these cases.
 *
 * When calling `hweno_impl`, `modified_neighbor_solution_buffer` should contain
 * one DataVector for each neighboring element (i.e. for each entry in
 * `neighbor_data`).
 */
template <typename Tag, size_t VolumeDim, typename PackagedData>
void hweno_impl(
    const gsl::not_null<std::unordered_map<
        std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, DataVector,
        boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>*>
        modified_neighbor_solution_buffer,
    const gsl::not_null<typename Tag::type*> tensor,
    const double neighbor_linear_weight, const Mesh<VolumeDim>& mesh,
    const Element<VolumeDim>& element,
    const std::unordered_map<
        std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>, PackagedData,
        boost::hash<std::pair<Direction<VolumeDim>, ElementId<VolumeDim>>>>&
        neighbor_data) noexcept {
  ASSERT(
      modified_neighbor_solution_buffer->size() == neighbor_data.size(),
      "modified_neighbor_solution_buffer->size() = "
          << modified_neighbor_solution_buffer->size()
          << "\nneighbor_data.size() = " << neighbor_data.size()
          << "\nmodified_neighbor_solution_buffer was incorrectly initialized "
             "before calling hweno_impl.");

  alg::for_each(
      neighbor_data, [&element, &mesh](const auto& neighbor_and_data) noexcept {
        ASSERT(Weno_detail::check_element_has_one_similar_neighbor_in_direction(
                   element, neighbor_and_data.first.first),
               "Found some amount of h-refinement; this is not supported");
        ASSERT(neighbor_and_data.second.mesh == mesh,
               "Found some amount of p-refinement; this is not supported");
      });

  for (size_t tensor_index = 0; tensor_index < tensor->size(); ++tensor_index) {
    const auto& tensor_component = (*tensor)[tensor_index];
    for (const auto& neighbor_and_data : neighbor_data) {
      const auto& primary_neighbor = neighbor_and_data.first;
      const auto neighbors_to_exclude =
          secondary_neighbors_to_exclude_from_fit<Tag>(
              mean_value(tensor_component, mesh), tensor_index, neighbor_data,
              primary_neighbor);

      DataVector& buffer =
          modified_neighbor_solution_buffer->at(primary_neighbor);
      solve_constrained_fit<Tag>(make_not_null(&buffer), tensor_component,
                                 tensor_index, mesh, element, neighbor_data,
                                 primary_neighbor, neighbors_to_exclude);
    }

    // Sum local and modified neighbor polynomials for the WENO reconstruction
    Weno_detail::reconstruct_from_weighted_sum(
        make_not_null(&((*tensor)[tensor_index])), neighbor_linear_weight,
        Weno_detail::DerivativeWeight::PowTwoEllOverEllFactorial, mesh,
        *modified_neighbor_solution_buffer);
  }
}

}  // namespace Limiters::Weno_detail

max_stars_repo_name: Ambrou/spectre

max_issues_repo_name: Ambrou/spectre

max_forks_repo_name: Ambrou/spectre


text: #include <boost/move/algo/adaptive_merge.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino