# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the Lesser GNU Public Licence, v2.1 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#

"""Streamplots (3D) --- :mod:`MDAnalysis.visualization.streamlines_3D`
===================================================================

:Authors: Tyler Reddy and Matthieu Chavent
:Year: 2014
:Copyright: Lesser GNU Public License v2.1+


The :func:`generate_streamlines_3d` function can generate a 3D flow field from
a MD trajectory, for instance, lipid molecules in a virus capsid. It can make
use of multiple cores to perform the analyis in parallel (using
:mod:`multiprocessing`).

.. rubric: References

.. footbibliography::

See Also
--------
MDAnalysis.visualization.streamlines : streamplots in 2D


.. autofunction:: generate_streamlines_3d

"""
import multiprocessing

import numpy as np
import numpy.testing
import scipy
import scipy.spatial.distance

import MDAnalysis


def determine_container_limits(
    topology_file_path, trajectory_file_path, buffer_value
):
    """Calculate the extent of the atom coordinates + buffer.

    A function for the parent process which should take the input trajectory
    and calculate the limits of the container for the system and return these
    limits.

    Parameters
    ----------
    topology_file_path : str
        topology file name
    trajectory_file_path : str
        trajectory file name
    buffer_value : float
        buffer value (padding) in +/- {x, y, z}
    """
    universe_object = MDAnalysis.Universe(
        topology_file_path, trajectory_file_path
    )
    all_atom_selection = universe_object.select_atoms(
        "all"
    )  # select all particles
    all_atom_coordinate_array = all_atom_selection.positions
    x_min, x_max, y_min, y_max, z_min, z_max = [
        all_atom_coordinate_array[..., 0].min(),
        all_atom_coordinate_array[..., 0].max(),
        all_atom_coordinate_array[..., 1].min(),
        all_atom_coordinate_array[..., 1].max(),
        all_atom_coordinate_array[..., 2].min(),
        all_atom_coordinate_array[..., 2].max(),
    ]
    tuple_of_limits = (
        x_min - buffer_value,
        x_max + buffer_value,
        y_min - buffer_value,
        y_max + buffer_value,
        z_min - buffer_value,
        z_max + buffer_value,
    )  # using buffer_value to catch particles near edges
    return tuple_of_limits


def produce_grid(tuple_of_limits, grid_spacing):
    """Produce a 3D grid for the simulation system.

    The partitioning is based on the tuple of Cartesian Coordinate limits
    calculated in an earlier step.

    Parameters
    ----------
    tuple_of_limits : tuple
        ``x_min, x_max, y_min, y_max, z_min, z_max``
    grid_spacing : float
        grid size in all directions in ångström

    Returns
    -------
    grid : array
        ``numpy.mgrid[x_min:x_max:grid_spacing, y_min:y_max:grid_spacing, z_min:z_max:grid_spacing]``

    """
    x_min, x_max, y_min, y_max, z_min, z_max = tuple_of_limits
    grid = np.mgrid[
        x_min:x_max:grid_spacing,
        y_min:y_max:grid_spacing,
        z_min:z_max:grid_spacing,
    ]
    return grid


def split_grid(grid, num_cores):
    """Split the grid into blocks of vertices.

    Take the overall `grid` for the system and split it into lists of cube
    vertices that can be distributed to each core.

    Parameters
    ----------
    grid : numpy.array
        3D array
    num_cores : int
        number of partitions to generate

    Returns
    -------
    list_dictionaries_for_cores : list of dict
    total_cubes : int
    num_sheets : int
    delta_array_shape : tuple

    """
    # unpack the x,y,z mgrid arrays
    x, y, z = grid
    num_z_values = z.shape[-1]
    num_sheets = z.shape[0]
    delta_array_shape = tuple(
        [n - 1 for n in x.shape]
    )  # the final target shape for return delta arrays is n-1 in each dimension

    ordered_list_per_sheet_x_values = []
    for (
        x_sheet
    ) in (
        x
    ):  # each x_sheet should have shape (25,23) and the same x value in each element
        array_all_x_values_current_sheet = x_sheet.flatten()
        ordered_list_per_sheet_x_values.append(
            array_all_x_values_current_sheet
        )
    ordered_list_per_sheet_y_values = []
    for y_columns in y:
        array_all_y_values_current_sheet = y_columns.flatten()
        ordered_list_per_sheet_y_values.append(
            array_all_y_values_current_sheet
        )
    ordered_list_per_sheet_z_values = []
    for z_slices in z:
        array_all_z_values_current_sheet = z_slices.flatten()
        ordered_list_per_sheet_z_values.append(
            array_all_z_values_current_sheet
        )

    ordered_list_cartesian_coordinates_per_sheet = []
    for x_sheet_coords, y_sheet_coords, z_sheet_coords in zip(
        ordered_list_per_sheet_x_values,
        ordered_list_per_sheet_y_values,
        ordered_list_per_sheet_z_values,
    ):
        ordered_list_cartesian_coordinates_per_sheet.append(
            list(zip(x_sheet_coords, y_sheet_coords, z_sheet_coords))
        )
    array_ordered_cartesian_coords_per_sheet = np.array(
        ordered_list_cartesian_coordinates_per_sheet
    )
    # now I'm going to want to build cubes in an ordered fashion, and in such a way that I can track the index /
    # centroid of each cube for domain decomposition / reconstruction and mayavi mlab.flow() input
    # cubes will be formed from N - 1 base sheets combined with subsequent sheets
    current_base_sheet = 0
    dictionary_cubes_centroids_indices = {}
    cube_counter = 0
    while current_base_sheet < num_sheets - 1:
        current_base_sheet_array = array_ordered_cartesian_coords_per_sheet[
            current_base_sheet
        ]
        current_top_sheet_array = array_ordered_cartesian_coords_per_sheet[
            current_base_sheet + 1
        ]  # the points of the sheet 'to the right' in the grid
        current_index = 0
        while current_index < current_base_sheet_array.shape[0] - num_z_values:
            # iterate through all the indices in each of the sheet arrays (careful to avoid extra
            # points not needed for cubes)
            column_z_level = (
                0  # start at the bottom of a given 4-point column and work up
            )
            while column_z_level < num_z_values - 1:
                current_list_cube_vertices = []
                first_two_vertices_base_sheet = current_base_sheet_array[
                    current_index : current_index + 2, ...
                ].tolist()
                first_two_vertices_top_sheet = current_top_sheet_array[
                    current_index : current_index + 2, ...
                ].tolist()
                next_two_vertices_base_sheet = current_base_sheet_array[
                    current_index
                    + num_z_values : 2
                    + num_z_values
                    + current_index,
                    ...,
                ].tolist()
                next_two_vertices_top_sheet = current_top_sheet_array[
                    current_index
                    + num_z_values : 2
                    + num_z_values
                    + current_index,
                    ...,
                ].tolist()
                for vertex_set in [
                    first_two_vertices_base_sheet,
                    first_two_vertices_top_sheet,
                    next_two_vertices_base_sheet,
                    next_two_vertices_top_sheet,
                ]:
                    current_list_cube_vertices.extend(vertex_set)
                vertex_array = np.array(current_list_cube_vertices)
                assert vertex_array.shape == (
                    8,
                    3,
                ), "vertex_array has incorrect shape"
                cube_centroid = np.average(
                    np.array(current_list_cube_vertices), axis=0
                )
                dictionary_cubes_centroids_indices[cube_counter] = {
                    "centroid": cube_centroid,
                    "vertex_list": current_list_cube_vertices,
                }
                cube_counter += 1
                current_index += 1
                column_z_level += 1
                if (
                    column_z_level == num_z_values - 1
                ):  # the loop will break but I should also increment the
                    # current_index
                    current_index += 1
        current_base_sheet += 1
    total_cubes = len(dictionary_cubes_centroids_indices)

    # produce an array of pseudo cube indices (actually the dictionary keys which are cube numbers in string format):
    pseudo_cube_indices = np.arange(0, total_cubes)
    sublist_of_cube_indices_per_core = np.array_split(
        pseudo_cube_indices, num_cores
    )
    # now, the split of pseudoindices seems to work well, and the above sublist_of_cube_indices_per_core is a list of
    # arrays of cube numbers / keys in the original dictionary
    # now I think I'll try to produce a list of dictionaries that each contain their assigned cubes based on the above
    #  per core split
    list_dictionaries_for_cores = []
    subdictionary_counter = 0
    for array_cube_indices in sublist_of_cube_indices_per_core:
        current_core_dictionary = {}
        items_to_pop = len(array_cube_indices)
        items_popped = 0
        while items_popped < items_to_pop:
            key, value = dictionary_cubes_centroids_indices.popitem()
            current_core_dictionary.update({key: value})
            items_popped += 1
        list_dictionaries_for_cores.append(current_core_dictionary)
        subdictionary_counter += 1
    return (
        list_dictionaries_for_cores,
        total_cubes,
        num_sheets,
        delta_array_shape,
    )


def per_core_work(
    start_frame_coord_array,
    end_frame_coord_array,
    dictionary_cube_data_this_core,
    MDA_selection,
    start_frame,
    end_frame,
):
    """Run the analysis on one core.

    The code to perform on a given core given the dictionary of cube data.
    """
    list_previous_frame_centroids = []
    list_previous_frame_indices = []
    # define some utility functions for trajectory iteration:

    def point_in_cube(
        array_point_coordinates, list_cube_vertices, cube_centroid
    ):
        """Determine if an array of coordinates are within a cube."""
        # the simulation particle point can't be more than half the cube side length away from the cube centroid in
        # any given dimension:
        array_cube_vertices = np.array(list_cube_vertices)
        cube_half_side_length = (
            scipy.spatial.distance.pdist(
                array_cube_vertices, "euclidean"
            ).min()
            / 2.0
        )
        array_cube_vertex_distances_from_centroid = (
            scipy.spatial.distance.cdist(
                array_cube_vertices, cube_centroid[np.newaxis, :]
            )
        )
        np.testing.assert_allclose(
            array_cube_vertex_distances_from_centroid.min(),
            array_cube_vertex_distances_from_centroid.max(),
            rtol=0,
            atol=1.5e-4,
            err_msg="not all cube vertex to centroid distances are the same, "
            "so not a true cube",
        )
        absolute_delta_coords = np.absolute(
            np.subtract(array_point_coordinates, cube_centroid)
        )
        absolute_delta_x_coords = absolute_delta_coords[..., 0]
        indices_delta_x_acceptable = np.where(
            absolute_delta_x_coords <= cube_half_side_length
        )
        absolute_delta_y_coords = absolute_delta_coords[..., 1]
        indices_delta_y_acceptable = np.where(
            absolute_delta_y_coords <= cube_half_side_length
        )
        absolute_delta_z_coords = absolute_delta_coords[..., 2]
        indices_delta_z_acceptable = np.where(
            absolute_delta_z_coords <= cube_half_side_length
        )
        intersection_xy_acceptable_arrays = np.intersect1d(
            indices_delta_x_acceptable[0], indices_delta_y_acceptable[0]
        )
        overall_indices_points_in_current_cube = np.intersect1d(
            intersection_xy_acceptable_arrays, indices_delta_z_acceptable[0]
        )
        return overall_indices_points_in_current_cube

    def update_dictionary_point_in_cube_start_frame(
        array_simulation_particle_coordinates, dictionary_cube_data_this_core
    ):
        """Basically update the cube dictionary objects assigned to this core to contain a new key/value pair
        corresponding to the indices of the relevant particles that fall within a given cube. Also, for a given cube,
        store a key/value pair for the centroid of the particles that fall within the cube.
        """
        cube_counter = 0
        for key, cube in dictionary_cube_data_this_core.items():
            index_list_in_cube = point_in_cube(
                array_simulation_particle_coordinates,
                cube["vertex_list"],
                cube["centroid"],
            )
            cube["start_frame_index_list_in_cube"] = index_list_in_cube
            if (
                len(index_list_in_cube) > 0
            ):  # if there's at least one particle in this cube
                centroid_particles_in_cube = np.average(
                    array_simulation_particle_coordinates[index_list_in_cube],
                    axis=0,
                )
                cube["centroid_of_particles_first_frame"] = (
                    centroid_particles_in_cube
                )
            else:  # empty cube
                cube["centroid_of_particles_first_frame"] = None
            cube_counter += 1

    def update_dictionary_end_frame(
        array_simulation_particle_coordinates, dictionary_cube_data_this_core
    ):
        """Update the cube dictionary objects again as appropriate for the second and final frame."""
        cube_counter = 0
        for key, cube in dictionary_cube_data_this_core.items():
            # if there were no particles in the cube in the first frame, then set dx,dy,dz each to 0
            if cube["centroid_of_particles_first_frame"] is None:
                cube["dx"] = 0
                cube["dy"] = 0
                cube["dz"] = 0
            else:  # there was at least one particle in the starting cube so we can get dx,dy,dz centroid values
                new_coordinate_array_for_particles_starting_in_this_cube = (
                    array_simulation_particle_coordinates[
                        cube["start_frame_index_list_in_cube"]
                    ]
                )
                new_centroid_for_particles_starting_in_this_cube = np.average(
                    new_coordinate_array_for_particles_starting_in_this_cube,
                    axis=0,
                )
                cube["centroid_of_paticles_final_frame"] = (
                    new_centroid_for_particles_starting_in_this_cube
                )
                delta_centroid_array_this_cube = (
                    new_centroid_for_particles_starting_in_this_cube
                    - cube["centroid_of_particles_first_frame"]
                )
                cube["dx"] = delta_centroid_array_this_cube[0]
                cube["dy"] = delta_centroid_array_this_cube[1]
                cube["dz"] = delta_centroid_array_this_cube[2]
            cube_counter += 1

    # now that the parent process is dealing with the universe object & grabbing required coordinates, each child
    # process only needs to take the coordinate arrays & perform the operations with its assigned cubes (no more file
    #  opening and trajectory iteration on each core--which I'm hoping will substantially reduce the physical memory
    # footprint of my 3D streamplot code)
    update_dictionary_point_in_cube_start_frame(
        start_frame_coord_array, dictionary_cube_data_this_core
    )
    update_dictionary_end_frame(
        end_frame_coord_array, dictionary_cube_data_this_core
    )
    return dictionary_cube_data_this_core


def produce_coordinate_arrays_single_process(
    topology_file_path,
    trajectory_file_path,
    MDA_selection,
    start_frame,
    end_frame,
):
    """Generate coordinate arrays.

    To reduce memory footprint produce only a single MDA selection and get
    desired coordinate arrays; can later send these coordinate arrays to all
    child processes rather than having each child process open a trajectory and
    waste memory.

    """
    universe_object = MDAnalysis.Universe(
        topology_file_path, trajectory_file_path
    )
    relevant_particles = universe_object.select_atoms(MDA_selection)
    # pull out coordinate arrays from desired frames:
    for ts in universe_object.trajectory:
        if ts.frame > end_frame:
            break  # stop here
        if ts.frame == start_frame:
            start_frame_relevant_particle_coordinate_array_xyz = (
                relevant_particles.positions
            )
        elif ts.frame == end_frame:
            end_frame_relevant_particle_coordinate_array_xyz = (
                relevant_particles.positions
            )
        else:
            continue
    return (
        start_frame_relevant_particle_coordinate_array_xyz,
        end_frame_relevant_particle_coordinate_array_xyz,
    )


def generate_streamlines_3d(
    topology_file_path,
    trajectory_file_path,
    grid_spacing,
    MDA_selection,
    start_frame,
    end_frame,
    xmin,
    xmax,
    ymin,
    ymax,
    zmin,
    zmax,
    maximum_delta_magnitude=2.0,
    num_cores="maximum",
):
    r"""Produce the x, y and z components of a 3D streamplot data set.

    Parameters
    ----------
    topology_file_path : str
            Absolute path to the topology file
    trajectory_file_path : str
            Absolute path to the trajectory file. It will normally be desirable
            to filter the trajectory with a tool such as GROMACS
            :program:`g_filter` (see :footcite:p:`Chavent2014`)
    grid_spacing : float
            The spacing between grid lines (angstroms)
    MDA_selection : str
            MDAnalysis selection string
    start_frame : int
            First frame number to parse
    end_frame : int
            Last frame number to parse
    xmin : float
            Minimum coordinate boundary for x-axis (angstroms)
    xmax : float
            Maximum coordinate boundary for x-axis (angstroms)
    ymin : float
            Minimum coordinate boundary for y-axis (angstroms)
    ymax : float
            Maximum coordinate boundary for y-axis (angstroms)
    maximum_delta_magnitude : float
            Absolute value of the largest displacement tolerated for the
            centroid of a group of particles ( angstroms). Values above this
            displacement will not count in the streamplot (treated as
            excessively large displacements crossing the periodic boundary)
    num_cores : int or 'maximum' (optional)
            The number of cores to use. (Default 'maximum' uses all available
            cores)

    Returns
    -------
    dx_array : array of floats
            An array object containing the displacements in the x direction
    dy_array : array of floats
            An array object containing the displacements in the y direction
    dz_array : array of floats
            An array object containing the displacements in the z direction

    Examples
    --------
    Generate 3D streamlines and visualize in `mayavi`_::

        import numpy as np

        import MDAnalysis
        import MDAnalysis.visualization.streamlines_3D

        import mayavi, mayavi.mlab

        # assign coordinate system limits and grid spacing:
        x_lower,x_upper = -8.73, 1225.96
        y_lower,y_upper = -12.58, 1224.34
        z_lower,z_upper = -300, 300
        grid_spacing_value = 20

        x1, y1, z1 = MDAnalysis.visualization.streamlines_3D.generate_streamlines_3d(
                        'testing.gro', 'testing_filtered.xtc',
                         xmin=x_lower, xmax=x_upper,
                         ymin=y_lower, ymax=y_upper,
                         zmin=z_lower, zmax=z_upper,
                         grid_spacing=grid_spacing_value, MDA_selection = 'name PO4',
                         start_frame=2, end_frame=3, num_cores='maximum')

        x, y, z = np.mgrid[x_lower:x_upper:x1.shape[0]*1j,
                          y_lower:y_upper:y1.shape[1]*1j,
                          z_lower:z_upper:z1.shape[2]*1j]

        # plot with mayavi:
        fig = mayavi.mlab.figure(bgcolor=(1.0, 1.0, 1.0), size=(800, 800), fgcolor=(0, 0, 0))
        for z_value in np.arange(z_lower, z_upper, grid_spacing_value):
            st = mayavi.mlab.flow(x, y, z, x1, y1, z1, line_width=1,
                                  seedtype='plane', integration_direction='both')
            st.streamline_type = 'tube'
            st.tube_filter.radius = 2
            st.seed.widget.origin = np.array([ x_lower,  y_upper,   z_value])
            st.seed.widget.point1 = np.array([ x_upper, y_upper,  z_value])
            st.seed.widget.point2 = np.array([ x_lower, y_lower,  z_value])
            st.seed.widget.resolution = int(x1.shape[0])
            st.seed.widget.enabled = False
        mayavi.mlab.axes(extent = [0, 1200, 0, 1200, -300, 300])
        fig.scene.z_plus_view()
        mayavi.mlab.savefig('test_streamplot_3D.png')
        # more compelling examples can be produced for vesicles and other spherical systems

    .. image:: test_streamplot_3D.png

    See Also
    --------
    MDAnalysis.visualization.streamlines.generate_streamlines


    .. _mayavi: http://docs.enthought.com/mayavi/mayavi/
    """
    # work out the number of cores to use:
    if num_cores == "maximum":
        num_cores = multiprocessing.cpu_count()  # use all available cores
    else:
        num_cores = num_cores  # use the value specified by the user
        # assert isinstance(num_cores,(int,long)), "The number of specified cores must (of course) be an integer."
    np.seterr(all="warn", over="raise")
    parent_cube_dictionary = {}  # collect all data from child processes here

    def log_result_to_parent(process_dict):
        parent_cube_dictionary.update(process_dict)

    # step 1: produce tuple of cartesian coordinate limits for the first frame
    # tuple_of_limits = determine_container_limits(topology_file_path = topology_file_path,trajectory_file_path =
    # trajectory_file_path,buffer_value=buffer_value)
    tuple_of_limits = (xmin, xmax, ymin, ymax, zmin, zmax)
    # step 2: produce a suitable grid (will assume that grid size / container size does not vary during simulation--or
    #  at least not beyond the buffer limit, such that this grid can be used for all subsequent frames)
    grid = produce_grid(
        tuple_of_limits=tuple_of_limits, grid_spacing=grid_spacing
    )
    # step 3: split the grid into a dictionary of cube information that can be sent to each core for processing:
    list_dictionaries_for_cores, total_cubes, num_sheets, delta_array_shape = (
        split_grid(grid=grid, num_cores=num_cores)
    )
    # step 3b: produce required coordinate arrays on a single core to avoid making a universe object on each core:
    start_frame_coord_array, end_frame_coord_array = (
        produce_coordinate_arrays_single_process(
            topology_file_path,
            trajectory_file_path,
            MDA_selection,
            start_frame,
            end_frame,
        )
    )
    # step 4: per process work using the above grid data split
    pool = multiprocessing.Pool(num_cores)
    for sub_dictionary_of_cube_data in list_dictionaries_for_cores:
        pool.apply_async(
            per_core_work,
            args=(
                start_frame_coord_array,
                end_frame_coord_array,
                sub_dictionary_of_cube_data,
                MDA_selection,
                start_frame,
                end_frame,
            ),
            callback=log_result_to_parent,
        )
    pool.close()
    pool.join()
    # so, at this stage the parent process now has a single dictionary with all the cube objects updated from all
    # available cores
    # the 3D streamplot (i.e, mayavi flow() function) will require separate 3D np arrays for dx,dy,dz
    # the shape of each 3D array will unfortunately have to match the mgrid data structure (bit of a pain): (
    # num_sheets - 1, num_sheets - 1, cubes_per_column)
    cubes_per_sheet = int(float(total_cubes) / float(num_sheets - 1))
    # produce dummy zero arrays for dx,dy,dz of the appropriate shape:
    dx_array = np.zeros(delta_array_shape)
    dy_array = np.zeros(delta_array_shape)
    dz_array = np.zeros(delta_array_shape)
    # now use the parent cube dictionary to correctly substitute in dx,dy,dz values
    current_sheet = 0  # which is also the current row
    y_index_current_sheet = 0  # sub row
    z_index_current_column = 0  # column
    total_cubes_current_sheet = 0
    for cube_number in range(0, total_cubes):
        dx_array[
            current_sheet, y_index_current_sheet, z_index_current_column
        ] = parent_cube_dictionary[cube_number]["dx"]
        dy_array[
            current_sheet, y_index_current_sheet, z_index_current_column
        ] = parent_cube_dictionary[cube_number]["dy"]
        dz_array[
            current_sheet, y_index_current_sheet, z_index_current_column
        ] = parent_cube_dictionary[cube_number]["dz"]
        z_index_current_column += 1
        total_cubes_current_sheet += 1
        if z_index_current_column == delta_array_shape[2]:
            # done building current y-column so iterate y value and reset z
            z_index_current_column = 0
            y_index_current_sheet += 1
            if (
                y_index_current_sheet == delta_array_shape[1]
            ):  # current sheet is complete
                current_sheet += 1
                y_index_current_sheet = 0  # restart for new sheet
                z_index_current_column = 0
                total_cubes_current_sheet = 0
    # now set velocity component values greater than a certain cutoff to 0,
    # because they tend to reflect spurious values (i.e., PBC jumping)
    dx_array[abs(dx_array) >= maximum_delta_magnitude] = 1.0
    dy_array[abs(dy_array) >= maximum_delta_magnitude] = 1.0
    dz_array[abs(dz_array) >= maximum_delta_magnitude] = 1.0
    return (dx_array, dy_array, dz_array)