% function w = ETAM_train_temporal(ptn, alpha, nextpattern)
% Inputs:
% ptn: a 2-D matrix with size neurons * patterns.
% alpha: the learning rate for each step, default = 0.005.
% nextpattern: 1 * patterns vector indicates the next state transition.
% Outputs:
% w: the weight using for ETAM with size neurons * (neurons + 1), the +1 
%    part is for constant threshold.

function w = ETAM_train_temporal(ptn, alpha, nextpattern)

	[ptn_h, ptn_w] = size(ptn);
	
    outofcube = ptn_h * ptn_h;
    if (nargin < 2)
        alpha = 0.005;
    end
    
    if (nargin < 3)
        nextpattern = 1:ptn_w;
    end
	
    w = ptn(:, nextpattern) * ptn';
    % normalized each row before use.
    for i = 1:ptn_h
        w(i,:) = w(i, :) / norm(w(i,:));
    end

    % the last column is for constant term.
	w = [w zeros(ptn_h, 1)];
	% the last row is for minus constant.
	ptn_1 = [ptn; -ones(1, ptn_w)];
	
	for i = 1:ptn_h
        %Cp = find(ptn(i, :) == 1);
        %Cn = find(ptn(i, :) == -1);
        Cp = find(ptn(i, nextpattern) == 1);
        Cn = find(ptn(i, nextpattern) == -1);
        
        dist = w(i, :) * ptn_1;
        [dp indexp] = min(dist(Cp));
        [dn indexn] = max(dist(Cn));
        
        if isempty(Cp)
            w(i, end) = outofcube;
            continue;
        elseif isempty(Cn)
            w(i, end) = -outofcube;
            continue;
        else
            pre_dm = -outofcube;
            DONE = 0;
            while (((dp-dn) / 2 > pre_dm) || (~DONE))
                DONE = 1;
                % shift the hyperplane to the midline of patterns.
                w(i, end) = w(i, end) + (dp + dn) / 2;
                % rotate the hyperplane to increase the distances.
                w1 = w(i, 1:(end - 1)); % record previous w in case of recovering.
                w(i, 1:(end - 1)) = w(i, 1:(end - 1)) + alpha * (ptn(:, Cp(indexp))' * ptn(i, nextpattern(Cp(indexp))) + ptn(:, Cn(indexn))' * ptn(i, nextpattern(Cn(indexn))));
                                
                % normalize w except constant term.
                w(i, 1:(end - 1)) = w(i, 1:(end - 1)) / norm(w(i, 1:(end - 1)));
                pre_dm = (dp - dn) / 2;
                                
                % trying for pattern recovery.
                sum = w(i, :) * ptn_1;
                
                %sum = sign(sum);
                %sum(find(sum == 0)) = 1;
                sum = (sum >= 0) * 2 - 1;
                
                if any(sum ~= ptn(i,nextpattern))
                %if any(sum ~= ptn(i,:))
                    DONE = 0;
                end
                
                dist = w(i, :) * ptn_1;
                [dp indexp] = min(dist(Cp));
                [dn indexn] = max(dist(Cn)); 
            end
            w(i, 1:(end - 1)) = w1;    % recovery 
        end 
    end
end