res=[]; res50=[]; 
global X; global Y; 
rand('seed',1); randn('seed',1);


load spec51_matrix.txt
K=spec51_matrix;  clear spec51_matrix;
k2=sqrt(diag(K))+1e-12; k2=k2*k2'; K=K ./ k2; clear k2; % Normalize matrix
K=spec; clear K;  % need to save copy of spec kernel
load blast_matrix.txt; hom=sparse(blast_matrix(:,1:end-1));  clear blast_matrix;
Ys=load('allseq_membership-noheaders.table'); % load data splits 
name='spectrum_kernel_profile[blast]'; % name of file to save
global K;
 

for i=1:size(Ys,2) %% loop over data splits

Y=Ys(:,i);  %% allocate correct train and test set labels for this data split
l=sum(Y==1)+sum(2==Y); lt=sum(Y==3)+sum(4==Y); % test & test set sizes
[trn]=find(Y==1 | Y==2); [tst]=find(Y==3 | Y==4);
Y(trn(1:l))=sign((Y(trn(1:l))==1)-0.5); Y(tst(1:lt))=sign((Y(tst(1:lt))==3)-0.5);
%%--------------- pull in domain homologs -----------
load blast_homologs addpts; f=addpts{i};  % already cached these computations
trn=[trn ;f']; trn=unique(trn); Y(f)=1;
%% ------------------------------------------------
trnhom = hom;
trnhom(:,tst) = 0;
trnhom = max(trnhom, sparse(diag(ones(7329,1))));
numhom = sum(trnhom,2);
trnhom = trnhom ./ repmat(numhom,1,7329);
unlab = ones(1,7329);
unlab(1,[trn,tst]) = 0;  %% zero out the rows for the unlabelled points
trnhom(find(unlab==1),:) = 0;
sparse(trnhom);
clear numhom; clear unlab;
K = sparse(trnhom * spec * trnhom');
clear trnhom;


X=Y'; Y=Y'; d=data_global(X,Y); dtrn=get(d,trn); dtst=get(d,tst);  
a=svm; a.optimizer='svmlight'; a.child=kernel('custom_fast'); 
a.algorithm.use_signed_output=0;
a.balanced_ridge=0.02*median(diag(K));
[tr a]=train(a,dtrn); r=test(a,dtst);
[dummy ind]=sort(-r.X); rr=r.Y(ind,:); rr(rr==-1)=0; 
 
clear K;

res=[res roc(rr)]
res50=[res50 roc(rr,50)]
save(name,'res','res50');

end 








