from __future__ import absolute_import, division, print_function
from cctbx.uctbx import unit_cell
from cctbx import miller, crystal, statistics
from cctbx.array_family import flex
from iotbx import mtz
from libtbx.utils import Sorry
import math, os
import numpy as np
from copy import deepcopy
from six.moves import cPickle as pickle
from collections import Counter
from .mod_merge_data import merge_data_handler
from .mod_mx import mx_handler
from .mod_leastsqr import good_unit_cell
from six.moves import range


class intensities_scaler(object):
    """
    Author      : Uervirojnangkoorn, M.
    Created     : 7/13/2014
    Merge equivalent reflections and report intensity and refinement statistics.
    """

    def __init__(self):
        """Constructor."""
        self.CONST_SE_MIN_WEIGHT = 0.17
        self.CONST_SE_MAX_WEIGHT = 1.0
        self.CONST_SIG_I_FACTOR = 1.5

    def write_stat_pickle(self, iparams, stat_dict):
        fname = iparams.run_no + "/stats/pickle_" + str(os.getpid()) + ".stat"
        if os.path.isfile(fname):
            pickle_stat = pickle.load(open(fname, "rb"))
            for key in stat_dict:
                if key in pickle_stat:
                    pickle_stat[key].append(stat_dict[key][0])
                else:
                    pickle_stat[key] = stat_dict[key]
            pickle.dump(pickle_stat, open(fname, "wb"))
        else:
            pickle.dump(stat_dict, open(fname, "wb"))

    def calc_avg_I_cpp(self, prep_output, iparams, avg_mode):
        group_no, group_id_list, miller_index, miller_indices_ori, I, sigI, G, B, p_set, rs_set, wavelength_set, sin_theta_over_lambda_sq, SE, uc_mean, wavelength_mean, pickle_filename_set, txt_out = (
            prep_output
        )
        from prime import Average_Mode, averaging_engine

        if avg_mode == "average":
            avg_mode_cpp = Average_Mode.Average
        elif avg_mode == "weighted":
            avg_mode_cpp = Average_Mode.Weighted
        elif avg_mode == "final":
            avg_mode_cpp = Average_Mode.Final
        else:
            raise Sorry("Bad averaging mode selected: %s" % avg_mode)
        sigma_max = iparams.sigma_rejection
        engine = averaging_engine(
            group_no,
            group_id_list,
            miller_index,
            miller_indices_ori,
            I,
            sigI,
            G,
            B,
            p_set,
            rs_set,
            wavelength_set,
            sin_theta_over_lambda_sq,
            SE,
            pickle_filename_set,
        )
        engine.avg_mode = avg_mode_cpp
        engine.sigma_max = sigma_max
        engine.flag_volume_correction = iparams.flag_volume_correction
        engine.n_rejection_cycle = iparams.n_rejection_cycle
        engine.flag_output_verbose = iparams.flag_output_verbose
        results = engine.calc_avg_I()
        mdh = merge_data_handler()
        mdh.extend_data(
            results.miller_index,
            results.I_avg,
            results.sigI_avg,
            (results.r_meas_top, results.r_meas_btm, results.multiplicity),
            (
                results.I_avg_even,
                results.I_avg_odd,
                results.I_avg_even_h,
                results.I_avg_odd_h,
                results.I_avg_even_k,
                results.I_avg_odd_k,
                results.I_avg_even_l,
                results.I_avg_odd_l,
            ),
            uc_mean,
            wavelength_mean,
        )
        return mdh, results.txt_obs_out, results.txt_reject_out

    def calc_mean_unit_cell(self, results):
        uc_array = [list(pres.uc_params) for pres in results if pres is not None]
        return np.mean(uc_array, 0), np.median(uc_array, 0), np.std(uc_array, 0)

    def calc_mean_postref_parameters(self, results):
        params_array = [
            [
                pres.G,
                pres.B,
                pres.ry,
                pres.rz,
                pres.re,
                pres.r0,
                pres.voigt_nu,
                pres.rotx,
                pres.roty,
                pres.R_final,
                pres.R_xy_final,
                pres.SE,
            ]
            for pres in results
            if (
                pres is not None
                and not math.isnan(pres.G)
                and not math.isnan(pres.B)
                and not math.isnan(pres.ry)
                and not math.isnan(pres.rz)
                and not math.isnan(pres.re)
                and not math.isnan(pres.r0)
                and not math.isnan(pres.voigt_nu)
                and not math.isnan(pres.rotx)
                and not math.isnan(pres.roty)
                and not math.isnan(pres.R_final)
                and not math.isnan(pres.R_xy_final)
                and not math.isnan(pres.SE)
            )
        ]
        return (
            np.mean(params_array, 0),
            np.median(params_array, 0),
            np.std(params_array, 0),
        )

    def prepare_output(self, results, iparams, avg_mode):
        if avg_mode == "average":
            cc_thres = 0
        else:
            cc_thres = iparams.frame_accept_min_cc
        std_filter = iparams.sigma_rejection
        if iparams.flag_weak_anomalous:
            if avg_mode == "final":
                target_anomalous_flag = iparams.target_anomalous_flag
            else:
                target_anomalous_flag = False
        else:
            target_anomalous_flag = iparams.target_anomalous_flag
        pr_params_mean, pr_params_med, pr_params_std = self.calc_mean_postref_parameters(
            results
        )
        G_mean, B_mean, ry_mean, rz_mean, re_mean, r0_mean, voigt_nu_mean, rotx_mean, roty_mean, R_mean, R_xy_mean, SE_mean = (
            pr_params_mean
        )
        G_med, B_med, ry_med, rz_med, re_med, r0_med, voigt_nu_med, rotx_med, roty_med, R_med, R_xy_med, SE_med = (
            pr_params_med
        )
        G_std, B_std, ry_std, rz_std, re_std, r0_std, voigt_nu_std, rotx_std, roty_std, R_std, R_xy_std, SE_std = (
            pr_params_std
        )
        # prepare data for merging
        miller_indices_all = flex.miller_index()
        miller_indices_ori_all = flex.miller_index()
        I_all = flex.double()
        sigI_all = flex.double()
        G_all = flex.double()
        B_all = flex.double()
        p_all = flex.double()
        rx_all = flex.double()
        rs_all = flex.double()
        rh_all = flex.double()
        SE_all = flex.double()
        sin_sq_all = flex.double()
        wavelength_all = flex.double()
        detector_distance_set = flex.double()
        R_init_all = flex.double()
        R_final_all = flex.double()
        R_xy_init_all = flex.double()
        R_xy_final_all = flex.double()
        pickle_filename_all = flex.std_string()
        filtered_results = []
        cn_good_frame, cn_bad_frame_SE, cn_bad_frame_uc, cn_bad_frame_cc, cn_bad_frame_G, cn_bad_frame_re = (
            0,
            0,
            0,
            0,
            0,
            0,
        )
        crystal_orientation_dict = {}
        for pres in results:
            if pres is not None:
                pickle_filepath = pres.pickle_filename.split("/")
                img_filename = pickle_filepath[len(pickle_filepath) - 1]
                flag_pres_ok = True
                # check SE, CC, UC, G, B, gamma_e
                if math.isnan(pres.G):
                    flag_pres_ok = False
                if math.isnan(pres.SE) or np.isinf(pres.SE):
                    flag_pres_ok = False
                if flag_pres_ok and SE_std > 0:
                    if abs(pres.SE - SE_med) / SE_std > std_filter:
                        flag_pres_ok = False
                        cn_bad_frame_SE += 1
                if flag_pres_ok and pres.CC_final < cc_thres:
                    flag_pres_ok = False
                    cn_bad_frame_cc += 1
                if flag_pres_ok:
                    if G_std > 0:
                        if abs(pres.G - G_med) / G_std > std_filter:
                            flag_pres_ok = False
                            cn_bad_frame_G += 1
                if flag_pres_ok:
                    if re_std > 0:
                        if abs(pres.re - re_med) / re_std > std_filter:
                            flag_pres_ok = False
                            cn_bad_frame_re += 1
                if flag_pres_ok and not good_unit_cell(
                    pres.uc_params, iparams, iparams.merge.uc_tolerance
                ):
                    flag_pres_ok = False
                    cn_bad_frame_uc += 1
                data_size = pres.observations.size()
                if flag_pres_ok:
                    cn_good_frame += 1
                    filtered_results.append(pres)
                    R_init_all.append(pres.R_init)
                    R_final_all.append(pres.R_final)
                    R_xy_init_all.append(pres.R_xy_init)
                    R_xy_final_all.append(pres.R_xy_final)
                    miller_indices_all.extend(pres.observations.indices())
                    miller_indices_ori_all.extend(pres.observations_original.indices())
                    I_all.extend(pres.observations.data())
                    sigI_all.extend(pres.observations.sigmas())
                    G_all.extend(flex.double([pres.G] * data_size))
                    B_all.extend(flex.double([pres.B] * data_size))
                    p_all.extend(pres.partiality)
                    rs_all.extend(pres.rs_set)
                    rh_all.extend(pres.rh_set)
                    sin_sq_all.extend(
                        pres.observations.two_theta(wavelength=pres.wavelength)
                        .sin_theta_over_lambda_sq()
                        .data()
                    )
                    SE_all.extend(flex.double([pres.SE] * data_size))
                    wavelength_all.extend(flex.double([pres.wavelength] * data_size))
                    detector_distance_set.append(pres.detector_distance_mm)
                    pickle_filename_all.extend(
                        flex.std_string([pres.pickle_filename] * data_size)
                    )
                    crystal_orientation_dict[
                        pres.pickle_filename
                    ] = pres.crystal_orientation
        # plot stats
        self.plot_stats(filtered_results, iparams)
        # write out updated crystal orientation as a pickle file
        if not iparams.flag_hush:
            pickle.dump(
                crystal_orientation_dict,
                open(iparams.run_no + "/" + "crystal.o", "wb"),
                pickle.HIGHEST_PROTOCOL,
            )
        # calculate average unit cell
        uc_mean, uc_med, uc_std = self.calc_mean_unit_cell(filtered_results)
        unit_cell_mean = unit_cell(tuple(uc_mean))
        # recalculate stats for pr parameters
        pr_params_mean, pr_params_med, pr_params_std = self.calc_mean_postref_parameters(
            filtered_results
        )
        G_mean, B_mean, ry_mean, rz_mean, re_mean, r0_mean, voigt_nu_mean, rotx_mean, roty_mean, R_mean, R_xy_mean, SE_mean = (
            pr_params_mean
        )
        G_med, B_med, ry_med, rz_med, re_med, r0_med, voigt_nu_med, rotx_med, roty_med, R_med, R_xy_med, SE_med = (
            pr_params_med
        )
        G_std, B_std, ry_std, rz_std, re_std, r0_std, voigt_nu_std, rotx_std, roty_std, R_std, R_xy_std, SE_std = (
            pr_params_std
        )
        # from all observations merge them
        crystal_symmetry = crystal.symmetry(
            unit_cell=tuple(uc_mean), space_group_symbol=str(iparams.target_space_group)
        )
        miller_set_all = miller.set(
            crystal_symmetry=crystal_symmetry,
            indices=miller_indices_all,
            anomalous_flag=target_anomalous_flag,
        )
        miller_array_all = miller_set_all.array(
            data=I_all, sigmas=sigI_all
        ).set_observation_type_xray_intensity()
        # sort reflections according to asymmetric-unit symmetry hkl
        perm = miller_array_all.sort_permutation(by_value="packed_indices")
        miller_indices_all_sort = miller_array_all.indices().select(perm)
        miller_indices_ori_all_sort = miller_indices_ori_all.select(perm)
        I_obs_all_sort = miller_array_all.data().select(perm)
        sigI_obs_all_sort = miller_array_all.sigmas().select(perm)
        G_all_sort = G_all.select(perm)
        B_all_sort = B_all.select(perm)
        p_all_sort = p_all.select(perm)
        rs_all_sort = rs_all.select(perm)
        wavelength_all_sort = wavelength_all.select(perm)
        sin_sq_all_sort = sin_sq_all.select(perm)
        SE_all_sort = SE_all.select(perm)
        pickle_filename_all_sort = pickle_filename_all.select(perm)
        miller_array_uniq = (
            miller_array_all.merge_equivalents()
            .array()
            .complete_array(d_min=iparams.merge.d_min, d_max=iparams.merge.d_max)
        )
        matches_uniq = miller.match_multi_indices(
            miller_indices_unique=miller_array_uniq.indices(),
            miller_indices=miller_indices_all_sort,
        )
        pair_0 = flex.int([pair[0] for pair in matches_uniq.pairs()])
        pair_1 = flex.int([pair[1] for pair in matches_uniq.pairs()])
        group_id_list = flex.int(
            [pair_0[pair_1[i]] for i in range(len(matches_uniq.pairs()))]
        )
        tally = Counter()
        for elem in group_id_list:
            tally[elem] += 1
        cn_group = len(tally)
        # preparte txt out stat
        txt_out = "Summary of refinement and merging\n"
        txt_out += " No. good frames:          %12.0f\n" % (cn_good_frame)
        txt_out += " No. bad cc frames:        %12.0f\n" % (cn_bad_frame_cc)
        txt_out += " No. bad G frames) :       %12.0f\n" % (cn_bad_frame_G)
        txt_out += " No. bad unit cell frames: %12.0f\n" % (cn_bad_frame_uc)
        txt_out += " No. bad gamma_e frames:   %12.0f\n" % (cn_bad_frame_re)
        txt_out += " No. bad SE:               %12.0f\n" % (cn_bad_frame_SE)
        txt_out += " No. observations:         %12.0f\n" % (len(I_obs_all_sort))
        txt_out += "Mean target value (BEFORE: Mean Median (Std.))\n"
        txt_out += " post-refinement:          %12.2f %12.2f (%9.2f)\n" % (
            np.mean(R_init_all),
            np.median(R_init_all),
            np.std(R_init_all),
        )
        txt_out += " (x,y) restraints:         %12.2f %12.2f (%9.2f)\n" % (
            np.mean(R_xy_init_all),
            np.median(R_xy_init_all),
            np.std(R_xy_init_all),
        )
        txt_out += "Mean target value (AFTER: Mean Median (Std.))\n"
        txt_out += " post-refinement:          %12.2f %12.2f (%9.2f)\n" % (
            np.mean(R_final_all),
            np.median(R_final_all),
            np.std(R_final_all),
        )
        txt_out += " (x,y) restraints:         %12.2f %12.2f (%9.2f)\n" % (
            np.mean(R_xy_final_all),
            np.median(R_xy_final_all),
            np.std(R_xy_final_all),
        )
        txt_out += " SE:                       %12.2f %12.2f (%9.2f)\n" % (
            SE_mean,
            SE_med,
            SE_std,
        )
        txt_out += " G:                        %12.3e %12.3e (%9.2e)\n" % (
            G_mean,
            G_med,
            G_std,
        )
        txt_out += " B:                        %12.2f %12.2f (%9.2f)\n" % (
            B_mean,
            B_med,
            B_std,
        )
        txt_out += " Rot.x:                    %12.2f %12.2f (%9.2f)\n" % (
            rotx_mean * 180 / math.pi,
            rotx_med * 180 / math.pi,
            rotx_std * 180 / math.pi,
        )
        txt_out += " Rot.y:                    %12.2f %12.2f (%9.2f)\n" % (
            roty_mean * 180 / math.pi,
            roty_med * 180 / math.pi,
            roty_std * 180 / math.pi,
        )
        txt_out += " gamma_y:                  %12.5f %12.5f (%9.5f)\n" % (
            ry_mean,
            ry_med,
            ry_std,
        )
        txt_out += " gamma_z:                  %12.5f %12.5f (%9.5f)\n" % (
            rz_mean,
            rz_med,
            rz_std,
        )
        txt_out += " gamma_0:                  %12.5f %12.5f (%9.5f)\n" % (
            r0_mean,
            r0_med,
            r0_std,
        )
        txt_out += " gamma_e:                  %12.5f %12.5f (%9.5f)\n" % (
            re_mean,
            re_med,
            re_std,
        )
        txt_out += " voigt_nu:                 %12.5f %12.5f (%9.5f)\n" % (
            voigt_nu_mean,
            voigt_nu_med,
            voigt_nu_std,
        )
        txt_out += " unit cell\n"
        txt_out += "   a:                      %12.2f %12.2f (%9.2f)\n" % (
            uc_mean[0],
            uc_med[0],
            uc_std[0],
        )
        txt_out += "   b:                      %12.2f %12.2f (%9.2f)\n" % (
            uc_mean[1],
            uc_med[1],
            uc_std[1],
        )
        txt_out += "   c:                      %12.2f %12.2f (%9.2f)\n" % (
            uc_mean[2],
            uc_med[2],
            uc_std[2],
        )
        txt_out += "   alpha:                  %12.2f %12.2f (%9.2f)\n" % (
            uc_mean[3],
            uc_med[3],
            uc_std[3],
        )
        txt_out += "   beta:                   %12.2f %12.2f (%9.2f)\n" % (
            uc_mean[4],
            uc_med[4],
            uc_std[4],
        )
        txt_out += "   gamma:                  %12.2f %12.2f (%9.2f)\n" % (
            uc_mean[5],
            uc_med[5],
            uc_std[5],
        )
        txt_out += "Parmeters from integration (not-refined)\n"
        txt_out += "  Wavelength:              %12.5f %12.5f (%9.5f)\n" % (
            np.mean(wavelength_all),
            np.median(wavelength_all),
            np.std(wavelength_all),
        )
        txt_out += "  Detector distance:       %12.5f %12.5f (%9.5f)\n" % (
            np.mean(detector_distance_set),
            np.median(detector_distance_set),
            np.std(detector_distance_set),
        )
        txt_out += "* (standard deviation)\n"
        # write out stat. pickle
        if not iparams.flag_hush:
            stat_dict = {
                "n_frames_good": [cn_good_frame],
                "n_frames_bad_cc": [cn_bad_frame_cc],
                "n_frames_bad_G": [cn_bad_frame_G],
                "n_frames_bad_uc": [cn_bad_frame_uc],
                "n_frames_bad_gamma_e": [cn_bad_frame_re],
                "n_frames_bad_SE": [cn_bad_frame_SE],
                "n_observations": [len(I_obs_all_sort)],
                "R_start": [np.mean(R_init_all)],
                "R_end": [np.mean(R_final_all)],
                "R_xy_start": [np.mean(R_xy_init_all)],
                "R_xy_end": [np.mean(R_xy_final_all)],
                "mean_gamma_y": [ry_mean],
                "std_gamma_y": [ry_std],
                "mean_gamma_z": [rz_mean],
                "std_gamma_z": [rz_std],
                "mean_gamma_0": [r0_mean],
                "std_gamma_0": [r0_std],
                "mean_gamma_e": [re_mean],
                "std_gamma_e": [re_std],
                "mean_voigt_nu": [voigt_nu_mean],
                "std_voigt_nu": [voigt_nu_std],
                "mean_a": [uc_mean[0]],
                "std_a": [uc_std[0]],
                "mean_b": [uc_mean[1]],
                "std_b": [uc_std[1]],
                "mean_c": [uc_mean[2]],
                "std_c": [uc_std[2]],
                "mean_alpha": [uc_mean[3]],
                "std_alpha": [uc_std[3]],
                "mean_beta": [uc_mean[4]],
                "std_beta": [uc_std[4]],
                "mean_gamma": [uc_mean[5]],
                "std_gamma": [uc_std[5]],
            }
            self.write_stat_pickle(iparams, stat_dict)
        return (
            cn_group,
            group_id_list,
            miller_indices_all_sort,
            miller_indices_ori_all_sort,
            I_obs_all_sort,
            sigI_obs_all_sort,
            G_all_sort,
            B_all_sort,
            p_all_sort,
            rs_all_sort,
            wavelength_all_sort,
            sin_sq_all_sort,
            SE_all_sort,
            uc_mean,
            np.mean(wavelength_all),
            pickle_filename_all_sort,
            txt_out,
        )

    def write_output(self, mdh, iparams, output_mtz_file_prefix, avg_mode):
        if iparams.flag_weak_anomalous:
            if avg_mode == "final":
                target_anomalous_flag = iparams.target_anomalous_flag
            else:
                target_anomalous_flag = False
        else:
            target_anomalous_flag = iparams.target_anomalous_flag
        uc_mean = mdh.uc_mean
        wavelength_mean = mdh.wavelength_mean
        # output mtz file and report binning stat
        miller_set_merge = crystal.symmetry(
            unit_cell=unit_cell(tuple(uc_mean)),
            space_group_symbol=str(iparams.target_space_group),
        ).build_miller_set(
            anomalous_flag=target_anomalous_flag, d_min=iparams.merge.d_min
        )
        mdh.generate_miller_array_from_miller_set(
            miller_set_merge, target_anomalous_flag
        )
        miller_array_complete = miller_set_merge.array()
        fake_data = flex.double([1.0] * len(miller_array_complete.indices()))
        miller_array_template_asu = miller_array_complete.customized_copy(
            data=fake_data, sigmas=fake_data
        ).resolution_filter(d_min=iparams.merge.d_min, d_max=iparams.merge.d_max)
        n_refl_all = mdh.get_size()
        # do another resolution filter here
        i_sel_res = mdh.miller_array_merge.resolution_filter_selection(
            d_min=iparams.merge.d_min, d_max=iparams.merge.d_max
        )
        mdh.reduce_by_selection(i_sel_res)
        n_refl_out_resolutions = n_refl_all - mdh.get_size()
        # remove outliers
        sequences = flex.int(range(mdh.get_size()))
        good_sequences = []
        for i_rejection in range(iparams.n_rejection_cycle):
            binner_merge = mdh.miller_array_merge.setup_binner(n_bins=200)
            for i_bin in range(1, 201):
                i_binner = binner_merge.bin_indices() == i_bin
                I_obs_bin = mdh.miller_array_merge.data().select(i_binner)
                sequences_bin = sequences.select(i_binner)
                if len(I_obs_bin) > 0:
                    I_obs_bin = mdh.miller_array_merge.data().select(i_binner)
                    try:
                        i_filter = (
                            flex.abs(
                                (I_obs_bin - np.median(I_obs_bin)) / np.std(I_obs_bin)
                            )
                            < 10
                        )
                    except Exception as e:
                        print(
                            "Warning: outlier rejection by bins failed because of floating point."
                        )
                        print(e)
                        i_filter = flex.bool([True] * len(I_obs_bin))
                    good_sequences.extend(list(sequences_bin.select(i_filter)))
        mdh.reduce_by_selection(flex.size_t(good_sequences))
        n_refl_outliers = n_refl_all - n_refl_out_resolutions - mdh.get_size()
        # get iso if given.
        mxh = mx_handler()
        flag_hklisoin_found, miller_array_iso = mxh.get_miller_array_from_reflection_file(
            iparams.hklisoin
        )
        # write output files
        if output_mtz_file_prefix != "":
            # write as mtz file
            miller_array_merge_unique = (
                mdh.miller_array_merge.merge_equivalents().array()
            )
            info = miller.array_info(wavelength=wavelength_mean)
            miller_array_merge_unique.set_info(info)
            mtz_dataset_merge = miller_array_merge_unique.as_mtz_dataset(
                column_root_label="IOBS"
            )
            mtz_dataset_merge.mtz_object().write(
                file_name=output_mtz_file_prefix + "_merge.mtz"
            )
            # write as cns file
            f_cns = open(output_mtz_file_prefix + "_merge.hkl", "w")
            miller_array_merge_unique.export_as_cns_hkl(file_object=f_cns)
            f_cns.close()
        # calculate merging stat table
        if True:
            # calculate isotropic B-factor
            try:
                mxh = mx_handler()
                asu_contents = mxh.get_asu_contents(iparams.n_residues)
                observations_as_f = mdh.miller_array_merge.as_amplitude_array()
                observations_as_f.setup_binner(auto_binning=True)
                wp = statistics.wilson_plot(
                    observations_as_f, asu_contents, e_statistics=True
                )
                B_merged = wp.wilson_b
            except Exception as e:
                B_merged = 0
                print(
                    "Warning: b-factor calculation in mod_util failed. Reset b-factor to 0"
                )
                print(e)
            # report binning stats
            txt_out = "\n"
            txt_out += "Isotropic B-factor:  %7.2f\n" % (B_merged)
            txt_out += "No. of reflections\n"
            txt_out += " all:                %7.0f\n" % (n_refl_all)
            txt_out += " outside resolution: %7.0f\n" % (n_refl_out_resolutions)
            txt_out += " outliers:           %7.0f\n" % (n_refl_outliers)
            txt_out += " total left:         %7.0f\n" % (mdh.get_size())
            txt_out += "Summary for " + output_mtz_file_prefix + "_merge.mtz\n"
            txt_out += "Bin Resolution Range     Completeness      <N_obs> |Rmerge  Rsplit   CC1/2   N_ind |CCiso   N_ind|CCanoma  N_ind| <I/sigI>   <I>    <sigI>    <I**2>\n"
            txt_out += "--------------------------------------------------------------------------------------------------------------------------------------------------\n"
            # for stat pickle
            sp_res, sp_complete, sp_n_obs, sp_cc12, sp_cc12_anom, sp_rmerge, sp_i_o_sigi, sp_isqr = (
                [],
                [],
                [],
                [],
                [],
                [],
                [],
                [],
            )
            # binning
            binner_template_asu = miller_array_template_asu.setup_binner(
                n_bins=iparams.n_bins
            )
            binner_template_asu_indices = binner_template_asu.bin_indices()
            # for stats on axis cones
            mdh_astar = deepcopy(mdh)
            mdh_bstar = deepcopy(mdh)
            mdh_cstar = deepcopy(mdh)
            mdh_astar.reduce_to_cone_on_axis((1, 0, 0), iparams.percent_cone_fraction)
            mdh_bstar.reduce_to_cone_on_axis((0, 1, 0), iparams.percent_cone_fraction)
            mdh_cstar.reduce_to_cone_on_axis((0, 0, 1), iparams.percent_cone_fraction)
            # prepare text out for axis cones
            txt_out_cone = "Summary of CC1/2 on three crystal axes\n"
            txt_out_cone += "Bin Resolution Range           CC1/2                      <I>                          N_refl           \n"
            txt_out_cone += "                        a*      b*      c*  |      a*        b*       c*    |    a*      b*     c*      \n"
            txt_out_cone += "---------------------------------------------------------------------------------------------------------\n"
            for i in range(1, iparams.n_bins + 1):
                i_binner = binner_template_asu_indices == i
                miller_indices_template_bin = miller_array_template_asu.indices().select(
                    i_binner
                )
                # for all reflections
                mdh_bin = deepcopy(mdh)
                mdh_bin.reduce_by_miller_index(miller_indices_template_bin)
                cc12, n_refl_cc12 = mdh_bin.get_cc12()
                cciso, n_refl_cciso = mdh_bin.get_cciso(miller_array_iso)
                cc_anom_acentric, n_refl_anom_acentric = mdh_bin.get_cc_anom()
                completeness = (
                    mdh_bin.get_size() / len(miller_indices_template_bin)
                ) * 100
                multiplicity = mdh_bin.get_multiplicity()
                txt_out += (
                    "%02d %7.2f - %7.2f %5.1f %6.0f / %6.0f %7.2f %7.2f %7.2f %7.2f %6.0f %7.2f %6.0f %7.2f %6.0f %8.2f %10.1f %8.1f %6.2f\n"
                    % (
                        i,
                        binner_template_asu.bin_d_range(i)[0],
                        binner_template_asu.bin_d_range(i)[1],
                        completeness,
                        mdh_bin.get_size(),
                        len(miller_indices_template_bin),
                        multiplicity,
                        mdh_bin.get_r_meas() * 100,
                        mdh_bin.get_r_split() * 100,
                        cc12 * 100,
                        n_refl_cc12,
                        cciso * 100,
                        n_refl_cciso,
                        cc_anom_acentric,
                        n_refl_anom_acentric,
                        mdh_bin.get_mean_IoversigI(),
                        mdh_bin.get_mean_I(),
                        mdh_bin.get_mean_sigI(),
                        mdh_bin.get_second_moment(),
                    )
                )
                # for reflections on cones
                mdh_astar_bin = deepcopy(mdh_astar)
                mdh_astar_bin.reduce_by_miller_index(miller_indices_template_bin)
                cc12_astar, n_refl_cc12_astar = mdh_astar_bin.get_cc12()
                mdh_bstar_bin = deepcopy(mdh_bstar)
                mdh_bstar_bin.reduce_by_miller_index(miller_indices_template_bin)
                cc12_bstar, n_refl_cc12_bstar = mdh_bstar_bin.get_cc12()
                mdh_cstar_bin = deepcopy(mdh_cstar)
                mdh_cstar_bin.reduce_by_miller_index(miller_indices_template_bin)
                cc12_cstar, n_refl_cc12_cstar = mdh_cstar_bin.get_cc12()
                txt_out_cone += (
                    "%02d %7.2f - %7.2f %7.2f %7.2f %7.2f %10.1f %10.1f %10.1f %6.0f %6.0f %6.0f\n"
                    % (
                        i,
                        binner_template_asu.bin_d_range(i)[0],
                        binner_template_asu.bin_d_range(i)[1],
                        cc12_astar * 100,
                        cc12_bstar * 100,
                        cc12_cstar * 100,
                        mdh_astar_bin.get_mean_I(),
                        mdh_bstar_bin.get_mean_I(),
                        mdh_cstar_bin.get_mean_I(),
                        n_refl_cc12_astar,
                        n_refl_cc12_bstar,
                        n_refl_cc12_cstar,
                    )
                )
                # for stat pickle
                sp_res.append(binner_template_asu.bin_d_range(i)[1])
                sp_complete.append(completeness)
                sp_n_obs.append(multiplicity)
                sp_cc12.append(cc12)
                sp_cc12_anom.append(cc_anom_acentric)
                sp_rmerge.append(mdh_bin.get_r_meas() * 100)
                sp_i_o_sigi.append(mdh_bin.get_mean_IoversigI())
                sp_isqr.append(mdh.get_second_moment())
            # txt out total for all reflections
            cc12, n_refl_cc12 = mdh.get_cc12()
            cciso, n_refl_cciso = mdh.get_cciso(miller_array_iso)
            cc_anom_acentric, n_refl_anom_acentric = mdh.get_cc_anom()
            txt_out += "--------------------------------------------------------------------------------------------------------------------------------------------------\n"
            txt_out += (
                "        TOTAL        %5.1f %6.0f / %6.0f %7.2f %7.2f %7.2f %7.2f %6.0f %7.2f %6.0f %7.2f %6.0f %8.2f %10.1f %8.1f %6.2f\n"
                % (
                    (mdh.get_size() / miller_array_template_asu.size()) * 100,
                    mdh.get_size(),
                    miller_array_template_asu.size(),
                    mdh.get_multiplicity(),
                    mdh.get_r_meas() * 100,
                    mdh.get_r_split() * 100,
                    cc12 * 100,
                    n_refl_cc12,
                    cciso * 100,
                    n_refl_cciso,
                    cc_anom_acentric,
                    n_refl_anom_acentric,
                    mdh.get_mean_IoversigI(),
                    mdh.get_mean_I(),
                    mdh.get_mean_sigI(),
                    mdh.get_second_moment(),
                )
            )
            txt_out += "--------------------------------------------------------------------------------------------------------------------------------------------------\n"
            txt_out += "\n"
            # txt out total for reflections on cones
            cc12_astar, n_refl_cc12_astar = mdh_astar.get_cc12()
            cc12_bstar, n_refl_cc12_bstar = mdh_bstar.get_cc12()
            cc12_cstar, n_refl_cc12_cstar = mdh_cstar.get_cc12()
            txt_out_cone += "----------------------------------------------------------------------------------------------------------\n"
            txt_out_cone += (
                "       total         %7.2f %7.2f %7.2f %10.1f %10.1f %10.1f %6.0f %6.0f %6.0f\n"
                % (
                    cc12_astar * 100,
                    cc12_bstar * 100,
                    cc12_cstar * 100,
                    mdh_astar.get_mean_I(),
                    mdh_bstar.get_mean_I(),
                    mdh_cstar.get_mean_I(),
                    n_refl_cc12_astar,
                    n_refl_cc12_bstar,
                    n_refl_cc12_cstar,
                )
            )
            txt_out_cone += "----------------------------------------------------------------------------------------------------------\n"
            txt_out_cone += "\n"
            txt_out_table1 = "Table1 (" + avg_mode + ")\n"
            txt_out_table1 += (
                "  Space group: "
                + str(mdh.miller_array_merge.space_group_info())
                + "\n"
            )
            txt_out_table1 += (
                "  Cell dimensions: %6.2f, %6.2f, %6.2f, %6.2f, %6.2f, %6.2f\n"
                % tuple(mdh.uc_mean)
            )
            txt_out_table1 += "  Resolution (A): %6.2f - %6.2f (%6.2f - %6.2f)\n" % (
                mdh.miller_array_merge.d_max_min()[0],
                mdh.miller_array_merge.d_max_min()[1],
                sp_res[-2],
                sp_res[-1],
            )
            txt_out_table1 += "  Rmerge: %6.2f (%6.2f)\n" % (
                mdh.get_r_meas() * 100,
                sp_rmerge[-1],
            )
            txt_out_table1 += "  CC1/2: %6.2f (%6.2f)\n" % (
                mdh.get_cc12()[0] * 100,
                sp_cc12[-1],
            )
            txt_out_table1 += "  I/sigI: %6.2f (%6.2f)\n" % (
                mdh.get_mean_IoversigI(),
                sp_i_o_sigi[-1],
            )
            txt_out_table1 += "  Completeness (%%): %6.2f (%6.2f)\n" % (
                (mdh.get_size() / miller_array_template_asu.size()) * 100,
                sp_complete[-1],
            )
            txt_out_table1 += "  Redundancy: %6.2f (%6.2f)\n" % (
                mdh.get_multiplicity(),
                sp_n_obs[-1],
            )
            # save data for stat. pickle in stat_dict
            if not iparams.flag_hush:
                stat_dict = {
                    "binned_resolution": [sp_res],
                    "binned_completeness": [sp_complete],
                    "binned_n_obs": [sp_n_obs],
                    "binned_cc12": [sp_cc12],
                    "binned_cc12_anom": [sp_cc12_anom],
                    "binned_rmerge": [sp_rmerge],
                    "binned_i_o_sigi": [sp_i_o_sigi],
                    "binned_isqr": [sp_isqr],
                    "total_res_max": [mdh.miller_array_merge.d_max_min()[0]],
                    "total_res_min": [mdh.miller_array_merge.d_max_min()[1]],
                    "total_completeness": [
                        (mdh.get_size() / miller_array_template_asu.size()) * 100
                    ],
                    "total_n_obs": [mdh.get_multiplicity()],
                    "total_cc12": [mdh.get_cc12()[0] * 100],
                    "total_rmerge": [mdh.get_r_meas() * 100],
                    "total_i_o_sigi": [mdh.get_mean_IoversigI()],
                    "space_group_info": [mdh.miller_array_merge.space_group_info()],
                }
                self.write_stat_pickle(iparams, stat_dict)
            txt_out += txt_out_cone + txt_out_table1
        return mdh, txt_out

    def plot_stats(self, results, iparams):
        # retrieve stats from results and plot them
        if iparams.flag_plot or iparams.flag_output_verbose:
            # for plotting set n_bins = 5 to avoid empty bin
            n_bins_plot = 5
            # get expected f^2
            try:
                mxh = mx_handler()
                asu_contents = mxh.get_asu_contents(iparams.n_residues)
                observations_as_f = results[0].observations.as_amplitude_array()
                binner_template_asu = observations_as_f.setup_binner(n_bins=n_bins_plot)
                wp = statistics.wilson_plot(
                    observations_as_f, asu_contents, e_statistics=True
                )
                expected_f_sq = wp.expected_f_sq
                mean_stol_sq = wp.mean_stol_sq
            except Exception:
                expected_f_sq = flex.double([0] * n_bins_plot)
                mean_stol_sq = flex.double(range(n_bins_plot))
                print("Warning: Wilson plot calculation in plot stats failed.")
            # setup list
            params_array = np.array(
                [
                    [
                        pres.R_init,
                        pres.R_final,
                        pres.R_xy_init,
                        pres.R_xy_final,
                        pres.G,
                        pres.B,
                        pres.rotx * 180 / math.pi,
                        pres.roty * 180 / math.pi,
                        pres.ry,
                        pres.rz,
                        pres.r0,
                        pres.re,
                        pres.voigt_nu,
                        pres.uc_params[0],
                        pres.uc_params[1],
                        pres.uc_params[2],
                        pres.uc_params[3],
                        pres.uc_params[4],
                        pres.uc_params[5],
                        pres.CC_final,
                        pres.pickle_filename,
                    ]
                    for pres in results
                ]
            )
            params = [
                "Rinit",
                "Rfinal",
                "Rxyinit",
                "Rxyfinal",
                "G",
                "B",
                "rot_x",
                "rot_y",
                "gamma_y",
                "gamma_z",
                "gamma_0",
                "gamma_e",
                "voigtnu",
                "a",
                "b",
                "c",
                "alpha",
                "beta",
                "gamma",
                "CC",
                "Filename",
            ]
        # keep parameter history if verbose is selected
        if iparams.flag_output_verbose:
            fileseq_list = flex.int()
            for file_in in os.listdir(iparams.run_no):
                if file_in.endswith(".paramhist"):
                    file_split = file_in.split(".")
                    fileseq_list.append(int(file_split[0]))
            if len(fileseq_list) == 0:
                new_fileseq = 0
            else:
                new_fileseq = flex.max(fileseq_list) + 1
            newfile_name = str(new_fileseq) + ".paramhist"
            txt_out_verbose = "\n".join([" ".join(p) for p in params_array])
            f = open(iparams.run_no + "/" + newfile_name, "w")
            f.write(txt_out_verbose)
            f.close()
        # plotting
        if iparams.flag_plot:
            try:
                import matplotlib.pyplot as plt
            except Exception as e:
                print("Warning: error importing matplotlib.pyplot")
                print(e)
                return
            n_rows = 3
            n_cols = int(math.ceil(len(params) / n_rows))
            num_bins = 10
            for i in range(len(params) - 1):
                tmp_params = params_array[:, i].astype(np.float)
                plt.subplot(n_rows, n_cols, i + 1)
                plt.hist(tmp_params, num_bins, normed=0, facecolor="green", alpha=0.5)
                plt.ylabel("Frequencies")
                plt.title(
                    params[i]
                    + "\nmu %5.1f med %5.1f sigma %5.1f"
                    % (np.mean(tmp_params), np.median(tmp_params), np.std(tmp_params))
                )
            plt.show()

    def combine_pre_merge(self, result, iparams):
        mi_all = flex.miller_index()
        mio_all = flex.miller_index()
        I_all = flex.double()
        sigI_all = flex.double()
        G_all = flex.double()
        B_all = flex.double()
        p_all = flex.double()
        rs_all = flex.double()
        wavelength_all = flex.double()
        sin_all = flex.double()
        SE_all = flex.double()
        uc_mean_set = []
        wavelength_mean_set = []
        pickle_filename_all = flex.std_string()
        for res in result:
            for prep_output in res:
                _, _, mi, mio, I, sigI, G, B, p, rs, wavelength, sin, SE, uc_mean, wavelength_mean, pickle_filename_set, txt_out = (
                    prep_output
                )
                mi_all.extend(mi)
                mio_all.extend(mio)
                I_all.extend(I)
                sigI_all.extend(sigI)
                G_all.extend(G)
                B_all.extend(B)
                p_all.extend(p)
                rs_all.extend(rs)
                wavelength_all.extend(wavelength)
                sin_all.extend(sin)
                SE_all.extend(SE)
                uc_mean_set.extend(uc_mean)
                wavelength_mean_set.append(wavelength_mean)
                pickle_filename_all.extend(pickle_filename_set)
        uc_mean = np.mean(np.array(uc_mean_set).reshape(-1, 6), axis=0)
        wavelength_mean = np.mean(wavelength_mean_set)
        ms_template = crystal.symmetry(
            unit_cell=tuple(uc_mean), space_group_symbol=str(iparams.target_space_group)
        ).build_miller_set(
            anomalous_flag=iparams.target_anomalous_flag, d_min=iparams.merge.d_min
        )
        ma_all = ms_template.array().customized_copy(
            indices=mi_all, data=I_all, sigmas=sigI_all
        )
        # sort reflections according to asymmetric-unit symmetry hkl
        perm = ma_all.sort_permutation(by_value="packed_indices")
        mi_all_sort = mi_all.select(perm)
        mio_all_sort = mio_all.select(perm)
        I_all_sort = I_all.select(perm)
        sigI_all_sort = sigI_all.select(perm)
        G_all_sort = G_all.select(perm)
        B_all_sort = B_all.select(perm)
        p_all_sort = p_all.select(perm)
        rs_all_sort = rs_all.select(perm)
        wavelength_all_sort = wavelength_all.select(perm)
        sin_all_sort = sin_all.select(perm)
        SE_all_sort = SE_all.select(perm)
        pickle_filename_all_sort = pickle_filename_all.select(perm)
        ma_uniq = (
            ma_all.merge_equivalents()
            .array()
            .complete_array(d_min=iparams.merge.d_min, d_max=iparams.merge.d_max)
        )
        matches_uniq = miller.match_multi_indices(
            miller_indices_unique=ma_uniq.indices(), miller_indices=mi_all_sort
        )
        pair_0 = flex.int([pair[0] for pair in matches_uniq.pairs()])
        pair_1 = flex.int([pair[1] for pair in matches_uniq.pairs()])
        group_id_list = flex.int(
            [pair_0[pair_1[i]] for i in range(len(matches_uniq.pairs()))]
        )
        tally = Counter()
        for elem in group_id_list:
            tally[elem] += 1
        cn_group = len(tally)
        return (
            cn_group,
            group_id_list,
            mi_all_sort,
            mio_all_sort,
            I_all_sort,
            sigI_all_sort,
            G_all_sort,
            B_all_sort,
            p_all_sort,
            rs_all_sort,
            wavelength_all_sort,
            sin_all_sort,
            SE_all_sort,
            uc_mean,
            wavelength_mean,
            pickle_filename_all_sort,
            "",
        )
