Random gradient descent See the variance

function c = solver(a)
% SOLVER:  Takes the amino acid chain, a, as input and returns the
%   corresponding configuration sequence, c.  The configuration sequence, c,
%   consists of the directions to the next successive amino acid:
%   Direction: 1 = east, i = north, -1 = west, -i = south
% Simple gradient descent with random seed at every iteration
% Slow because it spends more than a third of its time scoring, but elegant and very easy to understand: no special cases built in
% By Alex Backer

minimprov=1; % should be cost (in time) of one link loop
looptime=0.03;
maxtime=179; % sec

k1 = 1;
k2 = 0.142857247;
k3 = 0.05;

t0=cputime;
c = [ones(length(a)-1,1)];
currscore=score(c,a,k1,k2,k3,t0);
cv=[-1 i 1 -i];

lena=length(a);

improv=minimprov; % init
            minimprov=k2* exp (k3*looptime)/k1;
while improv>=minimprov & cputime-t0<maxtime-looptime,
    improv=0;
    for link=randperm(lena-1),
        loopstart=clock;
        for newconf=cv(randperm(4)),
            if newconf~=c(link), % use only new configurations to compare to current
                if link==1 | newconf~=-c(link-1), % Do not check config where next aa goes in prev aa's pos
                    cnew=c;
                    cnew(link)=newconf;
                    score=score(cnew,a,k1,k2,k3,t0);
                    if score>-1 & score<currscore,
                        c=cnew;
                        improv=currscore-score; % score decreases over iterations
                        currscore=score;
                    end
                end
            end
        end
        if link==1,
            looptime=cputime-loopstart;
            minimprov=k2* exp (k3*looptime)/k1;
        end
    end
end
            

function results=score(c,a,k1,k2,k3,t0)
    if size(c,2)==1
        p = [0;cumsum(c)];
    else
        p = [0 cumsum(c)];
    end

   % Check the solution
    
    if length(p)~=length(unique(p))
        %('The returned vector intersects itself')
        results=-1;
    else,
    % Score it
    % Calculate the distance matrix
    p = p(logical(a));
    [x,y] = meshgrid(1:length(p));
    dist = abs(p(x) - p(y));

    results = sum(sum(dist))/2;
    %score = k1 * results + k2* exp (k3*etime(clock,t0));

    end