mfbox_kernelica.m

toolbox for spm 5 for data, model free analysis

function [comp,psi]=mfbox_kernelica(source,ics,method)
% Copyright by Peter Gruber 
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

error(nargchk(1,3,nargin))     % check no. of input args

dim = size(source,1);

if (nargin<2), ics = dim; end
if (nargin<3), method = 'geneticgausspearson'; end

method = lower(method);

switch method
    case 'geneticgausspearson'
        borderBase = [2.6,4];
        borderSlope = [0,1];
        maxNumIterations = 100;
        epsilon = 0.001;
        [optcols,target,kernelmat] = basisopt(source,10*ics,ics,100*ics,'gauss');
        optcols = optcols(size(optcols,1),:);
        psi = projection(source,optcols,kernelmat,'gauss');
        comp = mfbox_fasticapearson(psi',epsilon,maxNumIterations,borderBase, ...
            borderSlope,ics,1,dim,'off');
    case 'geneticpolypearson'
        borderBase = [2.6,4];
        borderSlope = [0,1];
        maxNumIterations = 100;
        epsilon = 0.001;
        [optcols,target,kernelmat] = basisopt(source,10*ics,ics,100*ics,'poly',1,5);
        optcols = optcols(size(optcols,1),:);
        psi = projection(source,optcols,kernelmat,'poly',1,5);
        comp = mfbox_fasticapearson(psi',epsilon,maxNumIterations,borderBase, ...
            borderSlope,ics,1,dim,'off');
    case 'geneticgaussjade'
        [optcols,target,kernelmat] = basisopt(source,10*ics,ics,100*ics,'gauss');
        optcols = optcols(size(optcols,1),:);
        psi = projection(source,optcols,kernelmat,'gauss');
        [A,comp] = mfbox_jade(psi',ics);
    otherwise
        [A,comp] = mfbox_jade(source',ics);
end

return


function [optcols,target,kernel_opt]=basisopt(X,ind_num,ind_size,maxgen,kernel,c,d,verbose)

if (nargin<6), c = 1; end
if (nargin<7), d = 1; end
if (nargin<8), verbose = 0; end

ggap = 0.9;

if (verbose>0), warning('Number of individuals %d, Generation Gap: %f',ind_num,ggap); end

[rows_X,cols_X] = size(X);

chrom = initpop(ind_num, ind_size, cols_X);
ObjV = condobj(X,chrom,kernel,c,d);

if (verbose>0), warning('Minimale Konditionszahl des Ausgangchromosomensatzes: %f',min(ObjV)); end

gen = 0;

target = zeros(1,maxgen);
optcols = zeros(maxgen,size(chrom,2));
while (gen<maxgen)
    fitnv = ranking(ObjV);
    SelCh = select('sus',chrom,fitnv,ggap);
    SelCh = sort(SelCh,2);
    SelCh = recombin('xovsp',SelCh,0.7);
    SelCh = mut(SelCh,0.1,(cols_X+1)*ones(1,ind_size));
    [rows_SelCh,cols_SelCh] = size(SelCh);
    SelCh = SelCh+(SelCh==zeros(rows_SelCh,cols_SelCh).*round(cols_X*rand(rows_SelCh,cols_SelCh)));
    ObjVSel = condobj(X,SelCh,kernel,c,d);
    [chrom,ObjV] = reins(chrom,SelCh,1,1,ObjV,ObjVSel);
    gen = gen+1;
    [target(gen),index] = min(ObjV);
    optcols(gen,:) = chrom(index,:);

    if (verbose>0), warning('current generation %d',gen); end
end

[ObjV,kernel_opt] = condobj(X,optcols(gen,:),kernel,c,d);
if (verbose>1), plot(target); end


function psi=projection(X,basis,kernelm,kernel,c,d)

if (nargin<5), c = 1; end
if (nargin<6), d = 1; end

[rows_X,cols_X] = size(X);
[rows_basis,dim_basis] = size(basis);

kbasis = X(:,basis);

if (strcmp(kernel,'poly'))
    Xvdot = (kbasis'*X+c).^d;

    psi = kernelm^(-1/2)*Xvdot;
end

if (strcmp(kernel,'gauss'))
    Xvdot = [];
    for k = 1:dim_basis
        D = kbasis(:,k)*ones(1,cols_X);
        B = D-X;
        B = B.*B;
        diflength = sum(B,1);
        Xvdot = cat(1,Xvdot,-diflength/c);
    end
    Xvdot = exp(Xvdot);

    psi = kernelm^(-1/2)*Xvdot;
end


function chrom=initpop(ind_num,ind_size,max_value)

chrom = zeros(ind_num,ind_size);
for k = 1:ind_num
    permutation = (randperm(max_value));
    chrom(k,:) = permutation(1:ind_size);
end

chrom = sort(chrom,2);


function [Condnumbers,kernel_opt]=condobj(X,chrom,kernel,c,d)

[rows_chrom,cols_chrom] = size(chrom);

Condnumbers = zeros(rows_chrom,1);
if (strcmp(kernel,'poly'))
    cond_min = inf;
    for k = 1:rows_chrom
        sel_matrix = X(:,chrom(k,:));
        kernel_matrix = ((sel_matrix'*sel_matrix)+c).^d;
        Condnumbers(k) = cond(kernel_matrix);
        if (Condnumbers(k)<cond_min)
            cond_min = Condnumbers(k);
            kernel_opt = kernel_matrix;
        end
    end
end

if (strcmp(kernel,'gauss'))
    cond_min = inf;
    for k = 1:rows_chrom
        sel_matrix = X(:,chrom(k,:));
        EXPONENT = [];
        for l = 1:cols_chrom;
            d = sel_matrix(:,l);
            D = d*ones(1,cols_chrom);
            B = sel_matrix-D;
            B = B.*B;
            diflength = (sum(B,1))';
            EXPONENT = cat(2,EXPONENT,-diflength/c);
        end
        kernel_matrix = exp(EXPONENT);
        Condnumbers(k) = cond(kernel_matrix);
        if (Condnumbers(k)<cond_min)
            cond_min = Condnumbers(k,1);
            kernel_opt = kernel_matrix;
        end
    end
end


function [chrom,objvch]=reins(chrom,SelCh,subpop,InsOpt,objvch,ObjVSel)

NindP = size(chrom,1);
NindO = size(SelCh,1);

if (nargin == 2), subpop = 1; end
if (nargin > 2),
    if isempty(subpop), subpop = 1;
    elseif isnan(subpop), subpop = 1;
    end
end

Nind = NindP/subpop;
Nsel = NindO/subpop;

Isobjvch = 0; IsObjVSel = 0;
if (nargin>4),
    Isobjvch = 1;
end

if (nargin<4), insr = 1.0; sel = 0; end
if (nargin>=4),
    if (isempty(InsOpt)), insr = 1.0; sel = 0;
    elseif (isnan(InsOpt)), insr = 1.0; sel = 0;
    else
        insr = NaN; sel = NaN;
        if (length(InsOpt) >= 1), sel = InsOpt(1); end
        if (length(InsOpt) >= 2), insr = InsOpt(2); end
        if (isnan(sel)), sel = 0; end
        if (isnan(insr)), insr =1.0; end
    end
end

if (insr==0), return; end
NIns = min(max(floor(insr*Nsel+.5),1),Nind);

for irun=1:subpop,
    if (sel==1)
        [dummy,ChIx] = sort(-objvch((irun-1)*Nind+1:irun*Nind));
    else
        [dummy,ChIx] = sort(rand(Nind,1));
    end
    PopIx = ChIx((1:NIns)')+ (irun-1)*Nind;
    if (NIns<Nsel)
        [dummy,OffIx] = sort(ObjVSel((irun-1)*Nsel+1:irun*Nsel));
    else
        OffIx = (1:NIns)';
    end
    SelIx = OffIx((1:NIns)')+(irun-1)*Nsel;
    chrom(PopIx,:) = SelCh(SelIx,:);
    if (Isobjvch == 1 && IsObjVSel == 1), objvch(PopIx) = ObjVSel(SelIx); end
end


function newchrom=mut(oldchrom,Pm,baseV)

[Nind,Lind] = size(oldchrom);

if (nargin<2), Pm = 0.7/Lind; end
if (isnan(Pm)), Pm = 0.7/Lind; end

if (nargin<3), baseV = crtbase(Lind);  end
if (isnan(baseV)), baseV = crtbase(Lind);  end
if (isempty(baseV)), baseV = crtbase(Lind);  end

baseM = baseV(ones(Nind,1),:) ;

newchrom = rem(oldchrom+(rand(Nind,Lind)<Pm).*ceil(rand(Nind,Lind).*(baseM-1)),baseM);


function newchrom=recombin(REC_F,chrom,RecOpt,subpop)

Nind = size(chrom,1);

if (nargin<4), subpop = 1; end
if (nargin>3),
    if isempty(subpop), subpop = 1;
    elseif isnan(subpop), subpop = 1;
    end
end

Nind = Nind/subpop;

if (nargin<3), RecOpt = 0.7; end
if (nargin>2),
    if (isempty(RecOpt)), RecOpt = 0.7;
    elseif (isnan(RecOpt)), RecOpt = 0.7;
    end
end

newchrom = cell(1,subpop);
for irun = 1:subpop,
    chromSub = chrom((irun-1)*Nind+1:irun*Nind,:);
    newchromSub = feval(REC_F,chromSub,RecOpt);
    newchrom{irun} = newchromSub;
end
newchrom = cat(1,newchrom{:});


function fitnv=ranking(ObjV,RFun,subpop)

Nind = size(ObjV,1);

if (nargin<2), RFun = []; end
if (nargin>1), if isnan(RFun), RFun = []; end, end
if (numel(RFun) == 2)
    if (RFun(2) == 1), NonLin = 1;
    elseif (RFun(2) == 0), NonLin = 0;
    end
    RFun = RFun(1);
    if (isnan(RFun)), RFun = 2; end
end

if (nargin<3), subpop = 1; end
if (nargin>2)
    if isempty(subpop), subpop = 1;
    elseif isnan(subpop), subpop = 1;
    end
end

Nind = Nind/subpop;

if (isempty(RFun))
    RFun = 2*(0:Nind-1)'/(Nind-1);
elseif (numel(RFun) == 1)
    if NonLin == 1,
        Root1 = roots([RFun(1)-Nind,(RFun(1)*ones(1,Nind-1))]);
        RFun = (abs(Root1(1)) * ones(Nind,1)) .^ (0:Nind-1)';
        RFun = RFun / sum(RFun) * Nind;
    else
        RFun = 2-RFun + 2*(RFun-1)*[0:Nind-1]'/(Nind-1);
    end
end;

fitnv = [];

for irun = 1:subpop,
    ObjVSub = ObjV((irun-1)*Nind+1:irun*Nind);
    NaNix = isnan(ObjVSub);
    Validix = find(~NaNix);
    [ans,ix] = sort(-ObjVSub(Validix));

    ix = [find(NaNix);Validix(ix)];
    Sorted = ObjVSub(ix);

    i = 1;
    fitnvSub = zeros(Nind,1);
    for j = [find(Sorted(1:Nind-1)~=Sorted(2:Nind));Nind]',
        fitnvSub(i:j) = sum(RFun(i:j)) * ones(j-i+1,1) / (j-i+1);
        i =j+1;
    end

    [ans,uix] = sort(ix);
    fitnvSub = fitnvSub(uix);

    fitnv = [fitnv;fitnvSub];
end


function SelCh=select(sel_F,chrom,fitnv,ggap,subpop)

NindCh = size(chrom,1);

if (nargin<5), subpop = 1; end
if (nargin>4),
    if isempty(subpop), subpop = 1;
    elseif isnan(subpop), subpop = 1;
    end
end

Nind = NindCh/subpop;

if (nargin<4), ggap = 1; end
if (nargin>3),
    if isempty(ggap), ggap = 1;
    elseif isnan(ggap), ggap = 1;
    end
end

NSel = max(floor(Nind*ggap+.5),2);

SelCh = zeros(subpop,size(chrom,2));
for irun = 1:subpop,
    fitnvSub = fitnv((irun-1)*Nind+1:irun*Nind);
    ChrIx = feval(sel_F,fitnvSub,NSel)+(irun-1)*Nind;
    SelCh(irun,:) = chrom(ChrIx,:);
end


function basevec=crtbase(Lind,base)

[ml,LenL] = size(Lind);
if (nargin<2)
    base = 2*ones(LenL,1);
end
mb = size(base,1);

if (~ (ml>1 || mb>1)) && (~ ((LenL>1 && LenB>1 && LenL~=LenB) || ...
        (LenL==1 && LenB>1 ))) && (LenB==1 && LenL>1)
    base = base*ones(LenL,1);
end

basevec = [];
for i = 1:LenL
    basevec = [basevec,base(i)*ones(Lind(i),1)'];
end


function newchrom=xovsp(oldchrom,xovr)

newchrom = xovmp(oldchrom,xovr,1,0);


function NewChrIx=sus(fitnv,Nsel)

[Nind,ans] = size(fitnv);

cumfit = cumsum(fitnv);
trials = cumfit(Nind)/Nsel*(rand+(0:Nsel-1)');
Mf = cumfit(:,ones(1,Nsel));
Mt = trials(:,ones(1,Nind))';
[NewChrIx,ans] = find(Mt<Mf & [zeros(1,Nsel);Mf(1:Nind-1,:)]<=Mt);

[ans,shuf] = sort(rand(Nsel, 1));
NewChrIx = NewChrIx(shuf);


function newchrom=xovmp(oldchrom,Px,Npt,Rs)

[Nind,Lind] = size(oldchrom);

if (Lind<2), newchrom = oldchrom; return; end

if (nargin<4), Rs = 0; end
if (nargin<3), Npt = 0; Rs = 0; end
if (nargin<2), Px = 0.7; Npt = 0; Rs = 0; end
if (isnan(Px)), Px = 0.7; end
if (isnan(Npt)), Npt = 0; end
if (isnan(Rs)), Rs = 0; end
if (isempty(Px)), Px = 0.7; end
if (isempty(Npt)), Npt = 0; end
if (isempty(Rs)), Rs = 0; end

Xops = floor(Nind/2);
DoCross = rand(Xops,1)<Px;
odd = 1:2:Nind-1;
even = 2:2:Nind;

Mask = ~Rs|(oldchrom(odd,:)~=oldchrom(even,:));
Mask = cumsum(Mask,2);

xsites(:,1) = Mask(:,Lind);
if (Npt>=2)
    xsites(:,1) = ceil(xsites(:,1).*rand(Xops,1));
end
xsites(:,2) = rem(xsites+ceil((Mask(:,Lind)-1).*rand(Xops,1)).*DoCross-1,Mask(:,Lind))+1;

Mask = (xsites(:,ones(1,Lind))<Mask) == ...
    (xsites(:,2*ones(1,Lind))<Mask);

if (~Npt)
    shuff = rand(Lind,Xops);
    [ans,shuff] = sort(shuff);
    for i = 1:Xops
        oldchrom(odd(i),:) = oldchrom(odd(i),shuff(:,i));
        oldchrom(even(i),:) = oldchrom(even(i),shuff(:,i));
    end
end

newchrom(odd,:) = (oldchrom(odd,:).* Mask)+(oldchrom(even,:).*(~Mask));
newchrom(even,:) = (oldchrom(odd,:).*(~Mask))+(oldchrom(even,:).*Mask);

if (rem(Nind,2))
    newchrom(Nind,:) = oldchrom(Nind,:);
end

if (~Npt)
    [ans,unshuff] = sort(shuff);
    for i = 1:Xops
        newchrom(odd(i),:) = newchrom(odd(i),unshuff(:,i));
        newchrom(even(i),:) = newchrom(even(i),unshuff(:,i));
    end
end