function [TP,FP,AP,NsLvl]=testAlgoritmu(data,algoritmus,plotting,handles)

if (~exist('handles', 'var'))
    handles = createDefaultHandles();
end

% TP = true positives
% FP = false positives
% AP = actual positives in signal(numel data.corrSpikes)
% NsLvl = estimated noise level
%
% structure "data" must contain variables "signal" and "corrSpikes",
% where corrSpikes are numbers of samples where spikes should be

%spikeLocations; % numbers of samples, where are spikes detected

%% detect spikes
% Output of this stage = spikeLocations

switch algoritmus
    	case 'WaveclusDefault'
            handles=createDefaultHandles();
            
            [spikesTimestamps, ~, ~, ~]=spdetAmpDetect(data.signal, handles);
            fs = handles.spdet.sr;
            spikeLocations = spikesTimestamps*(fs/1e3);
            
        case 'WaveclusBest'
            handles.spdet.stdmin = 3.96;
            handles.spdet.stdmax = 30;
            handles.spdet.detect_fmin = 600;
            handles.spdet.detect_fmax = 2800; 
            handles.spdet.detection = 'both';
            
            [spikesTimestamps, ~, ~, ~]=spdetAmpDetect(data.signal, handles);
            fs = handles.spdet.sr;
            spikeLocations = spikesTimestamps*(fs/1e3);
            
        case 'NEOdetect'
            handles=createDefaultHandles();
            handles.spdet.stdmin = 4;
            handles.spdet.stdmax = 30;
            handles.spdet.detect_fmin = 770;
            handles.spdet.detect_fmax = 2130; 
            
            spikeLocations=NEOdetect(data.signal, handles);
            %fs = handles.spdet.sr;
            %spikeLocations = spikesTimestamps*(fs/1e3);
            
        case 'OsortWavelet'
            params=[];
            params.bandPass=[300 3000];
            params.nrNoiseTraces=0;
            params.prewhiten=0;
            params.samplingFreq=25000;
            params.limit=100;
            params.extractionThreshold = 0;
            params.alignMethod=3; %1 pos, 2 neg, 3 mixed
            params.peakAlignMethod=4; %1 findPeak, 2 none, 3 power, 4 mteo
            params.detectionMethod=5;
            dp3=[];
            dp3.scalesRange=[0.2 1.0];
            dp3.waveletName='bior1.5';
            params.detectionParams=dp3;
            
            [~, ~,~, spikeTimestamps, ~,~,~] = detectArtificialSpikes( data.signal, params );
            spikeLocations = spikeTimestamps;
            
        case 'OsortPower'
            params=[];
            params.bandPass=[300 3000];
            params.nrNoiseTraces=0;
            params.prewhiten=0;
            params.samplingFreq=25000;
            params.limit=100;
            params.extractionThreshold = 4;
            params.alignMethod=3; %1 pos, 2 neg, 3 mixed
            params.peakAlignMethod=1; %1 findPeak, 2 none, 3 power, 4 mteo
            params.detectionMethod=1; %1 power, 5 WDM
            dp1=[];
            dp1.kernelSize=18;
            params.detectionParams=dp1;
            
            [~, ~,~, spikeTimestamps, ~,~,~] = detectArtificialSpikes( data.signal, params );
            spikeLocations = spikeTimestamps;
            
    	case 'muj'          
            spikeLocations=najdiSpiky(data.signal,handles);
        case 'fast'   
            handles.fast = 1;
            spikeLocations=findSpikes(data.signal,handles);
        case 'prec'   
            handles.fast = 0;
            spikeLocations=findSpikes(data.signal,handles);
    otherwise
        error('Unknown algorithm name!')
end;

%% compute result
corrSpikesArray = data.corrSpikes;
falsePositives = zeros(numel(spikeLocations));

AP = numel(data.corrSpikes);
TP = 0;
FP = 0;

% true positives
for i = 1:numel(spikeLocations)
    found = find(abs((spikeLocations(i) - corrSpikesArray))<50);
    if found>0
        TP = TP+1;
        corrSpikesArray(found(1)) = -5;
    else 
        FP = FP+1;
        falsePositives(FP) = spikeLocations(i);
    end
end
undetectedSpikes = corrSpikesArray(corrSpikesArray>0);
falsePositives = falsePositives(falsePositives>0);

% Compute signal-to-noise ratio
handles=createDefaultHandles();
NsLvl = estimateNoiseLevel(data.signal, data.corrSpikes, handles);


%% plot result
if plotting>0
height = max(data.signal);
figure;
time=(1:numel(data.signal));
plot(time,data.signal);
hold on
% plot all detected spikes
scatter(spikeLocations,ones(size(spikeLocations))*height,'g');
% plot actual spikes (i.e. the ground truth)
scatter(data.corrSpikes,ones(size(data.corrSpikes))*height,'kx');
% plot undetected spikes
scatter(undetectedSpikes,ones(size(undetectedSpikes))*height*1.025,'r');
% plot false detected spikes
scatter(falsePositives,ones(size(falsePositives))*height*1.05,'b');
hold off;
legend('signal','detected','actual','undetected','false positive')
grid on
end
