def V_Multi(m, n):
    x,y,z=(m[0],m[1],m[2])
    a,b,c=(n[0],n[1],n[2])
    
    p,q,r = (a-x,b-y,c-z)
    
    return 90-degrees(acos( (x*p+y*q+z*r)/ (sqrt(x**2+y**2+z**2)*sqrt(p**2+q**2+r**2))))

def get_stderr_fit(f,Xdata,popt,pcov):
    
   Y= f(Xdata, popt)
   listdY=[]
   for i in xrange(len(popt)):
       p=popt[i]
       dp= abs(p)/1e6+1e-20
       popt[i]+=dp
       Yi= f(Xdata, popt)
       dY= (Yi-Y)/dp
       listdY.append(dY)
       popt[i]-=dp
   listdY= scipy.array(listdY)
   #list dy is the d in the derivation. it has N X M
   #pcov is N X N
   
   left= scipy.dot(listdY.T,pcov)
   right=scipy.dot(left,listdY)
   
   sigma2y= right.diagonal()
   #sigma2y is a standard function of fit
   mean_sigma2y= scipy.mean(right.diagonal())
   
   M= Xdata.shape[0]
   N= len(popt)
   avg_stddev_data=scipy.sqrt(M*mean_sigma2y/N)
   sigmay= scipy.sqrt(sigma2y)
   return sigmay,avg_stddev_data

def Rz_to_coshucosv(R,z,delta=1.):
    """
    NAME:

       Rz_to_coshucosv

    PURPOSE:

       calculate prolate confocal cosh(u) and cos(v) coordinates from R,z, and delta

    INPUT:

       R - radius

       z - height

       delta= focus

    OUTPUT:

       (cosh(u),cos(v))

    HISTORY:

       2012-11-27 - Written - Bovy (IAS)

    """
    d12= (z+delta)**2.+R**2.
    d22= (z-delta)**2.+R**2.
    coshu= 0.5/delta*(sc.sqrt(d12)+sc.sqrt(d22))
    cosv=  0.5/delta*(sc.sqrt(d12)-sc.sqrt(d22))
    return (coshu,cosv)

    def survival_function(loss_ratio, **kwargs):
        """
            Static method that prepares the calculation parameters
            to be passed to stats.lognorm.sf

            :param loss_ratio: current loss ratio
            :type loss_ratio: float

            :param kwargs: convenience dictionary
            :type kwargs: :py:class:`dict` with the following
                keys:
                    **vf** - vulnerability function as provided by
                            :py:class:`openquake.shapes.VulnerabilityFunction`
                    **col** - matrix column number
        """
        vuln_function = kwargs.get('vf')
        position = kwargs.get('col')

        vf_loss_ratio = vuln_function.loss_ratios[position]

        stddev = vuln_function.covs[position] * vf_loss_ratio

        variance = stddev ** 2.0

        sigma = sqrt(log((variance / vf_loss_ratio ** 2.0) + 1.0))
        mu = exp(log(vf_loss_ratio ** 2.0 /
            sqrt(variance + vf_loss_ratio ** 2.0)))

        return stats.lognorm.sf(loss_ratio, sigma, scale=mu)

    def specular_incidence(self, pol = 'TE'):
        """ Return a vector of plane wave amplitudes corresponding
            to specular incidence in the specified polarisation.

            i.e. all elements are 0 except the zeroth order.
        """
        # Element corresponding to 0th order, TE
        spec_TE = self.specular_order
        # Element corresponding to 0th order, TM
        spec_TM = self.specular_order + self.structure.num_pw_per_pol
        tot_num_pw = self.structure.num_pw_per_pol * 2

        inc_amp = np.mat(np.zeros(tot_num_pw, dtype='complex128')).T
        if   'TE' == pol:
            inc_amp[spec_TE] = 1
        elif 'TM' == pol:
            inc_amp[spec_TM] = 1
        elif 'R Circ' == pol:
            inc_amp[spec_TE] = 1/sqrt(2.)
            inc_amp[spec_TM] = +1j/sqrt(2.)
        elif 'L Circ' == pol:
            inc_amp[spec_TE] = 1/sqrt(2.)
            inc_amp[spec_TM] = -1j/sqrt(2.)
        else:
            raise NotImplementedError, \
            "Must select from the currently implemented polarisations; \
             TE, TM, R Circ, L Circ."

        return inc_amp

 def _getScalesDiag(self,termx=0):
     """
     Internal function for parameter initialization
     Uses 2 term single trait model to get covar params for initialization
     
     Args:
         termx:      non-noise term terms that is used for initialization 
     """
     assert self.P>1, 'VarianceDecomposition:: diagonal init_method allowed only for multi trait models' 
     assert self.noisPos!=None, 'VarianceDecomposition:: noise term has to be set'
     assert termx<self.n_randEffs-1, 'VarianceDecomposition:: termx>=n_randEffs-1'
     assert self.trait_covar_type[self.noisPos] not in ['lowrank','block','fixed'], 'VarianceDecomposition:: diagonal initializaiton not posible for such a parametrization'
     assert self.trait_covar_type[termx] not in ['lowrank','block','fixed'], 'VarianceDecimposition:: diagonal initializaiton not posible for such a parametrization'
     scales = []
     res = self._getH2singleTrait(self.vd.getTerm(termx).getK())
     scaleg = sp.sqrt(res['varg'].mean())
     scalen = sp.sqrt(res['varn'].mean())
     for term_i in range(self.n_randEffs):
         if term_i==termx:
             _scales = scaleg*self.diag[term_i]
         elif term_i==self.noisPos:
             _scales = scalen*self.diag[term_i]
         else:
             _scales = 0.*self.diag[term_i]
         if self.jitter[term_i]>0:
             _scales = sp.concatenate((_scales,sp.array([sp.sqrt(self.jitter[term_i])])))
         scales.append(_scales)
     return sp.concatenate(scales)

def test_periodogram_csd():
    """Test corner cases of  periodogram_csd"""

    arsig1, _, _ = utils.ar_generator(N=1024)
    arsig2, _, _ = utils.ar_generator(N=1024)

    tseries = np.vstack([arsig1, arsig2])
    
    Sk = np.fft.fft(tseries)

    f1, c1 = tsa.periodogram_csd(tseries)
    f2, c2 = tsa.periodogram_csd(tseries, Sk=Sk)
    npt.assert_equal(c1, c2)

    # Check that providing a complex signal does the right thing
    # (i.e. two-sided spectrum): 
    N = 1024 
    r, _, _ = utils.ar_generator(N=N)
    c, _, _ = utils.ar_generator(N=N)
    arsig1 = r + c * scipy.sqrt(-1)

    r, _, _ = utils.ar_generator(N=N)
    c, _, _ = utils.ar_generator(N=N)
    arsig2 = r + c * scipy.sqrt(-1)

    tseries = np.vstack([arsig1, arsig2])

    f, c = tsa.periodogram_csd(tseries)
    npt.assert_equal(f.shape[0], N) # Should be N, not the one-sided N/2 + 1

def drazin(A, tol):
    CB = A.copy()

    Bs = []
    Cs = []
    k = 1

    while not (sp.absolute(CB) < tol).all() and sp.absolute(la.det(CB)) < tol:
        U, s, Vh = la.svd(CB)
        S = sp.diag(s)
        S = S * (S > tol)
        r = sp.count_nonzero(S)
        B = sp.dot(U, sp.sqrt(S))
        C = sp.dot(sp.sqrt(S), Vh)
        B = B[:, 0:r]
        Bs.append(B)
        C = C[0:r, :]
        Cs.append(C)
        CB = sp.dot(C, B)
        k += 1

    D = sp.eye(A.shape[0])
    for B in Bs:
        D = sp.dot(D, B)
    if (sp.absolute(CB) < tol).all():
        D = sp.dot(D, CB)
    else:
        D = sp.dot(D, np.linalg.matrix_power(CB, -(k + 1)))
    for C in reversed(Cs):
        D = sp.dot(D, C)
    return D

def ndot_product(features1,
                 features2 = None):

    """
    generates kernel based on normalized dot product

    input :
    features1 : vectors representing the rows in the matrix
    features2 : vectors representing the columns in the matrix

    output :
    out : similarity matrix

    """

    features1.shape = features1.shape[0], -1
    features1 = features1/sp.sqrt((features1**2.).sum(1))[:, None]

    if features2 is None:
        features2 = features1
    else:
        features2.shape = features2.shape[0], -1
        features2 = features2/sp.sqrt((features2**2.).sum(1))[:, None]

    out = sp.dot(features1, features2.T)

    return out

    def _genBgTerm_fromXX(self,vTot,vCommon,XX,a=None,c=None):
        """
        generate background term from SNPs

        Args:
            vTot: variance of Yc+Yi
            vCommon: variance of Yc
            XX: kinship matrix
            a: common scales, it can be set for debugging purposes
            c: indipendent scales, it can be set for debugging purposes
        """
        vSpecific = vTot-vCommon

        SP.random.seed(0)
        if c==None: c = SP.randn(self.P)
        XX += 1e-3 * SP.eye(XX.shape[0])
        L = LA.cholesky(XX,lower=True)

        # common effect
        R = self.genWeights(self.N,self.P)
        A = self.genTraitEffect()
        if a is not None: A[0,:] = a
        Yc = SP.dot(L,SP.dot(R,A))
        Yc*= SP.sqrt(vCommon)/SP.sqrt(Yc.var(0).mean())

        # specific effect
        R = SP.randn(self.N,self.P)
        Yi = SP.dot(L,SP.dot(R,SP.diag(c)))
        Yi*= SP.sqrt(vSpecific)/SP.sqrt(Yi.var(0).mean())

        return Yc, Yi

def _sampling_matrix(hessian, cutoff=0, temperature=1, step_scale=1):    
    # basically need SVD of hessian - singular values and eigenvectors
    # hessian = u * diag(singVals) * vh
    u, sing_vals, vh = scipy.linalg.svd(0.5 * hessian)

    # scroll through the singular values and find the ones whose inverses will
    # be huge and set them to zero also, load up the array of singular values 
    # that we store
    # cutoff = (1.0/_.singVals[0])*1.0e03
    # double cutoff = _.singVals[0]*1.0e-02
    cutoff_sing_val = cutoff * max(sing_vals)

    D = 1.0/scipy.maximum(sing_vals, cutoff_sing_val)

    ## now fill in the sampling matrix ("square root" of the Hessian)
    ## note that sqrt(D[i]) is taken here whereas Kevin took sqrt(D[j])
    ## this is because vh is the transpose of his PT -JJW
    samp_mat = scipy.transpose(vh) * scipy.sqrt(D)

    # Divide the sampling matrix by an additional factor such
    # that the expected quadratic increase in cost will be about 1.
    cutoff_vals = scipy.compress(sing_vals < cutoff_sing_val, sing_vals)
    if len(cutoff_vals):
        scale = scipy.sqrt(len(sing_vals) - len(cutoff_vals)
                           + sum(cutoff_vals)/cutoff_sing_val)
    else:
        scale = scipy.sqrt(len(sing_vals))

    samp_mat /= scale
    samp_mat *= step_scale
    samp_mat *= scipy.sqrt(temperature)

    return samp_mat

 def test_gets_thermal_with_correlated(self):
     """Checks that the part of the freq_modes code that compensates the
     thermal for mode subtraction works."""
     self.data *= sp.sqrt(self.bw * 2)  # Makes thermal unity.
     # Need to add something correlated so the modes arn't just the
     # channels.
     correlated_overf = noise_power.generate_overf_noise(1, 
                                 -2, 0.5, self.dt, self.data.shape[0])
     correlated_overf += (rand.normal(size=(self.data.shape[0],))
                          * sp.sqrt((self.bw * 2) * 0.3))
     self.data += correlated_overf[:,None,None,None] / sp.sqrt(self.nf)
     Blocks = self.make_blocks()
     # Mask a channel out completly.
     for Data in Blocks:
         Data.data[:,:,:,3] = ma.masked
     model = 'freq_modes_over_f_4'  # Take out 20% of the thermal power.
     parameters = mn.measure_noise_parameters(Blocks,
                                             [model])
     right_ans = sp.ones(self.nf)
     right_ans[3] = T_infinity**2
     for p in parameters.itervalues():
         pars = p[model]
         thermal = pars['thermal']
         self.assertTrue(sp.allclose(thermal,
                                     right_ans, rtol=0.3))
         mean_thermal = sp.mean(thermal[right_ans==1])
         self.assertTrue(sp.allclose(mean_thermal, 1, rtol=0.05))
         self.assertTrue(sp.allclose(pars['over_f_mode_0']['thermal'],
                                     0.3, atol=0.1))

def ZYFF(Te, EIJ):
    """Computes `ZY` and `FF`, used in other functions.
    
    If `EIJ` is a scalar, the output has the same shape as `Te`. If `EIJ` is an
    array, the output has shape `EIJ.shape` + `Te.shape`. This should keep the
    output broadcastable with `Te`.
    
    Parameters
    ----------
    Te : array of float
        Electron temperature. Shape is arbitrary.
    EIJ : scalar float or array of float
        Energy difference.
    """
    # Expand the dimensions of EIJ to produce the desired output shape:
    Te = scipy.asarray(Te, dtype=float)
    EIJ = scipy.asarray(EIJ, dtype=float)
    for n in xrange(Te.ndim):
        EIJ = scipy.expand_dims(EIJ, axis=-1)
    
    ZY = EIJ / (1e3 * Te)
    
    FF = scipy.zeros_like(ZY)
    mask = (ZY >= 1.5)
    FF[mask] = scipy.log((ZY[mask] + 1) / ZY[mask]) - (0.36 + 0.03 * scipy.sqrt(ZY[mask] + 0.01)) / (ZY[mask] + 1)**2
    mask = ~mask
    FF[mask] = scipy.log((ZY[mask] + 1) / ZY[mask]) - (0.36 + 0.03 / scipy.sqrt(ZY[mask] + 0.01)) / (ZY[mask] + 1)**2
    
    return ZY, FF

def heart(scale,ndim,time):
    
    
    percent = 1.05 + 0.5*np.random.rand()
    
    real_heart = int(time*percent)
    
    ratio = 0.5*np.random.rand()  #xy ratio
    x = scipy.linspace(-2,2,real_heart/2)
    y1 = scipy.sqrt(1-(abs(x)-1)**2)
    y2 = -3*scipy.sqrt(1-(abs(x[::-1])/2)**0.5)
    
    Y = np.concatenate([y1,y2])
    X  = ratio*np.concatenate([x,x[::-1]])
    
    shift = np.random.randint(0,real_heart)
    Y = np.roll(Y,shift)[:time]
    X = np.roll(X,shift)[:time]
    traj = np.array([X,Y,np.zeros_like(Y)]).T
    
    alpha = 2*3.14 * np.random.rand()
    
    if ndim == 2:
        traj = traj[::,:2]
                    
   
    
    noise = diffusive(scale/800.,ndim,time+1,epsilon=1e-7)
    return scale*random_rot(traj,alpha,ndim) + noise[:-1]-noise[1:]

def get_stderr_fit(f,Xdata,popt,pcov):
	Y=f(Xdata,popt)
	listdY=[]
	for i in xrange(len(popt)):
		p=popt[i]
		dp=abs(p)/1e6+1e-20
		popt[i]+=dp
		Yi=f(Xdata,popt)
		dY=(Yi-Y)/dp
		listdY.append(dY)
		popt[i]-=dp
	listdY=scipy.array(listdY)
	#listdY is an array with N rows and M columns, N=len(popt), M=len(xdata[0])
	#pcov is an array with N rows and N columns
	left=scipy.dot(listdY.T,pcov) 
	#left is an array of M rows and N columns
	right=scipy.dot(left,listdY)
	#right is an array of M rows and M columns
	sigma2y=right.diagonal()
	#sigma2y is standard error of fit and function  of X
	mean_sigma2y=scipy.mean(right.diagonal())
	M=Xdata.shape[1];print M
	N=len(popt);print N
	avg_stddev_data=scipy.sqrt(M*mean_sigma2y/N)
	#this is because if exp error is constant at sig_dat,then mean_sigma2y=N/M*sig_dat**2
	sigmay=scipy.sqrt(sigma2y)
	return sigmay,avg_stddev_data

def generate_test_data_w_sum_stats(h2=0.5, n=100000, n_sample=100, m=50000, model='gaussian', 
                                         p=1.0, conseq_r2=0, m_ld_chunk_size=100):
    """
    Generate 
    """
    #Get LD sample matrix
    D_sample = genotypes.get_sample_D(200,conseq_r2=conseq_r2,m=m_ld_chunk_size)
    
    #Simulate beta_hats
    ret_dict = simulate_beta_hats(h2=h2, n=n, n_sample=n_sample, m=m, model=model, p=p, 
                                    conseq_r2=conseq_r2, m_ld_chunk_size=m_ld_chunk_size, D_sample=D_sample)
    
    #Simulate test genotypes
    test_snps = genotypes.simulate_genotypes_w_ld(n_sample=n_sample, m=m, conseq_r2=conseq_r2, 
                                                  m_ld_chunk_size=m_ld_chunk_size)
    ret_dict['test_snps'] = test_snps
    
    #Simulate test phenotypes
    phen_noise = stats.norm.rvs(0, sp.sqrt(1.0 - h2), size=n_sample) 
    phen_noise = sp.sqrt((1.0 - h2) / sp.var(phen_noise)) * phen_noise
    genetic_part = sp.dot(test_snps.T, ret_dict['betas'])
    genetic_part = sp.sqrt(h2 / sp.var(genetic_part)) * genetic_part
    test_phen = genetic_part + phen_noise
    ret_dict['test_phen'] = test_phen
    return ret_dict

def simulate_genotypes_w_ld(n, m, n_samples, m_ld_chunk_size=100, r2=0.9, validation = False):
	val_set_gen = []
	val_set_D = []
	for i in xrange(n_samples):
		snps = sp.zeros((m, n), dtype='float32')
		num_chunks = m / m_ld_chunk_size
		for chunk in xrange(num_chunks):
			X = sp.zeros((m_ld_chunk_size, n), dtype='float32')
			X[0] = stats.norm.rvs(size=n)
			for j in xrange(1, m_ld_chunk_size):
				X[j] = sp.sqrt(r2) * X[j - 1] + sp.sqrt(1 - r2) * stats.norm.rvs(size=n)
			start_i = chunk * m_ld_chunk_size
			stop_i = start_i + m_ld_chunk_size
			snps[start_i:stop_i] = X
		snps_means = sp.mean(snps, axis=1)
		snps_stds = sp.std(snps, axis=1)
		snps_means.shape = (m, 1)
		snps_stds.shape = (m, 1)
		snps = (snps - snps_means) / snps_stds
		val_set_gen.append(snps)
		if validation:
			val_set_D = 0
		else:
			val_set_D.append(get_sample_D(n=n, m=m, num_sim=100, r2=r2))
		if sp.isnan(val_set_gen).any():
			return simulate_genotypes_w_ld(n, m, n_samples, m_ld_chunk_size=100, r2=0.9)
	return val_set_gen, val_set_D

def pix2sky(header,x,y):
	hdr_info = parse_header(header)
	x0 = x-hdr_info[1][0]+1.	# Plus 1 python->image
	y0 = y-hdr_info[1][1]+1.
	x0 = x0.astype(scipy.float64)
	y0 = y0.astype(scipy.float64)
	x = hdr_info[2][0,0]*x0 + hdr_info[2][0,1]*y0
	y = hdr_info[2][1,0]*x0 + hdr_info[2][1,1]*y0
	if hdr_info[3]=="DEC":
		a = x.copy()
		x = y.copy()
		y = a.copy()
		ra0 = hdr_info[0][1]
		dec0 = hdr_info[0][0]/raddeg
	else:
		ra0 = hdr_info[0][0]
		dec0 = hdr_info[0][1]/raddeg
	if hdr_info[5]=="TAN":
		r_theta = scipy.sqrt(x*x+y*y)/raddeg
		theta = arctan(1./r_theta)
		phi = arctan2(x,-1.*y)
	elif hdr_info[5]=="SIN":
		r_theta = scipy.sqrt(x*x+y*y)/raddeg
		theta = arccos(r_theta)
		phi = artan2(x,-1.*y)
	ra = ra0 + raddeg*arctan2(-1.*cos(theta)*sin(phi-pi),
				   sin(theta)*cos(dec0)-cos(theta)*sin(dec0)*cos(phi-pi))
	dec = raddeg*arcsin(sin(theta)*sin(dec0)+cos(theta)*cos(dec0)*cos(phi-pi))

	return ra,dec

def compute_continuous_prob_value(parameters, distribution, rvs):
    mean = float(parameters[0])
    stddev = float(parameters[1]) / 100 * mean
    A = float(parameters[2])
    B = float(parameters[3])
    result = float("-inf")

    if rvs is None:
        while result <= A or result > B:

            if distribution == "normal":
                rvs = stats.norm.rvs
                result = rvs(mean, stddev)

            elif distribution == "lognormal":
                variance = stddev ** 2.0
                mu = log(mean ** 2.0 / sqrt(variance + mean ** 2.0))
                sigma = sqrt(log((variance / mean ** 2.0) + 1.0))
                rvs = stats.lognorm.rvs
                result = rvs(sigma, scale=scipy.exp(mu))

            elif distribution == "gamma":
                betha = (stddev) ** 2 / mean
                alpha = mean / betha
                rvs = stats.gamma.rvs
                result = rvs(alpha, scale=betha)
    else:
        result = 1

    return result

def relative_bin_deviation(h1, h2): # 79 us @array, 104 us @list \w 100 bins
    r"""
    Calculate the bin-wise deviation between two histograms.
    
    The relative bin deviation between two histograms :math:`H` and :math:`H'` of size
    :math:`m` is defined as:
    
    .. math::
    
        d_{rbd}(H, H') = \sum_{m=1}^M
            \frac{
                \sqrt{(H_m - H'_m)^2}
              }{
                \frac{1}{2}
                \left(
                    \sqrt{H_m^2} +
                    \sqrt{{H'}_m^2}
                \right)
              }
    
    *Attributes:*

    - semimetric (triangle equation satisfied?)
    
    *Attributes for normalized histograms:*

    - :math:`d(H, H')\in[0, \infty)`
    - :math:`d(H, H) = 0`
    - :math:`d(H, H') = d(H', H)`
    
    *Attributes for not-normalized histograms:*

    - :math:`d(H, H')\in[0, \infty)`
    - :math:`d(H, H) = 0`
    - :math:`d(H, H') = d(H', H)`
    
    *Attributes for not-equal histograms:*

    - not applicable 
    
    Parameters
    ----------
    h1 : sequence
        The first histogram.
    h2 : sequence
        The second histogram, same bins as ``h1``.
    
    Returns
    -------
    relative_bin_deviation : float
        Relative bin deviation between the two histograms.
    """
    h1, h2 = __prepare_histogram(h1, h2)
    numerator = scipy.sqrt(scipy.square(h1 - h2))
    denominator = (scipy.sqrt(scipy.square(h1)) + scipy.sqrt(scipy.square(h2))) / 2.
    old_err_state = scipy.seterr(invalid='ignore') # divide through zero only occurs when the bin is zero in both histograms, in which case the division is 0/0 and leads to (and should lead to) 0
    result = numerator / denominator
    scipy.seterr(**old_err_state)
    result[scipy.isnan(result)] = 0 # faster than scipy.nan_to_num, which checks for +inf and -inf also
    return scipy.sum(result)