/*
 *  calcNNSS.cpp
 *  
 *
 *  Created by Arne Pommerening on 20/01/2013.
 *  Copyright 2013 Philodendron International. All rights reserved.
 *
 */

#include <Rcpp.h>
#include <cmath>
#include <vector>


using namespace Rcpp;
using namespace std;


/**
 * Calculates the tree specific mingling variable.
 */
// [[Rcpp::export]]
NumericVector calcMingling(DataFrame xData, int mi) {
	vector<double> species = as< vector<double> > (xData["Species"]);
	NumericMatrix neighbour = as< NumericMatrix > (xData["neighbour"]);
	int n = species.size();
	NumericVector mingling(n);
	for(int i = 0; i < n; i++) {
		double sum = 0;
		for(int j = 0; j < mi; j++) {
          if(species[i] != species[neighbour(i, j) - 1])
			  sum++;
		}
		mingling[i] = sum / mi;
	}
    return mingling;
}

/**
 * Differentiation algorithm.
 */
double calcDiffComparison2(double value1, double value2) {
	double td = value1 / value2;
	if(td > 1)
		td = 1 / td;
	return 1 - td;
}

/**
 * Calculates the tree specific diameter differentiation variable (1st neighbour).
 */
// [[Rcpp::export]]
NumericVector calcDiff1(DataFrame xData) {
	vector<double> dbh = as< vector<double> > (xData["dbh"]);
	NumericMatrix neighbour = as< NumericMatrix > (xData["neighbour"]);
	int n = dbh.size();
	NumericVector td1(n);
	for (int i = 0; i < n; i++)
		td1[i] = calcDiffComparison2(dbh[i], dbh[neighbour(i, 0) - 1]);
	return td1;
}	

/**
 * Determines the differentiation class.
 */
// [[Rcpp::export]]
int tdToClass4(double td) {
	int xclass = 0;
	td += 1E-8;
	if (td < 0.3)
		xclass = 0;
	else if ((td >= 0.3) && (td < 0.5))
		xclass = 1;
	else if ((td >= 0.5) && (td < 0.7))
		xclass = 2;
	else
		xclass = 3;
	return xclass;
}

/**
 * Calculates the tree specific dominance variable.
 */
// [[Rcpp::export]]
NumericVector calcDominance(DataFrame xData, int mi) {
	vector<double> dbh = as< vector<double> > (xData["dbh"]);
	NumericMatrix neighbour = as< NumericMatrix > (xData["neighbour"]);
	int n = dbh.size();
	NumericVector dom(n);
	for(int i = 0; i < n; i++) {
		double sum = 0;
		for(int j = 0; j < mi; j++) {
			if(dbh[i] > dbh[neighbour(i, j) - 1])
				sum++;
		}
		dom[i] = sum / mi;
	}
    return dom;
}

/**
 * Calculates the tree specific diameter variogram index (1st neighbour).
 */
// [[Rcpp::export]]
NumericVector calcVarIndex(DataFrame xData, double xvar) {
	vector<double> dbh = as< vector<double> > (xData["dbh"]);
	NumericMatrix neighbour = as< NumericMatrix > (xData["neighbour"]);
	int n = dbh.size();
	NumericVector varInd(n);
	for (int i = 0; i < n; i++) {
		varInd[i] = 0.5 * pow(dbh[i] - dbh[neighbour(i, 0) - 1], 2);
	    varInd[i] /= xvar;
	}
	return varInd;
}

double translateX(double xmax, double x1, double y1, double x2, double y2) {
	double xx2, dx, dy, r, w;
	dx = x2 - x1;
	dy = y2 - y1;
	r = dx * dx + dy * dy;
	xx2 = x2;
	if(dx > xmax * 0.5) {
		w = pow(dx - xmax, 2);
		if(w + pow(dy, 2) < r)
			xx2 = x2 - xmax;
	}
	if(dx < -xmax * 0.5) {
		w = pow(dx + xmax, 2);
		if(w + pow(dy, 2) < r)
			xx2 = x2 + xmax;
	}
	return xx2;
}

/**
 * Returns translated y2. Modified last on 28.02.2007.
 */
double translateY(double ymax, double x1, double y1, double x2, double y2) {
	double yy2, dx, dy, r, w;
	dx = x2 - x1;
	dy = y2 - y1;
	r = dx * dx + dy * dy;
	yy2 = y2;
	if(dy > ymax * 0.5) {
		w = pow(dy - ymax, 2);
		if(pow(dx, 2) + w < r)
			yy2 = y2 - ymax;
	}
	if(dy < -ymax * 0.5) {
		w = pow(dy + ymax, 2);
		if(pow(dx, 2) + w < r)
			yy2 = y2 + ymax;
	}
	return yy2;
}

/**
 * Calculates the tree specific mean directional index.
 */
// [[Rcpp::export]]
NumericVector calcMDI(DataFrame xData, int mi, int edge, double xmax, double ymax) {
	vector<double> x = as< vector<double> > (xData["x"]);
	vector<double> y = as< vector<double> > (xData["y"]);
	NumericMatrix neighbour = as< NumericMatrix > (xData["neighbour"]);
	int n = x.size();
	NumericVector mdi(n);
	for (int i = 0; i < n; i++) {
	  double ex = 0;
	  double ey = 0;
	  for (int j = 0; j < mi; j++) {
			int neighbourindex = 0;
			double xx2 = x[neighbour(i, j) - 1];
			double yy2 = y[neighbour(i, j) - 1];
			if(edge == 1) {
				xx2 = translateX(xmax, x[i], y[i], x[j], y[j]);
				yy2 = translateY(ymax, x[i], y[i], x[j], y[j]);
			}
			double exx = x[i] - xx2;
			double eyy = y[i] - yy2;
			double dist = sqrt(exx * exx + eyy * eyy);
			ex += exx / dist;
			ey += eyy / dist;
	  }
	  mdi[i] = sqrt(ex * ex + ey * ey);
	}
	return mdi;
}