function [pMaxEnt, Lambda, paramLocalConst, paramDebug]=EMK_mod2(s,varargin)

%% Call internal test function when no arguments supplied
if nargin==0
    help EMK
    test_EMK
    return
end

showCost=0;
lowBound=1e-50;

lBound=min(-5, 3*min(s));
rBound=max(5, 3*max(s));
numIntBin=2^11;
intBin=linspace(lBound,rBound,numIntBin);
global rIntMat dIntBin;
rIntMat=[];
dIntBin=intBin(2)-intBin(1);
%% Gather Options

% build default options structure
Params=struct( ...
    'typeGlobalConst',[0 0 0 1], ... % list of types of global constraints
    'typeLocalConst',2, ... % 1:Parzen window    2:Gaussian kernel
    'numLocalConst',Inf, ... % number of local constraints
    'paramLocalConst', [],  ... % Parameters for local constraints, must has same number of rows as numLocalConst.
    ... % Parzen window:[min max]  Gaussian kernel:[mean std]
    'greedySearch',1 ... % greedy search for better result
    );

% load in user supplied options
Params=getopt(Params,varargin{:});

[K T]=size(s);  %sample size
typeLocalConst=Params.typeLocalConst;
numLocalConst=Params.numLocalConst;
paramLocalConst=Params.paramLocalConst;
numGlobalConst=4;          %number of global constriants
typeGlobalConst=Params.typeGlobalConst;               %
M=numGlobalConst; %number of the total constraints
greedySearch=Params.greedySearch;

if greedySearch==1
%     numBin=round((log10(T)-0.6)*160);
    numBin=round((log10(T))*80);
%     numBin=2^10;
else
    numBin=round((log10(T)-1.4)*40);
end

r(1).func=@(x) ones(size(x));
r(1).diff=@(x) zeros(size(x));
r(1).int=@(x) x;
rIntMat(1,:)=r(1).func(intBin);

r(2).func=@(x) x;
r(2).diff=@(x) ones(size(x));
r(2).int=@(x) (x.^2)/2;
rIntMat(2,:)=r(2).func(intBin);

r(3).func=@(x) x.^2;
r(3).diff=@(x) 2*x;
r(3).int=@(x) (x.^3)/3;
rIntMat(3,:)=r(3).func(intBin);

%The default is x/1+x^2
r(4).func=@(x) x./(1+x.^2);
r(4).diff=@(x) (1-x.^2)./((1+x.^2).^2);
r(4).int=@(x) 0.5*log(1+x.^2);
rIntMat(4,:)=r(4).func(intBin);


%Four types of constrain functions
rG(1).func=@(x) x.^4;
rG(1).diff=@(x) 4*x.^3;
rG(1).int=@(x) (x.^5)/5;

rG(2).func=@(x) abs(x)./(1+abs(x));
rG(2).diff=@(x) sign(x)./((1+abs(x)).^2);
rG(2).int=@(x) sign(x).*abs(x).*log(1+abs(x))+(abs(x)+1).*log(abs(x)+1)-abs(x);

rG(3).func=@(x) x.*abs(x)./(10+abs(x));
rG(3).diff=@(x) abs(x).*(20+abs(x))./((10+abs(x)).^2);
rG(3).int=@(x) 0.5*(-20*abs(x)+x.^2+200*log(10+abs(x)));

rG(4).func=@(x) x./(1+x.^2);
rG(4).diff=@(x) (1-x.^2)./((1+x.^2).^2);
rG(4).int=@(x) 0.5*log(1+x.^2);

stabFlag=zeros(1,size(typeGlobalConst,2));
entGlobalConst=Inf*ones(1,size(typeGlobalConst,2));
for i=1:size(typeGlobalConst,2)
    if typeGlobalConst(i)==0
        continue;
    end
    %Choose the fourth constraint according to typeGlobalConst
    r(4).func=rG(i).func;
    r(4).diff=rG(i).diff;
    r(4).int=rG(i).int;
    rIntMat(4,:)=r(4).func(intBin);
    
    p=@(x,Lambda) exp(-1+Lambda(1)+Lambda(2:end)'*rFuncs(r, size(Lambda,1), x));
    
    [LambdaTildeLE, Alpha(:,i)] = maxEntropyDensityEstLE( s, numGlobalConst, r, p, lBound, rBound );
    if showCost
        showDebugInfo(['Global', num2str(i),' LE'], LambdaTildeLE, Alpha(:,i), 1, []);
    end
    [LambdaGlobal(:,i), stabFlag(i), fConstCostTilde] = maxEntropyDensityEstNI...
        ( numGlobalConst, r, p, LambdaTildeLE, Alpha(:,i), i, lBound, rBound );
    entGlobalConst(i)=1-LambdaGlobal(:,i)'*Alpha(:,i);
    if showCost
        showDebugInfo(['Global', num2str(i)], LambdaGlobal(:,i), Alpha(:,i), stabFlag(i), fConstCostTilde);
    end
end
% if sum(stabFlag)==0
%     disp('Global constraints cannot converge1!');
%     typeGlobalConst=0;
%     Lambda=LambdaTildeLE;
%     LambdaTilde{1}=LambdaTildeLE;
%     return;
% end
minEnt=min(entGlobalConst(stabFlag==1));
typeGlobalConst=find(entGlobalConst==minEnt);
%Choose the fourth constraint according to typeGlobalConst
r(4).func=rG(typeGlobalConst).func;
r(4).diff=rG(typeGlobalConst).diff;
r(4).int=rG(typeGlobalConst).int;
rIntMat(4,:)=r(4).func(intBin);

p=@(x,Lambda) exp(-1+Lambda(1)+Lambda(2:end)'*rFuncs(r, size(Lambda,1), x));
LambdaTilde{1}=LambdaGlobal(:,typeGlobalConst);
Alpha=Alpha(:,typeGlobalConst);

%Local constraints
if greedySearch==1
    [~,kernelP,kernelT,~]=kde(s, numBin);
    kernelP=kernelP';
    kernelP(kernelP<lowBound)=lowBound;
    binWidth=kernelT(2)-kernelT(1);
    numBin=size(kernelT,2);
else
    [kernelP, kernelT]=hist(s,numBin);
    binWidth=kernelT(2)-kernelT(1);
    kernelP=kernelP/binWidth/T;
    kernelP(kernelP<lowBound)=lowBound;
end

Lambda=LambdaTilde{1};
logLikelihood(1)=-sum(log(p(s,Lambda)))/T;
Penalty(1)=(numGlobalConst)*K*log(T)/2/T;
MDL(1)=logLikelihood(1)+Penalty(1);

KernelGMM=4;
kernelPos=[];
bandWidth=[];
i=1;
while numLocalConst>0
    if size(paramLocalConst,1)>=i
        kernelPos(i)=paramLocalConst(i,1);
        bandWidth(i)=paramLocalConst(i,2);
        greedyKernelPos=paramLocalConst(i,1);
        greedyBandWidth=paramLocalConst(i,2);
        greedyDegree=1;
    else
        pdfTilde = @(x) p(x,Lambda);
        pdfP=pdfTilde(kernelT);
        pdfP(pdfP<lowBound)=lowBound;
%         diffP=kernelP-pdfP;            %Difference in terms of absolute value
%         diffP=abs(kernelP*log(kernelP/pdfP)).*(kernelP-pdfP); %Difference in terms of KL divergence
        diffP=(abs(kernelP.*log(kernelP./pdfP))+abs(pdfP.*log(pdfP./kernelP))).*(kernelP-pdfP); %This is what Gengshen used as a default diffP
        [volDiff, peakPos, peakWidth]=findDiffPDF(diffP, K, numBin);
        [volDiff, idxVolDiff]=sort(abs(volDiff), 'descend');
        peakPos=peakPos(idxVolDiff);
        peakWidth=peakWidth(idxVolDiff);
        if greedySearch==1
            greedyDegree=min(3,length(peakPos));
            greedyKernelPos=kernelT(peakPos(1:greedyDegree));
            greedyBandWidth=binWidth*(1+peakWidth(1:greedyDegree));
            
            numKernel=min(KernelGMM,length(peakPos));
            model.mu=[kernelT(peakPos(1:numKernel))];
            for idxKernel=1:numKernel
                model.Sigma(:,:,idxKernel)=(binWidth*(1+peakWidth(idxKernel)))^2;
            end
            model.weight = ones(1,numKernel)/numKernel;
            [model, ~] = localParamEstEM(s, log(pdfTilde(s)'), model);
            tmpKernelPos=model.mu(1:end-1);
            tmpBandWidth=sqrt(model.Sigma(1:end-1));
            tmpWeight=model.weight(1:end-1);
            [tmpWeight, idxWeight]=sort(abs(tmpWeight), 'descend');
            tmpKernelPos=tmpKernelPos(idxWeight);
            tmpBandWidth=tmpBandWidth(idxWeight);
            greedyDegree=min(2,length(tmpWeight));
            greedyKernelPos=[greedyKernelPos tmpKernelPos(1:greedyDegree)];
            greedyBandWidth=[greedyBandWidth squeeze(tmpBandWidth(:,:,1:greedyDegree))'];
        else
            greedyDegree=min(1,length(peakPos));
            greedyKernelPos=kernelT(peakPos(1:greedyDegree));
            greedyBandWidth=binWidth*(1+peakWidth(1:greedyDegree));
        end
        greedyKernelPos=greedyKernelPos((greedyBandWidth>0.1)&(greedyBandWidth<0.9));
        greedyBandWidth=greedyBandWidth((greedyBandWidth>0.1)&(greedyBandWidth<0.9));
        greedyDegree=length(greedyKernelPos);
    end
%     figure; plot(kernelT,diffP*100, 'r', kernelT, kernelP-pdfP, 'b', kernelT, kernelP.*log(kernelP./pdfP)*10, 'g');
    
    M=M+1;
    LambdaOld=Lambda;
    greedylogLikelihood=[];
    greedyLambda=[];
    for greedyIdx=1:greedyDegree
        kernelPos(i)=greedyKernelPos(greedyIdx);
        bandWidth(i)=greedyBandWidth(greedyIdx);
        if prod(abs(kernelPos(1:i-1)-kernelPos(i))+abs(bandWidth(1:i-1)-bandWidth(i)))==0
            %disp('Picked an existing kernel!');
            greedylogLikelihood(greedyIdx)=inf;
            continue;
        end
        if typeLocalConst==1
            %step function
            r(numGlobalConst+i).func=@(x) prod(double((x-(kernelPos(i)-bandWidth(i))).*((kernelPos(i)+bandWidth(i))-x)>=0),1);
        else
            %Gaussian kernel
            r(numGlobalConst+i).func=@(x) exp(-(x-kernelPos(i)).^2/2/(bandWidth(i).^2));
        end
        rIntMat(numGlobalConst+i,:)=r(numGlobalConst+i).func(intBin);
        p=@(x,Lambda) exp(-1+Lambda(1)+Lambda(2:end)'*rFuncs(r, size(Lambda,1), x));
        Lambda(M)=0;
        greedyAlpha(greedyIdx)=sum(r(numGlobalConst+i).func(s))/T;
        Alpha(M)=greedyAlpha(greedyIdx);
        [greedyLambda(:,greedyIdx), stabFlag, fConstCost] = maxEntropyDensityEstNI...
            ( M, r, p, Lambda, Alpha, typeGlobalConst, lBound, rBound );
        if stabFlag==1
            greedylogLikelihood(greedyIdx)=-sum(log(p(s,greedyLambda(:,greedyIdx))))/T;
        else
            greedylogLikelihood(greedyIdx)=inf;
        end
    end
    if prod(size(greedyLambda))==0
        %disp('Picked an existing kernel!');
        break;
    else
        [tmpLogLikelihood, idxLikelihood]=min(greedylogLikelihood);
        if tmpLogLikelihood==inf
            disp('Cannot find better solution!');
            break;
        end
        logLikelihood(i+1)=tmpLogLikelihood;
        Alpha(M)=greedyAlpha(idxLikelihood);
        Lambda=greedyLambda(:,idxLikelihood);
        kernelPos(i)=greedyKernelPos(idxLikelihood);
        bandWidth(i)=greedyBandWidth(idxLikelihood);
        if typeLocalConst==1
            %step function
            r(numGlobalConst+i).func=@(x) prod(double((x-(kernelPos(i)-bandWidth(i))).*((kernelPos(i)+bandWidth(i))-x)>=0),1);
        else
            %Gaussian kernel
            r(numGlobalConst+i).func=@(x) exp(-(x-kernelPos(i)).^2/2/(bandWidth(i).^2));
        end
        rIntMat(numGlobalConst+i,:)=r(numGlobalConst+i).func(intBin);
        p=@(x,Lambda) exp(-1+Lambda(1)+Lambda(2:end)'*rFuncs(r, size(Lambda,1), x));
    end
    
    
    
    if showCost
        showDebugInfo(['Local', num2str(i)], Lambda, Alpha, stabFlag, fConstCost);
    end
    
    if stabFlag==0
        disp('Global constraints cannot converge2!');
        Lambda=LambdaOld;
        M=M-1;
        break;
    end
    
    Penalty(i+1)=(numGlobalConst+3*i)*K*log(T)/2/T;
    MDL(i+1)=logLikelihood(i+1)+0.5*Penalty(i+1); % For all my experiments I had 0.5.  
    if MDL(i+1)>MDL(i)
        %disp('Order selected by ITC!');
        Lambda=LambdaOld;
        M=M-1;
        break;
    end
    numLocalConst=numLocalConst-1;
    i=i+1;
    LambdaTilde{i}=Lambda;
end

paramLocalConst={[kernelPos' bandWidth']};
if showCost
    [numGlobalConst typeGlobalConst numLocalConst typeLocalConst]
end

pMaxEnt=@(x) p(x,Lambda);
M=length(Lambda);

paramDebug={p, LambdaTilde, M, typeGlobalConst, typeLocalConst};
clear rIntMat dIntBin;
i=1;


function test_EMK
clc
clear
close all
dbstop if error;
dbstop if warning;
%% built-in test function for MaxEntDist
disp('Executing internal test function for MaxEntDist.')

%Loading functions:
dataPath=[pwd, '/Results/Latest/'];
%dirPath=[dataPath(1:strfind(pwd,'ICA')-1) 'ICA/'];
%addpath(genpath([dirPath 'utility']));
addpath(genpath([pwd]));

%% Simulate input data
T=1000;  %sample size

% beta=2;  %2:Gaussian
% mu=2;
% alpha=0.35;
% paramGGD=[alpha beta -mu; 1-alpha beta mu];        %alpha beta mu
% paramGGD=[1 1.2 1];        %alpha beta mu
% paramGGD=[0.2 1 -3; 0.5 2 0; 0.3 1 3];        %alpha beta mu  3 peaks
% paramGGD=[0.1 1 -4; 0.2 1 -2; 0.45 2 0; 0.25 1 3];        %alpha beta mu  4 peaks

numMixGGD=ceil(rand*5);
Alpha=rand(numMixGGD,1);
Alpha=Alpha/sum(Alpha);
Beta=2.^(rand(numMixGGD,1)*4-1);      %2:Gaussian
Mu=rand(numMixGGD,1)*14-7;
paramGGD=[Alpha, Beta, Mu];

s = genMixGGD( paramGGD, 1, 1, T );

% load('USC_SIPI/lena512.mat');
imagTest1=imread('USC_SIPI/7.1.09.tiff');
imagTest2=imread('USC_SIPI/boat.512.tiff');
% imagTest=imread('USC_SIPI/4.2.07.tiff');
% imagTest=imread('USC_SIPI/UMBC.jpg');
s1=double(reshape(imagTest1, 1, prod(size(imagTest1))));
s1=s1-mean(s1);
s1=s1/std(s1);
%imshow(imagTest1);
s2=double(reshape(imagTest2, 1, prod(size(imagTest2))));
s2=s2-mean(s2);
s2=s2/std(s2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%profile on;
[pMaxEnt1, Lambda1, paramLocalConst1, paramDebug1]=EMK(s1);
%profile off;

p1=paramDebug1{1};
LambdaTilde1=paramDebug1{2};
M1=paramDebug1{3};
typeGlobalConst1=paramDebug1{4};
typeLocalConst1=paramDebug1{5};

pMixGGD=getTrueMixGGDPDF(paramGGD);

%pdfTilde = @(x) p(x,LambdaTilde{1});
for i=1:length(LambdaTilde1)
    pdfTilde1{i}=@(x) p(x,LambdaTilde1{i});
end
isLocal=prod(size(Lambda1))-prod(size(LambdaTilde1{1}));
if isLocal==0
    pdf1 = pdfTilde1{1};
else
%     pdf = @(x) p(x,Lambda);
    pdf1 = pMaxEnt1;
end
plotDetails=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%profile on;
[pMaxEnt2, Lambda2, paramLocalConst2, paramDebug2]=EMK(s2);
%profile off;

p2=paramDebug2{1};
LambdaTilde2=paramDebug2{2};
M2=paramDebug2{3};
typeGlobalConst2=paramDebug2{4};
typeLocalConst2=paramDebug2{5};

pMixGGD=getTrueMixGGDPDF(paramGGD);

%pdfTilde = @(x) p(x,LambdaTilde{1});
for i=1:length(LambdaTilde2)
    pdfTilde2{i}=@(x) p(x,LambdaTilde2{i});
end
isLocal=prod(size(Lambda2))-prod(size(LambdaTilde2{1}));
if isLocal==0
    pdf2 = pdfTilde2{1};
else
%     pdf = @(x) p(x,Lambda);
    pdf2 = pMaxEnt2;
end

plotPDFs( s1, s2, pMixGGD, pdfTilde1, pdfTilde2, pdf1, pdf2, isLocal, M1, M2, paramLocalConst1, paramLocalConst2, plotDetails, imagTest1, imagTest2 );

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

fileName=[dataPath, 'EMK_M', num2str(M1), 'global', num2str(typeGlobalConst1), 'local', num2str(typeLocalConst1),...
    '_Peak', num2str(size(paramGGD,1)),'TestImages'];
print('-depsc', fileName);
eps2pdf([fileName, '.eps'], '/usr/local/bin/gs');
%eps2pdf([fileName, '.eps'], 'C:\Program Files\gs\gs9.04\bin\gswin64c.exe');
hgsave(fileName);
i=1;



function [Lambda, Alpha] = maxEntropyDensityEstLE( s, M, r, p, lBound, rBound )
T=size(s,2);

B=zeros(M-1,M-1);
for i=1:M-1
    for j=1:M-1
        B(i,j)=(r(i+1).int(s)*(r(j+1).diff(s)'))/T;
    end;
end;

A=zeros(M-1,1);
for i=1:M-1
    A(i)=sum(r(i+1).func(s))/T;
end

Lambda=-inv(B)*A;
Lambda=[0;Lambda];

PDF = @(x) p(x,Lambda);
probInt=quadgk(PDF,lBound,rBound);
Lambda(1)=-log(probInt);
Alpha=[1;A];


function [LambdaNew, stabFlag, fConstCost] = maxEntropyDensityEstNI( M, r, p, Lambda, Alpha, typeGlobalConst, lBound, rBound )
global rIntMat dIntBin;

maxIter=100;
tolerance = 1e-4;
stabTol = 4e-1;

fConstCost=Inf;
LambdaOld=Lambda;
LambdaNew=Lambda;
if typeGlobalConst == 1
    checkIdx=4;
else
    checkIdx=3;
end;
iterStable=1;
LambdaStable=[];
for iter = 1 : maxIter
    if LambdaOld(checkIdx)>0
        disp('Wrong lambda!!!');
        stabFlag=0;
        return;
    end
    PDF=@(x) p(x,LambdaOld);
    fConst=zeros(M,1);
    J=zeros(M,M);
    
    exp_term = exp(-1+LambdaOld'*rIntMat);
    for i=1:M
        fConst(i)=exp_term*(rIntMat(i,:)')*dIntBin-Alpha(i);
%         f=@(x) r(i).func(x).*PDF(x);
%         fConst(i)=quadgk(f,lBound,rBound)-Alpha(i);
    end
    
    %for debug
    fConstCost(iter)=sum(abs(fConst)/sum(abs(Alpha)));
    
    if fConstCost(iter)>4*fConstCost(iterStable)
        beta1=beta1/2;
        LambdaOld=LambdaStable;
    else
        beta1=1;
        iterStable=iter;
        LambdaStable=LambdaOld;
        
        for i=1:M
            for j=i:M
                J(i,j)= exp_term *(rIntMat(i,:).*rIntMat(j,:))';
                J(j,i)=J(i,j);
%                 f=@(x) r(i).func(x).*r(j).func(x).*PDF(x);
%                 J(i,j)=quadgk(f,lBound,rBound);
%                 J(j,i)=J(i,j);
            end;
        end;
        J = J * dIntBin;

        if cond(J)>1e15
            stabFlag=0;
            return;
        end
        LambdaDelta=J\fConst;
    end
    beta=1;
    if LambdaOld(checkIdx)>LambdaDelta(checkIdx)
        beta=0.99*LambdaOld(checkIdx)/LambdaDelta(checkIdx);
    end
    
    LambdaNew=LambdaOld-beta1*beta*LambdaDelta;
    %     if (sum(abs(LambdaNew-LambdaOld))/(LambdaOld'*LambdaOld))<tolerance
    if abs(1-((LambdaNew'*LambdaOld))/(LambdaOld'*LambdaOld))<tolerance
%         stabFlag=fConstCost(iter)<stabTol;
        break;
    end
    LambdaOld=LambdaNew;
end
% stabFlag=fConstCost(iter)<stabTol;
stabFlag=1;


function rFuncX=rFuncs(r, M, x)
rFuncX=[];
for i=2:M
    rFuncX=[rFuncX; r(i).func(x)];
end

function [volDiff, peakPos, peakWidth]=findDiffPDF(diffP, K, numBin)

abs_diffP = abs(diffP);
% find local peaks (logical)
peaksIdxMD=imregionalmax(abs_diffP);
% sort vector of peaks
peaksDiff1D=sort(abs_diffP(peaksIdxMD),'descend');
% matrix of local peaks (values)
peaksDiffMD=peaksIdxMD.*abs_diffP;

% iterate through each of of the peaks in descending order
for i=1:length(peaksDiff1D)
    % logicall identify current peak in vector
    peakIdx=(peaksDiffMD==peaksDiff1D(i));
    volDiffOld=0;
    peakIdxNew=peakIdx;
    % value of current peak
    volDiffNew=diffP(peakIdxNew);
    lag=0;
%     cmdIdx='Idx1';
%     for j=2:K
%         cmdIdx=[cmdIdx, ',Idx', num2str(j)];
%     end
%     cmdGetIdx=['[', cmdIdx, '] = ind2sub(size(peakIdx), find(peakIdx==1));'];
%     eval(cmdGetIdx);
%     eval(['peakIdxArray=[', cmdIdx, '];']);
    Idx1 = find(peakIdx==1);
    peakIdxArray = Idx1;
    while sign(diffP(peakIdx))*(volDiffNew-volDiffOld)>=0
        volDiffOld=volDiffNew;
        peakIdxOld=peakIdxNew;
        
        lag=lag+1;
        peakIdxArrayBound=[peakIdxArray-lag;peakIdxArray+lag];
        peakIdxArrayBoundTest=(peakIdxArrayBound>0).*(numBin>=peakIdxArrayBound);
        if prod(reshape(peakIdxArrayBoundTest,prod(size(peakIdxArrayBoundTest)),1))==0
            break;
        end
        
%         cmdIdx='Idx1-lag:Idx1+lag';
%         for j=2:K
%             cmdIdx=[cmdIdx, ',Idx',num2str(j),'-lag:Idx',num2str(j),'+lag'];
%         end
%         cmdGetIdx=['peakIdxNew(', cmdIdx, ')=1;'];
%         eval(cmdGetIdx);
        
        peakIdxNew_new = peakIdxNew;
        range = Idx1-lag:Idx1+lag;
        peakIdxNew(range) = ones(length(range),1);
        
        volDiffNew=sum(diffP(peakIdxNew));

        indexes = Idx1-lag:Idx1+lag;
        diffP(indexes);

    end
    volDiff(i)=volDiffOld;
    peakPos(i,:)=peakIdxArray;
    peakWidth(i)=lag-1;
end


function x = genMixGau( paramGGD, b, a, T )
N=size(paramGGD, 1);

Beta=paramGGD(:,2);
Mu=paramGGD(:,3);
Alpha=sqrt(gamma(1./Beta)./gamma(3./Beta));

Mix=zeros(N,T);
MixT=rand(1,T);
lowThe=0;
highThe=0;
for i=1:N
    highThe=lowThe+paramGGD(i,1);
    Mix(i,:)=(MixT>=lowThe).*(MixT<highThe);
    lowThe=highThe;
end;

s=zeros(N,T);
for i=1:N
    s(i,:)=randn(1, T);
    s(i,:)=(s(i,:)-mean(s(i,:)))/std(s(i,:))+Mu(i);
end;

s=s.*Mix;
x=sum(s,1);

x=filter(b, a, x);
x=(x-mean(x))/std(x);

% figure;
% hist(x,100);


function plotPDFs( s, pMixGGD, pTilde, p, isLocal, M, paramLocalConst, plotDetails )
%function plotPDFs( s1, s2, pMixGGD, pdfTilde1, pdfTilde2, p1, p2, isLocal, M1, M2, paramLocalConst1,paramLocalConst2, plotDetails,imagTest1, imagTest2 )

MyFontSize = 21;
MyMarkerSize = 16;
MyLineWidth = 4.2;
MyLineWidthThin=1.2;

figure;
T=size(s1,2);
NumBin=40;
subplot(2,2,1)
imshow(imagTest1);

subplot(2,2,2)
imshow(imagTest2);

%for histgram begin
[n1 t1]=hist(s1,NumBin);
dt1=t1(2)-t1(1);
n1=n1/T/dt1;
subplot(2,2,3)
bar(t1,n1,'b');
%for KME begin
hold on;
t=linspace(-3,3,2^14);
subplot(2,2,3)
plot(t, p1(t), 'r', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);xlim([-3 3])
%for KME end
%for histgram end
ylabel('p(x)','FontSize',MyFontSize,'FontName','Arial');
xlabel('x','FontSize',MyFontSize,'FontName','Arial');
set(gca,'FontSize',MyFontSize);
set(gca,'FontName','Arial');


%for KME end
[n2 t2]=hist(s2,NumBin);
dt2=t2(2)-t2(1);
n2=n2/T/dt2;
subplot(2,2,4)
bar(t2,n2,'b');
%for KME begin
hold on;
t=linspace(-3,3,2^14);
subplot(2,2,4)
plot(t, p2(t), 'r', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);xlim([-3 3])
%ylabel('p(x)','FontSize',MyFontSize,'FontName','Arial');
xlabel('x','FontSize',MyFontSize,'FontName','Arial');
set(gca,'FontSize',MyFontSize);
set(gca,'FontName','Arial');

if plotDetails
    numDetails=length(pTilde)-1;
    for i=2:numDetails
        MyColor=[(i-1)/numDetails 1-(i-1)/numDetails 0];
        hold on;
        plot(paramLocalConst{1}(i-1,1), 0,'*','MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth,'color',MyColor);
        hold on;
        plot(t, pTilde{i}(t), '--','MarkerSize', MyMarkerSize, 'linewidth', MyLineWidthThin,'color',MyColor);
    end
    hold on;
    if numDetails>=1
        plot(paramLocalConst{1}(end,1), 0,'r*','MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
    end
end

legendPos='NorthEast';
if abs(t(1))>abs(t(end))
    legendPos='NorthWest';
end
% legend('True PDF', 'GMM', 'KDE-DIFF', 'EMN', 'EMK', 'Location', legendPos)


i=1;



function plotHistPDFs( s, pMixGGD, pTilde, p, isLocal, M, paramLocalConst, plotDetails )

MyFontSize = 21;
MyMarkerSize = 16;
MyLineWidth = 2.8;
MyLineWidthThin=1.2;

figure;
T=size(s,2);
NumBin=200;

[pHist t]=hist(s,NumBin);
pHist=pHist/(t(2)-t(1))/sum(pHist);
%t=linspace(-5,5,NumBin);

[n]=ksdensity(s, t);
    
plot(t, pHist, 'k', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
hold on;
% bar(t,n);
hold on;
[mu, alpha, beta] = ggmle(s);
plot(t, ggpdf(t, mu, alpha, beta), 'g', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
hold on;
plot(t, n, 'b', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
hold on;
plot(t, pTilde{1}(t), 'm', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
if isLocal>0
    hold on;
    plot(t, p(t), 'r', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
end
if plotDetails
    for i=2:length(pTilde)
        MyColor=[1 0 1-(i-1)/length(pTilde)];
        hold on;
        plot(paramLocalConst{1}(i-1,1), 0,'*','MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth,'color',MyColor);
        hold on;
        plot(t, pTilde{i}(t), '--','MarkerSize', MyMarkerSize, 'linewidth', MyLineWidthThin,'color',MyColor);
    end
    hold on;
    plot(paramLocalConst{1}(end,1), 0,'r*','MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
end

legendPos='NorthEast';
if abs(t(1))>abs(t(end))
    legendPos='NorthWest';
end
legend('Hist', 'Param GGD', 'Nonparam Gau', 'MaxEnt wo local', 'MaxEnt w Gau', 'Location', legendPos);
ylabel('p(x)','FontSize',MyFontSize,'FontName','Arial');
xlabel('x','FontSize',MyFontSize,'FontName','Arial');
set(gca,'FontSize',MyFontSize);
set(gca,'FontName','Arial');

title(['fMRI M=', num2str(M)]);

i=1;


function properties = getopt(properties,varargin)
%GETOPT - Process paired optional arguments as 'prop1',val1,'prop2',val2,...
%
%   getopt(properties,varargin) returns a modified properties structure,
%   given an initial properties structure, and a list of paired arguments.
%   Each argumnet pair should be of the form property_name,val where
%   property_name is the name of one of the field in properties, and val is
%   the value to be assigned to that structure field.
%
%   No validation of the values is performed.
%%
% EXAMPLE:
%   properties = struct('zoom',1.0,'aspect',1.0,'gamma',1.0,'file',[],'bg',[]);
%   properties = getopt(properties,'aspect',0.76,'file','mydata.dat')
% would return:
%   properties =
%         zoom: 1
%       aspect: 0.7600
%        gamma: 1
%         file: 'mydata.dat'
%           bg: []
%
% Typical usage in a function:
%   properties = getopt(properties,varargin{:})

% Function from
% http://mathforum.org/epigone/comp.soft-sys.matlab/sloasmirsmon/bp0ndp$crq5@cui1.lmms.lmco.com

% dgleich
% 2003-11-19
% Added ability to pass a cell array of properties

if ~isempty(varargin) && (iscell(varargin{1}))
    varargin = varargin{1};
end;

% Process the properties (optional input arguments)
prop_names = fieldnames(properties);
TargetField = [];
for ii=1:length(varargin)
    arg = varargin{ii};
    if isempty(TargetField)
        if ~ischar(arg)
            error('Property names must be character strings');
        end
        %f = find(strcmp(prop_names, arg));
        if isempty(find(strcmp(prop_names, arg),1)) %length(f) == 0
            error('%s ',['invalid property ''',arg,'''; must be one of:'],prop_names{:});
        end
        TargetField = arg;
    else
        properties.(TargetField) = arg;
        TargetField = '';
    end
end
if ~isempty(TargetField)
    error('Property names and values must be specified in pairs.');
end


function showDebugInfo(titleStr, Lambda, Alpha, stabFlag, fConstCost)
disp(titleStr);
disp(['Lambda:', num2str(reshape(Lambda, 1, prod(size(Lambda))))]);
disp(['Alpha:', num2str(reshape(Alpha, 1, prod(size(Alpha))))]);
disp(['Stable:', num2str(stabFlag)]);

if prod(size(fConstCost))~=0
    figure;
    MyFontSize = 12;
    MyMarkerSize = 10;
    MyLineWidth = 1.3;
    plot(fConstCost, '-b', 'MarkerSize', MyMarkerSize, 'linewidth', MyLineWidth);
    title(titleStr,'FontSize',MyFontSize,'FontName','Arial')
    ylabel('Cost','FontSize',MyFontSize,'FontName','Arial');
    xlabel('Number of iterations','FontSize',MyFontSize,'FontName','Arial');
    set(gca,'FontSize',MyFontSize);
    set(gca,'FontName','Arial');
    set(gca,'YScale','log');
end
