############################################################
# make_unique
############################################################
function make_unique!(s)
	count=1
	sp=fsortperm(s)
	current=s[sp[1]]
	for index in 2:length(s)
		if s[sp[index]]==current
			s[sp[index]]=s[sp[index]]*"."*string(count)
			count+=1
		else
			current=s[sp[index]]
			count=1
		end
	end
end

function make_unique(s)
	r=similar(s)
	count=1
	sp=fsortperm(s)
	r[sp[1]]=s[sp[1]]
	for index in 2:length(s)
		if s[sp[index]]==s[sp[index-1]]
			r[sp[index]]=s[sp[index]]*"."*string(count)
			count+=1
		else
			r[sp[index]]=s[sp[index]]
			count=1
		end
	end
	r
end

############################################################
# read_csv
############################################################
function read_csv(path::String)
    df = GZip.open(path, "r") do io
        CSV.read(io, DataFrame, header = false)
    end
    return(df)
end

############################################################
# read_mtx
############################################################
function read_mtx(filename::String)
    _parseint(x) = parse(Int, x)    
    GZip.open(filename, "r") do mmfile
        # Read first line
        firstline = chomp(readline(mmfile))
        tokens = split(firstline)
        if length(tokens) != 5
            throw(string("Not enough words on first line: ", firstline))
        end
        if tokens[1] != "%%MatrixMarket"
            throw(string("Not a valid MatrixMarket header:", firstline))
        end
        (head1, rep, field, symm) = map(lowercase, tokens[2:5])
        if head1 != "matrix"
            throw("Unknown MatrixMarket data type: $head1 (only \"matrix\" is supported)")
        end

        eltypes = field == "real" ? Float32 :
                 field == "integer" ? Int32 :
                 throw(ParseError("Unsupported field $field (only real and complex are supported)"))

        symlabel = symm == "general" ? identity :
                   throw("Unknown matrix symmetry: $symm (only general, symmetric, skew-symmetric and hermitian are supported)")

        # Skip all comments and empty lines
        ll   = readline(mmfile)
        while length(chomp(ll))==0 || (length(ll) > 0 && ll[1] == '%')
            ll = readline(mmfile)
        end
        # Read matrix dimensions (and number of entries) from first non-comment line
        dd = map(_parseint, split(ll))
        if length(dd) < (rep == "coordinate" ? 3 : 2)
            throw(ParseError(string("Could not read in matrix dimensions from line: ", ll)))
        end
        rows = dd[1]
        cols = dd[2]
        dataf = CSV.read(mmfile, DataFrame, header = false, type = eltypes, delim = " ");
        return symlabel(sparse(dataf.Column1, dataf.Column2, dataf.Column3, rows, cols))
    end
end

############################################################
# sum_same_genes
############################################################
function sum_same_genes(raw_matrix::Array{Int32,2}, features::Vector{String}, barcodes::Vector{String})::NamedArray
    if (features |> unique |> length) == (features |> length)
        FinalMat = NamedArray(raw_matrix)
        setnames!(FinalMat, features, 1)
        setnames!(FinalMat, barcodes, 2)
    else
            # Count occurence of each genes
            countmap_features =  countmap(features)
            #Separate non unique genes with unique ones
            non_unique_features = (collect(keys(countmap_features))[values(countmap_features) .> 1])
            unique_features = collect(keys(countmap_features))[values(countmap_features) .== 1]
            #split matrix with those genes
            sub_X_unique_features = raw_matrix[in.(df_features.Column2, (unique_features,)),:]
            features1 = df_features.Column2[in.(df_features.Column2, (unique_features,))]
            sub_X_non_unique_features::SparseMatrixCSC{Int32,Int32} = reduce(vcat, map(x -> sum(raw_matrix[df_features.Column2 .== x, :], dims = 1), non_unique_features))
            FinalMat = vcat(sub_X_unique_features, sub_X_non_unique_features)
            ordered_features = append!(features1, non_unique_features)
            FinalMat = NamedArray(FinalMat)
            setnames!(FinalMat, ordered_features, 1)
            setnames!(FinalMat, barcodes[1], 2)
    end
    return(FinalMat)
end       

############################################################
# read_10X
############################################################
function read_10X(path::String; features_col::Int = 2, barcodes_col::Int = 1, unique_features::Bool = true)    
    # Check directory files
    dir_file = readdir(path)

    # Get files 
    barc_path = dir_file[occursin.(r"^barcode", dir_file)]
    feat_path = dir_file[occursin.(r"^gene|^feature", dir_file)]
    rmat_path = dir_file[occursin.(r"^matrix", dir_file)] 
    
    # Check if file exists and get full path
    isempty(barc_path) ? throw("$dir_file doesn't contain barcode file") : barc_path = join([path, barc_path[1]])
    isempty(feat_path) ? throw("$dir_file doesn't contain feature file") : feat_path = join([path, feat_path[1]])
    isempty(rmat_path) ? throw("$dir_file doesn't contain matrix file") : rmat_path = join([path, rmat_path[1]])

    # Read different components of teh 10X folder
    println("\nReading Matrix... \n")
    rmat = read_mtx(rmat_path);
    println("\nReading Features... \n")
    feat_df = read_csv(feat_path);
    println("\nReading Cell Barcodes... \n")
    barc_df = read_csv(barc_path);
    
    # Construct Named Arrays
    println("\nConstructing Named Arrays... \n")
    feat = feat_df[!, features_col]
    barc = barc_df[!, barcodes_col]
    feat_df_end = ncol(feat_df)
    if  ncol(feat_df) > 2
        list_assay = Dict{String,NamedArray{Int32,2}}()
        type_vec =  feat_df[!, feat_df_end]
        type =  type_vec |> unique
        for i in type
            feat_ind = feat_df[!,feat_df_end] .== i
            feat_split = feat_df[feat_ind, features_col];
            rmat_split = rmat[feat_ind,:];
                if unique_features
                    list_assay[i] = NamedArray(rmat_split, (make_unique(feat_split), barc))
                else
                    list_assay[i] = sum_same_genes(rmat_split, make_unique(feat_split), barc)
                end      
        end
    else
        if unique_features
            list_assay = NamedArray(rmat, (make_unique(feat), barc))
        else    
            list_assay = sum_same_genes(rmat, make_unique(feat), barc)
        end    
    end
    println("\nDONE! \n")
    return(list_assay)
end

############################################################
# init_pissaro
############################################################
function init_pissaro(matrix::AbstractArray; min_cells::Int = 5, min_features::Int = 0, assay_name::String = "RNA")::Pissaro
    # Assay
    nzm =  matrix.array .!= 0
    n_cell = Vector{Int32}((sum(nzm .!= 0, dims = 2))[:,1])
    n_feat = Vector{Int32}((sum(nzm .!= 0, dims = 1))[1,:])
    min_cell_vec = n_cell .>= min_cells
    min_feat_vec = n_feat .>= min_features
    matrix = matrix[min_cell_vec, min_feat_vec]
    features_data = DataFrame(features = names(matrix,1), nCells = n_cell[min_cell_vec])
    cell_data = DataFrame(cells = names(matrix, 2), nFeatures = n_feat[min_feat_vec], nCount = Vector{Int32}(sum(matrix.array, dims = 1)[1,:]))
    sp_matrix::SparseMatrixCSC{Float32,Int32} = spzeros(size(matrix)[1], size(matrix)[2]) 
    empty_mat = NamedArray(sp_matrix, (names(matrix,1), names(matrix,2)))
    dict_assay = Dict(assay_name => Assays(matrix, deepcopy(empty_mat), NamedArray{Float32}(1,1), features_data, Array{String,1}(), assay_name))
    dim_red = Dict{String,Reductions}()
    pissaro_object = Pissaro(dict_assay, dim_red, cell_data, assay_name, "")
    return(pissaro_object)
end

############################################################
# FastExpMean
############################################################
function FastExpMean(matrix::AbstractArray)
    a = mean(expm1.(matrix), dims = 2)
    res = log1p.(a)[:,1]
    return res
end   
############################################################
# FastLogVMR
############################################################
function FastLogVMR(matrix::AbstractArray)
    rmat = expm1.(matrix)
    a = mean(rmat, dims = 2)
    b = var(rmat , dims = 2)
    res = log1p.(a ./ b)[:,1]
    return res
end