Bootstrap approach for error determination

[jadball]

Hi Jon,

I remember us discussing the bootstrap approach for strain and position error determination some time ago - I seem to recall you mentioning that there was an old attempt at this somewhere in the ImageD11 codebase, but I couldn't find it here or in your fork. Can you remember where it was? If not, no worries, I can try to get something started myself and i'll contribute it back upstream.

Best wishes,

James

[jonwright]

Hi James - I can't find it right now but maybe there is something in the fork of @younes-elhachi ? I seem to remember he had a version in the past too. Pull requests are welcome - I hope I have not lost one! Best Jon

[jadball]

@jonwright I had a quick look but couldn't find anything. I'm working on the following implementation:

1. Do a grid index as normal with grid_index_parallel
2. Run makemap on the output map against all the peaks, with a tolerance sequence (say, down to 0.01)
3. Take all the peaks, and generate 10 different peaks files, each with a different 10% of peaks removed at random (ensuring the same peaks aren't removed from multiple columnfiles)
4. For each of the output peaks files, run makemap with 0.01 tolerance, thereby generating 10 final grain maps
5. Average out the 10 final grain maps, and use them to get a standard deviation

Unfortunately we run into similar problems as issue #153 - how do we sensibly average the grains together and get the errors? Translation is easy, but deformed unit cells are harder. My initial idea was to get the eps_grain for each grain, average them to get an average deformation. For the orientation, I am averaging the Rodrigues vectors and then converting back to a mean U matrix (to ensure the matrix remains symmetric). Then we can get the B matrix from the mean eps_grain, and the mean U matrix, and get a mean UBI. However I am running into problems:

b_array = np.array([a_grain.B for a_grain in grain_list])
mean_b = np.mean(b_array, axis=0)

reference_unit_cell = [3.6223, 3.6223, 3.6223, 90, 90, 90]

eps_array = np.array([a_grain.eps_grain(dzero_cell=reference_unit_cell) for a_grain in grain_list])
eps_mean = np.mean(eps_array, axis=0)
eps_stdev = np.std(eps_array, axis=0)

mean_b_from_epsilon = xfab.tools.epsilon_to_b(eps_mean, reference_unit_cell)

print(mean_b)
print(mean_b_from_epsilon)

mean_b gives me:

[[ 2.76068509e-01 -8.51720013e-06 -5.47933649e-07]
 [ 0.00000000e+00  2.76074043e-01  6.03700136e-06]
 [ 0.00000000e+00  0.00000000e+00  2.76070064e-01]]

but mean_b_from_epsilon gives me:

[[ 1.73458959e+00 -5.35039591e-05 -3.43484184e-06]
 [ 0.00000000e+00  1.73462436e+00  3.79229873e-05]
 [ 0.00000000e+00  0.00000000e+00  1.73459937e+00]]

I think mean_b doesn't make sense geometrically as you're averaging out the lengths in reciprocal space and you're not maintaining consistency between the side lengths and the angles, hence my idea to use the mean eps. However, the mean_b_from_epsilon seems very wrong. Do you have any pointers?

Thanks in advance!

[jonwright]

Is it just the two pi that xfab uses when returning the B matrix?

[jadball]

...indeed it is! Thank you. I'm going to leave this issue open until I have a pull request for you - just testing at the moment.

[jadball]

Hi @jonwright - I'm making progress but have run into a bit of a roadblock. My current idea is just to take the U of each grain, rotate it into the fundamental zone, then make a new UBI from it - that way hopefully the UBIs can be directly averaged to get a new UBI. My problem comes in determining the strain error. I'm using your deformation gradient tensor method to get the strain from the UBI, but I am unsure how I would go about propagating the UBI errors through to get an error in eps. The SVD in particular is throwing me off - do you have any ideas? My backup plan is to get the strains from each grain's UBI and find the standard deviation of those, although I'd prefer it if the strain error was directly derived from the error of each UBI matrix element, if that makes sense:

def bootstrap_mean_grain_and_errors(grain_list: List[id11_grain], phase):
    trans_array = np.array([grain.translation for grain in grain_list])
    new_trans = np.mean(trans_array, axis=0)
    trans_stdev = np.std(trans_array, axis=0)

    ubis_fz = np.empty((len(grain_list), 3, 3))
    for i, grain in enumerate(grain_list):
        # From pymicro
        # Take the orientation from the grain, move it to the fundamental zone
        u_fz, symm_used = move_rotation_to_FZ(grain.U, phase.symmetry)
        # Put the new funamental zone orientation back into a UBI
        ubis_fz[i] = np.linalg.inv(np.matmul(u_fz, grain.B))

    new_ubi = np.mean(ubis_fz, axis=0)
    ubi_stdev = np.std(ubis_fz, axis=0)

    new_grain = id11_grain(ubi=new_ubi, translation=new_trans)
    return new_grain, trans_stdev, ubi_stdev

[jonwright]

Don't you just need to compute a strain for each UBI and then take the mean and stddev of those?

For grains on the edge of the fundamental zone you can get problems as you have it. You might want to transform them all to be as close as possible to the first one on the list. If that one is near the edge of the FZ you want the others to reference to this and not pop up on the other side of a boundary...

[jadball]

Hi Jon,

Yes, that was my first idea - however, I'm getting inconsistent results when I compare the mean strain across all grains with the strain from the mean UBI:

    trans_array = np.array([grain.translation for grain in grain_list])
    new_trans = np.mean(trans_array, axis=0)
    trans_stdev = np.std(trans_array, axis=0)

    # Get the first grain orientation
    first_grain_U = grain_list[0].U
    grain_list[0].rotated_U = first_grain_U

    num_symms = phase.symmetry.symmetry_operators().shape[0]

    for grain in grain_list[1:]:
        angle_array = np.empty(shape=num_symms)
        rotated_grain_u_matrices = np.empty(shape=(num_symms, 3, 3))
        for i, symm_op in enumerate(phase.symmetry.symmetry_operators()):
            # Rotate this grain by the symmetry op
            o_i = np.dot(symm_op, grain.U)
            # Save it to a list for this grain
            rotated_grain_u_matrices[i] = o_i
            # Work out the angle between the first grain U and this rotated U
            delta = np.dot(first_grain_U, o_i.T)
            cw = (np.trace(delta) - 1.0) / 2.0
            if cw > 1. and cw - 1. < 10 * np.finfo("float32").eps:
                cw = 1.
            angle_array[i] = np.arccos(cw)  # returns radians

        # Get the smallest angle
        index = np.argmin(angle_array)
        closest_grain_rotation = rotated_grain_u_matrices[index]
        grain.rotated_U = closest_grain_rotation

    ubis_fz = np.empty((len(grain_list), 3, 3))
    eps_symms = np.empty((len(grain_list), 3, 3))
    eps_s_symms = np.empty((len(grain_list), 3, 3))
    for i, grain in enumerate(grain_list):
        grain_ubi_rotated = np.linalg.inv(np.matmul(grain.rotated_U, grain.B))
        eps_symms[i] = grain.eps_grain_matrix(dzero_cell=phase.reference_unit_cell)
        eps_s_symms[i] = grain.eps_sample_matrix(dzero_cell=phase.reference_unit_cell)
        ubis_fz[i] = grain_ubi_rotated

    new_ubi = np.mean(ubis_fz, axis=0)
    new_grain = id11_grain(ubi=new_ubi, translation=new_trans)

    mean_eps = np.mean(eps_symms, axis=0)
    grain_eps = new_grain.eps_grain_matrix(dzero_cell=phase.reference_unit_cell)
    print(mean_eps)
    print(grain_eps)

mean_eps:

[[-5.60662627e-17  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00 -9.60342916e-17 -2.02172675e-17]
 [ 0.00000000e+00 -2.01684493e-17 -7.99360578e-17]]

grain_eps:

[[-4.01143655e-04 -9.71995429e-07 -3.33813398e-06]
 [-9.71995429e-07 -4.23778181e-04 -6.03768441e-07]
 [-3.33813398e-06 -6.03768441e-07 -4.10767546e-04]]

[jonwright]

Can you share an example of grain_list somewhere? Is it a list corresponding to the same grain in repeated bootstrap fits, or all of the grains in a sample?

The code in ImageD11.sym_u works a bit differently for symmetry. There was some underlying problem with the definition of U in a strained grain that was avoided by not using U. I wonder if this is related.

[jadball]

Sure thing - I'm generating grain_list from some normally distributed orientations and translations - basically some test data that's meant to look like repeated bootstrap fits. Here's the grains exported into a map: https://gist.github.com/jadball/53b0f07fde83856a3311e05bba9c0cf5

Here's the functions I'm using to generate the normally-distributed UBIs:

def norm_scipy_random_rotations(n_rots, sigma):
    normally_distributed_rotations = [None] * n_rots

    for i in range(n_rots):
        quat = np.array([1, np.random.normal(loc=0, scale=sigma), np.random.normal(loc=0, scale=sigma),
                         np.random.normal(loc=0, scale=sigma)])
        norm_quaternion = quat / np.linalg.norm(quat)
        norm_quaternion_scipy_rot = scipy.Rotation.from_quat(norm_quaternion)
        normally_distributed_rotations[i] = norm_quaternion_scipy_rot

    return normally_distributed_rotations

def generate_random_ubis(n_ubis, u_sigma, phase):
    reference_unit_cell = np.array([3, 3, 3, 90, 90, 90])
    random_rotations = norm_scipy_random_rotations(n_ubis, u_sigma)
    random_rotation_matrices = np.array([rot.as_matrix() for rot in random_rotations])
    random_ubis = [xfab.tools.u_to_ubi(orien, reference_unit_cell) for orien in random_rotation_matrices]
    return random_ubis

And here's the full function i'm using to try and merge the observations and get errors:

def bootstrap_mean_grain_and_errors(grain_list: List[id11_grain], phase):
    trans_array = np.array([grain.translation for grain in grain_list])
    new_trans = np.mean(trans_array, axis=0)
    trans_stdev = np.std(trans_array, axis=0)

    # Get the first grain orientation
    first_grain_U = grain_list[0].U
    grain_list[0].rotated_U = first_grain_U

    # This is just all the cubic symmetry operators in an array:
    num_symms = phase.symmetry.symmetry_operators().shape[0]

    for grain in grain_list[1:]:
        angle_array = np.empty(shape=num_symms)
        rotated_grain_u_matrices = np.empty(shape=(num_symms, 3, 3))
        for i, symm_op in enumerate(phase.symmetry.symmetry_operators()):
            # Rotate this grain by the symmetry op
            o_i = np.dot(symm_op, grain.U)
            # Save it to a list for this grain
            rotated_grain_u_matrices[i] = o_i
            # Work out the angle between the first grain U and this rotated U
            delta = np.dot(first_grain_U, o_i.T)
            cw = (np.trace(delta) - 1.0) / 2.0
            if cw > 1. and cw - 1. < 10 * np.finfo("float32").eps:
                cw = 1.
            angle_array[i] = np.arccos(cw)  # returns radians

        # Get the smallest angle
        index = np.argmin(angle_array)
        closest_grain_rotation = rotated_grain_u_matrices[index]
        grain.rotated_U = closest_grain_rotation

    ubis_fz = np.empty((len(grain_list), 3, 3))
    eps_symms = np.empty((len(grain_list), 3, 3))
    eps_s_symms = np.empty((len(grain_list), 3, 3))
    for i, grain in enumerate(grain_list):
        grain_ubi_rotated = np.linalg.inv(np.matmul(grain.rotated_U, grain.B))
        eps_symms[i] = grain.eps_grain_matrix(dzero_cell=phase.reference_unit_cell)
        eps_s_symms[i] = grain.eps_sample_matrix(dzero_cell=phase.reference_unit_cell)
        ubis_fz[i] = grain_ubi_rotated

    new_ubi = np.mean(ubis_fz, axis=0)
    new_grain = id11_grain(ubi=new_ubi, translation=new_trans)

    eps_error = np.std(eps_symms, axis=0)
    eps_s_error = np.std(eps_s_symms, axis=0)

    return new_grain, trans_stdev, eps_error, eps_s_error

[jonwright]

Am I right in thinking that you expect zero strain for all of these grains? I will try to take a look later. This looks like all strains would be zero, just a variety of orientations:

random_ubis = [xfab.tools.u_to_ubi(orien, reference_unit_cell) for orien in random_rotation_matrices]

[jadball]

@jonwright yes that's correct - I'm considering them all strain-free for now while we troubleshoot this

[jonwright]

I think you are seeing an effect like bond libration. The lattice vectors (ubi) are rotated onto a spherical surface. The mean vector is inside this surface. So this mean shows a volume contraction. The amount of contraction should be proportional to the angular spread.

It seems the angular spread is about 2.5 degrees in your simulation? This is big enough to give the cross talk as I think you are beyond the small perturbation (linear) regime?

Real data should be within about 0.1 degrees if all went well?

[jadball]

@jonwright that makes sense, many thanks - decreasing the sigma in the norm_scipy_random_rotations function has decreased the angular spread and reduced the new strain value. It's good to know the constraint re: angular spread - what fraction of the peaks per grain would you recommend keeping? at the moment, for each sample I am using a fraction of 0.5, with around 100 samples per grain in total.

[jonwright]

This question is probably the reason we are really supposed to replace the peaks that are removed with computed peaks. Then you can guess that you don't want to replace "too many" of them.

If errors are only due to noise, when you have 2 times the number of peaks, the noise might be about sqrt(2) smaller. Doubtful that this is true in practice. It depends on how many peaks there were. You want to keep more observations than parameters, and there are 12 parameters and 1 or 2 good observations per peak (tth, eta, but omega is usually noisy). With >200 peaks you might get away with using half of them. With only 10 peaks then I would get very nervous trying to fit against 5. Does it depend on how many rings you measured perhaps?

My interest in this approach was to find out whether the resulting distribution is Gaussian. If some of the spots are overlapped, or assigned to the wrong grain, then I expect the fit to change a lot when they are removed. Usually, these spots show up as outliers anyway.

Doing 10 fits, each one excluding a different 10% of the data, seemed like a reasonable approach to me. Putting in replacement peaks might be quite easy to do. But it didn't seem worth the hassle for a 10% effect. Replacement stops the noise from inflating when you are running out of data and aims to get at a more 'standard error' by keeping the same number of observations. But running out of data usually means the grain has other problems anyway...