function [beta alpha b outer_itc] = QMKL_release(K, Q, truth)
%function [beta gap alpha b outer_itc] = QMKL(K, Q, truth)

N = numel(truth);
M = size(K,1);

iQ = pinv(Q);


%A few constants
epsilon = 10^-7;
C = 100;

%Initial guess for beta
beta = ones(M,1) * (1/M)^(1/2);


%Flags for determining which solution method to use.
CVX_Beta = false;
NEWTON = true && ~CVX_Beta;
GSL = false && ~(NEWTON || CVX_Beta);
SDP =true && ~(GSL || NEWTON || CVX_Beta);

NORMALIZE = true;
PLOT = false;
VERB = false;
extra_VERB = false;

w_gap = 1;

outer_itc = 1;
while w_gap >= epsilon && outer_itc < 50
    %Solve SVM for alpha(Q)
    KM = reshape(beta' * reshape(K, M, []), N, N);    
    alpha = SK_SVM(KM, truth, C);
    
    %Sm(m) = alpha' squeeze(K(m,:,:)) * alpha
    Sm = (alpha' * reshape(alpha' * reshape(K, M * N, [])', M, N)')'/2;

    
    %Objective function with previous beta and new alpha
    %St = beta' * Sm - sum(abs(alpha));
    St = beta' * iQ * beta - sum(abs(alpha));

    %Re-compute beta.
    if PLOT
        figure(2); plot(Sm, 'r+'); title('S');
        figure(3); plot(beta' * Q, 'b+'); title('beta * Q');
    end

    old_beta = beta;
    
    %CVX solution pretty much agrees with CPLEX
    %Multiplying by 10^ some power is for numerical reasons.
    if CVX_Beta
        cvx_quiet(true);
        cvx_begin
        variable beta(M);

        minimize -10^0 * transpose(Sm) * beta
        subject to
        beta >= 0;
        10^0 * transpose(beta) * Q * beta <= 10^0;
        cvx_end
    end
    
    if NEWTON
        converged = false;
        Wm = Sm .* old_beta.^2;
        Wm(old_beta <= 0) = 0;
                
        itc = 1;
        while ~converged && itc <35
            d = 2 * diag(Wm./(beta.^3));
            d(isnan(d)) = 0;
            
            H = Q + d;
            g = Q * beta - Wm./(beta.^2);
            g(isnan(g)) = 0;
            
            beta = beta - pinv(H)*g;
            
            if numel(find(isnan(beta))) > 0
                numel(find(isnan(H)))
                numel(find(isnan(g)))
            end
            
            if extra_VERB
%                 [beta old_beta Sm Wm]
%                 [Q * beta       Wm./(beta.^2)     inv(H)*g   beta]
                disp(sprintf('negative beta: %d\tt: %.4f\t||g||: %e\t||inv(H)g||: %e', numel(find(beta < 0)), t, norm(g), norm(inv(H)*g)));
                pause;
            end
            
            converged = norm(g) < 10^-12;
            itc = itc + 1;
        end
        itc
        
        beta(beta < 0) = 0;
        beta(isnan(beta)) = 0;
    end
    
    
    
    if GSL
        Wm = Sm .* old_beta.^2;
        Wm(old_beta <= 0) = 0;
        
        beta = QMKL_GSL(Wm, Q, iQ);
        
        if extra_VERB
            figure(3), plot(beta, 'g+');
            title('beta');

            %H = Q + diag(Wm./(beta.^3));
            %g = Q * beta - Wm./(beta.^2);
            %figure(4), plot(g'*inv(H), 'g+');
            %title('Newton step direction');
            %lambda = sum(Wm./beta);
            lambda = sqrt((Wm./beta.^2)' * iQ * (Wm./beta.^2));
            figure(4), plot(Wm./(beta.^2), lambda*Q*beta, 'b+');

            lambda
            figure(5), plot(Q * beta, 'b+');

            pause;
        end
        
        beta(beta < 0) = 0;
    end
    
    
    if SDP
        Wm = Sm .* old_beta.^2;
        Wm(old_beta <= 0) = 0;
        
        cvx_quiet(true);
        cvx_begin
        variable nu(M,1);
        variable beta(M,1);
        variable Z(M,M);

        minimize sum(Wm .* nu)
        subject to
        for i=1:M
            [nu(i) 1; 1 beta(i)] == semidefinite(2);
        end
        
%         trace(Q * Z) <= 1;
%         [1 transpose(beta); beta Z] == semidefinite(M + 1);
        transpose(beta) * Q * beta <= 1;
        cvx_end
        
        
        beta = double(beta);
        beta(beta < 0) = 0;
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Normalize if necessary
    if NORMALIZE %&& abs(sqrt(beta' * Q * beta) - 1) > 10^-4
        if VERB
            format long g
            disp(sprintf('Normalized: %f', sqrt(beta' * Q * beta)));
            disp(sprintf('Mean: %f Max: %f Min: %f', mean(beta), max(beta), min(beta)));
            format short g
        end
        
%         if sqrt(beta' * Q * beta) > 10^-4
%             beta = beta / sqrt(beta' * Q * beta);
%         end
        
        beta = beta / norm(beta);
    end
    
   
    if PLOT
        figure(4); plot(old_beta, 'g+'); title('old beta');
        figure(5); plot(beta, 'g+'); title('new beta');
        %figure(7); plot(alpha, 'g+'); title('alpha');
        pause;
    end


    %This is what they do in the Shogun source code, so I'm following it
    %here.
    %St_new = beta' * Sm - sum(abs(alpha));
    St_new = beta' * iQ * beta - sum(abs(alpha));
    w_gap = abs(1 - (St_new/St));

    if VERB
        disp(sprintf('St: %f\t\tSt_new: %f\t\tAlpha: %f\t\tMax(beta): %f', St, St_new, sum(abs(alpha)), max(beta)));
    end
%    gap = [gap w_gap];
    outer_itc = outer_itc + 1;
end

outer_itc

KM = reshape(beta' * reshape(K, M, []), N, N);
[alpha b] = SK_SVM(KM, truth, C);
end