function a=find_score_coef(param,IHD)

%a=zeros(4,1); % optimal polynomial coefficients


% % % param=[0.9 1 1.1 1.05 1.5 2]';
% % % IHD=[0 0 0 1 1 1]';

% p1=randn(100,1)+5;
% p2=randn(100,1)+2;
% param=[p1;p2];
% IHD=[ones(100,1);zeros(100,1)];

x1=param(IHD==1 & ~isnan(param)); % died
x2=param(IHD==0 & ~isnan(param)); % survved

N1=length(x1);
N2=length(x2);

X1=[x1.^3 x1.^2 x1.^1 ones(N1,1)];
X2=[x2.^3 x2.^2 x2.^1 ones(N2,1)];


X1=[x1.^2 x1.^1 ones(N1,1)];
X2=[x2.^2 x2.^1 ones(N2,1)];

s1=(sum(X1,1))';
s2=(sum(X2,1))';

S1=repmat(s1',N1,1);
S2=repmat(s2',N2,1);


% f     = @(x) (1-x(1)).^2 + 100*(x(2)-x(1).^2).^2;
% gradf = @(x) [-2*(1-x(1))-200*x(1).*(x(2)-x(1).^2);200*(x(2)-x(1).^2)];
 a0   = [0;0;.01;0];
 a0   = [.01; 0.1;0];
 
 
 
a= cg (X1,X2,a0,s1,s2,S1,S2);

figure(1)
scatter(x1,ones(size(x1)),'r*')
hold on
scatter(x2,zeros(size(x2)),'bo')

p1_trans=X1*a
p2_trans=X2*a
figure(2)
scatter(p1_trans,ones(size(p1_trans)),'r*')
hold on
scatter(p2_trans,ones(size(p2_trans)),'bo')
end