package yi.math;

/**
 * Class Matrix is to implement operations on two-dimensional matrices of type
 * double.
 * Copied from Matrix.java from Selim Aksoy
 * @author Selim Aksoy
 * @version 1.0 (Date: 3/26/2001)
 */
public class DoubleMatrix implements java.io.Serializable{

        public int height, width;
        public double[][] data;

        public DoubleMatrix( int matHeight, int matWidth ) {
                height = matHeight;
                width = matWidth;
                data = new double[ matHeight ][ matWidth ];
        }

        public DoubleMatrix( int matHeight, int matWidth, double[][] data) {
                height = matHeight;
                width = matWidth;
                this.data = data;
        }

        /**
         * add the elements of two matrices
         * @param inMatrix1
         * @param inMatrix2
         * @return
         */
        public static DoubleMatrix add( DoubleMatrix inMatrix1, DoubleMatrix inMatrix2 ) {

                int r, c;
                DoubleMatrix outMatrix = new DoubleMatrix( inMatrix1.height, inMatrix1.width );

                for( r = 0; r < inMatrix1.height; r++ )
                        for( c = 0; c < inMatrix1.width; c++ )
                                outMatrix.data[ r ][ c ] = inMatrix1.data[ r ][ c ] +
                                        inMatrix2.data[ r ][ c ];

                return outMatrix;

        }

        /**
         * subtract the elements of two matrices
         * @param inMatrix1
         * @param inMatrix2
         * @return
         */
        public static DoubleMatrix subtract( DoubleMatrix inMatrix1, DoubleMatrix inMatrix2 ) {

                int r, c;
                DoubleMatrix outMatrix = new DoubleMatrix( inMatrix1.height, inMatrix1.width );

                for( r = 0; r < inMatrix1.height; r++ )
                        for( c = 0; c < inMatrix1.width; c++ )
                                outMatrix.data[ r ][ c ] = inMatrix1.data[ r ][ c ] -
                                        inMatrix2.data[ r ][ c ];

                return outMatrix;

        }

        /**
         * multiply the elements of two matrices (element-wise multiplication)
         * @param inMatrix1
         * @param inMatrix2
         * @return
         */
        public static DoubleMatrix multiply( DoubleMatrix inMatrix1, DoubleMatrix inMatrix2 ) {

                int r, c;
                DoubleMatrix outMatrix = new DoubleMatrix( inMatrix1.height, inMatrix1.width );

                for( r = 0; r < inMatrix1.height; r++ )
                        for( c = 0; c < inMatrix1.width; c++ )
                                outMatrix.data[ r ][ c ] = inMatrix1.data[ r ][ c ] *
                                        inMatrix2.data[ r ][ c ];

                return outMatrix;

        }

        /**
         * copy the elements of one matrix into another
         * @param inMatrix
         * @param outRowSize
         * @param outColSize
         * @param targetHeight
         * @param targetWidth
         * @param startHeight
         * @param startWidth
         * @param endHeight
         * @param endWidth
         * @return
         */
        public static DoubleMatrix copy( DoubleMatrix inMatrix,
                int outRowSize, int outColSize, int targetHeight, int targetWidth,
                int startHeight, int startWidth, int endHeight, int endWidth ) {

                int r, c, h, w;
                DoubleMatrix outMatrix;
                outMatrix = new DoubleMatrix( outRowSize, outColSize );

                if( ( targetHeight >= 0 ) && ( targetHeight < outRowSize ) &&
                        ( targetWidth >= 0 ) && ( targetWidth < outColSize ) ) {
                        if( ( startHeight >= 0 ) && ( startHeight < inMatrix.height ) &&
                                ( startWidth >= 0 ) && ( startWidth < inMatrix.width ) ) {
                                h = endHeight - startHeight + 1;
                                w = endWidth - startWidth + 1;
                                if( ( h >= 1 ) && ( w >= 1 ) &&
                                        ( ( targetHeight + h - 1 ) < outRowSize ) &&
                                        ( ( targetWidth + w - 1 ) < outColSize ) ) {
                                        for( r = 0; r < h; r++ ) {
                                                for( c = 0; c < w; c++ ) {
                                                        outMatrix.data[ r + targetHeight ][
                                                                c + targetWidth ] =
                                                                inMatrix.data[ r + startHeight ][
                                                                c + startWidth ];
                                                }
                                        }
                                }
                        }
                } else {
                        System.out.println( "Matrix dimension error" );
                        System.exit( 1 );
                }

                return outMatrix;

        }

        /**
         * circularly shift the elements of a matrix by an offset
         * @param inMatrix
         * @param offset
         * @return
         */
        public static DoubleMatrix shift( DoubleMatrix inMatrix, int offset ) {

                int r, c;
                DoubleMatrix outMatrix = new DoubleMatrix( inMatrix.height, inMatrix.width );

                for( r = offset; r < inMatrix.height; r++ )
                        for( c = offset; c < inMatrix.width; c++ )
                                outMatrix.data[ r - offset ][ c - offset ] =
                                        inMatrix.data[ r ][ c ];

                for( r = offset; r < inMatrix.height; r++ )
                        for( c = 0; c < offset; c++ )
                                outMatrix.data[ r - offset ][ inMatrix.width + c - offset ] =
                                        inMatrix.data[ r ][ c ];

                for( r = 0; r < offset; r++ )
                        for( c = offset; c < inMatrix.width; c++ )
                                outMatrix.data[ inMatrix.height + r - offset ][ c - offset ] =
                                        inMatrix.data[ r ][ c ];

                for( r = 0; r < offset; r++ )
                        for( c = 0; c < offset; c++ )
                                outMatrix.data[ inMatrix.height + r - offset ][ inMatrix.width + c - offset ] =
                                        inMatrix.data[ r ][ c ];

                return outMatrix;

        }

        /**
         * 2-D Fast Fourier Transform (FFT)
         * @param Output_real
         * @param Output_imag
         * @param Input_real
         * @param Input_imag
         */
        public static void fft2( DoubleMatrix Output_real, DoubleMatrix Output_imag,
                DoubleMatrix Input_real, DoubleMatrix Input_imag ) {

                fftFour2(Output_real.data, Output_imag.data, Input_real.data,
                        Input_imag.data, Input_real.height, Input_real.width, 1);

        }

        /**
         * 2-D Inverse Fast Fourier Transform (IFFT)
         * @param Output_real
         * @param Output_imag
         * @param Input_real
         * @param Input_imag
         */
        public static void ifft2( DoubleMatrix Output_real, DoubleMatrix Output_imag,
                DoubleMatrix Input_real, DoubleMatrix Input_imag ) {

                int i, j;
                double NN;

                fftFour2(Output_real.data, Output_imag.data, Input_real.data,
                        Input_imag.data, Input_real.height, Input_real.width, -1);

                NN = (double) (Input_real.height*Input_real.width);

                for (i=0;i<Input_real.height;i++)
                        for (j=0;j<Input_real.width;j++) {
                                Output_real.data[i][j] /= NN;
                                Output_imag.data[i][j] /= NN;
                        }

        }

        /**
         * Methods below were directly copied from Ma's C code. They
         * implement the "butterfly" algorithms of FFT (were actually
         * copied from the book called Numerical Recipes).
         *
         * 2-D fourier transform of data with real part stored in
         * "rdata" and imaginary part in "idata" with size "rs" x
         * "cs". The result is in "fftr" and "ffti". The isign is
         * "isign" =  1 forward, and "isign" = -1 inverse
         *
         * @param fftr
         * @param ffti
         * @param rdata
         * @param idata
         * @param rs
         * @param cs
         * @param isign
         */
        private static void fftFour2( double[][] fftr, double[][] ffti,
                double[][] rdata, double[][] idata, int rs, int cs, int isign) {

                double[][] T = new double[ 2 * rs ][ cs ];
                double[] tmp1 = new double[ 2 * cs ];
                double[] tmp2 = new double[ 2 * rs ];
                int i, j;

                for (i=0;i<rs;i++) {
                        for (j=0;j<cs;j++) {
                                tmp1[j*2] = rdata[i][j];
                                tmp1[j*2+1] = idata[i][j];
                        }
                        fftFour1(tmp1, cs, isign);
                        for (j=0;j<cs;j++) {
                                T[i*2][j] = tmp1[j*2];
                                T[i*2+1][j] = tmp1[j*2+1];
                        }
                }

                for (i=0;i<cs;i++) {
                        for (j=0;j<rs;j++) {
                                tmp2[j*2] = T[j*2][i];
                                tmp2[j*2+1] = T[j*2+1][i];
                        }
                        fftFour1(tmp2,rs,isign);
                        for (j=0;j<rs;j++) {
                                fftr[j][i] = tmp2[j*2];
                                ffti[j][i] = tmp2[j*2+1];
                        }
                }

        }

        private static void fftFour1( double[] data, int nn, int isign ) {

                int n, mmax, m, j, istep, i;
                double wtemp, wr, wpr, wpi, wi, theta;
                double tempr, tempi;

                n = nn << 1;
                j = 1;
                for (i=1;i<n;i+=2) {
                        if (j > i) {
                                tempr=data[j-1];data[j-1]=data[i-1];data[i-1]=tempr;
                                tempr=data[j];data[j]=data[i];data[i]=tempr;
                        }
                        m = n >> 1;
                        while (m >= 2 && j > m) {
                                j -= m;
                                m >>= 1;
                        }
                        j += m;
                }
                mmax = 2;
                while (n > mmax) {
                        istep = 2*mmax;
                        theta = 6.28318530717959/(isign*mmax);
                        wtemp = Math.sin(0.5*theta);
                        wpr = -2.0*wtemp*wtemp;
                        wpi = Math.sin(theta);
                        wr = 1.0;
                        wi = 0.0;
                        for (m=1;m<mmax;m+=2) {
                                for (i=m;i<=n;i+=istep) {
                                        j = i+mmax;
                                        tempr = wr*data[j-1]-wi*data[j];
                                        tempi = wr*data[j]+wi*data[j-1];
                                        data[j-1] = data[i-1]-tempr;
                                        data[j] = data[i]-tempi;
                                        data[i-1] += tempr;
                                        data[i] += tempi;
                                }
                                wr = (wtemp=wr)*wpr-wi*wpi+wr;
                                wi = wi*wpr+wtemp*wpi+wi;
                        }
                        mmax = istep;
                }

        }

}