# using SuperResolution

mutable struct EvalSR
    _selectIDX::Array{Tuple{Int64,Int64},1}
    _modus::Int64
    wModel::Array{SuperResolution.ModelSR,1}
    RefModel::Array{SuperResolution.ModelSR,1}
    distKL::Array{SuperResolution.DistanceSR,1}
    wkldSR::Array{SuperResolution.clusWKLD,1}
    gapSR::Array{SuperResolution.clusGAP,1}

    max_mean_idx::Array{Int64,1}
    max_median_idx::Array{Int64,1}
    max_max_idx::Array{Int64,1}
    max_minDelta_idx::Array{Int64,1}

    diff_WK_max_idx::Array{Int64,1}
    diff_GAP_max_idx::Array{Int64,1}

    N_observations::Array{Int64,1}
    K_ground_truth::Array{Int64,1}

    K_KL_estimated_by_mean::Array{Int64,1}
    K_KL_estimated_by_median::Array{Int64,1}
    K_KL_estimated_by_max::Array{Int64,1}
    K_KL_estimated_by_minDelta::Array{Int64,1}
    K_KL_estimated_by_95quantile::Array{Int64,1}

    K_GAP_estimated::Array{Int64,1}

    corrDetecRate_KL_by_mean::Float64
    corrDetecRate_KL_by_median::Float64
    corrDetecRate_KL_by_max::Float64
    corrDetecRate_KL_by_minDelta::Float64
    corrDetecRate_KL_by_95quantile::Float64

    corrDetecRate_GAP::Float64

    VALID_DF::DataFrame
    error_DF::DataFrame

    locID_DF::Array{DataFrame,1}
    evo_DF::Array{DataFrame,1}

    function EvalSR(SEG::SuperResolution.SegmentSR,RESTR_SEG::Array{Tuple{Int64,Int64},1},modINPUT=5)#RESTR_SEG::Int64)
        obj=new()
        obj._selectIDX=RESTR_SEG
        obj._modus=modINPUT
        print("*** modINPUT = $(modINPUT) ***\n")
        NUM_SEG=length(RESTR_SEG)
        #obj._selectIDX=RESTR_SEG[1:RESTR_SEG]
        obj.wModel=Array{SuperResolution.ModelSR,1}(undef,NUM_SEG)
        obj.RefModel=Array{SuperResolution.ModelSR,1}(undef,NUM_SEG)
        obj.distKL=Array{SuperResolution.DistanceSR,1}(undef,NUM_SEG)
        obj.wkldSR=Array{SuperResolution.clusWKLD,1}(undef,NUM_SEG)
        obj.gapSR=Array{SuperResolution.clusGAP,1}(undef,NUM_SEG)

        obj.diff_WK_max_idx=Array{Int64,1}(undef,NUM_SEG)
        obj.diff_GAP_max_idx=Array{Int64,1}(undef,NUM_SEG)

        obj.K_ground_truth=Array{Int64,1}(undef,NUM_SEG)
        obj.N_observations=Array{Int64,1}(undef,NUM_SEG)

        select_mode=[   :K_KL_estimated_by_mean,                    #   modINPUT=1
                        :K_KL_estimated_by_median,                  #   modINPUT=2
                        :K_KL_estimated_by_max,                     #   modINPUT=3
                        :K_KL_estimated_by_minDelta,                #   modINPUT=4
                        :K_KL_estimated_by_95quantile]              #   modINPUT=5

        @time(for i=1:NUM_SEG
                print("\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n")
                print(">>> Starting Segment $(i)/$(NUM_SEG) with index $(SEG.idx[i])\n")
                print("\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n")
                print("(wModel:)\n")
                obj.wModel[i]=SuperResolution.ModelSR(SEG,RESTR_SEG[i],1)     #for segment create all models  WKLD
                # print("(RefModel:)\n")
                #obj.RefModel[i]=SuperResolution.ModelSR(SEG,RESTR_SEG[i],2)   #for segment create all models KLD
                # obj.RefModel[i]=obj.wModel[i]   #temporary workaround to save computation time

                obj.wModel[i].wkld[1:(end-1)].=obj.wModel[i].wkld[2:end]
                obj.wModel[i].weight[1:(end-1)].=obj.wModel[i].weight[2:end]
                obj.wModel[i].geom[1:(end-1)].=obj.wModel[i].geom[2:end]

                obj.wkldSR[i]=SuperResolution.clusWKLD(obj.wModel[i],1)    #cluster WKLD
                # obj.gapSR[i]=SuperResolution.clusGAP(obj.wModel[i],obj.RefModel[i],1)  #Gap statistic
        end)
        obj.K_ground_truth=map(i->SEG.K[CartesianIndex(RESTR_SEG[i])],1:NUM_SEG)
        obj.N_observations=map(i->SEG.N[CartesianIndex(RESTR_SEG[i])],1:NUM_SEG)

        obj.max_mean_idx=map(i->obj.wkldSR[i].max_mean_idx[1],1:length(obj.wkldSR))
        obj.max_median_idx=map(i->obj.wkldSR[i].max_median_idx[1],1:length(obj.wkldSR))
        obj.max_max_idx=map(i->obj.wkldSR[i].max_max_idx[1],1:length(obj.wkldSR))
        obj.max_minDelta_idx=map(i->obj.wkldSR[i].max_minDelta_idx[1],1:length(obj.wkldSR))
        # obj.diff_WK_max_idx=map(i->obj.gapSR[i].diff_WK_max_idx[1],1:length(obj.gapSR))

        #obj.K_KL_estimated_by_mean=map(i->obj.wModel[i].K[CartesianIndex(RESTR_SEG[i])],1:NUM_SEG)
        obj.K_KL_estimated_by_mean=map(i->size(obj.wModel[i].wkld[obj.max_mean_idx[i]],1)-1,1:NUM_SEG)
        obj.K_KL_estimated_by_median=map(i->size(obj.wModel[i].wkld[obj.max_median_idx[i]],1)-1,1:NUM_SEG)
        obj.K_KL_estimated_by_max=map(i->size(obj.wModel[i].wkld[obj.max_max_idx[i]],1)-1,1:NUM_SEG)
        obj.K_KL_estimated_by_minDelta=map(i->size(obj.wModel[i].wkld[obj.max_minDelta_idx[i]],1)-1,1:NUM_SEG)
        obj.K_KL_estimated_by_95quantile=map(i->length( obj.wkldSR[i].wkldMinDelta[ .!(obj.wkldSR[i].wkldMinDelta .< quantile(obj.wkldSR[i].wkldMinDelta, 0.95))] )+1,1:length(obj.wkldSR))  # +1 to account for indexing

        # obj.K_KL_estimated_by_95quantile=map(i->length( obj.wkldSR[i].wkldMin[ .!(obj.wkldSR[i].wkldMin .< quantile(obj.wkldSR[i].wkldMin, 0.95))] )+1,1:length(obj.wkldSR))  # +2 to account for indexing and +1 to account for floor quantile to closest existing wlkd value lower than the quantile
        # obj.K_GAP_estimated=map(i->size(obj.wModel[i].wkld[obj.diff_WK_max_idx[i]],1),1:NUM_SEG)
        #obj.corrDetecRate_GAP=sum(abs.(obj.K_ground_truth-obj.diff_WK_max_idx))/sum(obj.K_ground_truth)
        obj.corrDetecRate_KL_by_mean=mean(exp.((-abs.(obj.K_ground_truth-obj.K_KL_estimated_by_mean))./obj.K_ground_truth))
        obj.corrDetecRate_KL_by_median=mean(exp.((-abs.(obj.K_ground_truth-obj.K_KL_estimated_by_median))./obj.K_ground_truth))
        obj.corrDetecRate_KL_by_max=mean(exp.((-abs.(obj.K_ground_truth-obj.K_KL_estimated_by_max))./obj.K_ground_truth))
        obj.corrDetecRate_KL_by_minDelta=mean(exp.((-abs.(obj.K_ground_truth-obj.K_KL_estimated_by_minDelta))./obj.K_ground_truth))
        obj.corrDetecRate_KL_by_95quantile=mean(exp.((-abs.(obj.K_ground_truth-obj.K_KL_estimated_by_95quantile))./obj.K_ground_truth))
        # obj.corrDetecRate_GAP=mean(exp.((-abs.(obj.K_ground_truth-obj.K_GAP_estimated))./obj.K_ground_truth))
        obj.VALID_DF=DataFrame( idx=obj._selectIDX,
                                N_observations=obj.N_observations,
                                K_ground_truth=obj.K_ground_truth,
                                K_KL_estimated_by_mean=obj.K_KL_estimated_by_mean,
                                K_KL_estimated_by_median=obj.K_KL_estimated_by_median,
                                K_KL_estimated_by_max=obj.K_KL_estimated_by_max,
                                K_KL_estimated_by_minDelta=obj.K_KL_estimated_by_minDelta,
                                K_KL_estimated_by_95quantile=obj.K_KL_estimated_by_95quantile,
                                corrDetecRate_KL_by_mean=obj.corrDetecRate_KL_by_mean,
                                corrDetecRate_KL_by_median=obj.corrDetecRate_KL_by_median,
                                corrDetecRate_KL_by_max=obj.corrDetecRate_KL_by_max,
                                corrDetecRate_KL_by_minDelta=obj.corrDetecRate_KL_by_minDelta,
                                corrDetecRate_KL_by_95quantile=obj.corrDetecRate_KL_by_95quantile,
                                dev_by_mean=(obj.K_ground_truth-obj.K_KL_estimated_by_mean),
                                dev_by_median=(obj.K_ground_truth-obj.K_KL_estimated_by_median),
                                dev_by_max=(obj.K_ground_truth-obj.K_KL_estimated_by_max),
                                dev_by_minDelta=(obj.K_ground_truth-obj.K_KL_estimated_by_minDelta),
                                dev_by_95quantile=(obj.K_ground_truth-obj.K_KL_estimated_by_95quantile) )

        print("\n<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n")
        print("<<< K_ground_truth $(obj.K_ground_truth)\n")
        print("<<< N_observations $(obj.N_observations)\n")
        print("------------------------------------------------\n")
        # print("<<< K_KL_estimated_by_mean $(obj.K_KL_estimated_by_mean)\n")
        # print("<<< KL(mean) correct detection rate $(obj.corrDetecRate_KL_by_mean)\n")
        # print("<<< deviation from ground truth $(obj.K_ground_truth-obj.K_KL_estimated_by_mean)\n")
        # print("------------------------------------------------\n")
        print("<<< K_KL_estimated_by_median $(obj.K_KL_estimated_by_median)\n")
        print("<<< KL(median) correct detection rate $(obj.corrDetecRate_KL_by_median)\n")
        print("<<< deviation from ground truth $(obj.K_ground_truth-obj.K_KL_estimated_by_median)\n")
        print("------------------------------------------------\n")
        # print("<<< K_KL_estimated_by_max $(obj.K_KL_estimated_by_max)\n")
        # print("<<< KL(max) correct detection rate $(obj.corrDetecRate_KL_by_max)\n")
        # print("<<< deviation from ground truth $(obj.K_ground_truth-obj.K_KL_estimated_by_max)\n")
        # print("------------------------------------------------\n")
        print("<<< K_KL_estimated_by_minDelta $(obj.K_KL_estimated_by_minDelta)\n")
        print("<<< KL(minDelta) correct detection rate $(obj.corrDetecRate_KL_by_minDelta)\n")
        print("<<< deviation from ground truth $(obj.K_ground_truth-obj.K_KL_estimated_by_minDelta)\n")
        print("------------------------------------------------\n")
        print("<<< K_KL_estimated_by_95quantile $(obj.K_KL_estimated_by_95quantile)\n")
        print("<<< KL(95quantile) correct detection rate $(obj.corrDetecRate_KL_by_95quantile)\n")
        print("<<< deviation from ground truth $(obj.K_ground_truth-obj.K_KL_estimated_by_95quantile)\n")
        print("------------------------------------------------\n")
        # print("<<< K_GAP_estimated $(obj.K_GAP_estimated)\n")
        # print("<<< GAP correct detection rate $(obj.corrDetecRate_GAP)\n")
        # print("<<< deviation from ground truth $(obj.K_ground_truth-obj.K_GAP_estimated)\n")
        # print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n")
        obj.locID_DF=map(k->DataFrame(  mol_ID=obj.wModel[k].mol_ID,
                                        est_ID=obj.wModel[k].member[end-getproperty(obj, select_mode[modINPUT])[k]+1],
                                        x=obj.wModel[k].loc[1][:,1],
                                        y=obj.wModel[k].loc[1][:,2],
                                        # xC=vcat(map(i->repeat([obj.wModel[k].loc[end-obj.K_KL_estimated_by_95quantile[k]+1][i,:][1]],sum(obj.wModel[k].member[end-obj.K_KL_estimated_by_95quantile[k]+1].==unique(obj.wModel[k].member[end-obj.K_KL_estimated_by_95quantile[k]+1])[i])),1:length(unique(obj.wModel[k].member[end-obj.K_KL_estimated_by_95quantile[k]+1])))...),
                                        # yC=vcat(map(i->repeat([obj.wModel[k].loc[end-obj.K_KL_estimated_by_95quantile[k]+1][i,:][2]],sum(obj.wModel[k].member[end-obj.K_KL_estimated_by_95quantile[k]+1].==unique(obj.wModel[k].member[end-obj.K_KL_estimated_by_95quantile[k]+1])[i])),1:length(unique(obj.wModel[k].member[end-obj.K_KL_estimated_by_95quantile[k]+1])))...),
                                        x_clus_ID=vcat(map(i->obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,1][(1:length(obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,1]))[obj.wModel[k].clus_ID[end-getproperty(obj, select_mode[modINPUT])[k]+1].==i]],obj.wModel[k].member[end-getproperty(obj, select_mode[modINPUT])[k]+1])...),
                                        y_clus_ID=vcat(map(i->obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,2][(1:length(obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,2]))[obj.wModel[k].clus_ID[end-getproperty(obj, select_mode[modINPUT])[k]+1].==i]],obj.wModel[k].member[end-getproperty(obj, select_mode[modINPUT])[k]+1])...),
                                        x_mol_ID=vcat(map(i->mean(obj.wModel[k].loc[1][:,1][obj.wModel[k].mol_ID.==i]),obj.wModel[k].mol_ID)...),
                                        y_mol_ID=vcat(map(i->mean(obj.wModel[k].loc[1][:,2][obj.wModel[k].mol_ID.==i]),obj.wModel[k].mol_ID)...),
                                        intra_dist_clus_ID=map(j->pairwise(Euclidean(), transpose(hcat( obj.wModel[k].loc[1][:,1],obj.wModel[k].loc[1][:,2])),
                                                                                        transpose(hcat( vcat(map(i->obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,1][(1:length(obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,1]))[obj.wModel[k].clus_ID[end-getproperty(obj, select_mode[modINPUT])[k]+1].==i]],obj.wModel[k].member[end-getproperty(obj, select_mode[modINPUT])[k]+1])...),
                                                                                                        vcat(map(i->obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,2][(1:length(obj.wModel[k].loc[end-getproperty(obj, select_mode[modINPUT])[k]+1][:,2]))[obj.wModel[k].clus_ID[end-getproperty(obj, select_mode[modINPUT])[k]+1].==i]],obj.wModel[k].member[end-getproperty(obj, select_mode[modINPUT])[k]+1])...))),
                                                        dims=2)[j,j],1:length(obj.wModel[k].loc[1][:,1])),
                                        ),1:length(obj._selectIDX))

        obj.evo_DF=map(k->DataFrame(    K=map(i->length(unique(obj.wModel[k].member[i])),1:length(obj.wModel[k].member)),
                                        median_intra_dist_clus_ID=map(l->
                                                                    median(map(j->
                                                                        pairwise(Euclidean(), transpose(hcat( obj.wModel[k].loc[1][:,1],obj.wModel[k].loc[1][:,2])),
                                                                                        transpose(hcat( vcat(map(i->obj.wModel[k].loc[l][:,1][(1:length(obj.wModel[k].loc[l][:,1]))[obj.wModel[k].clus_ID[l].==i]],obj.wModel[k].member[l])...),
                                                                                                        vcat(map(i->obj.wModel[k].loc[l][:,2][(1:length(obj.wModel[k].loc[l][:,2]))[obj.wModel[k].clus_ID[l].==i]],obj.wModel[k].member[l])...))),
                                                                                dims=2)[j,j],
                                                                            1:length(obj.wModel[k].loc[1][:,1]))),
                                                                    1:length(obj.wModel[k].loc))
                                        ),1:length(obj._selectIDX))

        obj.error_DF=DataFrame(     idx=obj._selectIDX,
                                    clus_vs_mol_MSE=map(k->(1/length(obj.locID_DF[k].x_clus_ID)) * sum( (obj.locID_DF[k].x_clus_ID-obj.locID_DF[k].x_mol_ID).^2 + (obj.locID_DF[k].y_clus_ID-obj.locID_DF[k].y_mol_ID).^2 ),1:length(obj.locID_DF)),
                                    clus_vs_mol_MAE=map(k->(1/length(obj.locID_DF[k].x_clus_ID)) * sum(sqrt.( (obj.locID_DF[k].x_clus_ID-obj.locID_DF[k].x_mol_ID).^2 + (obj.locID_DF[k].y_clus_ID-obj.locID_DF[k].y_mol_ID).^2 )),1:length(obj.locID_DF))
        )
        print("------------------------------------------------\n")
        print("<<< mean clus_vs_mol_MSE $(mean(obj.error_DF.clus_vs_mol_MSE))\n")
        print("<<< median clus_vs_mol_MSE $(median(obj.error_DF.clus_vs_mol_MSE))\n")
        print("<<< sample standard deviation $(Statistics.std(obj.error_DF.clus_vs_mol_MSE))\n")
        print("------------------------------------------------\n")
        return obj
    end
end
export EvalSR

# map(i-> EVAL.wModel[1].loc[end-EVAL.K_KL_estimated_by_95quantile[1]+1][:,1][collect(1:length(EVAL.wModel[1].loc[end-EVAL.K_KL_estimated_by_95quantile[1]+1][:,1]))[EVAL.wModel[1].clus_ID[end-EVAL.K_KL_estimated_by_95quantile[1]+1].==i]],EVAL.wModel[1].member[end-EVAL.K_KL_estimated_by_95quantile[1]+1])
