"""
Author      : Uervirojnangkoorn, M.
Created     : 7/13/2014
Description :

The leastsqr_handler class refines post-refinement parameters for each frame m in:

J = sum(I_hi - (G0 * exp(-2B*sin_theta_over_lambda_sq) * p * I_h))^2

where I_hi is observed intensitiy, G0 and B are scale factors, and p
is the partiality function of these parameters:
rot_x, rot_y: crystal orientation angles around x- and y- axes
gamma_y, gamma_z: parameters to model anisotropy of crystal mosaicity
gamma_e: spectral dispersion
unit-cell parameters: a,b,c,alpha,beta,gamma.

The class implements Leveberg-Marquardt algorithms to scale the refined parameters
using the lamda updates. The unit-cell parameters are refined with restraints based
on the 7 crystal systems (6 conditions).
"""
from __future__ import absolute_import, division, print_function
import math
from cctbx.array_family import flex
from scitbx.matrix import sqr
from cctbx.uctbx import unit_cell
from cctbx.crystal_orientation import crystal_orientation
from .mod_lbfgs import lbfgs_handler
from .mod_lbfgs_partiality import lbfgs_partiality_handler
from .mod_partiality import partiality_handler
from six.moves import range


def coefficient_of_determination(y, y_model):
    mean_y = flex.mean(y)
    r_sqr = flex.sum((y_model - mean_y) ** 2) / flex.sum((y - mean_y) ** 2)
    return r_sqr


def standard_error_of_the_estimate(y, y_model, n_params):
    s = math.sqrt(flex.sum((y - y_model) ** 2) / (len(y) - n_params))
    return s


def good_unit_cell(uc_params, iparams, uc_tol, target_unit_cell=None):
    if iparams is None:
        expected_unit_cell = target_unit_cell
    else:
        expected_unit_cell = iparams.target_unit_cell
    flag_good_uc = False
    if (
        abs(uc_params[0] - expected_unit_cell.parameters()[0])
        <= (uc_tol * expected_unit_cell.parameters()[0] / 100)
        and abs(uc_params[1] - expected_unit_cell.parameters()[1])
        <= (uc_tol * expected_unit_cell.parameters()[1] / 100)
        and abs(uc_params[2] - expected_unit_cell.parameters()[2])
        <= (uc_tol * expected_unit_cell.parameters()[2] / 100)
        and abs(uc_params[3] - expected_unit_cell.parameters()[3])
        <= (uc_tol * expected_unit_cell.parameters()[3] / 100)
        and abs(uc_params[4] - expected_unit_cell.parameters()[4])
        <= (uc_tol * expected_unit_cell.parameters()[4] / 100)
        and abs(uc_params[5] - expected_unit_cell.parameters()[5])
        <= (uc_tol * expected_unit_cell.parameters()[5] / 100)
    ):
        flag_good_uc = True
    return flag_good_uc


class leastsqr_handler(object):
    """A wrapper class for least-squares refinement."""

    def __init__(self):
        """Intialitze parameters."""

    def get_filtered_data(
        self,
        filter_mode,
        filter_params,
        observations_in,
        alpha_angle_in,
        spot_pred_x_mm_in,
        spot_pred_y_mm_in,
        I_ref_in,
        partiality_in=False,
    ):
        if filter_mode == "resolution":
            i_sel = observations_in.resolution_filter_selection(
                d_min=filter_params[0], d_max=filter_params[1]
            )
        elif filter_mode == "sigma":
            i_sel = (observations_in.data() / observations_in.sigmas()) > filter_params[
                0
            ]
        elif filter_mode == "partiality":
            i_sel = partiality_in > filter_params[0]
        observations_out = observations_in.select(i_sel)
        alpha_angle_out = alpha_angle_in.select(i_sel)
        spot_pred_x_mm_out = spot_pred_x_mm_in.select(i_sel)
        spot_pred_y_mm_out = spot_pred_y_mm_in.select(i_sel)
        I_ref_out = I_ref_in.select(i_sel)
        return (
            observations_out,
            alpha_angle_out,
            spot_pred_x_mm_out,
            spot_pred_y_mm_out,
            I_ref_out,
        )

    def optimize_scalefactors(
        self,
        I_r_flex,
        observations_original,
        wavelength,
        crystal_init_orientation,
        alpha_angle,
        spot_pred_x_mm,
        spot_pred_y_mm,
        iparams,
        pres_in,
        observations_non_polar,
        detector_distance_mm,
        const_params,
    ):
        ph = partiality_handler()
        pr_d_min = iparams.postref.scale.d_min
        pr_d_max = iparams.postref.scale.d_max
        pr_sigma_min = iparams.postref.scale.sigma_min
        pr_partiality_min = iparams.postref.scale.partiality_min
        # filter by resolution
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "resolution",
            [pr_d_min, pr_d_max],
            observations_original,
            alpha_angle,
            spot_pred_x_mm,
            spot_pred_y_mm,
            I_r_flex,
        )
        # filter by sigma
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "sigma",
            [pr_sigma_min],
            observations_original_sel,
            alpha_angle_sel,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            I_ref_sel,
        )
        I_r_true = I_ref_sel[:]
        I_o_true = observations_original_sel.data()[:]
        if pres_in is not None:
            G, B, b0 = pres_in.G, pres_in.B, pres_in.B
        else:
            G, B, b0 = (1, 0, 0)
        refine_mode = "scale_factor"
        xinp = flex.double([G, B])
        args = (
            I_r_true,
            observations_original_sel,
            wavelength,
            alpha_angle_sel,
            crystal_init_orientation,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            detector_distance_mm,
            refine_mode,
            const_params,
            b0,
            None,
            iparams,
        )
        lh = lbfgs_handler(current_x=xinp, args=args)
        G_fin, B_fin = (lh.x[0], lh.x[1])
        rotx, roty, ry, rz, r0, re, voigt_nu, a, b, c, alpha, beta, gamma = const_params
        two_theta = observations_original.two_theta(wavelength=wavelength)
        sin_theta_over_lambda_sq = two_theta.sin_theta_over_lambda_sq().data()
        uc = unit_cell((a, b, c, alpha, beta, gamma))
        ph = partiality_handler()
        partiality_init, delta_xy_init, rs_init, dummy = ph.calc_partiality_anisotropy_set(
            uc,
            rotx,
            roty,
            observations_original.indices(),
            ry,
            rz,
            r0,
            re,
            voigt_nu,
            two_theta.data(),
            alpha_angle,
            wavelength,
            crystal_init_orientation,
            spot_pred_x_mm,
            spot_pred_y_mm,
            detector_distance_mm,
            iparams.partiality_model,
            iparams.flag_beam_divergence,
        )
        I_o_init = ph.calc_full_refl(
            observations_original.data(),
            sin_theta_over_lambda_sq,
            G,
            B,
            partiality_init,
            rs_init,
        )
        I_o_fin = ph.calc_full_refl(
            observations_original.data(),
            sin_theta_over_lambda_sq,
            G_fin,
            B_fin,
            partiality_init,
            rs_init,
        )
        SE_of_the_estimate = standard_error_of_the_estimate(I_r_flex, I_o_fin, 2)
        R_sq = coefficient_of_determination(I_r_flex, I_o_fin) * 100
        CC_init = flex.linear_correlation(I_r_flex, I_o_init).coefficient()
        CC_final = flex.linear_correlation(I_r_flex, I_o_fin).coefficient()
        err_init = (I_r_flex - I_o_init) / observations_original.sigmas()
        R_init = math.sqrt(flex.sum(err_init ** 2))
        err_final = (I_r_flex - I_o_fin) / observations_original.sigmas()
        R_final = math.sqrt(flex.sum(err_final ** 2))
        R_xy_init = 0
        R_xy_final = 0
        CC_iso_init = 0
        CC_iso_final = 0
        return (
            flex.double(list(lh.x)),
            (
                SE_of_the_estimate,
                R_sq,
                CC_init,
                CC_final,
                R_init,
                R_final,
                R_xy_init,
                R_xy_final,
                CC_iso_init,
                CC_iso_final,
            ),
        )

    def prepare_data_microcycle(
        self,
        refine_mode,
        iparams,
        observations_original,
        alpha_angle,
        spot_pred_x_mm,
        spot_pred_y_mm,
        I_r_flex,
        init_params,
        crystal_init_orientation,
        wavelength,
        detector_distance_mm,
    ):
        # prepare data
        if refine_mode == "crystal_orientation":
            pr_d_min = iparams.postref.crystal_orientation.d_min
            pr_d_max = iparams.postref.crystal_orientation.d_max
            pr_sigma_min = iparams.postref.crystal_orientation.sigma_min
            pr_partiality_min = iparams.postref.crystal_orientation.partiality_min
            pr_uc_tol = iparams.postref.unit_cell.uc_tolerance
        elif refine_mode == "reflecting_range":
            pr_d_min = iparams.postref.reflecting_range.d_min
            pr_d_max = iparams.postref.reflecting_range.d_max
            pr_sigma_min = iparams.postref.reflecting_range.sigma_min
            pr_partiality_min = iparams.postref.reflecting_range.partiality_min
            pr_uc_tol = iparams.postref.unit_cell.uc_tolerance
        elif refine_mode == "unit_cell":
            pr_d_min = iparams.postref.unit_cell.d_min
            pr_d_max = iparams.postref.unit_cell.d_max
            pr_sigma_min = iparams.postref.unit_cell.sigma_min
            pr_partiality_min = iparams.postref.unit_cell.partiality_min
            pr_uc_tol = iparams.postref.unit_cell.uc_tolerance
        elif refine_mode == "allparams":
            pr_d_min = iparams.postref.allparams.d_min
            pr_d_max = iparams.postref.allparams.d_max
            pr_sigma_min = iparams.postref.allparams.sigma_min
            pr_partiality_min = iparams.postref.allparams.partiality_min
            pr_uc_tol = iparams.postref.unit_cell.uc_tolerance
        # filter by resolution
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "resolution",
            [pr_d_min, pr_d_max],
            observations_original,
            alpha_angle,
            spot_pred_x_mm,
            spot_pred_y_mm,
            I_r_flex,
        )
        # filter by sigma
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "sigma",
            [pr_sigma_min],
            observations_original_sel,
            alpha_angle_sel,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            I_ref_sel,
        )
        # extract refined parameters
        G, B, rotx, roty, ry, rz, r0, re, voigt_nu, a, b, c, alpha, beta, gamma = (
            init_params
        )
        # filter by partiality
        two_theta = observations_original_sel.two_theta(wavelength=wavelength).data()
        uc = unit_cell((a, b, c, alpha, beta, gamma))
        ph = partiality_handler()
        partiality_init, delta_xy_init, rs_init, dummy = ph.calc_partiality_anisotropy_set(
            uc,
            rotx,
            roty,
            observations_original_sel.indices(),
            ry,
            rz,
            r0,
            re,
            voigt_nu,
            two_theta,
            alpha_angle_sel,
            wavelength,
            crystal_init_orientation,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            detector_distance_mm,
            iparams.partiality_model,
            iparams.flag_beam_divergence,
        )
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "partiality",
            [pr_partiality_min],
            observations_original_sel,
            alpha_angle_sel,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            I_ref_sel,
            partiality_in=partiality_init,
        )
        return (
            observations_original_sel,
            alpha_angle_sel,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            I_ref_sel,
        )

    def optimize(
        self,
        I_r_flex,
        observations_original,
        wavelength,
        crystal_init_orientation,
        alpha_angle,
        spot_pred_x_mm,
        spot_pred_y_mm,
        iparams,
        pres_in,
        observations_non_polar,
        detector_distance_mm,
    ):
        ph = partiality_handler()
        lph = lbfgs_partiality_handler()
        if iparams.postref.allparams.flag_on:
            refine_steps = ["allparams"]
        else:
            refine_steps = ["crystal_orientation"]
            if iparams.postref.reflecting_range.flag_on:
                refine_steps.append("reflecting_range")
            if iparams.postref.unit_cell.flag_on:
                refine_steps.append("unit_cell")
        # get miller array iso, if given.
        miller_array_iso = None
        # prepare data
        pr_d_min = iparams.postref.allparams.d_min
        pr_d_max = iparams.postref.allparams.d_max
        pr_sigma_min = iparams.postref.allparams.sigma_min
        pr_partiality_min = iparams.postref.allparams.partiality_min
        pr_uc_tol = iparams.postref.allparams.uc_tolerance
        cs = observations_original.crystal_symmetry().space_group().crystal_system()
        # filter by resolution
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "resolution",
            [pr_d_min, pr_d_max],
            observations_original,
            alpha_angle,
            spot_pred_x_mm,
            spot_pred_y_mm,
            I_r_flex,
        )
        # filter by sigma
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "sigma",
            [pr_sigma_min],
            observations_original_sel,
            alpha_angle_sel,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            I_ref_sel,
        )
        # initialize values only in the first sub cycle and the first refine step.
        spot_radius = ph.calc_spot_radius(
            sqr(crystal_init_orientation.reciprocal_matrix()),
            observations_original_sel.indices(),
            wavelength,
        )
        if pres_in is None:
            ry, rz, r0, re, voigt_nu, rotx, roty = (
                0,
                0,
                spot_radius,
                iparams.gamma_e,
                iparams.voigt_nu,
                0.0,
                0.0,
            )
            # apply constrain on the unit cell using crystal system
            uc_scale_inp = lph.prep_input(
                observations_original.unit_cell().parameters(), cs
            )
            uc_scale_constrained = lph.prep_output(uc_scale_inp, cs)
            a, b, c, alpha, beta, gamma = uc_scale_constrained
            const_params_scale = (
                rotx,
                roty,
                ry,
                rz,
                r0,
                re,
                voigt_nu,
                a,
                b,
                c,
                alpha,
                beta,
                gamma,
            )
            xopt_scalefactors, stats = self.optimize_scalefactors(
                I_r_flex,
                observations_original,
                wavelength,
                crystal_init_orientation,
                alpha_angle,
                spot_pred_x_mm,
                spot_pred_y_mm,
                iparams,
                pres_in,
                observations_non_polar,
                detector_distance_mm,
                const_params_scale,
            )
            G, B = xopt_scalefactors
        else:
            G, B, ry, rz, r0, re, voigt_nu, rotx, roty = (
                pres_in.G,
                pres_in.B,
                pres_in.ry,
                pres_in.rz,
                pres_in.r0,
                pres_in.re,
                pres_in.voigt_nu,
                0.0,
                0.0,
            )
            a, b, c, alpha, beta, gamma = pres_in.unit_cell.parameters()
            crystal_init_orientation = pres_in.crystal_orientation
        # filter by partiality
        two_theta = observations_original_sel.two_theta(wavelength=wavelength).data()
        uc = unit_cell((a, b, c, alpha, beta, gamma))
        partiality_init, delta_xy_init, rs_init, dummy = ph.calc_partiality_anisotropy_set(
            uc,
            rotx,
            roty,
            observations_original_sel.indices(),
            ry,
            rz,
            r0,
            re,
            voigt_nu,
            two_theta,
            alpha_angle_sel,
            wavelength,
            crystal_init_orientation,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            detector_distance_mm,
            iparams.partiality_model,
            iparams.flag_beam_divergence,
        )
        observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.get_filtered_data(
            "partiality",
            [pr_partiality_min],
            observations_original_sel,
            alpha_angle_sel,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            I_ref_sel,
            partiality_in=partiality_init,
        )
        I_r_true = I_ref_sel[:]
        I_o_true = observations_original_sel.data()[:]
        # calculate initial residual_xy error
        const_params_uc = (G, B, rotx, roty, ry, rz, r0, re, voigt_nu)
        xinp_uc = lph.prep_input((a, b, c, alpha, beta, gamma), cs)
        args_uc = (
            I_r_true,
            observations_original_sel,
            wavelength,
            alpha_angle_sel,
            crystal_init_orientation,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            detector_distance_mm,
            "unit_cell",
            const_params_uc,
            B,
            miller_array_iso,
            iparams,
        )
        uc_params_err = lph.func(xinp_uc, args_uc)
        init_residual_xy_err = flex.sum(uc_params_err ** 2)
        # calculate initial residual_pr error
        const_params_all = (G, B)
        xinp_all = flex.double([rotx, roty, ry, rz, r0, re, voigt_nu])
        xinp_all.extend(lph.prep_input((a, b, c, alpha, beta, gamma), cs))
        args_all = (
            I_r_true,
            observations_original_sel,
            wavelength,
            alpha_angle_sel,
            crystal_init_orientation,
            spot_pred_x_mm_sel,
            spot_pred_y_mm_sel,
            detector_distance_mm,
            "allparams",
            const_params_all,
            B,
            miller_array_iso,
            iparams,
        )
        all_params_err = lph.func(xinp_all, args_all)
        init_residual_err = flex.sum(all_params_err ** 2)
        # keep in list
        t_pr_list = [init_residual_err]
        t_xy_list = [init_residual_xy_err]
        refined_params_hist = [
            (G, B, rotx, roty, ry, rz, r0, re, voigt_nu, a, b, c, alpha, beta, gamma)
        ]
        txt_out = ""
        for i_sub_cycle in range(iparams.n_postref_sub_cycle):
            for j_refine_step in range(len(refine_steps)):
                refine_mode = refine_steps[j_refine_step]
                # prepare data
                init_params = (
                    G,
                    B,
                    rotx,
                    roty,
                    ry,
                    rz,
                    r0,
                    re,
                    voigt_nu,
                    a,
                    b,
                    c,
                    alpha,
                    beta,
                    gamma,
                )
                observations_original_sel, alpha_angle_sel, spot_pred_x_mm_sel, spot_pred_y_mm_sel, I_ref_sel = self.prepare_data_microcycle(
                    refine_mode,
                    iparams,
                    observations_original,
                    alpha_angle,
                    spot_pred_x_mm,
                    spot_pred_y_mm,
                    I_r_flex,
                    init_params,
                    crystal_init_orientation,
                    wavelength,
                    detector_distance_mm,
                )
                I_r_true = I_ref_sel[:]
                I_o_true = observations_original_sel.data()
                if refine_mode == "crystal_orientation":
                    xinp = flex.double([rotx, roty])
                    const_params = (
                        G,
                        B,
                        ry,
                        rz,
                        r0,
                        re,
                        voigt_nu,
                        a,
                        b,
                        c,
                        alpha,
                        beta,
                        gamma,
                    )
                elif refine_mode == "reflecting_range":
                    xinp = flex.double([ry, rz, r0, re, voigt_nu])
                    const_params = (G, B, rotx, roty, a, b, c, alpha, beta, gamma)
                elif refine_mode == "unit_cell":
                    xinp = lph.prep_input((a, b, c, alpha, beta, gamma), cs)
                    const_params = (G, B, rotx, roty, ry, rz, r0, re, voigt_nu)
                elif refine_mode == "allparams":
                    xinp = flex.double([rotx, roty, ry, rz, r0, re, voigt_nu])
                    xinp.extend(lph.prep_input((a, b, c, alpha, beta, gamma), cs))
                    const_params = (G, B)
                args = (
                    I_r_true,
                    observations_original_sel,
                    wavelength,
                    alpha_angle_sel,
                    crystal_init_orientation,
                    spot_pred_x_mm_sel,
                    spot_pred_y_mm_sel,
                    detector_distance_mm,
                    refine_mode,
                    const_params,
                    B,
                    miller_array_iso,
                    iparams,
                )
                lh = lbfgs_handler(current_x=xinp, args=args)
                xopt = flex.double(list(lh.x))
                if (
                    refine_mode == "crystal_orientation"
                    or refine_mode == "reflecting_range"
                    or refine_mode == "allparams"
                ):
                    current_residual_err = lh.f
                    # calculate residual_xy_error (for refine_mode = SF, CO, RR, and all params)
                    xinp_uc = lph.prep_input((a, b, c, alpha, beta, gamma), cs)
                    if refine_mode == "crystal_orientation":
                        rotx, roty = xopt
                    elif refine_mode == "reflecting_range":
                        ry, rz, r0, re, voigt_nu = xopt
                    elif refine_mode == "allparams":
                        rotx, roty, ry, rz, r0, re, voigt_nu = xopt[:7]
                        xinp_uc = xopt[7:]
                        a, b, c, alpha, beta, gamma = lph.prep_output(xinp_uc, cs)
                    const_params_uc = (G, B, rotx, roty, ry, rz, r0, re, voigt_nu)
                    xinp_uc = lph.prep_input((a, b, c, alpha, beta, gamma), cs)
                    args_uc = (
                        I_r_true,
                        observations_original_sel,
                        wavelength,
                        alpha_angle_sel,
                        crystal_init_orientation,
                        spot_pred_x_mm_sel,
                        spot_pred_y_mm_sel,
                        detector_distance_mm,
                        "unit_cell",
                        const_params_uc,
                        B,
                        miller_array_iso,
                        iparams,
                    )
                    uc_params_err = lph.func(xinp_uc, args_uc)
                    current_residual_xy_err = flex.sum(uc_params_err ** 2)
                elif refine_mode == "unit_cell":
                    current_residual_xy_err = lh.f
                    xopt_uc = lph.prep_output(xopt, cs)
                    a, b, c, alpha, beta, gamma = xopt_uc
                    # check the unit-cell with the reference intensity
                    xinp = flex.double([rotx, roty, ry, rz, r0, re, voigt_nu])
                    xinp.extend(lph.prep_input((a, b, c, alpha, beta, gamma), cs))
                    const_params_all = (G, B)
                    args_all = (
                        I_r_true,
                        observations_original_sel,
                        wavelength,
                        alpha_angle_sel,
                        crystal_init_orientation,
                        spot_pred_x_mm_sel,
                        spot_pred_y_mm_sel,
                        detector_distance_mm,
                        "allparams",
                        const_params_all,
                        B,
                        miller_array_iso,
                        iparams,
                    )
                    all_params_err = lph.func(xinp_all, args_all)
                    current_residual_err = flex.sum(all_params_err ** 2)
                flag_success = False
                if refine_mode == "allparams":
                    # if allparams refinement, only check the post-refine target function
                    if current_residual_err < (
                        t_pr_list[len(t_pr_list) - 1]
                        + (
                            t_pr_list[len(t_pr_list) - 1]
                            * iparams.postref.residual_threshold
                            / 100
                        )
                    ):
                        t_pr_list.append(current_residual_err)
                        t_xy_list.append(current_residual_xy_err)
                        refined_params_hist.append(
                            (
                                G,
                                B,
                                rotx,
                                roty,
                                ry,
                                rz,
                                r0,
                                re,
                                voigt_nu,
                                a,
                                b,
                                c,
                                alpha,
                                beta,
                                gamma,
                            )
                        )
                        flag_success = True
                else:
                    if current_residual_err < (
                        t_pr_list[len(t_pr_list) - 1]
                        + (
                            t_pr_list[len(t_pr_list) - 1]
                            * iparams.postref.residual_threshold
                            / 100
                        )
                    ):
                        if current_residual_xy_err < (
                            t_xy_list[len(t_xy_list) - 1]
                            + (
                                t_xy_list[len(t_xy_list) - 1]
                                * iparams.postref.residual_threshold_xy
                                / 100
                            )
                        ):
                            t_pr_list.append(current_residual_err)
                            t_xy_list.append(current_residual_xy_err)
                            refined_params_hist.append(
                                (
                                    G,
                                    B,
                                    rotx,
                                    roty,
                                    ry,
                                    rz,
                                    r0,
                                    re,
                                    voigt_nu,
                                    a,
                                    b,
                                    c,
                                    alpha,
                                    beta,
                                    gamma,
                                )
                            )
                            flag_success = True
                if flag_success is False:
                    G, B, rotx, roty, ry, rz, r0, re, voigt_nu, a, b, c, alpha, beta, gamma = refined_params_hist[
                        len(refined_params_hist) - 1
                    ]
                tmp_txt_out = (
                    refine_mode
                    + " %3.0f %6.4f %6.4f %6.4f %6.4f %10.8f %10.8f %10.8f %10.8f %10.8f %6.3f %6.3f %.4g %6.3f\n"
                    % (
                        i_sub_cycle,
                        G,
                        B,
                        rotx * 180 / math.pi,
                        roty * 180 / math.pi,
                        ry,
                        rz,
                        r0,
                        re,
                        voigt_nu,
                        a,
                        c,
                        t_pr_list[len(t_pr_list) - 1],
                        t_xy_list[len(t_pr_list) - 1],
                    )
                )
                txt_out += tmp_txt_out
        # apply the refined parameters on the full (original) reflection set
        two_theta = observations_original.two_theta(wavelength=wavelength).data()
        sin_theta_over_lambda_sq = (
            observations_original.two_theta(wavelength=wavelength)
            .sin_theta_over_lambda_sq()
            .data()
        )
        if pres_in is None:
            partiality_init, delta_xy_init, rs_init, rh_init = ph.calc_partiality_anisotropy_set(
                observations_original.unit_cell(),
                0.0,
                0.0,
                observations_original.indices(),
                0,
                0,
                spot_radius,
                iparams.gamma_e,
                iparams.voigt_nu,
                two_theta,
                alpha_angle,
                wavelength,
                crystal_init_orientation,
                spot_pred_x_mm,
                spot_pred_y_mm,
                detector_distance_mm,
                iparams.partiality_model,
                iparams.flag_beam_divergence,
            )
            I_o_init = ph.calc_full_refl(
                observations_original.data(),
                sin_theta_over_lambda_sq,
                1,
                0,
                partiality_init,
                rs_init,
            )
        else:
            partiality_init, delta_xy_init, rs_init, rh_init = ph.calc_partiality_anisotropy_set(
                pres_in.unit_cell,
                0.0,
                0.0,
                observations_original.indices(),
                pres_in.ry,
                pres_in.rz,
                pres_in.r0,
                pres_in.re,
                pres_in.voigt_nu,
                two_theta,
                alpha_angle,
                wavelength,
                crystal_init_orientation,
                spot_pred_x_mm,
                spot_pred_y_mm,
                detector_distance_mm,
                iparams.partiality_model,
                iparams.flag_beam_divergence,
            )
            I_o_init = ph.calc_full_refl(
                observations_original.data(),
                sin_theta_over_lambda_sq,
                pres_in.G,
                pres_in.B,
                partiality_init,
                rs_init,
            )
        partiality_fin, delta_xy_fin, rs_fin, rh_fin = ph.calc_partiality_anisotropy_set(
            unit_cell((a, b, c, alpha, beta, gamma)),
            rotx,
            roty,
            observations_original.indices(),
            ry,
            rz,
            r0,
            re,
            voigt_nu,
            two_theta,
            alpha_angle,
            wavelength,
            crystal_init_orientation,
            spot_pred_x_mm,
            spot_pred_y_mm,
            detector_distance_mm,
            iparams.partiality_model,
            iparams.flag_beam_divergence,
        )
        I_o_fin = ph.calc_full_refl(
            observations_original.data(),
            sin_theta_over_lambda_sq,
            G,
            B,
            partiality_fin,
            rs_fin,
        )
        SE_of_the_estimate = standard_error_of_the_estimate(I_r_flex, I_o_fin, 13)
        R_sq = coefficient_of_determination(I_r_flex, I_o_fin) * 100
        CC_init = flex.linear_correlation(I_r_flex, I_o_init).coefficient()
        CC_final = flex.linear_correlation(I_r_flex, I_o_fin).coefficient()
        err_init = (I_r_flex - I_o_init) / observations_original.sigmas()
        R_init = math.sqrt(flex.sum(err_init ** 2))
        err_final = (I_r_flex - I_o_fin) / observations_original.sigmas()
        R_final = math.sqrt(flex.sum(err_final ** 2))
        R_xy_init = math.sqrt(flex.sum(delta_xy_init ** 2))
        R_xy_final = math.sqrt(flex.sum(delta_xy_fin ** 2))
        if R_init < R_final or re > (iparams.gamma_e * 3):
            CC_final = CC_init
            R_final = R_init
            R_xy_final = R_xy_init
            if pres_in is None:
                G, B, r0, ry, rz, re, rotx, roty = (
                    1.0,
                    0.0,
                    spot_radius,
                    0.0,
                    0.0,
                    iparams.gamma_e,
                    0.0,
                    0.0,
                )
                a, b, c, alpha, beta, gamma = (
                    observations_original.unit_cell().parameters()
                )
            else:
                G, B, r0, ry, rz, re, rotx, roty = (
                    pres_in.G,
                    pres_in.B,
                    pres_in.r0,
                    pres_in.ry,
                    pres_in.rz,
                    pres_in.re,
                    0.0,
                    0.0,
                )
                a, b, c, alpha, beta, gamma = pres_in.unit_cell.parameters()
                crystal_init_orientation = pres_in.crystal_orientation
        # calculate CCiso if hklisoin is given
        CC_iso_init, CC_iso_final = (0, 0)
        if iparams.hklisoin is not None:
            if miller_array_iso is not None:
                from cctbx import miller

                matches = miller.match_multi_indices(
                    miller_indices_unique=miller_array_iso.indices(),
                    miller_indices=observations_non_polar.indices(),
                )
                I_iso_match = flex.double(
                    [miller_array_iso.data()[pair[0]] for pair in matches.pairs()]
                )
                I_o_init_match = flex.double(
                    [I_o_init[pair[1]] for pair in matches.pairs()]
                )
                I_o_fin_match = flex.double(
                    [I_o_fin[pair[1]] for pair in matches.pairs()]
                )
                CC_iso_init = flex.linear_correlation(
                    I_iso_match, I_o_init_match
                ).coefficient()
                CC_iso_final = flex.linear_correlation(
                    I_iso_match, I_o_fin_match
                ).coefficient()
        xopt = (G, B, rotx, roty, ry, rz, r0, re, voigt_nu, a, b, c, alpha, beta, gamma)
        return (
            xopt,
            (
                SE_of_the_estimate,
                R_sq,
                CC_init,
                CC_final,
                R_init,
                R_final,
                R_xy_init,
                R_xy_final,
                CC_iso_init,
                CC_iso_final,
            ),
            len(I_ref_sel),
        )
