<jupyter_start><jupyter_text>0. load data<jupyter_code>save_dir = "../temp"
if not os.path.exists(save_dir):
    os.makedirs(save_dir)
data_dir = "../data"
file_name = ['lfp_bred_1.mat', 'lfp_bred_2.mat', 'lfp_bred_3.mat']
rec_time = [0, 0.5]
data = [
    np.concatenate([io.loadmat('%s/%s'%(data_dir, fname))['v4_bred']
                    for fname in file_name], 0)[:,:,:100:2].T,
    np.concatenate([io.loadmat('%s/%s'%(data_dir, fname))['pfc_bred']
                    for fname in file_name], 0)[:,:,:100:2].T]
data[0].shape, data[1].shape
dims = [data[0].shape[1], data[1].shape[1]]
num_time = data[0].shape[0]
obs_trial = data[0].shape[2]<jupyter_output><empty_output><jupyter_text>load LaDynS run on data<jupyter_code>result_exp = np.load('result_exp.npz')
Sigma_exp = result_exp['Sigma_est']
weight_exp = result_exp['weight_est']<jupyter_output><empty_output><jupyter_text>cross precision<jupyter_code>num_blobs = 3
times = np.array([[1, 1, 12], [21, 18, 6], [32, 37, 9]])
H0 = np.zeros((times.shape[0], num_time, num_time))
for i, (time_0, time_1, duration) in enumerate(times):
    H0[i, time_0:time_0+duration, time_1:time_1+duration] \
    = np.abs(np.arange(duration)-np.arange(duration)[:,None]) < 2
import matplotlib.gridspec as gridspec
plt.figure(figsize=(4,3))
gridspec.GridSpec(1,4); plt.subplot2grid((1,4), (0,0), colspan=3, rowspan=1)
ldn.imshow(np.sum(H0, 0), time=[0,50], identity=True, vmax=1)
plt.ylabel(r'series 1 at time $t$'); plt.xlabel(r'series 2 at time $s$')
ax = plt.subplot2grid((1,4), (0,3))
ax.axis('off')<jupyter_output><empty_output><jupyter_text>latent direction<jupyter_code>mu = [np.mean(dat, -1) for dat in data]
covs = [
    np.matmul(dat - m[...,None], dat - m[...,None], axes=[(1,2),(2,1),(1,2)])/obs_trial
    for dat, m in zip(data, mu)]
weights_exp = [np.array(weight_exp[:num_time]),
               np.array(weight_exp[num_time:])]
betas_exp = [(cov @ w[...,None])[...,0] for cov ,w in zip(covs, weights_exp)]<jupyter_output><empty_output><jupyter_text>latent precision<jupyter_code>sighat = Sigma_exp.copy()
sighat_1 = sighat[:num_time,:num_time]
sighat_2 = sighat[num_time:,num_time:]
Y_mm_1 = np.concatenate([[0], np.cumsum(np.sqrt(-2*np.log(
    sighat_1[np.arange(1,num_time), np.arange(0,num_time-1)])))])
Y_mm_2 = np.concatenate([[0], np.cumsum(np.sqrt(-2*np.log(
    sighat_2[np.arange(1,num_time), np.arange(0,num_time-1)])))])
sigmm_1 = np.exp(-(np.square(Y_mm_1[:,None] - Y_mm_1[None,:])/2))
sigmm_2 = np.exp(-(np.square(Y_mm_2[:,None] - Y_mm_2[None,:])/2))
precs_auto = [linalg.inv(S + 0.5 * np.eye(num_time))
              for S in [sigmm_1, sigmm_2]]<jupyter_output><empty_output><jupyter_text>simulation function<jupyter_code>def generate_param(snr):
    prec_latent = np.array(np.block([
        [precs_auto[0], -snr*np.sum(H0,0)],
        [-snr*np.sum(H0,0).T, precs_auto[1]]
    ]))
    cov_latent = linalg.inv(prec_latent)
    cov_latent = cov_latent + np.real(np.min(linalg.eig(Sigma_exp)[0]) - np.min(linalg.eig(cov_latent)[0]))*np.eye(2*num_time)
    sd_latent = np.sqrt(np.diag(cov_latent))
    corr_latent = cov_latent / sd_latent / sd_latent[:,None]
    chol_latent = linalg.cholesky(corr_latent)
    
    return chol_latent
def generate_sim(num_trial, data, chol_latent):
    obs_trial = data[0].shape[2]
    
    dummy = [dat[:,:,np.random.choice(obs_trial, num_trial, replace=False)]
             for dat, dim in zip(data, dims)]
    ldummy = [np.squeeze(np.matmul(dir_i[:,None,:], dumm - np.mean(dumm, -1, keepdims=True)))
              for dir_i, dumm in zip(weights_exp, dummy)]

    # cov_ldummy = np.cov(ldummy[0], ldummy[1])
    # sd_ldummy = np.sqrt(np.diag(cov_ldummy))

    latent_variables = (chol_latent.T @ np.random.normal(0, 1, (2*num_time, num_trial)))

    pops = [dummy[0] + betas_exp[0][:,:,None] * (latent_variables[:num_time,None,:]-ldummy[0][:,None,:]),
            dummy[1] + betas_exp[1][:,:,None] * (latent_variables[num_time:,None,:]-ldummy[1][:,None,:])]
        
    return pops<jupyter_output><empty_output><jupyter_text>1. simulation  preliminary parameters<jupyter_code>offset_cross=10
offset_auto=10
lambdas= np.vstack([
    np.exp(np.linspace(np.log(0.1), np.log(0.001), 11)), #np.full((11,), 0.06), ## lambdas_cross
    np.zeros(11), ## lambdas_auto
    np.full((11,), 0.5), ## lambdas_diag
    np.zeros(11)]) ## lambdas_ridge
time_graph_cross = (np.abs(np.subtract(*np.where(
    np.full([num_time,num_time],True)))
       .reshape([num_time,num_time])) < offset_cross)
time_graph_auto = (np.abs(np.subtract(*np.where(
    np.full([num_time,num_time],True)))
       .reshape([num_time,num_time])) < offset_auto)
snrs = np.array([0.002, 0.004, 0.006, 0.008, 0.010])<jupyter_output><empty_output><jupyter_text>iterate simulation<jupyter_code>num_cv = 10
num_sim = 60
num_bst = 200
num_trial = 3000
for iter_sim in range(num_sim):
    for iter_snr, snr in enumerate(snrs):   
        if os.path.exists("%s/result_sim_%.1d_%.3d.npz"%(save_dir, iter_snr, iter_sim)):
            continue
        
        while True:
            with np.errstate(all='raise'):
                try:
                    start_sim = time.time()

                    # data simulation
                    chol_latent = generate_param(snr)
                    pops = generate_sim(num_trial, data, chol_latent)

                    # choosing tuning parameter
                    Omegas_cv = np.zeros([num_cv, lambdas.shape[-1], 2*num_time, 2*num_time])
                    weights_cv = []
                    for iter_cv in np.arange(num_cv):
                        start_cv = time.time()

                        # random sample
                        trial_sample = np.random.choice(num_trial, num_trial, False)
                        pops_bst = [pops[0], pops[1][:,:,trial_sample]]

                        # estimate for each parameter set
                        for iter_lambda, (lambda_cross, lambda_auto, lambda_diag, lambda_ridge) \
                        in enumerate(zip(*lambdas)): 
                            if iter_cv > 0:
                                init_arg = {'Omega_init': Omegas_cv[0,iter_lambda],
                                            'weight_init': weights_cv[iter_lambda]}
                            elif iter_lambda > 0:
                                init_arg = {'Omega_init': Omegas_cv[0,iter_lambda-1],
                                            'weight_init': weights_cv[iter_lambda-1]}
                            else:
                                init_arg = dict()
                                
                            Omega_lambda, Sigma_lambda, latent_lambda, weight_lambda \
                            = ldn.fit(pops_bst, lambda_diag, lambda_cross, offset_cross,
                                       lambda_auto, offset_auto, ths_glasso=1e-6, ths_lasso=1e-6,
                                       **init_arg)

                            # save result
                            Omegas_cv[iter_cv, iter_lambda] = Omega_lambda
                            if iter_cv == 0:
                                weights_cv.append(weight_lambda)

                        sys.__stdout__.write("%d-%d-%d-th cv finished, lapse: %.3fsec.\n"
                                             %(iter_snr+1, iter_sim+1, iter_cv+1, time.time()-start_cv))
                        sys.__stdout__.flush()
                
                except FloatingPointError:
                    sys.__stdout__.write("%d-%d-%d-th cv encountered FloatingPointError, lapse:%.3fsec.\n"
                                         %(iter_snr+1, iter_sim+1, iter_cv+1, time.time()-start_cv))
                    sys.__stdout__.write("Retrying...\n")
                    sys.__stdout__.flush()
                    continue
                
                sys.__stdout__.write("%d-%d-th simulation cv finished, lapse:%.3fsec.\n"
                                     %(iter_snr+1, iter_sim+1, time.time()-start_sim))
                sys.__stdout__.flush()
                break
                     
        # result tuning parameter
        fdrs_cv = np.sum(Omegas_cv[:,:,:num_time,num_time:]!=0,(0,-2,-1))/np.sum(time_graph_cross)/num_cv
        id_cv = np.argmax(fdrs_cv[fdrs_cv < 0.05])
        lambda_cross, lambda_auto, lambda_diag, lambda_ridge = lambdas.T[id_cv]
        
        # original estimation
        start_est = time.time()
        
        Omega_est, Sigma_est, latent_est, weight_est \
        = ldn.fit(pops, lambda_diag, lambda_cross, offset_cross,
                  lambda_auto, offset_auto, ths_glasso=1e-6, ths_lasso=1e-6)
        
        sys.__stdout__.write("%d-%d-th estimation finished, lapse: %.3fsec.\n"
                             %(iter_snr+1, iter_sim+1, time.time()-start_est))
        sys.__stdout__.flush()
        
        # permutation test
        Omegas_bst = np.zeros([num_bst, 2*num_time, 2*num_time])
        Sigmas_bst = np.zeros([num_bst, 2*num_time, 2*num_time])
        
        init_arg = {'Omega_init': Omegas_cv[0,id_cv],
                    'weight_init': weights_cv[id_cv]}
        for iter_bst in range(num_bst):
            while True:
                with np.errstate(all='raise'):
                    try:
                        start_bst = time.time()

                        # random sample
                        trial_sample = np.random.choice(num_trial, num_trial, True)
                        pops_bst = [pops[0], pops[1][:,:,trial_sample]]

                        # bootstrap
                        Omega_bst, Sigma_bst, latent_bst, weight_bst \
                        = ldn.fit(pops_bst, lambda_diag, lambda_cross, offset_cross,
                                  lambda_auto, offset_auto, ths_glasso=1e-6, ths_lasso=1e-6,
                                  **init_arg)
                        
                    except FloatingPointError:
                        sys.__stdout__.write("%d-%d-%d-th bootstrap encountered FloatingPointError, lapse:%.3fsec.\n"
                                             %(iter_snr+1, iter_sim+1, iter_bst+1, time.time()-start_bst))
                        sys.__stdout__.write("Retrying...\n")
                        sys.__stdout__.flush()
                        continue

                    # save result
                    Omegas_bst[iter_bst] = Omega_bst
                    Sigmas_bst[iter_bst] = Sigma_bst

                    sys.__stdout__.write("%d-%d-%d-th bootstrap finished, lapse: %.3fsec.\n"
                                  %(iter_snr+1, iter_sim+1, iter_bst+1, time.time()-start_bst))
                    sys.__stdout__.flush()
                    break
        
        np.savez("%s/result_sim_%.1d_%.3d"%(save_dir, iter_snr, iter_sim),
                 Omega_est=Omega_est, Sigma_est=Sigma_est,
                 Omegas_perm=Omegas_bst, Sigmas_perm=Sigmas_bst)

        sys.__stdout__.write("%d-th simulation finished, lapse: %.3fsec.\n"
                              %(iter_sim+1, time.time()-start_sim))
        sys.__stdout__.flush()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1440909090909091, 'fraction_numerical': 0.023545454545454546, 'mean_word_length': 2.088433464345873, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '413383', 'license': 'MIT', 'n_tokens_mistral': 4020, 'n_tokens_neox': 3779, 'n_words': 558, 'path': 'example/3_2_iterate_simulation.ipynb', 'repo_name': 'Qi-Xin/ladyns', 'size': '11000'}