function DATA=MD_classification(IHD,time_series_names,MD_DATA_D,MD_DATA_DD,MD_DATA,DESCRIPTORS,des_names,R,MD_VARI)
%for par_idx=1:num_params
%X=[MD_DATA_D MD_DATA_DD MD_DATA];
min_numdata1=min(sum(~isnan(MD_DATA_D),2));
min_numdata2=min(sum(~isnan(MD_DATA_DD),2));
MD_DATA_D=MD_DATA_D(:,1:min_numdata1);
MD_DATA_DD=MD_DATA_DD(:,1:min_numdata2);
X=[MD_DATA_D MD_DATA_DD MD_DATA MD_VARI];
%X=MD_VARI;
AGE_a=DESCRIPTORS(:,strcmp(des_names,'Age'));
Gender_a=DESCRIPTORS(:,strcmp(des_names,'Gender'));
ICUTYPE_a=[DESCRIPTORS(:,strcmp(des_names,'ICUType'))];
ICUTYPE_t=ICUTYPE_a;
ICUTYPE_a(ICUTYPE_t==1)=2;
ICUTYPE_a(ICUTYPE_t==2)=4;
ICUTYPE_a(ICUTYPE_t==3)=1;
ICUTYPE_a(ICUTYPE_t==4)=3;
X=[AGE_a ICUTYPE_a MD_DATA MD_VARI MD_DATA_D MD_DATA_DD];
%0.2708
%X=[R(:,~strcmp(time_series_names,'MechVent'))];
%0.277
% X=[R(:,strcmp(time_series_names,'TroponinT'))];
% X=[R(:,strcmp(time_series_names,'TroponinI')) X];
% X=[R(:,strcmp(time_series_names,'RespRate')) X];
% X=[R(:,strcmp(time_series_names,'Mg')) X];
% X=[R(:,strcmp(time_series_names,'Lactate')) X];
% X=[R(:,strcmp(time_series_names,'K')) X];
% X=[R(:,strcmp(time_series_names,'FiO2')) X];
% X=[R(:,strcmp(time_series_names,'Bilirubin')) X];
% X=[R(:,strcmp(time_series_names,'AST')) X];
% X=[R(:,strcmp(time_series_names,'ALT')) X];
% X=[R(:,strcmp(time_series_names,'ALP')) X];
% X=[R(:,strcmp(time_series_names,'Albumin')) X];
%
%
% max_N=size(MD_DATA_D,2);
% %X=X(:,1:max_N);
%
% %X_orig1=MD_DATA_D;
%
% scores=zeros(max_N,1);
%
% for idx=10:max_N
%
% X=[AGE_a ICUTYPE_a MD_DATA MD_VARI MD_DATA_D(:,1:idx) MD_DATA_DD(:,1:idx)];
% %[best_th0,max_score_0,D] = opt_th(X,IHD);
% %
% % [best_th1,max_score_1,D0] = opt_th(X(Gender_a==0,:),IHD(Gender_a==0));
% % [best_th2,max_score_2,D1] = opt_th(X(Gender_a==1,:),IHD(Gender_a==1));
%
% [best_th0,max_score_0,D] = opt_th_valid(X,IHD);
% % [best_th1,max_score_1,D0] = opt_th_valid(X(Gender_a==0,:),IHD(Gender_a==0));
% % [best_th2,max_score_2,D1] = opt_th_valid(X(Gender_a==1,:),IHD(Gender_a==1));
%
% scores(idx,:)=[max_score_0];
% end
%
% figure(1)
% plot(scores(:,1),'g')
% % hold on
% % plot(scores(:,2),'r')
% % plot(scores(:,3),'b')
% pause(.1)
[best_th0,max_score_0,D] = opt_th(X,IHD);
[best_th1,max_score_1,D0] = opt_th(X(Gender_a==0,:),IHD(Gender_a==0));
[best_th2,max_score_2,D1] = opt_th(X(Gender_a==1,:),IHD(Gender_a==1));
DATA(:,1)=0;
DATA(Gender_a==0,2)=D0(:,2);
DATA(Gender_a==1,2)=D1(:,2);
DATA(Gender_a==0,3)=D0(:,3);
DATA(Gender_a==1,3)=D1(:,3);
function [best_th,max_score1,BEST_DATA]=opt_th(X_,IHD_)
num_params=size(time_series_names,1);
DATA=zeros(size(IHD_,1),3);
B = mnrfit(X_,IHD_+1);
PHAT = mnrval(B,X_);
max_score1=0;
max_score2=0;
best_th=0;
for class_th=.1:.01:0.5
DATA(:,1)=str2double('0000');
DATA(:,2)=PHAT(:,2);
DATA(:,3)=PHAT(:,2)> class_th;
DATA(DATA(:,2)<0.01,2)=0.01;
DATA(DATA(:,2)>0.99,2)=0.99;
% if(~isempty(results))
% Calculate sensitivity (Se) and positive predictivity (PPV)
TP=sum(DATA(IHD_==1,3));
FN=sum(~DATA(IHD_==1,3));
FP=sum(DATA(IHD_==0,3));
Se=TP/(TP+FN);
PPV=TP/(TP+FP);
show=0; % if show is 1, the decile graph will be displayed by lemeshow()
H=lemeshow([IHD_ DATA(:,2)],show);
% Use the title of figure to display the results
title(['H= ' num2str(H) ' Se= ' num2str(Se) ' PPV= ' num2str(PPV) '. ' num2str(class_th) ])
% The event 1 score is the smaller of Se and PPV.
score1 = min(Se, PPV);
if score1>max_score1
max_score1=score1;
best_th=class_th;
max_score2=H;
BEST_DATA=DATA;
% display(['Unofficial Event 1 score: ' num2str(score1)]);
end
% end
end
% survivals1=sum(X==1 & IHD==0);
% deaths1=sum(X==1 & IHD==1);
% survivals2=sum(X==0 & IHD==0);
% deaths2=sum(X==0 & IHD==1);
%
% disp([' non NaNs:' num2str(deaths1/(deaths1+survivals1)) ])
% disp([' NaNs:' num2str(deaths2/(deaths2+survivals2)) ])
%
max_score1
end
function [best_th,max_score1,BEST_DATA]=opt_th_valid(X_,IHD_)
% num_params=size(time_series_names,1);
n=size(X_,1);
s=zeros(n,1);
s(1:round(n/2))=1;
s=boolean(s);
B = mnrfit(X_(s,:),IHD_(s)+1);
PHAT = mnrval(B,X_(~s,:));
BEST_DATA=[];
DATA_=zeros(size(PHAT,1),3);
max_score1=0;
max_score2=0;
best_th=0;
for class_th=.1:.01:0.5
DATA_(:,1)=str2double('0000');
DATA_(:,2)=PHAT(:,2);
DATA_(:,3)=PHAT(:,2)> class_th;
DATA_(DATA_(:,2)<0.01,2)=0.01;
DATA_(DATA_(:,2)>0.99,2)=0.99;
% if(~isempty(results))
% Calculate sensitivity (Se) and positive predictivity (PPV)
TP=sum(DATA_(IHD_(~s)==1,3));
FN=sum(~DATA_(IHD_(~s)==1,3));
FP=sum(DATA_(IHD_(~s)==0,3));
Se=TP/(TP+FN);
PPV=TP/(TP+FP);
% show=0; % if show is 1, the decile graph will be displayed by lemeshow()
% H=lemeshow([IHD_ DATA_(:,2)],show);
% Use the title of figure to display the results
% title(['H= ' num2str(H) ' Se= ' num2str(Se) ' PPV= ' num2str(PPV) '. ' num2str(class_th) ])
% The event 1 score is the smaller of Se and PPV.
score1 = min(Se, PPV);
if score1>max_score1
max_score1=score1;
best_th=class_th;
% max_score2=H;
BEST_DATA=DATA_;
% display(['Unofficial Event 1 score: ' num2str(score1)]);
end
% end
end
% survivals1=sum(X==1 & IHD==0);
% deaths1=sum(X==1 & IHD==1);
% survivals2=sum(X==0 & IHD==0);
% deaths2=sum(X==0 & IHD==1);
%
% disp([' non NaNs:' num2str(deaths1/(deaths1+survivals1)) ])
% disp([' NaNs:' num2str(deaths2/(deaths2+survivals2)) ])
%
%
% max_score1
end
end