function [bestSurrogate, spectralDiffBest] = siaaft_loop_1d_experimental(fourierCoeff, sortedValues, counterThresshold, iterativeCounterAmplitude, iterativeDenomAmplitude, algorithmVersion)
% INPUT: 
% fourierCoeff:              The 1 dimensional Fourier coefficients that describe the
%                            structure and implicitely pass the size of the matrix.
% sortedValues:              A vector with all the wanted amplitudes (e.g. LWC of LWP
%                            values) sorted in acending order.
% counterThresshold:         The number of unsuccesful iterations before
%                            the algorithm stops.
% iterativeCounterAmplitude: The fraction of amplitude adaptations is give
%                            by iterativeCounterAmplitude/iterativeDenomAmplitude.
%                            For some algorithm versions 2/4 can be a bit
%                            different that 1/2. Please look in the code.
% iterativeDenomAmplitude:   See variable iterativeCounterAmplitude above.
% algorithmVersion:          There are 4 different algorithm versions. 
%                            1: Partially stochastic version (default)
%                            2: Deterministic version
%                            3: Fully stochastic version (where the exact
%                               fraction of amplitude adaptations is
%                               performed).
%                            4: Fully stochastic version (more efficient,
%                               but not exactly the right fraction of
%                               amplitude adaptations.
% 
% OUTPUT:
% bestSurrogate:    The 1D SIAAFT surrogate time series
% spectralDiffBest: The RMS difference between the specified Fourier
%                   spectrum (from the template) and the spectrum of the
%                   surrogate.

% input checking
if (nargin < 3)
    counterThresshold = 1e4; % Number of iterative steps with no approvement before the algorithm stops.
end
if (nargin < 5)
    iterativeCounterAmplitude = 1;
    iterativeDenomAmplitude   = 5;
end
if (nargin < 6)
    algorithmVersion = 1;
end

% Settings
verbose = 0;   % Boolean indicating whether to write to the command window or not.
makePlots = 0; % Best used together during debugging (together with verbose).

% Initialise function
noValues = length(fourierCoeff);
% amplitudeDiff = 1;
spectralDiff = 1;
spectralDiffBest = 100;
standardDeviation = std(sortedValues);
noAdaptedValues = floor(noValues*iterativeCounterAmplitude/iterativeDenomAmplitude);

indices = init_regular_indices(iterativeCounterAmplitude, iterativeDenomAmplitude, noValues);
noIndices = size(indices,2);

% Set limits for plots.
if makePlots
    xmin = 1;
    xmax = noValues;
    ymin = sortedValues(1);
    ymax = sortedValues(end);
end

% The method starts with a randomized uncorrelated surrogate time series
% with the amplitudes (PDF) of sorted_values.
[dummy,index] = sort(rand(size(sortedValues)));
surrogate(index) = sortedValues;

% Main intative loop
counter = 1; % Counts the number of iterations in every stage.
number = 1;  % The index to the vector to use from the matrix indices (for algorithmVersion 1 and 2).
% counterSpec = 1;
% maxCounter = 1e4;

for stageCounter=1:2
    if (stageCounter == 2)
        % In this second stage all amplitudes and coefficients are used
        % in the iteration, just like in the original iaaft algorithm.
        disp('Second stage')
        iterativeCounterAmplitude = iterativeDenomAmplitude;
        noAdaptedValues = noValues;        
        % use the best surrogate as the new one.
        surrogate = bestSurrogate;   % The best surrogate of the first phase is the starting surrogate of the second phase
        spectralDiffBest = 100;      % make the second phase use the old surrogate as its best one up to now.
        counter = 1;
        indices = (1:noValues)';     % use all indices, not just a fraction, in the second stage.
        noIndices = 1;
        number = 1;
        counterThresshold = 100; % In the second round with the IAAFT algorithm, it doesn't help to iterate after convergence has stopped.
    end
    
    bestCounter = 0; % This counter counts the number of iteration, and is set to zero each time a new better fitting surrogate is found.
    while ( bestCounter < counterThresshold )
        % addapt the power spectrum
        x=ifft(surrogate);
        phase = angle(x);
        difference = sqrt( mean( (abs(x)-abs(fourierCoeff)).^2  ) );
        spectralDiff = difference/standardDeviation;
% 	    spectralDiffg(counterSpec) = spectralDiff;
%         counterSpec = counterSpec + 1;

        % If the calculation of the power spectrum reveals that the
        % surrogate time series from the previous iteration with perfect
        % amplitudes had the best power spectrum up to now, save this
        % surrogate and its accuracy and reset bestCounter.
        if ( (spectralDiff < spectralDiffBest) & (counter > 1) )
            fprintf(1, 'SpectralDiffBest: %e, bestCounter:, %d\n', spectralDiff, bestCounter);
            spectralDiffBest  = spectralDiff;
%             amplitudeDiffBest = amplitudeDiff;
            bestSurrogate     = surrogate;
            bestCounter = 0;
%            phaseCounter;
        else
            bestCounter = bestCounter + 1;
        end
           
        x = fourierCoeff .* exp(i*phase);
        surrogate = fft(x);
	
        % Plot surrogate after spectral adaptation.
        if ( verbose ), spectralDiff, end        
        if ( makePlots )
            subplot(3,1,1);
            plot(real(surrogate))
            title('Surrogate after spectal adaptation')
            axis([xmin xmax ymin ymax])        
            drawnow
        end
            
        % Adapt the amplitude distribution
        [dummy, index] = sort(real(surrogate)); % We need only the indices. The first value of index points to the highest value, etc.
 
        switch algorithmVersion
            case 1
                % random selection of noIndices vectors with regularly ordered indices
                number = ceil(rand(1)*noIndices);
                surrogate(index(indices(:,number)))=sortedValues(indices(:,number));	                        
            case 2
                % deterministic selection of noIndices vectors with regularly ordered indices
                number = increment(number, 1, 1, noIndices);
                surrogate(index(indices(:,number)))=sortedValues(indices(:,number));	        
            case 3
                % Fully random indices, with exactly the right fraction of amplitude adaptations.
                if (iterativeCounterAmplitude == iterativeDenomAmplitude)
                    surrogate(index) = sortedValues;
                else
                    [dummy, indices] = sort(rand(size(surrogate)));
                    indices = indices(1:noAdaptedValues);
                    surrogate(index(indices))=sortedValues(indices);
                end
            case 4
                % Fully random indices, more efficient version of algorithmVersion 3.
                if (iterativeCounterAmplitude == iterativeDenomAmplitude)
                    surrogate(index) = sortedValues;
                else                
                    indices = ceil(rand(noAdaptedValues,1)*noValues);
                    surrogate(index(indices))=sortedValues(indices);
                end
        end
        
        % Plot surrogate after amplitude adaptation.        
        if ( verbose ), amplitudeDiff, end       
        if ( makePlots )
            subplot(3,1,2);
            plot(real(surrogate))
            title('Surrogate after amplitude adaptation')
            axis([xmin xmax ymin ymax])        
            drawnow
        end
      
        counter=counter+1;
%        bestCounter;
        drawnow % Allows escaping the while loop with ^c.
	end % while not conveged   
end % for the 2 stages.

surrogate = double(surrogate);