function [MDL0ErrorDim2, mleSmallError, MDL0ErrorDim, mleBigError, NML0Error1, NML0Error, MDL1dim, NML, NML2, NML0Errorpart, NML0Errorpart1, MDL2dim] = NMLexampleAnalysis(P,x,n)

%Clintin P. Davis-Stober
%University of Missouri at Columbia, MO
%This function accepts a set of points, P, and a data vector, d, of
%relative frequencies and returns the unique best (if exists), by MDL criteria, set of points
%and their mixing weights.
%ALGORITHM CURRENT SAMPLES ONCE FROM THE UNIT CUBE!


tic

load FullSpace.mat;

p = size(P,2);
n1 = size(P,1);
MDL0dim = zeros(n1,1);
MDL0dimError = zeros(n1,1);
MDL1LML = zeros(n1*(n1-1)/2,1);
mle1dim = zeros(n1*(n1-1)/2, p);
numSam = 961;

Samples = unidrnd(n(1),numSam,p);


%Compute MDL for the high error 0-dimensional representations

NumDev = zeros(n1,1);
thet = zeros(n1,1);
NML0Error = zeros(n1,1);
mleBigError = zeros(n1,p);
for i = 1:n1
    NumDev(i) = sum(abs(n'.*P(i,:) - x'));
    thet(i) = NumDev(i)/sum(n);
    if thet(i) >= .20
        thet(i) =.20;
    elseif thet(i) <= .01
         thet(i) = .01;
    end
    
    MDL0dimError(i) = NegLike(P(i,:)*(1-2*thet(i)) + ones(1,p)*thet(i),x,n,p);
    mleBigError(i,:) = P(i,:)*(1-2*thet(i)) + ones(1,p)*thet(i);
    NML0Errorpart = zeros(numSam,1);
    theta = zeros(961,1);
    for j = 1:961
        NumDev2 = sum(abs(n'.*P(i,:) - U(j,:)));
        theta(j) = NumDev2/sum(n);
        if theta(j) >=.20
            theta(j) =.20;
        elseif theta(j) <= .01
         theta(j) = .01;    
        end
        NML0Errorpart(j) = exp(NegLike(P(i,:)*(1-2*theta(j)) + ones(1,p)*theta(j), U(j,:)',n,p));
    end
    NML0Error(i) = sum(NML0Errorpart);
        
end

    
MDL0ErrorDim = MDL0dimError - log(NML0Error);
    

%Compute MDL for the low error 0-dimensional representations

MDL0dimError1 = zeros(n1,1);
NumDev = zeros(n1,1);
NML0Error1 = zeros(n1,1);
thet1 = zeros(n1,1);
%NML0Error = zeros(n1,1);
mleSmallError = zeros(n1,p);
for i = 1:n1
    NumDev(i) = sum(abs(n'.*P(i,:) - x'));
    thet1(i) = NumDev(i)/sum(n);
    if thet1(i) >= .03
        thet1(i) =.03;
    elseif thet1(i) <= .01
         thet1(i) = .01;
    end
    MDL0dimError1(i) = NegLike(P(i,:)*(1-2*thet1(i)) + ones(1,p)*thet1(i),x,n,p);
    mleSmallError(i,:) = P(i,:)*(1-2*thet1(i)) + ones(1,p)*thet1(i);
    NML0Errorpart1 = zeros(numSam,1);
    theta1 = zeros(961,1);
    for j = 1:961
        NumDev2 = sum(abs(n'.*P(i,:) - U(j,:)));
        theta1(j) = NumDev2/sum(n);
        if theta1(j) >=.03
            theta1(j) =.03;
        elseif theta1(j) <= .01
         theta1(j) = .01;
        end
        NML0Errorpart1(j) = exp(NegLike(P(i,:)*(1-2*theta1(j)) + ones(1,p)*theta1(j), U(j,:)',n,p));
    end
    NML0Error1(i) = sum(NML0Errorpart1);
        
end
    
MDL0ErrorDim2 = MDL0dimError1 - log(NML0Error1);
 

%


%Compute MDL for all 1-dimensional representations

cc = 1:n1;
c = nchoosek(cc,2);

tol=1e-10;
 options=optimset('Algorithm','interior-point', 'TolX', tol, 'tolcon', .01);

 NML = zeros((n1*(n1-1))/2,1);
 y = c;
 %Samples = unidrnd(n(1),numSam,p);
 
parfor i = 1:(n1*(n1-1)/2)
    [mle, fval,j,z, x0, A, b] = maxLike([P(c(i,1),:); P(y(i,2),:)],x,n,p);
    MDL1LML(i) = -fval; 
    mle1dim(i,:) = mle;
    %Samples = unidrnd(n(1),numSam,p);
    NMLpart = zeros(961,1);
    for j=1:961
        [mle,fval,exitflag,output,lambda]=fmincon(@(theta) -NegLike(theta,U(j,:)',n,p), ...
    x0, [], [],A,b,ones(p,1)*.0001, ones(p,1)*.9999, [], options);
        NMLpart(j) = exp(-fval); 
    end
    NML(i) = sum(NMLpart);    

end
    
MDL1dim = MDL1LML - log(NML);



%Compute NML for the three-state mixture case



xnew = .3333*P(1,:) + .3333*P(2,:) + .3333*P(3,:);
[mle2,fval,exitflag,output,lambda]=fmincon(@(thetag) -NegLike(thetag,x,n,p), ...
    xnew,[-1 0; 1 -1; 0 1], [0;0;1],[],[],ones(p,1)*.0001, ones(p,1)*.9999, [], options);

 MDL2LML = exp(-fval);
 NML2part = zeros(numSam,1);
 
 parfor j=1:961
     [mle,fval,exitflag,output,lambda]=fmincon(@(thetag) -NegLike(thetag, U(j,:)',n,p), ...
    xnew, [-1 0; 1 -1; 0 1], [0;0;1],[],[],ones(p,1)*.0001, ones(p,1)*.9999, [], options);
      NML2part(j) = exp(-fval);
 end
NML2 = sum(NML2part);    

MDL2dim = log(MDL2LML) - log(NML2);
 



toc

end


function [fnegL] = NegLike(theta, x, n, p)
negll = zeros(p,1);

for i=1:p
    negll(i) = LogLike(theta(i), x(i), n(i));     
end
fnegL = sum(negll);
end


function [Neg] = LogLike(theta, x,n)

Neg = x*log(theta) + (n-x)*log(1-theta);


end

function [A,b,x0,j,z] = liney(P)


p = size(P,2);

A = eye(p,p);

j = P(1,:);
k = P(2,:);
z = k-j;
Q = zeros(p,2);

b = zeros(p,1);

for i = 1:p
    Q(i,1) = i;
    Q(i,2) = z(i);
end

Q = sortrows(Q,2);
if Q(p,2) == 0
    Q = sortrows(Q,-2);
end

w = zeros(1,2);
w = w + Q(p,:);

for i = 1:p
   b(i) = j(i) - (j(w(1,1))*z(i))/(z(w(1,1)));
end

for i = 1:p
    A(i,w(1,1)) = -(z(i))/z(w(1,1));
end
A(w(1,1),w(1,1)) = 0;

x0 = j + .5*z;

end

function [mle, fval,j,z, x0, A, b] = maxLike(P, x,n,p)

tol=1e-10;
 options=optimset('Algorithm','interior-point', 'TolX', tol, 'tolcon', .01);


[A,b,x0,j,z] = liney([P(1,:); P(2,:)]);

 
f = rand(size(x0)); % assumes x0 is the initial point
xnew = linprog(f,[],[],A,b,zeros(p,1),ones(p,1));


[mle,fval,exitflag,output,lambda]=fmincon(@(theta) -NegLike(theta,x,n,p), ...
    x0, [], [],A,b,ones(p,1)*.0001, ones(p,1)*.9999, [], options);


end



