#include <Rcpp.h>
#include <iostream>
using namespace Rcpp;


//int foo(int a, int b){
//  return a+b;}

// [[Rcpp::export]]
NumericVector inclusionprobabilitiesCpp(NumericVector a, double n){
  int nnull    = 0;  //Count probabilities = 0
  int nneg     = 0;  //Count probabilities < 0
  double sumA  = 0;  //Sum of elements of a

  for(int i=0;i<a.size();i++){
    if(a[i]==0){
      nnull++;}
    if(a[i]<0){
      nneg++;
      a[i] = 0;}
    sumA += a[i];
  }
  NumericVector pik1;

  if(nnull == a.size()){
    pik1    = a;}
  else{ // **
    pik1    = (n/sumA)*a;
    NumericVector tempPik = pik1;
    int pikSize           = 0;
    int pikCursor         = 0;

    for(int i=0;i<pik1.size();i++){
      if(pik1[i] <= 0){
	pik1[i] = 0;
      }
      else{
	pikSize++;}
    }
    //create, fill pik
    NumericVector pik(pikSize);

    for(int i=0;i<pik1.size();i++){
      if(pik1[i] != 0){
	pik[pikCursor] = pik1[i];
	pikCursor++;}
    }

    LogicalVector list1 = pik1 > 0;
    LogicalVector list2 = pik >= 1;
    int l               = 0;

    //find l, the number of elements of pik >= 1
    for(int i=0;i<list2.size();i++){
      if(list2[i]==true){
	l++;}
    }//end find l
    if(l>0){ //**
      int l1 = 0;
      NumericVector tempX(pik.size());
      while(l != l1){
	int xSize = pik.size();
	tempX = pik;
	//Where values are true, pik >=1 and we aren't modifying it yet
	// Instead, here we set other values x to (n-l)*x/(sum x)
	for(int i=0;i<pik.size();i++){
	  if(list2[i] == true){
	    tempX[i] = -1;
	    xSize--;
	  }
	  //std::cout<< tempX[i]<<" ";
	}
	NumericVector x(xSize);
	int xCursor = 0;
	double sumX = 0;
        for(int i=0;i<pik.size();i++){
	  if(list2[i] == false){

	    x[xCursor] = tempX[i];
	    sumX += tempX[i];
	    xCursor++;
	  }
	}// end define x
        x = x/sumX;
	xCursor = 0;
	for(int i=0;i<pik.size();i++){
	  if(list2[i]==false){
	    pik[i] = (n-l)*x[xCursor];
	    xCursor++;
	  }
	  else{
	    pik[i] = 1;
	  }
	}
	l1 = l;
	list2  = pik >=1;
	l = 0; //reset l
	//find l, the number of elements of pik >= 1
	for(int i=0;i<list2.size();i++){
	  if(pik[i]>=1){
	    l++;}
	}//end find l
      } //end while

      //pik1[list1]=pik
      pikCursor = 0;
      for(int i=0;i<pik1.size();i++){
	if(list1[i]==true){
	  pik1[i]=pik[pikCursor];
	  pikCursor++;
	}
      }
    } //end if **
  } //end else **
  return pik1; //a NumericVector
}

// [[Rcpp::export]]
NumericVector UPtilleCpp(NumericVector pik, double eps=0.000001){
  std::cout<<"oh hai";
  if(any(is_na(pik))){
    std::cout<<"there are missing values in the pik vector \n";
    return pik;}

    double n_double = 0;
    for(int i=0;i<pik.size();i++){
      n_double += pik[i];
    }
    int n = int(floor(n_double+.5)); //round n_double to nearest int
    //std::cout<<n_double <<" "<<n;

    LogicalVector list1 = (pik > eps);
    LogicalVector list2 = (pik < (1 - eps));
    LogicalVector list  = list1*list2;
    int pikbLength      = 0;

    for(int i =0;i<list.size();i++){
    //  std::cout<<list[i]<<" ";
      pikbLength+=list[i];
    }

    //Define pikb
    NumericVector pikb(pikbLength);
    int pikbCursor = 0;

    for(int i=0;i<list.size();i++){
      if(list[i]==true){
	pikb[pikbCursor] = pik[i];
	pikbCursor++;
      }
    } //end define pikb

    int N = pikbLength;
    NumericVector s = pik;

    if(N<1){
      std::cout<<"the pik vector has all elements outside of the range [eps,1-eps]";
      return pik;
    }
    else{ // **

      n = 0;
      for(int i=0;i<pikb.size();i++){
	n += pikb[i];
      }
      NumericVector sb(N,1.0);
      NumericVector b(N,1.0);

      for(int i=0;i<N-n;i++){ //**
	NumericVector a = inclusionprobabilitiesCpp(pikb,N-i);
//	int xx = foo(1,2);
      }//end for **
    }//end else **

  return pik;
}