computeEntropy.m

function [iOptRec,time] = computeEntropy(M,Ts,iTRecs,xyRec,show,PreiOptRec)
%COMPUTEENTROPY Compute a full entropy optimal design
%   This function computes the optimal receiver design for a given model
%   and number of receivers using the entropy as a quality measure.
%   Inputs:
%   M - model parameters
%   Ts - Time data as generated by 'compute_traveltimes3D'
%   iTRecs - index of receivers in Ts
%   xyRec - X and Y coordinates of receivers
%   show - boolean to toggle graph outputs

% initialize vectors
xOptRec = nan(1,size(iTRecs,1));
iOptRec = zeros(1, M.nRecMax);
useRec = false(1,size(iTRecs,1));
usedT = [];

% If receivers are precomputed, use those.
if nargin > 5
    iOptRec(1:length(PreiOptRec)) = PreiOptRec;
    nRecMin = length(PreiOptRec)+1;
    for i = 1:length(PreiOptRec)
        usedT = [usedT reshape(Ts(iTRecs(:,iOptRec(i)), :),[],1)];
    end
else
    nRecMin = 1;
end

% initialize figure for showing criterion map
if show
    p = numSubplots(M.nRecMax);
    figure(2)
    clf
end

% Loop until all nRecMax receivers are placed in a sequential fashion
tBegin = tic;
for i = nRecMin:M.nRecMax
    
    tic
    entropies = nan(size(xyRec,2),1);
    
    % Loop over all possible receiver locations and calculate the entropy
    parfor it = 1:size(xyRec,2)
        
        % Time points to use
        A = [reshape(Ts(iTRecs(:,it), :),[],1) usedT];
        
        % Calculate entropy
        entropies(it) = kdpee(A);
        
    end
    toc
    
    % Find optimal receiver
    [maxEntropy, iOptRec(i)] = max(entropies);
    xOptRec(i) = xyRec(iOptRec(i));
    useRec(iOptRec(i)) = true;
    
    % Plot the map of calculated entropy values
    if show
        plotMetric(entropies,xyRec,iOptRec,p,i,'ENT')
    end
    
    % Store the data points generated by selected receivers
    usedT = [usedT reshape(Ts(iTRecs(:,iOptRec(i)), :),[],1)];
    
    % Output status
    fprintf('Entropy in iteration %d is %0.4f \r',i,maxEntropy)
    fprintf('%d%% done. \r',round(i/M.nRecMax*100))
    
end
time = toc(tBegin);

% Output final status
fprintf('Optimal design is calculated. \r')
fprintf('Assigned %d unique receiver locations. \r', sum(useRec))
end