function [rho_R, rho_G, rho_B]=RGB_components_optim(T)
% Solves min sum(c={r,g,b})(sum(i=1:35)(zc_i+1 - zc_i)^2)
% s.t. T exp(zc) = s, c = {r,g,b}, s = {[1;0;0], [0;1;0], [0;0;1]}
% rho_R + rho_G + rho_B < 1, 
% where zc=ln(rho_c), 
% using Matlab's constrained optimization function fmincon.

% T is 3x36 matrix converting reflectance to D65-referenced linear rgb: T*rho=rgb.
% Each rho is a 36x1 vector of reconstructed reflectance values, in the range (0-1].

% For more information, see:
% http://scottburns.us/fast-rgb-to-spectrum-conversion-for-reflectances/
% Written by Scott Allen Burns, September 2018.
% Licensed under a Creative Commons Attribution-ShareAlike 4.0 International
% License (http://creativecommons.org/licenses/by-sa/4.0/).

options=optimoptions(@fmincon,'Algorithm','sqp', 'MaxFunEvals',500000, 'MaxIter',3000,...
    'ConstraintTolerance',1.0e-8, 'OptimalityTolerance',1.0e-8);
z_soln=fmincon(@my_obj, zeros(108,1), [],[],[],[],[],[], @my_con, options);
rho_R=exp(z_soln(1:36));
rho_G=exp(z_soln(37:72));
rho_B=exp(z_soln(73:108));

% nested functions for objective function and constraints

    function f=my_obj(z)
        % minimize norm of differences of successive log reflectances
        diff=z(1:35)-z(2:36);
        f=norm(diff);
        diff=z(37:71)-z(38:72);
        f=f+norm(diff);
        diff=z(73:107)-z(74:108);
        f=f+norm(diff);
    end

    function [c,ceq]=my_con(z)
        % require rgb of components to match
        ceq=zeros(9,1);
        current_rhoR=exp(z(1:36));
        current_rgbR=T*current_rhoR;
        ceq(1:3)=current_rgbR-[1;0;0];
        current_rhoG=exp(z(37:72));
        current_rgbG=T*current_rhoG;
        ceq(4:6)=current_rgbG-[0;1;0];
        current_rhoB=exp(z(73:108));
        current_rgbB=T*current_rhoB;
        ceq(7:9)=current_rgbB-[0;0;1];
        % upper limit on sum of components
        c = current_rhoR + current_rhoG + current_rhoB - 1;
    end
end