function list_feature_l=PNG_list_features(F,degree)
  list_feature_l=zeros(0,degree);
  list_feature_l=[zeros(F,degree-1) (1:F)'];
  if degree<=1 return; end
  counter=size(list_feature_l,1);
  for f1=1:F
    for f2=1:F
      if f2<f1 continue, end
      counter=counter+1;
      list_feature_l(counter,end-1:end)=[f1 f2];
    end
  end
  if degree<=2 return; end
  for f1=1:F
    for f2=1:F
      for f3=1:F
        if f2<f1 continue, end
        if f3<f2 continue, end
        counter=counter+1;
        list_feature_l(counter,end-2:end)=[f1 f2 f3];
      end
    end
  end
  if degree>3 error('script has to be modified') end
end


function list_feature_l=PNG_list_features2(F,degree)
%%more general function
  list_feature_l=zeros(1,degree);
  for d=1:degree
    new=zeros(0,degree);
    for f=1:F
      for l=1:size(list_feature_l,1)
        ind_non_zero=find(list_feature_l(l,:)~=0,1,'first');
        if isempty(ind_non_zero)
          new=[new;[zeros(1,degree-1) f]];
        else 
          ind_write=ind_non_zero-1;
          if ~(f<=list_feature_l(l,ind_non_zero)) continue, end
          new=[new;[zeros(1,ind_non_zero-2) f list_feature_l(l,ind_non_zero:end)]];
        end
      end
    end
    list_feature_l=[list_feature_l;new];
  end
end


function Xp=PNG_poly(X,degree)
  list_feature_l=PNG_list_features(F,degree);
  for l=1:size(list_feature_l,1)
    Xp(:,l)=1;
    for f_=1:size(list_feature_l)
      if list_feature_l(f_)==0 continue, end
      f=list_feature_l(f_);
      Xp(:,l)=Xp(:,l).*X(:,f);
    end
  end
end

S=load('training.mat'); 
S.X=Xp;
save('training.mat',S);

ind=randperm(size(X,1));
subset=X(ind(1:10),:);

function Xf=PNG_rbf(X,subset,gamma)
  Xf=zeros(size(X,1),size(subset,2));
  for v_=1:length(subset)
    x=X(subset(v_),:); 
    for n=1:size(X,1)
      Xf(n,v_)=kernel(X(n,:),x,gamma);
    end
  end
end

function val=kernel(x,y,gamma)
  val=exp(-gamma*sum((x-y).^2));
end

S=load('training.mat');
info.X=S.X; 
info.Y=S.Y; 
info.k=3;

function y_hat=PNG_predict(info,x)
  tab=zeros(size(info.X,1),1);
  for n=1:size(info.X,1)
    xn=info.X(n,:);
    tab(n)=dist(xn,x);
  end
  [~,order]=sort(tab);
  class_l=S.Y(order(1:k));
  classes=unique(class_l);
  counter=zeros(length(classes),1);
  for c_=1:length(classes)
    counter(c_)=sum(class_l==classes(c_));
  end
  [~,ind_c]=max(counter);
  y_hat=classes(ind_c);
end


function d=dist(xa,xb)
  d=sqrt(sum((x-y).^2));
end


function PNG_set_validation(alpha)
  S=load('training.mat');
  N=size(S.X,1);  
  ind=randperm(N);
  ind_validation=ind(1:floor(alpha*N));
  ind_training=ind(1+floor(alpha*N):end);
  X=S.X(ind_validation,:);
  Y=S.Y(ind_validation,:);
  metadata=S.medata;
  save('validation.mat','X','Y','metadata');
  X=S.X(ind_training,:);
  Y=S.Y(ind_training,:);
  metadata=S.medata;
  save('training.mat','X','Y','metadata');
end


function info=PNG_train1(beta)
  S=load('training.mat');
  d_best=[]; OA_best=-Inf; info_best=[]; 
  for d=1:3
    PNG_set_validation(beta);
    S=load('training.mat'); Xt=PNG_poly(X,degree); S.X=Xp; save('training.mat',S);
    S=load('validation.mat'); Xt=PNG_poly(X,degree); S.X=Xp; save('validation.mat',S);
    info=PNE_train2();
    score=PNC_score1(info,'validation.mat',PNE_predict1);
    if score.overall_accuracy>OA_best 
      d_best=d; 
      OA_best=score.overall_accuracy. 
      info_best=info;
    end
  end
  info=info_best; 
  info.d=d_best;
end