max/ourMELONS
1
0
Fork 0
Code Issues 1 Pull requests Projects Releases Packages Wiki Activity Actions
ourMELONS/matlab/linkage/linkageMix_fixK.m

1360 lines
48 KiB
Mathematica
Raw Normal View History

Added source Matlab code for reference
2019-12-16 16:47:21 +01:00
function [logml, npops, partitionSummary] = linkageMix_fixK(c,npops,nruns)
% Greedy search algorithm with fixed number of classes for linkage
% clustering.
global POP_LOGML; global PARTITION;
global CQ_COUNTS; global SP_COUNTS; %These counts are for populations
global CQ_SUMCOUNTS; global SP_SUMCOUNTS; %not for individuals
global LOGDIFF;
clearGlobalVars;
noalle = c.noalle;
adjprior = c.adjprior; %priorTerm = c.priorTerm;
rowsFromInd = c.rowsFromInd;
counts_cq = c.counts_cq; adjprior_cq = c.adjprior_cq;
counts_sp = c.counts_sp; adjprior_sp = c.adjprior_sp;
if isfield(c,'dist')
dist = c.dist; Z = c.Z;
end
clear c;
ninds = size(counts_cq,3);
if nargin < 2
npopstext = [];
if ninds>20
default = 20;
else
default = floor(ninds/2);
end
teksti = {'Number of populations:', ...
'Number of runs:'};
def = {num2str(default), '1'};
npopstextExtra = inputdlg(teksti ,...
'Input parameters for the computation algorithm',1,def);
if isempty(npopstextExtra) % cancel has been pressed
dispCancel
logml = 1; partitionSummary=1; npops=1;
return
end
npopstext = npopstextExtra{1};
nrunstext = npopstextExtra{2};
clear teksti npopstextExtra;
if isempty(npopstext)
return
else
npopsTable = str2num(npopstext);
npops = npopsTable(1);
if npops==1
logml = 1; partitionSummary=1; npops=1;
return
end
nrunsTable = str2num(nrunstext);
nruns = nrunsTable(1);
end
end
logmlBest = -1e50;
partitionSummary = -1e50*ones(100,2); % 100 best partitions (npops and logml)
partitionSummary(:,1) = zeros(100,1);
worstLogml = -1e50; worstIndex = 1;
for run = 1:nruns
dispLine;
disp(['Run ' num2str(run) '/' num2str(nruns) ...
', maximum number of populations ' num2str(npops) '.']);
initialPartition = admixture_initialization(npops, Z);
PARTITION = initialPartition;
[cq_counts, cq_sumcounts] = initialCounts(counts_cq);
% clear counts_cq;
CQ_COUNTS = cq_counts; clear cq_counts;
CQ_SUMCOUNTS = cq_sumcounts; clear cq_sumcounts;
[sp_counts, sp_sumcounts] = initialCounts(counts_sp);
% clear counts_sp;
SP_COUNTS = sp_counts; clear sp_counts;
SP_SUMCOUNTS = sp_sumcounts; clear sp_sumcounts;
logml = computeLogml(adjprior_cq, adjprior_sp);
POP_LOGML = computePopulationLogml(1:npops,adjprior_cq, adjprior_sp);
if logml>worstLogml
[partitionSummary, added] = addToSummary(logml, partitionSummary, worstIndex);
if (added==1)
[worstLogml, worstIndex] = min(partitionSummary(:,2));
end
end
%%%%
% Finding the best partition with the greedy search algorithm
%%%%
nRoundTypes = 7;
tested = zeros(nRoundTypes,1);
roundTypes = [1 1];
ready = 0; phase = 1;
ninds = length(PARTITION); % number of individuals
LOGDIFF = repmat(-Inf,ninds,npops);
disp(' ');
disp(['Mixture analysis started with initial ' num2str(npops) ' populations.']);
while ready ~= 1
changesMade = 0;
disp(['Performing steps: ' num2str(roundTypes)]);
for n = 1:length(roundTypes)
round = roundTypes(n);
% pack;
if tested(round) == 1
elseif round==1 % Moving one individual to another population
inds = randperm(ninds); % random order
changesMadeNow = 0;
for ind = inds
i1 = PARTITION(ind);
if length(find(PARTITION==i1))>1
indCqCounts = uint16(counts_cq(:,:,ind));
indSpCounts = uint16(counts_sp(:,:,ind));
changesInLogml = computeChanges(ind, adjprior_cq, ...
adjprior_sp, indCqCounts, indSpCounts);
[maxChange, i2] = max(changesInLogml);
if (i1~=i2 && maxChange>1e-5)
% Individual is moved
changesMade = 1;
if changesMadeNow == 0
disp('Action 1');
changesMadeNow = 1;
tested = zeros(nRoundTypes,1);
end
updateGlobalVariables(ind, i2, indCqCounts, ...
indSpCounts, adjprior_cq, adjprior_sp);
logml = logml+maxChange;
if logml>worstLogml
[partitionSummary, added] = addToSummary(logml, partitionSummary, worstIndex);
if (added==1)
[worstLogml, worstIndex] = min(partitionSummary(:,2));
end
end
end
end
end
if changesMadeNow == 0
tested(round) = 1;
end
elseif round==2 % Merging two populations and splitting the result
maxChange = -1e50;
poplogml = POP_LOGML;
partition = PARTITION;
cq_counts = CQ_COUNTS;
sp_counts = SP_COUNTS;
cq_sumcounts = CQ_SUMCOUNTS;
sp_sumcounts = SP_SUMCOUNTS;
logdiff = LOGDIFF;
% Two populations are merged first
for pop = 1:npops
changesInLogml = computeChanges2(pop, adjprior_cq, adjprior_sp);
changesInLogml(pop)=-1e50;
[biggest, index] = max(changesInLogml);
if biggest>maxChange
maxChange = biggest;
i1 = pop;
i2 = index;
end
end
totalChange = maxChange;
updateGlobalVariables2(i1, i2, adjprior_cq, adjprior_sp);
% A population is split in two
emptyPop = i1;
maxChange = -1e50;
for pop = 1:npops
inds2 = find(PARTITION==pop);
ninds2 = length(inds2);
if ninds2>1
% Computing the distance between individuals inds2
dist2 = laskeOsaDist(inds2, dist, ninds);
Z2 = computeLinkage(dist2');
npops2 = 2;
T2 = cluster_own(Z2, npops2);
changesInLogml = computeChanges3(T2, inds2, pop, ...
counts_cq, counts_sp, adjprior_cq, adjprior_sp);
biggest = changesInLogml(1,emptyPop);
if biggest > maxChange
maxChange = biggest;
movingInds = inds2(logical(T2==1));
end
end
end
indCqCounts = uint16(sum(counts_cq(:,:,movingInds),3));
indSpCounts = uint16(sum(counts_sp(:,:,movingInds),3));
updateGlobalVariables3(movingInds, emptyPop,indCqCounts, ...
indSpCounts, adjprior_cq, adjprior_sp);
totalChange = totalChange + maxChange;
% Individuals are moved between populations
inds = randperm(ninds); % random order
for ind = inds
i1 = PARTITION(ind);
if length(find(PARTITION==i1))>1
indCqCounts = uint16(counts_cq(:,:,ind));
indSpCounts = uint16(counts_sp(:,:,ind));
changesInLogml = computeChanges(ind, adjprior_cq, ...
adjprior_sp, indCqCounts, indSpCounts);
[maxChange, i2] = max(changesInLogml);
if (i1~=i2 && maxChange>1e-5)
updateGlobalVariables(ind, i2, indCqCounts, ...
indSpCounts, adjprior_cq, adjprior_sp);
totalChange = totalChange + maxChange;
end
end
end
if totalChange > 1e-5
disp('Action 2');
logml = logml + totalChange;
tested = zeros(nRoundTypes,1);
if logml>worstLogml
[partitionSummary, added] = addToSummary(logml, partitionSummary, worstIndex);
if (added==1)
[worstLogml, worstIndex] = min(partitionSummary(:,2));
end
end
else
PARTITION = partition;
POP_LOGML = poplogml;
CQ_COUNTS = cq_counts;
SP_COUNTS = sp_counts;
CQ_SUMCOUNTS = cq_sumcounts;
SP_SUMCOUNTS = sp_sumcounts;
LOGDIFF = logdiff;
tested(round) = 1;
end
elseif round==3 || round==4 % Splitting population into smaller groups
maxChange = 0;
for pop = 1:npops
inds2 = find(PARTITION==pop);
ninds2 = length(inds2);
if ninds2>5
% Computing the distance between individuals inds2
dist2 = laskeOsaDist(inds2, dist, ninds);
Z2 = computeLinkage(dist2');
% Number of groups:
if round==3
npops2 = max(min(20, floor(ninds2/2)),2);
elseif round==4
npops2 = 2;
end
T2 = cluster_own(Z2, npops2);
changesInLogml = computeChanges3(T2, inds2, pop, ...
counts_cq, counts_sp, adjprior_cq, adjprior_sp);
[biggest, index] = max(changesInLogml(1:end));
if biggest > maxChange
maxChange = biggest;
movingGroup = rem(index,npops2); % The group, which is moved
if movingGroup==0, movingGroup = npops2; end
movingInds = inds2(logical(T2==movingGroup));
i2 = ceil(index/npops2); % pop where movingGroup would be moved
end
end
end
if maxChange>1e-5
changesMade = 1;
tested = zeros(nRoundTypes,1);
if round==3
disp('Action 3');
else
disp('Action 4');
end
indCqCounts = uint16(sum(counts_cq(:,:,movingInds),3));
indSpCounts = uint16(sum(counts_sp(:,:,movingInds),3));
updateGlobalVariables3(movingInds, i2,indCqCounts, ...
indSpCounts, adjprior_cq, adjprior_sp);
logml = logml + maxChange;
if logml>worstLogml
[partitionSummary, added] = addToSummary(logml, partitionSummary, worstIndex);
if (added==1)
[worstLogml, worstIndex] = min(partitionSummary(:,2));
end
end
else
tested(round) = 1;
end
elseif round == 5 || round == 6
%Moving individuals out of population until positive change
%in logml has occured
pop=0;
changesMadeNow = 0;
%Saving old values
poplogml = POP_LOGML;
partition = PARTITION;
cq_counts = CQ_COUNTS;
sp_counts = SP_COUNTS;
cq_sumcounts = CQ_SUMCOUNTS;
sp_sumcounts = SP_SUMCOUNTS;
logdiff = LOGDIFF;
while (pop < npops && changesMadeNow == 0)
pop = pop+1;
totalChangeInLogml = 0;
inds = find(PARTITION==pop);
if round == 5
%Random order
aputaulu = [inds rand(length(inds),1)];
aputaulu = sortrows(aputaulu,2);
inds = aputaulu(:,1)';
elseif round == 6
inds = returnInOrder(inds, pop, counts_cq, counts_sp, ...
adjprior_cq, adjprior_sp);
end
i=0;
while (length(inds)>0 && i<length(inds) - 1)
i = i+1;
ind = inds(i);
indCqCounts = uint16(counts_cq(:,:,ind));
indSpCounts = uint16(counts_sp(:,:,ind));
changesInLogml = computeChanges(ind, adjprior_cq, ...
adjprior_sp, indCqCounts, indSpCounts);
changesInLogml(pop) = -1e50; % Varmasti ei suurin!!!
[maxChange, i2] = max(changesInLogml);
updateGlobalVariables(ind, i2, indCqCounts, ...
indSpCounts, adjprior_cq, adjprior_sp);
totalChangeInLogml = totalChangeInLogml+maxChange;
logml = logml+maxChange;
if round == 6
% Stop immediatly when change in logml is
% positive
if totalChangeInLogml > 1e-5
i=length(inds);
end
end
end
if totalChangeInLogml>1e-5
if round == 5
disp('Action 5');
elseif round == 6
disp('Action 6');
end
tested = zeros(nRoundTypes,1);
changesMadeNow=1;
changesMade = 1;
if logml>worstLogml
[partitionSummary, added] = addToSummary(logml, partitionSummary, worstIndex);
if (added==1)
[worstLogml, worstIndex] = min(partitionSummary(:,2));
end
end
else
% No better partition was found, restoring the old
% values
PARTITION = partition;
POP_LOGML = poplogml;
CQ_COUNTS = cq_counts;
SP_COUNTS = sp_counts;
CQ_SUMCOUNTS = cq_sumcounts;
SP_SUMCOUNTS = sp_sumcounts;
LOGDIFF = logdiff;
logml = logml - totalChangeInLogml;
end
end
clear partition; clear poplogml;
if changesMadeNow == 0
tested(round) = 1;
end
elseif round == 7
emptyPop = npops + 1;
j = 0;
pops = randperm(npops);
changesMadeNow = 0;
while (j < npops)
j = j +1;
pop = pops(j);
inds2 = find(PARTITION == pop);
ninds2 = length(inds2);
if ninds2 > 5
partition = PARTITION;
cq_sumcounts = CQ_SUMCOUNTS;
cq_counts = CQ_COUNTS;
sp_sumcounts = SP_SUMCOUNTS;
sp_counts = SP_COUNTS;
poplogml = POP_LOGML;
logdiff = LOGDIFF;
% A new population is created temporarily
npops = npops + 1;
POP_LOGML(npops) = 0;
CQ_COUNTS(:,:,npops) = zeros(size(CQ_COUNTS(:,:,1)));
CQ_SUMCOUNTS(npops,:) = zeros(size(CQ_SUMCOUNTS(1,:)));
SP_COUNTS(:,:,npops) = zeros(size(SP_COUNTS(:,:,1)));
SP_SUMCOUNTS(npops,:) = zeros(size(SP_SUMCOUNTS(1,:)));
dist2 = laskeOsaDist(inds2, dist, ninds);
Z2 = computeLinkage(dist2');
T2 = cluster_own(Z2, 2);
% movingInds = inds2(find(T2 == 1));
movingInds = inds2(logical(T2 == 1));
changesInLogml = computeChanges3(T2, inds2, pop, ...
counts_cq, counts_sp, adjprior_cq, adjprior_sp);
totalChangeInLogml = changesInLogml(1, emptyPop);
indCqCounts = uint16(sum(counts_cq(:,:,movingInds),3));
indSpCounts = uint16(sum(counts_sp(:,:,movingInds),3));
updateGlobalVariables3(movingInds, emptyPop,indCqCounts, ...
indSpCounts, adjprior_cq, adjprior_sp);
% Individuals are moved between populations
inds = randperm(ninds); % random order
for ind = inds
i1 = PARTITION(ind);
if length(find(PARTITION==i1))>1
indCqCounts = uint16(counts_cq(:,:,ind));
indSpCounts = uint16(counts_sp(:,:,ind));
changesInLogml = computeChanges(ind, adjprior_cq, ...
adjprior_sp, indCqCounts, indSpCounts);
[maxChange, i2] = max(changesInLogml);
if (i1~=i2 && maxChange>1e-5)
updateGlobalVariables(ind, i2, indCqCounts, ...
indSpCounts, adjprior_cq, adjprior_sp);
totalChangeInLogml = totalChangeInLogml + maxChange;
end
end
end
% Two populations are merged
if length(find(any(CQ_SUMCOUNTS,2))) == npops
maxChange = -1e50;
for pop = 1:npops
changesInLogml = computeChanges2(pop, adjprior_cq, adjprior_sp);
changesInLogml(pop)=-1e50;
[biggest, index] = max(changesInLogml);
if biggest>maxChange
maxChange = biggest;
i1 = pop;
i2 = index;
end
end
totalChangeInLogml = totalChangeInLogml + maxChange;
updateGlobalVariables2(i1, i2, adjprior_cq, adjprior_sp);
end
if totalChangeInLogml > 1e-5
changesMade = 1;
logml = logml + totalChangeInLogml;
npops = removeEmptyPops; % The temporary population is removed
POP_LOGML = computePopulationLogml(1:npops, adjprior_cq, adjprior_sp);
if logml>worstLogml
[partitionSummary, added] = addToSummary(logml, partitionSummary, worstIndex);
if (added==1)
[worstLogml, worstIndex] = min(partitionSummary(:,2));
end
end
if changesMadeNow == 0
disp('Action 7');
changesMadeNow = 1;
end
changesMadeNow = 1;
tested = zeros(nRoundTypes, 1);
j = npops;
else
% No better partition was found, restoring the old
% values
PARTITION = partition;
POP_LOGML = poplogml;
CQ_COUNTS = cq_counts;
SP_COUNTS = sp_counts;
CQ_SUMCOUNTS = cq_sumcounts;
SP_SUMCOUNTS = sp_sumcounts;
LOGDIFF = logdiff;
npops = npops-1;
end
end
end
if changesMadeNow == 0
tested(round) = 1;
end
end
end
if changesMade == 0
if phase==1
phase = 2;
elseif phase==2
phase = 3;
elseif phase==3
phase = 4;
elseif phase==4;
phase = 5;
elseif phase==5
ready = 1;
end
else
changesMade = 0;
end
if ready==0
if phase==1
roundTypes=[1];
elseif phase==2
roundTypes=[2];
elseif phase==3
roundTypes=[5 5 7];
elseif phase==4
roundTypes=[4 3 1 1];
elseif phase==5
roundTypes=[6 2 7 3 4 1];
end
end
end
% Saving results
npops = removeEmptyPops;
POP_LOGML = computePopulationLogml(1:npops, adjprior_cq, adjprior_sp);
disp(['Found partition with ' num2str(npops) ' populations.']);
disp(['Log(ml) = ' num2str(logml)]);
disp(' ');
if logml>logmlBest
% Updating the best found partition
logmlBest = logml;
npopsBest = npops;
partitionBest = PARTITION;
cq_countsBest = CQ_COUNTS;
sp_countsBest = SP_COUNTS;
cq_sumcountsBest = CQ_SUMCOUNTS;
sp_sumcountsBest = SP_SUMCOUNTS;
pop_logmlBest = POP_LOGML;
logdiffbest = LOGDIFF;
end
end
logml = logmlBest;
npops = npopsBest;
PARTITION = partitionBest;
CQ_COUNTS = cq_countsBest;
SP_COUNTS = sp_countsBest;
CQ_SUMCOUNTS = cq_sumcountsBest;
SP_SUMCOUNTS = sp_sumcountsBest;
POP_LOGML = pop_logmlBest;
LOGDIFF = logdiffbest;
%--------------------------------------------------------------------------
% The next three functions are for computing the initial partition
% according to the distance between the individuals
function initial_partition=admixture_initialization(nclusters,Z)
T=cluster_own(Z,nclusters);
initial_partition=T;
%--------------------------------------------------------------------------
function T = cluster_own(Z,nclust)
% true=logical(1);
% false=logical(0);
maxclust = nclust;
% Start of algorithm
m = size(Z,1)+1;
T = zeros(m,1);
% maximum number of clusters based on inconsistency
if m <= maxclust
T = (1:m)';
elseif maxclust==1
T = ones(m,1);
else
clsnum = 1;
for k = (m-maxclust+1):(m-1)
i = Z(k,1); % left tree
if i <= m % original node, no leafs
T(i) = clsnum;
clsnum = clsnum + 1;
elseif i < (2*m-maxclust+1) % created before cutoff, search down the tree
T = clusternum(Z, T, i-m, clsnum);
clsnum = clsnum + 1;
end
i = Z(k,2); % right tree
if i <= m % original node, no leafs
T(i) = clsnum;
clsnum = clsnum + 1;
elseif i < (2*m-maxclust+1) % created before cutoff, search down the tree
T = clusternum(Z, T, i-m, clsnum);
clsnum = clsnum + 1;
end
end
end
function T = clusternum(X, T, k, c)
m = size(X,1)+1;
while(~isempty(k))
% Get the children of nodes at this level
children = X(k,1:2);
children = children(:);
% Assign this node number to leaf children
t = (children<=m);
T(children(t)) = c;
% Move to next level
k = children(~t) - m;
end
%--------------------------------------------------------------------------
function dist2 = laskeOsaDist(inds2, dist, ninds)
% Muodostaa dist vektorista osavektorin, joka sis<EFBFBD>lt<EFBFBD><EFBFBD> yksil<EFBFBD>iden inds2
% v<EFBFBD>liset et<EFBFBD>isyydet. ninds=kaikkien yksil<EFBFBD>iden lukum<EFBFBD><EFBFBD>r<EFBFBD>.
ninds2 = length(inds2);
apu = zeros(nchoosek(ninds2,2),2);
rivi = 1;
for i=1:ninds2-1
for j=i+1:ninds2
apu(rivi, 1) = inds2(i);
apu(rivi, 2) = inds2(j);
rivi = rivi+1;
end
end
apu = (apu(:,1)-1).*ninds - apu(:,1) ./ 2 .* (apu(:,1)-1) + (apu(:,2)-apu(:,1));
dist2 = dist(apu);
%--------------------------------------------------------------------------
function Z = computeLinkage(Y, method)
[k, n] = size(Y);
m = (1+sqrt(1+8*n))/2;
if k ~= 1 || m ~= fix(m)
error('The first input has to match the output of the PDIST function in size.');
end
if nargin == 1 % set default switch to be 'co'
method = 'co';
end
method = lower(method(1:2)); % simplify the switch string.
% monotonic = 1;
Z = zeros(m-1,3); % allocate the output matrix.
N = zeros(1,2*m-1);
N(1:m) = 1;
n = m; % since m is changing, we need to save m in n.
R = 1:n;
for s = 1:(n-1)
X = Y;
[v, k] = min(X);
i = floor(m+1/2-sqrt(m^2-m+1/4-2*(k-1)));
j = k - (i-1)*(m-i/2)+i;
Z(s,:) = [R(i) R(j) v]; % update one more row to the output matrix A
I1 = 1:(i-1); I2 = (i+1):(j-1); I3 = (j+1):m; % these are temp variables.
U = [I1 I2 I3];
I = [I1.*(m-(I1+1)/2)-m+i i*(m-(i+1)/2)-m+I2 i*(m-(i+1)/2)-m+I3];
J = [I1.*(m-(I1+1)/2)-m+j I2.*(m-(I2+1)/2)-m+j j*(m-(j+1)/2)-m+I3];
switch method
case 'si' %single linkage
Y(I) = min(Y(I),Y(J));
case 'av' % average linkage
Y(I) = Y(I) + Y(J);
case 'co' %complete linkage
Y(I) = max(Y(I),Y(J));
case 'ce' % centroid linkage
K = N(R(i))+N(R(j));
Y(I) = (N(R(i)).*Y(I)+N(R(j)).*Y(J)-(N(R(i)).*N(R(j))*v^2)./K)./K;
case 'wa'
Y(I) = ((N(R(U))+N(R(i))).*Y(I) + (N(R(U))+N(R(j))).*Y(J) - ...
N(R(U))*v)./(N(R(i))+N(R(j))+N(R(U)));
end
J = [J i*(m-(i+1)/2)-m+j];
Y(J) = []; % no need for the cluster information about j.
% update m, N, R
m = m-1;
N(n+s) = N(R(i)) + N(R(j));
R(i) = n+s;
R(j:(n-1))=R((j+1):n);
end
%--------------------------------------------------------------------------
function changes = computeChanges(ind, adjprior_cq, adjprior_sp, ...
indCqCounts, indSpCounts)
% Computes changes in log-marginal likelihood if individual ind is
% moved to another population
%
% Input:
% ind - the individual to be moved
% adjprior_cq & _sp - adjpriors for cliques and separators
% indCqCounts, indSpCounts - counts for individual ind
%
% Output:
% changes - table of size 1*npops. changes(i) = difference in logml if
% ind is move to population i.
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
global PARTITION; global POP_LOGML;
global LOGDIFF;
npops = size(CQ_COUNTS,3);
changes = LOGDIFF(ind,:);
i1 = PARTITION(ind);
i1_logml = POP_LOGML(i1);
changes(i1) = 0;
sumCq = uint16(sum(indCqCounts,1));
sumSp = uint16(sum(indSpCounts,1));
CQ_COUNTS(:,:,i1) = CQ_COUNTS(:,:,i1)-indCqCounts;
CQ_SUMCOUNTS(i1,:) = CQ_SUMCOUNTS(i1,:)-sumCq;
SP_COUNTS(:,:,i1) = SP_COUNTS(:,:,i1)-indSpCounts;
SP_SUMCOUNTS(i1,:) = SP_SUMCOUNTS(i1,:)-sumSp;
new_i1_logml = computePopulationLogml(i1, adjprior_cq, adjprior_sp);
CQ_COUNTS(:,:,i1) = CQ_COUNTS(:,:,i1)+indCqCounts;
CQ_SUMCOUNTS(i1,:) = CQ_SUMCOUNTS(i1,:)+sumCq;
SP_COUNTS(:,:,i1) = SP_COUNTS(:,:,i1)+indSpCounts;
SP_SUMCOUNTS(i1,:) = SP_SUMCOUNTS(i1,:)+sumSp;
i2 = find(changes==-Inf);
i2 = setdiff(i2,i1);
i2_logml = POP_LOGML(i2);
ni2 = length(i2);
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)+repmat(indCqCounts, [1 1 ni2]);
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)+repmat(sumCq,[ni2 1]);
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)+repmat(indSpCounts, [1 1 ni2]);
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:) + repmat(sumSp,[ni2 1]);
new_i2_logml = computePopulationLogml(i2, adjprior_cq, adjprior_sp);
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)-repmat(indCqCounts, [1 1 ni2]);
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)-repmat(sumCq,[ni2 1]);
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)-repmat(indSpCounts, [1 1 ni2]);
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:) - repmat(sumSp,[ni2 1]);
% a = repmat(sumSp,[npops-1 1]);
changes(i2) = new_i1_logml - i1_logml ...
+ new_i2_logml - i2_logml;
LOGDIFF(ind,:) = changes;
%------------------------------------------------------------------------------------
function changes = computeChanges2(i1, adjprior_cq, adjprior_sp)
% Computes changes in log marginal likelihood if population i1 is combined
% with another population
%
% Input:
% i1 - the population to be combined
% adjprior_cq & _sp - adjpriors for cliques and separators
%
% Output:
% changes - table of size 1*npops. changes(i) = difference in logml if
% i1 is combined with population i.
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
global POP_LOGML;
npops = size(CQ_COUNTS,3);
changes = zeros(npops,1);
i1_logml = POP_LOGML(i1);
indCqCounts = CQ_COUNTS(:,:,i1);
indSpCounts = SP_COUNTS(:,:,i1);
sumCq = uint16(sum(indCqCounts,1));
sumSp = uint16(sum(indSpCounts,1));
new_i1_logml = 0;
i2 = [1:i1-1 , i1+1:npops];
i2_logml = POP_LOGML(i2);
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)+repmat(indCqCounts, [1 1 npops-1]);
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)+repmat(sumCq,[npops-1 1]);
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)+repmat(indSpCounts, [1 1 npops-1]);
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)+ repmat(sumSp,[npops-1 1]);
% a = repmat(sumSp,[npops-1 1]);
% if ~any(sumSp)
% a(:,[1:size(a,2)])=[];
% end
% SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)+ a ;
new_i2_logml = computePopulationLogml(i2, adjprior_cq, adjprior_sp);
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)-repmat(indCqCounts, [1 1 npops-1]);
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)-repmat(sumCq,[npops-1 1]);
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)-repmat(indSpCounts, [1 1 npops-1]);
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)- repmat(sumSp,[npops-1 1]);
changes(i2) = new_i1_logml - i1_logml ...
+ new_i2_logml - i2_logml;
%------------------------------------------------------------------------------------
function changes = computeChanges3(T2, inds2, i1, counts_cq, counts_sp, ...
adjprior_cq, adjprior_sp)
% Computes changes in log marginal likelihood if subpopulation of i2 is
% moved to another population
%
% Input:
% T2 - partition of inds2 to subpopulations
% inds2 - individuals in population i1
% i2
% counts_cq, counts_sp - counts for individuals
%
% Output:
% changes - table of size length(unique(T2))*npops.
% changes(i,j) = difference in logml if subpopulation inds2(find(T2==i)) of
% i2 is moved to population j
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
global POP_LOGML;
npops = size(CQ_COUNTS,3);
npops2 = length(unique(T2));
changes = zeros(npops2,npops);
%cq_counts = CQ_COUNTS;
%sp_counts = SP_COUNTS;
%cq_sumcounts = CQ_SUMCOUNTS;
%sp_sumcounts = SP_SUMCOUNTS;
i1_logml = POP_LOGML(i1);
for pop2 = 1:npops2
% inds = inds2(find(T2==pop2));
inds = inds2(logical(T2==pop2));
ninds = length(inds);
if ninds>0
indCqCounts = uint16(sum(counts_cq(:,:,inds),3));
indSpCounts = uint16(sum(counts_sp(:,:,inds),3));
sumCq = uint16(sum(indCqCounts,1));
sumSp = uint16(sum(indSpCounts,1));
CQ_COUNTS(:,:,i1) = CQ_COUNTS(:,:,i1)-indCqCounts;
CQ_SUMCOUNTS(i1,:) = CQ_SUMCOUNTS(i1,:)-sumCq;
SP_COUNTS(:,:,i1) = SP_COUNTS(:,:,i1)-indSpCounts;
SP_SUMCOUNTS(i1,:) = SP_SUMCOUNTS(i1,:)-sumSp;
new_i1_logml = computePopulationLogml(i1, adjprior_cq, adjprior_sp);
CQ_COUNTS(:,:,i1) = CQ_COUNTS(:,:,i1)+indCqCounts;
CQ_SUMCOUNTS(i1,:) = CQ_SUMCOUNTS(i1,:)+sumCq;
SP_COUNTS(:,:,i1) = SP_COUNTS(:,:,i1)+indSpCounts;
SP_SUMCOUNTS(i1,:) = SP_SUMCOUNTS(i1,:)+sumSp;
i2 = [1:i1-1 , i1+1:npops];
i2_logml = POP_LOGML(i2)';
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)+repmat(indCqCounts, [1 1 npops-1]);
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)+repmat(sumCq,[npops-1 1]);
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)+repmat(indSpCounts, [1 1 npops-1]);
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)+ repmat(sumSp,[npops-1 1]);
new_i2_logml = computePopulationLogml(i2, adjprior_cq, adjprior_sp)';
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)-repmat(indCqCounts, [1 1 npops-1]);
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)-repmat(sumCq,[npops-1 1]);
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)-repmat(indSpCounts, [1 1 npops-1]);
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)- repmat(sumSp,[npops-1 1]);
changes(pop2,i2) = new_i1_logml - i1_logml ...
+ new_i2_logml - i2_logml;
end
end
%--------------------------------------------------------------------------
function changes = computeChanges5(inds, i1, i2, counts_cq, counts_sp, ...
adjprior_cq, adjprior_sp)
% Computes change in logml if individual of inds is moved between
% populations i1 and i2
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
global POP_LOGML; global PARTITION;
ninds = length(inds);
changes = zeros(ninds,1);
i1_logml = POP_LOGML(i1);
i2_logml = POP_LOGML(i2);
for i = 1:ninds
ind = inds(i);
if PARTITION(ind)==i1
pop1 = i1; %from
pop2 = i2; %to
else
pop1 = i2;
pop2 = i1;
end
indCqCounts = uint16(counts_cq(:,:,ind));
indSpCounts = uint16(counts_sp(:,:,ind));
sumCq = uint16(sum(indCqCounts,1));
sumSp = uint16(sum(indSpCounts,1));
CQ_COUNTS(:,:,pop1) = CQ_COUNTS(:,:,pop1)-indCqCounts;
CQ_SUMCOUNTS(pop1,:) = CQ_SUMCOUNTS(pop1,:)-sumCq;
SP_COUNTS(:,:,pop1) = SP_COUNTS(:,:,pop1)-indSpCounts;
SP_SUMCOUNTS(pop1,:) = SP_SUMCOUNTS(pop1,:) - sumSp;
CQ_COUNTS(:,:,pop2) = CQ_COUNTS(:,:,pop2)+indCqCounts;
CQ_SUMCOUNTS(pop2,:) = CQ_SUMCOUNTS(pop2,:)+sumCq;
SP_COUNTS(:,:,pop2) = SP_COUNTS(:,:,pop2)+indSpCounts;
SP_SUMCOUNTS(pop2,:) = SP_SUMCOUNTS(pop2,:) + sumSp;
new_logmls = computePopulationLogml([i1 i2], adjprior_cq, adjprior_sp);
changes(i) = sum(new_logmls);
CQ_COUNTS(:,:,pop1) = CQ_COUNTS(:,:,pop1)+indCqCounts;
CQ_SUMCOUNTS(pop1,:) = CQ_SUMCOUNTS(pop1,:)+sumCq;
SP_COUNTS(:,:,pop1) = SP_COUNTS(:,:,pop1)+indSpCounts;
SP_SUMCOUNTS(pop1,:) = SP_SUMCOUNTS(pop1,:)+sumSp;
CQ_COUNTS(:,:,pop2) = CQ_COUNTS(:,:,pop2)-indCqCounts;
CQ_SUMCOUNTS(pop2,:) = CQ_SUMCOUNTS(pop2,:)-sumCq;
SP_COUNTS(:,:,pop2) = SP_COUNTS(:,:,pop2)-indSpCounts;
SP_SUMCOUNTS(pop2,:) = SP_SUMCOUNTS(pop2,:)-sumSp;
end
changes = changes - i1_logml - i2_logml;
%-------------------------------------------------------------------------------------
function updateGlobalVariables(ind, i2, indCqCounts, indSpCounts, ...
adjprior_cq, adjprior_sp)
% Updates global variables when individual ind is moved to population i2
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
global PARTITION; global POP_LOGML;
global LOGDIFF;
i1 = PARTITION(ind);
PARTITION(ind)=i2;
sumCq = uint16(sum(indCqCounts,1));
sumSp = uint16(sum(indSpCounts,1));
CQ_COUNTS(:,:,i1) = CQ_COUNTS(:,:,i1)-indCqCounts;
CQ_SUMCOUNTS(i1,:) = CQ_SUMCOUNTS(i1,:)-sumCq;
SP_COUNTS(:,:,i1) = SP_COUNTS(:,:,i1)-indSpCounts;
SP_SUMCOUNTS(i1,:) = SP_SUMCOUNTS(i1,:)-sumSp;
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)+indCqCounts;
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)+sumCq;
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)+indSpCounts;
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)+sumSp;
POP_LOGML([i1 i2]) = computePopulationLogml([i1 i2], adjprior_cq, adjprior_sp);
LOGDIFF(:,[i1 i2]) = -Inf;
inx = [find(PARTITION==i1); find(PARTITION==i2)];
LOGDIFF(inx,:) = -Inf;
%---------------------------------------------------------------------------------
function updateGlobalVariables2(i1, i2, adjprior_cq, adjprior_sp)
% Updates global variables when all individuals from population i1 are moved
% to population i2
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
global PARTITION; global POP_LOGML;
global LOGDIFF;
% inds = find(PARTITION==i1);
% PARTITION(inds) = i2;
PARTITION(logical(PARTITION==i1)) = i2;
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)+CQ_COUNTS(:,:,i1);
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)+CQ_SUMCOUNTS(i1,:);
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)+SP_COUNTS(:,:,i1);
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)+SP_SUMCOUNTS(i1,:);
CQ_COUNTS(:,:,i1) = 0;
CQ_SUMCOUNTS(i1,:) = 0;
SP_COUNTS(:,:,i1) = 0;
SP_SUMCOUNTS(i1,:) = 0;
POP_LOGML(i1) = 0;
POP_LOGML(i2) = computePopulationLogml(i2, adjprior_cq, adjprior_sp);
LOGDIFF(:,[i1 i2]) = -Inf;
inx = [find(PARTITION==i1); find(PARTITION==i2)];
LOGDIFF(inx,:) = -Inf;
%------------------------------------------------------------------------------------
function updateGlobalVariables3(muuttuvat, i2, indCqCounts, indSpCounts, ...
adjprior_cq, adjprior_sp)
% Updates global variables when individuals muuttuvat are moved to
% population i2
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
global PARTITION; global POP_LOGML;
global LOGDIFF;
i1 = PARTITION(muuttuvat(1));
PARTITION(muuttuvat) = i2;
sumCq = uint16(sum(indCqCounts,1));
sumSp = uint16(sum(indSpCounts,1));
CQ_COUNTS(:,:,i1) = CQ_COUNTS(:,:,i1)-indCqCounts;
CQ_SUMCOUNTS(i1,:) = CQ_SUMCOUNTS(i1,:)-sumCq;
SP_COUNTS(:,:,i1) = SP_COUNTS(:,:,i1)-indSpCounts;
SP_SUMCOUNTS(i1,:) = SP_SUMCOUNTS(i1,:)-sumSp;
CQ_COUNTS(:,:,i2) = CQ_COUNTS(:,:,i2)+indCqCounts;
CQ_SUMCOUNTS(i2,:) = CQ_SUMCOUNTS(i2,:)+sumCq;
SP_COUNTS(:,:,i2) = SP_COUNTS(:,:,i2)+indSpCounts;
SP_SUMCOUNTS(i2,:) = SP_SUMCOUNTS(i2,:)+sumSp;
POP_LOGML([i1 i2]) = computePopulationLogml([i1 i2], adjprior_cq, adjprior_sp);
LOGDIFF(:,[i1 i2]) = -Inf;
inx = [find(PARTITION==i1); find(PARTITION==i2)];
LOGDIFF(inx,:) = -Inf;
%----------------------------------------------------------------------
function inds = returnInOrder(inds, pop, counts_cq, counts_sp, ...
adjprior_cq, adjprior_sp)
% Returns individuals inds in order according to the change in the logml if
% they are moved out of the population pop
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
ninds = length(inds);
apuTaulu = [inds, zeros(ninds,1)];
for i=1:ninds
ind = inds(i);
indCqCounts = uint16(counts_cq(:,:,ind));
indSpCounts = uint16(counts_sp(:,:,ind));
sumCq = uint16(sum(indCqCounts,1));
sumSp = uint16(sum(indSpCounts,1));
CQ_COUNTS(:,:,pop) = CQ_COUNTS(:,:,pop)-indCqCounts;
CQ_SUMCOUNTS(pop,:) = CQ_SUMCOUNTS(pop,:)-sumCq;
SP_COUNTS(:,:,pop) = SP_COUNTS(:,:,pop)-indSpCounts;
SP_SUMCOUNTS(pop,:) = SP_SUMCOUNTS(pop,:)-sumSp;
apuTaulu(i, 2) = computePopulationLogml(pop, adjprior_cq, adjprior_sp);
CQ_COUNTS(:,:,pop) = CQ_COUNTS(:,:,pop)+indCqCounts;
CQ_SUMCOUNTS(pop,:) = CQ_SUMCOUNTS(pop,:)+sumCq;
SP_COUNTS(:,:,pop) = SP_COUNTS(:,:,pop)+indSpCounts;
SP_SUMCOUNTS(pop,:) = SP_SUMCOUNTS(pop,:)+sumSp;
end
apuTaulu = sortrows(apuTaulu,2);
inds = apuTaulu(ninds:-1:1,1);
%--------------------------------------------------------------------------
function clearGlobalVars
global CQ_COUNTS; CQ_COUNTS = [];
global CQ_SUMCOUNTS; CQ_SUMCOUNTS = [];
global SP_COUNTS; SP_COUNTS = [];
global SP_SUMCOUNTS; SP_SUMCOUNTS = [];
global PARTITION; PARTITION = [];
global POP_LOGML; POP_LOGML = [];
global LOGDIFF; LOGDIFF = [];
%--------------------------------------------------------------------------
function npops = removeEmptyPops
% Removes empty pops from all global COUNTS variables.
% Updates PARTITION and npops
global CQ_COUNTS;
global CQ_SUMCOUNTS;
global SP_COUNTS;
global SP_SUMCOUNTS;
global PARTITION;
global LOGDIFF;
notEmpty = find(any(CQ_SUMCOUNTS,2));
CQ_COUNTS = CQ_COUNTS(:,:,notEmpty);
CQ_SUMCOUNTS = CQ_SUMCOUNTS(notEmpty,:);
SP_COUNTS = SP_COUNTS(:,:,notEmpty);
SP_SUMCOUNTS = SP_SUMCOUNTS(notEmpty,:);
LOGDIFF = LOGDIFF(:,notEmpty);
for n=1:length(notEmpty)
% apu = find(PARTITION==notEmpty(n));
% PARTITION(apu)=n;
PARTITION(logical(PARTITION==notEmpty(n))) = n;
end
npops = length(notEmpty);
%--------------------------------------------------------------------------
function [partitionSummary, added] = addToSummary(logml, partitionSummary, worstIndex)
% Tiedet<EFBFBD><EFBFBD>n, ett<EFBFBD> annettu logml on isompi kuin huonoin arvo
% partitionSummary taulukossa. Jos partitionSummary:ss<EFBFBD> ei viel<EFBFBD> ole
% annettua logml arvoa, niin lis<EFBFBD>t<EFBFBD><EFBFBD>n worstIndex:in kohtaan uusi logml ja
% nykyist<EFBFBD> partitiota vastaava nclusters:in arvo. Muutoin ei tehd<EFBFBD> mit<EFBFBD><EFBFBD>n.
global PARTITION;
apu = isempty(find(abs(partitionSummary(:,2)-logml)<1e-5,1));
if apu
% Nyt l<EFBFBD>ydetty partitio ei ole viel<EFBFBD> kirjattuna summaryyn.
npops = length(unique(PARTITION));
partitionSummary(worstIndex,1) = npops;
partitionSummary(worstIndex,2) = logml;
added = 1;
else
added = 0;
end
%--------------------------------------------------------------------------
function [counts, sumcounts] = initialCounts(ind_counts)
global PARTITION;
pops = unique(PARTITION);
npops = max(pops);
counts = zeros(size(ind_counts,1), size(ind_counts,2), npops,'uint16');
sumcounts = zeros(npops, size(ind_counts,2),'uint16');
for i = 1:npops
inds = find(PARTITION == i);
counts(:,:,i) = sum(ind_counts(:,:,inds), 3);
sumcounts(i,:) = sum(counts(:,:,i),1);
end
%--------------------------------------------------------------------------
function logml = computeLogml(adjprior_cq, adjprior_sp)
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
cq_counts = double(CQ_COUNTS);
cq_sumcounts = double(CQ_SUMCOUNTS);
sp_counts = double(SP_COUNTS);
sp_sumcounts = double(SP_SUMCOUNTS);
npops = size(CQ_COUNTS, 3);
cq_logml = sum(sum(sum(gammaln(cq_counts+repmat(adjprior_cq,[1 1 npops]))))) ...
- npops*sum(sum(gammaln(adjprior_cq))) - ...
sum(sum(gammaln(1+cq_sumcounts)));
sp_logml = sum(sum(sum(gammaln(sp_counts+repmat(adjprior_sp,[1 1 npops]))))) ...
- npops*sum(sum(gammaln(adjprior_sp))) - ...
sum(sum(gammaln(1+sp_sumcounts)));
logml = cq_logml - sp_logml;
clear cq_counts cq_sumcounts sp_counts sp_sumcounts;
%--------------------------------------------------------------------------
function popLogml = computePopulationLogml(pops, adjprior_cq, adjprior_sp)
% Palauttaa length(pops)*1 taulukon, jossa on laskettu korikohtaiset
% logml:t koreille, jotka on m<EFBFBD><EFBFBD>ritelty pops-muuttujalla.
global CQ_COUNTS; global CQ_SUMCOUNTS;
global SP_COUNTS; global SP_SUMCOUNTS;
cq_counts = double(CQ_COUNTS);
cq_sumcounts = double(CQ_SUMCOUNTS);
sp_counts = double(SP_COUNTS);
sp_sumcounts = double(SP_SUMCOUNTS);
nall_cq = size(CQ_COUNTS,1);
nall_sp = size(SP_COUNTS, 1);
ncliq = size(CQ_COUNTS,2);
nsep = size(SP_COUNTS, 2);
z = length(pops);
popLogml_cq = ...
squeeze(sum(sum(reshape(...
gammaln(repmat(adjprior_cq,[1 1 z]) + cq_counts(:,:,pops)) ...
,[nall_cq ncliq z]),1),2)) - sum(gammaln(1+cq_sumcounts(pops,:)),2) - ...
sum(sum(gammaln(adjprior_cq)));
popLogml_sp = ...
squeeze(sum(sum(reshape(...
gammaln(repmat(adjprior_sp,[1 1 z]) + sp_counts(:,:,pops)) ...
,[nall_sp nsep z]),1),2)) - sum(gammaln(1+sp_sumcounts(pops,:)),2) - ...
sum(sum(gammaln(adjprior_sp)));
popLogml = popLogml_cq - popLogml_sp;
clear cq_counts cq_sumcounts sp_counts sp_sumcounts;
%-------------------------------------------------------------------
%--------------------------------------------------------------
function newline = takeLine(description,width)
%Returns one line from the description: line ends to the first
%space after width:th mark.
% newLine = description(1:width);
n = width+1;
while ~isspace(description(n)) && n<length(description)
n = n+1;
end;
newline = description(1:n);
function dispLine
disp('---------------------------------------------------');
function dispCancel
disp('** CANCELLED');
function num2 = omaRound(num)
% Py<EFBFBD>rist<EFBFBD><EFBFBD> luvun num 1 desimaalin tarkkuuteen
num = num*10;
num = round(num);
num2 = num/10;
%---------------------------------------------------------
function mjono = logml2String(logml)
% Palauttaa logml:n string-esityksen.
mjono = ' ';
if abs(logml)<10000
%Ei tarvita e-muotoa
mjono(7) = palautaYks(abs(logml),-1);
mjono(6) = '.';
mjono(5) = palautaYks(abs(logml),0);
mjono(4) = palautaYks(abs(logml),1);
mjono(3) = palautaYks(abs(logml),2);
mjono(2) = palautaYks(abs(logml),3);
pointer = 2;
while mjono(pointer)=='0' && pointer<7
mjono(pointer) = ' ';
pointer=pointer+1;
end
if logml<0
mjono(pointer-1) = '-';
end
else
suurinYks = 4;
while abs(logml)/(10^(suurinYks+1)) >= 1
suurinYks = suurinYks+1;
end
if suurinYks<10
mjono(7) = num2str(suurinYks);
mjono(6) = 'e';
mjono(5) = palautaYks(abs(logml),suurinYks-1);
mjono(4) = '.';
mjono(3) = palautaYks(abs(logml),suurinYks);
if logml<0
mjono(2) = '-';
end
elseif suurinYks>=10
mjono(6:7) = num2str(suurinYks);
mjono(5) = 'e';
mjono(4) = palautaYks(abs(logml),suurinYks-1);
mjono(3) = '.';
mjono(2) = palautaYks(abs(logml),suurinYks);
if logml<0
mjono(1) = '-';
end
end
end
function digit = palautaYks(num,yks)
% palauttaa luvun num 10^yks termin kertoimen
% string:in<EFBFBD>
% yks t<EFBFBD>ytyy olla kokonaisluku, joka on
% v<EFBFBD>hint<EFBFBD><EFBFBD>n -1:n suuruinen. Pienemmill<EFBFBD>
% luvuilla tapahtuu jokin py<EFBFBD>ristysvirhe.
if yks>=0
digit = rem(num, 10^(yks+1));
digit = floor(digit/(10^yks));
else
digit = num*10;
digit = floor(rem(digit,10));
end
digit = num2str(digit);
%--------------------------------------------------------------------------
function [emptyPop, pops] = findEmptyPop(npops)
% Palauttaa ensimm<EFBFBD>isen tyhj<EFBFBD>n populaation indeksin. Jos tyhji<EFBFBD>
% populaatioita ei ole, palauttaa -1:n.
global PARTITION;
pops = unique(PARTITION)';
if (length(pops) ==npops)
emptyPop = -1;
else
popDiff = diff([0 pops npops+1]);
emptyPop = min(find(popDiff > 1));
end
%--------------------------------------------------------------------------
Reference in a new issue Copy permalink