/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-    */
/* ex: set filetype=cpp softtabstop=4 shiftwidth=4 tabstop=4 cindent expandtab: */

/*
  $Id: vctDynamicNArrayLoopEngines.h,v 1.4 2007/04/26 19:33:57 anton Exp $

  Author(s):	Daniel Li
  Created on:	2006-07-05

  (C) Copyright 2006-2007 Johns Hopkins University (JHU), All Rights
  Reserved.

--- begin cisst license - do not edit ---

This software is provided "as is" under an open source license, with
no warranty.  The complete license can be found in license.txt and
http://www.cisst.org/cisst/license.txt.

--- end cisst license ---
*/


#ifndef _vctDynamicNArrayLoopEngines_h
#define _vctDynamicNArrayLoopEngines_h

/*!
  \file
  \brief Declaration of vctDynamicNArrayLoopEngines
 */


#include <cisstCommon/cmnPortability.h>
#include <cisstCommon/cmnThrow.h>

#include <cisstVector/vctFixedSizeVector.h>
#include <cisstVector/vctContainerTraits.h>


/*!
  \brief Container class for the dynamic nArray engines.

  \sa SoNi SoNiNi SoNiSi NoNiNi NoNiSi NoSiNi NioSi NioNi NoNi Nio NioSiNi MinAndMax
*/
template <unsigned int _dimension>
class vctDynamicNArrayLoopEngines
{
public:
    /* define types */
    typedef unsigned int size_type;
    typedef int stride_type;
    typedef int difference_type;
    typedef unsigned int index_type;

    VCT_NARRAY_TRAITS_TYPEDEFS(_dimension);

	/*! Helper function to throw an exception whenever dimensions mismatch.
      This enforces that a standard message is sent. */
    inline static void ThrowDimensionMismatchException(void) throw(std::runtime_error)
    {
        cmnThrow(std::runtime_error("vctDynamicNArrayLoopEngines: Dimensions of nArrays don't match"));
    }


    /*! Helper function to throw an exception whenever sizes mismatch.
      This enforces that a standard message is sent. */
    inline static void ThrowSizeMismatchException(void) throw(std::runtime_error)
    {
        cmnThrow(std::runtime_error("vctDynamicNArrayLoopEngines: Sizes of nArrays don't match"));
    }


    /*! Helper function to calculate the strides to next dimension. */
    inline static void CalculateSTND(nstride_type & stnd,
                                     const nsize_type & sizes,
                                     const nstride_type & strides)
    {
        // set up iterators
        typename nsize_type::const_iterator sizesIter = sizes.begin();
        typename nstride_type::const_iterator stridesIter = strides.begin();
        typename nstride_type::iterator stndIter = stnd.begin();
        const typename nstride_type::const_iterator stndIterEnd = stnd.end();

        *stndIter = 0;
        ++sizesIter;
        ++stridesIter;
        ++stndIter;

        stride_type skippedStrides;
        for (; stndIter != stndIterEnd;
             ++sizesIter, ++stridesIter, ++stndIter)
        {
            skippedStrides = static_cast<stride_type>(*sizesIter) * (*stridesIter);
            *stndIter = *(stridesIter-1) - skippedStrides;
        }
    }


    /*! Helper function to calculate the pointer offsets to the next dimension. */
    inline static void CalculateOTND(nstride_type & otnd,
                                     const nstride_type & strides,
                                     const nstride_type & stnd)
    {
        // set up iterators
        stride_type previousOTND = *(strides.rbegin());
        typename nstride_type::const_reverse_iterator stndIter = stnd.rbegin();
        typename nstride_type::reverse_iterator otndIter = otnd.rbegin();
        const typename nstride_type::const_reverse_iterator otnd_rend = otnd.rend();

        *otndIter = previousOTND;
        ++otndIter;

        for (; otndIter != otnd_rend;
             ++otndIter, ++stndIter)
        {
            *otndIter = *stndIter + previousOTND;
            previousOTND = *otndIter;
        }
    }


    /*! Helper function to calculate the initial targets. */
    template <class _elementType>
    inline static void InitializeTargets(vctFixedSizeVector<const _elementType *, _dimension> & targets,
                                         const nsize_type & sizes,
                                         const nstride_type & strides,
                                         const _elementType * basePtr)
    {
        typedef _elementType value_type;

        // set up iterators
        typename nsize_type::const_iterator sizesIter = sizes.begin();
        typename nstride_type::const_iterator stridesIter = strides.begin();
        typename vctFixedSizeVector<const value_type *, _dimension>::iterator targetsIter = targets.begin();
        const typename vctFixedSizeVector<const value_type *, _dimension>::const_iterator targetsIterEnd = targets.end();

        stride_type offset;
        for (; targetsIter != targetsIterEnd;
             ++targetsIter, ++sizesIter, ++stridesIter)
        {
            offset = static_cast<stride_type>(*sizesIter) * (*stridesIter);
            *targetsIter = basePtr + offset;
        }
    }


    /*! Helper function to synchronize the given nArray's current pointer
      with the master nArray's current pointer. */
    template <class _elementType>
    inline static void SyncCurrentPointer(const _elementType * & currentPointer,
                                          const nstride_type & otnd,
                                          difference_type numberOfWrappedDimensions)
    {
        const typename nstride_type::const_reverse_iterator otndBegin = otnd.rbegin();
        currentPointer += otndBegin[numberOfWrappedDimensions];
    }


    /*! Helper function to increment the current pointer and any necessary
      target pointers. */
    template <class _elementType, class _pointerType>
    inline static dimension_type IncrementPointers(vctFixedSizeVector<const _elementType *, _dimension> & targets,
                                                    _pointerType & currentPointer,
                                                    const nstride_type & strides,
                                                    const nstride_type & stnd)
    {
        typedef _elementType value_type;

        // set up iterators
        typename vctFixedSizeVector<const value_type *, _dimension>::reverse_iterator targetsIter = targets.rbegin();
        const typename vctFixedSizeVector<const value_type *, _dimension>::const_reverse_iterator targets_rbeg = targets.rbegin();
        // typename vctFixedSizeVector<const value_type *, _dimension>::reverse_iterator targets_innerIter;
        typename nstride_type::const_reverse_iterator stridesIter = strides.rbegin();
        // const typename nstride_type::const_reverse_iterator strides_rend = strides.rend();
        typename nstride_type::const_reverse_iterator stndIter = stnd.rbegin();
        // typename nstride_type::const_reverse_iterator stnd_innerIter;
        dimension_type numberOfWrappedDimensions = 0;

        //* Below is Daniel's code
#if 0
        // update the current pointer
        currentPointer += *stridesIter;

        // update all necessary targets
        while (currentPointer == *targetsIter) {
            currentPointer += *stndIter;
            ++stridesIter;
            if (stridesIter == strides_rend) {
                ++numberOfWrappedDimensions;
                break;
            }
            *targetsIter += *stridesIter;

            targets_innerIter = targetsIter;
            stnd_innerIter = stndIter + 1;
            if (targets_innerIter != targets_rbeg) {
                do {
                    --targets_innerIter;
                    --stnd_innerIter;

                    *targets_innerIter += *stnd_innerIter;
                } while (targets_innerIter != targets_rbeg);
            }

            ++targetsIter;
            ++stndIter;
            ++numberOfWrappedDimensions;
        }
#else
        //* Below is Ofri's code
        // Consider the following structure.
        //
        // initial target setup :    // performed before engine loop is begun, outside of this function
        //       currPtr = base;
        //       for each k
        //           target[k] = currPtr + size[k] * stride[k];  // This should be invariant after every dimension wrap
        //           stnd[k] = stride[k-1] - size[k] * stride[k];   // This is constant through the run
        //
        // wrap(nw)  :      // nw is the number of dimensions wrapped
        //       currPtr += sum_{i =(d-nw)}^{d-1}( stnd[i] );  // obsolete if currPtr already determined; see below
        //       for k = (d-nw) to (d-1)
        //           target[k] = currPtr + size[k] * stride[k];
        //
        // check_wrap :     // this updates currentPtr (see above) and finds out nw (number of wrapped dimensions)
        //       k = d-1;
        //       nw = 0;
        //       currPtr += stride[k];
        //       while (currPtr == target[k]) {
        //           currPtr += stnd[k];
        //           --k;
        //           ++nw;
        //       }
        //       return nw;
        //
        // NOTE:  The invariant condition above and the initialization of stnd leads to
        //       size[k] * stride[k] = stride[k-1] - stnd[k];

        currentPointer += *stridesIter;
        while (currentPointer == *targetsIter) {
            currentPointer += *stndIter;
            ++targetsIter;
            ++stndIter;
            ++numberOfWrappedDimensions;
            if (numberOfWrappedDimensions == _dimension)
                return numberOfWrappedDimensions;
        }
        
        if (numberOfWrappedDimensions == 0)
            return numberOfWrappedDimensions;
        
        stridesIter += numberOfWrappedDimensions;
        int targetOffset;
        do {
            --targetsIter;
            --stndIter;
            targetOffset = *stridesIter - *stndIter;
            *targetsIter = currentPointer + targetOffset;
            --stridesIter;
        } while (targetsIter != targets_rbeg);
#endif
        return numberOfWrappedDimensions;
    }


    template <class _incrementalOperationType, class _elementOperationType>
    class SoNi
    {
    public:
        typedef typename _incrementalOperationType::OutputType OutputType;

        template <class _inputNArrayType>
        static OutputType Run(const _inputNArrayType & inputNArray)
        {
            typedef typename _inputNArrayType::value_type value_type;

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            vctFixedSizeVector<const value_type *, _dimension> inputTargets;
            const value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = inputNArray.dimension();

            CalculateSTND(inputSTND, inputSizes, inputStrides);
            InitializeTargets(inputTargets, inputSizes, inputStrides, inputPointer);

            OutputType incrementalResult = _incrementalOperationType::NeutralElement();

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                incrementalResult = _incrementalOperationType::Operate(
                    incrementalResult, _elementOperationType::Operate(*inputPointer) );

                numberOfWrappedDimensions =
                    IncrementPointers(inputTargets, inputPointer, inputStrides, inputSTND);
            }
            return incrementalResult;
        }   // Run method

    };  // SoNi class


    template <class _incrementalOperationType, class _elementOperationType>
    class SoNiNi
    {
    public:
        typedef typename _incrementalOperationType::OutputType OutputType;

        template <class _input1NArrayType, class _input2NArrayType>
        static OutputType Run(const _input1NArrayType & input1NArray,
                              const _input2NArrayType & input2NArray)
        {
            typedef typename _input1NArrayType::value_type value_type1;
            typedef typename _input2NArrayType::value_type value_type2;

            // declare all variables used for input1NArray
            const nsize_type & input1Sizes = input1NArray.sizes();
            const nstride_type & input1Strides = input1NArray.strides();
            nstride_type input1STND;
            vctFixedSizeVector<const value_type1 *, _dimension> input1Targets;
            const value_type1 * input1Pointer = input1NArray.Pointer();

            // declare all variables used for input2NArray
            const nsize_type & input2Sizes = input2NArray.sizes();
            const nstride_type & input2Strides = input2NArray.strides();
            nstride_type input2STND;
            nstride_type input2OTND;
            const value_type2 * input2Pointer = input2NArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = input1NArray.dimension();

            CalculateSTND(input1STND, input1Sizes, input1Strides);
            CalculateSTND(input2STND, input2Sizes, input2Strides);
            CalculateOTND(input2OTND, input2Strides, input2STND);
            InitializeTargets(input1Targets, input1Sizes, input1Strides, input1Pointer);

            OutputType incrementalResult = _incrementalOperationType::NeutralElement();

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                incrementalResult = _incrementalOperationType::Operate(
                    incrementalResult, _elementOperationType::Operate(*input1Pointer, *input2Pointer) );

                numberOfWrappedDimensions =
                    IncrementPointers(input1Targets, input1Pointer, input1Strides, input1STND);

                SyncCurrentPointer(input2Pointer, input2OTND, numberOfWrappedDimensions);
            }

            return incrementalResult;
        }   // Run method

    };  // SoNiNi class


    template <class _incrementalOperationType, class _elementOperationType>
    class SoNiSi
    {
    public:
        typedef typename _incrementalOperationType::OutputType OutputType;

        template <class _inputNArrayType, class _inputScalarType>
        static OutputType Run(const _inputNArrayType & inputNArray,
                              const _inputScalarType inputScalar)
        {
            typedef typename _inputNArrayType::value_type value_type;

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            vctFixedSizeVector<const value_type *, _dimension> inputTargets;
            const value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = inputNArray.dimension();

            CalculateSTND(inputSTND, inputSizes, inputStrides);
            InitializeTargets(inputTargets, inputSizes, inputStrides, inputPointer);

            OutputType incrementalResult = _incrementalOperationType::NeutralElement();

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                incrementalResult = _incrementalOperationType::Operate(
                    incrementalResult, _elementOperationType::Operate(*inputPointer, inputScalar) );

                numberOfWrappedDimensions =
                    IncrementPointers(inputTargets, inputPointer, inputStrides, inputSTND);
            }

            return incrementalResult;
        }   // Run method

    };  // SoNiSi class


    template <class _elementOperationType>
    class NoNiNi
    {
    public:
        template <class _outputNArrayType, class _input1NArrayType, class _input2NArrayType>
        static void Run(_outputNArrayType & outputNArray,
                        const _input1NArrayType & input1NArray,
                        const _input2NArrayType & input2NArray)
        {
            typedef typename _outputNArrayType::value_type output_value_type;
            typedef typename _input1NArrayType::value_type input1_value_type;
            typedef typename _input2NArrayType::value_type input2_value_type;

            // declare all variables used for outputNArray
            const nsize_type & outputSizes = outputNArray.sizes();
            const nstride_type & outputStrides = outputNArray.strides();
            nstride_type outputSTND;
            vctFixedSizeVector<const output_value_type *, _dimension> outputTargets;
            output_value_type * outputPointer = outputNArray.Pointer();

            // declare all variables used for input1NArray
            const nsize_type & input1Sizes = input1NArray.sizes();
            const nstride_type & input1Strides = input1NArray.strides();
            nstride_type input1STND;
            nstride_type input1OTND;
            const input1_value_type * input1Pointer = input1NArray.Pointer();

            // declare all variables used for input2NArray
            const nsize_type & input2Sizes = input2NArray.sizes();
            const nstride_type & input2Strides = input2NArray.strides();
            nstride_type input2STND;
            nstride_type input2OTND;
            const input2_value_type * input2Pointer = input2NArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = outputNArray.dimension();

            CalculateSTND(outputSTND, outputSizes, outputStrides);
            CalculateSTND(input1STND, input1Sizes, input1Strides);
            CalculateSTND(input2STND, input2Sizes, input2Strides);
            CalculateOTND(input1OTND, input1Strides, input1STND);
            CalculateOTND(input2OTND, input2Strides, input2STND);
            InitializeTargets(outputTargets, outputSizes, outputStrides, outputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                *outputPointer = _elementOperationType::Operate(*input1Pointer, *input2Pointer);

                numberOfWrappedDimensions =
                    IncrementPointers(outputTargets, outputPointer, outputStrides, outputSTND);

                SyncCurrentPointer(input1Pointer, input1OTND, numberOfWrappedDimensions);
                SyncCurrentPointer(input2Pointer, input2OTND, numberOfWrappedDimensions);
            }
        }   // Run method

    };  // NoNiNi class


    template <class _elementOperationType>
    class NoNiSi
    {
    public:
        template <class _outputNArrayType, class _inputNArrayType, class _inputScalarType>
        static void Run(_outputNArrayType & outputNArray,
                        const _inputNArrayType & inputNArray,
                        const _inputScalarType inputScalar)
        {
            typedef typename _outputNArrayType::value_type output_value_type;
            typedef typename _inputNArrayType::value_type input_value_type;

            // declare all variables used for outputNArray
            const nsize_type & outputSizes = outputNArray.sizes();
            const nstride_type & outputStrides = outputNArray.strides();
            nstride_type outputSTND;
            vctFixedSizeVector<const output_value_type *, _dimension> outputTargets;
            output_value_type * outputPointer = outputNArray.Pointer();

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            nstride_type inputOTND;
            const input_value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = outputNArray.dimension();

            CalculateSTND(outputSTND, outputSizes, outputStrides);
            CalculateSTND(inputSTND, inputSizes, inputStrides);
            CalculateOTND(inputOTND, inputStrides, inputSTND);
            InitializeTargets(outputTargets, outputSizes, outputStrides, outputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                *outputPointer = _elementOperationType::Operate(*inputPointer, inputScalar);

                numberOfWrappedDimensions =
                    IncrementPointers(outputTargets, outputPointer, outputStrides, outputSTND);

                SyncCurrentPointer(inputPointer, inputOTND, numberOfWrappedDimensions);
            }
        }   // Run method

    };  // NoNiSi class


    template <class _elementOperationType>
    class NoSiNi
    {
    public:
        template <class _outputNArrayType, class _inputScalarType, class _inputNArrayType>
        static void Run(_outputNArrayType & outputNArray,
                        const _inputScalarType inputScalar,
                        const _inputNArrayType & inputNArray)
        {
            typedef typename _outputNArrayType::value_type output_value_type;
            typedef typename _inputNArrayType::value_type input_value_type;

            // declare all variables used for outputNArray
            const nsize_type & outputSizes = outputNArray.sizes();
            const nstride_type & outputStrides = outputNArray.strides();
            nstride_type outputSTND;
            vctFixedSizeVector<const output_value_type *, _dimension> outputTargets;
            output_value_type * outputPointer = outputNArray.Pointer();

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            nstride_type inputOTND;
            const input_value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = outputNArray.dimension();

            CalculateSTND(outputSTND, outputSizes, outputStrides);
            CalculateSTND(inputSTND, inputSizes, inputStrides);
            CalculateOTND(inputOTND, inputStrides, inputSTND);
            InitializeTargets(outputTargets, outputSizes, outputStrides, outputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                *outputPointer = _elementOperationType::Operate(inputScalar, *inputPointer);

                numberOfWrappedDimensions =
                    IncrementPointers(outputTargets, outputPointer, outputStrides, outputSTND);

                SyncCurrentPointer(inputPointer, inputOTND, numberOfWrappedDimensions);
            }
        }   // Run method

    };  // NoSiNi class


    template <class _elementOperationType>
    class NioSi
    {
    public:
        template <class _inputOutputNArrayType, class _inputScalarType>
        static void Run(_inputOutputNArrayType & inputOutputNArray,
                        const _inputScalarType inputScalar)
        {
            typedef typename _inputOutputNArrayType::value_type value_type;

            const nsize_type & inputOutputSizes = inputOutputNArray.sizes();
            const nstride_type & inputOutputStrides = inputOutputNArray.strides();
            nstride_type inputOutputSTND;
            vctFixedSizeVector<const value_type *, _dimension> inputOutputTargets;
            value_type * inputOutputPointer = inputOutputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = inputOutputNArray.dimension();

            CalculateSTND(inputOutputSTND, inputOutputSizes, inputOutputStrides);
            InitializeTargets(inputOutputTargets, inputOutputSizes, inputOutputStrides, inputOutputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                _elementOperationType::Operate(*inputOutputPointer, inputScalar);

                numberOfWrappedDimensions =
                    IncrementPointers(inputOutputTargets, inputOutputPointer, inputOutputStrides, inputOutputSTND);
            }
        }   // Run method

    };  // NioSi class


    template <class _elementOperationType>
    class NioNi
    {
    public:
        template <class _inputOutputNArrayType, class _inputNArrayType>
        static void Run(_inputOutputNArrayType & inputOutputNArray,
                        const _inputNArrayType & inputNArray)
        {
            typedef typename _inputOutputNArrayType::value_type input_output_value_type;
            typedef typename _inputNArrayType::value_type input_value_type;

            // declare all variables used for inputOutputNArray
            const nsize_type & inputOutputSizes = inputOutputNArray.sizes();
            const nstride_type & inputOutputStrides = inputOutputNArray.strides();
            nstride_type inputOutputSTND;
            vctFixedSizeVector<const input_output_value_type *, _dimension> inputOutputTargets;
            input_output_value_type * inputOutputPointer = inputOutputNArray.Pointer();

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            nstride_type inputOTND;
            const input_value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = inputOutputNArray.dimension();

            CalculateSTND(inputOutputSTND, inputOutputSizes, inputOutputStrides);
            CalculateSTND(inputSTND, inputSizes, inputStrides);
            CalculateOTND(inputOTND, inputStrides, inputSTND);
            InitializeTargets(inputOutputTargets, inputOutputSizes, inputOutputStrides, inputOutputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                _elementOperationType::Operate(*inputOutputPointer, *inputPointer);

                numberOfWrappedDimensions =
                    IncrementPointers(inputOutputTargets, inputOutputPointer, inputOutputStrides, inputOutputSTND);

                SyncCurrentPointer(inputPointer, inputOTND, numberOfWrappedDimensions);
            }
        }   // Run method

    };  // NioNi class


    template <class _elementOperationType>
    class NoNi
    {
    public:
        template <class _outputNArrayType, class _inputNArrayType>
        static inline void Run(_outputNArrayType & outputNArray,
                               const _inputNArrayType & inputNArray)
        {
            typedef typename _outputNArrayType::value_type output_value_type;
            typedef typename _inputNArrayType::value_type input_value_type;

            // declare all variables used for outputNArray
            const nsize_type & outputSizes = outputNArray.sizes();
            const nstride_type & outputStrides = outputNArray.strides();
            nstride_type outputSTND;
            vctFixedSizeVector<const output_value_type *, _dimension> outputTargets;
            output_value_type * outputPointer = outputNArray.Pointer();

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            nstride_type inputOTND;
            const input_value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = outputNArray.dimension();

            CalculateSTND(outputSTND, outputSizes, outputStrides);
            CalculateSTND(inputSTND, inputSizes, inputStrides);
            CalculateOTND(inputOTND, inputStrides, inputSTND);
            InitializeTargets(outputTargets, outputSizes, outputStrides, outputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                *outputPointer = _elementOperationType::Operate(*inputPointer);

                numberOfWrappedDimensions =
                    IncrementPointers(outputTargets, outputPointer, outputStrides, outputSTND);

                SyncCurrentPointer(inputPointer, inputOTND, numberOfWrappedDimensions);
            }
        }   // Run method

    };  // NoNi class


    template <class _elementOperationType>
    class Nio
    {
    public:
        template <class _inputOutputNArrayType>
        static void Run(_inputOutputNArrayType & inputOutputNArray)
        {
            typedef typename _inputOutputNArrayType::value_type value_type;

            const nsize_type & inputOutputSizes = inputOutputNArray.sizes();
            const nstride_type & inputOutputStrides = inputOutputNArray.strides();
            nstride_type inputOutputSTND;
            vctFixedSizeVector<const value_type *, _dimension> inputOutputTargets;
            value_type * inputOutputPointer = inputOutputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = inputOutputNArray.dimension();

            CalculateSTND(inputOutputSTND, inputOutputSizes, inputOutputStrides);
            InitializeTargets(inputOutputTargets, inputOutputSizes, inputOutputStrides, inputOutputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                _elementOperationType::Operate(*inputOutputPointer);

                numberOfWrappedDimensions =
                    IncrementPointers(inputOutputTargets, inputOutputPointer, inputOutputStrides, inputOutputSTND);
            }
        }   // Run method

    };  // Nio class


    template <class _inputOutputElementOperationType, class _scalarNArrayElementOperationType>
    class NioSiNi
    {
    public:
        template <class _inputOutputNArrayType, class _inputScalarType, class _inputNArrayType>
        static void Run(_inputOutputNArrayType & inputOutputNArray,
                        const _inputScalarType inputScalar,
                        const _inputNArrayType & inputNArray)
        {
            typedef typename _inputOutputNArrayType::value_type input_output_value_type;
            typedef typename _inputNArrayType::value_type input_value_type;

            // declare all variables used for inputOutputNArray
            const nsize_type & inputOutputSizes = inputOutputNArray.sizes();
            const nstride_type & inputOutputStrides = inputOutputNArray.strides();
            nstride_type inputOutputSTND;
            vctFixedSizeVector<const input_output_value_type *, _dimension> inputOutputTargets;
            input_output_value_type * inputOutputPointer = inputOutputNArray.Pointer();

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            nstride_type inputOTND;
            const input_value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = inputOutputNArray.dimension();

            CalculateSTND(inputOutputSTND, inputOutputSizes, inputOutputStrides);
            CalculateSTND(inputSTND, inputSizes, inputStrides);
            CalculateOTND(inputOTND, inputStrides, inputSTND);
            InitializeTargets(inputOutputTargets, inputOutputSizes, inputOutputStrides, inputOutputPointer);

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                _inputOutputElementOperationType::Operate(
                    *inputOutputPointer, _scalarNArrayElementOperationType::Operate(inputScalar, *inputPointer) );

                numberOfWrappedDimensions =
                    IncrementPointers(inputOutputTargets, inputOutputPointer, inputOutputStrides, inputOutputSTND);

                SyncCurrentPointer(inputPointer, inputOTND, numberOfWrappedDimensions);
            }
        }   // Run method

    };  // NioSiNi class


    class MinAndMax
    {
    public:
        template <class _inputNArrayType>
        static void Run(const _inputNArrayType & inputNArray,
                        typename _inputNArrayType::value_type & minValue,
                        typename _inputNArrayType::value_type & maxValue)
        {
            typedef typename _inputNArrayType::value_type value_type;

            // declare all variables used for inputNArray
            const nsize_type & inputSizes = inputNArray.sizes();
            const nstride_type & inputStrides = inputNArray.strides();
            nstride_type inputSTND;
            vctFixedSizeVector<const value_type *, _dimension> inputTargets;
            const value_type * inputPointer = inputNArray.Pointer();

            dimension_type numberOfWrappedDimensions = 0;
            const dimension_type maxWrappedDimensions = inputNArray.dimension();

            CalculateSTND(inputSTND, inputSizes, inputStrides);
            InitializeTargets(inputTargets, inputSizes, inputStrides, inputPointer);

            value_type minElement, maxElement, inputElement;
            minElement = maxElement = *inputPointer;

            while (numberOfWrappedDimensions != maxWrappedDimensions)
            {
                inputElement = *inputPointer;

                if (inputElement < minElement)
                {
                    minElement = inputElement;
                }
                else if (inputElement > maxElement)
                {
                    maxElement = inputElement;
                }

                numberOfWrappedDimensions =
                    IncrementPointers(inputTargets, inputPointer, inputStrides, inputSTND);
            }

            minValue = minElement;
            maxValue = maxElement;
        }   // Run method

    };  // MinAndMax class
};  // vctDynamicNArrayLoopEngines


#endif  // _vctDynamicNArrayLoopEngines_h


// ****************************************************************************
//                              Change History
// ****************************************************************************
//
// $Log: vctDynamicNArrayLoopEngines.h,v $
// Revision 1.4  2007/04/26 19:33:57  anton
// All files in libraries: Applied new license text, separate copyright and
// updated dates, added standard header where missing.
//
// Revision 1.3  2007/04/02 01:25:54  anton
// vctDynamicNArrayLoopEngines.h: Commented unused variables in new engines
// implementation.
//
// Revision 1.2  2007/03/08 04:29:25  anton
// cisstVector NArray: Now templated by dimension (see #258).  Added template
// specialization for slices of dimension 1 NArray to return value_type.  Used
// Ofri's implementation of IncrementPointer() in loop engines (previous one did
// not handle more than dimension 3).  Also renamed DimSize, DimStride and
// Dimension to size(d), stride(d), and dimensioni (for STL flavor).
//
// Revision 1.1  2006/12/09 04:13:25  anton
// cisstVector: Imported code for NArray from Daniel Li with minor cleanup
// for license, dos2unix, a few renamed methods and port to gcc.
//
//
// ****************************************************************************