notebook.community

Edit and run 





In [ ]:

Generate data for illustration #1





In [ ]:



    


%matplotlib inline
from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt

from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression

import pyximport
pyximport.install()
import os

absolute_code_path = '/home/mackelab/Desktop/Projects/Stitching/code/pyLDS_dev/'
os.chdir(absolute_code_path +'pylds')

from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats

from scipy.io import savemat # store results for comparison with Matlab code   


def update(model):
    model.EM_step()
    return model.log_likelihood()                    

#########################
#  set some parameters  #
#########################

n = 3
p = 9
T = 10000

num_sets = 10;

# tentative observation scheme:
obs_scheme = {'sub_pops': (np.arange(0, p//2 + 1), np.arange(p//2, p)),
              'obs_pops': np.array((0,1)),
              'obs_time': np.array((T//2,T))
             }

for idx_d in range(num_sets):
    
    print('generating set #', idx_d)
    pars_true, _ = gen_pars(n, p, u_dim=0, 
                 pars_in=None, 
                 obs_scheme=None,
                 gen_A='full', lts=np.linspace(0.95, 0.98, n),
                 gen_B='random', 
                 gen_Q='identity', 
                 gen_mu0='random', 
                 gen_V0='stable', 
                 gen_C='random', 
                 gen_d='scaled', 
                 gen_R='fraction',
                 diag_R_flag=True,
                 x=None, y=None, u=None)
    
    means_C =  np.sum(pars_true['C'],axis=1)
    pars_true['C'] -= means_C.reshape(p,1)
    norms_C = np.sum(pars_true['C']*pars_true['C'],axis=1)
    pars_true['C'] /= np.sqrt(norms_C.reshape(p,1))


    ###################
    #  generate data  #
    ###################

    truemodel = LDS(
        dynamics_distn=AutoRegression(A=pars_true['A'].copy(),sigma=pars_true['Q'].copy()),
        emission_distn=Regression_diag(A=np.hstack((pars_true['C'].copy(), pars_true['d'].copy().reshape(p,1))),
                                       sigma=pars_true['R'].copy(), affine=True),
                    )
    truemodel.mu_init = pars_true['mu0'].copy()
    truemodel.sigma_init = pars_true['V0'].copy()

    data, stateseq = truemodel.generate(T)

    plt.figure(figsize=(10,10))
    plt.subplot(1,2,1)
    plt.imshow(np.cov(data.T), interpolation='none')
    plt.subplot(1,2,2)
    plt.imshow(np.corrcoef(data.T), interpolation='none')
    plt.colorbar()
    plt.show()
    
    save_file = '../../../results/cosyne_poster/illustration_1/data/LDS_save_idx'+str(idx_d) 
    save_file_m = {'x': truemodel.states_list[0].stateseq, 
                   'y': truemodel.states_list[0].data,
                   'u' : [], 
                   'T' : truemodel.states_list[0].T, 
                   'Trial': len(truemodel.states_list), 
                   'truePars':pars_true,
                   'obsScheme' : obs_scheme}

    savemat(save_file,save_file_m) # does the actual saving
    
    np.savez(save_file, x=truemodel.states_list[0].stateseq,
                        y= truemodel.states_list[0].data,
                        T=truemodel.states_list[0].T, 
                        Trial=len(truemodel.states_list), 
                        truePars=pars_true,
                        sub_pops=obs_scheme['sub_pops'],            
                        obs_time=obs_scheme['obs_time'],            
                        obs_pops=obs_scheme['obs_pops'])



    











In [ ]:



    


means_C =  np.sum(pars_true['C'],axis=1)
pars_true['C'] -= means_C.reshape(p,1)
norms_C = np.sum(pars_true['C']*pars_true['C'],axis=1)
pars_true['C'] /= np.mean(np.sqrt(norms_C.reshape(p,1)))

Generate data for illustration #2





In [ ]:



    


from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt

from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression

import pyximport
pyximport.install()
import os

absolute_code_path = '/home/mackelab/Desktop/Projects/Stitching/code/pyLDS_dev/'
os.chdir(absolute_code_path +'pylds')

from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats

from scipy.io import savemat # store results for comparison with Matlab code   


np.random.seed(1 * 42 + 2)

if False:

    save_file = '../../../results/cosyne_poster/illustration_2/data/LDS_save_idx'+str(idx_d) 
    from scipy.io import savemat # store results for comparison with Matlab code   
    from scipy.linalg import solve_discrete_lyapunov as dtlyap # solve discrete-time Lyapunov equation
    save_file_m = {'x': truemodel.states_list[0].stateseq, 
                   'y': truemodel.states_list[0].data,
                   'u' : [], 
                   'T' : truemodel.states_list[0].T, 
                   'Trial': len(truemodel.states_list), 
                   'truePars':pars_true,
                   'obsScheme' : obs_scheme}

    savemat(save_file,save_file_m) # does the actual saving
    
    np.savez(save_file, x=truemodel.states_list[0].stateseq,
                        y= truemodel.states_list[0].data,
                        T=truemodel.states_list[0].T, 
                        Trial=len(truemodel.states_list), 
                        truePars=pars_true,
                        sub_pops=obs_scheme['sub_pops'],            
                        obs_time=obs_scheme['obs_time'],            
                        obs_pops=obs_scheme['obs_pops'])

Generate data for simulation #1





In [ ]:



    


from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt

from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression

import pyximport
pyximport.install()
import os

absolute_code_path = '/home/mackelab/Desktop/Projects/Stitching/code/pyLDS_dev/'
os.chdir(absolute_code_path +'pylds')

from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats

from scipy.io import savemat # store results for comparison with Matlab code   

np.random.seed(2 * 42 + 1)

def update(model):
    model.EM_step()
    return model.log_likelihood()                    

#########################
#  set some parameters  #
#########################


ns = np.sort(np.hstack([np.arange(1,10)+1,np.arange(1,10)+1]))
num_sets = ns.size;

p = 30
T = 10000

for idx_d in range(num_sets):
    
    n = ns[idx_d]
    
    if np.mod(idx_d,2)==0:
        overlap = ns[idx_d]
    else:
        overlap = 2
        
    # tentative observation scheme:
    obs_scheme = {'sub_pops': (np.arange(0,p//2+np.ceil(overlap/2.).astype(np.int64)), 
                               np.arange(p//2-np.floor(overlap/2.).astype(np.int64),p)),
              'obs_pops': np.array((0,1)),
              'obs_time': np.array((T//2,T))
             }
    
    print('generating set #', idx_d)
    pars_true, _ = gen_pars(n, p, u_dim=0, 
                 pars_in=None, 
                 obs_scheme=None,
                 gen_A='full', lts=np.linspace(0.95, 0.98, n),
                 gen_B='random', 
                 gen_Q='identity', 
                 gen_mu0='random', 
                 gen_V0='stable', 
                 gen_C='random', 
                 gen_d='scaled', 
                 gen_R='fraction',
                 diag_R_flag=True,
                 x=None, y=None, u=None)

    ###################
    #  generate data  #
    ###################

    truemodel = LDS(
        dynamics_distn=AutoRegression(A=pars_true['A'].copy(),sigma=pars_true['Q'].copy()),
        emission_distn=Regression_diag(A=np.hstack((pars_true['C'].copy(), pars_true['d'].copy().reshape(p,1))),
                                       sigma=pars_true['R'].copy(), affine=True),
                    )
    truemodel.mu_init = pars_true['mu0'].copy()
    truemodel.sigma_init = pars_true['V0'].copy()

    data, stateseq = truemodel.generate(T)


    save_file = '../../../results/cosyne_poster/simulation_1/data/LDS_save_n'+str(n)+'_idx'+str(np.mod(idx_d,2)) 
    save_file_m = {'x': truemodel.states_list[0].stateseq, 
                   'y': truemodel.states_list[0].data,
                   'u' : [], 
                   'T' : truemodel.states_list[0].T, 
                   'Trial': len(truemodel.states_list), 
                   'truePars':pars_true,
                   'obsScheme' : obs_scheme}

    savemat(save_file,save_file_m) # does the actual saving
    
    np.savez(save_file, x=truemodel.states_list[0].stateseq,
                        y= truemodel.states_list[0].data,
                        T=truemodel.states_list[0].T, 
                        Trial=len(truemodel.states_list), 
                        truePars=pars_true,
                        sub_pops=obs_scheme['sub_pops'],            
                        obs_time=obs_scheme['obs_time'],            
                        obs_pops=obs_scheme['obs_pops'])

Generate data for simulation #2





In [ ]:



    


from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt

from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression

import pyximport
pyximport.install()
import os

absolute_code_path = '/home/mackelab/Desktop/Projects/Stitching/code/pyLDS_dev/'
os.chdir(absolute_code_path +'pylds')

from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats

from scipy.io import savemat # store results for comparison with Matlab code   

np.random.seed(2 * 42 + 2)

def update(model):
    model.EM_step()
    return model.log_likelihood()                    

#########################
#  set some parameters  #
#########################


num_sets = 10;

n = 5
p = 50
T = 5000

overlap = 2

# tentative observation scheme:
sub_pop1 = np.arange( 0,16)
sub_pop2 = np.arange(14,30)
sub_pop3 = np.hstack([sub_pop1[0],np.arange(sub_pop2[-1], p)])
obs_scheme = {'sub_pops': (sub_pop1, sub_pop2, sub_pop3),
              'obs_pops': np.array((0,1,2)),
              'obs_time': np.array((T//3,T//2,T))
             }

for idx_d in range(num_sets):
    
    
    print('generating set #', idx_d)
    pars_true, _ = gen_pars(n, p, u_dim=0, 
                 pars_in=None, 
                 obs_scheme=None,
                 gen_A='full', lts=np.linspace(0.95, 0.98, n),
                 gen_B='random', 
                 gen_Q='identity', 
                 gen_mu0='random', 
                 gen_V0='stable', 
                 gen_C='random', 
                 gen_d='scaled', 
                 gen_R='fraction',
                 diag_R_flag=True,
                 x=None, y=None, u=None)

    ###################
    #  generate data  #
    ###################

    truemodel = LDS(
        dynamics_distn=AutoRegression(A=pars_true['A'].copy(),sigma=pars_true['Q'].copy()),
        emission_distn=Regression_diag(A=np.hstack((pars_true['C'].copy(), pars_true['d'].copy().reshape(p,1))),
                                       sigma=pars_true['R'].copy(), affine=True),
                    )
    truemodel.mu_init = pars_true['mu0'].copy()
    truemodel.sigma_init = pars_true['V0'].copy()

    data, stateseq = truemodel.generate(T)


    save_file = '../../../results/cosyne_poster/simulation_2/data/LDS_save_idx'+str(idx_d) 
    save_file_m = {'x': truemodel.states_list[0].stateseq, 
                   'y': truemodel.states_list[0].data,
                   'u' : [], 
                   'T' : truemodel.states_list[0].T, 
                   'Trial': len(truemodel.states_list), 
                   'truePars':pars_true,
                   'obsScheme' : obs_scheme}

    savemat(save_file,save_file_m) # does the actual saving
    
    np.savez(save_file, x=truemodel.states_list[0].stateseq,
                        y= truemodel.states_list[0].data,
                        T=truemodel.states_list[0].T, 
                        Trial=len(truemodel.states_list), 
                        truePars=pars_true,
                        sub_pops=obs_scheme['sub_pops'],            
                        obs_time=obs_scheme['obs_time'],            
                        obs_pops=obs_scheme['obs_pops'])

Generate data for experiment #1

generated by Matlab neural network simulation!

Generate data for experiment #2

generated by Matlab neural network simulation!

Content source: mackelab/pyLDS_dev

Similar notebooks:

notebook.community | gallery | about