% Script to find very different TCV Equilibria using clustering techniques
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% load a large mat file containing outputs from mdb (TCV matlab database)
load liutcv.mat

ip   = reshape(IP  ,121,[]);
time = reshape(TIME,121,[]);
ft   = reshape(FT  ,121,[]);
ipmax = max(ip); ipmin = min(ip); k = find(abs(ipmin)>abs(ipmax)); ipmax(k) = ipmin(k);
ip0   = ip(1,:);
k = abs(ip) > abs(ipmax*.9);
k = logical(sum(k) - cumsum(k)) & ~isnan(ft) & abs(ip) > 10e3;
% dipdt = (ip([2 2:end],:)-ip([1 1:end-1],:))./(time([2 2:end],:)-time([1 1:end-1],:));
% k = cumsum(dipdt.*sign(ipmax) < -1e5) == 0;
ks = abs(ipmax)<1e6 & abs(ip0)<15e3 & abs(ipmax)>50e3;
k(:,~ks) = false;
k = k(:);
IP   = IP (k  );
IEF  = IEF(k,:);
SHOT = SHOT(k,:);
TIME = TIME(k,:);
n = 100;
IPlim = linspace(-500e3,500e3,n+1);
s = zeros(n,2);
SAMP = [];
for k = 1:length(IPlim)-1
 l = find(IP >= IPlim(k) & IP < IPlim(k+1));
 if length(l) > 1
  y = linkage(IEF(l,:));
  t = cluster(y,'cutoff',250,'criterion','distance');
  s(k,:) = [max(t) length(l)];
  for kt = 1:max(t)
   k = l(t == kt);
   x = IEF(k,:);
   [~,ks] = min(sum((x-mean(x,1)).^2,2));
   SAMP = [SAMP k(ks)];
  end
 end
end
l = SAMP;
SAMP = [];
y = linkage(IEF(l,:));
% t = cluster(y,'cutoff'  , 250,'criterion','distance');
t = cluster(y,'maxclust',5000,'criterion','distance');
for kt = 1:max(t)
 k = l(t == kt);
 x = IEF(k,:);
 [~,ks] = min(sum((x-mean(x,1)).^2,2));
 SAMP = [SAMP k(ks)];
end

x = sortrows([SHOT(SAMP) TIME(SAMP)]);
shot = x(:,1);
time = x(:,2);

%% save as .mat
save liutcvsamp shot time

%% save as .dat
fid = fopen('liutcvselsamp.dat','w');
for ii=1:numel(shot)
    fprintf(fid,'%d  %5.3f\n',shot(ii),time(ii));
end
fclose(fid);