using DataStructures;
using SparseArrays;
using Serialization;
using DelimitedFiles;

include("Kernel_integral_Names.jl");
include("Pairwise_distance.jl");
#include("knn.jl");
#include("kd_tree.jl");
#include("Data_tools.jl");
#include("Data_structures/int_BinaryTree.jl");

"""
B = Data[offset+1 : offset+N, : ] is an ``N×dim`` array. Here the metric incoporates Q delays, so the effective space is ``R^{Q'dim}``, where ``Q' = |Q|+1`` . This creates a NXN matrix K such that K[i,j] = exp( d_Q^2( B_i, B_j) ).
"""
function Gaussian_kernel( Data::Array{Float64,2}, offset::Int64, N::Int64, Q::Int64, epsilon::Float64, numproc::Int64; divide_by_Q::Bool=true ) :: Tuple
    #------------------System contants------------------------------
    zero_thresh = 1.0E-10;
    delete_temp_files = true;
    normalize_dimension_ranges = true;
    #---------------------------------------------------------------
    pairwise_param = Get_Pairwise_dist_param( Data, Q, epsilon; normalize_dimension_ranges = normalize_dimension_ranges, divide_by_Q = divide_by_Q, offset = offset, N=N);
    K = Gaussian_kernel( Data, offset, N, pairwise_param, numproc);
    return (pairwise_param, K);
end

"""
B = Data[offset+1 : offset+N, : ] is an ``N×dim`` array. Here the metric incoporates Q delays, so the effective space is ``R^{Q'dim}``, where ``Q' = |Q|+1`` . This creates a NXN matrix K such that K[i,j] = exp( d_Q^2( B_i, B_j) ).
"""
function Gaussian_kernel( Data::Array{Float64,2}, offset::Int64, N::Int64, pairwise_param :: Pairwise_dist_param, numproc::Int64) ::SparseMatrixCSC{Float64,Int64}
    #------------------System contants------------------------------
    delete_temp_files = true;
    #---------------------------------------------------------------
    total_count = Pairwise_dist_UpperDiag(Data, offset, N, pairwise_param, numproc = numproc);
    GC.gc();
    counts = readdlm(Gaussian_kernel_counts_name(), Int64); counts = vec(counts); print("\nAbout to combine now");
    K = Combine_PairDist_into_GaussKer(counts, offset, pairwise_param.epsilon, numproc ); print("\nCombined");
    if(delete_temp_files)
      rm( Gaussian_kernel_counts_name(), force=true );
      Threads.@threads for i=1:numproc
        rm( Pairwise_dist_tempname(i), force=true );
      end
    end
    GC.gc();
    return K;
end
#=============================================================#
function Combine_PairDist_into_GaussKer(counts::Array{Int64,1}, offset::Int64, epsilon::Float64, numproc::Int64) :: SparseMatrixCSC{Float64,Int64}
  total_count = sum(counts);
  x_index= zeros(Int64, total_count);
  y_index= zeros(Int64, total_count);
  value=zeros(Float64, total_count);
  #counter=0;
  #.....................Update limits...........................
  index_limits = Array{Int64,2}(undef,numproc, 2);
  index_limits[1,1] = 1; index_limits[1,2] = counts[1];
  for i=2:numproc
    index_limits[i,1] = index_limits[i-1,2]+1;
    index_limits[i,2] = index_limits[i-1,2]+counts[i];
  end
  #.............................................................
  #end_index=0;
  Threads.@threads for i=1:numproc #
    #start_index=end_index+1;    end_index += counts[i];
    Combine_PairDist_into_GaussKer_thread( x_index, y_index, value, index_limits, i, offset, epsilon);
  end
  print("\nDistances collected in arrays");#sleep(2);
  K = sparse(x_index,y_index,value); print("\nArray collected");
  #......................................................
  #...................Clear memory.......................
  index_limits=1.0; x_index=1.0; y_index=1.0; value=1.0;
  GC.gc();
  #......................................................
  return K;
end
#=============================================================#
function Combine_PairDist_into_GaussKer_thread( x_index::Array{Int64,1}, y_index::Array{Int64,1}, value::Array{Float64,1}, index_limits::Array{Int64,2}, threadid::Int64, offset::Int64, epsilon :: Real  )
  filename = Pairwise_dist_tempname(threadid); #print("\nReading : ", filename);display(index_limits[threadid,:]);
  open(filename) do filereader
    for i = index_limits[threadid,1] : index_limits[threadid,2]
      index = parse( Int64, readline(filereader) ); #if(i== 190151243) print("\n",index);end
      x_index[i] = index - offset;
      index = parse( Int64, readline(filereader) ); #if(i== 190151243) print("\n",index);end
      y_index[i] = index - offset;
      d=parse( Float64, readline(filereader) ); #if(i== 190151243) print("\n",d);end
      value[i] = exp(-d/epsilon);
    end
  end
  #print("...finished");
end
