clear all;clc;  
time = cputime;

%% parameters
n=800;
m=100;
windows = 6;
K = 50;
w = 1;

%% preparing
input = importdata( 'train_20171215.txt' );
dd = input.data;
date = unique(dd(:,1));
data = zeros(length(date),3);
for i=1:length(date)
    temp = find( dd(:,1)==i );
    data(i,1) = i;
    data(i,3) = sum( dd(temp,4) );
    data(i,2) = mean( dd(temp,2) );
end

%% step1
for i=1:n
    train(i,1:windows-1) = data(2+i:i+windows,3) ;
    train(i,windows) = ( data(i+windows+1,3)-data(i+windows,3) );%/data(i+windows,3);
end
for i=1:m
    verification(i,1:windows-1) = data(2+n+i:n+i+windows,3);
    verification(i,windows) = ( data(i+windows+n+1,3)-data(n+i+windows,3) );%/data(n+i+windows,3);
end
[n,p] = size(train);
X = train(:,1:end-1);
label = train(:,end);
[m,p] = size(verification);
Xv = verification(:,1:end-1);
real = verification(:,end);

load RR.mat;
ss =1;
for K=20:10:90
    for kk=1:10
%% step2: k-means++ clustering / Recurrent Phase
[ C,index ] = ClusteringOP( X,K,w );
[K,p] = size(C);

%% step3: metric matrices
for i=1:K
    D(:,i) = vSOPD(C(i,:),X,w);
    Dv(:,i) = vSOPD(C(i,:),Xv,w);
end
temp = max(max(D));
D = D/temp;
Dv = Dv/temp;
save predata D Dv label real C;

%% step4: Harmonic Phase
d = 1;
x_input_domain = [0,1;0,1e+06];
Ehaha = @(gamma)OA(gamma,1);
%[gamma, fx, X_sample, F_sample, Result] = IMGPO_default_run ( Ehaha, d, x_input_domain,200,1,1,1,1);
[gamma, fx] = IMGPO_default_run ( Ehaha, d, x_input_domain,100,1,1,1,1);

%% step5: Prediction
load predata.mat;
[ alpha,result(1,:),estimate ] = OptimizationAlpha( D,label,gamma );  %·µ»¹²âÆÀ
[ result(2,:),prediction ] = forecasting( Dv, gamma, alpha, real );   %Ô¤²â
rr(kk,:) = result(2,:);
    end
    R1(ss,:) = mean( rr );
    R2(ss,:) = var( rr );
    R3(ss,:) = max( rr );
    R4(ss,:) = min( rr );
    ss=ss+1;
end
time = cputime - time

%save K50 train verification result gamma K C windows;
%save RR R1 R2 R3 R4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ dist ]=vSOPD( X,Y,w ) 

% input an n-dimensional row vector X and an k*m-dimensional matrix Y
% w denotes the warping window
% output the SOPD results, which is an k-dimensional column vector

if w == 1
    dist = pdist2( X,Y )';
else
    n=length( X );  
    [k,m]=size( Y );   
    XX=repmat( X,k,1);
    D=10000*ones(k,n,m);                              %D=zeros(k,n,m);  
    D(:,1,1)=abs ( XX(:,1)-Y(:,1) );                     % starting point
    for j=2:w
        D(:,1,j)=abs( XX(:,1)-Y(:,j) ) + D(:,1,j-1);    % Y axis 
    end
    for i=2:w
        D(:,i,1)=abs( XX(:,i)-Y(:,1) ) + D(:,i-1,1);    % X axis
    end
    for i=2:n
        for j=max(2,i-w):1:min(m,i+w)
            tmp(1,:)=D(:,i-1,j-1);
            tmp(2,:)=D(:,i,j-1);
            tmp(3,:)=D(:,i-1,j);
            D(:,i,j)=abs( XX(:,i)-Y(:,j) ) + min(tmp)';
        end
    end
    dist=D(:,n,m);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ dist ]=wSOPD( X,Y,w ) 

% input an n-dimensional row vector X and an m-dimensional row vector Y
% w denotes the warping window
% output the corresponding SOPD result

if w == 1
    dist = pdist2(X,Y);
else
    n=length( X );
    m=length( Y );
    D=100000*ones(n,m);                       % D=zeros(n,m);  
    D(1,1)=abs ( X(1)-Y(1) );                     %  starting point
    for j=2:w
        D(1,j)=abs( X(1)-Y(j) ) + D(1,j-1);    % Y axis  
    end
    for i=2:w
        D(i,1)=abs( X(i)-Y(1) ) + D(i-1,1);    % X axis
    end
    for i=2:n
        for j=max(2,i-w):1:min(m,i+w)
            tmp(1,:)=D(i-1,j-1);
            tmp(2,:)=D(i,j-1);
            tmp(3,:)=D(i-1,j);
            D(i,j)=abs( X(i)-Y(j) ) + min(tmp)';
        end
    end
    dist=D(n,m);
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ C,index ] = ClusteringOP( X,K,w )
% order-subpermutation clustering
Cstart = kstart( X,K,w );                      % kmeans++ initialization
Cre = krestart( X,Cstart,w );               % re-initialization
[ C,index ] = kclustering( X,Cre,w );   % kmeans++ clustering
end

function Cstart = kstart(X, K, w)
[n, p] = size(X);
S = RandStream.getGlobalStream;
% Select the first seed by sampling uniformly at random
index = zeros(1,K);
[C(:,1),index(1)] = datasample(S,[1:n],1,2);
minDist = inf(n,1);
% Select the rest of the seeds by a probabilistic model
for ii = 2:K
    for j=1:n
        mini(j,:) = min( wSOPD( X(j,:),X(index(1:ii-1),:),w ) );
    end
    minDist = min( minDist,mini );
    denominator = sum(minDist);
    sampleProbability = minDist/denominator;
    [C(:,ii), index(ii)] = datasample(S,[1:n],1,2,'Replace',false,...
         'Weights',sampleProbability);        
end
Cstart = C;
end

function [Cre] = krestart( X,Cstart,w )
[n,p] = size(X);
K = length(Cstart);
jndex = setdiff( [1:n],Cstart );
C = X(Cstart,:);
for i=1:length(jndex)
    d = vSOPD( X(jndex(i),:),C,w );
    up = find( d==min(d) );
    if length(up) > 1
        up = up(1);
    end 
    C(up,:) = mean( [C(up,:); X(jndex(i),:)] ); 
end
Cre=C;
end

function [ C,index ] = kclustering( X,Cre,w )
[n,p]=size(X);
[K,p]=size(Cre);
index = cell(K,1);
for i=1:n
    d = vSOPD( X(i,:),Cre,w );
    up = find( d==min(d) );
    index{up(1),1}(end+1) = i;
end
for i=1:K
    cc = index{i,1};
    C(i,:) = mean( X(cc,:) );
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ result,prediction ] = forecasting( Dv, gamma, alpha, real )
    [m,K]=size(Dv);
    phi = zeros(m,K);  
    for i=1:m
        temp = find( Dv(i,:) < gamma );
        phi(i,temp) = exp( -Dv(i,temp) );
    end
    prediction = phi * alpha;
    result(1) = sum( ( prediction - real ).^2 ) / m;
    result(2) = length( find( prediction.*real >0 ) ) / length(real);
    temp1 = find(real>0);
    result(3) = length( find( prediction(temp1)>0 ) )/length(temp1); 
    temp2 = find(real<0);
    result(4) = length( find( prediction(temp2)<0 ) )/length(temp2);
    result(5) = length(temp1)/length(temp2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ result ] = OA(gamma,oa) 
% MSE is the target of this optimization function
load predata.mat;
[n,K] = size(D);
phi = newsample( D,gamma );
alpha = (phi'*phi)^-1 * phi' * label;
estimate = phi * alpha;
EMMSE = sum( (  estimate - label ).^2 )/n;
accuracy = length (  estimate.* label >0 )/ length(label);
if oa==1
    result = EMMSE;
else
    result = accuracy;
end
end








