function [alfa, beta, gama, SROCC, LCC, MSE] = ParOptim_ALM (X, Y, K)
% (c) ZHANG Fan, Oct. 2010 ~ Apr. 2013
% Parameter optimization of Additive Log-logisitc Model
% ALM: q = 1/(1+a*sum(exp(gama)*Am.^alfa.*Ph.^beta)^b)
% Am                 - Amplitude feature vectors for each scale of an image
% Ph                 - Phase feature vectors for each scale of an image
% alfa, beta, & gama - parameters of models for each scale
% a & b              - parameters of each database for monotonic mapping
% q                  - normalized image quality
% Given features Am & Ph and MOS (ground truth of q) of multiple databases,
% this code solves the parameters leading to maximum likelihood.
% Formatted input
% X                  - features{DBs, scales}(patches,images, feature types)
% Y                  - q{DBs}(images)
% K                  - number of patches per DB per scale {DBs, scales}
% System requirment: 8G+ memory, 64-bit OS

% Reference:
% F. Zhang, W. Lin, Z. Chen and K. N. Ngan, ¡§Additive log-logistic model 
% for video quality assessment,¡¨IEEE Transaction on Image Processing. 
% volumn 22, issue 4, pp. 1536-1547, 2013.

global dbselect; % Index of selected databases
dbselect = [1:8 12:14];
global distype;  % Index of selected distortion types for each database
distype = {1:5,1:5,1:2,1:15,1:6,1:2,1:5,1:3,[],[],1:10,1:2,1:6,1:2};
ts = 3;          % The number of image scales to be considered
global X; global K; global Y;
X = cell(length(dbselect),ts); K = cell(length(dbselect),ts); Y = cell(length(dbselect),1); % 14 DB, 5 scales

isaconfig = {'804'};t=1; % 804 - a total of four 8-dimensinal ISA subspaces
%isaconfig = {'404','604','804','1004','1204','1404','1604'};
%isaconfig = {'232','321','414','512','610','709','808','907','106','125','164','213','322','641'};
databasename = {'LIVE','IVC','Toyama','TID2008','A57','WIQ','CSIQ','LAR',...
                '','','Enrico','BA','FourierSB','Meerwald'};

%% Format the input (features and DMOS)
if exist('input.mat','file')
    for scale = 1 : ts
        for dbind = 1:length(dbselect)
            load(['.\Features\' isaconfig{t} '\' databasename{dbselect(dbind)} num2str(scale) '.mat'],'T','R1','subject');%T: Ph. R1: Am. Subject: MOS
            ind = ismember(subject(:,2),distype{dbselect(dbind)});
            X{dbind,scale}(:,:,1) = R1(:,ind); % Am,Amplitude
            X{dbind,scale}(:,:,2) = T(:,ind);  % Ph,Phase
            K{dbind,scale} = subject(ind,4)';  % MOS/DMOS, ground truth of q
            % Normalized to DMOS
            if ismember(dbselect(dbind),[2 3 4 8 12 13 14])
                subject(ind,1) = 1 - subject(ind,1);
            end
            if scale == 1
                Y{dbind}(:,1) = subject(ind,1); % DMOS
            end
            clear T; clear R1; clear subject;
        end
    end
    % save -v7.3 input.mat X K Y;
else
    load input.mat;
end

%% parameter initialization
alfa = [2.1 .6 1.2 1.6 1.8]; 
beta = [1 1 1 1 1 1]; 
gama = log([1 1000 1000 5 10]); 
a = [6.761441,17.84420,18.13421,5.215124,16.04316,5.546129,9.111989,14.57309,0,0,0,11.66329,11.34763,14.18893];
b = [-0.776,-2.184,-2.542,-0.6059,-1.922,-0.692,-1.036,-1.838,1,1,1,-1.410,-1.691,-1.994];
x0 = [alfa(1:ts) beta(1:ts) gama(2:ts) a(dbselect) b(dbselect)];

%% parameter opmtimization
options = optimset('GradObj','on','Display','iter','Hessian','on');
[x,fval] = fminunc(@Proposal,x0,options);
alfa = x(1:ts);
beta = x(ts+1:2*ts);
gama = [0 x(2*ts+1:3*ts-1)];

%% performance calculation
[SROCC,LCC,MSE]=Vis(x);
end


function [f,g,H] = Proposal(x)
global X; global K; global Y; global dbselect;
f = 0;
g = zeros(length(x),1,'double'); 
H = zeros(length(x),'double');   
ts = size(X,2); % total number of scales
if sum(x(1:ts*2)<0) > 0 % alfa & beta should be non-negative
    f = 1e8;
else

for dbind = 1 : length(dbselect)
    N = size(X{dbind,1},2);
    d = zeros(ts,N,'single'); dlogAm = zeros(ts,N,'single'); dlogPh = zeros(ts,N,'single');
    dlog2Am = zeros(ts,N,'single'); dlog2Ph = zeros(ts,N,'single'); dlogAmlogPh = zeros(ts,N,'single');
    for scale = 1 : ts
        alfa = x(scale);
        beta = x(scale+ts);
        if scale <=3
            ap = abs(X{dbind,scale}(:,:,1)).^alfa.*X{dbind,scale}(:,:,2).^beta;
        else % keep the sign of Am
            ap = sign(X{dbind,scale}(:,:,1)).*abs(X{dbind,scale}(:,:,1)).^alfa.*X{dbind,scale}(:,:,2).^beta;
        end
        logap = log(abs(X{dbind,scale}(:,:,1:2))+eps);
        % sum(Am^alfa.*Ph^beta)
        d(scale,:) = sum(ap)./K{dbind,scale};
        % ~dD/dalfa: sum(Am^alfa.*Ph^beta.*logAm)
        dlogAm(scale,:) = sum(ap.*logap(:,:,1))./K{dbind,scale};
        % ~dD/dbeta: sum(Am^alfa.*Ph^beta.*logPh)
        dlogPh(scale,:) = sum(ap.*logap(:,:,2))./K{dbind,scale};
        % ~d2D/d2alfa: sum(Am^alfa.*Ph^beta.*logAm^2)
        dlog2Am(scale,:) = sum(ap.*logap(:,:,1).^2)./K{dbind,scale};
        % ~d2D/d2beta: sum(Am^alfa.*Ph^beta.*logPh^2)
        dlog2Ph(scale,:) = sum(ap.*logap(:,:,2).^2)./K{dbind,scale};
        % ~d2D/(dalfa dbeta): sum(Am^alfa.*Ph^beta.*logAm.*logPh)
        dlogAmlogPh(scale,:) = sum(ap.*logap(:,:,1).*logap(:,:,2))./K{dbind,scale};
        % Already normalized by K
    end
    gama = exp([0;x(2*ts+1:3*ts-1)']);
    D = sum(bsxfun(@times,d(1:ts,:),gama)); 
    if sum(D(:)<0)>0
        f = -1e8;
        break;
    end
    a = exp(x(dbind+3*ts-1));
    b = x(dbind+3*ts-1+length(dbselect));
    q = 1./(1 + a * D.^b);
    dq = (q - Y{dbind}(:,1)');             % Y{dbind}(:,1) - DMOS 
    dLdD = b*dq./D;                        % dlogL/dD
    d2Ld2D = (b*(q-1).*q-dq)./(D.*D);      % d2logL/d2D
    q1qD = (q-1).*q./D;                    % (q-1).*q./D
    bq1qDdqD = (b*(q-1).*q.*log(D)+dq)./D; % [b*(q-1).*q.*log(D)+dq]./D
    f = f + sum(Y{dbind}(:,1)'.*log(q./(1-q+eps))+log(1-q+eps));
    for scale = 1 : ts
        g(scale) = g(scale) + gama(scale)*sum(dLdD.*dlogAm(scale,:));      % dL/dalfa
        g(scale+ts) = g(scale+ts) + gama(scale)*sum(dLdD.*dlogPh(scale,:));% dL/dbeta
        if scale >1  % gama(1)=1 constantly
            g(scale+2*ts-1) = g(scale+2*ts-1) + gama(scale)*sum(dLdD.*d(scale,:)); % dL/dloggama
        end
        %  d2L/d2D * dD/dai * dD/daj  
        for j = 1 : ts
            H(scale,j)=H(scale,j)+gama(scale)*gama(j)*sum(d2Ld2D.*dlogAm(scale,:).*dlogAm(j,:));
            H(scale,j+ts)=H(scale,j+ts)+gama(scale)*gama(j)*sum(d2Ld2D.*dlogAm(scale,:).*dlogPh(j,:));
            H(scale+ts,j+ts)=H(scale+ts,j+ts)+gama(scale)*gama(j)*sum(d2Ld2D.*dlogPh(scale,:).*dlogPh(j,:));
            if j > 1 % gama(1)=1 constantly
                H(scale,j+2*ts-1)=H(scale,j+2*ts-1)+gama(scale)*gama(j)*sum(d2Ld2D.*dlogAm(scale,:).*d(j,:));
                H(scale+ts,j+2*ts-1)=H(scale+ts,j+2*ts-1)+gama(scale)*gama(j)*sum(d2Ld2D.*dlogPh(scale,:).*d(j,:));
                H(scale+2*ts-1,j+2*ts-1)=H(scale+2*ts-1,j+2*ts-1)+gama(scale)*gama(j)*sum(d2Ld2D.*d(scale,:).*d(j,:));
            end
        end
        % Note that mostly d2D/daidaj=0, except the same scale.
        H(scale,scale)=H(scale,scale)+gama(scale)*sum(dLdD.*dlog2Am(scale,:));
        H(scale,scale+ts)=H(scale,scale+ts)+gama(scale)*sum(dLdD.*dlogAmlogPh(scale,:));
        H(scale+ts,scale+ts)=H(scale+ts,scale+ts)+gama(scale)*sum(dLdD.*dlog2Ph(scale,:));
        H(scale,dbind+3*ts-1) = b*gama(scale)*sum(q1qD.*dlogAm(scale,:));
        H(scale+ts,dbind+3*ts-1) = b*gama(scale)*sum(q1qD.*dlogPh(scale,:));
        H(scale,dbind+3*ts-1+length(dbselect)) = gama(scale)*sum(bq1qDdqD.*dlogAm(scale,:));
        H(scale+ts,dbind+3*ts-1+length(dbselect)) = gama(scale)*sum(bq1qDdqD.*dlogPh(scale,:));
        if scale > 1
            H(scale,scale+2*ts-1)=H(scale,scale+2*ts-1)+gama(scale)*sum(dLdD.*dlogAm(scale,:));
            H(scale+ts,scale+2*ts-1)=H(scale+ts,scale+2*ts-1)+gama(scale)*sum(dLdD.*dlogPh(scale,:));
            H(scale+2*ts-1,scale+2*ts-1)=H(scale+2*ts-1,scale+2*ts-1)+gama(scale)*sum(dLdD.*d(scale,:));
            H(scale+2*ts-1,dbind+3*ts-1) = b*gama(scale)*sum(q1qD.*d(scale,:));
            H(scale+2*ts-1,dbind+3*ts-1+length(dbselect)) = gama(scale)*sum(bq1qDdqD.*d(scale,:));
        end
    end
    g(dbind+3*ts-1) = sum(dq);                         % dL/dloga
    g(dbind+3*ts-1+length(dbselect)) = sum(dq.*log(D));% dL/db
    H(dbind+3*ts-1,dbind+3*ts-1) = sum((q-1).*q);      % d2L/d2a
    H(dbind+3*ts-1,dbind+3*ts-1+length(dbselect)) = sum((q-1).*q.*log(D)); % d2L/(dadb)
    H(dbind+3*ts-1+length(dbselect),dbind+3*ts-1+length(dbselect)) = sum((q-1).*q.*log(D).^2); % d2L/d2b
    H = tril(H',-1)+triu(H);
end
f = double(-f); g = -g; H = -H; % Due to maximization
end
end


function [SROCC,LCC,MSE]=Vis(x)
global X; global K; global Y;global dbselect;
ts = size(X,2); 
LCC = zeros(1,length(dbselect)); 
SROCC = zeros(1,length(dbselect));
MSE = zeros(1,length(dbselect));
for dbind = 1:length(dbselect) 
    d = zeros(ts,size(X{dbind,1},2),'single');
    for scale = 1 : ts
        d(scale,:) = sum(abs(X{dbind,scale}(:,:,1)).^x(scale).*X{dbind,scale}(:,:,2).^x(scale+ts))./K{dbind,scale};
    end
    for scale = 4 : ts
        d(scale,:) = sum(sign(X{dbind,scale}(:,:,1)).*abs(X{dbind,scale}(:,:,1)).^x(scale).*X{dbind,scale}(:,:,2).^x(scale+ts))./K{dbind,scale};
    end
    D = sum(bsxfun(@times,d(1:ts,:),exp([0;x(2*ts+1:3*ts-1)']))); 
    q = 1./(1 + exp(x(dbind+3*ts-1)) * D.^x(dbind+3*ts-1+length(dbselect)));
    
    SROCC(dbind) = corr(Y{dbind}(:,1),q','type','spearman');
    % Adaptive monotonic mapping for each DB
    b = glmfit(log(D'),Y{dbind}(:,1),'binomial','link','logit');
    q = glmval(b,log(D'),'logit'); q = q';
    LCC(dbind) = corr(Y{dbind}(:,1),q','type','pearson');
    MSE(dbind) = sum((Y{dbind}(:,1)-q').^2)/length(q);
end
end