function [betamat,Emat]= gs_select(X,k1,k,m,n,r1,r2,mu0,V0,...
    betadraw,v0,S0,Edraw,f,g,nf,ng,l0,r0,lpos,rpos,Tdraw,y)


R = r1+r2;
betamat=zeros(k,r2); 
Emat=zeros((m*(m+1)/2)-1,r2); %collects unique elements of E as a vector
% recall first element is one, so we can skip a row; 

ff=(1:2:(2*n-1))'; % indexes odd rows only, i.e. X1 stuff, needed for draws of ydraw below
gg=(2:2:2*n)'; % indexes even rows only, i.e. X2 stuff, needed for draws of ydraw below 


ydraw=[Tdraw y]; %convenient for draws of beta and E

thoucount=1000;
% start main loop
i=1;
for i=1:R
   
   %%%%%%%%%%%%%%%%%%%%%%
   % draw beta
   %%%%%%%%%%%%%%%%%%%%%
   
   int1=0; 
   int2=0;
   j=1;
   for j=1:n
       Xi=X{j}; 
       int1=int1+(Xi'/Edraw)*Xi;
       int2=int2+(Xi'/Edraw)*ydraw(j,:)';
   end
   V1=inv(inv(V0)+int1);
   mu1=V1*(inv(V0)*mu0+int2);
   betadraw=mvnrnd(mu1,V1)';
   if i>r1
       betamat(:,i-r1)=betadraw;
   end
 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % draw E
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  v1=v0+n;
  yr= reshape(ydraw',2*n,1);%stacking of yi pairs - 2xn by 1
  Xr=cell2mat(X); %2xn by k
  int3=yr-Xr*betadraw;
  int3=reshape(int3,m,n)';%n x 2
  S1=S0+(int3'*int3); 
  Edraw=iw_sig11(v1,inv(S1),1);
  if i>r1
       Emat(:,i-r1)=[Edraw(1,2);Edraw(2,2)];
   end
 
  %%%%%%%%%%%%%%%%%%%
  % draw ystar
  %%%%%%%%%%%%%%%%%%%
  sig11=Edraw(1,1); %"1", of course
  sig12=Edraw(2,1);
  sig22=Edraw(2,2);
  rho=sig12/(sqrt(sig11*sig22));
  beta1=betadraw(1:k1);
  beta2=betadraw(k1+1:end);
  X1=Xr;
  X2=Xr;
  X1(gg,:)=[]; % keep odd rows only
  X1(:,k1+1:end)=[]; % get rid of trailing zeros to get n by k1
  X2(ff,:)=[]; % keep even rows only
  X2(:,1:k1)=[]; % get rid of leading zeros to get n by k2
    
  %Case 1: Ti=1
  %%%%%%%%%%%%%%%%%%%
  % STEP 1: draw Tstar for T=1 cases, conditional on y(observed)
  mu=X1*beta1+(sig12/sig22)*(y-X2*beta2);%conditional mean
  mu(f)=[]; % eliminate zero amounts, so now mu is length g
  int=tnormrnd(mu,sqrt(sig11*(1-rho^2)),lpos,rpos);
  Tdraw(g,1)=int;
  
  %Case 2: Ti=0
  %%%%%%%%%%%%%%%%%%%
  % STEP 1: draw Tstar for T=1 cases, conditional on y(observed)
  mu=X1*beta1+(sig12/sig22)*(y-X2*beta2);%conditional mean
  mu(g)=[]; 
  int=tnormrnd(mu,sqrt(sig11*(1-rho^2)),l0,r0);
  Tdraw(f,1)=int;
  
  
  ydraw=[Tdraw y]; %update ydraw
    
  if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
      
  
        
    
end


