

%%% need to install matlab toolboxes:
%%%   FastICA http://research.ics.aalto.fi/ica/fastica
%%%   FSLNets (downloaded later than 26 August 2016) http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSLNets
%%%   htmlTableToCell http://uk.mathworks.com/matlabcentral/fileexchange/34968-htmltabletocell
%%%   nearestSPD http://uk.mathworks.com/matlabcentral/fileexchange/42885-nearestspd
%%%   FDR https://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/software/fdr
addpath ~/matlab ~/matlab/FSLNets ~/matlab/FastICA_25

%%% Import big text file with all imaging IDPs (except rfMRI IDPs), obtained from UK Biobank.
%%% This should be arranged with one column per IDP (starting with subject ID),
%%% using exactly the same IDP ordering as found in the IDPinfo.txt file found in
%%% http://www.fmrib.ox.ac.uk/ukbiobank/gwaspaper
%%% The top line should list the IDP names: a horizontal version of the first column in IDPinfo.txt
%%% The second line should be the first subject's IDP values, etc.
IDPs=importdata('IDPs.txt');
ALL_IDPs=IDPs.data;
N=size(ALL_IDPs,1);

%%% load rfMRI partial correlation netmats and node-amplitudes (i.e. from separate files than the above)
%%% these files should come (from UK Biobank) with one line per subject, no header line, and
%%% for each line, start with the subject ID and then have the unwrapped netmat or node amplitudes.
d=25; 
poop=load(sprintf('10k_rfMRI_d%d_partialcorr_v1.txt',d));
grot=zeros(N,size(poop,2)-1)/0;  for i = 1:size(poop,1), j=find(ALL_IDPs(:,1)==poop(i,1)); if size(j,1)>0, grot(j,:)=poop(i,2:end); end; end;
NET25=grot;
poop=load(sprintf('10k_rfMRI_d%d_NodeAmplitudes_v1.txt',d));
grot=zeros(N,size(poop,2)-1)/0;  for i = 1:size(poop,1), j=find(ALL_IDPs(:,1)==poop(i,1)); if size(j,1)>0, grot(j,:)=poop(i,2:end); end; end;
NODEamps25=grot;
d=100;
poop=load(sprintf('10k_rfMRI_d%d_partialcorr_v1.txt',d));
grot=zeros(N,size(poop,2)-1)/0;  for i = 1:size(poop,1), j=find(ALL_IDPs(:,1)==poop(i,1)); if size(j,1)>0, grot(j,:)=poop(i,2:end); end; end;
NET100=grot;
poop=load(sprintf('10k_rfMRI_d%d_NodeAmplitudes_v1.txt',d));
grot=zeros(N,size(poop,2)-1)/0;  for i = 1:size(poop,1), j=find(ALL_IDPs(:,1)==poop(i,1)); if size(j,1)>0, grot(j,:)=poop(i,2:end); end; end;
NODEamps100=grot;

%%% get the list of IDP sub-modalities
! cat IDPinfo.txt | awk '{print $1}' > grot.txt
po=fopen(sprintf('grot.txt'),'r');  IDPnames=textscan(po,'%s');  fclose(po);  IDPnames0=IDPnames{1};
! cat new_IDPs_1_info.txt | awk '{print $1}' > grot.txt
po=fopen(sprintf('grot.txt'),'r');  IDPnames=textscan(po,'%s');  fclose(po);  IDPnames=IDPnames{1}; IDPnames1=IDPnames(2:end);
! cat new_IDPs_2_info.txt | awk '{print $1}' > grot.txt
po=fopen(sprintf('grot.txt'),'r');  IDPnames=textscan(po,'%s');  fclose(po);  IDPnames=IDPnames{1}; IDPnames2=IDPnames(2:end);
IDPnames=[ IDPnames0(1:43); IDPnames2; IDPnames1; IDPnames0(44:end) ];
j=length(IDPnames)+1;
for i=1:size(NODEamps25,2),   IDPnames{j}=sprintf('rfMRI amplitudes (ICA25 node %d)',i);   j=j+1;  end;
for i=1:size(NODEamps100,2),  IDPnames{j}=sprintf('rfMRI amplitudes (ICA100 node %d)',i);  j=j+1;  end;
for i=1:size(NET25,2),        IDPnames{j}=sprintf('rfMRI connectivity (ICA25 edge %d)',i);        j=j+1;  end;
for i=1:size(NET100,2),       IDPnames{j}=sprintf('rfMRI connectivity (ICA100 edge %d)',i);       j=j+1;  end;

%%% OK for now let's just code it by hand
IDP_modality_types(1:18)=0;       % QC
IDP_modality_types(19:182)=1;     % T1
IDP_modality_types(183)=2;        % T2
IDP_modality_types(184:197)=3;    % SWI
IDP_modality_types(198:213)=4;    % tfMRI
IDP_modality_types(214:888)=5;    % dMRI
IDP_modality_types(889:964)=6;    % rfMRI amplitudes
IDP_modality_types(965:size(IDPnames,1))=7; % rfMRI netmats


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%  FreeSurfer dude

% ! cat aseg_data.txt | head -n 1 > FS_IDPs.txt
% ! cat aparc_data.txt | head -n 1 >> FS_IDPs.txt
! cat aseg_data.txt | tail -n +2 > ~/grot
ASEG=load('~/grot');
! cat aparc_data.txt | tail -n +2 > ~/grot
SURF=load('~/grot');
! /bin/rm ~/grot
FSnames = dataread('file', 'FS_IDPs.txt', '%s', 'delimiter', '\n');
FS=zeros(size(ALL_IDPs,1),size(ASEG,2)+size(SURF,2)-2)/0;
[~,grotI,grotJ]=intersect(ASEG(:,1),ALL_IDPs(:,1));
FS(grotJ,1:size(ASEG,2)-1)=ASEG(grotI,2:end);
[~,grotI,grotJ]=intersect(SURF(:,1),ALL_IDPs(:,1));
FS(grotJ,size(ASEG,2):end)=SURF(grotI,2:end);
clear ASEG SURF;

FS =       FS(:,[91:122 130:156 157:267 305:416 1618:1728 1766:1877]); % main FS variable sub-selection
FSnames=FSnames([91:122 130:156 157:267 305:416 1618:1728 1766:1877]);

grot = load('fs_indices.2017.10.13.txt'); % secondary FS variable sub-selection
FSnames=FSnames(grot);
FS=FS(:,grot);

po=fopen(sprintf('FSnames.txt'),'w');
for i=1:length(FSnames)
%for i=1:65
  grot=FSnames{i};
  grot=strrep(grot,'_',' ');
  if (i<60)
    grot=strrep(grot,'-',' ');
    grot=strrep(grot,'Inf Lat Vent','Inferior Lateral Ventricle');
    grot=strrep(grot,' CC ',' Corpus Callosum (5mm slab) ');
    grot=strrep(grot,'BrainSegVol','Brain Segmentation Volume');
    grot=strrep(grot,'lhC','left-hemisphere C');
    grot=strrep(grot,'rhC','right-hemisphere C');
  end
  grot=strrep(grot,' ant ',' anterior ');
  grot=strrep(grot,' inf ',' inferior ');
  grot=strrep(grot,' sup ',' superior ');
  grot=strrep(grot,' lat ',' lateral ');
  grot=strrep(grot,' lh ',' left-hemisphere ');
  grot=strrep(grot,' rh ',' right-hemisphere ');
  grot=strrep(grot,' S ',' sulcus ');
  grot=strrep(grot,' G ',' gyrus ');
  grot=strrep(grot,' G&S ',' gyrus&sulcus ');
  grot=strrep(grot,'Lat Fis','lateral-fissure');
  grot=strrep(grot,'med-olfact','medial-olfactory');
  grot=strrep(grot,'Ins lg','insular long');
  grot=strrep(grot,'a2009s','Destrieux atlas');
  grot=strrep(grot,'DKTatlas','Desikan-Killiany-Tourville atlas');
  disp(grot);
  fprintf(po,'%d %s   FreeSurfer %s\n',i,FSnames{i},grot);
end
fclose(po);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% get full list of subjects and vars
%%% load non-IDP data obtained from UK Biobank (including info from the supplied HTML doc files).
%%% The data text files should be arranged to have one line per subject (no header line),
%%% and with one column for each of the rows in the main initial table in the HTML.
%%% Hence the first column corresponds to variable number "0" in the HTML.
%%% The data text file should be preprocessed to only have numerical (and NaN) values, and be space-separated.
%%% Hence dates and/or times of importance should be converted to pure numerical form.
varsSUBS=[]; varsVARS=[];
%%% The kk list of numbers below corresponds to if you have multiple distinct sets of variables obtained
%%% from UK Biobank (ie multiple pairs of HTML and data text files).
kk=[ 8252 9268 ]; % put in multiple entries to load complementary downloads
for k = kk
  varsORIG=load(sprintf('SMS/ukb%d_Oct2015.txt',k));
  varsSUBS=sort(unique([varsSUBS;varsORIG(:,1)]));
  clear table;  table.idTableBy.plaintextPreceedingTable='Ukb2docs';  varsHTML=htmlTableToCell(sprintf('SMS/ukb%d.html',k),table);
  varsVARS=sort(unique([ varsVARS {varsHTML{3:end,2}}]));
  % [k length(varsSUBS) length(varsVARS)]
end
%%% incrementally fill out full vars matrix and headers lists, overwriting earlier versions with later ones
vars=zeros(length(varsSUBS),length(varsVARS))/0;  clear varsHeader; clear varsSource;
for k = kk
  varsORIG=load(sprintf('SMS/ukb%d_Oct2015.txt',k));
  clear table;  table.idTableBy.plaintextPreceedingTable='Ukb2docs';  varsHTML=htmlTableToCell(sprintf('SMS/ukb%d.html',k),table);
  [~,grotSa,grotSb]=intersect(varsSUBS,varsORIG(:,1));
  [~,grotVa,grotVb]=intersect(varsVARS,{varsHTML{3:end,2}});
  vars(grotSa,grotVa)=varsORIG(grotSb,1+grotVb);
  varsHeader(grotVa')={varsHTML{2+grotVb,5}};
  varsSource(grotVa')=k;
end

%%% fill out empty descriptions (occurs for repeat visits), and insert ". " before "Uses"
for i=size(varsHeader,2):-1:1
  if prod(size(varsHeader{i}))==0
    j=i-1;
    while prod(size(varsHeader{j}))==0
      j=j-1;
    end
    varsHeader{i}=[varsHeader{j} sprintf(' visit-%d',1+i-j)];
  end
  j=regexp(varsHeader{i},'Uses');
  if prod(size(j))>0
    grot=varsHeader{i};  varsHeader{i} = [ grot(1:j-1) '. ' grot(j:end) ];
  end
end

%%% discard non-imaged subjects
varsORIG=vars; vars=zeros(N,size(varsORIG,2))/0;
for i=1:N
  j=find(ALL_IDPs(i,1)==varsSUBS);
  if prod(size(j))==1
    vars(i,:)=varsORIG(j,:);
  else
    ALL_IDPs(i,1)
  end
end
%clear varsORIG;   % can clear this once we're sure we don't want the full subject list again

%%% identify remaining bad vars on the basis of too much missing data, etc
badvars=[];
for i=1:size(vars,2)
  Y=vars(:,i);
  if sum(abs(mod(Y(~isnan(Y)),1))) == 0  % is this an integer or categorical variable?
    Y(Y<0)=NaN;  % I *think* that in this case negative numbers should be treated as NaN
  end
  grotKEEP=~isnan(Y);  
  grot=(Y(grotKEEP)-median(Y(grotKEEP))).^2; grot=max(grot/mean(grot));   grot=0;  % do we have extreme outliers?  currently ignored
  grotSTD=std(Y(grotKEEP)-mean(Y(grotKEEP))); % for huge constant values this isn't zero using juse std(Y)
  if (sum(grotKEEP)>50) & (grotSTD>0) & (max(sum(nets_class_vectomat(Y(grotKEEP))))/length(Y(grotKEEP))<0.95) & (grot<100)
      % the 3rd thing above is:  is the size of the largest equal-values-group too large?
    vars(:,i)=Y;
  else
    %[i sum(grotKEEP) std(Y(grotKEEP)) max(sum(nets_class_vectomat(Y(grotKEEP))))/length(Y(grotKEEP)) grot]
    badvars=[badvars i];
  end
end
badvars=unique(badvars);  varskeep=setdiff([1:size(vars,2)],badvars);

%%% identify "bad" vars on the basis of redundancy
for i=varskeep
  for j=varskeep
    if j>i & corr(vars(:,i),vars(:,j),'rows','pairwise')>0.99 & corr(isnan(vars(:,i)),isnan(vars(:,j)))>0.2
      if sum(~isnan(vars(:,i))) > sum(~isnan(vars(:,j)))
        k=j;
      else
        k=i;
      end
      %[i j sum(~isnan(vars(:,i))) sum(~isnan(vars(:,j))) k]
      badvars=[badvars k];
    end
  end
end
badvars=unique(badvars);  varskeep=setdiff([1:size(vars,2)],badvars);

%%%%%%%%%%%%%%%%%%%% categories
% 1  age, sex, brain MRI protocol Phase
% 2  genetics
% 3  early life factors
% 10 lifestyle and environment - general
% 11 lifestyle and environment - exercise and work
% 12 lifestyle and environment - food and drink
% 13 lifestyle and environment - alcohol
% 14 lifestyle and environment - tobacco
% 20 physical measures - general
% 21 physical measures - bone density and sizes
% 22 physical measures - cardiac & blood vessels
% 23 hearing test
% 24 eye test
% 25 physical activity measures
% 30 blood assays
% 31 brain IDPs
% 32 cognitive phenotypes
% 50 health and medical history, health outcomes
% 51 mental health self-report
% 99 misc, ignored

vt{1}=[31 34 22200 25780];
vt{2}=[21000 22000:22125 22201:22325];
vt{3}=[52 129 130 1677 1687 1697 1737 1767 1777 1787 20022];
vt{10}=[3:6 189 680 699 709 728 738 1031 1797 1807 1835 1845 1873 1883 2139 2149 2375 2385 2395 2405 2267 ...
        2277 2714:10:2834 2946 3526 3536 3546 3581 3591 3659 3669 3700 3710 3720 3829 3839 3849 3872 3882 3912 3942 3972 3982 4501 4674 4825 4836 5057 ...
        6138 6142 6139:6141 6145:6146 6160 10016 10105 10114 10721 10722 10740 10749 10860 10877 10886 ...
        20074:20075 20110:20113 20118:20119 20121 22501 22599 22606 22700 22702 22704];
vt{11}=[1001 1011 767 777 796 806 816 826 845 864 874 884 894 904 914 924 943 971 981 991 1021 1050:10:1220 2624 2634 ...
        3426 3637 3647 6143 6162 6164 10953 10962 10971 22604 22605 22607:22615 22620 22630 22631 22640:22655 104900 104910 104920];
vt{12}=[1289:10:1389 1408:10:1548 2654 3680 6144 10007 10767 10776 10855 10912 20084:20094 20098:20109 100001:100009 100011:100019 100021:100025 ...
        100010:10:100560 100760:10:104670];
vt{13}=[1558:10:1628 2664 3731 3859 4407 4418 4429 4440 4451 4462 5364 20095:20097 20117 100580:10:100740];
vt{14}=[1239:10:1279 2644 2867:10:2907 2926 2936 3159 3436:10:3506 5959 6157 6158 6183 6194 20116 22506:22508];
vt{20}=[46:50 1707 1717 1727 1747 1757 2306 3062:3065 3088 3089 3160 10691 10694:10696 20015 21001 21002 22400:22410 22427 23098:23130 23244:23289];
vt{21}=[77 78 3083:3085 3143:3144 3148 4100:4101 4104:4106 4119:4120 4124:4125 4138:4141 4143:4146 23200:23243 23290:23320];
vt{22}=[93:95 102 4079 4080 4194 4195 4196 5983 5984 5986 5992 5993 6019 6020 6022 6024 6032:6034 6038 6039 ...
        12673:12687 12697 12698 12702 21021 22420:22426];
vt{23}=[4229 4230 4232 4233 4236 4240:4244 4268:4270 4275:4277 10793 20019 20021];
vt{24}=[5084:5089 5096:5119 5132:5135 5149 5155:5164 5181 5198 5221 5237 5251 5254:5257 5262:5265 5276 5292 5306 5324:5325 5327 6070 6072];
vt{25}=[90002:90003 90011:90013 90015 90019:90083 90086:90089 90091:90150 90159:90177 90179 90182:90195];
vt{30}=[74 30000:10:30300 30104 30112 30114 30172 30174 30242 30252 30254];
vt{31}=[25000:25746 25761:25768 25781:25920];
vt{32}=[111 398:404 630 4250:4256 4258:4259 4282:4283 4285:4286 4288:4292 4957 4968 4979 4990 5001 5012 5556 5699 5779 5790 5866 ...
        10137:10141 10144 10146:10147 10241 10609 10612 20016 20018 20023 20082 20128:20157 20159 20169:2:20191 20195 20196 20198 20200 ...
        20229 20230 20240 20242 20244:20248];
vt{50}=[84 87 92 134:137 2178 2188 2207:10:2257 2296 2316 2335:10:2365 2415 2443:10:2473 2492 2674 ...
        2684 2704 2844 2956:10:2986 3005 3079 3140 ...
        3393 3404 3414 3571 3606 3616 3627 3741 3751 3761 3773 3786 3799 3809 3894 3992 ...
        4012 4022 4041 4056 4067 4689 4700 4717 4728 4792 4803 4814 5408 5419 5430 5441 5452 5463 5474 5485 5496 5507 ...
        5518 5529 5540 5610 5832 5843 5855 5877 5890 5901 5912 5923 5934 5945 6119 6147 6148 6149 ...
        6150 6151 6152 6153 6154 6155 6159 6177 6179 6205 20001:20011 22126:22181 22502:22505 22616 22618 22619 40001:41253 42000:42013];
vt{51}=[1920:10:2110 4526 4537 4548 4559 4570 4581 4598 4609 4620 4631 4653 5375 5386 5663 5674 6156];
vt{99}=[19 21 53 54 68 120 132 757 1647 2129 3137 12139 12140 12141 12187 12188 12223 12224 12253 12254 12623 12624 12651 12652 ...
        12663 12664 12671 12695 12699 12700 12704 12706 12848 12851 12854 ...
        20012:20014 20024 20049 20053:20054 20061:20062 20077 20078 20083 20115 20158 20201:20227 20249:20253 21003 22499 22500 ...
        22600:22603 22617 22660:22664 25747:25753 40000 105010 105030 110005 110006];

vt_use=[ 3 10:14 20:25 30 32 51 ];
vt_not_use=[ 99 1 2 31 50 ];

%%% check if there are still existing vars not listed by hand
grotI=[];  for i=varskeep
  ii=varsVARS(i);  ii=ii{1};  ii(strfind(ii,'-'):end)=[];  grotI=unique([grotI str2num(ii)]);
end
grotII=[];  for i=[vt_use vt_not_use]
  grotII=unique([grotII vt{i}]);
end
setdiff(grotI,grotII)  % should be empty if we already list above all relevant vars.

%%% sanity check of any given vt (variable type)
for i=vt{20}
  j=nets_cellfind(varsVARS,sprintf('%d-',i));
  if ~isempty(j), disp( sprintf('%d %s',i,varsHeader{j(1)} ) ); end;
end

%%% identify vars of special interest
dosI=nets_cellfind(varsVARS,'53-2.0'); scan_date=vars(:,dosI); sum(isnan(scan_date))   
  scan_date(isnan(scan_date))=max(scan_date)+0.1; % sometimes needed with missing dates
yobI=nets_cellfind(varsVARS,'34-0.0'); mobI=nets_cellfind(varsVARS,'52-0.0'); birth_date=vars(:,yobI)+(vars(:,mobI)-0.5)/12;
age=scan_date - birth_date;
sexI=nets_cellfind(varsVARS,'31-0.0');  sex=vars(:,sexI);
badvars=[ badvars dosI yobI mobI sexI ];

%%% discard some more vars
for ii = vt_not_use
  ii
  for i= vt{ii}
    i
    j=nets_cellfind(varsVARS,sprintf('%d-',i));
    if ~isempty(j)
      %disp([i j])
      %disp( sprintf('%d %s',i,varsHeader{j(1)} ) )
      badvars=[badvars j];
    end
  end
end
badvars=unique(badvars);  varskeep=setdiff([1:size(vars,2)],badvars);  varskeepOLD=varskeep;

%%% reorder kept vars according to nice categorisations
varskeepNEW=[]; varskeepVT=[];
for ii=vt_use
  for i = vt{ii}
    j=nets_cellfind(varsVARS,sprintf('%d-',i));
    for jj=j
      if ~isempty(find(varskeep==jj))
        varskeepNEW=[varskeepNEW jj];
        varskeepVT=[varskeepVT ii];
      end
    end
  end
end
varskeep=varskeepNEW;
for i=1:length(varskeep)
   disp( sprintf('%d %d %s',i,varskeepVT(i),varsHeader{varskeep(i)} ) )
end

%%% sanitise categorical variables to make them more quantitative
varsK=vars(:,varskeep);
[~,grotDC,~]=intersect(varskeep,nets_cellfind(varsHeader,' data-coding ')); clear poop; poop2=[];
for i=grotDC'
  Y=varsK(:,i);
  if sum(abs(mod(Y(~isnan(Y)),1))) == 0  % is this an integer or categorical variable?
    Y(Y<0)=NaN;  % I *think* that in this case negative numbers should be treated as NaN
    poop(i)=sum(Y==555)/sum(~isnan(Y));
    if (poop(i)>0.0025)
      grot=varsHeader{varskeep(i)};
      grot=grot(strfind(grot,' data-coding '):end);
      grot=grot(14:end);
      grot=str2num(grot(1:min(strfind(grot,' '))-1));
      poop2=[poop2 grot];
      Y(Y==555)=0.5;
      Y(Y==444)=0.25;
      Y(Y==111)=NaN; % "varied" (WTF?)
      Y(Y>50)= max(Y(Y<50))+1;
      unique(Y(~isnan(Y)))'
      varsK(:,i)=Y;
    end
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% prepare confounds and main working matrices

%%% table position
grot=load('10129subjects_HeadTablePos_03August17.txt'); % see recent showcase variables 25756-9 to get these
TablePos=zeros(N,size(grot,2)-2)/0;
for i=1:N                                                                                                                                                           
  j=find(ALL_IDPs(i,1)==grot(:,1));                                                                                                                                 
  if ~isempty(j)                                                                                                                                                    
    TablePos(i,:)=grot(j,2:end-1);     %%% = [ FOV-X FOV-Y FOV-Z table-Z ]
  end;                                                                                                                                                              
end                                                                                                                                                                 

%%%% create date-drift confounds - beta
%%%  maybe should remove other confounds (including bed position) before doing this

%%% median-based standardise data, then remove outliers
%%%% (so grotR is a temporary copy of all IDPs *including* QC-IDPs)
grotR=[ALL_IDPs(:,2:end) NODEamps25 NODEamps100 NET25 NET100];
grotR=grotR-repmat(nanmedian(grotR),size(grotR,1),1); grotR=grotR./repmat(nanmedian(abs(grotR)),size(grotR,1),1); grotR(grotR>5)=5; grotR(grotR<-5)=-5;

%%%% impute missing data in grotR using low-rank-with-shrinkage imputation
%%%% (again, this is all just for the temporary version of the IDPs used to create the date confounds)
[impute_replace,grot_impute_noreplace,grotpcaU,grotpcaS,grotpcaV]=nets_impute(grotR,1000);

%%%% spline-based smoothing of each IDP separately, on a regular date grid
%%%% SmoothingParam is set to give a reasonably large amount of temporal smoothing (but not too large)
%%%% thin-plate-spline with decent end-effect-treatment would be better
grotR=impute_replace; grotX=[min(scan_date):1/365:max(scan_date)]'; grotRs=zeros(length(grotX),size(grotR,2));
for i=1:size(grotR,2)
  i
  grot=[scan_date grotR(:,i)];  grot=grot(sum(isnan(grot),2)==0,:);
  grotf=fit(grot(:,1),grot(:,2),'smoothingspline','SmoothingParam',0.999);
  grotRs(:,i)=feval(grotf,grotX);
end

%%%% PCA to get top (smoothed) temporal eigenvectors
[grotUs,grotSs,grotVs]=nets_svds(nets_demean(grotRs),0);
grotUs=grotUs(:,1:10); % top 10 confounds (less would probably be OK)

%%%% sample back onto original dates, ie for each subject separately
confDATE=zeros(length(scan_date),size(grotUs,2));
for i=1:length(scan_date)
  grot=abs(grotX-scan_date(i)); % for this subject, find the best date-point in the regular date grid
  grot=find(grot==min(grot));
  confDATE(i,:)=grotUs(grot(1),:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%% make confounds

%%% in general we are including two versions of most confound variables, to make the final confound
%%% matrix as "powerful" as possible;
%%%   - one version with MAD-based outlier cleanup
%%%   - one version with inormal Gaussianisation

% first define some common working variables, all demeaned
a=nets_demean(age); a2=nets_demean(a.^2); s=nets_demean(sex);
q=nets_demean(TablePos(:,1:4));  % confounds relating to bed position in scanner (x,y,z,table)

%%% age-related confounds
confA1=[a a2 a.*s a2.*s];   % note - no MAD-based outlier removal for age confounds - doesn't work and isn't needed.
confA2=inormal(confA1);     % fully "Gaussianised" version of the same

%%% age-unrelated confounds
%%% ALL_IDPs(:,[12 13 18]) = rfMRI head motion  /  tfMRI head motion  /  head size scaling
grot=[ALL_IDPs(:,[12 13 18]) ALL_IDPs(:,[12 13]).^2 q q.^2 confDATE];
grot=grot-repmat(nanmedian(grot),size(grot,1),1); grot=grot./repmat(nanmedian(abs(grot)),size(grot,1),1); grot(abs(grot)>5)=NaN;
confB1=grot;

confB2=inormal([ALL_IDPs(:,[12 13 18]) q q.^2 confDATE ]);

%%% create final full confound matrix, with standardisation (simple Z transform) and set missing data = 0
conf = nets_normalise([ confA1 confA2 confB1 confB2 s ]);  conf(isnan(conf))=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%  create deconfounded and de-outliered versions of the IDPs
%%%%  ALL_IDPs_i = deoutliered (but not deconfounded)
%%%%  ALL_IDPs_i_deconf = deoutliered and fully deconfounded
%%%%  ALL_IDPs_i_deconfnoage = deoutliered and deconfounded (but no age-related deconfounding)
%%%%  _i_ is Gaussianised version;  _m_ is MAD de-outliering

grot = [ ALL_IDPs(:,19:end) NODEamps25 NODEamps100 NET25 NET100 ]; % all IDPs except QC-IDPs
ALL_IDPs_i = inormal(grot);
grot=grot-repmat(nanmedian(grot),size(grot,1),1); grot=grot./repmat(nanmedian(abs(grot)),size(grot,1),1); grot(abs(grot)>5)=NaN;
ALL_IDPs_m=nets_normalise(grot);
ALL_IDPs_i_deconf      = nets_unconfound(ALL_IDPs_i,conf);
ALL_IDPs_i_deconfnoage = nets_unconfound(ALL_IDPs_i,conf(:,(size(confA1,2)+size(confA2,2)+1):end));
ALL_IDPs_m_deconf      = nets_unconfound(ALL_IDPs_m,conf);
ALL_IDPs_m_deconfnoage = nets_unconfound(ALL_IDPs_m,conf(:,(size(confA1,2)+size(confA2,2)+1):end));

%%%%%%%%%%%%%%%%%%%%%%%%%%% FreeSurfer  dude

grot=FS;
FS_i = inormal(grot);
grot=grot-repmat(nanmedian(grot),size(grot,1),1); grot=grot./repmat(nanmedian(abs(grot)),size(grot,1),1); grot(abs(grot)>5)=NaN;
FS_m=nets_normalise(grot);
FS_i_deconf      = nets_unconfound(FS_i,conf);
FS_i_deconfnoage = nets_unconfound(FS_i,conf(:,(size(confA1,2)+size(confA2,2)+1):end));
FS_m_deconf      = nets_unconfound(FS_m,conf);
FS_m_deconfnoage = nets_unconfound(FS_m,conf(:,(size(confA1,2)+size(confA2,2)+1):end));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

grot=varsK;
vars_i = inormal(grot);
grot=grot-repmat(nanmedian(grot),size(grot,1),1); grot=grot./repmat(nanmedian(abs(grot)),size(grot,1),1); grot(abs(grot)>5)=NaN;
vars_m=nets_normalise(grot);
vars_i_deconf      = nets_unconfound(vars_i,conf);
vars_i_deconfnoage = nets_unconfound(vars_i,conf(:,(size(confA1,2)+size(confA2,2)+1):end));
vars_m_deconf      = nets_unconfound(vars_m,conf);
vars_m_deconfnoage = nets_unconfound(vars_m,conf(:,(size(confA1,2)+size(confA2,2)+1):end));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% genetic (population) confounds (including data subsetting for white-British, etc)

%%%% find which subjects were used in the main genetics stuff for white-British?
KG_IDs=load('ukb810_8428_ids.txt'); % see paper for definition, to sort this for your own ID list
%%KG_IDs=load('ukb810_930_ids.txt'); % initial modest replication sample
[~,KG,~]=intersect(ALL_IDPs(:,1),KG_IDs);

%%%% new GPCs (top 40 PCA components and associated IDs)
! cat ukb_sqc_v2.txt | cut -d ' ' -f26-65 > ~/grotgrot
GPCs=load('~/grotgrot');
! cat ukb810_cal_chr1_v2_s488374.fam | awk '{print $1}' > ~/grotgrot
GPCsID=load('~/grotgrot');
grot=zeros(size(ALL_IDPs,1),40)/0;  [~,grotI,grotJ]=intersect(GPCsID,ALL_IDPs(:,1));  grot(grotJ,:)=GPCs(grotI,:);  GPCs=grot; % reduce to imaged subjects
GPCs=GPCs(KG,:); 
GPCs=nets_normalise(GPCs); GPCs(isnan(GPCs))=0;  % small number of NaNs, so zero-centre things and replace NaNs with 0.

ALL_IDPs_i_KG_deconf_GPCdeconf=nets_unconfound(ALL_IDPs_i(KG,:),[conf(KG,:) GPCs]); % deconfounded white-British-subsetted IDPs
FS_i_KG_deconf_GPCdeconf=nets_unconfound(FS_i(KG,:),[conf(KG,:) GPCs]);
vars_i_KG_deconf_GPCdeconf=nets_unconfound(vars_i(KG,:),[conf(KG,:) GPCs]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%% body confounds

BV=vars_i_KG_deconf_GPCdeconf;

[~,grot,~]=intersect(varskeep,nets_cellfind(varsVARS,'4079-')); [~,~,grot,~,grotV]=nets_impute(BV(:,grot),2); BPd=grot(:,1)*sign(grotV(1,1));
[~,grot,~]=intersect(varskeep,nets_cellfind(varsVARS,'4080-')); [~,~,grot,~,grotV]=nets_impute(BV(:,grot),2); BPs=grot(:,1)*sign(grotV(1,1));

[~,grot,~]=intersect(varskeep,nets_cellfind(varsVARS,'50-')); [~,~,grot,~,grotV]=nets_impute(BV(:,grot),2); Height=grot(:,1)*sign(grotV(1,1));
[~,grot,~]=intersect(varskeep,nets_cellfind(varsVARS,'21002-')); Weight=nets_impute(BV(:,grot),2);

% get Head measures
[~,grot,~]=intersect(varskeep,nets_cellfind(varsVARS,'23306-')); BONE_HeadBMC=BV(:,grot);
[~,grot,~]=intersect(varskeep,nets_cellfind(varsVARS,'23226-')); BONE_HeadBMD=BV(:,grot);
grotB=~isnan( BONE_HeadBMC .* BONE_HeadBMD);
%%%% then run the above code with J=21 to see the list of the top 20 bone things that correlate with IDPs
J=21; grot=find(varskeepVT==J); [groty,groti]=sort(max(-log10(unicorrPi_deconf(:,grot))));  varsHeader(varskeep(grot(groti)))'
%%% run PCA
[grotX,~,grotU,grotS,grotV]=nets_impute(BV(grotB,grot(groti(end:-1:end-20))),10);
% then view grotV and take the components which are most relevant to head, and (eg) correlate with predage
BONE_PC1=0*BONE_HeadBMC/0; BONE_PC1(grotB)=grotU(:,1)*sign(grotV(1,1));
BONE_PC2=0*BONE_HeadBMC/0; BONE_PC2(grotB)=grotU(:,2)*sign(grotV(1,2));

BodyConf=[BPd BPs Height Weight BONE_HeadBMC BONE_HeadBMD BONE_PC1 BONE_PC2];
BodyConf1=nets_impute(nets_normalise(BodyConf),6);   % impute missing data
BodyConf2=nets_normalise(BodyConf(:,5:8)); BodyConf2(isnan(BodyConf2))=0;  % zero missing data
BodyConf3=nets_normalise([BodyConf1 BodyConf1.^2 BodyConf2 BodyConf2.^2]);  % all possible combinations and powers....

ALL_IDPs_i_KG_deconf_GPCdeconf_BODYdeconf =nets_unconfound(ALL_IDPs_i(KG,:), [conf(KG,:) GPCs BodyConf3]);
FS_i_KG_deconf_GPCdeconf_BODYdeconf       =nets_unconfound(FS_i(KG,:),       [conf(KG,:) GPCs BodyConf3]);
vars_i_KG_deconf_GPCdeconf_BODYdeconf     =nets_unconfound(vars_i(KG,:),     [conf(KG,:) GPCs BodyConf3]);


po=fopen(sprintf('BodyConfounds09Aug2017.txt'),'w');
fprintf(po,'ID BloodPressureDiastolic BloodPressureSystolic Height Weight HeadBMC HeadBMD BonePC1 BonePC2\n');
for i=1:size(ALL_IDPs,1)
  fprintf(po,'%d %f %f %f %f %f %f %f %f \n',ALL_IDPs(i,1),BPd(i),BPs(i),Height(i),Weight(i),BONE_HeadBMC(i),BONE_HeadBMD(i),BONE_PC1(i),BONE_PC2(i));
end
fclose(po);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ICA on netmat edges

grot=ALL_IDPs_i_KG_deconf_GPCdeconf(:,(end+1-size(NET25,2)-size(NET100,2)):end);
NK=~isnan(sum(grot,2)); NKN=sum(NK)
NETs_i_deconf_NK=grot(NK,:);
% split-half randomly or by sex
grotN1=1:2:NKN; grotN2=2:2:NKN;
[NpcaU,NpcaS,NpcaV]   =nets_svds(NETs_i_deconf_NK,0);
[NpcaU1,NpcaS1,NpcaV1]=nets_svds(NETs_i_deconf_NK(grotN1,:),0);
[NpcaU2,NpcaS2,NpcaV2]=nets_svds(NETs_i_deconf_NK(grotN2,:),0);

%%% first-pass: compare a range of dimensionalities
%for J=20:-1:10, grot2=[]; figure; for JJ=1:5
%%% second-pass: find best ICA run for chosen dimensionality
J=14; JJJ=6; grot2=[]; for JJ=1:1000 % J=dimensionality input to ICA, JJJ=how many we aim to keep

[icaSd,icaAd,icaWd]   =fastica(NpcaV(:,1:J)','approach','symm','g','tanh','epsilon',1e-13,'maxNumIterations',3000,'lastEig',J);
[icaSd1,icaAd1,icaWd1]=fastica(NpcaV1(:,1:J)','approach','symm','g','tanh','epsilon',1e-13,'maxNumIterations',3000,'lastEig',J);
[icaSd2,icaAd2,icaWd2]=fastica(NpcaV2(:,1:J)','approach','symm','g','tanh','epsilon',1e-13,'maxNumIterations',3000,'lastEig',J);

[icaSdb,~,icaSd1b,~,icaSd2b,~] = ssorder2(icaSd',icaAd',icaSd1',icaAd1',icaSd2',icaAd2');
grot=corr([icaSdb icaSd1b icaSd2b]); % figure; imagesc(grot,[-1 1]); figure;plot(diag(grot(J+1:2*J,2*J+1:end))');

%%% first-pass
%grot2=[grot2;diag(grot(J+1:2*J,2*J+1:end))'];end;plot(grot2'+randn(size(grot2'))*0.01);hold on;plot(mean(grot2)','k','LineWidth',2); title(sprintf('J=%d',J));end;
%%% second-pass
grot2=[grot2 mean(diag(grot(J+1:J+JJJ,2*J+1:2*J+JJJ)))]
if grot2(end)>.929, break; end; end; icaSdb=icaSdb(:,1:JJJ)';

grot = nets_normalise(NETs_i_deconf_NK * icaSdb');  
icaU=zeros(size(KG,1),size(grot,2))/0; icaU(NK,:)=grot;
NETICA_i_KG_deconf_GPCdeconf=icaU;
NETICA_i_KG_deconf_GPCdeconf_BODYdeconf=nets_unconfound(icaU, [conf(KG,:) GPCs BodyConf3]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

