Source code for

# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2022 by the BurnMan team, released under the GNU
# GPL v2 or later.

from __future__ import absolute_import

import numpy as np
from sympy import Matrix
from ..classes.combinedmineral import CombinedMineral
from ..classes.solution import Solution
from ..utils.chemistry import site_occupancies_to_strings

def _decompose_3D_matrix(W):
    Decomposes a 3D matrix W_ijk where E = W_ijk p_i p_j p_k
    into a subregular form where
    E = G_i p_i + WB_ij (1 - p_j + p_i) / 2 + WT_ijk p_i p_j p_k,
    and i < j < k.

    W : 3D numpy array

    new_endmember_excesses : 1D numpy array
        The array G_i

    new_binary_matrix : 2D numpy array
        The upper triangular matrix WB_ij

    new_ternary_terms : list of lists of length 4
        A list where each item is in the form [i, j, k, WT_ijk]

    n_mbrs = len(W)
    # New endmember components
    # W_iii needs to be copied, otherwise just a view onto W
    new_endmember_excesses = np.copy(np.einsum('iii->i', W))

    # Removal of endmember components from 3D representation
    W -= (np.einsum('i, j, k->ijk',
                    new_endmember_excesses, np.ones(n_mbrs),
          + np.einsum('i, j, k->ijk',
                      np.ones(n_mbrs), new_endmember_excesses,
          + np.einsum('i, j, k->ijk',
                      np.ones(n_mbrs), np.ones(n_mbrs),

    # Transformed 2D components
    # (i=j, i=k, j=k)
    new_binary_matrix = (np.einsum('jki, jk -> ij', W, np.identity(n_mbrs))
                         + np.einsum('jik, jk -> ij', W, np.identity(n_mbrs))
                         + np.einsum('ijk, jk -> ij', W,

    # Wb is the 3D matrix corresponding to the terms in the binary matrix,
    # such that the two following print statements produce the same answer
    # for a given array of endmember proportions
    Wb = (np.einsum('ijk, ij->ijk', W, np.identity(n_mbrs))
          + np.einsum('ijk, jk->ijk', W, np.identity(n_mbrs))
          + np.einsum('ijk, ik->ijk', W, np.identity(n_mbrs)))

    # Remove binary component from 3D representation
    # The extra terms are needed because the binary term in the formulation
    # of a subregular solution model given by
    # Helffrich and Wood includes ternary components (the sum_k X_k part)..
    W -= Wb + (np.einsum('ij, k', new_binary_matrix, np.ones(n_mbrs))
               - np.einsum('ij, ik->ijk',
                           new_binary_matrix, np.identity(n_mbrs))
               - np.einsum('ij, jk->ijk',
                           new_binary_matrix, np.identity(n_mbrs)))/2.

    # Find the 3D components Wijk by adding the elements at
    # the six equivalent positions in the matrix
    new_ternary_terms = []
    for i in range(n_mbrs):
        for j in range(i+1, n_mbrs):
            for k in range(j+1, n_mbrs):
                val = (W[i, j, k] + W[j, k, i]
                       + W[k, i, j] + W[k, j, i]
                       + W[j, i, k] + W[i, k, j]).round(decimals=12)
                if np.abs(val) > 1.e-12:
                    new_ternary_terms.append([i, j, k, val])

    return (new_endmember_excesses, new_binary_matrix, new_ternary_terms)

def _subregular_matrix_conversion(new_basis, binary_matrix,
                                  ternary_terms=None, endmember_excesses=None):
    Converts the arrays reguired to describe a subregular solution
    from one endmember basis to another.

    The excess nonconfigurational energies of the subregular solution model
    are described as follows:
    E = G_i p_i + WB_ij (1 - p_j + p_i) / 2 + WT_ijk p_i p_j p_k,
    and i < j < k.

    new_basis : 2D numpy array
        The new endmember basis, given as amounts of the old endmembers.

    binary_matrix : 2D numpy array
        The upper triangular matrix WB_ij

    ternary_terms : list of lists of length 4
        A list where each item is in the form [i, j, k, WT_ijk]

    endmember_excesses : 1D numpy array
        The array G_i

    new_endmember_excesses : 1D numpy array
        The array G_i

    new_binary_matrix : 2D numpy array
        The upper triangular matrix WB_ij

    new_ternary_terms : list of lists of length 4
        A list where each item is in the form [i, j, k, WT_ijk]
    n_mbrs = len(binary_matrix)
    # Compact 3D representation of original interactions
    W = (np.einsum('i, jk -> ijk', np.ones(n_mbrs), binary_matrix)
         + np.einsum('ij, jk -> ijk', binary_matrix, np.identity(n_mbrs))
         - np.einsum('ij, ik -> ijk', binary_matrix, np.identity(n_mbrs)))/2.

    # Add endmember components to 3D representation
    if endmember_excesses is not None:
        W += (np.einsum('i, j, k->ijk', endmember_excesses,
                        np.ones(n_mbrs), np.ones(n_mbrs))
              + np.einsum('j, i, k->ijk', endmember_excesses,
                          np.ones(n_mbrs), np.ones(n_mbrs))
              + np.einsum('k, i, j->ijk', endmember_excesses,
                          np.ones(n_mbrs), np.ones(n_mbrs)))/3.

    # Add ternary values to 3D representation
    if ternary_terms is not None:
        for i, j, k, val in ternary_terms:
            W[i, j, k] += val

    # Transformation to new 3D representation
    A = new_basis.T
    Wn = np.einsum('il, jm, kn, ijk -> lmn', A, A, A, W)

    new_endmember_excesses, new_binary_terms, new_ternary_terms = _decompose_3D_matrix(

    return (new_endmember_excesses, new_binary_terms, new_ternary_terms)

def complete_basis(basis):
    Creates a full basis by filling remaining rows with
    rows of the identity matrix with row indices not
    in the column pivot list of the basis RREF

    n, m = basis.shape
    if n < m:
        pivots = list(Matrix(basis).rref()[1])
        return np.concatenate((basis,
                               np.identity(m)[[i for i in range(m)
                                               if i not in pivots], :]),
        return basis

[docs]def transform_solution_to_new_basis(solution, new_basis, n_mbrs=None, solution_name=None, endmember_names=None, molar_fractions=None): """ Transforms a solution model from one endmember basis to another. Returns a new Solution object. Parameters ---------- solution : :class:`burnman.Solution` object The original solution object. new_basis : 2D numpy array The new endmember basis, given as amounts of the old endmembers. n_mbrs : float (optional) The number of endmembers in the new solution (defaults to the length of new_basis) solution_name : string (optional) A name corresponding to the new solution endmember_names : list of strings (optional) A list corresponding to the names of the new endmembers. molar_fractions : numpy array (optional) Fractions of the new endmembers in the new solution. Returns ------- solution : :class:`burnman.Solution` object The transformed solution """ new_basis = np.array(new_basis) if n_mbrs is None: n_mbrs, n_all_mbrs = new_basis.shape else: _, n_all_mbrs = new_basis.shape if solution_name is None: name = 'child solution' else: name = solution_name solution_type = solution.solution_type if solution_type == 'ideal': ESV_modifiers = [[0., 0., 0.] for v in new_basis] elif (solution_type == 'asymmetric' or solution_type == 'symmetric'): A = complete_basis(new_basis).T old_alphas = solution.solution_model.alphas alphas = np.einsum('i, ij', solution.solution_model.alphas, A) inv_diag_alphas = np.diag(1./alphas) B = np.einsum('ij, jk, kl->il', np.diag(old_alphas), A, inv_diag_alphas) alphas = list(alphas[0:n_mbrs]) Qe = np.einsum('ik, ij, kl->jl', solution.solution_model.We, B, B) Qs = np.einsum('ik, ij, kl->jl', solution.solution_model.Ws, B, B) Qv = np.einsum('ik, ij, kl->jl', solution.solution_model.Wv, B, B) def new_interactions(Q, n_mbrs): return [[float((Q[i, j] + Q[j, i] - Q[i, i] - Q[j, j]) * (alphas[i] + alphas[j])/2.) for j in range(i+1, n_mbrs)] for i in range(n_mbrs-1)] energy_interaction = new_interactions(Qe, n_mbrs) entropy_interaction = new_interactions(Qs, n_mbrs) volume_interaction = new_interactions(Qv, n_mbrs) ESV_modifiers = [[Qe[i, i] * alphas[i], Qs[i, i] * alphas[i], Qv[i, i] * alphas[i]] for i in range(n_mbrs)] elif solution_type == 'subregular': full_basis = complete_basis(new_basis) def new_interactions(W, n_mbrs): return [[[W[i, j], W[j, i]] for j in range(i+1, n_mbrs)] for i in range(n_mbrs-1)] # N.B. initial endmember_excesses are zero Emod, We, ternary_e = _subregular_matrix_conversion(full_basis, solution.solution_model.We, solution.solution_model.ternary_terms_e) Smod, Ws, ternary_s = _subregular_matrix_conversion(full_basis, solution.solution_model.Ws, solution.solution_model.ternary_terms_s) Vmod, Wv, ternary_v = _subregular_matrix_conversion(full_basis, solution.solution_model.Wv, solution.solution_model.ternary_terms_v) energy_interaction = new_interactions(We, n_mbrs) entropy_interaction = new_interactions(Ws, n_mbrs) volume_interaction = new_interactions(Wv, n_mbrs) ESV_modifiers = [[Emod[i], Smod[i], Vmod[i]] for i in range(n_mbrs)] else: raise Exception('The function to change basis for the ' '{0} solution model has not yet been ' 'implemented.'.format(solution_type)) # Create site formulae new_occupancies = np.array(new_basis).dot( solution.solution_model.endmember_occupancies) new_multiplicities = np.array(new_basis).dot( solution.solution_model.site_multiplicities) site_formulae = site_occupancies_to_strings(solution.solution_model.sites, new_multiplicities, new_occupancies) # Create endmembers endmembers = [] for i, vector in enumerate(new_basis): nonzero_indices = np.nonzero(vector)[0] if len(nonzero_indices) == 1: endmembers.append([solution.endmembers[nonzero_indices[0]][0], site_formulae[i]]) else: mbr = CombinedMineral([solution.endmembers[idx][0] for idx in nonzero_indices], [vector[idx] for idx in nonzero_indices], ESV_modifiers[i]) mbr.params['formula'] = {key: value for (key, value) in mbr.params['formula'].items() if value > 1.e-12} endmembers.append([mbr, site_formulae[i]]) if endmember_names is not None: for i in range(n_mbrs): endmembers[i][0].params['name'] = endmember_names[i] endmembers[i][0].name = endmember_names[i] if n_mbrs == 1: endmembers[0][0].name = name endmembers[0][0].parent = solution endmembers[0][0].basis = new_basis return endmembers[0][0] else: new_solution = Solution(name=name, solution_type=solution_type, endmembers=endmembers, energy_interaction=energy_interaction, volume_interaction=volume_interaction, entropy_interaction=entropy_interaction, alphas=alphas, molar_fractions=molar_fractions) new_solution.parent = solution new_solution.basis = new_basis return new_solution