//
//  PartialPairMarkConnection.cpp
//  
//
//  Created by Arne Pommerening on 10/25/13.
//  Copyright 2013 Philodendron International. All rights reserved.
//

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


/**
 * Estimates the partial pair and mark connection functions.
 */
// [[Rcpp::export]]
DataFrame calcPartialPairMarkConnection(DataFrame TreeList, double xmax, double ymax, NumericVector speciesFreq, NumericVector speciesCodes, int m, double dd, double h, int kernel) {
	vector<double> x = as< vector<double> > (TreeList["x"]);
    vector<double> y = as< vector<double> > (TreeList["y"]);
    vector<int> species = as< vector<int> > (TreeList["species"]);
    // cout << m << endl;
    const double pi = 4.0 * atan(1.0);
    int n = species.size();
    double la = n / (xmax * ymax);
    double ql = la * la;
	NumericVector r(m), g(m), x11(m), x12(m), x22(m), x33(m), x13(m), x23(m), g11(m), g12(m), g22(m), g33(m), g13(m), g23(m);
    NumericVector ll(speciesFreq.size());
    for(int i = 0; i < speciesFreq.size(); i++) {
        ll[i] = speciesFreq[i] / (xmax * ymax);
    }        
    for(int i = 0; i < m; i++) {
		r[i] = 0;
        g[i] = 0;
        x11[i] = 0;
        x12[i] = 0;
        x22[i] = 0;
        x33[i] = 0;
        x13[i] = 0;
        x23[i] = 0;
        g11[i] = 0;
        g12[i] = 0;
        g22[i] = 0;
        g33[i] = 0;
        g13[i] = 0;
        g23[i] = 0;
    }
    double rz = m * dd + h;
    for(int i = 1; i < n; i++) {
        double xx = x[i];
        double yy = y[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;
							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;
                                if(speciesFreq.size() > 0)
                                    if((species[i] == speciesCodes[0]) && (species[j] == speciesCodes[0])) 
                                        x11[l] += rh;
                                if(speciesFreq.size() > 1){
                                    if(((species[i] == speciesCodes[0]) && (species[j] == speciesCodes[1])) || ((species[i] == speciesCodes[1]) && (species[j] == speciesCodes[0])) ) 
                                        x12[l] += rh;
                                    if((species[i] == speciesCodes[1]) && (species[j] == speciesCodes[1])) 
                                        x22[l] += rh;
                                }
                                if(speciesFreq.size() > 2){
                                    if((species[i] == speciesCodes[2]) && (species[j] == speciesCodes[2])) 
                                        x33[l] += rh;
                                    if(((species[i] == speciesCodes[0]) && (species[j] == speciesCodes[2])) || ((species[i] == speciesCodes[2]) && (species[j] == speciesCodes[0])) )                                     
                                        x13[l] += rh;
                                    if(((species[i] == speciesCodes[1]) && (species[j] == speciesCodes[2])) || ((species[i] == speciesCodes[2]) && (species[j] == speciesCodes[1])) )                                     
                                        x23[l] += rh;
                                }    
                            }
                        }
                    }
                }
            }    
        }
    }
    for(int l = 0; l < m; l++) {
		r[l] = (l + 1) * dd; 
        if(g[l] < 1E-20) {
            if(speciesFreq.size() > 0) {
                x11[l] = -1;
                g11[l] = -1;
            }
            if(speciesFreq.size() > 1) {
                x12[l] = -1;
                x22[l] = -1;
                g12[l] = -1;
                g22[l] = -1;
            }
            if(speciesFreq.size() > 2) {
                x13[l] = -1;
                x23[l] = -1;
                x33[l] = -1;
                g13[l] = -1;
                g23[l] = -1;
                g33[l] = -1;
            }
        }
        else {
            if(speciesFreq.size() > 0) {
                if(ll[0] > 0)
                    g11[l] = x11[l] / pow(ll[0], 2);
                x11[l] /= g[l];
            }
            if(speciesFreq.size() > 1) {
                if((ll[0] > 0) && (ll[1] > 0))
                    g12[l] = x12[l] / (2 * ll[0] * ll[1]);
                if(ll[1] > 0)
                    g22[l] = x22[l] / pow(ll[1], 2);  
                x12[l] /= g[l];
                x22[l] /= g[l];
            }
            if(speciesFreq.size() > 2) {
                if((ll[0] > 0) && (ll[2] > 0))
                    g13[l] = x13[l] / (2 * ll[0] * ll[2]);
                if((ll[1] > 0) && (ll[2] > 0))
                    g23[l] = x23[l] / (2 * ll[1] * ll[2]);
                if(ll[2] > 0)
                    g33[l] = x33[l] / pow(ll[2], 2);  
                x13[l] /= g[l];
                x23[l] /= g[l];
                x33[l] /= g[l];
            }
        }
    }
    DataFrame xresult;
    if(speciesFreq.size() == 1)
        xresult = DataFrame::create(Named("r") = r, Named("g11") = g11, Named("p11") = x11);
    if(speciesFreq.size() == 2)
        xresult = DataFrame::create(Named("r") = r, Named("g11") = g11, Named("g12") = g12, Named("g22") = g22, Named("p11") = x11, Named("p12") = x12, Named("p22") = x22);
    if(speciesFreq.size() == 3)
        xresult = DataFrame::create(Named("r") = r, Named("g11") = g11, Named("g12") = g12, Named("g13") = g13, Named("g22") = g22, Named("g23") = g23, Named("g33") = g33,Named("p11") = x11, Named("p12") = x12, Named("p13") = x13, Named("p22") = x22, Named("p23") = x23, Named("p33") = x33); 
    return xresult;
}

/**
 * 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);
}