function varargout = predictability(varargin)
% PREDICTABILITY MATLAB code for predictability.fig
%      PREDICTABILITY, by itself, creates a new PREDICTABILITY or raises the existing
%      singleton*.
%
%      H = PREDICTABILITY returns the handle to a new PREDICTABILITY or the handle to
%      the existing singleton*.
%
%      PREDICTABILITY('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in PREDICTABILITY.M with the given input arguments.
%
%      PREDICTABILITY('Property','Value',...) creates a new PREDICTABILITY or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before predictability_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to predictability_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help predictability

% Last Modified by GUIDE v2.5 02-Jul-2012 18:46:00

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @predictability_OpeningFcn, ...
                   'gui_OutputFcn',  @predictability_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before predictability is made visible.
function predictability_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to predictability (see VARARGIN)

% Choose default command line output for predictability
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes predictability wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = predictability_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

%-------------------- Standard / Empty callbacks -------------------

function FileMenu_Callback(hObject, eventdata, handles)

function OpenMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to OpenMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
file = uigetfile('*.fig');
if ~isequal(file, 0)
    open(file);
end
function PrintMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to PrintMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
printdlg(handles.figure1)

function CloseMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to CloseMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
selection = questdlg(['Close ' get(handles.figure1,'Name') '?'],...
                     ['Close ' get(handles.figure1,'Name') '...'],...
                     'Yes','No','Yes');
if strcmp(selection,'No')
    return;
end

delete(handles.figure1)

function editGamma1_Callback(hObject, eventdata, handles)
function editGamma2_Callback(hObject, eventdata, handles)
function editGamma3_Callback(hObject, eventdata, handles)
function editGamma4_Callback(hObject, eventdata, handles)
function editZ1_Callback(hObject, eventdata, handles)
function editZ2_Callback(hObject, eventdata, handles)
function editZ3_Callback(hObject, eventdata, handles)
function editZ4_Callback(hObject, eventdata, handles)
function editdZ1_Callback(hObject, eventdata, handles)
function editdZ2_Callback(hObject, eventdata, handles)
function editdZ3_Callback(hObject, eventdata, handles)
function editdZ4_Callback(hObject, eventdata, handles)
function editTimeStep_Callback(hObject, eventdata, handles)
function editDuration_Callback(hObject, eventdata, handles)
function editEnsembleSize_Callback(hObject, eventdata, handles)
function editEnsembleDuration_Callback(hObject, eventdata, handles)

%----------------------- CreateFcn -----------------------
% all CreateFcn do the same thing => generic call

function generic_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function editGamma1_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editGamma2_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editGamma3_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editGamma4_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editZ1_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editZ2_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editZ3_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editZ4_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editdZ1_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editdZ2_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editdZ3_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editdZ4_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editTimeStep_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editDuration_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editEnsembleSize_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);
function editEnsembleDuration_CreateFcn(hObject, eventdata, handles)
generic_CreateFcn(hObject, eventdata, handles);

%-------------------------- CallBacks ----------------------

function pushbuttonRunSimulation_Callback(hObject, eventdata, handles)
gam1 = ReadEdit(handles.editGamma1);
gam2 = ReadEdit(handles.editGamma2);
gam3 = ReadEdit(handles.editGamma3);
gam4 = ReadEdit(handles.editGamma4);

T  = ReadEdit(handles.editDuration);
dT = ReadEdit(handles.editTimeStep);

z1 = ReadEdit(handles.editZ1); z1=z1(1)+i*z1(2);
z2 = ReadEdit(handles.editZ2); z2=z2(1)+i*z2(2);
z3 = ReadEdit(handles.editZ3); z3=z3(1)+i*z3(2);
z4 = ReadEdit(handles.editZ4); z4=z4(1)+i*z4(2);
z = point_vortices([gam1 ; gam2 ; gam3 ; gam4], [z1 ; z2 ; z3; z4], T,dT);
handles.simu=z;

dz1 = ReadEdit(handles.editdZ1); z1=z1+dz1(1)+i*dz1(2);
dz2 = ReadEdit(handles.editdZ2); z2=z2+dz2(1)+i*dz2(2);
dz3 = ReadEdit(handles.editdZ3); z3=z3+dz3(1)+i*dz3(2);
dz4 = ReadEdit(handles.editdZ4); z4=z4+dz4(1)+i*dz4(2);
z = point_vortices([gam1 ; gam2 ; gam3 ; gam4], [z1 ; z2 ; z3; z4], T,dT);
handles.simu2=z;
handles.dt=dT;

figure;
plot(real(z'),-imag(z')); axis image
xlabel('x'); ylabel('y');
guidata(hObject, handles);

function pushbuttonEnsembleSimulate_Callback(hObject, eventdata, handles)
% hObject    handle to pushbuttonEnsembleSimulate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
gam1 = ReadEdit(handles.editGamma1);
gam2 = ReadEdit(handles.editGamma2);
gam3 = ReadEdit(handles.editGamma3);
gam4 = ReadEdit(handles.editGamma4);

T  = ReadEdit(handles.editEnsembleDuration);
dT = ReadEdit(handles.editTimeStep);
P  = ReadEdit(handles.editEnsembleSize);

z1 = perturb(handles.editZ1, handles.editdZ1,P);
z2 = perturb(handles.editZ2, handles.editdZ2,P);
z3 = perturb(handles.editZ3, handles.editdZ3,P);
z4 = perturb(handles.editZ4, handles.editdZ4,P);

ensemble = ensemble_vortices([gam1 ; gam2 ; gam3 ; gam4], [z1 z2 z3 z4], T,dT);
handles.ensemble = ensemble;
handles.dt=dT;
guidata(hObject, handles);

function pushbuttonAnimate_Callback(hObject, eventdata, handles)
% hObject    handle to pushbuttonAnimate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
z=handles.simu; dt=handles.dt;
x=real(z(:)); xmin=min(x); xmax=max(x);
y=imag(z(:)); ymin=min(y); ymax=max(y);
h=figure; 
for ii=20:5:size(z,2)
    zi=z(:,ii-19:ii);
    plot(real(zi'),-imag(zi')); 
    axis equal; axis([xmin xmax ymin ymax]);
    title(sprintf('t=%d',floor(dt*ii)));
    xlabel('x'); ylabel('y');
    getframe(h);
end

function pushbuttonEnsembleAnimate_Callback(hObject, eventdata, handles)
z=handles.ensemble; dt=handles.dt;
x=real(z(:)); xmin=min(x); xmax=max(x);
y=imag(z(:)); ymin=min(y); ymax=max(y);
h=figure;
for ii=20:5:size(z,3)
    zi=z(:,1,ii);
    plot(real(zi'),-imag(zi'),'r.'); 
    axis equal; axis([xmin xmax ymin ymax]);
    title(sprintf('t=%d',floor(dt*ii)));
    xlabel('x'); ylabel('y');
    getframe(h);
end

function pushPDF_Callback(hObject, eventdata, handles)
% hObject    handle to pushPDF (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
figure; 

subplot(1,2,1);
%z=handles.simu; x=real(z(4,:)); y=imag(z(4,:));
z=handles.simu; x=real(z(1,:)); y=imag(z(1,:));
[n1 xx yy]=hist2(x,y);
pcolor(xx,yy,n1); shading flat; axis image;
colormap(1-gray); colorbar
xlabel('x'); ylabel('y');

subplot(1,2,2);
%z=handles.simu2; x=real(z(4,:)); y=imag(z(4,:));
z=handles.simu2; x=real(z(1,:)); y=imag(z(1,:));
[n2 xx yy]=hist2(x,y);
pcolor(xx,yy,n2); shading flat; axis image;
colormap(1-gray); colorbar
xlabel('x'); ylabel('y');

figure; plot(n1(:),n2(:),'k.',n1(:),n1(:),'r--');
xlabel('Density (unperturbed simulation)');
ylabel('Density (perturbed initial condition)');

function pushbuttonEnsemblePDF_Callback(hObject, eventdata, handles)
%z=handles.ensemble; z=squeeze(z(:,4,:));
z=handles.ensemble; z=squeeze(z(:,1,:));
N=floor(size(z,1)/2); z1=z(1:N,:); z2=z(N+1:end,:);
z1=z1(:); z2=z2(:);
whos z z1 z2
figure; 

subplot(1,2,1);
[n1 xx yy]=hist2(real(z1),imag(z1));
pcolor(xx,yy,n1); shading flat; axis image;
colormap(1-gray); colorbar

subplot(1,2,2);
[n2 xx yy]=hist2(real(z2),imag(z2));
pcolor(xx,yy,n2); shading flat; axis image;
colormap(1-gray); colorbar

figure; plot(n1(:),n2(:),'k.');
xlabel('Density (first ensemble)');
ylabel('Density (second ensemble)');

function pushbuttonSpread_Callback(hObject, eventdata, handles)
z=real(handles.simu); z1=squeeze(z(1,:));
z=real(handles.simu2); z2=squeeze(z(1,:));
dz=abs(z1-z2); z=[z1 ; z2];

dt = handles.dt; t=(1:numel(dz))*dt;
h=figure;
subplot(3,1,1); plot(t,z1,'k',t,z2,'r--');
xlabel('time'); ylabel('x_1'); legend('control','perturbed');
subplot(3,1,2); plot(t,dz,'b');
xlabel('time'); ylabel('\delta x_1');
subplot(3,1,3); semilogy(t,dz,'b');
xlabel('time'); ylabel('\delta x_1');

function pushbuttonEnsembleSpread_Callback(hObject, eventdata, handles)
z=real(handles.ensemble); z=squeeze(z(:,1,:)); dz=std(z);
dt = handles.dt; t=(1:numel(dz))*dt;
h=figure;
subplot(3,1,1); plot(t,z);
xlabel('time'); ylabel('x_1');
subplot(3,1,2); plot(t,dz);
xlabel('time'); ylabel('\delta x_1');
subplot(3,1,3); semilogy(t,dz);
xlabel('time'); ylabel('\delta x_1');

%------------------------ Functions ----------------------------------

function val = ReadEdit(edit)
val = eval(get(edit, 'String'));
function z=point_vortices(gam,z0,T,dt)
N=ceil(T/dt); dt=T/N; M=numel(gam);
zcur = z0; z=zeros(M,N);
elapsed=0;
for ii=1:N
    tic;
    % inner iteration solving for znext
    % znext = zcur + dt*f((zcur+znext)/2)
    znext=zcur;
    for jj=1:10
        zmid = conj(znext+zcur)/2; % dz/dt = i/conj(z)
        dz = 0*zmid; % dz/dt
        for kk=1:M
            dzk = 0; % d/dt z(kk)
            for ll=1:M
                if(kk ~= ll)
                    dzk = dzk + i*gam(ll)/(zmid(kk)-zmid(ll));
                end
            end
            dz(kk)=dzk;
        end
        znext = zcur+dt*dz;
    end
    z(:,ii) = zcur;
    zcur = znext;
    elapsed=elapsed+toc;
    if(elapsed>0.5)
        elapsed = elapsed-0.5;
        fprintf('Outer iteration %d/%d\n', ii,N);
    end
end

function [n xx yy]=hist2(x,y)
xx=min(x)+(max(x)-min(x))*(0:.05:1);
yy=min(y)+(max(y)-min(y))*(0:.05:1);
n=zeros(numel(yy)-1,numel(xx)-1);
for ii=2:numel(yy)
    for jj=2:numel(xx)
        n(ii-1,jj-1) = sum((x>=xx(jj-1)) & (x<xx(jj)) & (y>=yy(ii-1)) & (y<yy(ii)) );
    end
end
xx=(xx(2:end)+xx(1:end-1))/2; 
yy=(yy(2:end)+yy(1:end-1))/2; 



function z=perturb(editZ, editdZ, P)
z=ReadEdit(editZ); dz=ReadEdit(editdZ);
x = z(1)+dz(1)*randn(P,1);
y = z(2)+dz(2)*randn(P,1);
z=x+i*y;
function z=ensemble_vortices(gam,z0,T,dt)
% z0 = PxM array

N=ceil(T/dt); dt=T/N; M=numel(gam); P=size(z0,1);
zcur = z0; z=zeros(P,M,N);
elapsed=0;
for ii=1:N
    tic;
    % inner iteration solving for znext
    % znext = zcur + dt*f((zcur+znext)/2)
    znext=zcur; % PxM
    for jj=1:10
        zmid = conj(znext+zcur)/2; % dz/dt = i/conj(z)
        dz = zeros(P,M); % dz/dt
        for kk=1:M
            dzk = zeros(P,1); % d/dt z(kk)
            for ll=1:M
                if(kk ~= ll)
                    dzk = dzk + i*gam(ll)./(zmid(:,kk)-zmid(:,ll));
                end
            end
            dz(:,kk)=dzk;
        end
        znext = zcur+dt*dz;
    end
    z(:,:,ii) = zcur;
    zcur = znext;
    elapsed=elapsed+toc;
    if(elapsed>0.5)
        elapsed = elapsed-0.5;
        fprintf('Outer iteration %d/%d\n', ii,N);
    end
end


