rand('state',100);
n=200;
d = 10;
noiselevel = 0.6;
x = rand(n,d);
for i = 1:(d-2)
   x(:,i) = (x(:,i) + x(:,i+1) + x(:,i+2))/3;
end
ftrue = 4 * x(:,1) /3 + pi * sin(pi * x(:,3)) + 8 * x(:,6).^5 + 2/(2.71828 - 1) * exp(x(:,8)) - 5;
y = ftrue + randn(n,1) * noiselevel;

K3darray = zeros(n,n,d);
for i = 1:d
   K3darray(:,:,i) = newsplinematrix(x(:,i),n);
end

folds = 5;
splitsize = floor(n/folds) * ones(folds, 1);
splitsize(folds) = n - (folds - 1) * floor(n/folds);
p = randperm(n);

theta = ones(d,1);
Kth0 = zeros(n,n);
for i = 1:d
   Kth0 = Kth0 + theta(i) * K3darray(:,:, i);
end

lambda0 = tunela02(Kth0, y, p, folds, splitsize);
lambda0

m = 3 * ceil(log2(8 * d));

bigKth0 = [Kth0 + lambda0 * eye(n), ones(n,1); ones(1,n), 0];

cb0 = bigKth0 \ [y;0];
c0 = cb0(1:n);
b0 = cb0(n+1);
G0 = zeros(n,d);
for inneriter = 1:d
   G0(:, inneriter) = K3darray(:,:,inneriter) * c0;
end

solution = zeros(n+1, d+1, m);
fest = zeros(n, m);
for iter = 1:m
   iter
   M = 2^(iter/3 - 2);
   solution(:,:,iter) = solve1(K3darray, y, d, G0, c0, b0, lambda0, M);
   fest(:,iter) = solution(1:n, 1:d, iter) * (solution(n+1, 1:d, iter))' + solution(n+1, d+1, iter);
   sqrtmse(iter) = norm(fest(:,iter) - ftrue);
end


cv = zeros(m,folds);

for run = 1:folds
   tuneindex = p(((run - 1) * floor(n/folds) + 1) : ((run - 1) * floor(n/folds) + splitsize(run)));
   temp = zeros(n, 1);
   temp(tuneindex) = 1;
   ytune = y(logical(temp));
   ytrain = y(logical(1-temp));

   ntrain = length(ytrain);

   K3dtrain = zeros(ntrain,ntrain,d);
   for i = 1:d
      K3dtrain(:,:,i) = K3darray(logical(1-temp), logical(1-temp), i);
   end

   theta = ones(d,1);
   Kth0 = zeros(ntrain,ntrain);
   for i = 1:d
      Kth0 = Kth0 + theta(i) * K3dtrain(:,:, i);
   end

   bigKth0 = [Kth0 + lambda0 * eye(ntrain), ones(ntrain,1); ones(1,ntrain), 0];

   cb0 = bigKth0 \ [ytrain;0];
   c0 = cb0(1:ntrain);
   b0 = cb0(ntrain+1);
   G0 = zeros(ntrain,d);
   for inneriter = 1:d
      G0(:, inneriter) = K3dtrain(:,:,inneriter) * c0;
   end
   for iter = 1:m
      iter
      M = 2^(iter/3 - 2);
      solutiontrain = solve1(K3dtrain, ytrain, d, G0, c0, b0, lambda0, M);
      Kpred = zeros(n - ntrain, ntrain);
      for inneriter = 1:d
         Kpred = Kpred + solutiontrain(ntrain+1, inneriter) * K3darray(logical(temp), logical(1-temp),inneriter);
      end
      prediction = Kpred * solutiontrain(1:ntrain, d+1) + solutiontrain(ntrain+1, d+1);
      cv(iter, run) = (norm(ytune - prediction))^2;
   end
end

totalcv = cv * ones(folds, 1);
[c, bestiter] = min(sqrtmse);
[C, cviter] = min(totalcv);
cviter
bestiter

solution(n+1, 1:d, cviter)
solution(n+1, 1:d, bestiter)

festcomponents = solution(1:n, 1:d, cviter) * diag(solution(n+1, 1:d, cviter));
fbestcomponents = solution(1:n, 1:d, bestiter) * diag(solution(n+1, 1:d, bestiter));

subplot (4,4,9), scatter(x(:,1),festcomponents(:,1))
subplot (4,4,10), scatter(x(:,3),festcomponents(:,3))
subplot (4,4,11), scatter(x(:,6),festcomponents(:,6))
subplot (4,4,12), scatter(x(:,8),festcomponents(:,8))

subplot (4,4,13), scatter(x(:,1),fbestcomponents(:,1))
subplot (4,4,14), scatter(x(:,3),fbestcomponents(:,3))
subplot (4,4,15), scatter(x(:,6),fbestcomponents(:,6))
subplot (4,4,16), scatter(x(:,8),fbestcomponents(:,8))

