function out=hedgeMakeCloud2(data,xs)
%hedgeMakeCloud2 takes historical cloud data and generates a random
%sequence drawn from the same PDF. This is then shuffled until
%its autocorrelation is sufficiently similar to the autocorrelation of the 
%original sequence.
%=========================================================================
% Usage:
% out=hedgeMakeCloud2(data,xs)
% Inputs:
%       * data=vector of cloud cover fractions
%       * xs = bins for the PDF, in NOAA cloud cover units (i.e. 0:10)
% Output, out = new random cloud cover vector in NOAA cloud cover units

%length of cloud cover vector
numDays=length(data);

%calculate target autocorrelation
xcorr1=xcorr(data,data);

% get PDF of values in cloud data
ys=hist(data,xs);

%populate a vector with random values from that PDF
dataTemp=zeros(numDays,1);
for i=1:numDays
    dataTemp(i)=randarb(xs,ys);
end

% calculate autocorrelation of random vector as a baseline
score=xcorr(dataTemp,dataTemp);

%monte carlo search parameters:
scoreTemp=sqrt(sum((score-xcorr1).^2)); %current dist between autocorrelations
bestScore=scoreTemp;
closeEnough=1000;                       %what distance is good enough - this is empirical, for seasons of length 214 days.

%initialize tick counter for monte carlo
i=0;

%number of sequential days to move around at once
windowSize=2;

%shuffle until the score is good enough.
while scoreTemp>closeEnough
    %find two random points in the vector to swap.
    which1=ceil(rand()*(numDays-windowSize));
    which2=ceil(rand()*(numDays-windowSize));
    %switch the data (and subsequent window) at these points
    dataProposed=dataTemp;
    dataProposed(which1:which1+windowSize)=dataTemp(which2:which2+windowSize);
    dataProposed(which2:which2+windowSize)=dataTemp(which1:which1+windowSize);
    %update counter
    i=i+1;
    %get new distance score
    xcorrTemp=xcorr(dataProposed,dataProposed);
    scoreTemp=sqrt(sum((xcorrTemp-xcorr1).^2));
    % if it's an improvement, keep it, and update monte carlo parameters
    if scoreTemp < bestScore
        bestScore=scoreTemp;
%         disp(corr2(xcorrTemp,xcorr1));  %comment this line in to see correlation between autocorrelation vectors
        dataTemp=dataProposed;
    end
end

%output the final sequence
out=dataTemp;