text:
function [dGf0, dHf0, mf, aveHbound, aveZi, lambda, gpfnsp] = calcdGHT(dGzero, dHzero, zi, nH, pHr, is, temp, chi, Legendre, LegendreCHI, printLevel)
% Calculates the standard transformed Gibbs energy of a reactant
%
% Reproduces the function of T (in Kelvin), pHa (electrode pH), and ionic strength (is) that
% gives the standard transformed Gibbs energy of formation of a reactant
% (sum of species) and the standard transformed enthalpy of a reactant.
%
% Assuming p pseudoisomer species corresponding to one reactant
%
% Optional output dependent on multiple precision toolbox
%
% USAGE:
%
%    [dGf0, dHf0, mf, aveHbound, aveZi, lambda, gpfnsp] = calcdGHT(dGzero, dHzero, zi, nH, pHr, is, temp, chi, Legendre, LegendreCHI, printLevel)
%
% INPUTS:
%    dGzero:         `p x 1` standard Gibbs energy of formation at 298.15 K
%    zi:             `p x 1` electric charge
%    nH:             `p x 1` number of hydrogen atoms in each species
%    pHr:            real pH of 5 to 9 (see realpH.m)
%    is:             ionic strength 0 to 0.35 M
%    temp:           temperature 273.15 K to 313.15 K
%
% OPTIONAL INPUT:
%    dHzero:         `p x 1` standard enthalpy of formation at 298.15 K
%    chi:            electrical potential
%    Legendre:       {(1), 0} Legendre Transformation for specifc pHr?
%    LegendreCHI:    {(1), 0} Legendre Transformation for specifc electrical potential?
%
% OUTPUT:
%    dGf0:           reactant standard transformed Gibbs energy of formation
%    dHf0:           reactant standard transformed enthalpy of formation
%    mf:             mole fraction of each species within a pseudoisomer group
%    aveHbound:      average number of protons bound to a reactant
%    aveZi:          average charge of a reactant
%    lambda:         activity coefficient for each metabolite species
%    gpfnsp:         metabolite species standard transformed Gibbs energy of formation
%
% OPTIONAL OUTPUT:
%
%    dHf0:           standard transformed enthalpy of formation
%
% The values of the standard transformed Gibbs energy of formation
% and the standard transformed enthalpy of formation can be calculated
% temperature in the range 273.15 K to 313.15 K, using the van't Hoff equation.
% `pHr` in the range 5 to 9 (these correspond to the pH range of the species in Alberty's tables)
% ionic strength in the range 0 to 0.35 M.
% See Mathematica program `calcdGHT p289 Alberty 2003`
%
% The changes in the values of standard transformed Gibbs energy of formation
% and the standard transformed enthalpy of formation might be improved if
% knowlegde of standard molar heat capacity was available for each species
% See `p41 Alberty 2003`
%
% Multiple Precision Toolbox for MATLAB by Ben Barrowes (mptoolbox_1.1)
% http://www.mathworks.com/matlabcentral/fileexchange/6446
%
% .. Author: -Ronan M.T. Fleming

if ~exist('printLevel','var')
    printLevel=0;
end

electricalTerm = 0; % initialize electricalTerm

if pHr<5 || pHr>9
    if 0
        error('pHr out of applicable range, 5 - 9.');
    else
        warning('pHr out of applicable range, 5 - 9.');
    end
end
if temp<273.15 || temp>313.15
    error('temperature out of applicable range, 273.15 - 313.15');
end
if is<0 || is>0.35
    error('ionic strength out of applicable range, 0 - 0.35');
end

if ~exist('Legendre','var')
    Legendre=1;
end
if ~exist('LegendreCHI','var')
    LegendreCHI=1;
end

%check if multiple precision toolbox is properly installed
if strcmp(which('mp'),'') || 1 %TODO install
    if printLevel>0
        fprintf('%s\n','No multiple precision toolbox: NaN returned if exp(x) gets too large');
    end
    R=8.31451;
    %Energies are expressed in kJ mol^-1.*)
    R=R/1000;
    %standard temperature with a capital T
    T=298.15;

    %Faraday Constant (kJ/mol)
    F=96.48; %kJ/mol

    p=9.20483;
    q=10^3;
    r=1.284668;
    s=10^5;
    u=4.95199;
    v=10^8;
    % Eq. 3.7-3 p48 Alberty 2003
    % Around 298.15K
    % alpha = 1.10708 - 1.54508*temp/(10^3) +5.95584*temp^2/(10^6)
    % Eq. 3.7-4 p49 Alberty 2003
    % R*T*alpha, where alpha is the temperature dependent coeffficient, sometimes A,
    % in the extended Debye-Huckel equation
    gibbscoeff = (9.20483*temp)/10^3 - (1.284668*temp^2)/10^5 + (4.95199*temp^3)/10^8;

    %If standard enthalpy of formation is known, and independent of
    %temperature an adjustment for temperature can be made.
    %(calcdGHT p289 Alberty 2003)
    if isempty(dHzero)
        %dGzeroT = (dGzero*temp)/T + dHzero*(1 - temp/T);
        dGzeroT = dGzero;%(dGzero*temp)/T + dHzero*(1 - temp/T);
    else
        % van't Hoff equation
        dGzeroT = (dGzero*temp)/T + dHzero*(1 - temp/T);
    end

    %pHr
    if Legendre
        %Eq. 4.4-9/10 p67 Alberty 2003
        %note the use of culture temperature
        pHterm = nH*R*temp*log(10^-pHr);
        %Eq 4.4-10 Alberty 2003 with temp dependent gibbscoeff
        istermG = (gibbscoeff*(zi.^2 - nH)*is^0.5)/(1 + 1.6*is^0.5);
    else
        %no Legendre transformation for pHr
        pHterm = 0;
        %Eq 3.6-3 Alberty 2003 with temp dependent gibbscoeff
        istermG = (gibbscoeff*(zi.^2)*is^0.5)/(1 + 1.6*is^0.5); %omit the -nH if  no Legendre
    end

    if LegendreCHI
        if 0
            %By convention, we assume the chemical potential of a metabolite
            %includes an electrical potential term
            % u = u0 + RT*log(activity) + F*zi*chi;
            %The Legendre transformation for electrical potential is
            % u' = u -  F*zi*chi = u0 + RT*log(activity);
            %So the following line will negate the effect of an electrical
            %potential ONLY if it has previously been added.
            electricalTerm=-(F*(chi/1000))*zi;
            %eq 8.5-1 p148 Alberty 2003
        else
            %Imaginary Legendre Transformation for Electrical Potential, to
            %take account of the fact that we have not previously added an
            %electrical potential term to the standard Gibbs energy.
            electricalTerm=0;
            %The charge and change in chemical potential for multiphase
            %reactions is taken into account in
            %deltaG0concFluxConstraintBounds.m
        end
    end

    %standard transformed Gibbs energy of each species
    gpfnsp = dGzeroT - pHterm - istermG - electricalTerm;

    %isomer group thermodynamic standard transformed Gibbs energy of
    %formation for a reactant with more than one metabolite species
    if length(dGzero)==1
        dGf0=gpfnsp;
    else
        %need to approximate log(sum(exp(-gpfnsp/(R*temp))));
        dGf0 = -R*temp*maxstar(-gpfnsp/(R*temp));
    end

    %Mole fraction
    %see 3.5-12 p45 Alberty 2003
    mf=exp((dGf0-gpfnsp)/(R*temp));
%     fprintf('%s\n','Cannot calculate mole fractions without multiple
%     precision toolbox');

    %activity coefficient
    lambda=double(exp(-(gibbscoeff*(zi.^2)*is^0.5)/(1 + 1.6*is^0.5)/(R*temp)));

    %average number of H+ ions bound by a reactant
    aveHbound=mf'*nH;

    %average charge of a reactant
    aveZi=mf'*zi;
%     fprintf('%s\n',int2str(length(dGzero)));

    %isomer group thermodynamic standard transformed Enthalpy of
    %formation for a reactant
    %%%%%%% makes script faster to leave this out for now.
    dHzero=[];
    if ~isempty(dHzero)
        %make temperature a smaller variable
        t=temp;
        switch length(dGzero)
            case 1
                %corresponds to Simplify[-(t^2*D[dGf0/t, t])] in Albertys code for
                %one species reactant
                %see Mathematica file dHfn.nb
                A=dGzero(1);
                B=dHzero(1);
                C=zi(1);
                D=nH(1);
                dHf0 =(B*(1+1.6*is^0.5)*s*v + (C^2)*(is^0.5)*(t^2)*(2*s*t*u-r*v)+D*(is^0.5)*(t^2)*(-2*s*t*u+r*v)) / ((1+1.6*is^0.5)*s*v);

                if isnan(dHf0)
                    error('No multiple precision toolbox: NaN returned if exp(x) gets too large')
                end
            case 2
                %see Mathematica file dHfn.nb
                A=dGzero(1);
                B=dHzero(1);
                C=zi(1);
                D=nH(1);
                a=dGzero(2);
                b=dHzero(2);
                c=zi(2);
                d=nH(2);
                %translated to matlab from mathematica by hand
                dHf0 =((10^-pHr)^d*exp(-(b*((1/t)-(1/T))+a/T-(((c^2-d)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R)*(b*(1+1.6*is^0.5)*s*v+...
                    c^2*is^0.5*t^2*(2*s*t*u-r*v)+d*is^0.5*t^2*(-2*s*t*u+r*v))+...
                    (10^-pHr)^D*exp(-(B*((1/t)-(1/T))+A/T-(((C^2-D)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R)*(B*(1+1.6*is^0.5)*s*v+...
                    C^2*is^0.5*t^2*(2*s*t*u-r*v)+D*is^0.5*t^2*(-2*s*t*u+r*v)))/...
                    (((10^-pHr)^d*exp((b*((-1/t)+(1/T))-a/T+(((c^2-d)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R)+...
                    (10^-pHr)^D*exp((B*((-1/t)+(1/T))-A/T+(((C^2-D)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R))*(1+1.6*is^0.5)*s*v);
                %         Mathematica Expression to Matlab m-file Converter by Harri Ojanen, Rutgers University
                %         dHfn2=((10.^((-1).*pHr)).^d.*exp(R.^(-1).*(b.*((-1).*t.^(-1)+T.^(-1))+(-1).*a.* ...
                %             T.^(-1)+0.1E1.*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).* ...
                %             is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).* ...
                %             r.*v))))+(10.^((-1).*pHr)).^D.*exp(R.^(-1).*(B.*((-1).*t.^(-1)+T.^(-1))+( ...
                %             -1).*A.*T.^(-1)+0.1E1.*(C.^2+(-0.1E1).*D).*(0.1E1+0.16E1.*is.^0.5E0).^( ...
                %             -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                %             -0.1E1).*r.*v))))).^(-1).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).*s.^(-1).*v.^( ...
                %             -1).*((10.^((-1).*pHr)).^d.*exp((-1).*R.^(-1).*(b.*(t.^(-1)+(-1).*T.^(-1) ...
                %             )+a.*T.^(-1)+(-0.1E1).*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.*is.^0.5E0).^( ...
                %             -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                %             -0.1E1).*r.*v)))).*(b.*(1+0.16E1.*is.^0.5E0).*s.*v+c.^2.*is.^0.5E0.* ...
                %             t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+d.*is.^0.5E0.*t.^2.*((-0.2E1).*s.* ...
                %             t.*u+r.*v))+(10.^((-1).*pHr)).^D.*exp((-1).*R.^(-1).*(B.*(t.^(-1)+(-1).* ...
                %             T.^(-1))+A.*T.^(-1)+(-0.1E1).*(C.^2+(-0.1E1).*D).*(0.1E1+0.16E1.* ...
                %             is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                %             s.*t.*u+(-0.1E1).*r.*v)))).*(B.*(1+0.16E1.*is.^0.5E0).*s.*v+C.^2.* ...
                %             is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+D.*is.^0.5E0.*t.^2.*(( ...
                %             -0.2E1).*s.*t.*u+r.*v)));

                if isnan(dHf0)
                    %see Mathematica file dHfn.nb
                    error('No multiple precision toolbox: NaN returned if exp(x) gets too large')
                end
            case 3

                A=dGzero(1);
                B=dHzero(1);
                C=zi(1);
                D=nH(1);
                a=dGzero(2);
                b=dHzero(2);
                c=zi(2);
                d=nH(2);
                e=dGzero(3);
                f=dHzero(3);
                g=zi(3);
                h=nH(3);
                %see Mathematica file dHfn.nb
                %Mathematica Expression to Matlab m-file Converter by Harri Ojanen, Rutgers University
                dHf0 = ((10.^((-1).*pHr)).^d.*exp(R.^(-1).*(b.*((-1).*t.^(-1)+T.^(-1))+(-1).*a.* ...
                    T.^(-1)+0.1E1.*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).* ...
                    is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).* ...
                    r.*v))))+(10.^((-1).*pHr)).^D.*exp(R.^(-1).*(B.*((-1).*t.^(-1)+T.^(-1))+( ...
                    -1).*A.*T.^(-1)+0.1E1.*(C.^2+(-0.1E1).*D).*(0.1E1+0.16E1.*is.^0.5E0).^( ...
                    -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                    -0.1E1).*r.*v))))+(10.^((-1).*pHr)).^h.*exp(R.^(-1).*(f.*((-1).*t.^(-1)+ ...
                    T.^(-1))+(-1).*e.*T.^(-1)+0.1E1.*(g.^2+(-0.1E1).*h).*(0.1E1+0.16E1.* ...
                    is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                    s.*t.*u+(-0.1E1).*r.*v))))).^(-1).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).*s.^( ...
                    -1).*v.^(-1).*((10.^((-1).*pHr)).^d.*exp((-1).*R.^(-1).*(b.*(t.^(-1)+(-1) ...
                    .*T.^(-1))+a.*T.^(-1)+(-0.1E1).*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.* ...
                    is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                    s.*t.*u+(-0.1E1).*r.*v)))).*(b.*(1+0.16E1.*is.^0.5E0).*s.*v+c.^2.* ...
                    is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+d.*is.^0.5E0.*t.^2.*(( ...
                    -0.2E1).*s.*t.*u+r.*v))+(10.^((-1).*pHr)).^D.*exp((-1).*R.^(-1).*(B.*( ...
                    t.^(-1)+(-1).*T.^(-1))+A.*T.^(-1)+(-0.1E1).*(C.^2+(-0.1E1).*D).*(0.1E1+ ...
                    0.16E1.*is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.* ...
                    v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v)))).*(B.*(1+0.16E1.*is.^0.5E0).*s.*v+ ...
                    C.^2.*is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+D.*is.^0.5E0.* ...
                    t.^2.*((-0.2E1).*s.*t.*u+r.*v))+(10.^((-1).*pHr)).^h.*exp((-1).*R.^(-1).* ...
                    (f.*(t.^(-1)+(-1).*T.^(-1))+e.*T.^(-1)+(-0.1E1).*(g.^2+(-0.1E1).*h).*( ...
                    0.1E1+0.16E1.*is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*( ...
                    p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v)))).*(f.*(1+0.16E1.*is.^0.5E0).* ...
                    s.*v+g.^2.*is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+h.* ...
                    is.^0.5E0.*t.^2.*((-0.2E1).*s.*t.*u+r.*v)));

                %       dHf0 = ((10.^((-1).*pHr)).^d.*...
                %             exp(R.^(-1).*(b.*((-1).*t.^(-1)+T.^(-1))+(-1).*a.* ...
                %             T.^(-1)+0.1E1.*(c.^2+(-0.1E1).*d).*(1+1.6*is.^0.5).^(-1).* ...
                %             is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).* ...
                %             r.*v))))...
                %             +(10.^((-1).*pHr)).^D.*...
                %             exp(R.^(-1).*(B.*((-1).*t.^(-1)+T.^(-1))+( ...
                %             -1).*A.*T.^(-1)+0.1E1.*(C.^2+(-0.1E1).*D).*(1+1.6*is.^0.5).^( ...
                %             -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                %             -0.1E1).*r.*v))))...
                %             +(10.^((-1).*pHr)).^h.*...
                %             exp(R.^(-1).*(f.*((-1).*t.^(-1)+ ...
                %             T.^(-1))+(-1).*e.*T.^(-1)+0.1E1.*(g.^2+(-0.1E1).*h).*(0.1E1+0.16E1.* ...
                %             is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                %             s.*t.*u+(-0.1E1).*r.*v))))...
                %             ).^(-1).*(1+1.6*is.^0.5).^(-1).*s.^(-1).*v.^(-1).*...%%% end of denominator
                %             ((10.^((-1).*pHr)).^d.*...
                %             exp((-1).*R.^(-1).*(b.*(t.^(-1)+(-1) ...
                %             .*T.^(-1))+a.*T.^(-1)+(-0.1E1).*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.* ...
                %             is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                %             s.*t.*u+(-0.1E1).*r.*v))))...
                %             .*(b*(1+1.6*is^0.5)*s*v...
                %             +c^2*is.^0.5*t^2*(2.*s.*t.*u-r*v)+d.*is.^0.5*t^2*(-2*s*t*u+r*v))...
                %             +(10.^((-1).*pHr)).^D.*...
                %             exp((-1).*R.^(-1).*(B.*( ...
                %             t.^(-1)+(-1).*T.^(-1))+A.*T.^(-1)+(-0.1E1).*(C.^2+(-0.1E1).*D).*(0.1E1+ ...
                %             0.16E1.*is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.* ...
                %             v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v))))...
                %             .*(B.*(1+1.6*is.^0.5).*s.*v+ ...
                %             C.^2*is^0.5*t^2*(2*s*t*u-r*v)+D*is.^0.5*t^2*(-2*s*t*u+r*v))...
                %             +(10.^(-pHr)).^h.*...
                %             exp((-1).*R.^(-1).* ...
                %             (f.*(t.^(-1)+(-1).*T.^(-1))+e.*T.^(-1)+(-0.1E1).*(g.^2+(-0.1E1).*h).*( ...
                %             1+1.6*is.^0.5).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*( ...
                %             p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v))))...
                %             *(f*(1+1.6*is.^0.5)*s*v+g^2*is^0.5*t.^2*(2*s*t*u-r*v)+h*is^0.5*t^2*(-2*s*t*u+r*v)));
                %
                %             denominatorInv = ((10.^((-1).*pHr)).^d.*...
                %             exp(R.^(-1).*(b.*((-1).*t.^(-1)+T.^(-1))+(-1).*a.* ...
                %             T.^(-1)+0.1E1.*(c.^2+(-0.1E1).*d).*(1+1.6*is.^0.5).^(-1).* ...
                %             is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).* ...
                %             r.*v))))...
                %             +(10.^((-1).*pHr)).^D.*...
                %             exp(R.^(-1).*(B.*((-1).*t.^(-1)+T.^(-1))+( ...
                %             -1).*A.*T.^(-1)+0.1E1.*(C.^2+(-0.1E1).*D).*(1+1.6*is.^0.5).^( ...
                %             -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                %             -0.1E1).*r.*v))))...
                %             +(10.^((-1).*pHr)).^h.*...
                %             exp(R.^(-1).*(f.*((-1).*t.^(-1)+ ...
                %             T.^(-1))+(-1).*e.*T.^(-1)+0.1E1.*(g.^2+(-0.1E1).*h).*(0.1E1+0.16E1.* ...
                %             is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                %             s.*t.*u+(-0.1E1).*r.*v))))...
                %             ).^(-1).*(1+1.6*is.^0.5).^(-1).*s.^(-1).*v.^(-1);%%% end of denominator
                %
                %             dInv1=(10.^((-1).*pHr)).^d.*...
                %             exp(R.^(-1).*(b.*((-1).*t.^(-1)+T.^(-1))+(-1).*a.* ...
                %             T.^(-1)+0.1E1.*(c.^2+(-0.1E1).*d).*(1+1.6*is.^0.5).^(-1).* ...
                %             is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).* ...
                %             r.*v))));
                %
                %         dInv2=(10.^((-1).*pHr)).^D.*...
                %             exp(R.^(-1).*(B.*((-1).*t.^(-1)+T.^(-1))+( ...
                %             -1).*A.*T.^(-1)+0.1E1.*(C.^2+(-0.1E1).*D).*(1+1.6*is.^0.5).^( ...
                %             -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                %             -0.1E1).*r.*v))));
                %
                %         dInv3=(10.^((-1).*pHr)).^h.*...
                %             exp(R.^(-1).*(f.*((-1).*t.^(-1)+ ...
                %             T.^(-1))+(-1).*e.*T.^(-1)+0.1E1.*(g.^2+(-0.1E1).*h).*(0.1E1+0.16E1.* ...
                %             is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                %             s.*t.*u+(-0.1E1).*r.*v))));
                %
                %         dinvEnd=(1+1.6*is.^0.5).^(-1).*s.^(-1).*v.^(-1);
                %
                %         denominatorInv=(dInv1+dInv2+dInv3)^-1*dinvEnd;
                %
                %
                %
                %
                %             numerator=((10.^((-1).*pHr)).^d.*...
                %             exp((-1).*R.^(-1).*(b.*(t.^(-1)+(-1) ...
                %             .*T.^(-1))+a.*T.^(-1)+(-0.1E1).*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.* ...
                %             is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                %             s.*t.*u+(-0.1E1).*r.*v))))...
                %             .*(b*(1+1.6*is^0.5)*s*v...
                %             +c^2*is.^0.5*t^2*(2.*s.*t.*u-r*v)+d.*is.^0.5*t^2*(-2*s*t*u+r*v))...
                %             +(10.^((-1).*pHr)).^D.*...
                %             exp((-1).*R.^(-1).*(B.*( ...
                %             t.^(-1)+(-1).*T.^(-1))+A.*T.^(-1)+(-0.1E1).*(C.^2+(-0.1E1).*D).*(0.1E1+ ...
                %             0.16E1.*is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.* ...
                %             v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v))))...
                %             .*(B.*(1+1.6*is.^0.5).*s.*v+ ...
                %             C.^2*is^0.5*t^2*(2*s*t*u-r*v)+D*is.^0.5*t^2*(-2*s*t*u+r*v))...
                %             +(10.^(-pHr)).^h.*...
                %             exp((-1).*R.^(-1).* ...
                %             (f.*(t.^(-1)+(-1).*T.^(-1))+e.*T.^(-1)+(-0.1E1).*(g.^2+(-0.1E1).*h).*( ...
                %             1+1.6*is.^0.5).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*( ...
                %             p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v))))...
                %             *(f*(1+1.6*is.^0.5)*s*v+g^2*is^0.5*t.^2*(2*s*t*u-r*v)+h*is^0.5*t^2*(-2*s*t*u+r*v)));
                %             dHf0=denominatorInv*numerator;

                if isnan(dHf0)
                    %see Mathematica file dHfn.nb
                    error('No multiple precision toolbox: NaN returned if exp(x) gets too large')
                end
            otherwise
                dHf0 = NaN;
        end
        if isnan(dHf0)
            %see Mathematica file dHfn.nb
            error('Too many species in reactant so standard transformed enthalpy is NaN')
        end
    else
        %changed to NaN 20 July 2009
        dHf0=NaN;
    end
else
    fprintf('%s\n','Using multiple precision toolbox');
    dGzero=mp(dGzero);
    dHzero=mp(dHzero);
    zi=mp(zi);
    nH=mp(nH);
    pHr=mp(pHr);
    is=mp(is);
    chi=mp(chi);
    %culture temp
    temp=mp(temp);

    R=mp(8.31451);
    %Energies are expressed in kJ mol^-1.*)
    R=R/1000;
    %standard temperature with a capital T
    T=mp(298.15);
    %Faraday Constant (kJ/mol)
    F=96.48; %kJ/mol
    F=mp(F);

    p=mp(9.20483);
    q=mp(10^3);
    r=mp(1.284668);
    s=mp(10^5);
    u=mp(4.95199);
    v=mp(10^8);

    %RTalpha p 49 Alberty 2003
    %where alpha is the Debye-Huckel Constant
    gibbscoeff = (9.20483*temp)/10^3 - (1.284668*temp^2)/10^5 + (4.95199*temp^3)/10^8;

    %If standard enthalpy of formation is known, and independent of
    %temperature an adjustment for temperature can be made.
    %(calcdGHT p289 Alberty 2003)
    if isempty(dHzero) || double(temp)==T;
        %dGzeroT = (dGzero*temp)/T + dHzero*(1 - temp/T);
        dGzeroT = dGzero;%(dGzero*temp)/T + dHzero*(1 - temp/T);
    else
        dGzeroT = (dGzero*temp)/T + dHzero*(1 - temp/T);
    end

    %pHr
    if Legendre
        %Eq. 4.4-9/10 p67 Alberty 2003
        %note the use of culture temperature
        pHterm = nH*R*temp*log(10^-pHr);
        %Eq 4.4-10 Alberty 2003 with temp dependent gibbscoeff
        istermG = (gibbscoeff*(zi.^2 - nH)*is^0.5)/(1 + 1.6*is^0.5);
    else
        %no Legendre transformation for pHr
        pHterm = 0;
        %Eq 3.6-3 Alberty 2003 with temp dependent gibbscoeff
        istermG = (gibbscoeff*(zi.^2)*is^0.5)/(1 + 1.6*is^0.5); %omit the -nH if  no Legendre
    end

    if LegendreCHI
        if 0
            %By convention, we assume the chemical potential of a metabolite
            %includes an electrical potential term
            % u = u0 + RT*log(activity) + F*zi*chi;
            %The Legendre transformation for electrical potential is
            % u' = u -  F*zi*chi = u0 + RT*log(activity);
            %So the following line will negate the effect of an electrical
            %potential ONLY if it has previously been added.
            electricalTerm=-(F*(chi/1000))*zi;
            %eq 8.5-1 p148 Alberty 2003
        else
            %Imaginary Legendre Transformation for Electrical Potential, to
            %take account of the fact that we have not previously added an
            %electrical potential term to the standard Gibbs energy.
            electricalTerm=0;
            %The charge and change in chemical potential for multiphase
            %reactions is taken into account in
            %deltaG0concFluxConstraintBounds.m
        end
    end

    %standard transformed Gibbs energy of each species
    gpfnsp = dGzeroT - pHterm - istermG - electricalTerm;

    %partition function
    pf=sum(exp(-gpfnsp/(R*temp)));

    %mole fraction
    if length(dGzero)==1
        mf=1;
    else
        %mole fraction of each species if there is more than one
        lin_gpfnsp=exp(-gpfnsp/(R*temp));
        %cast back into a double
        mf=double(lin_gpfnsp/pf);
    end

    %activity coefficient
    lambda=double(exp(-(gibbscoeff*(zi.^2)*is^0.5)/(1 + 1.6*is^0.5)/(R*temp)));

    %average number of H+ ions bound by a reactant
    aveHbound=mf'*double(nH);

    %average number of H+ ions bound by a reactant
    aveZi=mf'*double(zi);

    %isomer group thermodynamic standard transformed Gibbs energy of
    %formation for a reactant with more than one metabolite species
    if length(dGzero)==1
        dGf0=gpfnsp;
    else
        dGf0 = -R*temp*log(pf);
    %     dGf0 = -R*temp*maxstar(-gpfnsp/(R*temp));
    end

    %isomer group thermodynamic standard transformed Enthalpy of
    %formation for a reactant
    %%%%%%% makes script faster to leave this out for now.
    dHzero=[];
    %%%%%%%
    if ~isempty(dHzero)
        %make temperature a smaller variable
        t=temp;
        switch length(dGzero)
            case 1
                %corresponds to Simplify[-(t^2*D[dGf0/t, t])] in Albertys code for
                %one species reactant
                %see Mathematica file dHfn.nb
                A=dGzero(1);
                B=dHzero(1);
                C=zi(1);
                D=nH(1);
                dHf0 =(B*(1+1.6*is^0.5)*s*v + (C^2)*(is^0.5)*(t^2)*(2*s*t*u-r*v)+D*(is^0.5)*(t^2)*(-2*s*t*u+r*v)) / ((1+1.6*is^0.5)*s*v);
            case 2
                %see Mathematica file dHfn.nb
                A=dGzero(1);
                B=dHzero(1);
                C=zi(1);
                D=nH(1);
                a=dGzero(2);
                b=dHzero(2);
                c=zi(2);
                d=nH(2);
                %translated to matlab from mathematica by hand
                dHf0 =((10^-pHr)^d*exp(-(b*((1/t)-(1/T))+a/T-(((c^2-d)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R)*(b*(1+1.6*is^0.5)*s*v+...
                    c^2*is^0.5*t^2*(2*s*t*u-r*v)+d*is^0.5*t^2*(-2*s*t*u+r*v))+...
                    (10^-pHr)^D*exp(-(B*((1/t)-(1/T))+A/T-(((C^2-D)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R)*(B*(1+1.6*is^0.5)*s*v+...
                    C^2*is^0.5*t^2*(2*s*t*u-r*v)+D*is^0.5*t^2*(-2*s*t*u+r*v)))/...
                    (((10^-pHr)^d*exp((b*((-1/t)+(1/T))-a/T+(((c^2-d)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R)+...
                    (10^-pHr)^D*exp((B*((-1/t)+(1/T))-A/T+(((C^2-D)*is^0.5*(p*s*v+q*t*(s*t*u-r*v)))/((1+1.6*is^0.5)*q*s*v)))/R))*(1+1.6*is^0.5)*s*v);
                %         Mathematica Expression to Matlab m-file Converter by Harri Ojanen, Rutgers University
                %         dHfn2=((10.^((-1).*pHr)).^d.*exp(R.^(-1).*(b.*((-1).*t.^(-1)+T.^(-1))+(-1).*a.* ...
                %             T.^(-1)+0.1E1.*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).* ...
                %             is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).* ...
                %             r.*v))))+(10.^((-1).*pHr)).^D.*exp(R.^(-1).*(B.*((-1).*t.^(-1)+T.^(-1))+( ...
                %             -1).*A.*T.^(-1)+0.1E1.*(C.^2+(-0.1E1).*D).*(0.1E1+0.16E1.*is.^0.5E0).^( ...
                %             -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                %             -0.1E1).*r.*v))))).^(-1).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).*s.^(-1).*v.^( ...
                %             -1).*((10.^((-1).*pHr)).^d.*exp((-1).*R.^(-1).*(b.*(t.^(-1)+(-1).*T.^(-1) ...
                %             )+a.*T.^(-1)+(-0.1E1).*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.*is.^0.5E0).^( ...
                %             -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                %             -0.1E1).*r.*v)))).*(b.*(1+0.16E1.*is.^0.5E0).*s.*v+c.^2.*is.^0.5E0.* ...
                %             t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+d.*is.^0.5E0.*t.^2.*((-0.2E1).*s.* ...
                %             t.*u+r.*v))+(10.^((-1).*pHr)).^D.*exp((-1).*R.^(-1).*(B.*(t.^(-1)+(-1).* ...
                %             T.^(-1))+A.*T.^(-1)+(-0.1E1).*(C.^2+(-0.1E1).*D).*(0.1E1+0.16E1.* ...
                %             is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                %             s.*t.*u+(-0.1E1).*r.*v)))).*(B.*(1+0.16E1.*is.^0.5E0).*s.*v+C.^2.* ...
                %             is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+D.*is.^0.5E0.*t.^2.*(( ...
                %             -0.2E1).*s.*t.*u+r.*v)));
            case 3
                A=dGzero(1);
                B=dHzero(1);
                C=zi(1);
                D=nH(1);
                a=dGzero(2);
                b=dHzero(2);
                c=zi(2);
                d=nH(2);
                e=dGzero(3);
                f=dHzero(3);
                g=zi(3);
                h=nH(3);
                %see Mathematica file dHfn.nb
                %Mathematica Expression to Matlab m-file Converter by Harri Ojanen, Rutgers University
                dHf0 = ((10.^((-1).*pHr)).^d.*exp(R.^(-1).*(b.*((-1).*t.^(-1)+T.^(-1))+(-1).*a.* ...
                    T.^(-1)+0.1E1.*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).* ...
                    is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).* ...
                    r.*v))))+(10.^((-1).*pHr)).^D.*exp(R.^(-1).*(B.*((-1).*t.^(-1)+T.^(-1))+( ...
                    -1).*A.*T.^(-1)+0.1E1.*(C.^2+(-0.1E1).*D).*(0.1E1+0.16E1.*is.^0.5E0).^( ...
                    -1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*(s.*t.*u+( ...
                    -0.1E1).*r.*v))))+(10.^((-1).*pHr)).^h.*exp(R.^(-1).*(f.*((-1).*t.^(-1)+ ...
                    T.^(-1))+(-1).*e.*T.^(-1)+0.1E1.*(g.^2+(-0.1E1).*h).*(0.1E1+0.16E1.* ...
                    is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                    s.*t.*u+(-0.1E1).*r.*v))))).^(-1).*(0.1E1+0.16E1.*is.^0.5E0).^(-1).*s.^( ...
                    -1).*v.^(-1).*((10.^((-1).*pHr)).^d.*exp((-1).*R.^(-1).*(b.*(t.^(-1)+(-1) ...
                    .*T.^(-1))+a.*T.^(-1)+(-0.1E1).*(c.^2+(-0.1E1).*d).*(0.1E1+0.16E1.* ...
                    is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.*v+q.*t.*( ...
                    s.*t.*u+(-0.1E1).*r.*v)))).*(b.*(1+0.16E1.*is.^0.5E0).*s.*v+c.^2.* ...
                    is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+d.*is.^0.5E0.*t.^2.*(( ...
                    -0.2E1).*s.*t.*u+r.*v))+(10.^((-1).*pHr)).^D.*exp((-1).*R.^(-1).*(B.*( ...
                    t.^(-1)+(-1).*T.^(-1))+A.*T.^(-1)+(-0.1E1).*(C.^2+(-0.1E1).*D).*(0.1E1+ ...
                    0.16E1.*is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*(p.*s.* ...
                    v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v)))).*(B.*(1+0.16E1.*is.^0.5E0).*s.*v+ ...
                    C.^2.*is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+D.*is.^0.5E0.* ...
                    t.^2.*((-0.2E1).*s.*t.*u+r.*v))+(10.^((-1).*pHr)).^h.*exp((-1).*R.^(-1).* ...
                    (f.*(t.^(-1)+(-1).*T.^(-1))+e.*T.^(-1)+(-0.1E1).*(g.^2+(-0.1E1).*h).*( ...
                    0.1E1+0.16E1.*is.^0.5E0).^(-1).*is.^0.5E0.*q.^(-1).*s.^(-1).*v.^(-1).*( ...
                    p.*s.*v+q.*t.*(s.*t.*u+(-0.1E1).*r.*v)))).*(f.*(1+0.16E1.*is.^0.5E0).* ...
                    s.*v+g.^2.*is.^0.5E0.*t.^2.*(0.2E1.*s.*t.*u+(-0.1E1).*r.*v)+h.* ...
                    is.^0.5E0.*t.^2.*((-0.2E1).*s.*t.*u+r.*v)));
            otherwise
                dHf0 = NaN;
        end
        if isnan(dHf0)
            %see Mathematica file dHfn.nb
%             error('Too many species in reactant so standard transformed enthalpy is NaN')
            fprintf('%s\n','Too many species in reactant so standard transformed enthalpy is NaN')
        end
        %cast back into normal matlab double
        dHf0=double(dHf0);
    else
        %changed to NaN 20 July 2009
        dHf0=NaN;
    end
    %cast back into normal matlab double
    dGf0=double(dGf0);
    gpfnsp=double(gpfnsp);
end

%
%This is a helper function to compute the log(sum(exp(v))) of a vector v
function y = maxstar(x, w, dim)
% maxstar   Log of a sum of exponentials.
%   For vectors, maxstar(x) is equivalent to log(sum(exp(x))).
%   For matrices, maxstar(x) is a row vector and maxstar operates on
%   each column of x. For N-D arrays, maxstar(x) operates along the
%   first non-singleton dimension.
%
%   maxstar(x,w) is the log of a weighted sum of exponentials,
%   equivalent to log(sum(w.*exp(x))). Vectors w and x must be
%   the same length. For matrix x, the weights w can be input as
%   a matrix the same size as x, or as a vector of the same length
%   as columns of x. Weights may be zero or negative, but the result
%   sum(w.*exp(x)) must be greater than zero.
%
%   maxstar(x, [], dim) operates along the dimension dim, and has
%   the same dimensions as the MATLAB function max(x, [], dim).
%
%   Note:
%   The max* function is described in Lin & Costello, Error Control
%   Coding, 2nd Edition, equation 12.127, in the two-argument form
%     max*(x1,x2) = max(x1,x2) + log(1 + exp(-abs(x1-x2))).
%   The function max* can be applied iteratively:
%     max*(x1,x2,x3) = max*(max*(x1,x2),x3).
%   Functions max(x) ~ max*(x), and min(x) ~ -max*(-x).
%
%   Algorithm:
%   The double precision MATLAB expresson log(sum(exp(x))) fails
%   if all(x < -745), or if any(x > 706). This is avoided using
%   m = max(x) in  max*(x) = m + log(sum(exp(x - m))).
%
%   Example: If x = [2 8 4
%                    7 3 9]
%
%   then maxstar(x,[],1) is [7.0067 8.0067 9.0067],
%
%   and  maxstar(x,[],2) is [8.0206
%                            9.1291].

% 2006-02-10   R. Dickson
% 2006-03-25   Implemented N-D array features following a suggestion
%              from John D'Errico.
%
%   Uses: max, log, exp, sum, shiftdim, repmat, size, zeros, ones,
%         length, isempty, error, nargin, find, reshape

if nargin < 1 || nargin > 3
    error('Wrong number of input arguments.');
end

[x, n] = shiftdim(x);
szx = size(x);

switch nargin
    case 1
        w = [];
        dim = 1;
    case 2
        dim = 1;
    case 3
        dim = dim - n;
end

if isempty(w)
    % replicate m = max(x) to get mm, with size(mm) == size(x)
    m = max(x,[],dim);
    szm = ones(size(szx));
    szm(dim) = szx(dim);
    mm = repmat(m,szm);
    y = m + log(sum(exp(x - mm), dim));
else
    w = shiftdim(w);
    szw = size(w);
    % protect the second condition with a short-circuit or
    if ~(length(szw) == length(szx)) || ~all(szw == szx)
        if size(w,1) == size(x,dim)
            % replicate w with repmat so size(w) == size(x)
            szw = ones(size(szx));
            szw(dim) = size(w,1);
            w = reshape(w, szw);
            szr = szx;
            szr(dim) = 1;
            w = repmat(w, szr);
        else
            error('Length of w must match size(x,dim).');
        end
    end

    % Move the weight into the exponent xw and find
    % m = max(xw) over terms with positive weights
    ipos = find(w>0);
    xw = -Inf*zeros(szx);
    xw(ipos) = x(ipos) + log(w(ipos));
    m = max(xw,[],dim);
    % replicate m with repmat so size(mm) == size(x)
    szm = ones(size(szx));
    szm(dim) = szx(dim);
    mm = repmat(m,szm);
    exwp = zeros(szx);
    exwp(ipos) = exp(xw(ipos)-mm(ipos));
    % check for terms with negative weights
    ineg = find(w<0);
    if ~isempty(ineg)
        exwn = zeros(szx);
        exwn(ineg) = exp(x(ineg) + log(-w(ineg)) - mm(ineg));
        y = m + log(sum(exwp, dim) - sum(exwn, dim));
    else
        y = m + log(sum(exwp, dim));
    end
end
%

author: opencobra

repo: cobratoolbox


text:
function [g1, g2] = lfmwhiteXrbfwhiteKernGradient(lfmKern, rbfKern, t1, varargin)

% LFMWHITEXRBFWHITEKERNGRADIENT Compute a cross gradient between an LFM-WHITE
% and an RBF-WHITE kernels.
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel between an
% LFM-WHITE and an RBF-WHITE kernels for the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFM-WHITE kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE kernel.
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the LFM-WHITE kernel, for
% ordering see lfmwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the RBF-WHITE kernel, for
% ordering see rbfwhiteKernExtractParam.
%
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel between an
% LFM-WHITE and an RBF-WHITE kernels for the multiple output kernel.
% ARG lfmKern : the kernel structure associated with the LFM-WHITE kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the LFM-WHITE kernel, for
% ordering see lfmwhiteKernExtractParam.
% RETURN g2 : gradient of the parameters of the RBF-WHITE kernel, for
% ordering see rbfwhiteKernExtractParam.
%
% SEEALSO : multiKernParamInit, multiKernCompute, lfmwhiteKernParamInit,
% rbfwhiteKernParamInit, lfmwhiteKernExtractParam, rbfwhiteKernExtractParam
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 5
    t2 = t1;
else
    t2 = varargin{1};
end
covGrad = varargin{end};

if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if lfmKern.variance ~= rbfKern.variance
  error('Kernels cannot be cross combined if they have different variances.')
end
if lfmKern.isStationary ~= rbfKern.isStationary
  error('Stationary and non-stationary kernels cannot be cross combined.')
end

g1 = zeros(1, lfmKern.nParams);
g2 = zeros(1, rbfKern.nParams);

T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
deltaT = T1 - T2;
indT = double(T1 >= T2);

% Parameters required in the computation of the kernel
isStationary = lfmKern.isStationary;
mass = lfmKern.mass;
spring = lfmKern.spring;
damper = lfmKern.damper;
variance = lfmKern.variance;
sensitivity = lfmKern.sensitivity;

invWidth = rbfKern.inverseWidth;

alpha = lfmKern.alpha;
omega = lfmKern.omega;
gamma = lfmKern.gamma;
gammaTilde = alpha - j*omega;

% Gradients of theta w.r.t. mass, spring and damper

gradThetaMass = [1 0 0];
gradThetaAlpha = [-damper/(2*(mass^2)) 0 1/(2*mass)];
gradThetaOmega = [(damper^2-2*mass*spring)/(2*(mass^2)*sqrt(4*mass*spring-damper^2)) ...
    1/(sqrt(4*mass*spring-damper^2)) -damper/(2*mass*sqrt(4*mass*spring-damper^2))];
gradThetaGamma = gradThetaAlpha + j*gradThetaOmega;
gradThetaGammaTilde = gradThetaAlpha - j*gradThetaOmega;

% Computation of the normalised kernel (i.e. variance = 1 and sensitivity = 1)
% and auxiliary functions and constants

lfmKern.variance = 1;
lfmKern.sensitivity = 1;
rbfKern.variance = 1;
K = lfmwhiteXrbfwhiteKernCompute(lfmKern, rbfKern, t1, t2);

c = 1 / (j*4*mass*omega);

varphiT1T2 = gamma * deltaT .* indT + 0.5 * invWidth * (deltaT.^2) .* (1-indT);
zT1T2 = sqrt(0.5*invWidth) * (-deltaT .* (1-indT) + gamma / invWidth);
varphiT1T2Tilde = gammaTilde * deltaT .* indT + 0.5 * invWidth * (deltaT.^2) .* (1-indT);
zT1T2Tilde = sqrt(0.5*invWidth) * (-deltaT .* (1-indT) + gammaTilde / invWidth);
if (isStationary == false)
    varphiT10 = gamma * T1;
    varphi0T2 = 0.5 * invWidth * (T2.^2);
    z0T2 = sqrt(0.5*invWidth) * (T2 + gamma/invWidth);
    varphiT10Tilde = gammaTilde * T1;
    z0T2Tilde = sqrt(0.5*invWidth) * (T2 + gammaTilde/invWidth);
    %gradThetaVarphiT10 = T1;
end

gradThetaVarphiT1T2 = deltaT .* indT;
%gradThetaZ = 1/sqrt(2*invWidth);

% Gradient w.r.t. the mass, spring and damper (lfmKern)
for i = 1:3
    gradThetaPsi = - exp(-varphiT1T2) ...
        .* (gradThetaGamma(i) * wofzHui(j*zT1T2) .* gradThetaVarphiT1T2 ...
        + sqrt(2/invWidth) * gradThetaGamma(i) * (1/sqrt(pi)-zT1T2.*wofzHui(j*zT1T2)));
    gradThetaPsiTilde = - exp(-varphiT1T2Tilde) ...
        .* (gradThetaGammaTilde(i) * wofzHui(j*zT1T2Tilde) .* gradThetaVarphiT1T2 ...
        + sqrt(2/invWidth) * gradThetaGammaTilde(i) * (1/sqrt(pi)-zT1T2Tilde.*wofzHui(j*zT1T2Tilde)));
    if (isStationary == false)
        gradThetaPsi = gradThetaPsi + exp(-(varphiT10 + varphi0T2)) ...
            .* (gradThetaGamma(i) * T1 .* wofzhui(j*z0T2) ...
            + sqrt(2/invWidth) * gradThetaGamma(i) * (1/sqrt(pi)-z0T2.*wofzHui(j*z0T2)));
        gradThetaPsiTilde = gradThetaPsiTilde + exp(-(varphiT10Tilde + varphi0T2)) ...
            .* (gradThetaGammaTilde(i) * T1 .* wofzhui(j*z0T2Tilde) ...
            + sqrt(2/invWidth) * gradThetaGammaTilde(i) * (1/sqrt(pi)-z0T2Tilde.*wofzHui(j*z0T2Tilde)));
    end
    g1(i) = sensitivity * variance * sum(sum((-(gradThetaMass(i)/mass + gradThetaOmega(i)/omega) ...
        * K + c * (gradThetaPsiTilde - gradThetaPsi)) .* covGrad));
end

% Gradient w.r.t. the inverse width (rbfKern)
gradInvWidthVarphiT1T2 = 0.5 * (deltaT.^2) .* (1-indT);
gradInvWidthZT1T2 = -(deltaT .* (1-indT) + gamma / invWidth) / sqrt(8*invWidth);
gradInvWidthZT1T2Tilde = -(deltaT .* (1-indT) + gammaTilde / invWidth) / sqrt(8*invWidth);
gradInvWidthPsi = - exp(-varphiT1T2) ...
    .* (wofzHui(j*zT1T2) .* gradInvWidthVarphiT1T2 ...
    + 2 * (1/sqrt(pi)-zT1T2.*wofzHui(j*zT1T2)) .* gradInvWidthZT1T2);
gradInvWidthPsiTilde = - exp(-varphiT1T2Tilde) ...
    .* (wofzHui(j*zT1T2Tilde) .* gradInvWidthVarphiT1T2 ...
    + 2 * (1/sqrt(pi)-zT1T2Tilde.*wofzHui(j*zT1T2Tilde)) .* gradInvWidthZT1T2Tilde);
if (isStationary == false)
    gradInvWidthVarphi0T2 = 0.5*(T2.^2);
    gradInvWidthZ0T2 = (T2 - gamma / invWidth) / sqrt(8*invWidth);
    gradInvWidthZ0T2Tilde = (T2 - gammaTilde / invWidth) / sqrt(8*invWidth);
    gradInvWidthPsi = gradInvWidthPsi + exp(-(varphiT10 + varphi0T2)) ...
        .* (gradInvWidthVarphi0T2 .* wofzhui(j*z0T2) ...
        + 2 * (1/sqrt(pi)-z0T2.*wofzHui(j*z0T2)) .* gradInvWidthZ0T2);
    gradInvWidthPsiTilde = gradInvWidthPsiTilde + exp(-(varphiT10Tilde + varphi0T2)) ...
        .* (gradInvWidthVarphi0T2 .* wofzhui(j*z0T2Tilde) ...
        + 2 * (1/sqrt(pi)-z0T2Tilde.*wofzHui(j*z0T2Tilde)) .* gradInvWidthZ0T2Tilde);
end
g2(1) = sensitivity * variance * sum(sum((c * (gradInvWidthPsiTilde - gradInvWidthPsi)) .* covGrad));

% Gradient w.r.t. sigma_r^2
g1(4) = sensitivity * sum(sum(K .* covGrad));
g2(2) = 0; % Otherwise it is counted twice

% Gradient w.r.t. sensitivity (only lfmKern)
g1(5) = variance * sum(sum(K .* covGrad));

% Ensure that the gradients are real
g1 = real(g1);
g2 = real(g2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% first casadi test for mpc fpr mobile robots
clear all
close all
clc

% CasADi v3.4.5
% addpath('C:\Users\mehre\OneDrive\Desktop\CasADi\casadi-windows-matlabR2016a-v3.4.5')
% CasADi v3.5.5
addpath('C:\Users\mehre\OneDrive\Desktop\CasADi\casadi-windows-matlabR2016a-v3.5.5')

import casadi.*

T = 0.2; %[s]
N = 10; % prediction horizon
rob_diam = 0.3;

v_max = 0.6; v_min = -v_max;
omega_max = pi/4; omega_min = -omega_max;

x = SX.sym('x'); y = SX.sym('y'); theta = SX.sym('theta');
states = [x;y;theta]; n_states = length(states);

v = SX.sym('v'); omega = SX.sym('omega');
controls = [v;omega]; n_controls = length(controls);
rhs = [v*cos(theta);v*sin(theta);omega]; % system r.h.s

f = Function('f',{states,controls},{rhs}); % nonlinear mapping function f(x,u)
U = SX.sym('U',n_controls,N); % Decision variables (controls)
P = SX.sym('P',n_states + n_states);
% parameters (which include the initial state and the reference state)

X = SX.sym('X',n_states,(N+1));
% A vector that represents the states over the optimization problem.

obj = 0; % Objective function
g = [];  % constraints vector

Q = zeros(3,3); Q(1,1) = 1;Q(2,2) = 5;Q(3,3) = 0.1; % weighing matrices (states)
R = zeros(2,2); R(1,1) = 0.5; R(2,2) = 0.05; % weighing matrices (controls)

st  = X(:,1); % initial state
g = [g;st-P(1:3)]; % initial condition constraints
for k = 1:N
    st = X(:,k);  con = U(:,k);
    obj = obj+(st-P(4:6))'*Q*(st-P(4:6)) + con'*R*con; % calculate obj
    st_next = X(:,k+1);
    f_value = f(st,con);
    st_next_euler = st+ (T*f_value);
    g = [g;st_next-st_next_euler]; % compute constraints
end
% make the decision variable one column  vector
OPT_variables = [reshape(X,3*(N+1),1);reshape(U,2*N,1)];

nlp_prob = struct('f', obj, 'x', OPT_variables, 'g', g, 'p', P);

opts = struct;
opts.ipopt.max_iter = 2000;
opts.ipopt.print_level =0;%0,3
opts.print_time = 0;
opts.ipopt.acceptable_tol =1e-8;
opts.ipopt.acceptable_obj_change_tol = 1e-6;

solver = nlpsol('solver', 'ipopt', nlp_prob,opts);

args = struct;

args.lbg(1:3*(N+1)) = 0;  % -1e-20  % Equality constraints
args.ubg(1:3*(N+1)) = 0;  % 1e-20   % Equality constraints

args.lbx(1:3:3*(N+1),1) = -2; %state x lower bound
args.ubx(1:3:3*(N+1),1) = 2; %state x upper bound
args.lbx(2:3:3*(N+1),1) = -2; %state y lower bound
args.ubx(2:3:3*(N+1),1) = 2; %state y upper bound
args.lbx(3:3:3*(N+1),1) = -inf; %state theta lower bound
args.ubx(3:3:3*(N+1),1) = inf; %state theta upper bound

args.lbx(3*(N+1)+1:2:3*(N+1)+2*N,1) = v_min; %v lower bound
args.ubx(3*(N+1)+1:2:3*(N+1)+2*N,1) = v_max; %v upper bound
args.lbx(3*(N+1)+2:2:3*(N+1)+2*N,1) = omega_min; %omega lower bound
args.ubx(3*(N+1)+2:2:3*(N+1)+2*N,1) = omega_max; %omega upper bound
%----------------------------------------------
% ALL OF THE ABOVE IS JUST A PROBLEM SET UP

% THE SIMULATION LOOP SHOULD START FROM HERE
%-------------------------------------------
t0 = 0;
x0 = [0 ; 0 ; 0.0];    % initial condition.
xs = [1.5 ; 1.5 ; 0.0]; % Reference posture.

xx(:,1) = x0; % xx contains the history of states
t(1) = t0;

u0 = zeros(N,2);        % two control inputs for each robot
X0 = repmat(x0,1,N+1)'; % initialization of the states decision variables

sim_tim = 20; % Maximum simulation time

% Start MPC
mpciter = 0;
xx1 = [];
u_cl=[];

% the main simulaton loop... it works as long as the error is greater
% than 10^-6 and the number of mpc steps is less than its maximum
% value.
tic
while(norm((x0-xs),2) > 1e-2 && mpciter < sim_tim / T)
    args.p   = [x0;xs]; % set the values of the parameters vector
    % initial value of the optimization variables
    args.x0  = [reshape(X0',3*(N+1),1);reshape(u0',2*N,1)];
    sol = solver('x0', args.x0, 'lbx', args.lbx, 'ubx', args.ubx,...
        'lbg', args.lbg, 'ubg', args.ubg,'p',args.p);
    u = reshape(full(sol.x(3*(N+1)+1:end))',2,N)'; % get controls only from the solution
    xx1(:,1:3,mpciter+1)= reshape(full(sol.x(1:3*(N+1)))',3,N+1)'; % get solution TRAJECTORY
    u_cl= [u_cl ; u(1,:)];
    t(mpciter+1) = t0;
    % Apply the control and shift the solution
    [t0, x0, u0] = shift(T, t0, x0, u,f);
    xx(:,mpciter+2) = x0;
    X0 = reshape(full(sol.x(1:3*(N+1)))',3,N+1)'; % get solution TRAJECTORY
    % Shift trajectory to initialize the next step
    X0 = [X0(2:end,:);X0(end,:)];
    mpciter
    mpciter = mpciter + 1;
end;
toc

ss_error = norm((x0-xs),2)
Draw_MPC_point_stabilization_v1 (t,xx,xx1,u_cl,xs,N,rob_diam)

author: MMehrez


text:
function calpak_test616 ( )

%*****************************************************************************80
%
%% CALPAK_TEST616 tests YEAR_LENGTH_ROMAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST616\n' );
  fprintf ( 1, '  For the Roman calendar:\n' );
  fprintf ( 1, '  YEAR_LENGTH_ROMAN determines the length of a year.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Year  Length\n' );
  fprintf ( 1, '\n' );

  for y = 96 : 100
    sy = y_to_s_roman ( y );
    fprintf ( 1, '  %10s  %3d\n', sy, year_length_roman ( y ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [zValues,zCosts,zGuesses,inConvHull,meanMax,exitFlags] = ...
    findTDistributedProjections_fmin(data,trainingData,...
                                        trainingEmbedding,parameters)
%findTDistributedProjections_fmin is called by runEmbbedings to find the
%optimal embeddings of a set of wavelet amplitudes into a previously
%defined t-SNE embedding

    readout = 5000;

    perplexity = parameters.perplexity;
    sigmaTolerance = parameters.sigmaTolerance;
    maxNeighbors = parameters.maxNeighbors;
    batchSize = parameters.embedding_batchSize;
    maxOptimIter = parameters.maxOptimIter;

    N = length(data(:,1));
    zValues = zeros(N,2);
    zGuesses = zeros(N,2);
    zCosts = zeros(N,1);
    batches = ceil(N/batchSize);
    inConvHull = false(N,1);
    meanMax = zeros(N,1);
    exitFlags = zeros(N,1);

    options = optimset('Display','off','maxiter',maxOptimIter);

    for j=1:batches
        fprintf(1,'\t Processing batch #%4i out of %4i\n',j,batches);
        idx = (1:batchSize) + (j-1)*batchSize;
        idx = idx(idx <= N);
        current_guesses = zeros(length(idx),2);
        current = zeros(length(idx),2);
        currentData = data(idx,:);
        tCosts = zeros(size(idx));
        current_poly = false(length(idx),1);

        D2 = findListKLDivergences(currentData,trainingData);
        current_meanMax = zeros(length(idx),1);

        parfor i=1:length(idx)

            if mod(i,readout) == 0
                fprintf(1,'\t\t Image #%5i\n',i);
            end

            [~,p] = returnCorrectSigma_sparse(D2(i,:),perplexity,sigmaTolerance,maxNeighbors);
            idx2 = p>0;
            z = trainingEmbedding(idx2,:);
            [~,maxIdx] = max(p);
            a = sum(bsxfun(@times,z,p(idx2)'));

            guesses = [a;trainingEmbedding(maxIdx,:)];

            b = zeros(2,2);
            c = zeros(2,1)
            flags = zeros(2,1);

            q = convhull(z);
            q = z(q,:);

            [b(1,:),c(1),flags(1)] = fminsearch(@(x)calculateKLCost(x,z,p(idx2)),guesses(1,:),options);
            [b(2,:),c(2),flags(2)] = fminsearch(@(x)calculateKLCost(x,z,p(idx2)),guesses(2,:),options);
            polyIn = inpolygon(b(:,1),b(:,2),q(:,1),q(:,2));

            if sum(polyIn) > 0
                pp = find(polyIn);
                [~,mI] = min(c(polyIn));
                mI = pp(mI);
                current_poly(i) = true;
            else
                [~,mI] = min(c);
                current_poly(i) = false;
            end

            exitFlags(i) = flags(mI);
            current_guesses(i,:) = guesses(mI,:);
            current(i,:) = b(mI,:);
            tCosts(i) = c(mI);
            current_meanMax(i) = mI;

        end

        zGuesses(idx,:) = current_guesses;
        zValues(idx,:) = current;
        zCosts(idx) = tCosts;
        inConvHull(idx) = current_poly;
        meanMax(idx) = current_meanMax;

    end

    zValues(~inConvHull,:) = zGuesses(~inConvHull,:);

author: gordonberman

repo: MotionMapper


text:
function varargout = process_mri_deface( varargin )
% PROCESS_MRI_DEFACE: Remove the facial features from an MRI.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2019-2023
%          Inspired from SPM12 function spm_deface

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Deface MRI volumes';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import anatomy'};
    sProcess.Index       = 41;
    sProcess.Description = 'https://surfer.nmr.mgh.harvard.edu/fswiki/mri_deface';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = 'NewSubject';
    % Option: Method
    sProcess.options.method_title.Comment = '
Defacing method:';
    sProcess.options.method_title.Type    = 'label';
    sProcess.options.method.Comment = {'FreeSurfer: mri_deface', 'BrainSuite: Remove the face with a pre-defined mask', 'SPM: Cut below a plane in MNI coordinates (a.x+b.y+c.z+d=0)'; ...
                                       'freesurfer', 'brainsuite', 'spm'};
    sProcess.options.method.Type    = 'radio_label';
    sProcess.options.method.Value   = 'spm';
    % SPM: MNI plane coordinates
    sProcess.options.mniplane.Comment = '      [a, b, c, d]: ';
    sProcess.options.mniplane.Type    = 'value';
    sProcess.options.mniplane.Value   = {[0, -11, 9.6, 1], 'list', 3};
    sProcess.options.mniplanedef.Comment = '      Cut face: a=0, b=-11, c=9.6, d=1
      Cut neck: a=0, b=0, c=1, d=0.115 (Z<-115)';
    sProcess.options.mniplanedef.Type    = 'label';
    % Option: Recompute head surface
    sProcess.options.defacehead.Comment = 'Recompute head surface';
    sProcess.options.defacehead.Type    = 'checkbox';
    sProcess.options.defacehead.Value   = 1;
    sProcess.options.defacehead.Group   = 'output';
    % Option: Overwrite
    sProcess.options.overwrite.Comment = 'Overwrite existing files';
    sProcess.options.overwrite.Type    = 'checkbox';
    sProcess.options.overwrite.Value   = 0;
    sProcess.options.overwrite.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};

    % ===== GET OPTIONS =====
    OPTIONS.isOverwrite = sProcess.options.overwrite.Value;
    OPTIONS.Method = sProcess.options.method.Value;
    OPTIONS.isInteractive = 0;
    OPTIONS.isDefaceHead = sProcess.options.defacehead.Value;
    if (length(sProcess.options.mniplane.Value{1}) == 4)
        OPTIONS.MNIplane = sProcess.options.mniplane.Value{1};
    else
        bst_report('Error', sProcess, [], 'Invalid MNI plane equation (must be 4 values: a b c d).');
        return;
    end
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return;
    end

    % ===== GET SUBJECT =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    if isempty(iSubject)
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" does not exist.']);
        return
    end
    % Check if a MRI is available for the subject
    if isempty(sSubject.Anatomy)
        bst_report('Error', sProcess, [], ['No MRI available for subject "' SubjectName '".']);
        return
    end

    % ===== DEFACE ALL VOLUMES =====
    % Deface volumes
    [DefacedMri, errMsg] = Compute(iSubject, OPTIONS);
    % Error handling
    if isempty(DefacedMri)
        bst_report('Error', sProcess, [], errMsg);
        return;
    elseif ~isempty(errMsg)
        bst_report('Warning', sProcess, [], errMsg);
    end

    OutputFiles = {'import'};
end

%% ===== EXTERNAL CALL =====
% USAGE:  [DefacedFiles, errMsg] = Compute(MriFiles, OPTIONS=[])
%         [DefacedFiles, errMsg] = Compute(iSubject, OPTIONS=[])
function [DefacedFiles, errMsg] = Compute(MriFiles, OPTIONS)
    % Parse inputs
    if ischar(MriFiles)
        MriFiles = {MriFiles};
    elseif isnumeric(MriFiles)
        iSubject = MriFiles;
        % Get volumes to deface
        sSubject = bst_get('Subject', iSubject);
        iNoAtlas = find(cellfun(@(c)isempty(strfind(c, '_volatlas')), {sSubject.Anatomy.FileName}));
        MriFiles = {sSubject.Anatomy(iNoAtlas).FileName};
        % Put the default volume first
        if ~isempty(sSubject.iAnatomy) && ismember(sSubject.iAnatomy, iNoAtlas)
            iOrder = [sSubject.iAnatomy, setdiff(1:length(iNoAtlas), sSubject.iAnatomy)];
            MriFiles = MriFiles(iOrder);
        end
    end
    if (nargin < 2) || isempty(OPTIONS)
        OPTIONS = struct();
    end
    % Ask method to user in interactive mode
    if (~isfield(OPTIONS, 'Method') || isempty(OPTIONS.Method)) && (~isfield(OPTIONS, 'isInteractive') || OPTIONS.isInteractive)
        % Windows: no options anyway (FreeSurfer's mri_deface not available)
        if ismember(bst_get('OsType',0), {'win32', 'win64'})
            OPTIONS.Method = 'spm';
        % Other OS: Ask for choice
        else
            res = java_dialog('combo', 'Select the defacing method:
', 'Deface volume', [], {'spm', 'freesurfer'}, 'freesurfer');
            if isempty(res)
                errMsg = 'Aborted by user.';
                return;
            end
            OPTIONS.Method = res;
        end
    end
    % Default options
    Def_OPTIONS = struct(...
        'Method',        'spm', ...
        'MNIplane',      [0, -11, 9.6, 1], ...
        'isOverwrite',   0, ...
        'isInteractive', 1, ...
        'isDefaceHead',  1);
    OPTIONS = struct_copy_fields(OPTIONS, Def_OPTIONS, 0);
    % Initialize returned variables
    DefacedFiles = {};
    errMsg = [];
    fileTag = ' | deface';

    % Progress bar
    isProgress = bst_progress('isVisible');
    if ~isProgress
        bst_progress('start', 'Deface MRI', 'Initialization...');
    end
    % Process each input file separately
    MriHead = [];
    for iFile = 1:length(MriFiles)
        bst_progress('text', 'Loading input MRI...');
        % Get subject index
        [sSubject, iSubject, iAnatomy] = bst_get('MriFile', MriFiles{iFile});
        % If MRI was already defaced: skip
        if (length(MriFiles) > 1) && ~isempty(strfind(sSubject.Anatomy(iAnatomy).Comment, fileTag))
            continue;
        end
        % Check if it is loaded in memory
        sMri = bst_memory('GetMri', MriFiles{iFile});
        % If not: load it from the file
        if isempty(sMri)
            sMri = in_mri_bst(MriFiles{iFile});
        end

        % Switch depending on the method
        switch (OPTIONS.Method)
            case 'spm'
                % Compute cut mask on first volume only
                if (iFile == 1)
                    [sMri, maskCut, errMsg] = CutMriPlane(sMri, OPTIONS.MNIplane);
                    % Error management (if MNI normalization failed)
                    if ~isempty(errMsg)
                        if ~isProgress
                            bst_progress('stop');
                        end
                        return;
                    end
                % Following volumes: just a apply the same cut
                else
                    sMri.Cube(maskCut) = 0;
                end

            case 'freesurfer'
                % Get path to mri_deface (download if necessary)
                [exePath, talFile, faceFile, errMsg] = InstallMriDeface(OPTIONS.isInteractive);
                if ~isempty(errMsg)
                    if ~isProgress
                        bst_progress('stop');
                    end
                    return;
                end
                disp(['BST> Deface: Using ' exePath]);
                % Save temporary MRI file
                bst_progress('text', 'Saving temporary nii file...');
                TmpDir = bst_get('BrainstormTmpDir', 0, 'deface');
                fileNii = bst_fullfile(TmpDir, 'orig.nii');
                out_mri_nii(sMri, fileNii, 'int16');
                % Change current folder to .brainstorm/tmp, because mri_deface saves a useless log in the current folder
                curdir = pwd;
                cd(TmpDir);
                % Call mri_deface
                bst_progress('text', 'Running mri_deface...');
                fileNiiDefaced = bst_fullfile(TmpDir, 'orig_defaced.nii');
                cmdDeface = [exePath ' ' fileNii ' ' talFile ' ' faceFile ' ' fileNiiDefaced];
                status = system(cmdDeface);
                if (status ~= 0)
                    errMsg = 'Error calling mri_deface (see console for details).';
                    if ~isProgress
                        bst_progress('stop');
                    end
                    return;
                end
                % Restore initial folder
                cd(curdir);
                % Read defaced file
                bst_progress('text', 'Reading defaced file...');
                sMriDefaced = in_mri_nii(fileNiiDefaced, 0, 0, 0);
                % Saves defaced volume
                sMri.Cube = sMriDefaced.Cube;
                % Delete the temporary files
                file_delete(TmpDir, 1, 1);

            case 'brainsuite'
                % Get ICBM152 defacing template
                sTemplate = bst_get('AnatomyDefaults', 'ICBM152');
                MaskFile = bst_fullfile(sTemplate.FilePath, 'facemask_300z.nii.gz');
                if isempty(MaskFile)
                    errMsg = ['Could not find face mask: ' MaskFile];
                    if ~isProgress
                        bst_progress('stop');
                    end
                    return;
                end
                % Load mask file
                sMask = in_mri(MaskFile, 'Nifti1');
                sMask.NCS.R = eye(3);
                sMask.NCS.T = [-99; -135; -184];
                % Reslice to the space of the MRI to deface
                sMaskReslice = mri_reslice(sMask, sMri, 'ncs', 'ncs');
                % Apply mask
                sMri.Cube(sMaskReslice.Cube == 0) = 0;

            otherwise
                errMsg = ['Invalid defacing method: ' OPTIONS.Method];
                if ~isProgress
                    bst_progress('stop');
                end
                return;
        end
        % Add comment tag
        sMri.Comment = [sMri.Comment, fileTag];

        % Remove initial file header
        sMri.Header = [];
        sMri.Histogram = [];
        % Remove file history (may contain information in the original file names)
        sMri.History = [];
        sMri = bst_history('add', sMri, 'process', ['process_mri_deface: ', OPTIONS.Method]);
        if strcmpi(OPTIONS.Method, 'spm')
            sMri = bst_history('add', sMri, 'process', sprintf('Cut below MNI plane: %1.3fx + %1.3fy + %1.3fz + %1.3f = 0', OPTIONS.MNIplane));
        end

        % Save defaced MRI
        bst_progress('text', 'Saving results to database...');
        if OPTIONS.isOverwrite
            % Update file structure
            bst_save(file_fullpath(MriFile), sMri, 'v6');
            DefacedFiles{end+1} = MriFile;
            % Unload from memory to force reloading
            bst_memory('UnloadMri', MriFile);
        % Add new file
        else
            DefacedFiles{end+1} = db_add(iSubject, sMri, 0);
            iAnatomy = length(sSubject.Anatomy) + 1;
        end
        % Update database registration
        sSubject.Anatomy(iAnatomy).FileName = DefacedFiles{end};
        sSubject.Anatomy(iAnatomy).Comment = sMri.Comment;
        bst_set('Subject', iSubject, sSubject);
        % Refresh tree
        panel_protocols('UpdateNode', 'Subject', iSubject);
        panel_protocols('SelectNode', [], 'subject', iSubject, -1);
        % If we recompute the head, it is on the defaced version of the first MRI in the list
        if (iFile == 1)
            MriHead = DefacedFiles{end};
        end
    end

    % ===== RECOMPUTE HEAD SURFACES =====
    if OPTIONS.isDefaceHead && ~isempty(MriHead)
        % Get subject
        [sSubject, iSubject] = bst_get('MriFile', MriHead);
        % Get all head surfaces
        [sSurfaces, iSurfaces] = bst_get('SurfaceFileByType', iSubject, 'Scalp');
        % Remove the BEM surfaces (too smooth to be problematic)
        if ~isempty(sSurfaces)
            iBem = find(~cellfun(@(c)isempty(strfind(c, '_bem_')), {sSurfaces.FileName}));
            if ~isempty(iBem)
                sSurfaces(iBem) = [];
                iSurfaces(iBem) = [];
            end
        end
        % If there are surfaces to delete
        if ~isempty(sSurfaces)
            % Delete head surfaces
            file_delete(file_fullpath({sSurfaces.FileName}), 1);
            % Update database structure
            HeadComment = [sSurfaces(1).Comment, fileTag];
            sSubject.Surface(iSurfaces) = [];
            bst_set('Subject', iSubject, sSubject);
            % Compute new head surface
            sSubject.Anatomy(sSubject.iAnatomy).FileName;
            tess_isohead(MriHead, 10000, 0, 2, HeadComment);
        end
    end

    % Close progress bar
    if ~isProgress
        bst_progress('stop');
    end
end

%% ===== CUT MRI PLANE =====
function [sMri, maskCut, errMsg] = CutMriPlane(sMri, MNIplane)
    % If the MNI normalization are not available: compute it now
    if (~isfield(sMri, 'NCS') || isempty(sMri.NCS) || ~isfield(sMri.NCS, 'R') || ~isfield(sMri.NCS, 'T') || isempty(sMri.NCS.R) || isempty(sMri.NCS.T))
        [sMri, errMsg] = bst_normalize_mni(sMri);
        if ~isempty(errMsg)
            maskCut = [];
            return;
        end
    else
        errMsg = '';
    end
    % Get MNI transformation
    vox2mni = cs_convert(sMri, 'voxel', 'mni');
    % Get cut plane in MRI coordinates
    cutPlane = MNIplane * vox2mni;
    % Get voxel indices under the MNI plane defined in input
    mriSize = size(sMri.Cube);
    [i,j,k] = ndgrid(1:mriSize(1), 1:mriSize(2), 1:mriSize(3));
    maskCut = (cutPlane(1)*i + cutPlane(2)*j + cutPlane(3)*k + cutPlane(4) < 0);
    % Checking for errors
    strPlane = sprintf('%1.3fx + %1.3fy + %1.3fz + %1.3f = 0', MNIplane);
    if (nnz(maskCut) == 0)
        errMsg = ['No voxels are located below the plane: ', strPlane];
        return;
    elseif (nnz(maskCut) == numel(maskCut))
        errMsg = ['All the voxels are located below the plane: ', strPlane];
        return;
    end
    % Replicate over multiple volumes
    if (length(mriSize) == 4) && (mriSize(4) > 1)
        maskCut = repmat(maskCut, [1 1 1 mriSize(4)]);
    end
    % Set to zero the voxels below the plane
    sMri.Cube(maskCut) = 0;
end

%% ===== INSTALL MRI_DEFACE =====
function [exePath, talFile, faceFile, errMsg] = InstallMriDeface(isInteractive)
    exePath = [];
    errMsg = [];
    talFile = [];
    faceFile = [];
    curdir = pwd;
    % Test if FreeSurfer is installed
    fsHome = getenv('FREESURFER_HOME');
    if ~isempty(fsHome)
        exePath = bst_fullfile(fsHome, 'bin', 'mri_deface');
        talFile = bst_fullfile(fsHome, 'average', 'talairach_mixed_with_skull.gca');
        faceFile = bst_fullfile(fsHome, 'average', 'face.gca');
        if file_exist(exePath) && file_exist(talFile) && file_exist(faceFile)
            disp(['BST> mri_deface found in FreeSurfer folder: ' fsHome]);
            return;
        end
        exePath = [];
    end
    % Executable mri_deface not found: Trying to install locally
    disp('BST> Deface: Variable FREESURFER_HOME not set. Trying to install mri_deface locally...');
    % Get download url
    osType = bst_get('OsType', 0);
    switch(osType)
        case 'linux32',  url = 'https://surfer.nmr.mgh.harvard.edu/pub/dist/mri_deface/mri_deface-v1.22-Linux.gz';
        case 'linux64',  url = 'https://surfer.nmr.mgh.harvard.edu/pub/dist/mri_deface/mri_deface-v1.22-Linux64.gz';
        case 'mac32',    url = 'https://surfer.nmr.mgh.harvard.edu/pub/dist/mri_deface/mri_deface-v1.22-MacOS-Leopard-intel.gz';
        case 'mac64',    url = 'https://surfer.nmr.mgh.harvard.edu/pub/dist/mri_deface/mri_deface-v1.22-MacOS-Leopard-intel.gz';
        otherwise,       errMsg = 'The program mri_deface is not available for your operating system.'; return;
    end
    talUrl = 'https://surfer.nmr.mgh.harvard.edu/pub/dist/mri_deface/talairach_mixed_with_skull.gca.gz';
    faceUrl = 'https://surfer.nmr.mgh.harvard.edu/pub/dist/mri_deface/face.gca.gz';

    % Local folder where to install mri_deface
    mriDefaceDir = bst_fullfile(bst_get('BrainstormUserDir'), 'mri_deface', osType);
    % If dir doesn't exist in user folder, try to look for it in the Brainstorm folder
    if ~isdir(mriDefaceDir)
        mriDefaceDirMaster = bst_fullfile(bst_get('BrainstormHomeDir'), 'mri_deface', osType);
        if isdir(mriDefaceDirMaster)
            mriDefaceDir = mriDefaceDirMaster;
        end
    end
    % URL file defines the current version
    urlFile = bst_fullfile(mriDefaceDir, 'url');

    % Read the previous download url information
    if isdir(mriDefaceDir) && file_exist(urlFile)
        fid = fopen(urlFile, 'r');
        prevUrl = fread(fid, [1 Inf], '*char');
        fclose(fid);
    else
        prevUrl = '';
    end

    % Local installation
    exePath = bst_fullfile(mriDefaceDir, 'mri_deface');
    talFile = bst_fullfile(mriDefaceDir, 'talairach_mixed_with_skull.gca');
    faceFile = bst_fullfile(mriDefaceDir, 'face.gca');
    % If binary file doesnt exist: download
    if ~isdir(mriDefaceDir) || ~file_exist(exePath) || ~strcmpi(prevUrl, url)
        % If folder exists: delete
        if isdir(mriDefaceDir)
            file_delete(mriDefaceDir, 1, 3);
        end
        % Create folder
        res = mkdir(mriDefaceDir);
        if ~res
            errMsg = ['Error: Cannot create folder' 10 mriDefaceDir];
            return
        end
        % Message
        if isInteractive
            isOk = java_dialog('confirm', ...
                ['FreeSurfer or mri_deface are not installed on your computer (or out-of-date).' 10 10 ...
                 'Download and the latest version of mri_deface?'], 'mri_deface');
            if ~isOk
                errMsg = 'Download aborted by user';
                return;
            end
        end
        % Download file
        errMsg1 = gui_brainstorm('DownloadFile', url, [exePath '.gz'], 'Download mri_deface');
        errMsg2 = gui_brainstorm('DownloadFile', talUrl, [talFile '.gz'], 'mri_deface template (1)');
        errMsg3 = gui_brainstorm('DownloadFile', faceUrl, [faceFile '.gz'], 'mri_deface template (2)');
        % If file was not downloaded correctly
        if ~isempty(errMsg1) || ~isempty(errMsg2) || ~isempty(errMsg3)
            errMsg = ['Impossible to download mri_deface:' 10 errMsg1];
            return;
        end
        % Display again progress bar
        bst_progress('text', 'Installing mri_deface...');
        % Unzip file
        cd(mriDefaceDir);
        system(['gunzip ' exePath '.gz']);
        system(['gunzip ' talFile '.gz']);
        system(['gunzip ' faceFile '.gz']);
        system(['chmod a+x ' exePath]);
        cd(curdir);
        % Save download URL in folder
        fid = fopen(urlFile, 'w');
        fwrite(fid, url);
        fclose(fid);
    end
    % If the executable is still not accessible
    if ~file_exist(exePath)
        errMsg = ['mri_convert could not be installed in: ' mriDefaceDir];
    end
end

author: brainstorm-tools


text:
classdef TestLSDDetector
    %TestLSDDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','tsukuba_l.png');
        kfields = {'angle', 'class_id', 'octave', 'pt', 'response', ...
            'size', 'startPoint', 'endPoint', 'startPointInOctave', ...
            'endPointInOctave', 'lineLength', 'numOfPixels'};
    end

    methods (Static)
        function test_detect_img
            obj = cv.LSDDetector();
            img = imread(TestLSDDetector.im);

            klines = obj.detect(img, 'Scale',2, 'NumOctaves',1);
            validateattributes(klines, {'struct'}, {'vector'});
            assert(all(ismember(TestLSDDetector.kfields, fieldnames(klines))));
        end

        function test_detect_imgset
            img = imread(TestLSDDetector.im);
            imgs = {img, img};
            obj = cv.LSDDetector();

            klines = obj.detect(imgs, 'Scale',2, 'NumOctaves',1);
            validateattributes(klines, {'cell'}, ...
                {'vector', 'numel',numel(imgs)});
            cellfun(@(kpt) validateattributes(kpt, {'struct'}, ...
                {'vector'}), klines);
            cellfun(@(kpt) assert(all(ismember(TestLSDDetector.kfields, ...
                fieldnames(kpt)))), klines);
        end

        function test_detect_mask
            img = imread(TestLSDDetector.im);
            mask = zeros(size(img), 'uint8');
            mask(:,1:end/2) = 255;  % only search left half of the image
            obj = cv.LSDDetector();

            klines = obj.detect(img, 'Mask',mask);
            validateattributes(klines, {'struct'}, {'vector'});
            assert(all(ismember(TestLSDDetector.kfields, fieldnames(klines))));
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function subpak_test04 ( )

%*****************************************************************************80
%
%% TEST04 tests AXIS_LIMITS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  xmin = 67.3;
  xmax = 114.7;
  ndivs = 6;

  [ pxmin, pxmax, pxdiv, nticks ] = axis_limits ( xmin, xmax, ndivs );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  AXIS_LIMITS adjusts plot limits to "nicer" values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Input XMIN =    %f\n', xmin );
  fprintf ( 1, '  Input XMAX =    %f\n', xmax );
  fprintf ( 1, '  Input NDIVS =   %d\n', ndivs );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Output PXMIN =  %f\n', pxmin );
  fprintf ( 1, '  Output PXMAX =  %f\n', pxmax );
  fprintf ( 1, '  Output PXDIV =  %f\n', pxdiv );
  fprintf ( 1, '  Output NTICKS = %d\n', nticks );

  xmin = -26.0;
  xmax = +26.0;
  ndivs = 10;

  [ pxmin, pxmax, pxdiv, nticks ] = axis_limits ( xmin, xmax, ndivs );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Input XMIN =    %f\n', xmin );
  fprintf ( 1, '  Input XMAX =    %f\n', xmax );
  fprintf ( 1, '  Input NDIVS =   %d\n', ndivs );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Output PXMIN =  %f\n', pxmin );
  fprintf ( 1, '  Output PXMAX =  %f\n', pxmax );
  fprintf ( 1, '  Output PXDIV =  %f\n', pxdiv );
  fprintf ( 1, '  Output NTICKS = %d\n', nticks );

  return
end


text:
function net = netunpak(net, w)
%NETUNPAK Separates weights vector into weight and bias matrices.
%
%	Description
%	NET = NETUNPAK(NET, W) takes an net network data structure NET and  a
%	weight vector W, and returns a network data structure identical to
%	the input network, except that the componenet weight matrices have
%	all been set to the corresponding elements of W.  If there is  a MASK
%	field in the NET data structure, then the weights in W are placed in
%	locations corresponding to non-zero entries in the mask (so W should
%	have the same length as the number of non-zero entries in the MASK).
%
%	See also
%	NETPAK, NETFWD, NETERR, NETGRAD
%

%	Copyright (c) Ian T Nabney (1996-2001)

unpakstr = [net.type, 'unpak'];

% Check if we are being passed a masked set of weights
if (isfield(net, 'mask'))
   if length(w) ~= size(find(net.mask), 1)
      error('Weight vector length does not match mask length')
   end
   % Do a full pack of all current network weights
   pakstr = [net.type, 'pak'];
   fullw = feval(pakstr, net);
   % Replace current weights with new ones
   fullw(logical(net.mask)) = w;
   w = fullw;
end

net = feval(unpakstr, net, w);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [strength, time, x, y, storedMacs] = parseWirelessData(fileName)

% PARSEWIRELESSDATA Load wireless strength data.
% FORMAT
% DESC loads wireless strength data as written by the University of
% Washington by Brian Ferris.
% ARG fileName : the file name to load the data from.
% RETURN strength : the strengths of the wireless access points
% read.
% RETURN time : the time at which the wireless access points were
% read.
% RETURN x : the 'ground truth' x location of the device at the
% time.
% RETURN y : the 'ground truth' y location of the device at the
% time.
% RETURN macs : the mac addresses of the devices that were read.
%
% SEEALSO : lvmLoadData
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
%
% MODIFICATIONS : Neil D. Lawrence, 2009

% DATASETS

fid = fopen(fileName);
if fid == -1
  error(['No such file name ' fileName])
end
readLine = fgets(fid);
storedMacs = {};
counter = 0;
data = [];
curTime = -1;
timeIndex = 1;
while readLine ~= -1
  counter = counter + 1;
  parts = stringSplit(readLine, ' ');
  index = -1;
  for i = 1:length(storedMacs)
    if strcmp(parts{2}, storedMacs{i})
      index = i;
    end
  end
  if index == -1
    index = length(storedMacs) + 1;
    storedMacs{index} = parts{2};
  end
  pointTime = str2num(parts{1});
  if curTime == -1
    curTime = pointTime;
  end
  if pointTime > curTime
    timeIndex = timeIndex + 1;
    curTime = pointTime;
  end
  data(counter, :) = [timeIndex index str2num(parts{1})  str2num(parts{3}) str2num(parts{4}) str2num(parts{5})];
  readLine = fgets(fid);
end
start = min(data(:, 3));
data(:, 3) = data(:, 3) - start + 1;
indMax = max(data(:, 2));
points = max(data(:, 1));

x = zeros(points, 1);
y = zeros(points, 1);
time = zeros(points, 1);
strength = repmat(-92, points, indMax);

for i = 1:size(data, 1);
  x(data(i, 1), 1) = data(i, 4);
  y(data(i, 1), 1) = data(i, 5);
  time(data(i, 1), 1) = data(i, 3);
  strength(data(i, 1), data(i, 2)) = data(i, 6);
end
[storedMacs, ind] = sort(storedMacs);
strength = strength(:, ind);


text:
% em_test_setup.m
%
% create sample image, system matrix, and sinograms for examples
% and testing of Poisson emission maximum likelihood (ML) algorithms
% creates: ig sg xtrue G proj ci ytrue ri yi
%
% Copyright Jan 1998, Jeff Fessler, University of Michigan

% true emission image
if ~isvar('xtrue'), printm 'xtrue'
	if ~isvar('ig')
		ig = image_geom('nx', 64, 'ny', 60, 'fov', 500);
	end
	xtrue = read_zubal_emis('nx', ig.nx, 'ny', ig.ny);
	mumap = read_zubal_attn('nx', ig.nx, 'ny', ig.ny);
	im plc 2 3
	im(1, xtrue, 'emission image'), cbar
	im(2, mumap, 'attenuation map'), cbar

	% reconstruction mask (which pixels do we estimate?)
	ig.mask = ig.circ(220, 180) > 0;
	im(3, ig.mask + xtrue, 'support mask + xtrue')
end

% system matrix G
if ~isvar('G'), printm 'system'
	sg = sino_geom('par', 'nb', ig.nx+2, 'na', ig.ny*3/2, ...
		'dr', 528 / (ig.nx+2));

	% simple strip-integral system model
	G = Gtomo2_strip(sg, ig, 'single', 1);

	if isvar('f.wtf') && has_mex_jf
		if exist(f.wtf), delete(f.wtf), end
		wtf_write(f.wtf, G, ig.nx, ig.ny, sg.nb, sg.na);
	end
	if isvar('f.wtr') && has_mex_jf && has_aspire
		if exist(f.wtr), delete(f.wtr), end
		os_run(sprintf('wt -chat 0 col2row %s %s', f.wtr, f.wtf))
	end
end

% noisy measurements
if ~isvar('yi'), printm 'data yi'
	proj = G * xtrue;
	li = G * mumap;
	printm('Maximum line integral = %g', max(li(:)))
	if ~isvar('f.count'), f.count = 1e5; end
	% detector efficiency variations per CTI 931 PET scanner
	ci = exp(0.3 * randn(size(proj)));
	ci = ci .* exp(-li);
	ci = f.count / sum(ci(:) .* proj(:)) * ci;
	ci = dsingle(ci);
	ytrue = ci .* proj;
	if ~isvar('f.randpercent')
		f.randpercent = 10;
	end
	ri = f.randpercent / 100 * mean(ytrue(:)) * sg.ones;
	ri = dsingle(ri);
	rng(0)
	yi = poisson(ytrue + ri);

	im(4, ytrue, 'ytrue: true projections'), cbar
	im(5, yi, 'yi: noisy projections'), cbar
	clear ytrue proj
end

% FBP reconstruction
if ~isvar('xfbp'), printm 'fbp'
	xfbp = em_fbp(sg, ig, yi, ci, ri);
	xfbp = max(xfbp, 0);
	im(6, xfbp, 'FBP Reconstruction'), cbar
end

% save to files if needed
if isvar('f.yi')
	fld_write(f.yi, yi, 'check', 0)
end
if isvar('f.ci')
	fld_write(f.ci, ci, 'check', 0)
end
if isvar('f.ri')
	fld_write(f.ri, ri, 'check', 0)
end
if isvar('f.mask')
	fld_write(f.mask, ig.mask, 'check', 0)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%% The Class embedding
%
%%
%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function calc_dist_matrix_shape_collection(path_shapes,path_distance_matrix,num_workers)
    d = dir([path_shapes,'*.mat']);

    if ~exist(path_distance_matrix, 'dir')
       mkdir(path_distance_matrix)
    end

    myCluster = parcluster;
    myCluster.NumWorkers = num_workers;
    parpool(num_workers);

    parfor i=1:numel(d)
        try
            S=load([path_shapes,d(i).name]);
            D = calc_dist_matrix(S); D = single(D);
            parsave([path_distance_matrix,d(i).name],D);
            display(i)
        catch
            display(d(i).name)
        end
    end

    delete(gcp('nocreate'))
end

author: OshriHalimi


text:
function [ftlb] = kJ2ftlb(kJ)
% Convert energy or work from kilojoules to foot-pounds.
% Chad A. Greene 2012
ftlb = kJ*737.56217557;

author: Sable

repo: mcbench-benchmarks


text:
classdef prtBrvMixtureMcmcTraining

    properties

         randnState = randn('seed');
         randState = rand('seed');
         startTime = now;
         endTime = [];

         membershipModels
         componentMemberships
         componentLogLikelihoods
         componentLogLikelihoodsInModel

         iterations = struct('logLikelihood',[]);
    end
    methods
        function self = prtBrvMixtureVbTraining(varargin)
            self = prtBrvMixtureMcmcTraining(self,varargin{:});
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function b = iseven(x)
%ISEVEN True for even numbers.
%
%   ISEVEN(X) returns 1's where the elements of X are even numbers and 0's
%   where they are not.
%
%   See also ISINT, ISODD.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-11-03 18:55:21 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 1, nargin));
   if ~isnumeric(x)
      error('Argument must be a numeric array.');
   end

   cls = class(x);                      % class of input argument
   if isempty(x)
      b = feval(cls, x);                % return empty array of same class
   else
      switch cls
         case 'double'
            b = ~mod(x, 2);
         case 'single'
            % "mod" is not defined for class "single"; so convert input to double, compare,
            % and convert back
            b = single(~mod(double(x), 2));
         case {'uint8', 'uint16', 'uint32', 'uint64'}
            b = ~bitand(x, 1);
         case {'int8', 'int16', 'int32', 'int64'}
            error('Not implemented for classes int8, int16, int32, and int64.');
         otherwise
            error('Argument is of unrecognized class.');
      end
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [pair_label] = generate_pair_label(label)
% this file converts the labels of a set of scale labels (from the smallest candidate to the largest one) into pair-wise
% labels.

num_scale = numel(label);
base_ind = ceil(num_scale/2);

pair_label = [2:base_ind, base_ind:num_scale-1;...
               1:base_ind-1, base_ind+1:num_scale];

end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function blas1_z_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests DZNRM2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  x = [ ...
     2.0 - 1.0 * i, ...
    -4.0 - 2.0 * i, ...
     3.0 + 1.0 * i, ...
     2.0 + 2.0 * i, ...
    -1.0 - 1.0 * i ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  DZNRM2 returns the Euclidean norm of a complex vector.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The vector X:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f\n', j, real ( x(j) ), imag ( x(j) ) );
  end

  incx = 1;
  norm = dznrm2 ( n, x, incx );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The L2 norm of X is %f\n', norm );

  return
end


text:
function startup()
 fprintf('Entering startup\n');
  curdir = fileparts(mfilename('fullpath'));
    if ~(isdeployed)
    fprintf('Adding paths\n');
      addpath(genpath(fullfile(curdir, 'utils')));
      addpath(genpath(fullfile(curdir, 'functions')));
      addpath(genpath(fullfile(curdir, 'bin')));
      addpath(genpath(fullfile(curdir, 'experiments')));
      addpath(genpath(fullfile(curdir, 'imdb')));
      addpath(genpath(fullfile(curdir, 'datasets/coco2014/')));
     fprintf('Starting caffe\n');
      caffe_path = fullfile(curdir, 'external', 'caffe', 'matlab');
      if exist(caffe_path, 'dir') == 0
        error('matcaffe is missing from external/caffe/matlab; See README.md');
      end
      addpath(genpath(caffe_path));
    fprintf('Done\n');
    end
    fprintf('mask_rcnn startup done\n');
end

author: jasjeetIM

repo: Mask-RCNN

save_path: github-repos/MATLAB/jasjeetIM-Mask-RCNN


text:
%Computes the position difference in meters between two coordinates defined
%with Llh.
%Inputs:
%ref=[Lat, lon, height]' of reference.
%height=[h1;h2;...;hn] where n is the number of coordinate
%hpos must be oneof the following: (each position must be on a new row)
%hpos=3rd row of Cen=[-cos(Lat)cos(lon) -cos(Lat)sin(lon) -sin(Lat)];
%or
%hpos=[Lat lon]

%Output: Metric difference between corrdinates defined in the nav frame of
%reference.

function metricdif=posdiff(hpos, height, ref)
E_SQR=0.00669437999014;
SM_AXIS=6378137;

%Convert ref into ecef and compute Cen(ref)
sL=sin(ref(1));
cL=cos(ref(1));
sl=sin(ref(2));
cl=cos(ref(2));
h=ref(3);

Re=SM_AXIS/(sqrt(1.0-E_SQR*sL*sL));
refECEF=[(Re+h)*cL*cl (Re+h)*cL*sl (Re*(1-E_SQR)+h)*sL];
CenRef=[-sL*cl -sL*sl cL;-sl cl 0;-cL*cl -cL*sl -sL];

%%%%Convert hpos to ecef
%Be sure that hpos and height is in correct order
if (size(height,2)>size(height,1))
    height=height';
    hpos=hpos';
end

%Compute the normal to the elipsoid in ECEF frame
if size(hpos,2)==2 %hpos is defined as lat/lon. Convert it to 3rd row of Cen
    normal=[-cos(hpos(:,1)).*cos(hpos(:,2)) -cos(hpos(:,1)).*sin(hpos(:,2)) -sin(hpos(:,1))];
elseif (size(hpos,2)==3) %hpos is defined
    normal=hpos;
end
Re=SM_AXIS./(1.0-E_SQR*normal(:,3).^2).^0.5;

%Compute ecef coordinates
posECEF=-[(Re+height).*normal(:,1) (Re+height).*normal(:,2) (Re*(1-E_SQR)+height).*normal(:,3)];

%Difference
metricdif=(posECEF-ones(size(posECEF,1),1)*refECEF)*CenRef';

repo: ekfmonoslam


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtDomRegularize2D.m                          |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 25.11.2018                                    |
%|  / 0 \  |   LAST MODIF : 25.11.2018                                    |
%| ( === ) |   SYNOPSIS   : Singular kernel regularization (in debug mode)|
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N   = 1e2;
gss = 3;

% Spherical mesh
circle = mshCircle(N,1);

% Square mesh
square = mshSquare(2*N,[3 3]);

% Graphical representation
figure
plot(circle)
hold on
plot(square,'w')
axis equal
alpha(0.5)

% Domain
sigma = dom(circle,gss);

% Finite elements
u = fem(circle,'P0');
v = fem(circle,'P1');

% Gren kernel
Gxy      = @(X,Y) femGreenKernel(X,Y,'[log(r)]',[]);
dyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[log(r)]1',[]);
dyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[log(r)]2',[]);
dyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[log(r)]3',[]);

% Single radiation P0
M  = integral(square.vtx,sigma,Gxy,u);
Mr = regularize(square.vtx,sigma,'[log(r)]',u);
norm(M+Mr,'inf')

% Double radiation P0
M  = integral(square.vtx,sigma,dyGxy,ntimes(u));
Mr = regularize(square.vtx,sigma,'grady[log(r)]',ntimes(u));
norm(M+Mr,'inf')

% Single radiation P1
M  = integral(square.vtx,sigma,Gxy,v);
Mr = regularize(square.vtx,sigma,'[log(r)]',v);
norm(M+Mr,'inf')

% Single layer P0
M  = integral(sigma,sigma,u,Gxy,u);
Mr = regularize(sigma,sigma,u,'[log(r)]',u);
norm(M+Mr,'inf')

% Double layer P0
M  = integral(sigma,sigma,u,dyGxy,ntimes(u));
Mr = regularize(sigma,sigma,u,'grady[log(r)]',ntimes(u));
norm(M+Mr,'inf')

% Single layer P1
M  = integral(sigma,sigma,u,Gxy,v);
Mr = regularize(sigma,sigma,u,'[log(r)]',v);
norm(M+Mr,'inf')

% Hypersingular P1
M  = integral(sigma,sigma,ntimes(u),Gxy,ntimes(v));
Mr = regularize(sigma,sigma,ntimes(u),'[log(r)]',ntimes(v));
norm(M+Mr,'inf')

M  = integral(sigma,sigma,nxgrad(u),Gxy,nxgrad(v));
Mr = regularize(sigma,sigma,nxgrad(u),'[log(r)]',nxgrad(v));
norm(M+Mr,'inf')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%  Make ANALYZE 7.5 data structure specified by a 3D or 4D matrix.
%  Optional parameters can also be included, such as: voxel_size,
%  origin, datatype, and description.
%
%  Once the ANALYZE structure is made, it can be saved into ANALYZE 7.5
%  format data file using "save_untouch_nii" command (for more detail,
%  type: help save_untouch_nii).
%
%  Usage: ana = make_ana(img, [voxel_size], [origin], [datatype], [description])
%
%  Where:
%
%	img:		a 3D matrix [x y z], or a 4D matrix with time
%			series [x y z t]. When image is in RGB format,
%			make sure that the size of 4th dimension is
%			always 3 (i.e. [R G B]). In that case, make
%			sure that you must specify RGB datatype to 128.
%
%	voxel_size (optional):	Voxel size in millimeter for each
%				dimension. Default is [1 1 1].
%
%	origin (optional):	The AC origin. Default is [0 0 0].
%
%	datatype (optional):	Storage data type:
%		2 - uint8,  4 - int16,  8 - int32,  16 - float32,
%		64 - float64,  128 - RGB24
%			Default will use the data type of 'img' matrix
%			For RGB image, you must specify it to 128.
%
%	description (optional):	Description of data. Default is ''.
%
%  e.g.:
%     origin = [33 44 13]; datatype = 64;
%     ana = make_ana(img, [], origin, datatype);    % default voxel_size
%
%  ANALYZE 7.5 format: http://www.rotman-baycrest.on.ca/~jimmy/ANALYZE75.pdf
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function ana = make_ana(varargin)

   ana.img = varargin{1};
   dims = size(ana.img);
   dims = [4 dims ones(1,8)];
   dims = dims(1:8);

   voxel_size = [0 ones(1,3) zeros(1,4)];
   origin = zeros(1,5);
   descrip = '';

   switch class(ana.img)
      case 'uint8'
         datatype = 2;
      case 'int16'
         datatype = 4;
      case 'int32'
         datatype = 8;
      case 'single'
         datatype = 16;
      case 'double'
         datatype = 64;
      otherwise
         error('Datatype is not supported by make_ana.');
   end

   if nargin > 1 & ~isempty(varargin{2})
      voxel_size(2:4) = double(varargin{2});
   end

   if nargin > 2 & ~isempty(varargin{3})
      origin(1:3) = double(varargin{3});
   end

   if nargin > 3 & ~isempty(varargin{4})
      datatype = double(varargin{4});

      if datatype == 128 | datatype == 511
         dims(5) = [];
         dims = [dims 1];
      end
   end

   if nargin > 4 & ~isempty(varargin{5})
      descrip = varargin{5};
   end

   if ndims(ana.img) > 4
      error('NIfTI only allows a maximum of 4 Dimension matrix.');
   end

   maxval = round(double(max(ana.img(:))));
   minval = round(double(min(ana.img(:))));

   ana.hdr = make_header(dims, voxel_size, origin, datatype, ...
	descrip, maxval, minval);
   ana.filetype = 0;
   ana.ext = [];
   ana.untouch = 1;

   switch ana.hdr.dime.datatype
   case 2
      ana.img = uint8(ana.img);
   case 4
      ana.img = int16(ana.img);
   case 8
      ana.img = int32(ana.img);
   case 16
      ana.img = single(ana.img);
   case 64
      ana.img = double(ana.img);
   case 128
      ana.img = uint8(ana.img);
   otherwise
      error('Datatype is not supported by make_ana.');
   end

   return;					% make_ana

%---------------------------------------------------------------------
function hdr = make_header(dims, voxel_size, origin, datatype, ...
	descrip, maxval, minval)

   hdr.hk   = header_key;
   hdr.dime = image_dimension(dims, voxel_size, datatype, maxval, minval);
   hdr.hist = data_history(origin, descrip);

   return;					% make_header

%---------------------------------------------------------------------
function hk = header_key

    hk.sizeof_hdr       = 348;			% must be 348!
    hk.data_type        = '';
    hk.db_name          = '';
    hk.extents          = 0;
    hk.session_error    = 0;
    hk.regular          = 'r';
    hk.hkey_un0         = '0';

    return;					% header_key

%---------------------------------------------------------------------
function dime = image_dimension(dims, voxel_size, datatype, maxval, minval)

   dime.dim = dims;
   dime.vox_units = 'mm';
   dime.cal_units = '';
   dime.unused1 = 0;
   dime.datatype = datatype;

   switch dime.datatype
   case   2,
      dime.bitpix = 8;  precision = 'uint8';
   case   4,
      dime.bitpix = 16; precision = 'int16';
   case   8,
      dime.bitpix = 32; precision = 'int32';
   case  16,
      dime.bitpix = 32; precision = 'float32';
   case  64,
      dime.bitpix = 64; precision = 'float64';
   case 128
      dime.bitpix = 24;  precision = 'uint8';
   otherwise
      error('Datatype is not supported by make_ana.');
   end

   dime.dim_un0 = 0;
   dime.pixdim = voxel_size;
   dime.vox_offset = 0;
   dime.roi_scale = 1;
   dime.funused1 = 0;
   dime.funused2 = 0;
   dime.cal_max = 0;
   dime.cal_min = 0;
   dime.compressed = 0;
   dime.verified = 0;
   dime.glmax = maxval;
   dime.glmin = minval;

   return;					% image_dimension

%---------------------------------------------------------------------
function hist = data_history(origin, descrip)

   hist.descrip = descrip;
   hist.aux_file = 'none';
   hist.orient = 0;
   hist.originator = origin;
   hist.generated = '';
   hist.scannum = '';
   hist.patient_id = '';
   hist.exp_date = '';
   hist.exp_time = '';
   hist.hist_un0 = '';
   hist.views = 0;
   hist.vols_added = 0;
   hist.start_field = 0;
   hist.field_skip = 0;
   hist.omax = 0;
   hist.omin = 0;
   hist.smax = 0;
   hist.smin = 0;

   return;					% data_history

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA