..
    : PyBEST: Pythonic Black-box Electronic Structure Tool
    : Copyright (C) 2016-- The PyBEST Development Team
    :
    : This file is part of PyBEST.
    :
    : PyBEST is free software; you can redistribute it and/or
    : modify it under the terms of the GNU General Public License
    : as published by the Free Software Foundation; either version 3
    : of the License, or (at your option) any later version.
    :
    : PyBEST is distributed in the hope that it will be useful,
    : but WITHOUT ANY WARRANTY; without even the implied warranty of
    : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    : GNU General Public License for more details.
    :
    : You should have received a copy of the GNU General Public License
    : along with this program; if not, see <http://www.gnu.org/licenses/>
    :
    : --

.. _user_linalg_objects:

Matrix classes
##############

Methods implemented in PyBEST operate on and return custom multidimensional
array classes that store the actual arrays as their private attributes 
(_array). The linalg module contains two types of matrix representations:
dense matrices and Cholesky matrices. Algebraic operations are implemented
as methods.

DenseNIndex (N = One, Two, Three, Four) classes
===============================================

In the dense representation, all the elements are stored in memory. The dense
matrix can be created using DenseLinalgFactory or directly by using the constructor

.. code-block:: python

    my_matrix = DenseThreeIndex(n1, n2, n3)

where arguments n1, n2, and n3 are three integers defining the dimensionality
of the three-index object. DenseNIndex objects are usually returned as
products of mathematical operations.

CholeskyFourIndex class
=======================

The Cholesky Decomposition representation operates on the four-index object A
that is internally stored in the form of two three-index object which fulfills
A = LL*. The interface is analogous to the DenseFourIndex class interface, so
CholeskyFourIndex instances are treated as four-index arrays.

LinalgFactory classes
#####################

The DenseNIndex and CholeskyFourIndex matrices can be created using linear
algebra factory instances - DenseLinalgFactory and CholeskyLinagsFactory.
See examples to compare different methods of creating new N-index
multidimensional arrays.

Using PyBEST classes to perform algebraic operations
######################################################

Creation of dense multidimensional arrays
=========================================

.. code-block:: python

    from pybest.linalg.dense import DenseLinalgFactory, DenseTwoIndex

    # Create 9x9 and 9x9x9 dense matrices using factory instance.
    lf = DenseLinalgFactory(9)
    matrix2d_0 = lf.create_two_index()
    matrix3d_0 = lf.create_three_index()

    # Create 3x9 dense matrix using direct initialization.
    matrix2d_1 = DenseTwoIndex(3, 9)

    # Create an empty matrix with the same shape as matrix2d_0.
    matrix2d_2 = matrix2d_0.new()

    # Create a copy of the matrix.
    matrix2d_3 = matrix2d_0.copy()


Filling matrix with numbers
===========================

.. code-block:: python

    # Fill matrix with random numbers.
    my_matrix.randomize()

    # Fill diagonal elements of the subblock [3:6, 3:6] (defined by begin and
    # end keywords) of the matrix.
    my_matrix.assign_diagonal(2, begin=3, end=6).

    # Set all elements of a matrix to 0
    my_matrix.assign(0)

    # Set all elements to 3 in a subblock [2:4,1:3] of my_matrix.
    my_matrix.assign(3, begin0=2, end0=4, begin1=1, end1=3)
    print(my_matrix.get_element[2,2]) # Prints 3
    print(my_matrix.get_element[0,0]) # Prints 0

    # Fill matrix with the content of another matrix.
    my_matrix.assign(other_matrix)

    # Fill a subblock (defined using begin0, end0, begin1, end1 keywords) of
    # my_matrix with content of a subblock of another matrix (defined by
    # begin2, end2, begin3, end3 keywords).
    my_matrix.assign(other_matrix, begin0=2, begin2=2)


General NIndexObject multiplication - contraction
=================================================

Evaluates the Einstein summation convention on the operands that are
CholeskyFourIndex or DenseNIndex (where N = One, Two, Three, Four) instances.

Arguments
---------

    :subscript: (string)
          specifies the subscripts for summation as comma separated list
          of subscript labels, for example: 'abcd,efcd->abfe'.

    :operands: (DenseNIndex)
          these are the other arrays for the operation. If out keyword is
          not specified and the result of operation is not a scalar value,
          the last operand is treated as out.

Keyword arguments
-----------------

    :out: (DenseNIndex)
          The product of operation is added to out.

    :factor: (float)
          The product of operation is multiplied by a factor before it is
          added to out.

    :clear: (boolean)
           The out is cleared before the product of operation is added to
           out if clear is set to True.

    :select: (string)
           Specifies the contraction algorithm: One of:

               * 'opt_einsum' - opt_einsum.contract,
               * 'einsum' - numpy.einsum function,
               * 'td' - operations with utilization of numpy.tensordot function,
               * 'blas' - operations written in c++,
               * None - first tries 'td' or 'blas' routine, then "einsum"

           'td' is usually much faster but it requires 
           more memory and it is available only for specific cases.

    :optimize: (False, True, 'optimal')
           Specifies contraction path for operation if select is set to
           'einsum'. For details, see numpy.einsum_path.

    :begin0, end0, begin1, ...: (int)
           The lower and upper bounds for summation.

.. code-block:: python

    # Matrix-vector multiplication.
    lf = DenseLinalgFactory(4)
    my_vector = lf.create_one_index()
    my_matrix = lf.create_two_index()
    my_result = my_matrix.contract("ab,b", my_vector)

    # Four-index object contracted with the subblock of the two-index object.
    # output_2 is a transposed output_1.
    output_1 = array_4d.contract("abcd,cd->ab", array_2d)
    output_2 = array_4d.contract("abcd,cd->ba", array_2d)

    # Multiply a result of matrix-matrix multiplication by factor -1 and add it
    # to another matrix.
    matrix_0.contract("ab,ba", matrix_1, out=matrix_2, factor=-1)

    # There is no difference in syntax if you contract DenseFourIndex or
    # CholeskyFourIndex object.
    output_1 = dense4ind_object.contract("abcd,cd->ab", array_2d)
    output_2 = cholesky4ind_object.contract("abcd,cd->ab", array_2d)


Keyword arguments begin0, end0, begin1, ... can be used to perform contraction
on a subblock of an N-dimentional array. In the following example, we perform
the operation on the subblock of the electron repulsion integrals that are
represented as CholeskyFourIndex object. Please note, that there is no need to
use all of the begin0, end0, begin1, ... keyword arguments since their default
values allows to include all indices.

.. code-block:: python

    # nocc = number of occupied orbitals
    nocc = 5

    # subblock_ooov = defines boundaries of the subblock of the array
    subblock_ooov = {"end0": nocc, "end1": nocc, "end2": nocc, "begin3": nocc}
    
    # eri = electron repulsion integrals (CholeskyFourIndex instance)
    exchange_ooov = eri.contract("abcd->abcd", **subblock_ooov)
    exchange_ooov.iadd_transpose((0, 2, 1, 3), factor=-2)


General NIndexObject slicing
============================

The syntax for the slicing operation is analogous as for the contract method.

Arguments
---------

    :subscripts: (string)
          specifies the subscripts for summation as comma-separated list
          of subscript labels, for example: 'abcd,efcd->abfe'.

Keyword arguments
-----------------

    :out: (DenseNIndex)
          The product of operation is added to out.

    :factor: (float)
          The product of operation is multiplied by a factor before it is
          added to out.

    :clear: (boolean)
           The out is cleared before the product of operation is added to
           out if clear is set to True.

    :begin0, end0, begin1, ...: (int)
           The lower and upper bounds for summation.

.. code-block:: python

    # Get a slice
    array_3d = array_4d.contract("abcb->abc")
    
    # Get a slice from a given subblock
    subblock = {"begin0": 1, "end0": 9, "begin3": 1, "end3": 9}
    array_2d = array_4d.contract("abba->ab", **subblock)
    
    # Get a slice and add it to array_out after setting its values to zero
    array_4d.contract("abcb->abc", out=array_out, clear=True)


Adding arrays
=============

The kinetic energy operator and nuclear attraction energy operator are
represented by DenseTwoIndex objects. A nonrelativistic one-body energy
operator can be obtained by summing these two operators.

.. code-block:: python

    kinetic_energy_operator = compute_kinetic(obasis)
    nuclear_attraction_energy_operator = compute_nuclear(obasis)
    
    one_body_energy_operator = kinetic_energy_operator.copy()
    one_body_energy_operator.iadd(nuclear_attraction_energy_operator)

Scaling by a scalar factor
==========================

.. code-block:: python

    any_nindex_object.iscale(0.5)