#include "mex.h"
#include "matrix.h"
#include "math.h"
#include "stdlib.h"

/*
 *   Copyright (C) 2009 Cesare Magri
 *   Version: 1.0.2
 */

/*
 *-------
 * LICENSE
 * -------
 * This software is distributed free under the condition that:
 *
 * 1. it shall not be incorporated in software that is subsequently sold;
 *
 * 2. the authorship of the software shall be acknowledged and the following
 *    article shall be properly cited in any publication that uses results
 *    generated by the software:
 *
 *      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
 *      toolbox for the fast information analysis of multiple-site LFP, EEG
 *      and spike train recordings. BMC Neuroscience 2009 10(1):81;
 *
 * 3.  this notice shall remain in place in each source file.
 *
 */

/* ----------
 * DISCLAIMER
 * ----------
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
    /* Input:                                                              */
    double *R;
    double *Nt;
    
    /* Dimensions of R:                                                    */
    mwSize Nc, maxNt, Ns;
    
    /* Indexes:                                                            */
    mwIndex c, s;
    mwIndex i, j;
    mwIndex currentIndx;
    mwIndex randomIndx;
    
    /* Output:                                                             */
    double *Rsh;
    
    /*
     * --------------------------------------------------------------------
     * The routine is based on the following simple algorithm for permuting
     * an array (i.e., any ordering of the array is equally probable or any
     * element has equal chance of being in any position):
     *
     *      For i=1 to n-1:
     *      - Let j=r(i+1)
     *      - Swap a[i] and a[j]
     *
     * where r[n] is a rand number generated between 0 and n-1.
     * --------------------------------------------------------------------
     */
        
    /* Reading input:                                                      */
    R  = mxGetPr(prhs[0]);
    Nt = mxGetPr(prhs[1]);
    
    /* Reading dimensions of R:                                            */
    Nc    = mxGetDimensions(prhs[0])[0];
    maxNt = mxGetDimensions(prhs[0])[1];
    /* If the users inputs a response matrix with a single stimulus then   */
    /* the number of dimensions for the response matrix will be 2 and      */
    /* asking for the number of elements in the third dimension will leak  */
    /* into a meaningless value. Thus first we need to check that there is */
    /* indeed a third dimension. We do this using MXGETNUMBEROFDIMENSIONS  */
    if(mxGetNumberOfDimensions(prhs[0])>2)
        Ns = mxGetDimensions(prhs[0])[2];
    else
        Ns = 1;
    
    /* Allocating output:                                                  */
    plhs[0] = mxCreateNumericArray(3, mxGetDimensions(prhs[0]), mxDOUBLE_CLASS, 0);
    Rsh = mxGetPr(plhs[0]);
    
    for(s=0;s<Ns;s++)
    {      
        /* For i = 1 to Nt[s]-1                                            */
        for(i=0;i<Nt[s];i++)
        {
            /* Generating random number between 0 and i                    */
            j = rand() % (i+1);

            for(c=0;c<Nc;c++)
            {
                currentIndx = c + i*Nc + s*Nc*maxNt;
                randomIndx  = c + j*Nc + s*Nc*maxNt;

                Rsh[currentIndx] = Rsh[randomIndx];
                Rsh[randomIndx]  = R[currentIndx];
            }
        }
    }
}