using Distributions
using NumericExtensions
using Devectorize
#using Permutations
#using GramSchmidt
#using IterativeSolvers
using Iterators
using Gadfly
using DataFrames

function F_func(X,Y)
		n1 = size(X,1)
		n2 = size(Y,1)
		p = size(X,2)
		mu1 = zeros(p)
		mu2 = zeros(p)
		f = 0.0
		mdiff = zeros(p)
		mean!(mu1,X,1)
		mean!(mu2,Y,1)
		# for i=1:p
		#    mdiff[i] = mu2[i] - mu1[i]
	 #    end
	    @devec mdiff = mu2 - mu1
		#@show mdiff
		spool = zeros(p,p)
		n0 = n1 + n2 - 2
		ndel  = 1/n1 + 1/n2
		spool = ((n1-1).*cov(X) + (n2-1).*cov(Y))/n0 ### covariance matrix
		#@show spool
		f =  ((n0 - p + 1)/(n0*p*ndel)).*mdiff'*inv(spool)*mdiff
		return f[1,1] #, ((n0 - p + 1)/(n0*p*ndel)), mdiff'*inv(spool)*mdiff
end
###################################
function BayesFac_Block_noloc(pam) 
	n1 = ifloor(pam[1])
	n2 = ifloor(pam[2])
	 p = ifloor(pam[3])
	 k = ifloor(pam[4])
	 m = ifloor(pam[5])
	kap0 = pam[6]
	pct = pam[7]
	mu = pam[8]
	Nsim = 100  ### number of data sets to simulate
    N = 2000      ### number of random projections
	p0 = ifloor(p*pct) ## percent of elements that are null
	mu0 = zeros(p)
	lb0 = 0.0
	cst0 = 0.0
	
	out = zeros(N)
	indx = zeros(p0)
	
	f = 0.0
    R = zeros(p,m) 
    g1 = zeros(n1,m)
    g2 = zeros(n2,m)
	g10 = zeros(n1,m)
    g20 = zeros(n2,m)
    C01 = ones(m,m)
    C02 = ones(m,m)
    C0 = ones(m,m)
    dmut = ones(m)
    
    m1 = zeros(m)
    m2 = zeros(m)
    res = zeros(Nsim,9)
	
    sigp = zeros(p,p)
     B = .85.*eye(25) +.15.*ones(25,25)
     num_b = ifloor(maximum(p)./25)
   for i in 1:num_b
     idx=(25.*(i-1)+1):(i.*25)
     sigp[idx,idx] = B
   end
    Y1 = zeros(p,n1)
    Y2 = zeros(p,n2)
    d1 = MvNormal(sigp)
    d2 = MvNormal(mu0,sigp)
	d0 = MvNormal(eye(p))
	
	## define constant 
	n_del0 = n1*n2/(n1+n2)
	n_del1 = 1/(1/n_del0 + 1/kap0)
	t0 = n_del0/kap0
	n = n1 + n2
	#c0 = n-m-1
	c1 = m/((1+t0)*(n-m-1))
	c3 = m/(n-m-1)
	c4 = .5*(n-1)
	c5 = -.5*m*log(1+t0)
	c6 = (n_del0 + kap0)/(n*n_del0*n_del0)
	#d0 = MvNormal(mu00,sigp)
	temp0=0.0
    for j in 1:Nsim                          ## simulate Nsim data sets
	 sample!([1:p],indx,replace=false)
	 rand!(d1,mu0)  ### Stores the true mean vector use sample in StatsBase package
	 mu0[indx] = 0
	  rand!(d1,Y1)
	  rand!(d2,Y2)
	
	    ## Standardization step
     scale!(vec(sqrt(var(Y1,2))),Y1) 
     scale!(vec(sqrt(var(Y2,2))),Y2) 
     
	   ######################################### End of initializations
	   for i in 1:N
	    rand!(d0,R) 
	    At_mul_B!(g1,Y1,R)
	    At_mul_B!(g2,Y2,R)
		f = F_func(g1,g2)
	    mean!(m1,g1,1)
	    mean!(m2,g2,1)
		broadcast!(-,g10,g1,mean(g1,1)) ## return g1 after centering the columns and changes g1 
        broadcast!(-,g20,g2,mean(g2,1)) ## return g2 after centering the columns and changes g2
		At_mul_B!(C01,g10,g10) ## Compute the g1'*g1 and stores it in C01
		At_mul_B!(C02,g20,g20)

	     @devec C0 = C01 + C02
	     @devec dmut = m2 - m1
	    
		lb0 =  k==0 ? 0.0 : c6.*dmut'*inv(C0 + n_del1.*At_mul_B(dmut',dmut'))*dmut
		out[i] = c5 - c4*(log(1 + c1*f) - log(1 + c3*f)) - log(1+t0) + log(lb0[1,1]+m) - log(m)
	   end
	  res[j,1] = n1*1.0 
	  res[j,2] = n2*1.0
	  res[j,3] = p*1.0
	  res[j,4] = k*1.0
	  res[j,5] = m*1.0
	  res[j,6] = kap0
	  res[j,7] = pct
	  res[j,8] = mu
	  res[j,9] = mean(out)
    end
     println("Done!")	
     return [mean(res,1)',mean(res[:,9] .> 1),mean(res[:,9] .> 3),mean(res[:,9] .> 5)]' ## averages over data sets 
end