an error which indicates a modification of the X array by reference.

fast_lmm_corr_error.py

from __future__ import division, print_function
import numpy as np
import pandas as pd
from pysnptools.snpreader import SnpData
from pysnptools.kernelreader import KernelData
from fastlmm.inference.lmm import LMM

def snp_from_pandas(df, pos = None):
    indiv_ids = zip(df.index.tolist(), ["."] * len(df.index))
    snp_ids = [("_".join(["%s"% y for y in x]) if type(x) in (tuple, list) else "%s"% x) for x in df.columns.tolist()]
    if pos is None:
        pos = [[1,1,0] for x in df.columns.tolist() ]
    elif pos is True:
        pos = [[1,x,x] for x in df.columns.tolist() ]
    elif type(pos) in (list, tuple):
        pass
    else:
        raise ValueError("unknown `pos` type")
#     pos = [[i%2+1,0,0] for i in xrange(len(df.columns))]
    out = SnpData(indiv_ids, snp_ids, df.as_matrix(), pos)
    #out.pos = df.columns.tolist() 
    return out
#########################################################
N = 90
p = 700
pu = 10
delta = 0.5

np.random.seed(20)


snpd =  snp_from_pandas( 
    pd.DataFrame( np.random.randint(0,2, size=(N, p)), 
                 index = ["ind%u" % n for n in range(N)],
                 columns = ["snp%u" % n for n in range(p)],
                ).astype(float) 
    )

K_gt = KernelData(iid=snpd.iid, val=np.corrcoef(snpd.val)).standardize()
beta = pd.Series(np.zeros(p),  index = ["snp%u" % n for n in range(p)])
beta.iloc[N//2] = 10
##
phu =  snp_from_pandas( 
    pd.DataFrame( np.random.randn(*(N, pu)), 
                 index = ["ind%u" % n for n in range(N)],
                 columns = ["ph%u" % n for n in range(pu)],
                ).astype(float) 
    )
K_phu = KernelData(iid=phu.iid, val=np.corrcoef(phu.val)).standardize()
##
rnad =  snp_from_pandas( 
    pd.DataFrame(snpd.val.dot(beta) + \
                 phu.val.dot( delta*np.random.randn( pu)) + \
                 np.random.randn(N), index = ["ind%u" % n for n in range(N)]) 
    )

def check_corr_nans(X):
    XX = X.T.dot(X)
    print(XX.shape)
    print(np.isinf(XX).any())
    print(np.isnan(XX).any())
    # [Sxx,Uxx]= LA.eigh(XX)
    [i,j] = np.where(np.isnan(XX))
    print("nan indices:", list(zip(i,j))[:3])

X = np.c_[np.ones([snpd.iid_count,1]), snpd.val]
check_corr_nans(X)

lmm2 = LMM()
# lmm2.setK(K0=K_ca.val, K1=K_gt.val)
lmm2.setK(K0=K_phu.val, K1=K_gt.val)
lmm2.setX2( X )
lmm2.sety(rnad.val.ravel())
try:
    res2 = lmm2.findA2(penalty = 5.0)
    print(res2)
except Exception as ee:
    print(ee)
    
check_corr_nans(X)