function [ComputerState,logdata] = qsearch( varargin)
% QSEARCH Perform a quantum search using Grover's algorithm.
%
% 
% 
% qsearch(n, M, ['opts']) performs a quantum search on an n-qubit database 
% looking for M randomly chosen solutions.  n must be an integer, and M 
% must be a scalar.  
%
% qsearch(n, ['opts']) is equivilant to qsearch(n, 1, ['opts']).
%
% qsearch(n, ANSWERS, ['opts'] ), with ANSWERS a vector of length > 1, 
% performs a quantum search on an n-qubit database looking for the values
% specified in answers.  This will not work if you want to search for only
% one pre-specified value. 
%
% qsearch(n, 1, ANSWER, ['opts'] ), with answer a scalar will perform a 
% quantum search on an n-qubit database looking for the value specified in 
% answer. M will be set to 1.
%
%
% 'opts' is an optional string that may be used to control the 
% methods used to conduct the search. 'opts' may contain up to five 
% characters, in any order, one for each of the five parameters that may be
% adjusted: matrix method, Hadamard method, logging method, warnings, and the number
% of iterations.  'opts' is not case sensitive.
% 
% Matrix method controls the type of matrix used to create the oracle and
% phase-shift matrixes.  Characters controlling the matrix method are:
%
% 's'  Use sparse double-precision matrices (the default). This uses 16*2^n
% bytes of memory to store the phase-shift and oracle matrices.
% 'd'  Use double-precision, full matrices. This will require 16*2^(2n)
% bytes of memory to store the phase-shift and oracle matrices. 
% 'i'  Use single-precision, full matrices. This will require 8*2^(2n)
% bytes of memory to store the phase-shift and oracle matrices. 
% 
% NOTE: MATLAB does not support sparse, single-precision matrices. 
%
%
% Hadamard method controls the method used to compute the hadamard
% transform. This can have a significant effect on runtime, since the
% Hadamard transform is used twice in each iteration. Characters
% controlling the Hadamard method are:
%
% 'p'  Pre-compute the necessary additions and subtractions to apply the 
% Hadamard transform on n qubits, and then apply those additions and 
% subtractions using the function prehadamard.m.  For large numbers of 
% qubits, this provides a very large speedup over the other two methods. 
% This is the default.
% 'f'  Use the built-in function FWHT from the signal processing toolkit to
% compute the Hadamard transform.
% 'm'  Use a 2^n by 2^n matrix to compute the Hadamard transform. This is
% generally is faster than using the built-in function, but slower than the 
% precomputed function.
%
%
% Logging method controls whether and how the intermediate states of the
% computer are stored.  Characters controlling the output method are:
%
% 'r'  At each iteration, store the amplitudes of every computational basis
% state as a row in the matrix OUTPUTS. 
% 'k'  At each iteration, write the amplitudes of every computational basis
% state as a line in the tab-separated file 'iterations.txt'.  This will
% significantly increase runtime if file I/O to the location of
% 'iterations.txt' is slow (i.e. the file is being transferred over a 
% network).
% 'n'  Do not log the amplitudes of each computational basis state. Only
% output the final state after running through itmax iterations.  This is
% the default.
%
% By default, the program prints to stderr the parameters for a search
% before executing that search, and issues several warnings designed to
% avoid accidential out-of memory errors and their effect on the
% performance of your entire computer. To disable this behavior, include
% the character 'u' in the control string.
%
% The number of iterations controls how many times the program executes the
% Grover operator. By default, the program calculates the number of
% iterations according to the following equation:
% theta=asin(2*sqrt(M*(N-M))/N);
% itmax=round(acos(sqrt(M/N))/theta);
%
% If any character is specified for the number of iterations, an additional
% integer argument ITMAXADJ must be given after the 'opts' string.
% Characters controlling the number of iterations are:
%
% 't'  Ignore the formula for computing itmax, and set itmax to the value
% of itmaxadj.  This can be useful in combination with some of the output
% methods to see the behavior of the Grover operator over a large number of
% iterations.
% 'a'  Adjust itmax by itmaxadj.  This sets itmax=itmax+itmaxadj. This can
% be useful in combination with some of the output methods to see the
% behavior of the Grover operator in the vicinity of the calculated
% solution.



% start by reading the input arguments, and handling them appropriately.


if nargin<1
    error('Not enough input arguments.')
elseif ischar(varargin{1})
    error('n must be an integer.')
elseif length(varargin{1})>1
    error('n must be a scalar')
end


n=varargin{1};
N=pow2(n);

% ---------------------------------------------------------------------
% Figure out whether the function we're using for the search has been
% specified on the command line, and if not, generate it.
% ---------------------------------------------------------------------

% if second argument not specified, or is an option string, assume we're 
% searching for only one solution
if nargin==1 || ischar(varargin{2})
    M=1;
    optstart=2;
    
    % if second argument is a vector, interpret it as the list of answers
elseif length(varargin{2})>1
    answers=varargin{2};
    optstart=3;
    
    % if second argument is 1 and third argument is a scalar, 
    % third argument is the answer to look for.
elseif nargin>2 && varargin{2}==1 && length(varargin{3})==1 ...
         && ~ischar(varargin{3})
    answers=varargin{3};
    optstart=4;
    
    % If we've reached this point, there are at least two arguments, the
    % first two arguments are both integers, and the second argument is
    % a scalar other than one.  Interpret the second argument as M.
else 
    M=varargin{2};
    optstart=3;
end

    % if we don't know what the answers are yet, we need to generate them.
if ~exist('answers','var')
    rng('shuffle');
    i=1;
    while i<=M
        answers(1,i)=randi(N,1)-1;
        % make sure we've added a new solution. If the new solution was already
        % there, repeat the loop with the same value of i to pick another.
        if size(answers)==size(unique(answers))
            i=i+1;
        end
    end
else
    answers=unique(answers);
    M=length(answers);
end

% sanity checks
if M>N 
    error('You''re asking me to find more than N solutions in N entries!')
end
if max(answers)>N
    error('You''re asking me to find a value that is too large!')
end

% determime the number of times to run the grover operator

theta=asin(2*sqrt(M*(N-M))/N);
itmax=round(acos(sqrt(M/N))/theta);

% -------------------------------------------------------------------
% Parse options string, and set appropriate options
% -------------------------------------------------------------------

% check if any options were set
if nargin>=optstart &&ischar(varargin{optstart})
    
    % check for the matrix method. Sparse is the default.
    if max(varargin{optstart}=='D') || max(varargin{optstart}=='d')
        matrixmethod='d';
    elseif max(varargin{optstart}=='I') || max(varargin{optstart}=='i')
        matrixmethod='i';
    else
        matrixmethod='s';
    end
    
    % check for the hadamard method. Precomputed is the default.
    if max(varargin{optstart}=='F') || max(varargin{optstart}=='f')
        hadamardmethod='f';
    elseif max(varargin{optstart}=='M') || max(varargin{optstart}=='m')
        hadamardmethod='m';
    else
        hadamardmethod='p';
    end
        
    % check for the loggingmethod
    if max(varargin{optstart}=='R') || max(varargin{optstart}=='r')
        loggingmethod='r';
    elseif max(varargin{optstart}=='K') || max(varargin{optstart}=='k')
        loggingmethod='k';
    else 
        loggingmethod='n';
    end    
    
    % should we allow pausing?
    if max(varargin{optstart}=='U') || max(varargin{optstart}=='u')
        pause off;
        warn=false;
    else
        pause on;
        warn=true;
    end
    
    % check to see if we need to adjust itmax
    if nargin>optstart && ~ischar(varargin{optstart+1}) ...
        && length(varargin{optstart+1})==1
            if max(varargin{optstart}=='t') || max(varargin{optstart}=='T')
                itmax=varargin{optstart+1};
            elseif max(varargin{optstart}=='a') || max(varargin{optstart}=='A')
                itmax=itmax+varargin{optstart+1};
            end
    end
else % if no options were set, use the defaults.
    matrixmethod='s';
    loggingmethod='n';
    hadamardmethod='p';
    warn=true;
    pause on;
    
end
        

       

    

% 

% 
% % switch loggingmethod
% %     case {'d','disk','disk space'}
% %         loggingmethod='disk space';
% %     case {'m','memory','mem','r','ram','RAM'}
% %         loggingmethod='memory';
% %     case {'n','none',[]}
% %         loggingmethod='none';
% %     otherwise
% %         error('I don''t know that output method. \nValid options are ''n'', ''none'', ''m'', ''memory'',''d'', ''disk''.')
% % end
% 



% Make a lookup table for translating the single-character options to
% human-readable ones. 

trans{'s'}='Sparse';
trans{'d'}='full, Double';
trans{'i'}='full, sIngle';

trans{'p'}='Pre-computed function';
trans{'f'}='built-in Function';
trans{'m'}='Matrix';

trans{'r'}='memoRy';
trans{'k'}='disK';
trans{'n'}='None';

trans{'u'}='do not paUse';

if warn
    % print to stderr
    fprintf(2,['Running quantum search using the following parameters:\n' ...
        'n=%.0f\t\tmatrixmethod=%s\nM=%.0f\t\thadamardmethod=%s\nitmax=%.0f\toutputmethod=%s\n ']...
        ,n,trans{matrixmethod},M,trans{hadamardmethod},itmax,trans{loggingmethod});

    % estimate storage needs, to avoid out of memory errors

    % RAM used by oracle and phase matrices
    switch matrixmethod
        case 's'
            ramused=2*24*N;
        case 'i'
            ramused=2*4*N^2;
        case 'd'
            ramused=2*8*N^2;
    end

    % RAM used by the hadamard matrix
    switch hadamardmethod
        case 'm'
            ramused=ramused+8*N^2;
    end
    
    % RAM used to store the computer state
    ramused=ramused+8*N;

    % RAM/disk space used to store the output data
    switch loggingmethod
        case 'r'
            ramused=ramused+8*N*itmax;
        case 'k'
            diskused=10*N*itmax;
            if diskused>pow2(500,20)
            warning(['This will require %g megabytes of disk space to store the output\n' ...
                'Are you sure you want to continue?\n' ...
                'Press enter to continue, or Ctrl+C to abort.'],pow2(diskused,-20))
                pause
            end
    end


    if ramused>pow2(2,30) && warn
        warning(['This will require %g gigabytes of memory\n' ...
            'Are you sure you want to continue?\n' ...
            'Press enter to continue, or Ctrl+C to abort.'],pow2(ramused,-30))
        pause
    end
        
% if n>14 && ~isequal(loggingmethod(1),'n') && warn
%     warning(['You''ve asked me to do a computation on %.0f qubits. \n' ...
%         'This will require more than %g gigabytes of %c just to store' ...
%         'the output and sparse matrices, and on the order of %g operations. \n \n' ...
%         'This is likely to take a long time.\n' ...
%         'Are you sure you want to continue?\n' ... 
%         'Press enter to continue, or Ctrl+C to abort.'],n,pow2(outsize,-30),loggingmethod,numcalcs)
%     pause
% end

end
% ----------------------------------------------------------------
% Beginning of the actual program
% ----------------------------------------------------------------



switch hadamardmethod
% If we're using matrices, generate a hadamard matrix of appropriate type
    case 'm'
        switch matrixmethod
            case 'i'
                H=single(hadamard(N)/sqrt(N));
            otherwise
                H=hadamard(N)/sqrt(N);
        end
% If using the custom hadamard function, compute which additions and 
% subtractions to do ahead of time, so we only have to do it once.      
    case 'p'
       [index1,index2]=hadsorts(n);
end


% Set up where we're going to store output, if we are going to store output.
switch loggingmethod
    case 'r'
        % pre-allocate memory for the output array. For large values of itmax or N,
        % this is significantly faster than building the array one row at a time.
        logdata(1:itmax,1:N)=0;
    case 'k'
        logfile=fopen('iterations.txt','w');
        logdata=[];
    case 'n'
        logdata=[];
    otherwise
        error('unknown loggingmethod');
end
% Set up the oracle and phase-shift matrices
switch matrixmethod
    case 'i'
        OracleMatrix=generate_single_oracle(n,answers);
        PhaseMatrix=generate_single_phase(n);
    case 's'
        OracleMatrix=generate_sparse_oracle(n,answers);
        PhaseMatrix=generate_sparse_phase(n);
    case 'd'
        OracleMatrix=generate_full_oracle(n,answers);
        PhaseMatrix=generate_full_phase(n);
    otherwise
        error('unknown matrixmethod');
end


% ComputerState is the current state vector of the computer
% initialize the ComputerState to |0>
ComputerState(1:N)=0;
ComputerState(1)=1;
% intial state has computer at all zeros. Amplitude of the all-zero state
% is 1
% oracle qubut is initialized to (|0>-|1>)\sqrt(2), and remains unchanged
% throughout the computation, so we don't need to simulate it.

% Do initial Hadamard transformation to create superposition
switch hadamardmethod
    case 'm'
        ComputerState=ComputerState*H;
    case 'p'
        ComputerState=prehadamard(ComputerState,index1,index2);
    case 'f'
        ComputerState=fwht(ComputerState,[],'hadamard');
end

% set up logging, if applicable
switch loggingmethod
    case 'r'
        logdata(1,:)=ComputerState;
    case 'k'
        for w=1:N  
            fprintf(logfile,'%d\t',ComputerState(w));
            fprintf(logfile,'\n');
        end
end

% do the computations
for i=1:itmax
    ComputerState=ComputerState*OracleMatrix;
    switch hadamardmethod
        case 'p' % use precomputed func for hadamard gates
            ComputerState=prehadamard(ComputerState,index1,index2);
            ComputerState=ComputerState*PhaseMatrix;
            ComputerState=prehadamard(ComputerState,index1,index2);
        case 'f' % use built-in func for hadamard gates
            ComputerState=fwht(ComputerState,[],'hadamard');
            ComputerState=ComputerState*PhaseMatrix;
            ComputerState=fwht(ComputerState,[],'hadamard');
        case 'm' % use matrices for hadamard gates
            ComputerState=ComputerState*H;
            ComputerState=ComputerState*PhaseMatrix;
            ComputerState=ComputerState*H;
    end
    switch loggingmethod
        case 'r' % log intermediate states to memory
            logdata(i+1,:)=ComputerState;
        case 'k' % log intermediate states to file
            for w=1:N  
                fprintf(logfile,'%d\t',ComputerState(w));
                fprintf(logfile,'\n');
            end
    end
end

            

