#!/usr/bin/env python
# -*- coding: utf-8 -*-
# 
# Copyright 2016 <+YOU OR YOUR COMPANY+>.
# 
# This 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, or (at your option)
# any later version.
# 
# This software 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 software; see the file COPYING.  If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street,
# Boston, MA 02110-1301, USA.
# 

import numpy
import scipy.linalg
import numpy.linalg
import math
from TracyWidom import * 
from gnuradio import gr


class MME(gr.sync_block):
    """
    docstring for block MME
    """
    def __init__(self,samples,slots,false,algo):
	self.samples = samples
	self.slots = slots
	self.false = false
	self.algo = algo
        gr.sync_block.__init__(self,
            name="MME",
            in_sig=[numpy.float32],
            out_sig=[numpy.float32])
    '''
    def R(self, x):
	x0 = x - numpy.mean(x)
	L = self.slots
	Ns = len(x0)
	lbd = numpy.empty(L)
	for l in xrange(L):
		if l > 0:
			xu = x0[:-l]
		else:
			xu = x0
			lbd[l] = numpy.dot(xu, x0[l:])/(Ns-l)

	return scipy.linalg.toeplitz(lbd)
    '''
    #compute the eigenvalues
    def lbd(self, x):
	#R = self.R(x)
        R = numpy.cov(x)      #compute the received sample covariance matrix
	lbd = numpy.linalg.eigvalsh(R)
	return numpy.abs(lbd)

    #compute the threshold
    def thr_MME(self,Ns,M,false):
	x = math.sqrt(Ns)+math.sqrt(M*Ns)
	y = math.sqrt(Ns)-math.sqrt(M*Ns)
	f = TracyWidom().cdfinv(1-false)
	thr1 = x*x/(y*y)*(1+x**(-2/3)/(M*Ns*Ns)**(1/6)*f)
        return thr1
    
    def work(self, input_items, output_items):
        in0 = input_items[0]
        out = output_items[0]
        Ns = self.samples      #number of samples
        M = self.slots         #number of slots
	false = self.false     #probability of false alarm
        # <+signal processing here+>
        n = len(in0)/Ns/M
        in1 = numpy.array(in0)
        in1 = in1[0:n*Ns*M:1]  #slice the array
        in1 = numpy.reshape(in1,(n*M,Ns)) # >>> x = np.arange(16.0).reshape(4, 4)
        in2 = numpy.split(in1,n)
        if(self.algo == 0):    #use MME detection
	    for i in range(n):
	        lbd = self.lbd(in2[i])
	        lbd.sort()
	        out[i] = lbd[-1]/lbd[0]
                if (out[i] >= self.thr_MME(Ns,M,false)):
		    out[i] = 1
                else: 
		    out[i] = 0
		    
        if(self.algo == 1):   #use EME detection
	    for i in range(n):
	        lbd = self.lbd(in2[i])
	        lbd.sort()
	        out[i] = numpy.sum(lbd**2)/lbd[0]
      
        return len(output_items[0])