//
//  MinglingMarkDifferentiation.cpp
//  
//
//  Created by Arne Pommerening on 1/29/13.
//  Copyright 2013 Philodendron International. All rights reserved.
//

#include <Rcpp.h>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace Rcpp;
using namespace std;


/**
 * Calculates expected differentiation (Pommerening et al., 2011).
 */
double calcExpectedDiff(NumericVector zSort) {
    int n = zSort.size();
	vector<double> help_r(n);
	help_r[0] = 0;
	for (int i = 1; i < n; i++)
		help_r[i] = help_r[i - 1] + zSort[i - 1];
	double TdExp = 0;
	for (int i = 0; i < n; i++) {
		double dummy = help_r[i] / zSort[i];
		TdExp += 2 * dummy;
	}
	TdExp = 1 - 1 / ((n - 1.0) * n) * TdExp;
	// cout << "TdExp: ";
	// cout << TdExp;
	// cout << endl;
	return TdExp;
}

/**
 * Estimates the mark mingling and mark differentiation function (Pommerening et al., 2011).
 */
// [[Rcpp::export]]
DataFrame calcMingDiffFunctions(DataFrame TreeList, double xmax, double ymax, int m, double dd, double h, NumericVector speciesFreq, NumericVector zSort, int kernel) {
	vector<double> x = as< vector<double> > (TreeList["x"]);
    vector<double> y = as< vector<double> > (TreeList["y"]);
    vector<double> dbh = as< vector<double> > (TreeList["dbh"]);
    vector<double> species = as< vector<double> > (TreeList["Species"]);
    const double pi = 4.0 * atan(1.0);
    int n = dbh.size();
    double la = n / (xmax * ymax);
    double ql = la * la;
	NumericVector g(m),  r(m), xm(m), ming(m), xl(m), xn(m), xk(m);
    for(int i = 0; i < m; i++) {
		r[i] = 0;
        g[i] = 0;
        xk[i] = 0;
        xm[i] = 0;
        ming[i] = 0;
        xl[i] = 0;
        xn[i] = 0;
    }
    double rz = m * dd + h;
    for(int i = 1; i < n; i++) {
        double xx = x[i];
        double yy = y[i];
        double zz = dbh[i];
        for(int j = 0; j < i; j++) {
            double cx = fabs(xx - x[j]);
            if(cx < rz) {
                double cy = fabs(yy - y[j]);
                if(cy < rz) {
                    double tz = sqrt(cx * cx + cy * cy);
                    if(tz < rz) {
                        double cc = 1 / ((xmax - cx) * (ymax - cy));
                        for(int l = 0; l < m; l++) {
                            double xr = (l + 1) * dd;
                            if(tz <= xr) {
								xk[l] += 2 * cc; 
								xl[l] += 2 * cc * (1 - min(zz, dbh[j]) / max(zz, dbh[j])); 
                                if(species[i] != species[j])
								    xn[l] += 2 * cc;
							}
							double t = fabs(xr - tz);
                            if(t <= h) {
                                double rh = 0;
                                if(kernel == 0)
                                    rh = (1 - pow(t / h, 2)) * 0.75 / h;
                                else
                                    rh = 0.5 / h;
                                double cr = pi * xr;
                                rh *= cc / cr;
                                g[l] += rh;
                                xm[l] += rh * (1 - min(zz, dbh[j]) / max(zz, dbh[j]));
                                if(species[i] != species[j])
								  ming[l] += rh;                                
                            }
                        }
                    }
                }
            }    
        }
    }
    double em = 0;
    for(int i = 0; i < speciesFreq.size(); i++)
        em += (speciesFreq[i] * (n - speciesFreq[i])) / (n * (n - 1));
    double ed = calcExpectedDiff(zSort);
    for(int l = 0; l < m; l++) {
		xk[l] = sqrt(xk[l] / (pi * ql));
        r[l] = (l + 1) * dd; 
        xl[l] = sqrt(xl[l] / (pi * ql * ed));
        xn[l] = sqrt(xn[l] / (pi * ql * em));
		if(g[l] < 1E-20) {
            xm[l] = -1;
            ming[l] = -1;
        }
        else {
            xm[l] /= (g[l] * ed);
            ming[l] /= (g[l] * em);
        }
    }
    return DataFrame::create(Named("r") = r, Named("L") = xk, Named ("tau") = xm, Named ("Lt") = xl, Named ("nu") = ming, Named ("Ln") = xn);
}

/**
 * Shifts all trees of species "mark" by the random measures zx and cy. Includes periodic boundary conditions.
 */
// [[Rcpp::export]]
DataFrame shiftSpeciesGroup(DataFrame TreeList, NumericVector marks, double zx, double zy, double xmax, double ymax) {
	vector<double> x = as< vector<double> > (TreeList["x"]);
    vector<double> y = as< vector<double> > (TreeList["y"]);
    vector<double> species = as< vector<double> > (TreeList["Species"]);
    vector<int> xmarks = as< vector<int> > (marks);
    int n = x.size();
    for(int i = 0; i < n; i++) {
		if(std::find(xmarks.begin(), xmarks.end(), species[i]) != xmarks.end()) {
			x[i] += zx;
			if(x[i] > xmax)
			   x[i] -= xmax;
			y[i] += zy;
			if(y[i] > ymax)
				y[i] -= ymax;
		}
	}
    return DataFrame::create(Named("x") = x, Named ("y") = y);
}

/**
 * Calculates the minimum distance between points in an observation window.
 */
// [[Rcpp::export]]
double calcMinDistance(NumericVector x, NumericVector y) {
	double dm = 1E50;
	double xx = 0;
    int n = x.size();
    for(int i = 1; i < n; i++) {
        for(int j = 0; j < i; j++) {
			xx = pow(x[i] - x[j], 2) + pow(y[i] - y[j], 2);
			dm = min(xx, dm);
		}
	}
	return sqrt(dm);
}