function install(varargin)
% INSTALL   Install script for NESTool.
%    INSTALL creates the NESTool folder and (optionally) 
%    the needed entries in the startup.m file.
% 
%    INSTALL PATH creates the NESTool folder 
%    in the specified PATH.
%    
%    This installation script was generated by using 
%    the MAKEINSTALL tool. For further information
%    visit http://matlab.pucicu.de

% Copyright (c) 2008-2012
% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
% http://www.pik-potsdam.de
%
% Copyright (c) 2001-2008
% Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
% THIS IS A GENERATED INSTALL-FILE, DO NOT EDIT!
% Generation date: 18-Mar-2014 18:34:46
% $Date: 2014/03/12 16:17:50 $
% $Revision: 3.30 $

install_file='';install_path='';installfile_info.date='';installfile_info.bytes=[];
time_stamp='';checksum='';checksum_file=''; instpaths = '';
errcode=0;

try
  warning('off')
  if isoctave
     disp(['  You are trying to install the toolbox in Octave.',10,'  Some compatibility issues might appear.'])
     warning('off','Octave:possible-matlab-short-circuit-operator')
     more off
  end
  if nargin
    install_path = varargin{1};
  end

  if exist('install.log','file') == 2, delete('install.log'), end
  %rehash
  disp('------------------------')
  disp('  INSTALLATION NESTool');
  disp('------------------------')
  install_file=[mfilename,'.m'];
  currentpath=pwd; time_stamp='time_stamp not yet obtained'; checksum='checksum not yet obtained';
  fid = 0;
  
%%%%%%% read the archive
%%%%%%% and look for checksum and date in archive
  errcode=90;
  disp('  Reading the archiv ')
  fid=fopen(install_file,'r'); 
  fseek(fid,0,'eof'); eofbyte=ftell(fid);
  fseek(fid,30910,'bof'); % location where the container starts
  while 1
     temp=fgetl(fid);
     startbyte=ftell(fid);
     if length(temp)>1
       if strcmpi(temp,'%<-- Header begins here -->')
          errcode=90.1;
          checksum=fgetl(fid);
          temp1=fgetl(fid);
          temp2=fgetl(fid);
       end
       if strcmpi(temp,'%<-- Header ends here -->')
          startbyte=ftell(fid);
          break
       end
     end
  end
  checksum(1:2)=[];
  fseek(fid,startbyte,'bof');
  errcode=90.2;
  A=fread(fid,eofbyte);
  errcode=90.3;
  checksum_file=dec2hex(sum((1:length(A))'.*A));
  if ~strcmpi(checksum_file,checksum)
    error(['The installation file is corrupt!',10,'Ensure that the archive container was ',...
           'not modified (check FTP/ ',10,'proxy/ firewall settings, anti-virus scanner for emails etc.)!'])
  else
    disp(['  Checksum test passed (', checksum,')'])
  end
  fclose(fid);

  disp(['  NESTool version ', temp2(3:end),''])
  time_stamp=temp1(3:end); disp(['  NESTool time stamp ', time_stamp,'']); 
  
  errcode=91;
  if isunix
    toolboxpath='NESTool';
  else
    toolboxpath='NESTool';
  end
  
  
%%%%%%% check for older versions
  
  p=path; i1=0;
  rem_old = '';
  
  while any([findstr([lower(toolboxpath),'demo'],lower(p)) findstr(lower(toolboxpath),lower(p)) findstr('',lower(p))]>i1) && ~strcmpi('N',rem_old)
    errcode=92;
    i1=[findstr([lower(toolboxpath),'demo'],lower(p)) findstr(lower(toolboxpath),lower(p)) findstr('',lower(p))];
    if ~isempty(i1)
      i1=i1(1);
      if isunix, i2=findstr(':',p); else i2=findstr(';',p); end
      i3=i2(i2>i1);                 % last index pathname
      if ~isempty(i3), i3=i3(1)-1; else i3=length(p); end
      i4=i2(i2<i1);                 % first index pathname
      if ~isempty(i4), i4=i4(end)+1; else i4=1; end
      oldtoolboxpath=p(i4:i3);
      if isempty(rem_old)
          disp(['  Old NESTool found in ', oldtoolboxpath,''])
          rem_old = input('> Delete old toolbox? Y/N [Y]: ','s');
      end
      if isempty(rem_old), rem_old = 'Y'; end
      if strcmpi('Y',rem_old)
%%%%%%% removing old entries in startup-file
        errcode=94;
        rmpath(oldtoolboxpath)
        err = savepath;
        if err, disp('  ** Warning: No write access to pathdef.m file!'), end
        if i4>1, p(i4-1:i3)=''; else p(i4:i3)=''; end
        startup_exist = exist('startup','file');
        if isoctave startup_exist = exist(fullfile('~','.octaverc'),'file'); end
        if startup_exist
             startupfile=which('startup');
             startuppath=startupfile(1:findstr('startup.m',startupfile)-1);
             if isoctave
                startuppath = ['~',filesep];
                startupfile = fullfile('~','.octaverc');
             end
             errcode=94.1;
             if ~isunix
               if isoctave
                   toolboxroot=fullfile(matlabroot);
               else
                   toolboxroot=fullfile(matlabroot,'toolbox');
               end
               curr_pwd = pwd; home_pwd = matlabroot; 
             else
               toolboxroot=startuppath;
               curr_pwd = pwd; cd ('~'); home_pwd = pwd; cd(curr_pwd);
             end
             fid = fopen(startupfile,'r');
             k = 1;
             while 1
                 tmp = fgetl(fid);
                 if ~ischar(tmp), break, end
                 instpaths{k} = tmp;
                 k = k + 1;
             end
             fclose(fid);
             k=1;
             while k <= length(instpaths)
               if findstr(oldtoolboxpath,strrep(instpaths{k},'~',home_pwd))
                 errcode=94.2;
                 instpaths(k)=[];
               else
                 k=k+1;
               end
             end
             fid=fopen(startupfile,'w');
             errcode=94.3;
             if fid < 0
               disp(['  ** Warning: Could not get access to ',startupfile,'.']);
               disp('  ** Could not remove toolbox from the startup.m file.');
               disp('  ** Ensure that you have write access!');
             else
               for i2=1:length(instpaths), 
                 fprintf(fid,'%s\n', char(instpaths{i2})); 
               end
               fclose(fid);
             end
        end
%%%%%%% removing old paths
        errcode=93;
        if exist(oldtoolboxpath,'dir') == 7
           disp(['  Change to ',oldtoolboxpath,''])
           cd(oldtoolboxpath)
           dirnames='';filenames='';
           temp='.:';
           errcode=93.1;
           while ~isempty(temp)
             [temp1 temp]=strtok(temp,':');
             if ~isempty(temp1)
               dirnames=[dirnames; {temp1}];
               x2=dir(temp1);
               for i=1:length(x2)
                 if ~x2(i).isdir, filenames=[filenames; {[temp1,'/', x2(i).name]}]; end
         	   if x2(i).isdir && ~strcmp(x2(i).name,'.') && ~strcmp(x2(i).name,'..'), temp=[temp,temp1,filesep,x2(i).name,':']; end
               end
             end
           end
           errcode=93.2;
           if isoctave, confirm_recursive_rmdir (false, 'local'); end
           dirnames=strrep(dirnames,['.',filesep],'');
           for i=1:length(dirnames),l(i)=length(dirnames{i}); end
           [i i4]=sort(l);
           dirnames=dirnames(fliplr(i4));
           errcode=93.3;
           for i=1:length(dirnames)
              if dirnames{i} == '.', continue, end
              delete([dirnames{i}, filesep,'*']),
              if exist('rmdir') == 5 && exist(dirnames{i}) == 7, rmdir(dirnames{i},'s'); else, delete(dirnames{i}), end
              disp(['  Removing files in ',char(dirnames{i}),''])
           end
           errcode=93.4;
           cd(currentpath)
           if exist('rmdir') == 5 && exist(oldtoolboxpath) == 7, rmdir(oldtoolboxpath,'s'); else, delete(oldtoolboxpath), end
           errcode=93.5;
           disp(['  Removing ',oldtoolboxpath,''])
        end
%%%%%%%
      end
    end
    p = path; i1 = 0;
  end
  clear p i i1 i2 i3 i4 temp* x2
  
%%%%%%% add entry into startpath in startup.m
  i=findstr(toolboxpath,path);
  startupPos = 0;
  startup_exist = exist('startup','file');
  if isoctave startup_exist = exist(fullfile('~','.octaverc'),'file'); end
  if startup_exist
        errcode=95.1;
        startupfile=which('startup');
        startuppath=startupfile(1:findstr('startup.m',startupfile)-1);
        if isoctave
           startuppath = ['~',filesep];
           startupfile = fullfile('~','.octaverc');
        end
  
        if ~isunix
           errcode=95.11;
           %toolboxroot=fullfile(matlabroot,'toolbox');
           if isoctave
               toolboxroot=matlabroot;
           else
               toolboxroot=strtok(userpath,';');
           end
        else
           errcode=95.12;
           toolboxroot=startuppath;
        end
        fid = fopen(startupfile,'r');
        k = 1;
        while 1
            tmp = fgetl(fid);
            if ~ischar(tmp), break, end
            instpaths{k} = tmp;
            k = k + 1;
        end
        fclose(fid);
  end
  if isempty(i)
     errcode=95;
     startupfilestr = 'startup.m';
     if isoctave startupfilestr = fullfile('~','.octaverc'); toolboxroot = pkg('prefix'); end

     %% check whether the default toolbox path exists
     err = 1;
     if ~isunix
           errcode=95.01;
           if isoctave
               toolboxroot=matlabroot;
           else
               if isempty(userpath), userpath('reset'), end
               toolboxroot=strtok(userpath,';');
           end
           if exist(toolboxroot,'file') ~= 7, err=mkdir(toolboxroot); end
           cd(toolboxroot)
           startupfile=fullfile(toolboxroot,'startup.m');
           if isoctave, startupfile=fullfile('~','.octaverc');end
           instpaths={''};
     else
           errcode=95.02;
           cd ~
           startuppath=[pwd,filesep];
           if isoctave
              status = fileattrib(pkg('prefix'));
              result = stat(pkg('prefix'));
              uid = geteuid;
              if status && result.uid ~= uid
                  disp(['  ** No writeable package folder found. Will define it as ~/octave.',10,'     If installation fails, you have to create it manually, set it by',10,'     calling ''pkg prefix ~/octave'' and re-run the installation!'])
                  err=mkdir('~/octave');
                  pkg prefix ~/octave;
                  toolboxroot = pkg('prefix');
              end

              if exist(pkg('prefix'),'file') ~= 7, err=mkdir(pkg('prefix')); end
              cd(pkg('prefix'))
  	          startupfile=fullfile(startuppath,'.octaverc');
           else
              if isempty(userpath), userpath('reset'), end
              up = textscan(userpath,'%s','delimiter',':'); up=up{1};
              startuppath = '';
              for k = 1:length(up)
                  if exist(up{k},'file') == 7
                      startuppath = up{k};
                  end
              end
              if isempty(startuppath)
                  err=mkdir(up{1});
                  startuppath = up{1};
              end
              cd(startuppath)
              startupfile=fullfile(startuppath,'startup.m');
           end
     end
     if ~err, error('Could not create toolbox path. Please check whether you have write access or whether there is another file of such a name blocking the creation of the path.'), end
           
     if exist(startupfilestr,'file')
        errcode=95.1;
     else
        errcode=95.2;
        if ~isunix
           errcode=95.21;
           if isoctave
               toolboxroot=matlabroot;
               startupfile=fullfile('~','.octaverc');
           else
               toolboxroot=strtok(userpath,';');
               startupfile=fullfile(toolboxroot,'startup.m');
           end
           cd(toolboxroot)
           instpaths={''};
        else
           errcode=95.22;
           cd ~
           startuppath=[pwd,filesep];
           if isoctave
              cd(pkg('prefix'))
              startupfile=fullfile(startuppath,'.octaverc');
           else
              up = textscan(userpath,'%s','delimiter',':'); up=up{1};
              startuppath = '';
              for k = 1:length(up)
                  if exist(up{k},'file') == 7
                      startuppath = up{k};
                  end
              end
              if isempty(startuppath)
                  err=mkdir(up{1});
                  startuppath = up{1};
              end
              cd(startuppath)
              startupfile=fullfile(startuppath,'startup.m');
           end
           
  	     toolboxroot=startuppath;
  	     instpaths={''};
        end
     end
    
     errcode=95.23;
     if ~isempty(install_path)
           switch ( exist(install_path,'dir') )
              case 0
                in = input(['> Create ', install_path, '? Y/N [Y]: '],'s');
                if isempty(in), in = 'Y'; end
                if strcmpi('Y',in)
                   err = mkdir(install_path);
                   if ~err
                     disp(['  ** Could not create ', install_path, '! Using ',startuppath,' as installation path.'])
                   else                
                     toolboxroot = install_path;
                   end                
                else 
                   disp(['  ** Do not create ', install_path, '! Using ',startuppath,' as installation path.'])
                end
              case 2
                disp(['  ** ', install_path, ' is not a directory! Using ',startuppath,' as installation path.'])
              case 7
                toolboxroot = install_path;
           end
     end

     errcode=95.3;
     TBfullpath=fullfile(toolboxroot,toolboxpath);
     if ~exist(TBfullpath,'dir'), mkdir(toolboxroot,toolboxpath); end

%%%%%%% resolve relative path (starting with ./ and ../) to absolute path

     isrelpath = findstr('./',TBfullpath);
     isrelpathDos = findstr('.\',TBfullpath);
     if ( ~isempty(isrelpath) && isrelpath == 1 ) || ( ~isempty(isrelpathDos) && isrelpathDos == 1 )
         TBfullpath = fullfile(pwd, TBfullpath(3:end));
     end

     isrelpath = findstr('../',TBfullpath);
     isrelpathDos = findstr('..\',TBfullpath);
     if ( ~isempty(isrelpath) && isrelpath == 1 ) || ( ~isempty(isrelpathDos) && isrelpathDos == 1 )
         TBfullpath = fullfile(pwd, TBfullpath(4:end));
     end
     

%%%%%%% ask where to add entry in startup file

     disp(['> In order to get permanent access, the toolbox should be added',10,'> to the top (default) or end (E) of your startup path.'])
     in = input('> Add toolbox permanently into your startup path (highly recommended)? Y/E/N [Y]: ','s');
     if isempty(in), in = 'Y'; end
     if strcmpi('Y',in)
       startupPos = '-begin';
       disp('  Adding Toolbox at the top of the startup.m file')
     elseif strcmpi('E',in)
       startupPos = '-end';
       disp('  Adding Toolbox at the end of the startup.m file')
     end

     if startupPos
         errcode=95.4;
         loc = ['addpath ''',TBfullpath,''' ', startupPos];
         if ~ismember(loc, instpaths)
             instpaths{end+1} = loc;
         end
     end

  else
     errcode=96;
     startupfile=which('startup');
     startuppath=startupfile(1:findstr('startup.m',startupfile)-1);
     if isoctave
         startuppath = ['~', filesep];
         startupfile = fullfile('~','.octaverc');
     end
     if ~isunix
        if isoctave
            toolboxroot=matlabroot;
        else
            toolboxroot=strtok(userpath,';');
        end
        %toolboxroot=fullfile(matlabroot,'toolbox');
     else
        toolboxroot=startuppath;
     end
    
        errcode=96.21;
        if ~isempty(install_path)
           switch ( exist(install_path,'dir') )
              case 0
                disp(['> Create ', install_path, '?'])
                in = input('> Create ', install_path, '? Y/N [Y]: ','s');
                if isempty(in), in = 'Y'; end
                if strcmpi('Y',in)
                   err = mkdir(install_path);
                   if ~err
                     disp(['  ** Could not create ', install_path, '! Using ',startuppath,' as installation path.'])
                   else                
                     toolboxroot = install_path;
                   end                
                else 
                   disp(['  ** Do not create ', install_path, '! Using ',startuppath,' as installation path.'])
                end
              case 2
                disp(['  ** ', install_path, ' is not a directory! Using ',startuppath,' as installation path.'])
              case 7
                toolboxroot = install_path;
           end
        end
        TBfullpath=fullfile(toolboxroot,toolboxpath);
        if ~exist(TBfullpath,'dir'), mkdir(toolboxroot,toolboxpath); end
  end
  
  
%%%%%%% read the archive
  errcode=97;
  fprintf('  Importing the archiv ')
  cd(currentpath)
  max_sc=30;
  scale=[sprintf('%3.0f',0),'% |',repmat(' ',1,max_sc),'|'];
  fprintf('%s',scale)
  b=''; c.file=''; c.data=''; bfile=''; folder=''; i2=0;
    
    % read ASCII
    errcode=97.1;
    i1=findstr(char(A'),'%<-- ASCII begins here');
    i2=findstr(char(A'),'%<-- ASCII ends here -->');
    
    for k=1:length(i1),
      wb=round(i1(k)*max_sc/eofbyte);
      errcode=97.11;
      scale=[sprintf('%3.0f',100*wb/max_sc),'% |',repmat('*',1,wb),repmat(' ',1,max_sc-wb),'|'];
      fprintf([repmat('\b',1,max_sc+7),'%s'],scale)
      B=A(i1(k):i2(k)-2);
     
      i4=find(B == 10);
      i3=[0;i4(1:end-1)]+1;
      filename=strrep(...
                strrep(char(B(i3(1):i4(1)-1)'),'%<-- ASCII begins here: __',''),...
  	      '__ -->','');
      errcode=['97.1',reshape(dec2hex(double(filename))',1,length(filename)*2)];
      filename = strrep(filename, '/', filesep);
      i6=findstr(filename, filesep);
      if i6>0
         folder=[folder; {filename(1:i6(end)-1)}];
      end
      c(k).file={filename};
      c(k).data=strrep(char(B(i4(1)+3:end)'),[char(10),'%@'],char(10));
    end
  
    % read binary
    errcode=97.2;
    i1=findstr(char(A'),'%<-- Binary begins here');
    i2=findstr(char(A'),'%<-- Binary ends here -->');
    for k=1:length(i1),
      wb=round(i1(k)*max_sc/eofbyte);
      errcode=97.21;
      scale=[sprintf('%3.0f',100*wb/max_sc),'% |',repmat('*',1,wb),repmat(' ',1,max_sc-wb),'|'];
      fprintf([repmat('\b',1,max_sc+7),'%s'],scale)
      B=A(i1(k):i2(k));
  
      i4=find(B == 10);
      i3=[0;i4(1:end-1)]+1;
      i5=findstr(char(B(i3(1):i4(1))'),'__');
      filename=char(B(i5(1)+2:i5(2)-1)');
      errcode=['97.2',reshape(dec2hex(double(filename))',1,length(filename)*2)];
      filename = strrep(filename, '/', filesep);
      i6=findstr(filename, filesep);
      if i6>0
      folder=[folder; {filename(1:i6(end)-1)}];
      end
      bfile=[bfile; {filename}];
      nbytes=str2double(char(B(i5(3)+2:i5(4)-1)'));
      if nbytes>=2
        temp=reshape(B(i4(1)+1:i4(1)+nbytes),2,nbytes/2);
        b=[b;{temp(2,:)}];
      else
        b=[b;{''}];
      end
  
    end
    wb=max_sc;
    scale=[sprintf('%3.0f',100*wb/max_sc),'% |',repmat('*',1,wb),repmat(' ',1,max_sc-wb),'|'];
    fprintf([repmat('\b',1,max_sc+7),'%s'],scale)
  fprintf('\n')
  clear temp* i i1 i2 i3 i4 i5 i6 A B
  
  errcode='97.3';
  if exist(TBfullpath,'dir') ~= 7, error(['Could not enter toolbox path',10,'  ', TBfullpath,10,'. Please check whether you have write access or whether there is another file of such a name blocking the creation of the path.']), end
  cd(TBfullpath)
  disp(['  Toolbox folder is ',TBfullpath,''])
  
%%%%%%% make sub-directories
  errcode='97.4';
  for i=1:length(folder);
     i1=folder{i}; i2='.'; olddir=pwd;
     if exist([TBfullpath, filesep, i1],'file') ~= 7
        while ~isempty(i2) && ~isempty(i1)
          cd(i2)
          [i2 i1]=strtok(i1,filesep);
          if exist([pwd, filesep, i2],'file') ~= 7 && exist([pwd, filesep, i2],'file')
             disp(['  ** Warning: Found ', i2, ' will be overwritten.'])
             errcode=['97.5',reshape(dec2hex(double(i2))',1,length(i2)*2)];
             delete(i2)
          end
          if ~exist([pwd, filesep, i2],'file')
            disp(['  Make directory ',pwd, filesep, i2])
            errcode=['97.6',reshape(dec2hex(double(i2))',1,length(i2)*2)];
            err = mkdir(i2);
            if err==0, disp(['  ** Warning: Could not create ',pwd, filesep, i2,'!',10,'     Installation will probably fail.']), end
            errcode=97.7;
            if ~(strcmpi(i2,'private') || strcmpi(i2(1),'+')) && any(i2 ~= '@')
                 % add entries in startup.m
                 if isempty(startupPos) startupPos = '-end'; end
                 loc = ['addpath ''',pwd, filesep, i2,''' ', startupPos];
                 eval(loc);
                 if ~any(ismember(loc, instpaths)) && all(startupPos~=0)
                     instpaths{end+1} = loc;
                 end
            end
          end
        end
        cd(olddir)
     end
  end
  
%%%%%%% make toolbox accessible in the current matlab session
  if ~startupPos, startupPos = '-end'; end
  addpath(TBfullpath,startupPos);


%%%%%%% write startup file
  if startupPos
    errcode=95.4;
    if ~exist('startupfile','var')
        if isoctave
            startupfile = fullfile('~','.octaverc');
        else
            up = textscan(userpath,'%s','delimiter',':'); up=up{1};
            startuppath = '';
            for k = 1:length(up)
                if exist(up{k},'file') == 7
                    startuppath = up{k};
                end
            end
            if isempty(startuppath)
                err=mkdir(up{1});
                startuppath = up{1};
            end
            startupfile = fullfile(startuppath,'startup.m');
        end
    end

    fid=fopen(startupfile,'w');
    if fid < 0
      disp(['  ** Warning: Could not get access to ',startupfile,'.']);
      disp('  ** Could not add toolbox into the startup.m file.');
      disp('  ** Ensure that you have write access!');
    else
      %for i2=1:length(instpaths), fprintf(fid,'%s\n', char(instpaths{i2})); disp(char(instpaths{i2}));end
      for i2=1:length(instpaths), fprintf(fid,'%s\n', char(instpaths{i2})); end
      fclose(fid);
    end
  end


%%%%%%% write the programme files
  errcode=98;
  num_errors=0;
  for i=1:length(c),
    disp(['  Creating ',char(c(i).file)])
    fid=fopen(char(c(i).file),'w');
    if fid < 0
       disp(['  ** Warning: Could not get access to ',char(c(i).file),'.']);
       disp('  ** Ensure that you have write access in this filesystem!');
       num_errors=num_errors+1;
       if num_errors == 2; 
         disp('Abort!')
         disp('Too much errors due to write access failure in this filesystem.')
         cd(currentpath), return
       end
    else
      if strcmpi(c(i).file,'info.xml')
        v=version;
        release=str2double(v(findstr(v,'(R')+2:findstr(v,')')-1));
        if release>12, area='toolbox'; icon_path='$toolbox/matlab/icons'; else area='matlab'; icon_path='$toolbox/matlab/general'; end
        i3=findstr(c(i).data,'<matlabrelease>'); i4=findstr(c(i).data,'</matlabrelease>');
        if ~isempty(i3) && i4>i3;
          c(i).data=strrep(c(i).data,c(i).data(i3:i4-1),['<matlabrelease>',num2str(release)]);
        end
        i3=findstr(c(i).data,'<area>'); i4=findstr(c(i).data,'</area>');
        if ~isempty(i3) && i4>i3;
          c(i).data=strrep(c(i).data,c(i).data(i3:i4-1),['<area>',area]);
        end
        c(i).data=strrep(c(i).data,'<icon>$toolbox/matlab/general',['<icon>',icon_path]);
        c(i).data=strrep(c(i).data,'<icon>$toolbox/matlab/icons',['<icon>',icon_path]);
      end
      fprintf(fid,'%s',char(c(i).data));
      fclose(fid);
    end
  end

%%%%%%% pcode the programme files
  for i=1:length(c),
    try
      [tPath tFile tExt]=fileparts(char(c(i).file));
      if strcmpi(tExt,'.m') && ~strcmpi(tFile,'Readme') && ~strcmpi(tFile,'Contents')
        if mislocked(char(c(i).file)), munlock(char(c(i).file)); clear(char(c(i).file)); end
        disp(['  Pcode ',char(c(i).file),''])
        if sum(c(i).file == '.') < 2
           pcode(char(c(i).file),'-inplace')
        end  
      end  
    catch
    end
  end

%%%%%%% write the binary files
  num_errors=0;
  for i=1:length(b),
    disp(['  Creating ',char(bfile(i)),''])
    fid=fopen(char(bfile(i)),'w');
    if fid < 0
       disp(['  ** Warning: Could not get access to ',char(bfile(i)),'.']);
       disp('  ** Ensure that you have write access in this filesystem!');
       num_errors=num_errors+1;
       if num_errors == 2; 
         disp('Abort!')
         disp('Too much errors due to write access failure in this filesystem.')
         cd(currentpath), return
       end
    else
      fwrite(fid,b{i}); 
      fclose(fid);
    end
  end
  tx=version; tx=strtok(tx,'.'); 
  if str2double(tx)>=5 && exist('rehash','builtin'), rehash, end
  if str2double(tx)>=6 && exist('rehash','builtin'), eval('rehash toolboxcache'), end
  
%%%%%%% removing installation file
  errcode=99;
  cd(currentpath)
  i = input('> Delete installation file? Y/N [Y]: ','s');
  if isempty(i), i = 'Y'; end
  
  if strcmpi('Y',i)
    disp('  Removing installation file')
    delete(install_file)
  end
  
  disp('  Installation finished!')
  
  if ~exist('rehash','builtin')
    disp('  ** Warning: Could not rehash your Matlab system.')
    disp('  ** Probably a restart of Matlab will be necessary in order')
    disp('  ** to get access to the installed toolbox.')
  end
  
  disp('------------------------')
  
  disp('For an overview type:')
  disp(['helpwin ',toolboxpath])
  warning('on')
  if isoctave
      more on
  end
  
  
%%%%%%% error handling

catch
  z2=whos;x_lasterr=lasterr;y_lastwarn=lastwarn;
  if ~strcmpi(x_lasterr,'Interrupt')
    if fid>-1, 
      try, z_ferror=ferror(fid); catch, z_ferror=''; end
    else
      z_ferror='File not found.'; 
    end
    installfile_info=dir([currentpath,filesep,install_file]);
    fid=fopen(fullfile(currentpath,'install.log'),'w');
    checksum_test=findstr(x_lasterr,'The installation file is corrupt!');
    if isempty(checksum_test),checksum_test=0; end
    if ~checksum_test
      fprintf(fid,'%s\n','A critical error has occurred. Please inform the distributor');
      fprintf(fid,'%s\n','of the toolbox, where the error occured and send us the entire');
      fprintf(fid,'%s\n','screen output of the installation, the following error');
      fprintf(fid,'%s\n','report, and the informations about the toolbox (distributor,');
      fprintf(fid,'%s\n','name, URL etc.). Provide a brief description of what you were');
      fprintf(fid,'%s\n','doing when this problem occurred.');
      fprintf(fid,'%s\n','E-mail or FAX this information to us at:');
      fprintf(fid,'%s\n','    E-mail:  marwan@pik-potsdam.de');
      fprintf(fid,'%s\n','       Fax:  ++49 +331 288 2640');
      fprintf(fid,'%s\n\n\n','Thank you for your assistance.');
      fprintf(fid,'%s\n',repmat('-',50,1));
      fprintf(fid,'%s\n',datestr(now,0));
      fprintf(fid,'%s\n',['Matlab ',char(version),' on ',computer]);
      fprintf(fid,'%s\n',repmat('-',50,1));
      fprintf(fid,'%s\n','Makeinstall Version ==> 3.30 ');
      fprintf(fid,'%s\n',['Install File ==> ',install_file,'/',installfile_info.date,'/',num2str(installfile_info.bytes)]);
      fprintf(fid,'%s\n',['Container ==> ',time_stamp,'/',checksum]);
      fprintf(fid,'%s\n\n',repmat('-',50,1));
      fprintf(fid,'%s\n',x_lasterr);
      fprintf(fid,'%s\n',y_lastwarn);
      fprintf(fid,'%s\n',z_ferror);
      fprintf(fid,'%s\n',[' errorcode ==> ',num2str(errcode)]);
      fprintf(fid,'%s\n',' workspace dump ==>');
      if ~isempty(z2), 
        fprintf(fid,'%s\n',['Name',char(9),'Size',char(9),'Bytes',char(9),'Class']);
        for j=1:length(z2);
          fprintf(fid,'%s',[z2(j).name,char(9),num2str(z2(j).size),char(9),num2str(z2(j).bytes),char(9),z2(j).class]);
          if ~strcmp(z2(j).class,'cell') && ~strcmp(z2(j).class,'struct')
            content=eval(z2(j).name);
            try, content=mat2str(content(1:min([size(content,1),500]),1:min([size(content,2),500])));end
            fprintf(fid,'\t%s',content(1:min([length(content),500])));
          elseif strcmp(z2(j).class,'cell')
            content=eval(z2(j).name);
            fprintf(fid,'\t');
            for j2=1:min([length(content),500])
              if isnumeric(content{j2})
                fprintf(fid,'{%s} ',content{j2}(1:end));
              elseif iscell(content{j2})
                fprintf(fid,'{%s} ',content{j2}{1:end});
              end
            end
          elseif strcmp(z2(j).class,'struct')
            content=fieldnames(eval(z2(j).name));
            content=char(content); content(:,end+1)=' '; content=content';
            fprintf(fid,'\t%s',content(:)');
          end
          fprintf(fid,'%s\n','');
        end
      end
    else
      fprintf(fid,'%s\n','Installation aborted due to a failed checksum test!');
      fprintf(fid,'%s\n',['Checksum should be:     ', checksum]);
      fprintf(fid,'%s\n\n',['Checksum of archive is: ', checksum_file]);
      fprintf(fid,'%s\n','Ensure that the installation file was not modified by any');
      fprintf(fid,'%s\n','other programme, as an anti-virus scanner for emails, a');
      fprintf(fid,'%s\n','mis-configured HTTP proxy or FTP programme.');
    end
    fclose(fid);
    disp('----------------------------');
    disp('       ERROR OCCURED ');
    disp('    during installation');
    disp('----------------------------');
    disp(x_lasterr);
    if errcode == 95.21
       disp('----------------------------');
       disp('   This error means that the toolbox directory could not');
       disp('   be created. Probably, there is a file of the same name.');
       disp('   or there are in-appropriate permission settings in its');
       disp('   parent folder.');
       disp('   Please check the output of the following command:');
       disp('      userpath(''reset'')');
       disp('   which might help to locate the problem.');
       disp('----------------------------');
    end
    if ~checksum_test
      disp(z_ferror);
      disp(['   errorcode is ',num2str(errcode)]);
      disp('----------------------------');
      disp('   A critical error has occurred. Please inform the distributor');
      disp('   of the toolbox, where the error occured and send us the entire');
      disp('   screen output of the installation, the error report report');
      disp('   and the informations about the toolbox (distributor, name,');
      disp('   URL etc.). For your convenience, this information has been')
      disp('   recorded in: ')
      disp(['   ',fullfile(currentpath,'install.log')]), disp(' ')
      disp('   Provide a brief description of what you were doing when ')
      disp('   this problem occurred.'), disp(' ')
      disp('   E-mail or FAX this information to us at:')
      disp('       E-mail:  marwan@pik-potsdam.de')
      disp('          Fax:  ++49 +331 288 2640'), disp(' ')
      disp('   Thank you for your assistance.')
    end
  end
  warning('on')
  cd(currentpath)
  if isoctave
      more on
  end
end

function flag = isoctave
% ISOCTAVE   Checks whether the code is running in Octave
%   ISOCTAVE is returning the value TRUE if executed within the
%   Octave environment, else it is returning FALSE (e.g. when
%   called within Matlab.

a = ver('Octave');

if ~isempty(a) && strfind(a(1).Name,'Octave')
    flag = true;
else
    flag = false;
end
% -------------------------------------------
% GENERATED ENTRIES - DO NOT MODIFY ANYTHING!
%<-- Header begins here -->
%@76ACC88690
%@18-Mar-2014 18:34:46
%@1.01 
%<-- Header ends here -->
%<-- ASCII begins here: __Contents.m__ -->
%@% NESTool
%@% Version 1.01 18-Mar-2014
%@%
%@%    mihist                -  computes the mutual information between variables X and Y.
%@%    panlayerchronensemble -  creates alternative timescales for layer counted data.
%@%    ar1sur                -  computes autocorrelated irregular surrogate time series.
%@%    car                   -  gives irregular time series for coupled autoregressive (AR1) processes.
%@%    deinstall_nestool     -    Removes NESTool.
%@%    eventsynchro          -  calculates event synchronization for irregular time series.
%@%    extractwindow         -  cuts (irregular) time series according to the boundaries.
%@%    gausssmooth           -  smoothes irregularly sampled time series with gaussian weights.
%@%    generate_t            -  generates arbitrary time axes with specified properties.
%@%    gmi                   -  Mutual Information estimates for non-equidistantly sampled time series.
%@%    imi                   -  Interpolated Cross-MI estimates for non-equidistantly sampled time series.
%@%    intersectt            -  Intersection of two time strings.
%@%    nest_scatter          -  produces a scatterplot of two irregular time series.
%@%    nexcf                 -  Cross-correlation estimates for non-equidistantly sampled time series.
%@%    pspek                 -  computes the power spectrum for a given autocorrelation function.
%@%    simgram               -  calculates windowed cross similarity between two time series.
%@%    similarity            -  computes statistical associations for irregular time series.
%@%    slidingwindow         -  gives overlapping partition windows for time series.
%@%    tar1                  -  simulates irregular time series of threshold autoregressive processes.
%@%    tauest_ls             -  estimates persistence of irregular time series with least squares.
%@%    test_nest             -  script for NESToolbox V.0.1 delta containing illustration test cases.
%@
%@% Generated at 06-Sep-2013 10:46:07 by MAKEINSTALL
%@% Modified at 18-Mar-2014 18:34:46 by MAKEINSTALL
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __MIhist.m__ -->
%@function Z = MIhist(x, y, nbins)
%@%MIHIST computes the mutual information between variables X and Y.
%@%   Z=MIHIST(X,Y,NBINS) computes the marginal and joint densities for the
%@%   variables X and Y from one (marginal) or two-dimensional (joint)
%@%   histograms with NBINS bins. The formula can be written as a sum over
%@%   joint p(x,y) and marginal densities p(x) and p(y)
%@%   MI(X,Y)=sum p(x,y)*log[ p(x,y)/(p(x)*p(y))],
%@%   or from the entropies MI(X,Y)=Ix+Iy-Ixy,
%@%   as the sum of the individual entropies minus the joint entropy.
%@%   Please refer to Cover&Thomas, 2010 for for the information-theoretic
%@%   background and Rehfeld,2013 for applications.
%@% 
%@%    REFERENCE:
%@%       gMI: K. Rehfeld, N. Marwan, S. Breitenbach, J. Kurths: Late
%@%       Holocene Asian Monsoon Dynamics from small but complex paleoclimate
%@%       networks, 41(1), 3–19 (2013).
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Norbert Marwan, Kira Rehfeld
%@% ver: 1.0  last rev: 2013-08-21
%@%
%@% see also: GMI,IMI,SIMILARITY
%@
%@
%@
%@% normalise the data and replace the values with integers
%@% in the range [1 nbins]
%@x = x - repmat(min(x), size(x,1), 1);
%@y = y - repmat(min(y), size(y,1), 1);
%@x = x ./ repmat(max(x) + eps, size(x,1), 1);
%@y = y ./ repmat(max(y) + eps, size(y,1), 1);
%@    
%@    %multiply by nbins-> frequencies
%@    x = floor(x * nbins) + 1;
%@    y = floor(y * nbins) + 1;
%@    
%@    % compute probabilities
%@    Z = zeros(1,size(x,2));
%@    for i = 1:size(x,2)
%@    
%@        Pxy = accumarray([x(:,i) y(:,i)] + 1, 1);
%@        Px = sum(Pxy,1);
%@        Py = sum(Pxy,2);
%@
%@        Pxy = Pxy / sum(Pxy(:));
%@        Px = Px / sum(Px(:));
%@        Py = Py / sum(Py(:));
%@
%@        % entropies
%@        Ix = -sum((Px(Px ~= 0)) .* log(Px(Px ~= 0)));
%@        Iy = -sum((Py(Py ~= 0)) .* log(Py(Py ~= 0)));
%@        Ixy = -sum(Pxy(Pxy ~= 0) .* log(Pxy(Pxy ~= 0)));                                                                    
%@        
%@        % mutual information
%@        Z(i) = Ix + Iy - Ixy;
%@        
%@    end    
%<-- ASCII ends here -->
%<-- ASCII begins here: __PANlayerchronensemble.m__ -->
%@function [ RecordEnsemble] = PANlayerchronensemble(proxymeas,timeorig,noEns,varargin)
%@%PANLAYERCHRONENSEMBLE creates alternative timescales for layer counted data.
%@%   [Rensemble] = PANlayerchronensemble(proxymeas,timeorig,noEns,'perc',1)
%@%   takes the original timescale timeorig (monotonically increasing regular
%@%   timescale) and the corresponding proxy measurements proxymeas and
%@%   creates noEns alternative timescales in which up to 1 percent counting
%@%   error is introduced. It is assumed here that, while counting, a
%@%   growth layer could have been missed while counting at any point in the
%@%   record, but double-counting may not occur.
%@%
%@%
%@%     EXAMPLE 
%@%     tx=1:500;tx=tx'; % create sample time vector
%@%     x=randn(size(tx));
%@%     % example for percentwise error (increasing with depth)
%@%     NoEns=100;
%@%     [LCrec1] = PANlayerchronensemble(x,tx,NoEns,'%',1);
%@%     % example for absolute counting error error (i.e. +/- 5 years)
%@%     [LCrec] = PANlayerchronensemble(x,tx,NoEns,'abs',1);
%@%     xmat=repmat(x,1,NoEns); % create 
%@%     subplot(2,1,1)
%@%     plot(LCrec(:,2:10),xmat(:,2:10),'k');hold on;
%@%     plot(LCrec1(:,2:10),xmat(:,2:10),'g');xlim([0 500])
%@%     ylabel('pseudoproxy measurement value');
%@%     xlabel('time')
%@%     subplot(2,1,2)
%@%     plot(1:length(tx),var(LCrec(:,2:end),0,2),'k'); hold on;
%@%     plot(1:length(tx),var(LCrec1(:,2:end),0,2),'g');
%@%     xlabel('time');ylabel('variance (time)');xlim([0 500])
%@%     legend('percentwise age error increase','absolute layer counting error')
%@%    
%@%
%@%    REFERENCE:
%@%
%@%       K. Rehfeld and J. Kurths: Similarity measures for irregular
%@%       and age uncertain time series, submitted to Clim. Past. Discuss,
%@%       2013.
%@%
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 0.9  last rev: 2013-08-21
%@
%@% Base data
%@TE=max(timeorig);
%@T0=min(timeorig);
%@TimeSp=TE-T0;
%@Nobs=length(proxymeas);
%@
%@
%@% find what kind of uncertainty is given
%@for i=1:2:length(varargin)-1    
%@    switch varargin{i}
%@        case {'layer', 'perc', 'p', '%'}
%@            % percentwise
%@            mxerr=round(TimeSp*varargin{i+1}/100);
%@            if mxerr<1,
%@                error('too small an error');
%@            end
%@            T0=repmat(min(timeorig),noEns,1);
%@        case {'abs', 'a', 'layered, abs'}
%@            % absolute error
%@            mxerr=varargin{i+1};
%@            % starting year could be wrong by pm 5 years
%@            T0err=randi(2*mxerr+1,noEns,1)-1-mxerr;
%@            T0=T0err+T0;
%@        otherwise % layer
%@            mxerr=TimeSp*varargin{i+1}/100;
%@    end
%@end
%@
%@% The error could be 0, but it could be also anything up to mxerr. --> For
%@% each of the noEns, find a (random) maximum 
%@E=randi(mxerr+1,noEns,1);% 0 is not treated by randi=> subtract to allow "perfect" chron
%@E=E-1;
%@
%@
%@for k=1:noEns,
%@    % devise timescale
%@    tries=0;
%@    compute=1;
%@    while compute
%@        tries=tries+1;
%@        tsk=T0(k):1:(TE+E(k)+(T0(k)-min(timeorig))); % potential times span from T0 to Tmax+Error (T0 already incorporates error)
%@        Nexcess=length(tsk)-Nobs;
%@        % assuming that each observation really represents one year (i.e., no
%@        % year was double-counted!), the time scale error results in "excess
%@        % years" that have to be removed. Each year is just as likely to NOT
%@        % have been recorded as any other.
%@        
%@        for ii=1:Nexcess,
%@            remind=randi(length(tsk),1,1);
%@            tsk(remind)=[];
%@        end
%@        % what to do if the time scale is too short?  => there has to be a year
%@        % for each observation!
%@        
%@        if length(tsk)==length(proxymeas) % ok, everything workt
%@            T(:,k)=tsk';
%@            clear tsk Nexcess
%@            compute=0; % exit while loop for this timescale iteration
%@        else % try again
%@            if tries==100,
%@                error('something is wrong here')
%@            else
%@                compute=1; % continue to try to get a time scale
%@                continue
%@                
%@            end
%@        end
%@    end
%@    
%@end
%@% plot(T)
%@% hold on;
%@% plot(timeorig,'k','LineWidth', 3)
%@
%@RecordEnsemble=[proxymeas T];
%@end
%<-- ASCII ends here -->
%<-- ASCII begins here: __README.txt__ -->
%@----------------------------------------------------------------------
%@NESToolbox
%@Copyright (c) 2014, Kira Rehfeld
%@http://www.tocsy.pik-potsdam.de/nest.php
%@
%@Authors: Kira Rehfeld <kira.rehfeld@awi.de>, with ideas and contributions
%@	 from Norbert Marwan and ideas from Jobst Heitzig, Bedartha Goswami
%@	 and Sebastian Breitenbach.
%@Version: 1.1 (18-03-2014)
%@----------------------------------------------------------------------
%@
%@
%@The NESToolbox is a collection of algorithms to perform similarity
%@estimation for  irregularly sampled time series as they arise for example
%@in the geosciences. It is implemented as a toolbox for the widely used
%@software MATLAB (see http://www.mathworks.com) and the freely available 
%@open-source software OCTAVE (see
%@http://www.http://www.gnu.org/software/octave/). 
%@
%@At the center of the toolbox are the functions for linear and nonlinear
%@similarity estimation for irregular time series, which are based on 
%@Gaussian-kernel weight functions. Cross-correlation estimation, as it is 
%@used in most standard time series analysis, cannot be performed for 
%@irregular time series directly, as the vector index cannot be used as a 
%@substitute for time differences. The conventional approach, interpolation 
%@of the time series to regular grid followed by the use of standard 
%@estimators, has bias side effects [1]. The function similarity.m makes 
%@alternative approaches such as the Gaussian-kernel-based cross correlation 
%@[1], the nonlinear Gaussian-kernel-based mutual information [2] or the 
%@Event Synchronization function [3] available under a single, unified 
%@syntax. Additional functions allow for the estimation of uncorrelated time 
%@series surrogates to test the significance of similarity estimates as well 
%@as nonlinear trends, power spectra and weighted scatterplots. The usage of 
%@these functions is illustrated in the script TEST_NEST.m, and a list of the 
%@available functions can be found in Contents.m.
%@ 
%@To install the toolbox, call the downloaded script install.m from the 
%@Matlab or Octave commandline. The toolbox will be installed in the user's 
%@toolbox folder. In Matlab this is usually a folder called matlab in the 
%@users home directory). In Octave, the standard folder in the home directory 
%@is called octave, and please execute the command savepath afterwards, in 
%@order to ensure that the installed folder will be set permanently. Then 
%@open the script TEST_NEST.m for usage examples and tests.
%@
%@
%@-----------------------------------------------------------------------
%@REQUIREMENTS
%@
%@The toolbox was tested using Matlab R2011b and R2012b with and without the 
%@optimization toolbox present. If the Optimization toolbox is not available, 
%@the parameter estimation in the function TAUEST_LS.m cannot be performed 
%@using nonlinear least squares. The function AR1Sur.m therefore computes the
%@lag-1 autocorrelation using the Gaussian-kernel based estimator instead. 
%@This also the case in Octave 3.2.4, where the toolbox also requires the 
%@additional, freely available, package octave-signal.
%@
%@
%@----------------------------------------------------------------------
%@REFERENCES
%@
%@        [1] K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@        correlation analysis techniques for irregularly sampled time series,
%@        Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@        [2] K. Rehfeld, N. Marwan, S. Breitenbach, J. Kurths: Late
%@        Holocene Asian Monsoon Dynamics from small but complex paleoclimate
%@        networks, Climate Dynamics, 41(1), 3–19 (2013).
%@        [3] K. Rehfeld and J. Kurths: Similarity estimators for irregular 
%@	and age-uncertain time series, Clim. Past, 10, 107-122 (2014).
%@
%@----------------------------------------------------------------------
%@LICENSE
%@
%@The NESToolbox is free software; you can redistribute it and/or
%@modify it under the terms of the GNU General Public License,
%@version 3, as published by the Free Software Foundation.
%@
%@This program is distributed in the hope that it will be useful, but
%@without any warranty; without even the implied warranty of
%@merchantability or fitness for a particular purpose.  See the GNU
%@General Public License for more details.
%@
%@A copy of the GNU General Public License, version 3, is available at
%@http://www.r-project.org/Licenses/GPL-3
%@----------------------------------------------------------------------
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __ar1sur.m__ -->
%@function [Xsur]=ar1sur(tx,x,Nsur)
%@% AR1SUR computes autocorrelated irregular surrogate time series.
%@% [Xsur]=AR1SUR(tx,x,Nsur) calculates a matrix of NSur surrogate time series, Xsur, that fits the
%@% time axis of the  Sx with time axis Tx of a given Vector X(tx). 
%@% The persistence time (or: autocorrelation) is estimated from a least
%@% squares fit of AR1-processes to the data (standard) or, if the
%@% Optimization toolbox is not available, using the autocorrelation function
%@% of the data (e.g. in Octave).
%@% 
%@%
%@%    REFERENCE:
%@%
%@%       ESF: K. Rehfeld and J. Kurths: Similarity measures for irregular
%@%       and age uncertain time series, submitted to Clim. Past. Discuss,
%@%       2013.
%@%
%@% see also: NEXCF,GMI,SIMILARITY,CAR,TAR1
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% (c) Alfred Wegener Institute for Polar and Marine Research 2014
%@% Kira Rehfeld
%@% ver: 1.1  last rev: 2014-01-25
%@
%@acf=0;
%@
%@if length(tx) ~= length(x); error('tx and x not same length');
%@end;
%@%dtau=median(diff(tx));
%@if exist('lsqnonlin','file') % use least-squares estimator if available
%@[Tau]=tauest_ls(tx,x,1);
%@else
%@    % use AR1-coefficient from autocorrelation function
%@    acf=true;
%@    Tau=NaN;
%@end
%@
%@
%@if Tau<min(diff(tx))||~isreal(Tau)||acf==true,
%@    % use ACF coefficients to estimate persistence
%@    lag=[1:2]*mean(diff(tx));
%@    
%@     Cxx=similarity(tx,x,tx,x,lag,'gXCF');
%@     if any(isnan(Cxx))% choose smaller decorrelation scale
%@     	Cxx=similarity(tx,x,tx,x,1:2,'gXCF');
%@     	end
%@     	
%@     	
%@     EstTau=-[1,2].*mean(diff(lag))./log(Cxx([1,2]))';
%@     Tau=mean(real(EstTau));
%@     est_phi=exp(-mean(diff(tx))/(Tau));
%@
%@end
%@
%@Xsur=zeros(length(tx),Nsur);
%@%stream=RandStream('mlfg6331_64','RandnAlg','Polar');
%@for k=1:Nsur
%@    x=randn(size(tx));
%@    % adjust standard deviation of the noise such that the
%@    % standard deviation of the resulting time series is 1:
%@    Eps=randn(length(tx),1);
%@    Eps=(Eps-mean(Eps))./std(Eps);
%@    Xsuri=Eps;
%@    % possible correction
%@    %%KK=sqrt((1-exp(-2*diff(tx)./TauKern))./(2./TauKern));
%@    KK=1;
%@    for i=2:length(tx)
%@        Xsuri(i)=exp(-(tx(i)-tx(i-1))./Tau).*Xsuri(i-1) +KK.*Eps(i);
%@    end
%@    
%@    Xsuri=(Xsuri-mean(Xsuri))./std(Xsuri) *std(x);
%@    
%@    Xsur(:,k)=Xsuri;
%@    clear Xsuri
%@    
%@end
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __car.m__ -->
%@function [x y]=car(tx,ty, phi, coupl_strength, lag, interp_method)
%@% CAR gives irregular time series for coupled autoregressive (AR1) processes.
%@%   [X,Y]=CAR(tx,ty,phi,coupl_strength,lag,interp_method) returns observations X and
%@%   Y for the observation time vectors tx and ty (monotonically increasing
%@%   column vectors). Phi (double value between 0 and 1) determines the
%@%   autocorrelation in X, coupl_strength is the coupling strength (double value
%@%   between 0 and 1 and lag (integer value) determines the lag at which the
%@%   time series are coupled. The autocorrelated time series are first
%@%   generated on a regular timescale at very high resolution and
%@%   subsequently re-sampled by interpolation with the method interp_method
%@%   ('linear','spline') to the timescales tx and ty.
%@%   
%@%   EXAMPLE:
%@%   % Definition of the processes (see Rehfeld et al., 2011 for more
%@%   % details)
%@%   % X_n=exp((t_n-t_n-1)/Tau *X_n-1+Eps_n=phi*X_{n-1}+Eps_n;
%@%   % Y_n=coupl_strength*X_(n-lag)+(1-coupl_strength)*Eps2_n.
%@%   % make sample time scales, here: regular timescales for tx and ty
%@%     tx=1:1:100;tx=tx';ty=tx;   
%@%     phi=0.1; 
%@%     % AR1 coefficient (here equivalent to true Cross-Correl. at
%@%     %lag of coupling
%@%     coupl_strength=0.9; % coupling strength
%@%     lag=2; % expected to be positive for causality, X always drives Y
%@%     interp_method='linear'
%@%     [x y]=car(tx,ty, phi, coupl_strength, lag, interp_method)
%@%
%@%
%@%    REFERENCE:
%@%       K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@%       correlation analysis techniques for irregularly sampled time series,
%@%       Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@%
%@%
%@% see also: NEXCF,GMI,EVENTSYNCHRO,SIMILARITY,TAR
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 1.0  last rev: 2013-08-22
%@
%@
%@% set random stream according to clock - does not yet work in Octave
%@%RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(1000*clock)))
%@
%@% set the resolution at which the AR1process is simulated at high
%@% resolution as 10 times the actual mean resolution
%@dtsur=min(mean(diff(tx)), mean(diff(ty)))./10;
%@% find beginning and end points and create regular timescale between them
%@t1=min(min(tx), min(ty))-lag; % subtract lag  to compensate transient effects;
%@t2=max(max(tx),max(ty));
%@t=(t1-lag):dtsur:t2;t=t';
%@% find persistence time with relation to the actual processes
%@fi=exp(-dtsur/(-1/log(phi)));
%@
%@%% First step: Create the driving process at high resolution (X)
%@Eps1=randn(length(t),1); % innovation term
%@
%@xsur=Eps1;% initialization
%@for i=2:length(t)
%@    xsur(i)=fi.*xsur(i-1) +sqrt(1-fi^2)*Eps1(i);
%@end
%@
%@xsur=(xsur-mean(xsur))./std(xsur);
%@
%@%% Second step: Couple the second timeseries
%@Eps2=randn(length(t),1);
%@ysur=Eps2;
%@
%@% determine appropriate lag of coupling for ysur/xsur
%@L=ceil(lag./dtsur);
%@% ysur is coupled to xsur
%@ysur(L+1:end)=coupl_strength*xsur(1:end-L)+sqrt(1-coupl_strength^2)*Eps2(L+1:end);
%@% normalize
%@ysur=(ysur-mean(ysur))./std(ysur);
%@
%@%% Third step, interpolate to (irregular) desired time scale
%@
%@y = interp1(t,ysur,ty,interp_method,'extrap');
%@x = interp1(t,xsur,tx,interp_method,'extrap');
%@y=(y-mean(y))./std(y);
%@x=(x-mean(x))./std(x);
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __deinstall.m__ -->
%@function deinstall
%@%DEINSTALL   Removes NESTool.
%@%    DEINSTALL removes all files of NESTool from
%@%    the filesystem and its entry from the Matlab
%@%    startup file.
%@%    
%@%    This installation script was generated by using 
%@%    the MAKEINSTALL tool. For further information
%@%    visit http://matlab.pucicu.de
%@
%@% Copyright (c) 2008-2012
%@% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
%@% http://www.pik-potsdam.de
%@%
%@% Copyright (c) 2002-2008
%@% Norbert Marwan, Potsdam University, Germany
%@% http://www.agnld.uni-potsdam.de
%@%
%@% Generation date: 18-Mar-2014 18:34:46
%@% $Date: 2014/03/12 16:17:50 $
%@% $Revision: 3.30 $
%@
%@error(nargchk(0,0,nargin));
%@
%@try
%@  if isoctave
%@      more off
%@  end
%@  fid = 0;
%@  warning('off')
%@  disp('------------------------')
%@  disp('    REMOVING NESTool    ')
%@  disp('------------------------')
%@  currentpath=pwd;
%@  oldtoolboxpath = fileparts(which(mfilename));
%@
%@  disp(['  NESTool found in ', oldtoolboxpath,''])
%@  i = input('> Delete NESTool? Y/N [Y]: ','s');
%@  if isempty(i), i = 'Y'; end
%@
%@  if strcmpi('Y',i)
%@%%%%%%% check for entries in startup
%@  
%@        p=path; i1=0; i = ''; number_warnings_pathdef = 0;
%@  
%@        while findstr(upper('NESTool'),upper(p)) > i1
%@           i1=findstr(upper('NESTool'),upper(p));
%@           if ~isempty(i1)
%@               i1=i1(end);
%@               if isunix, i2=findstr(':',p); else, i2=findstr(';',p); end
%@               i3=i2(i2>i1);                 % last index pathname
%@               if ~isempty(i3), i3=i3(1)-1; else, i3=length(p); end
%@               i4=i2(i2<i1);                 % first index pathname
%@               if ~isempty(i4), i4=i4(end)+1; else, i4=1; end
%@               rmtoolboxpath=p(i4:i3);
%@%%%%%%% removing entry in startup-file
%@               rmpath(rmtoolboxpath)
%@               err = savepath;
%@               if number_warnings_pathdef == 0 && err, disp('  ** Warning: No write access to pathdef.m file!'), number_warnings_pathdef = number_warnings_pathdef+1; end
%@               if i4>1, p(i4-1:i3)=''; else, p(i4:i3)=''; end
%@               startup_exist = exist('startup','file');
%@               if isoctave startup_exist = exist(fullfile('~','.octaverc'),'file'); end
%@               if startup_exist
%@                    startupfile=which('startup');
%@                    startuppath=startupfile(1:findstr('startup.m',startupfile)-1);
%@                    if isoctave
%@                        startuppath = ['~',filesep];
%@                        startupfile = fullfile('~','.octaverc');
%@                    end
%@                    fid = fopen(startupfile,'r');
%@                    k = 1;
%@                    while 1
%@                       tmp = fgetl(fid);
%@                       if ~ischar(tmp), break, end
%@                       instpaths{k} = tmp;
%@                       k = k + 1;
%@                    end
%@                    k=1;
%@                    while k <= length(instpaths)
%@                        if ~isempty(findstr(rmtoolboxpath,instpaths{k}))
%@                            disp(['  Removing startup entry ', instpaths{k}])
%@                            instpaths(k)=[];
%@                        end
%@                        k=k+1;
%@                    end
%@                    fid=fopen(startupfile,'w');
%@                    for i2=1:length(instpaths), 
%@                        fprintf(fid,'%s\n', char(instpaths{i2})); 
%@                    end
%@                    fclose(fid);
%@               end
%@           end
%@           p = path; i1 = 0;
%@       end
%@%%%%%%% removing old paths
%@        if exist(oldtoolboxpath,'dir') == 7
%@           if isoctave, confirm_recursive_rmdir (false, 'local'); end
%@           disp(['  Removing files in ',oldtoolboxpath,''])
%@           cd(oldtoolboxpath)
%@           dirnames='';filenames='';
%@           temp='.:';
%@           while ~isempty(temp)
%@               [temp1 temp]=strtok(temp,':');
%@               if ~isempty(temp1)
%@                   dirnames=[dirnames; {temp1}];
%@                   x2=dir(temp1);
%@                   for i=1:length(x2)
%@                       if ~x2(i).isdir, filenames=[filenames; {[temp1,'/', x2(i).name]}]; end
%@         	             if x2(i).isdir && ~strcmp(x2(i).name,'.') && ~strcmp(x2(i).name,'..'), temp=[temp,temp1,filesep,x2(i).name,':']; end
%@                   end
%@               end
%@           end
%@           dirnames = strrep(dirnames,['.',filesep],'');
%@           dirnames(strcmpi('.',dirnames)) = [];
%@           l = zeros(length(dirnames),1); for i=1:length(dirnames),l(i)=length(dirnames{i}); end
%@           [i i4]=sort(l); i4 = i4(:);
%@           dirnames=dirnames(flipud(i4));
%@           for i=1:length(dirnames)
%@              delete([dirnames{i}, filesep,'*'])
%@              if exist('rmdir') == 5 && exist(dirnames{i}) == 7, rmdir(dirnames{i},'s'); else, delete(dirnames{i}), end
%@              disp(['  Removing files in ',char(dirnames{i}),''])
%@           end
%@           if exist(currentpath), cd(currentpath), else, cd .., end
%@           if strcmpi(currentpath,oldtoolboxpath), cd .., end
%@           if exist('rmdir') == 5 && exist(oldtoolboxpath) == 7, rmdir(oldtoolboxpath,'s'); else, delete(oldtoolboxpath), end
%@           disp(['  Removing folder ',oldtoolboxpath,''])
%@        end
%@       disp(['  NESTool now removed.'])
%@  else
%@       disp(['  Nothing happened. Keep smiling.'])
%@  end
%@  tx=version; tx=strtok(tx,'.'); if str2double(tx)>=6 && exist('rehash','builtin'), rehash, end
%@  warning on
%@  if isoctave
%@      more on
%@  end
%@  if exist(currentpath,'dir') ~= 7, cd(fileparts(currentpath)), else, cd(currentpath), end
%@  
%@%%%%%%% error handling
%@
%@catch
%@  x=lasterr;y=lastwarn;
%@  if ~strcmpi(lasterr,'Interrupt')
%@    if fid>-1, 
%@      try, z=ferror(fid); catch, z='No error in the installation I/O process.'; end
%@    else
%@      z='File not found.'; 
%@    end
%@    fid=fopen('deinstall.log','w');
%@    fprintf(fid,'%s\n','A critical error has occurred. Please inform the distributor');
%@    fprintf(fid,'%s\n','of the toolbox, where the error occured and send us the entire');
%@    fprintf(fid,'%s\n','screen output of the installation, the following error');
%@    fprintf(fid,'%s\n','report, and the informations about the toolbox (distributor,');
%@    fprintf(fid,'%s\n','name, URL etc.). Provide a brief description of what you were');
%@    fprintf(fid,'%s\n','doing when this problem occurred.');
%@    fprintf(fid,'%s\n','E-mail or FAX this information to us at:');
%@    fprintf(fid,'%s\n','    E-mail:  marwan@pik-potsdam.de');
%@    fprintf(fid,'%s\n','       Fax:  ++49 +331 288 2640');
%@    fprintf(fid,'%s\n\n\n','Thank you for your assistance.');
%@    fprintf(fid,'%s\n',repmat('-',50,1));
%@    fprintf(fid,'%s\n',datestr(now,0));
%@    fprintf(fid,'%s\n',['Matlab ',char(version),' on ',computer]);
%@    fprintf(fid,'%s\n',repmat('-',50,1));
%@    fprintf(fid,'%s\n','NESTool');
%@    fprintf(fid,'%s\n',x);
%@    fprintf(fid,'%s\n',y);
%@    fprintf(fid,'%s\n',z);
%@    fclose(fid);
%@    disp('----------------------------');
%@    disp('       ERROR OCCURED ');
%@    disp('   during deinstallation');
%@    disp('----------------------------');
%@    disp(x);
%@    disp(z);
%@    disp('----------------------------');
%@    disp('   A critical error has occurred. Please inform the distributor');
%@    disp('   of the toolbox, where the error occured and send us the entire');
%@    disp('   screen output of the installation, the error report report');
%@    disp('   and the informations about the toolbox (distributor, name,');
%@    disp('   URL etc.). For your convenience, this information has been')
%@    disp('   recorded in: ')
%@    disp(['   ',fullfile(pwd,'deinstall.log')]), disp(' ')
%@    disp('   Provide a brief description of what you were doing when ')
%@    disp('   this problem occurred.'), disp(' ')
%@    disp('   E-mail or FAX this information to us at:')
%@    disp('       E-mail:  marwan@pik-potsdam.de')
%@    disp('          Fax:  ++49 +331 288 2640'), disp(' ')
%@    disp('   Thank you for your assistance.')
%@  end
%@  warning('on')
%@  if exist(currentpath,'dir') == 7, cd(fileparts(currentpath)), else, cd(currentpath), end
%@  if isoctave
%@      more on
%@  end
%@end
%@
%@function flag = isoctave
%@% ISOCTAVE   Checks whether the code is running in Octave
%@%   ISOCTAVE is returning the value TRUE if executed within the
%@%   Octave environment, else it is returning FALSE (e.g. when
%@%   called within Matlab.
%@
%@a = ver('Octave');
%@
%@if ~isempty(a) && strfind(a(1).Name,'Octave')
%@    flag = true;
%@else
%@    flag = false;
%@end
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __eventsynchro.m__ -->
%@function [Q q,varargout] = eventsynchro(tx,x,ty,y,lag,conflevel)
%@%EVENTSYNCHRO calculates event synchronization for irregular time series.
%@%
%@%   Event synchronization as described in Quiroga et al, 2002, assesses the
%@%   co-occurence of events in time series {tx,x} and {ty,y}. Events, in our
%@%   case, are defined as observation falling onto the margins of the
%@%   distributions, beyond the confidence levels. Q gives a measure of the
%@%   coupling strength, q a sense of the direction.
%@%
%@%       EXAMPLE
%@%       tx=[1:1000]';ty=tx;x=randn(size(tx));y=randn(size(ty));
%@%       subplot(3,1,1)
%@%       plot(x,y,'*')
%@%       y(x>0)=x(x>0);
%@%       hold on;plot(x,y,'r*');
%@%       xlabel('x'),ylabel('y');legend('uncorrelated', 'event-correlated')
%@%       conflevel=0.9; % events are below the 0.05 and .95fth level;
%@%       lag=0;% absolute lag in time
%@%       [Q,q]=eventsynchro(tx,x,ty,y,lag,0.9)
%@%       % evsx and evsy give the event times in x/y
%@%       [Q q evsx evsy] = eventsynchro(tx,x,ty,y,lag,conflevel);
%@%       
%@%       subplot(3,1,2)
%@%       plot(tx,x);hold on;plot(ty,y,'r');legend('x','y');xlabel('t'),ylabel('variable')
%@%       subplot(3,1,3)
%@%       stem(tx(evsx),1.05*ones(size(tx(evsx)))); hold on; stem(tx(evsx),0.95*ones(size(tx(evsx))),'r');
%@%
%@%       % and tau the temporal threshold to decide on coincidence/non-coincidence:
%@%       [Q q evsx evsy tau] = eventsynchro(tx,x,ty,y,lag,conflevel);
%@%
%@%    REFERENCE:
%@%
%@%       ESF: K. Rehfeld and J. Kurths: Similarity measures for irregular
%@%       and age uncertain time series, submitted to Clim. Past. Discuss,
%@%       2013.
%@%
%@% see also: NEXCF,GMI,SIMILARITY
%@
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 1.0  last rev: 2013-08-21
%@
%@
%@clo=(1-conflevel)/2;
%@chi=conflevel+clo;
%@evsx=find((x>=quantile(x,chi))|(x<quantile(x,clo)));
%@evsy=find((y>=quantile(y,chi))|(y<quantile(y,clo)));
%@
%@
%@for k=1:length(lag)
%@    
%@    txs=tx-lag(k);
%@    % fixed tau, not as in Quiroga, 2002
%@    tau=max(mean(diff(txs)),floor(abs(min(min(diff(txs(evsx))),min(diff(ty(evsy))))/2)));
%@    
%@    for i=1:length(evsx)
%@        for j=1:length(evsy)
%@            TD(i,j)=txs(evsx(i))-ty(evsy(j));
%@        end
%@    end
%@    
%@    %% Now we threshold the time difference matrix to get the coincidence matrices J
%@    TDx=zeros(size(TD));TDy=zeros(size(TD));
%@    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@    % !!!!!!!!!!!!!!!
%@    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@    
%@    % find events that occur in x before they occur in y
%@    TDx=double((TD>0)&(TD<tau));%Jxy=1 <=> 0<tdx<tau % Malik, 2011
%@    %TDx=double((TD>0)&(TD<=tau));%Jxy=1 <=> 0<tdx<=tau, as in Quiroga,2002
%@    % alternative definition
%@    TDy=double((-TD>0)&(-TD<tau));%Jxy=1 <=> 0<tdx<tau % Malik, 2011
%@    %TDy=double((-TD>0)&(-TD<=tau));%Jxy=1 <=> 0<tdx<=tau,  as in Quiroga,2002
%@    
%@    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@    % find those that occur simultaneously & set the matrix value to 1/2
%@    ind=find(TD==0);
%@    TDx(ind)=1/2;
%@    TDy(ind)=1/2;
%@        
%@    % to get c_xy and c_yx, just sum the matrices
%@    ctx=sum(sum(TDx));
%@    cty=sum(sum(TDy));
%@    
%@    %% Now calculate Q and q, the coupling strength and direction estimates
%@    
%@    denominator=sqrt(length(evsx)*length(evsy));
%@    
%@    Q=(ctx+cty)/denominator;
%@    q=(cty-ctx)/denominator;
%@    
%@    
%@end
%@
%@varargout{1}=evsx;
%@varargout{2}=evsy;
%@varargout{3}=tau;
%@
%@end
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __extractwindow.m__ -->
%@function [txw,xw] = extractwindow(tx,x,links,rechts,warng)
%@%EXTRACTWINDOW cuts (irregular) time series according to the boundaries.
%@%   [txw,xw] = extractwindow(tx,x,links,rechts,warng) cuts the time series
%@%   {tx,x} so that all observation times lie between the boundaries of
%@%   [links rechts], where these are the left (younger) and right (older)
%@%   limits. It is especially useful when dealing with irregular time series,
%@%   where the index bears no indication of the time.
%@%   EXAMPLE
%@%       % We take two time series from stochastic processes X and Y, {tx,x} and
%@%       % {ty,y}, which are defined as follows: 
%@%       % Sampling times 
%@%       tx = (1:1001)'; ty = (1:1001)'; 
%@%       % Original signals 
%@%       x = randn(1001,1);
%@%       x(2:end) = 0.5*x(1:end-1)+randn(1000,1); 
%@%       y = randn(1001,1);y_nl=y;
%@%       y(2:end) = 0.7*x(1:end-1)+randn(1000,1);
%@%       % We have, however, in reality, observed only 80% of the points in y, 
%@%       % yielding a time series {ty_miss,y_miss}:
%@%       rem = randi(length(ty), floor(length(ty)*0.2),1); 
%@%       ty_miss = ty; ty_miss(rem) = [];y_miss = y; y_miss(rem) = [];
%@%
%@%       % the time series are not equally long, but we want to compute a
%@%       % statistic for the overlapping portions
%@%       t1=min(min(tx),min(ty_miss));tEnd=max(max(tx),max(ty_miss));
%@%       [txw,xw]=extractwindow(tx,x,t1,tEnd,0)
%@%       [tyw,yw]=extractwindow(ty_miss,y_miss,t1,tEnd,0)
%@%
%@%       % we can now, for example, compute the Gaussian-kernel based cross
%@%       % correlation for the two time series:
%@%       lag=[-25:1:25]
%@%       Cxy = similarity(tx, x, ty_miss, y_miss,lag, 'gXCF');
%@%       plot(lag,Cxy), xlabel('Lag'), ylabel('Cross-correlation C_{xy}')
%@%
%@% see also: SIMILARITY,SLIDINGWINDOW
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 1.0  last rev: 2013-11-09
%@
%@if isempty(warng),
%@warng=1;
%@end
%@
%@txw=[];xw=[];
%@% length of the time series x and y
%@Nx = length(x);
%@
%@if (Nx ~= size(tx))
%@   error('Vectors for sampling times and observations have to be of same length.'), 
%@end
%@
%@if (rechts<=links)
%@    error('Borders are not correctly specified')
%@end
%@
%@% if there are no observations in the given time range
%@if  (rechts<=min(tx))|| (links>=max(tx)),
%@    if warng,
%@        warning('empty time series slot')
%@    end
%@    txw=[]; %return an empty time series
%@    xw=[];
%@else % obs in range
%@    mx=find(tx>links& tx<=rechts);
%@    txw=tx(mx);
%@    xw=x(mx,:);
%@end
%@
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@
%@
%@end
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __gausssmooth.m__ -->
%@function [Xsmooth,W] = gausssmooth(tx,x,h,varargin)
%@%GAUSSSMOOTH smoothes irregularly sampled time series with gaussian weights.
%@%   [Xsmooth] = GAUSSSMOOTH(tx,x,h) computes a smoothed version of
%@%   the time series observations x, where the amount of smoothing is
%@%   determined by the kernel width h. The first/last points should be taken with
%@%   caution as they are only smoothed onesidedly.
%@%
%@%   EXAMPLE:
%@%   tx = (1:1001)';
%@%   % Original signal is random with a trend:
%@%   x = randn(1001,1)+linspace(0,2,length(tx))';
%@%   plot(tx,x);xlabel('time index');xlim([0 1000]);
%@%   % use GAUSSSMOOTH to obtain the trend (low-pass filter):
%@%   x_trend=gausssmooth(tx,x,50)
%@%   hold on; plot(tx,x_trend,'r')
%@%   % or use it to extract the high-frequency variability by removing the
%@%   % trend component (high-pass filter)
%@%   x_fluct=x-x_trend;
%@%   plot(tx,x_fluct,'g')
%@%   legend('orig','trend','high-frequency')
%@%
%@%
%@%    REFERENCE:
%@%       K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@%       correlation analysis techniques for irregularly sampled time series,
%@%       Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% (c) Alfred Wegener Institute for Polar and Marine Research 2014-
%@% 
%@% Kira Rehfeld
%@% ver: 1.1  last rev: 2014-01-25 (vectorization)
%@
%@% varargin -> cut off
%@
%@
%@Sx=size(x);St=size(tx);
%@% match the sizes! Sx(1) should be equal to length(tx), i.e. St(1)
%@if ~(Sx(1)==St(1))
%@	error('Gausssmooth: time and observation vector differ in length!')
%@	end
%@	
%@if h==0, % if smoothing = zero, return original series
%@	Xsmooth=x;
%@	W=[];
%@else	
%@
%@W=zeros(length(tx));% make sure the matrices come columnwise!
%@
%@% computing the weights 
%@% row i should contain the weights for time point i in tx/x
%@for k=1:length(tx)
%@	%k-th row=kth observation
%@	W(:,k)=1/sqrt(2*pi*h^2).*(exp(-((tx(:)-tx(k)).^2)./(2*h^2)));%
%@end
%@
%@% matrix of output time series (columnwise)
%@weightsum=sum(W,2);
%@Xsmooth=W*x./repmat(weightsum,1,Sx(2));
%@	
%@ end
%@ 
%@end
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __generate_t.m__ -->
%@function [tx] = generate_t(dt,tmin,tmax,method,varargin)
%@%GENERATE_T generates arbitrary time axes with specified properties.
%@%     [tx] = GENERATE_T(dt,tmin,tmax,method) generates a timescale vector tx
%@%     from tmin to tmax with a mean time step of dt and according to the
%@%     method indicated.
%@%     The generated monotonically increasing vectors obey tmin<=tx<=tmax and
%@%     have a target mean time resolution of dt. The methods chosen may need
%@%     additional inputs, as indicated in the examples.
%@%     Standard input:
%@%     dt - desired sampling interval tmin, tmax, min and max of the new vector
%@%     time points method: any of 'linear', 'gamma', 'quasireg', 'poisson',
%@%     'trapez','drift', 'gdrift', 'sindrift', 'gskew'.
%@%     The first optional parameter, necessary for the irregular sampling
%@%     distributions gives the quality of the sampling, e.g the percentage of
%@%     missing points or the skewness of the distribution.
%@% 
%@%     Further optional arguments:
%@%     'mean' -> standardize sampling to dt
%@%     'NoNorm' -> explicitly do not standardize sampling (default)
%@% 
%@%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@%     %%
%@%     %Examples:
%@%     L=100; % length of the time series
%@%     dt=1; % average sampling rate of time series
%@%     %%% Gamma distributed sampling intervals
%@%     ty=generate_t(dt,1,L,'gamma', 0.1, 'mean');
%@%     subplot(3,2,1)
%@%     hist(diff(ty)); %distribution of time series sampling points
%@%     title('gamma distributed')
%@%     subplot(3,2,2)
%@%     plot(diff(ty))
%@%     title('sampling interval')
%@%     %%% missing values
%@%     ty=generate_t(dt,1,L,'missing', 0.9);
%@%     subplot(3,2,3)
%@%     hist(diff(ty)); %distribution of time series sampling points
%@%     title('missing values')
%@%     subplot(3,2,4)
%@%     plot(diff(ty))
%@%     title('sampling interval')
%@% 
%@%     %%% drifting sampling rate
%@%     ty=generate_t(dt,1,L,'gdrift', 0.1);
%@%     subplot(3,2,5)
%@%     hist(diff(ty)); %distribution of time series sampling points
%@%     title('gamma distributed, drifting mean sampling interval')
%@%     subplot(3,2,6)
%@%     plot(diff(ty))
%@%     title('sampling rate')
%@% 
%@%     %%
%@%     figure()
%@%     ty=generate_t(dt,1,L,'gdrift', 0.1, 'mean');
%@%     subplot(3,2,1)
%@%     hist(diff(ty)); %distribution of time series sampling points
%@%     title('gamma distributed, drifting mean')
%@%     subplot(3,2,2)
%@%     plot(diff(ty))
%@%     title('sampling interval')
%@%     %%%
%@%     ty=generate_t(dt,1,L,'sindrift', 10);
%@%     subplot(3,2,3)
%@%     hist(diff(ty)); %distribution of time series sampling points
%@%     title('gamma distributed, sinusoidally modulated')
%@%     subplot(3,2,4)
%@%     plot(diff(ty))
%@%     title('sampling interval')
%@%     %%%
%@%     ty=generate_t(dt,1,L,'drift', 0.1);
%@%     subplot(3,2,5)
%@%     hist(diff(ty)); %distribution of time series sampling points
%@%     title('drifting mean sampling interval')
%@%     subplot(3,2,6)
%@%     plot(diff(ty))
%@%     title('sampling rate')
%@
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@%    REFERENCE:
%@%       gXCF: K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@%       correlation analysis techniques for irregularly sampled time series,
%@%       Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@%       gMI: K. Rehfeld, N. Marwan, S. Breitenbach, J. Kurths: Late
%@%       Holocene Asian Monsoon Dynamics from small but complex paleoclimate
%@%       networks, 41(1), 3–19 (2013).
%@%
%@%
%@% see also: NEXCF,GMI,EVENTSYNCHRO
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 0.9  last rev: 2013-08-22
%@
%@%% Varargin distribution
%@ToMean=false;
%@
%@optargin = size(varargin,2);
%@stdargin = nargin - optargin;
%@if stdargin <4
%@    error('Too few arguments');
%@end
%@
%@%% Check of inputs
%@GM=any([ (strcmpi('gamma',method)) strcmpi('gdrift', method)]); %gamma
%@ToMean=any(strcmpi(varargin, 'mean'));
%@NoNorm=any(strcmpi(varargin, 'NoNorm'));
%@
%@if ~ToMean, NoNorm=true; end;
%@if  GM&&(isempty(varargin)||~isa(varargin{1}, 'numeric')),
%@    shape=1;
%@    scale=dt/shape;
%@end
%@
%@if GM&&(NoNorm&&isa(varargin{1}, 'numeric')),
%@    meanv=dt;
%@    shape=1/varargin{1};
%@    scale=meanv/shape;
%@end
%@
%@if optargin>2&&GM&&isa(varargin{1},'numeric')&&isa(varargin{2}, 'numeric'),
%@    
%@    shape=varargin{1};
%@    scale=varargin{2};
%@    
%@    sprintf('gamma %g-%g', shape,scale)
%@end
%@if GM&&isa(varargin{1}, 'numeric')&& ToMean,
%@    % sprintf('standardize to mean! =dt')
%@    meanv=dt;
%@    shape=1/varargin{1};
%@    scale=meanv/shape;
%@    
%@    %sprintf('gamma %g-%g', shape,scale)
%@    
%@end
%@
%@%% Switch to computation method
%@
%@switch lower(method)
%@    
%@    case {'linint','regular','linear'}
%@        temp=tmin:dt:tmax;
%@        NoNorm=true;
%@        
%@    case {'gamma'}
%@        L=2*ceil((tmax-tmin)/dt);
%@        R=gamrnd(shape,scale,L,1);
%@        R=R(R>eps);
%@        temp=cumsum(R);
%@        
%@    case {'quasireg'}
%@        %    r = 1 + 2.*randn(100,1);
%@        maxdev=0.0001;
%@        L=2*ceil((tmax-tmin)/dt);
%@        R=dt+maxdev*randn(L,1);
%@        temp=cumsum(R);
%@        NoNorm=false;
%@        
%@    case {'poisson'}
%@        lambda=dt;
%@        L=10*ceil((tmax-tmin)/dt);
%@        R=poissrnd(lambda,L,1);
%@        R=R(R>eps);
%@        temp=cumsum(R);
%@        ratio=median(diff(temp))/dt;
%@        tx=temp((tmin*ratio<temp)&(temp<tmax*ratio));
%@        
%@    case {'missing'}
%@        %construct data on regular grid, take out varargin{1}*100 percent of the points randomly
%@        tx=tmin:dt:tmax;
%@        perc=varargin{1};
%@        L=length(tx);
%@        Rem=randi(L,floor(L*perc),1);
%@        tx(Rem)=[];
%@        temp=tx;
%@        
%@    case{'trapez'}
%@        %tx=tmin:dt:tmax;
%@        dt1= 8/5*dt;
%@        dt2=0.5*dt1;
%@        part1=tmin:dt2:(tmax/2-tmax/8);
%@        part2=max(part1):dt1:(tmax/2+tmax/8);part2(1)=[];
%@        part3=max(part2):dt2:tmax;part3(1)=[];
%@        tx=[part1 part2 part3];
%@        %tx=[tmin:dt2:(tmax/2-tmax/8)  (tmax/2-tmax/8 +dt2):dt1:(tmax/2+tmax/8) [(tmax/2+tmax/8 +dt1):dt2:tmax]]';
%@        temp=tx;
%@        
%@    case{'drift'}
%@        %tx=tmin:dt:tmax;
%@        dt1= 0.5*dt;
%@        dt2=2*dt;
%@        dtx=linspace(dt1,dt2,ceil((tmax-tmin)/dt));
%@        tx=cumsum(dtx)';
%@        temp=tx;
%@        
%@    case{'gdrift'}
%@        %tx=tmin:dt:tmax;
%@        dt1= 0.5*dt;
%@        dt2=2*dt;
%@        dtx=linspace(dt1,dt2,ceil((tmax-tmin)/dt));
%@        ddtx=gamrnd(dt,1,1,length(dtx)).*dtx;
%@        tx=cumsum(ddtx)';
%@        temp=tx;
%@        tx=temp((tmin<temp)&(temp<tmax));
%@        NoNorm=true;
%@        temp=tx;
%@        
%@    case{'gskew'}
%@        %tx=tmin:dt:tmax;
%@        SK=linspace(0.01,3,ceil(tmax-tmin)/dt)';
%@        AL=4./(SK.^2);
%@        BET=dt./AL;
%@        ddtx=gamrnd(AL,BET,length(SK),1);
%@        tx=cumsum(ddtx)';
%@        temp=tx;
%@        tx=temp((tmin<temp)&(temp<tmax));
%@        NoNorm=true;
%@        temp=tx;
%@        
%@        
%@    case{'sindrift'}
%@        % sinusoidally varying sampling (e.g. in winter more sparse than in
%@        % summer) modulating gamma distributed sampling intervals
%@        if isa(varargin{1}, 'numeric')
%@            PERI=varargin{1};
%@        else
%@            PERI=1;
%@        end
%@        AL=4./(0.5^2);
%@        BET=dt/AL;
%@        R=gamrnd(AL,BET,5*ceil(tmax-tmin)/dt,1);
%@        R=R(R>eps);
%@        dtemp=cumsum(R);
%@        size(R)
%@        size(dtemp)
%@        temp=cumsum(R.*abs(sin(2*pi*dtemp/PERI)));
%@        ratio=mean(diff(temp))/dt;
%@        %tx=temp((tmin<temp)&(temp<tmax));
%@        tx=temp((ratio*tmin<temp)&(temp<ratio*tmax));
%@        ToMean=true;    
%@end
%@%% Make sure that the median sampling rate is set
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@
%@if ToMean,
%@    ratio=mean(diff(temp)/dt);
%@    tx=temp((tmin*ratio<temp)&(temp<tmax*ratio));
%@    tx=tx*dt./mean(diff(tx));
%@elseif NoNorm
%@    tx=temp((tmin<temp)&(temp<tmax));
%@    
%@else
%@    ratio=median(diff(temp))/dt;
%@    tx=temp((tmin*ratio<temp)&(temp<tmax*ratio));
%@    tx=tx*dt./median(diff(tx));
%@    
%@end
%@
%@S=size(tx);
%@if S(2)>S(1),tx=tx';end
%<-- ASCII ends here -->
%<-- ASCII begins here: __gmi.m__ -->
%@function varargout = gmi(tx,x,ty,y,lag,varargin)
%@%GMI Mutual Information estimates for non-equidistantly sampled time series.
%@%    C = GMI(TX, X, TY, Y,SIGMA) returns the cross-mutual information function estimates
%@%    between the time series X and Y, for which the sampling times are given in
%@%    the column vectors TX and TY. The lags for which the Correlation functions
%@%    are estimated are given in the vector LAG=-T:DT:T. 
%@%
%@%    C = GMI(TX, X, TY, Y, LAG, 'kernel',H) computes the correlation using a kernel
%@%    width H, given in units of DT. H=0.5 amounts to a kernel width of 0.5*DT.
%@%
%@%    C = GMI(TX, X, TY, Y, LAG) returns the cross-correlation function at lags
%@%    given by vector LAG.
%@%
%@%    [C,LAGS] = GMI(...) returns the gaussian cross-mi and the vector of
%@%    time lags (LAGS). Note that LAG refers to the absolute time
%@%    differences, not the vector index differences (e.g. give -5 0 5 not -1
%@%    0 1 if the average spacing of the time vector is 5).
%@%
%@%
%@%    EXAMPLE: 
%@%       % We take two time series from stochastic processes X and Y, {tx,x} and
%@%       % {ty,y}, which are defined as follows: 
%@%       % Sampling times 
%@%       tx = (1:1001)'; ty = (1:1001)'; 
%@%       % Original signals 
%@%       x = randn(1001,1);
%@%       x(2:end) = 0.5*x(1:end-1)+randn(1000,1); 
%@%       y = randn(1001,1);
%@%       y(2:end) = 0.7*x(1:end-1).^2+randn(1000,1);
%@%
%@%       % We have, however, in reality, observed only 50% of the points in y, 
%@%       % yielding a time series {ty_miss,y_miss}:
%@%       rem = randi(length(ty), floor(length(ty)*0.5),1); 
%@%       ty_miss = ty; ty_miss(rem) = [];
%@%       y_miss = y; y_miss(rem) = [];
%@%
%@%       % We can now calculate the ACF for the time series {ty_miss,y_miss}:
%@%       % First we create a lag vector,
%@%       lag = -50:50;
%@%       % and we define the kernel width we want to employ:
%@%       H = 0.25; 
%@%       % Now the ACF of {ty_miss;y_miss} is given by
%@%       Cy = gmi(ty_miss, y_miss, ty_miss, y_miss, lag,'kernel', H); 
%@%       plot(lag,Cy), xlabel('Lag'), ylabel('Auto-MI C_{y}')
%@%       % and the CCF, indicating in this example the lag and the 
%@%       % magnitude of the coupling, is given by 
%@%       Cxy = gmi(tx, x, ty, y, lag,'kernel', H);
%@%       plot(lag,Cxy), xlabel('Lag'), ylabel('Cross-MI C_{xy}')
%@%
%@%    REFERENCE:
%@%       K. Rehfeld, N. Marwan, S. Breitenbach, J. Kurths: Late
%@%       Holocene Asian Monsoon Dynamics from small but complex paleoclimate
%@%       networks, 41(1), 3–19 (10013).
%@%
%@%    See also SIMILARITY,NEMI,ESF.
%@%    
%@% (c) Potsdam Institute for Climate Impact Research 10011-10013
%@% Kira Rehfeld, Norbert Marwan
%@% ver: 0.9  last rev: 10013-08-21
%@
%@nbins=10;%default value, if it isn't supplied
%@SIGMA=0.5;
%@
%@if length(varargin)>1,
%@    for i=1:(length(varargin)-1)
%@        bool1=strcmpi('nbins', varargin{i});
%@        if bool1, nbins=varargin{i+1};end %number of bins for weight matrix
%@        bool2=strcmpi('kernel', varargin{i})||strcmpi('width', varargin{i});
%@        if bool2, SIGMA=varargin{i+1};end %weight matrix kernel width
%@    end
%@    
%@end
%@%%%%% Checks for time series length etc
%@if length(tx) ~= length(x); error('tx and X not same length');
%@end;
%@if length(ty) ~= length(y); error('ty and Y not same length');
%@end;
%@
%@
%@% make sure column vectors are given
%@S1=size(x);
%@S2=size(y);
%@S3=size(tx);
%@S4=size(ty);
%@if S1(2)>S1(1),x=x';end
%@if S2(2)>S2(1),y=y';end
%@if S3(2)>S3(1),tx=tx';end
%@if S4(2)>S4(1),ty=ty';end
%@
%@dt=max(median(diff(tx)),median(diff(ty)));
%@tx=tx./dt;
%@ty=ty./dt;
%@l=nan(size(lag));
%@c=l;
%@
%@% repeat computation independently for each lag
%@for k=1:length(lag) 
%@    % make gaussian weighted reconstruction where time points are available
%@    [X,Y] = nest_scatter(tx,x,ty,y,SIGMA,-lag(k),0);
%@    % compute histogram-based MI for these reconstructions
%@    c(k)=MIhist(X,Y,nbins);
%@    l(k)=lag(k);
%@    clear X Y
%@end
%@
%@%% output of results
%@if nargout == 1
%@    varargout = {c};
%@elseif nargout == 2
%@    varargout(1) = {c};
%@    varargout(2) = {l};
%@elseif nargout == 0
%@    c
%@end
%@
%@
%@
%@end
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __imi.m__ -->
%@function varargout= imi(tx,x,ty,y,lag,varargin)
%@%IMI Interpolated Cross-MI estimates for non-equidistantly sampled time series.
%@% See GMI for syntax - this is just a helper function for comparison in SIMILARITY.
%@% ver: 0.9   last rev: 2013-08-21
%@%
%@% see also: NEXCF, SIMILARITY
%@
%@optargin = size(varargin,2);
%@stdargin = nargin - optargin;
%@nbins=10;%default value, if it isn't supplied in the argument
%@
%@if length(varargin)>1,
%@    for i=1:optargin,
%@        bool1=strcmpi('nbins', varargin{i});
%@        if bool1, nbins=varargin{i+1};end %number of bins for weight matrix
%@    end
%@    
%@end
%@
%@%%%%% Checks
%@if length(tx) ~= length(x); error('tx and X not same length');
%@end;
%@if length(ty) ~= length(y); error('ty and Y not same length');
%@end;
%@
%@
%@
%@% dt=mean(diff(lag));
%@
%@
%@S1=size(x);
%@S2=size(y);
%@S3=size(tx);
%@S4=size(ty);
%@
%@if S1(2)>S1(1),x=x';end
%@if S2(2)>S2(1),y=y';end
%@if S3(2)>S3(1),tx=tx';end
%@if S4(2)>S4(1),ty=ty';end
%@
%@%dt=mean(diff(lag));
%@dt=max(median(diff(tx)),median(diff(ty)));
%@tx=tx./dt;
%@ty=ty./dt;
%@l=nan(size(lag));
%@c=l;
%@
%@
%@
%@
%@for k=1:length(lag)
%@    
%@    txs=tx-lag(k);
%@
%@    if isempty(txs), error('can not calculate, time series empty'),end
%@    t0=max(min(txs),min(ty));
%@    t1=min(max(txs),max(ty));
%@    N=floor((t1-t0)); % already standardized to unit step 1
%@    tint=linspace(t0,t1,N);
%@    xint=interp1(txs,x,tint,'linear', 'extrap');
%@    yint=interp1(ty,y,tint,'linear', 'extrap');
%@    yint=yint-mean(yint);
%@    xint=xint-mean(xint);
%@    
%@    X=xint';
%@    Y=yint';
%@    c(k)=MIhist(X,Y,nbins);
%@    
%@    l(k)=lag(k);
%@end
%@
%@S=size(l);
%@if S(2)>S(1), l=l';end
%@
%@S=size(c);
%@if S(2)>S(1), c=c';end
%@%% output of results
%@if nargout == 1
%@    varargout = {c};
%@elseif nargout == 2
%@    varargout(1) = {c};
%@    varargout(2) = {l};
%@    elseif nargout == 4
%@    varargout(1) = {c};
%@    varargout(2) = {l};
%@        varargout(3) = {X};
%@        varargout(4) = {Y};
%@
%@elseif nargout == 0
%@    c
%@end
%@
%@
%@
%@    
%<-- ASCII ends here -->
%<-- ASCII begins here: __info.xml__ -->
%@<productinfo xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://www.mathworks.com/namespace/info/v1/info.xsd">
%@<?xml-stylesheet type="text/xsl" href="http://www.mathworks.com/namespace/info/v1/info.xsd"?>
%@
%@<matlabrelease>NaN</matlabrelease>
%@<name>NESTool</name>
%@<type>matlab</type>
%@<icon>$toolbox/matlab/general/matlabicon.gif</icon>
%@
%@<list>
%@
%@<listitem>
%@<label>Help</label>
%@<callback>helpwin NESTool/</callback>
%@<icon>$toolbox/matlab/general/bookicon.gif</icon>
%@</listitem>
%@
%@</list>
%@
%@</productinfo>
%<-- ASCII ends here -->
%<-- ASCII begins here: __intersectt.m__ -->
%@function [i1 i2] = intersectt(t1,t2)
%@% INTERSECTT Intersection of two time strings.
%@%   [I1 I2] = INTERSECTT(T1,T2) for vectors T1 and T2, returns
%@%   the index vectors I1 and I2 for which the time periods T1
%@%   and T2 are overlapping.
%@%
%@%   Example:
%@%      t1 = 1:50:1000;
%@%      t2 = 532:63:3000;
%@%
%@%      [i1 i2] = intersectt(t1,t2);
%@%      % returns the overlapping time vectors 
%@%      t1(i1) = [551 601 651 701 751 801 851 901 951]
%@%      t2(i2) = [532 595 658 721 784 847 910]
%@
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Norbert Marwan
%@% ver: 1.0   last rev: 2013-09-05
%@
%@
%@% allow only vectors
%@t1 = t1(:);
%@t2 = t2(:);
%@
%@% ensure that both vectors are increasing
%@% (this is arbitrary decision, could also be decreasing) 
%@if t1(1) > t1(end) | t2(1) > t2(end)
%@    error('Time vector should be increasing')
%@end
%@
%@% find start and end in time series 1
%@start1 = find(t2(1) <= t1,1);
%@end1 = find(t2(end) > t1,1,'last');
%@i1 = start1:end1;
%@
%@% find start and end in time series 2
%@start2 = find(t1(1) <= t2,1);
%@end2 = find(t1(end) > t2,1,'last');
%@i2 = start2:end2;
%@
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __nest_scatter.m__ -->
%@function varargout = nest_scatter(tx,x,ty,y,SIGMAG,lag,plt)
%@% NEST_SCATTER produces a scatterplot of two irregular time series.
%@% The function is based on Gaussian weighted local reconstructions for time
%@% series observations, as used in Rehfeld et al., 2013 for Mutual
%@% Information estimation (see reference).
%@%
%@%
%@%       EXAMPLE
%@%       % We take two time series from stochastic processes X and Y, {tx,x} and
%@%       % {ty,y}, which are defined as follows:
%@%       % Sampling times
%@%       tx = (1:1001)'; ty = (1:1001)';
%@%       % Original signals
%@%       x = randn(1001,1);
%@%       x(2:end) = 0.5*x(1:end-1)+randn(1000,1);
%@%       y = randn(1001,1);
%@%       y(2:end) = 0.7*x(1:end-1)+randn(1000,1);
%@%
%@%       % We have, however, in reality, observed only 50% of the points in y,
%@%       % yielding a time series {ty_miss,y_miss}:
%@%       rem = randi(length(ty), floor(length(ty)*0.5),1);
%@%       ty_miss = ty; ty_miss(rem) = [];
%@%       y_miss = y; y_miss(rem) = [];
%@%       % a scatterplot of (x,y_miss) is not possible, as they have different
%@%   	% lengths and observation times. Therefore local reconstructions of
%@%       % y_miss have to be used, where SIGMA determines the width of the
%@%       % reconstruction kernel and lag is a potential time offset:
%@%       SIGMA=1; % kernel width in units of dt, where
%@%       % dt=max(median(diff(tx)),median(diff(ty))); if the average spacing
%@%       % of the time series is 10 years, a kernel width of 0.5 results in an
%@%       % effective "smoothing" with a gaussian kernel smoother of bandwidth
%@%       % 5 years.
%@%       lag=0;
%@%       plt=0;% do not plot the scatterplot directly   
%@%       % obtain the time series 
%@%       [X,Y,T]=nest_scatter(tx,x,ty_miss,y_miss,SIGMA,lag,0)
%@%
%@%       % make the scatterplot, where the size of the rings indicates the
%@%       % time difference between the observations/ the amount of weight.
%@%       [X,Y,T,h,Wx,Wy]=nest_scatter(tx,x,ty_miss,y_miss,SIGMA,lag,1)
%@%       % Here, h is the handle to the scatterplot figure, Wx and Wy are
%@%       the weights for the reconstructed time series {T,X} and {T,Y}
%@%
%@%       % Please note that this missing data example does not necessarily 
%@%       % present a suitable application but is rather for illustration. The
%@%       % approach may be more sensible for irregular time series, than for
%@%       % those which only miss a few observations.
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 1.0  last rev: 2013-09-04
%@
%@tx_orig=tx;
%@ty_orig=ty;
%@dt=max(median(diff(tx)),median(diff(ty)));
%@tx=tx./dt;
%@% standardize the time series observation times by the mean sampling rate
%@% for easier computation
%@ty=ty./dt;
%@h=[];
%@ind=[];
%@X=[];
%@Y=[];
%@Wx=[];
%@Wy=[];
%@txs=tx+lag; % allowing for a delay between the time series
%@
%@
%@width=3;
%@wx_max=1/sqrt(2*pi*SIGMAG^2)*exp(0); % maximum weight of kernel
%@ind=0;
%@for i=1:length(txs)
%@    % step through txs to find suitable timepoints for reconstruction in y
%@    %find the timepoints in y which are not more than width from the
%@    %respective timepoint in x, txs(i).
%@    indy=((ty>(txs(i)-width))&ty<(txs(i)+width));
%@    %find the candidates in x for comparison (some smoothing due to kernel)
%@    indx=((txs>(txs(i)-width))&txs<(txs(i)+width));
%@    %calculate weights
%@    if ~any(indy)||~any(indx),
%@         %disp('jump')
%@        continue,
%@        
%@    end
%@    ind=ind+1;
%@    wx=1/sqrt(2*pi*SIGMAG^2)*exp(-(txs(indx)-txs(i)).^2./(2*SIGMAG^2))./wx_max;
%@    wy=1/sqrt(2*pi*SIGMAG^2)*exp(-(ty(indy)-txs(i)).^2./(2*SIGMAG^2))./wx_max;
%@    Wx(ind)=sum(wx)/width;
%@    Wy(ind)=sum(wy)/width;
%@    %P(ind,1)=txs(i);
%@    X(ind,1)=sum(wx.*x(indx))./sum(wx);
%@    Y(ind,1)=sum(wy.*y(indy))./sum(wy);
%@    T(ind,1)=tx(i);
%@    clear wx wy
%@end
%@%%
%@
%@for j=1:length(ty)
%@    %find the candidates in y for comparison
%@    indx=((tx>(ty(j)-width))&(tx<(ty(j)+width)));
%@    %find the candidates in y for comparison
%@    indy=((ty>(ty(j)-width))&ty<(ty(j)+width));
%@    %calculate weights
%@    if ~any(indy)||~any(indx),
%@        continue,
%@    end
%@    ind=ind+1;
%@    wx=1/sqrt(2*pi*SIGMAG^2)*exp(-(txs(indx)-ty(j)).^2./(2*SIGMAG^2))./wx_max;
%@    wy=1/sqrt(2*pi*SIGMAG^2)*exp(-(ty(indy)-ty(j)).^2./(2*SIGMAG^2))./wx_max;
%@    Wx(ind)=sum(wx)/width;
%@    Wy(ind)=sum(wy)/width;
%@    
%@    X(ind,1)=sum(wx.*x(indx))./sum(wx);
%@    Y(ind,1)=sum(wy.*y(indy))./sum(wy);
%@    T(ind,1)=ty(j);
%@    clear wx wy
%@    
%@end
%@
%@
%@indnan=isnan(X)|isnan(Y);
%@X(indnan)=[];
%@Y(indnan)=[];
%@T(indnan)=[];
%@Wx(indnan)=[];
%@Wy(indnan)=[];
%@
%@[Tsort,indsrt]=sort(T); % new actual time series, sorted increasingly
%@X=X(indsrt); Y=Y(indsrt);T=T(indsrt); Wx=Wx(indsrt);Wy=Wy(indsrt);
%@clear Tsort
%@if any(diff(T)<eps)
%@     %disp('non-distinctiveness could cause errors, removing double observations')
%@     ind=find(diff(T(:,1))<eps);
%@     T(ind)=[];X(ind)=[];Y(ind)=[];
%@end
%@
%@
%@%% Plot the scatterplot!
%@%h=plot(X,Y,'*');
%@
%@if plt,
%@    mrksz=2;col=[0.5 0.5 0.5];
%@    lw=2;
%@    scaling=(Wx+Wy)/2;
%@    for k=1:length(X)
%@        %plot(X(k),Y(k),'o','MarkerSize',mrksz*scaling(k), 'MarkerFaceColor',col.*scaling(k), 'MarkerEdgeColor',col.*scaling(k))
%@        h(k)=plot(X(k),Y(k),'ko','MarkerSize',mrksz*scaling(k),'LineWidth',lw*scaling(k));
%@        hold on;
%@    end
%@    xlabel(['X[ t+(' sprintf('%g',lag) ')]'])
%@    ylabel('Y[t]')
%@end
%@
%@
%@%% assign output arguments
%@if nargout==2,
%@    varargout{1}=X;
%@    varargout{2}=Y;
%@elseif nargout==3
%@    varargout{1}=X;
%@    varargout{2}=Y;
%@    varargout{3}=T*dt;% re-scale time index to original units
%@elseif nargout==6
%@    varargout{1}=X;
%@    varargout{2}=Y;
%@    varargout{3}=T*dt;
%@    varargout{4}=h;
%@    varargout{5}=Wx;
%@    varargout{6}=Wy;
%@end
%@
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __nexcf.m__ -->
%@function varargout = nexcf(tx, x, ty, y,lag,h, verbose)
%@% NEXCF Cross-correlation estimates for non-equidistantly sampled time series.
%@%    C = NEXCF(TX, X, TY, Y) returns the cross-correlation function estimates
%@%    between the time series X and Y, for which the sampling times are given in
%@%    the column vectors TX and TY. The lags for which the Correlation functions
%@%    are estimated are given in the vector LAG=-T:DT:T. If this vector is not
%@%    provided, where T=RANGE([MAX(TX(1),TY(1)),MIN(TX(END),TY(END))])/2 and
%@%    DT=T/10. Unless specified, the default kernel width is 0.25*DT.
%@%
%@%    C = NEXCF(TX, X, TY, Y, LAG, H) computes the correlation using a kernel
%@%    width H, given in units of DT. H=0.25 amounts to a kernel width of 0.25*DT.
%@%
%@%    C = NEXCF(TX, X, TY, Y, LAG) returns the cross-correlation function at lags
%@%    given by vector LAG.
%@%
%@%    [C,LAGS] = NEXCF(...) returns the correlation function and the vector of
%@%    time lags (LAGS). Note that LAG refers to the absolute time
%@%    differences, not the vector index differences (e.g. give -5 0 5 not -1
%@%    0 1 if the average spacing of the time vector is 5).
%@%
%@%
%@%    EXAMPLE:
%@%       % We take two time series from stochastic processes X and Y, {tx,x} and
%@%       % {ty,y}, which are defined as follows:
%@%       % Sampling times
%@%       tx = (1:201)'; ty = (1:201)';
%@%       % Original signals
%@%       x = randn(201,1);
%@%       x(2:end) = 0.5*x(1:end-1)+randn(200,1);
%@%       y = randn(201,1);
%@%       y(2:end) = 0.7*x(1:end-1)+randn(200,1);
%@%
%@%       % We have, however, in reality, observed only 50% of the points in y,
%@%       % yielding a time series {ty_miss,y_miss}:
%@%       rem = randi(length(ty), floor(length(ty)*0.5),1);
%@%       ty_miss = ty; ty_miss(rem) = [];
%@%       y_miss = y; y_miss(rem) = [];
%@%
%@%       % We can now calculate the ACF for the time series {ty_miss,y_miss}:
%@%       % First we create a lag vector,
%@%       lag = -50:50;
%@%       % and we define the kernel width we want to employ:
%@%       H = 0.25;
%@%       % Now the ACF of {ty_miss;y_miss} is given by
%@%       Cy = nexcf(ty_miss, y_miss, ty_miss, y_miss, lag, H);
%@%       plot(lag,Cy), xlabel('Lag'), ylabel('Auto-correlation C_{y}')
%@%       % and the CCF, indicating in this example the lag and the
%@%       % magnitude of the coupling, is given by
%@%       Cxy = nexcf(tx, x, ty, y, lag, H);
%@%       plot(lag,Cxy), xlabel('Lag'), ylabel('Cross-correlation C_{xy}')
%@%
%@%    REFERENCE:
%@%       Ref.: K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@%       correlation analysis techniques for irregularly sampled time series,
%@%       Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@%
%@%    See also SIMILARITY,NEMI,ESF.
%@
%@% (c) Potsdam Institute for Climate Impact Research 2011-2013
%@% Kira Rehfeld, Norbert Marwan
%@% ver: 1.4   last rev: 2013-09-02
%@
%@
%@%% output check
%@if nargout > 2
%@    error('Too many output arguments')
%@end
%@
%@%% input check
%@error(nargchk(2, 7, nargin, 'struct'))
%@
%@verb = 0;
%@useCuda = 1; %default
%@pd=0;
%@
%@if nargin == 7
%@    verb = verbose;
%@end
%@
%@% check if pdist2 is available
%@if ~exist('pdist2','file')
%@    % octave doesn't have pdist2, use modified functions at the end
%@    pd=1;
%@end
%@
%@% check cuda capability
%@try
%@   if exist('gpuDeviceCount','file')
%@       gpu = gpuDeviceCount;
%@   else
%@       gpu = 0;
%@   end
%@catch
%@    gpu = 0;
%@end
%@
%@
%@
%@if useCuda && gpu && verb
%@    disp('Wow! CUDA support available! You can enjoy accelerated calculation!')
%@end
%@if ~useCuda && gpu && verb
%@    disp('CUDA support available but not used.')
%@end
%@
%@% length of the time series x and y
%@Nx = length(x); Ny = length(y);
%@
%@% check the orientation of the vectors (vector multiplication in pdist...)
%@tx=tx(:);ty=ty(:);x=x(:);y=y(:);
%@
%@
%@% check for differences in length of time vector and observation vector
%@if (Nx ~= length(tx)) || (Ny ~=length(ty)),
%@    error('Vectors for sampling times and observations have to be of same length.'),
%@end
%@
%@% check for overlap in time-axes
%@t_min = max(min(tx), min(ty));
%@t_max = min(max(tx), max(ty));
%@
%@if (t_max-t_min)<eps
%@    error('The time-series do not overlap.')
%@end
%@
%@% check if h is given
%@if nargin < 6
%@    h = 0.25;
%@end
%@
%@% check if lag is given
%@% (if not, we use some default lag)
%@if nargin < 5
%@    T = range([t_min,t_max])/2;
%@    DT = T/10;
%@    %DT=10*mean(diff(T));
%@    lag(:,1) = -T:DT:T;
%@end
%@
%@% check for ACF or CCF
%@if (length(tx)==length(ty)) && ~any(x-y)
%@    flag_acf = 1;
%@    showType = 'Calculate ACF';
%@else
%@    flag_acf = 0;
%@    showType = 'Calculate XCF';
%@end
%@if verb, disp(showType), end
%@
%@%% data preprocessing
%@% make the time index dimensionless to make l and tx/ty comparable
%@if numel(lag) == 1
%@    dtlag = 1;
%@else
%@    %dtlag = abs(mean(diff(lag)));
%@    dtlag=max(mean(diff(tx)),mean(diff(ty)));
%@end
%@tx = tx/dtlag;
%@ty = ty/dtlag;
%@normlag = lag/dtlag;
%@
%@%normalise the data
%@x = (x - mean(x)) / std(x);
%@y = (y - mean(y)) / std(y);
%@
%@
%@
%@
%@%% main part of the calculation
%@
%@% Gauss kernel
%@b = @(t,h2,h2pi)(exp(-(t.^2)./h2) / h2pi);
%@
%@if useCuda && gpu
%@    % initialise correlation vector
%@    C = gpuArray(zeros(length(lag),1)); % correlation vector
%@    
%@    % pairwise time-distances between all observation
%@    t_dist = gpuArray(pdist2(tx,ty, @(xi,yi)(xi-yi))); t_dist = t_dist(:);
%@    t_distx = gpuArray(pdist2(tx,tx, @(xi,yi)(xi-yi)));
%@    t_disty = gpuArray(pdist2(ty,ty, @(xi,yi)(xi-yi)));
%@    
%@    % prodcut matrix of data
%@    xy_dist = gpuArray(pdist2(x,y, @(xi,yi)(xi.*yi))); xy_dist = xy_dist(:);
%@    
%@else
%@    % initialise correlation vector
%@    C = zeros(length(lag),1); % correlation vector
%@    if ~pd
%@        % pairwise time-distances between all observation
%@        t_dist = pdist2(tx,ty, @(xi,yi)(xi-yi)); t_dist = t_dist(:);
%@        t_distx = pdist2(tx,tx, @(xi,yi)(xi-yi));
%@        t_disty = pdist2(ty,ty, @(xi,yi)(xi-yi));
%@        
%@        % prodcut matrix of data
%@        xy_dist = pdist2(x,y, @(xi,yi)(xi.*yi)); xy_dist = xy_dist(:);
%@    else
%@        % octave -> use my_pdist
%@        t_dist = pdist2Diff(tx,ty); t_dist = t_dist(:);
%@        t_distx = pdist2Diff(tx,tx);
%@        t_disty = pdist2Diff(ty,ty);
%@        
%@        % prodcut matrix of data
%@        xy_dist = pdist2Prod(x,y); xy_dist = xy_dist(:);
%@    end
%@end
%@
%@% lower and upper limit of time differences
%@min_dist = min(t_dist(:));
%@max_dist = max(t_dist(:));
%@
%@h2 = 2*h^2;
%@h2pi = sqrt(h2 * pi);
%@
%@if verb, fprintf('lags... %03d %%',0), end
%@for i = 1:length(lag)
%@    if verb
%@        percent = round((i/length(lag))*100);
%@        fprintf('\b\b\b\b\b\b %03d %%',percent);
%@    end
%@    
%@    % check whether lag is outside possible time difference range
%@    if (normlag(i) <= min_dist) || (normlag(i) >= max_dist)
%@        C(i) = NaN;
%@        continue
%@    end
%@    
%@    % new distance matrix corresponding to considered lag
%@    t_dist_d = t_dist+normlag(i);
%@    
%@    % count the number of observations falling inside the kernel
%@    kernel_content = (t_dist_d <= 5*h) & (t_dist_d >= -5*h);
%@    % check whether this number is smaller than 5
%@    % (if yes, we do not believe in the correlation value)
%@    if sum(kernel_content) < 5
%@        C(i) = NaN;
%@        warning([sprintf('Lag %0.3f outside confident time region. Consider smaller range for lags.',lag(i))]);
%@        continue
%@    end
%@    
%@    % matrix of the weights (= apply kernel)
%@    WL = b(t_dist_d,h2,h2pi);
%@    Weight = sum(WL);
%@    % sum of the product of kernel with products x*y is the correlation
%@    C_ = sum(xy_dist.*WL) ./ Weight;
%@    C(i) = C_;
%@end
%@
%@
%@%% normalise the correlation values to a range [-1:1] using the variances of the time series
%@%matrix of the weights (= Gaussian kernel) for lag = 0
%@%WL = (exp(-((t_dist(:)).^2)./(2*h^2))) / sqrt(2*pi*h^2);
%@WLx = (exp(-((t_distx(:)).^2)./h2)) / h2pi;
%@WLy = (exp(-((t_disty(:)).^2)./h2)) / h2pi;
%@
%@WeightX = sum(WLx);
%@WeightY = sum(WLy);
%@
%@% product matrix of data x and y
%@if ~pd
%@    x_dist = pdist2(x,x, @(xi,yi)(xi.*yi));
%@    y_dist = pdist2(y,y, @(xi,yi)(xi.*yi));
%@else
%@    x_dist = pdist2Prod(x,x);
%@    y_dist = pdist2Prod(y,y);
%@end
%@
%@% corrected variance as weighted mean, taking into account the irregularity
%@varX = sum( x_dist(:).*WLx) ./ WeightX;
%@varY = sum( y_dist(:).*WLy) ./ WeightY;
%@
%@%% normalise the correlation estimates
%@C = C/(sqrt(varX) * sqrt(varY));
%@
%@if any(C(:)>1),C=C/max(C(:));end
%@% numerically, it can happen that the variance is underestimated-> the
%@% correlation, as the ratio of covariance to variance overestimated. As a
%@% rule, we forbid the Correlation to be larger than 1.
%@%%
%@% transform gpuArray data to a standard Matlab array
%@if useCuda && gpu
%@    C = gather(C);
%@end
%@
%@%% output of results
%@if nargout == 1
%@    varargout = {C};
%@elseif nargout == 2
%@    varargout(1) = {C};
%@    varargout(2) = {lag};
%@elseif nargout == 0
%@    C
%@end
%@
%@%%%%%%% Define some functions that are not available in Octave (pdist2-
%@%%%%%%% replacements)
%@    function D = pdist2Prod(X,Y)
%@        n = length(X);
%@        m = length(Y);
%@        D = X * Y';
%@    end
%@    function D = pdist2Diff(X,Y)
%@        n = length(X);
%@        m = length(Y);
%@        D = ones(n,1) * Y' - X * ones(1,m) ;
%@    end
%@end
%<-- ASCII ends here -->
%<-- ASCII begins here: __pspek.m__ -->
%@function [freq, power,h]=pspek(acf, dt, varargin)
%@%PSPEK computes the power spectrum for a given autocorrelation function.
%@%   [F, P]=pspek(ACF, DT) computes the power spectrum P over the
%@%   frequencies F via the fast fourier transform of the autocorrelation
%@%   function ACF. DT is the spacing of the ACF lags. If a third argument,
%@%   PT, is given, this true/false value indicates whether a plot of the
%@%   spectrum is desired. Optional arguments need to be given in pairs, i.e.
%@%   {'ccf', 1} or {'nfft',1024}: [F, P]=pspek(ACF, DT,'nfft',1024).
%@%   
%@%   EXAMPLE:
%@%  
%@%     tx = (1:1001)';
%@%     % Original signals
%@%     x = randn(1001,1);
%@%     % signal is mix of sinusoid and noise, partly autocorrelated  
%@%     x(2:end) = sin(2*pi*tx(2:end)/50)+0.2*x(1:end-1)+0.1*randn(1000,1);
%@%     x=(x-mean(x))./std(x);
%@%     lag=-500:1:500; % lags for acf
%@%     [c]=similarity(tx,x,tx,x,lag,'gXCF'); % compute acf
%@%     plt=true;
%@%     [f,p]=pspek(c,mean(diff(lag)),plt);
%@%
%@% see also: SIMILARITY,NEXCF
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 1.0  last rev: 2013-10-14
%@
%@
%@
%@h=[];
%@cross=false;
%@NFFT = 2^nextpow2(length(acf));
%@pt=0; % do not plot by default
%@
%@% check optional input arguments
%@if length(varargin)>1,
%@    for i=1:length(varargin),
%@        bool1=strcmpi('ccf', varargin{i});
%@        if bool1, cross=varargin{i+1};end %cross spectrum 
%@         bool2=strcmpi('nfft', varargin{i});
%@        if bool2, NFFT=varargin{i+1};end %number of frequencies for fft
%@        
%@    end
%@elseif length(varargin)==1
%@    % plotting argument
%@    pt=varargin{1};
%@end
%@
%@if pt,
%@disp(sprintf('Computing FFT of ACF with %g frequencies in FFT',NFFT))
%@end
%@
%@
%@Fs = 1/dt;  % Sampling frequency
%@T = 1/Fs; % Sampling period
%@
%@Y = fft(acf,NFFT)/length(acf); % PSD
%@f = [Fs/2*linspace(0,1,NFFT/2+1)]'; % frequency vector
%@
%@
%@%reduction to positive frequencies, i.e. doubling of the contained power.
%@if ~cross
%@    power=abs(Y);%should be pos. semidefinite, but usually isn't
%@    power=(power(1:floor(NFFT/2)+1))';
%@    power=power./abs(sum(power));%Y(1)
%@    freq=f(1:end);
%@else
%@    Y(1)=[] ;
%@    %power=real(Y); 
%@    power=abs(Y);%should be pos. semidefinite, but usually isn't
%@    power=power./sum(power);
%@    
%@    freq=[-flipud(f(2:end));f(2:end-1)];
%@    freq=flipud(fftshift(freq));
%@    f= [Fs*linspace(0,1,NFFT)]';
%@end
%@
%@
%@if pt,
%@    h=plot(1./(freq),power);
%@    xlabel('period [unit]')
%@    ylabel('power')
%@end
%@
%@S=size(freq);
%@if S(2)>S(1),
%@    freq=freq';
%@end
%@S=size(power);
%@if S(2)>S(1),
%@    power=power';
%@end
%@
%@end
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __simgram.m__ -->
%@function [c_out, l_out, t_out, N_out] =simgram(tx,x,ty,y,varargin)
%@%SIMGRAM calculates windowed cross similarity between two time series.
%@%
%@%   C = SIMGRAM(tx,x,ty,y,'overlap',OVERLAP,'windw',windw,'lag',LAGS)
%@%   calculates the windowed cross similarity between the signals in vector
%@%   x and vector y with their corresponding time axis tx and ty,
%@%   respectively. The similarity is estimated for each lag shift in the
%@%   vector LAGS. SIMGRAM splits the signals into overlapping segments and
%@%   forms the columns of C with their cross similarity values for each lag
%@%   as specified by vector LAGS. Each column of C contains the cross
%@%   similarity function between the short-term, time-localized signals
%@%   {tx,x} and {ty,y}. Time increases linearly across the columns of C,
%@%   from left to right.  Lag increases linearly down the rows. If T is the
%@%   length of the signals (in time units), C is a matrix with Nlags rows
%@%   (Nlags is the length of the lag vector) and
%@%         k = fix((T-windw)/((1-overlap)*windw)+1)
%@%   columns.
%@%
%@%   Note: LAGS is a vector, typically symmetric, and given in the original
%@%   units of the time series tx and ty temporal spacing and WINDW has to
%@%   be given in the same time unit!
%@%   Different methods  for cross-similarity estimation are available and
%@%   chosen according to the argument 'method':
%@%       'gXCF': Gaussian-kernel cross-correlation (Rehfeld, NPG 2011) or NEXCF
%@%       'iXCF': Interpolation to a regular timescale for both time series, 
%@%               standard xcorr estimator (Rehfeld, 2011)
%@%        'gMI': Gaussian-kernel mutual information (Rehfeld, ClimDyn 2012)
%@%               or GMI
%@%        'iMI': Interpolation to regular timeseries, then use MIhist
%@%        'ESF': Event Synchronization function (Rehfeld, CPD 2013)
%@%
%@%
%@%   [C,L,T] = SIMGRAM(...) returns a column of lag L and one of time T
%@%   at which the similarity coefficients are computed. L has length equal
%@%   to the number of rows of C, T has length k.
%@%
%@%
%@%
%@%   C = SIMGRAM(tx,x,ty,y) calculates windowed cross similarity using
%@%   default settings; the defaults are LAG=-(10:1:10)*dt,where dt is the
%@%   common mean sampling rate, WINDW = floor(0.1*T), OVERLAP = 0 and
%@%   METHOD is 'gXCF', i.e. Gaussian-kernel-based cross correlation.
%@%
%@%   C = SIMGRAM(tx,x,ty,y,optargs) calculates windowed cross similarity
%@%   using optional arguments to control the computation. The arguments are
%@%   expected to come in pairs, see the example below.
%@%               'method' - and one of {'gXCF','iXCF','gMI', 'iMI','ESF'}
%@%               'overlap'- and a value between [0 and 1] 
%@%               'lag' - symmetric lag vector, e.g. [-1:1:1]
%@%               'windw' - window width in time units, should contain >80 points
%@%
%@%   SIMGRAM(tx,x,ty,y) with no output arguments plots the windowed
%@%   Gaussian-kernel cross correlation using the current figure.
%@%
%@%
%@%
%@%   Example
%@%       x = randn(1000,1);tx=1:length(x);ty=tx;
%@%       y = x + linspace(0,10,length(tx))'.*randn(length(x),1);
%@%       y=(y-mean(y))./std(y);
%@%       simgram(tx,x,ty,y,'method','gXCF','lag',[-10:1:10])
%@%       % or with MI:
%@%       simgram(tx,x,ty,y,'method','gMI','lag',[-10:1:10],'overlap',0.5)
%@%
%@%
%@%    REFERENCE:
%@%
%@%       gXCF: K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@%       correlation analysis techniques for irregularly sampled time series,
%@%       Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@%       gMI: K. Rehfeld, N. Marwan, S. Breitenbach, J. Kurths: Late
%@%       Holocene Asian Monsoon Dynamics from small but complex paleoclimate
%@%       networks, 41(1), 3???19 (2013).
%@%       ESF: K. Rehfeld and J. Kurths: Similarity measures for irregular
%@%       and age uncertain time series, submitted to Clim. Past. Discuss,
%@%       2013.
%@%
%@%
%@%   See also NEXCF,SIMILARITY,GMI
%@%
%@% Copyright (c) 2013
%@% Kira Rehfeld, Norbert Marwan
%@% Potsdam Institute for Climate Impact Research, Germany
%@% ver: 0.9c  last rev: 2013-12-19
%@% log: added detrending option
%@
%@
%@
%@error(nargchk(4,14,nargin))
%@verbose = 0;
%@kernelwidth = .25;
%@
%@
%@tx = tx(:); ty = ty(:);
%@x = x(:); y = y(:);
%@
%@
%@% check input and inital setting of parameters
%@
%@nx = length(x); ny = length(y);
%@% if length(tx) ~= nx | length(ty) ~= ny
%@%     error('Time vector and data vector have to be of equal length.')
%@% end
%@% if nx < ny    % zero-pad x if it has length less than y
%@%     x(ny) = 0; nx = ny;
%@% end
%@%
%@% if ny < nx    % zero-pad y if it has length less than x
%@%     y(nx) = 0;
%@% end
%@
%@
%@dt=max(median(diff(tx)),median(diff(ty))); % check approproate common sampling rate
%@lag = (-10:1:10)*dt;
%@overlap = 0;
%@method='gXCF';
%@disp=1; 
%@detr=0;
%@
%@% time start and end
%@%t1 = %min([tx;ty]);
%@t1=max([min(tx),min(ty)]);
%@t2=min([max(tx),max(ty)]);
%@
%@%t2 = max([tx;ty]);
%@windw = floor((t2-t1)/10);% default: 10 time windows spanning the overlapping period
%@if range([t1,t2]) < windw
%@    error('Window size has to be larger than the overlapping time span.')
%@end
%@
%@% cut time series to overlapping period to avoid empty windows[B
%@[tx,x]=extractwindow(tx,x,t1,t2,1);
%@[ty,y]=extractwindow(ty,y,t1,t2,1);
%@
%@if length(varargin)>1,
%@    for i=1:length(varargin),
%@        bool1=strcmpi('lag', varargin{i});
%@        
%@        if bool1, lag=varargin{i+1};end % lag vector
%@        bool4=strcmpi('windw', varargin{i})||strcmpi('width', varargin{i});
%@        if bool4, windw=varargin{i+1};
%@            if windw < 0, error('Requires positive time period for window length.'), end
%@            if length(windw) > 1, error('Requires WINDW to be a scalar.'), end
%@        end
%@        
%@        bool2=strcmpi('overlap',varargin{i});
%@        if bool2,overlap=varargin{i+1};
%@            if overlap < 0, error('Requires positive integer value for OVERLAP.'), end
%@            if length(overlap) > 1, error('Requires OVERLAP to be a scalar.'), end
%@            if overlap >= 1, error('Requires OVERLAP to be strictly less than 1. OVERLAP=0 means no overlap, OVERLAP=0.5 is equal to an overlap of 50% of the window width.'), end
%@        end
%@        
%@        bool3=strcmpi('method',varargin{i});
%@        if bool3,method=varargin{i+1};%gXCF,iXCF,gMI,iMI,ESF are the options
%@        end
%@        
%@        bool4=strcmpi('verbose',varargin{i});
%@        if bool4,verbose=varargin{i+1};% 0->no waitbar, 1-> waitbar
%@        end  
%@        
%@        bool5=strcmpi('detrending',varargin{i});
%@        if bool5,detr=varargin{i+1};% detrending in the individual windows
%@            %'detrending'
%@        end    
%@    end
%@    
%@end
%@
%@
%@
%@
%@% vector of time windows
%@[center, links, rechts] = slidingwindow(t1,t2,windw,overlap);
%@nw=length(links);% no of windows
%@C = zeros(length(lag), nw);
%@N = zeros(nw, 2);
%@
%@
%@if verbose, h = waitbar(0,'Compute cross similarity for each window'); end
%@%for t = tw, if verbose, waitbar(cnt/nw), end
%@for t=1:length(links)
%@    if verbose, waitbar(t/length(links)), end
%@    
%@    [txw,xw]=extractwindow(tx,x,links(t),rechts(t),1);
%@    [tyw,yw]=extractwindow(ty,y,links(t),rechts(t),1);
%@    
%@    
%@    if isempty(txw) || isempty(tyw)
%@        % if this time segment is empty-> warning!
%@        warning(sprintf('Non-overlapping segment found at %f-%f',t,t+windw))
%@        Cxy=NaN(length(lag),1);
%@    else
%@        % detrending (optional)
%@        if detr,
%@            xw=xw-gausssmooth(txw,xw,detr);
%@            yw=yw-gausssmooth(tyw,yw,detr);
%@%              figure()
%@%              plot(txw,xw,'r')
%@%              hold on;
%@%              plot(tyw,yw,'b')
%@        end
%@        %Cxy = nexcf(txw,xw,tyw,yw,lag,kernelwidth,0);
%@        Cxy = similarity(txw,xw,tyw,yw,lag,method);
%@        
%@    end
%@    C(:,t) = Cxy;
%@    N(t,1) = length(txw);
%@    N(t,2) = length(tyw);
%@    clear twx tyw xw yw Cxy
%@end
%@
%@if verbose, delete(h), end
%@
%@warning on
%@
%@if any(N(:)<50),
%@warning('less than 50 data points in a time window! re-set limits')
%@end
%@
%@% display and output result
%@if nargout == 0
%@    newplot
%@    % imagesc(tw+windw/2, lag, C)
%@    imagesc(center, lag, C)
%@    xlabel('Time'), ylabel('Lag'), axis xy
%@    title(['Windowed ' method], 'fontweight', 'bold')
%@    hcb=colorbar;
%@    get(hcb)
%@elseif nargout == 1,
%@    c_out = C;
%@elseif nargout == 2,
%@    c_out = C;
%@    l_out = lag;
%@elseif nargout == 3,
%@    c_out = C;
%@    l_out = lag;
%@    %t_out = tw;
%@    t_out = center;
%@elseif nargout == 4,
%@    c_out = C;
%@    l_out = lag;
%@    %t_out = tw;
%@    t_out = center;
%@    N_out = N;
%@end
%@
%@
%@function Y = normalize(X)
%@Y = X - repmat(mean(X), size(X,1), 1);
%@s = sqrt( sum(conj(Y) .* Y) / (size(Y,1) - 1) );
%@Y = Y ./ repmat(s, size(Y,1), 1);
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __similarity.m__ -->
%@function [C, varargout] = similarity(tx,x,ty,y,lag,method,varargin)
%@%SIMILARITY computes statistical associations for irregular time series.
%@%   C=SIMILARITY(TX,X,TY,Y,'gXCF') returns Gaussian kernel correlation
%@%   estimates for the statistical association between the time series
%@%   {TX,X} and {TY,Y}, for which the sampling times are given in the column
%@%   vectors TX and TY and the observed values in X and Y. The lags for
%@%   which the associations are estimated are given in the vector
%@%   LAG=-T:DT:T, where DT is the lag spacing. SIMILARITY acts as a wrapper
%@%   for different methods to estimate the interdependences. Different
%@%   methods are available and chosen according to the argument:
%@%       'gXCF': Gaussian-kernel cross-correlation (Rehfeld, NPG 10011) or NEXCF
%@%       'iXCF': Interpolation to a regular timescale for both time series, standard xcorr estimator (Rehfeld, 10011)
%@%        'gMI': Gaussian-kernel mutual information (Rehfeld, ClimDyn 10012) or NEMI
%@%        'iMI': Interpolation to regular timeseries, then use NEMI
%@%        'ESF': Event Synchronization function (Rehfeld, CPD 2013)
%@%        'LS': all of the above (to investigate the Linkstrength (Rehfeld, CPD 2013)
%@%
%@%   Without additional input standard parameters are chosen for each
%@%   method. Additional input can be, for example,
%@%   C=SIMILARITY(TX,X,TY,Y,'gXCF','kernel',0.5) to set a specific kernel
%@%   width,C=SIMILARITY(TX,X,TY,Y,'gMI','nbins',5) to set 5 instead of 10
%@%   (default) bins for mutual information estimation,(...,'biascorr',1)
%@%   enforces bias-correction with uncorrelated surrogates for mutual
%@%   information,(...,'conflevel',0.8) sets an 80% confidence threshold for
%@%   event synchronization.
%@%
%@%
%@%       % EXAMPLE
%@%       % We take two time series from stochastic processes X and Y, {tx,x} and
%@%       % {ty,y}, which are defined as follows:
%@%       % Sampling times
%@%       tx = (1:1001)'; ty = (1:1001)';
%@%       % Original signals
%@%       x = randn(1001,1);
%@%       x(2:end) = 0.5*x(1:end-1)+randn(1000,1);
%@%       y = randn(1001,1);y_nl=y;
%@%       y(2:end) = 0.7*x(1:end-1)+randn(1000,1);
%@%       y_nl(2:end) = 0.7*x(1:end-1).^2+randn(1000,1);
%@%       % We have, however, in reality, observed only 80% of the points in y,
%@%       % yielding a time series {ty_miss,y_miss}:
%@%       rem = randi(length(ty), floor(length(ty)*0.2),1);
%@%       ty_miss = ty; ty_miss(rem) = [];
%@%       y_miss = y; y_miss(rem) = [];
%@%       y_miss_nl=y_nl;y_miss_nl(rem)=[];
%@%
%@%       % We can now calculate the ACF for the time series {ty_miss,y_miss}:
%@%       % First we create a lag vector
%@%       lag = -50:50;
%@%       % Now the ACF of {ty_miss;y_miss} is given by
%@%       Cy = similarity(ty_miss, y_miss, ty_miss, y_miss, lag,'gXCF');
%@%       plot(lag,Cy), xlabel('Lag'), ylabel('Auto-correlation C_{y}')
%@%       % and the CCF, indicating in this example the lag and the
%@%       % magnitude of the coupling, is given by
%@%       Cxy = similarity(tx, x, ty_miss, y_miss, lag, 'gXCF');
%@%       plot(lag,Cxy), xlabel('Lag'), ylabel('Cross-correlation C_{xy}')
%@%       % For a nonlinear relationship, for example a quadratic dependence,
%@%       % the cross correlation will not give a strong indication of the
%@%       % coupling:
%@%       C_nl_xcf = similarity(tx, x, ty_miss, y_miss_nl, lag, 'gXCF');
%@%       plot(lag,C_nl_xcf), xlabel('Lag'),
%@%       ylabel('cross-dependence')
%@%       % and compare this to the output of the MI or ESF, where a peak
%@%       % indicates the lag of coupling (in this example, -1).
%@%       C_nl_gmi = similarity(tx, x, ty_miss, y_miss_nl, lag, 'gMI');
%@%       C_nl_esf = similarity(tx, x, ty_miss, y_miss_nl, lag, 'ESF');
%@%       hold on; plot(lag,C_nl_gmi,'r');
%@%       hold on; plot(lag,C_nl_esf,'g');
%@%       legend('gXCF','gMI','ESF')
%@%
%@%       % Link strength - all similarity measures
%@%        C_ls = similarity(tx, x, ty_miss, y_miss_nl, lag, 'LS');
%@%       % or as a subset, i.e. 'gXCF', 'gMI' and 'ES'
%@%       C_lsss = similarity(tx, x, ty_miss, y_miss_nl, lag, 'LS','subset',[1,3,5]);
%@%
%@%
%@%    REFERENCE:
%@%       gXCF: K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@%       correlation analysis techniques for irregularly sampled time series,
%@%       Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@%       gMI: K. Rehfeld, N. Marwan, S. Breitenbach, J. Kurths: Late
%@%       Holocene Asian Monsoon Dynamics from small but complex paleoclimate
%@%       networks, 41(1), 3???19 (2013).
%@%       ESF: K. Rehfeld and J. Kurths: Similarity measures for irregular
%@%       and age uncertain time series, submitted to Clim. Past. Discuss,
%@%       2013.
%@%
%@%
%@% see also: NEXCF,GMI,EVENTSYNCHRO
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 0.9b  last rev: 2014-01-17
%@
%@
%@
%@c=[];
%@
%@x=(x-mean(x))/std(x); % centralize & standardize time series
%@y=(y-mean(y))/std(y);
%@dt=max(median(diff(tx)),median(diff(ty))); % check approproate common sampling rate
%@tx=tx./dt;
%@ty=ty./dt;
%@
%@%*****%****%*****%****for bias correction of MI
%@%biascorr=1;
%@Nsur=100;
%@varargout{1}=NaN;
%@varargout{2}=NaN;
%@%*****%****%*****%**** pre-defined options
%@mitoccscale=1;% convert MI to CC scale ([0,1])
%@biascorr=0; % 0 if no bias correction for MI desired
%@nbins=10; % no of bins for MI estimation
%@SIGMA_xcf=0.25; % kernel width for gXCF
%@SIGMA_gmi=0.25;
%@conflevel=0.9;% threshold level for ESF
%@subset=1;
%@similaritymeas={'gXCF','iXCF','gMI','iMI','ESF'};
%@bool6=0;
%@
%@%*****%****%*****%**** check optional input parameters
%@if length(varargin)>1,
%@    for i=1:length(varargin),
%@        bool1=strcmpi('nbins', varargin{i});
%@        if bool1, nbins=varargin{i+1};end %number of bins for weight matrix
%@        bool4=strcmpi('kernel', varargin{i})||strcmpi('width', varargin{i});
%@        if bool4, SIGMA=varargin{i+1};SIGMA_xcf=SIGMA;SIGMA_gmi=SIGMA;end %weight matrix kernel width
%@        bool2=strcmpi('biascorr',varargin{i});
%@        if bool2,biascorr=varargin{i+1};end % correct bias using autocorrelated ts
%@        bool3=strcmpi('MItoCCscale',varargin{i});
%@        if bool2,mitoccscale=varargin{i+1};end % correct bias using autocorrelated ts
%@        bool5=strcmpi('conflevel', varargin{i});
%@        if bool5,confl=varargin{i+1};end;
%@        bool6=strcmpi('subset', varargin{i})&&strcmpi(method,'LS');
%@        if bool6,subset=varargin{i+1};end;
%@    end
%@    
%@end
%@
%@% if method='LS'and no additional argument --> use all methods available
%@if strcmpi(method,'LS')&&~bool6
%@    subset=1:5;
%@    bool6=1;
%@end
%@
%@%% Compute similarity
%@
%@for k=1:length(subset) % loop through all if 'Linkstrength' is desired
%@    
%@    if bool6
%@        method=similaritymeas{subset(k)}; % use method no k if linkstrength is computed
%@    end
%@    
%@    switch method
%@        case {'gXCF','GXCF','gxcf'}
%@            %gaussian kernel correlation, Rehfeld et al., 10011
%@            if isempty(varargin),
%@                SIGMA_xcf=0.25;
%@            end
%@            
%@            c=nexcf(tx,x,ty,y,lag,SIGMA_xcf,0);
%@            
%@            if (any(isnan(c))) % try larger kernel width
%@                warning('Gaussian kernel width insufficient -- use larger kernel')
%@                SIGMA_xcf=1;
%@                c=nexcf(tx,x,ty,y,lag,SIGMA_xcf);
%@            end
%@            
%@        case {'iXCF','ixcf','IXCF'}
%@            % interpolate the time series to a regular grid and quantify similarity
%@            % using standard matlab xcorr (based on FFT, expects regular sampling)
%@            if (all(diff(tx)==dt)&&all(diff(ty)==dt)),%if regularly sampled
%@                % find which of the time series is maybe shorter
%@                S=[length(tx) length(ty)];
%@                if S(2)>S(1)
%@                    % ty as reference for tint
%@                    tint=ty;
%@                    xint=interp1(tx,x,tint,'linear','extrap');
%@                    yint=y;
%@                else
%@                    % tx as reference
%@                    tint=tx;
%@                    yint=interp1(ty,y,tint,'linear','extrap');
%@                    xint=x;
%@                end
%@                xint=xint-mean(x);
%@                yint=yint-mean(y);
%@                
%@            else %interpolate to common sampling
%@                %   dt=median(diff(lag)); % alternate option, if diff(lag) is not
%@                %   equal to max(median(diff(tx)),median(diff(ty)));
%@                Lx=max(tx);Lx1=min(tx);
%@                Ly=max(ty);Ly1=min(ty);
%@                %tint=linspace(min(Lx1,Ly1),max(Lx,Ly),floor((max(Lx,Ly)-min(Lx1,Ly1)))./dt);
%@                % mean time difference has been standardized above!
%@                tint=linspace(min(Lx1,Ly1),max(Lx,Ly),floor((max(Lx,Ly)-min(Lx1,Ly1))));
%@
%@                xint=interp1(tx,x,tint,'linear','extrap');
%@                yint=interp1(ty,y,tint,'linear','extrap');
%@                yint=yint-mean(yint); xint=xint-mean(xint);
%@                if any(isnan(xint))||any(isnan(yint)), error('interpolation error!'),end
%@            end
%@            % compute covariance
%@            [c,dummy]=xcorr(xint,yint,(length(lag)-1)/2, 'unbiased');
%@            c=c./(std(xint)*std(yint)); % compute correlation
%@            c=c';%column vector
%@            clear dummy % [c,~] does not work in Octave
%@        case {'gMI','gmi','GMI'}
%@            if isempty(varargin)
%@                SIGMA_gmi=1;
%@            end
%@            c0=gmi(tx,x,ty,y,lag,'nbins',nbins,'kernel',SIGMA_gmi);
%@            if biascorr,
%@                %disp('correcting MI Bias')
%@                Xsur=AR1Sur(tx,x,Nsur);
%@                Ysur=AR1Sur(ty,y,Nsur);
%@                cB=NaN(length(c0),Nsur);
%@                for i=1:Nsur
%@                    cB=gmi(tx,Xsur(:,i),ty,Ysur(:,i),lag,'kernel',SIGMA_gmi,'nbins',nbins);
%@                end
%@                Bias=mean(cB,2);
%@                BiasStd=std(cB,0,2);
%@                c=c0-Bias;
%@                varargout{1}=Bias;
%@                varargout{2}=BiasStd;
%@            else
%@                c=c0;
%@            end
%@            
%@            if mitoccscale
%@                c=real(MItoCCscale(c));
%@            end
%@            
%@        case {'iMI','imi','IMI'}
%@            % MI with prior interpolation to mean data sampling
%@            c0=imi(tx,x,ty,y,lag,'nbins',nbins);
%@            [c0 l0 X Y]=imi(tx,x,ty,y,lag,'nbins',nbins);
%@            
%@            if biascorr,
%@                %disp('correcting MI Bias')
%@                Xsur=ar1sur(tx,x,Nsur);
%@                Ysur=ar1sur(ty,y,Nsur);
%@                cB=NaN(length(c0),Nsur);
%@                for i=1:Nsur
%@                    cB(:,i)=imi(tx,Xsur(:,i),ty,Ysur(:,i),lag,'nbins',nbins);
%@                end
%@                Bias=mean(cB,2);
%@                BiasStd=std(cB,0,2);
%@                c=c0-Bias;
%@                varargout{1}=Bias;
%@                varargout{2}=BiasStd;
%@            else
%@                c=c0;
%@            end
%@            if mitoccscale
%@                c=real(MItoCCscale(c));
%@            end
%@            
%@        case {'evsync', 'ES','ESF','esf'}
%@            if isempty(varargin),
%@                conflevel=0.9;
%@            end
%@            for l=1:length(lag)
%@                [Q(l) q(l)] = eventsynchro(tx,x,ty,y,lag(l),conflevel);
%@            end
%@            % Q gives an association strength
%@            % q gives a direction, if >1, x "drives" y, if <1, other way around
%@            c=Q;
%@            varargout{1}=q;
%@        otherwise
%@            error('no appropriate method given to similarity - check spelling')
%@    end
%@    
%@    
%@    S=size(c);
%@    
%@    if S(2)>S(1),
%@        c=c';
%@    end
%@    
%@    C(:,k)=c;
%@    clear c
%@end
%@
%@%%%% additional function
%@    function [rmi] = MItoCCscale(mi)
%@        % MITOCCSCALE transforms Mutual Information estimates to a [0 1] scale.
%@        %
%@        %   MItoCCscale transforms the estimated Mutual Information mi to a scale
%@        %   comparable to the cross correlation, under the assumption that (X,Y) are
%@        %   multivariate normal distributed with a correlation of rmi:
%@        %   mi=-1/2*log(1-rmi^2)
%@        rmi=sqrt(1-exp(-2.*mi));
%@    end
%@
%@end
%@
%@
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __slidingwindow.m__ -->
%@function [center, links, rechts] = slidingwindow(t_min,t_max,windw,overlap)
%@%SLIDINGWINDOW gives overlapping partition windows for time series.
%@%   [center, links, rechts] = slidingwindow(t_min,t_max,windw,overlap)
%@%   gives the window partitions of a time series from t_min to t_max,
%@%   provided the windows have the width windw(>0) and a given overlap
%@%   (0-1).
%@%   EXAMPLE:
%@%   t_min=0;t_max=1000;% start and end point
%@%   windw=300;% 100 year windows
%@%   overlap=0.8;% 80% overlap
%@%   [center, links, rechts] = slidingwindow(t_min,t_max,windw,overlap)
%@%   % produces 3 possible
%@%   % where center gives the centerpoints of the windows, vector links the left
%@%   % boundary and rechts the right boundaries.
%@%
%@% see also: EXTRACTWINDOW,SIMGRAM
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 1.0  last rev: 2013-08-22
%@
%@
%@if any([overlap<0,windw<0,t_max<t_min]),
%@    error('Check the inputs to function!')
%@end
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@% determine the total overlapping portion
%@T=t_max-t_min;
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@% calculate the number of sections of the given width and overlap
%@if overlap<eps,
%@    no_overl=true;
%@    %just start at first timepoint,
%@    lmax=floor((t_max-t_min)/windw);
%@else
%@    no_overl=false;
%@    un=(1-overlap).*windw;
%@    ob=floor(T-windw);
%@    lmax=floor(ob./un +1); %number of shifts with overlap and window width
%@end
%@
%@center=nan(lmax,1);
%@links=nan(lmax,1);
%@rechts=nan(lmax,1);
%@
%@% find the starting point, the first section
%@for l=1:lmax,
%@    if l==1 % find the starting point of the overlapping part
%@        links(1)= t_min;
%@        rechts(1)=links(1)+windw;
%@        
%@    else
%@        if ~no_overl,
%@            links(l)=windw*(1-overlap)+links(l-1);
%@            rechts(l)=links(l)+windw;
%@        else
%@            links(l)=rechts(l-1)+1;
%@            rechts(l)=links(l)+windw;
%@        end
%@    end
%@    center(l)=links(l)+abs(links(l)-rechts(l))/2;
%@    
%@end
%@
%@end
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __tar1.m__ -->
%@function [x y]=tar1(tx,ty, phi,KAPPA,lag, interp_method,thresh)
%@% TAR1 simulates irregular time series of threshold autoregressive processes.
%@%   [X,Y]=TAR1(tx,ty,phi,KAPPA,lag,interp_method) returns observations X and
%@%   Y for the observation time vectors tx and ty (monotonically increasing
%@%   column vectors). Phi (double value between 0 and 1) determines the
%@%   autocorrelation in X, KAPPA is the coupling strength above and below
%@%   the threshold (two double values between 0 and 1 and lag (integer
%@%   value) determines the lag at which the time series are coupled. The
%@%   autocorrelated time series are first generated on a regular timescale
%@%   at very high resolution and subsequently re-sampled by interpolation
%@%   with the method interp_method ('linear','spline') to the timescales tx
%@%   and ty.
%@%   The input vectors tx,ty must be column vectors!
%@%   KAPPA: vector with coupling parameters
%@%   THRESH: threshold (1 value for two kappas)
%@%
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@% Definition of the processes:
%@%  X_n=exp((t_n-t_n-1)/Tau *X_n-1+Eps_n=phi*X_{n-1}+Eps_n;
%@%  Y_n=KAPPA(g)*X_(n-lag)+(1-KAPPA(g)*Eps2_n.
%@% Where KAPPA(g)=kappa1 if X_(n-lag)>thresh, KAPPA(g)=kappa2 if <thresh
%@% ==> a one threshold AR model with two linear dependence regimes
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@% Example:
%@%     tx=1:1:100;tx=tx';ty=tx;
%@%     phi=0.1;
%@%     thresh=0;
%@%     KAPPA=[0.9 0.1]; % strong correl. for low values, lower for high
%@%     %values
%@%     lag=2; % expected to be positive for causality!
%@%     interp_method='linear'
%@%     [x y]=tar1(tx,ty, phi, KAPPA, lag, interp_method,thresh)
%@%     plot(tx,x);hold on; plot(ty,y,'r');legend('x', 'y');xlabel('time')
%@%
%@%
%@%   More information on Threshold-AR processes: 
%@%   http://www.ny.frb.org/research/staff_reports/sr87.pdf
%@%   And on this implementation with one threshold and two linear regimes: 
%@%   REFERENCE:
%@%       K. Rehfeld and J. Kurths: Similarity measures for irregular
%@%       and age uncertain time series, submitted to Clim. Past. Discuss,
%@%       2013.
%@%
%@% see also: NEXCF,GMI,EVENTSYNCHRO,SIMILARITY,TAR
%@%
%@% (c) Potsdam Institute for Climate Impact Research 2013
%@% Kira Rehfeld
%@% ver: 1.0  last rev: 2013-08-22
%@
%@
%@% set random stream according to clock - does not (yet) work in Octave
%@%RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(1000*clock)))
%@
%@% set the resolution at which the processes are simulated at 10 times the
%@% actually desired mean resolution
%@dtsur=min(mean(diff(tx)), mean(diff(ty)))./10;
%@
%@t1=min(min(tx), min(ty))-lag; % subtract lag  to compensate transient effect
%@t2=max(max(tx),max(ty));
%@% coupling strengths
%@coupl_strength=exp(log(KAPPA));
%@t=t1:dtsur:t2;
%@t=t';
%@% persistence time for X
%@fi=exp(-dtsur/(-1/log(phi)));
%@
%@%% First step: Create driving process at high resolution (X)
%@Eps1=randn(length(t),1);
%@%xsur=filter(1,d,Eps1');
%@
%@xsur=Eps1;
%@for i=2:length(t)
%@    xsur(i)=fi.*xsur(i-1) +sqrt(1-fi^2)*Eps1(i);
%@end
%@
%@xsur=(xsur-mean(xsur))./std(xsur);
%@
%@%% Second step: Couple the second timeseries
%@Eps2=randn(length(t),1);
%@ysur=Eps2;
%@
%@L=ceil(lag./dtsur);
%@
%@ for k=L+1:length(xsur)
%@     if xsur(k-L)<thresh% if below thresh, use kappa1 
%@         couplthresh(k)=coupl_strength(1);
%@     else % kappa2 if above threshold
%@         couplthresh(k)=coupl_strength(2);
%@
%@     end
%@     ysur(k)=couplthresh(k)*xsur(k-L) +sqrt(1-couplthresh(k)^2)*Eps2(k-L);
%@ end
%@
%@%ysur(L+1:end)=coupl_strength*(sin(xsur(1:end-L))+cos(ysur(1:end-L)))+sqrt(1-coupl_strength^2)*Eps2(L+1:end);
%@
%@
%@ysur=real((ysur-mean(ysur))./std(ysur));
%@
%@%% Third step, interpolate to (irregular) desired time scale
%@
%@y = interp1(t,ysur,ty,interp_method,'extrap');
%@x = interp1(t,xsur,tx,interp_method,'extrap');
%@y=(y-mean(y))./std(y);
%@x=(x-mean(x))./std(x);
%<-- ASCII ends here -->
%<-- ASCII begins here: __tauest_ls.m__ -->
%@function varargout = tauest_ls(tx,x,lambda0)
%@%TAUEST_LS estimates persistence of irregular time series with least squares.
%@% [Tau, fval] = TAUEST_LS(tx,x,lambda0) returns the Ohrnstein-Uhlenbeck
%@% exponential coeffient Tau from the time series {tx,x}, assuming it is a
%@% continuous-time autoregressive processes of first order, with
%@% reference to the temporal spacing lambda0.
%@%
%@% EXAMPLE:
%@%
%@%   % prepare an exemplatory time series
%@%   tx=1:1000;tx=tx';
%@%   x=randn(size(tx));
%@%   fi=0.7;% autocorrelation parameter
%@%   eps=randn(size(tx));
%@%   % make x autocorrelated
%@%   for t=2:length(tx)
%@%   x(t)=fi*x(t-1)+sqrt(1-fi^2)*eps(t);
%@%   end
%@%   % estimate persistence
%@%   lambda0=mean(diff(tx));% sampling time of the process
%@%   [Tau,fval]=tauest_ls(tx,x,lambda0)
%@%   % convert from persistence time to AR1 coefficient
%@%    est_fi=exp(-mean(diff(tx))/(Tau))
%@%
%@%    REFERENCE:
%@%       Ref.: K. Rehfeld, N. Marwan, J. Heitzig, J. Kurths: Comparison of
%@%       correlation analysis techniques for irregularly sampled time series,
%@%       Nonlinear Processes in Geophysics, 18(3), 389-404 (2011).
%@%
%@%    See also SIMILARITY,NEMI,ESF.
%@    
%@% (c) Potsdam Institute for Climate Impact Research 2011-2013
%@% Kira Rehfeld
%@% ver: 1.0   last rev: 2013-08-21
%@
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@% parameters for the least squares fit
%@%options=optimset('lsqnonlin');
%@%options.Display='off';
%@% options = optimset('Algorithm', {'levenberg-marquardt',.005});
%@%%
%@% set the persistence time into context with the sampling time
%@a0=exp(-lambda0);
%@
%@% least squares fit of AR1 processes
%@%[a,fval]=lsqnonlin(@nestedfun,a0,eps,(max(tx)-min(tx)),options);
%@if exist('lsqnonlin','file') % use least-squares estimator if available
%@	%[x,resnorm] = lsqnonlin(fun,x0,lb,ub,options)
%@[a,fval]=lsqnonlin(@nestedfun,a0,eps,(max(tx)-min(tx))); % tested without options?
%@else
%@	error('least squares fitting for octave not working yet');
%@%elseif exist('leasqr','file') % octave
%@%[fval,a,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
%@  %                 leasqr(x,y,p0,@nestedfun,{stol,niter,wt,dp,dFdp,options})
%@  %[fval,a]=leasqr(tx,x,a0,@nestedfunO)
%@end
%@
%@
%@
%@
%@% use leasqr for octave!
%@    function y=nestedfun(a)
%@        dt=diff(tx)';
%@        Xi=x(2:end);
%@        Xim1=x(1:end-1);
%@        Phi=a.^dt;
%@        % % technically, the variance of the noise term, K, depends on sampling and AR1 parameter.
%@        % % practically, the estimation doesn't converge if this is included.
%@        % % as a rather robust simplification, this is fixed here.
%@        %K=std(x)*sqrt((1-a.^(2*dt))/(-2*log(a)));
%@        K=1;
%@        y=(Xi(:)-Phi(:).*Xim1(:))./K;
%@        
%@    end
%@
%@    %function y=nestedfunO(a,tx)
%@        %dt=diff(tx)';
%@        %Xi=x(2:end);
%@        %Xim1=x(1:end-1);
%@        %Phi=a.^dt;
%@        %% % technically, the variance of the noise term, K, depends on sampling and AR1 parameter.
%@        %% % practically, the estimation doesn't converge if this is included.
%@        %% % as a rather robust simplification, this is fixed here.
%@        %%K=std(x)*sqrt((1-a.^(2*dt))/(-2*log(a)));
%@        %K=1;
%@        %y=(Xi(:)-Phi(:).*Xim1(:))./K;
%@        %
%@    %end
%@% Convert from AR1 coefficient to persistence time.
%@Tau=-1./log(a);
%@
%@%assign output
%@if nargout == 1
%@    varargout = {Tau};
%@elseif nargout == 2
%@    varargout(1) = {Tau};
%@    varargout(2) = {fval};
%@end
%@end
%@
%@
%@
%<-- ASCII ends here -->
%<-- ASCII begins here: __test_nest.m__ -->
%@% Test script for NESToolbox containing test cases for usage illustration.
%@% For example, run this script in Matlab or Octave to check if all functions are
%@% working.
%@
%@%% 1st step: create irregular time series for testing purpose
%@% Irregular and systematic sampling schemes can be simulated with
%@% generate_t. Hereby it is possible to control sampling resolution dt,
%@% minimum and maximum times and the method by which the irregularity is
%@% created. Gamma-distributed sampling times can be created, as well as
%@% missing value timescales.
%@dt=1;% inter-obs. time
%@Tmin=1;Tmax=1000;
%@% create the irregular time sampling
%@tx = generate_t(dt,Tmin,Tmax,'gamma',1,1); % last args control skewness of distribution
%@ty = generate_t(dt,Tmin,Tmax,'linear'); %regular timescale for comparison
%@% Original signals
%@phi=0.7;coupl=0.9;coupllag=3;
%@
%@%% Simulate linearly coupled processes and sample irregular time series
%@% Using the (irregular/ regular) time vectors above, coupled processes are
%@% simulated on these time scales by car (linear coupling) or
%@% tar1(nonlinear coupling).
%@
%@% for similarity estimation
%@maxlag=20;
%@
%@%%%%%%%%%%%%%%%%%%%%%%%%%%% car
%@[x,y]=car(tx,ty,phi,coupl,coupllag,'linear');
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@
%@% plot the results
%@figure()
%@
%@% plot time series
%@subplot(3,2,1)
%@plot(tx,x,'b');hold on;
%@plot(ty,y,'r');legend('x','y');xlabel('time');ylabel('a.u.')
%@title('time series: linear coupling')
%@
%@% plot sampling
%@subplot(3,2,2)
%@hist(diff(tx));hold on;
%@stem(1,Tmax,'r')
%@xlabel('sampling time interval')
%@ylabel('frequency')
%@title('sampling of time series')
%@
%@% plot scatterplot
%@subplot(3,2,3)
%@nest_scatter(tx,x,ty,y,1,3,1)
%@title('Weighted scatterplot')
%@
%@% plot cross-similarity
%@subplot(3,2,4)
%@lag=-maxlag:dt:maxlag;
%@[cxy]=similarity(tx,x,ty,y,lag,'gXCF');
%@[cxy_gmi]=similarity(tx,x,ty,y,lag,'gMI');
%@[cxy_esf]=similarity(tx,x,ty,y,lag,'esf');
%@plot(lag,cxy,'r');hold on;
%@plot(lag,cxy_gmi,'g');plot(lag,cxy_esf,'b');legend('gXCF','gMI','ESF');
%@title('Cross-correlation function with gXCF')
%@xlim([-maxlag maxlag]);xlabel('lag');ylabel('\rho')
%@
%@% compute and plot power spectrum 
%@subplot(3,2,5)
%@%use larger no of lags for better sampling of the spectrum 
%@maxlagPSD=fix(0.5*Tmax/dt);lagPSD=-maxlagPSD:dt:maxlagPSD;
%@[cxy_ps]=similarity(tx,x,tx,x,lagPSD,'gXCF');
%@plt=true;
%@[fy,py,h]=pspek(cxy_ps,mean(diff(lagPSD)),plt);
%@title('Power spectrum of X')
%@set(h,'Color','r')
%@xlim([0 500])
%@
%@
%@%% Simulate nonlinear processes on irregular time vector
%@% This threshold-AR model used in tar1 has two different regimes. If the
%@% value in X at time t-lag is above the threshold, the value in Y at time t
%@% is obtained from positive coupling, if it is below it is obtained from
%@% negative coupling. This is where the limits of cross correlation lie!
%@
%@%%%%%%%%%%%%%%%%%%%%%%%%%%% tar1
%@thresh=0;
%@tx2 = generate_t(dt,Tmin,Tmax,'gamma',1,2); % more irregular
%@ty2 = generate_t(dt,Tmin,Tmax,'linear');
%@[x2,y2]=tar1(tx2,ty2,phi,coupl*[-1 1],coupllag,'linear',thresh);
%@%%%%%%%%%%%%%%%%%%%%%%%%%%%
%@
%@% plot the results
%@figure()
%@
%@% plot time series
%@subplot(2,2,1)
%@plot(tx2,x2,'b');hold on;
%@plot(ty2,y2,'r');legend('x2','y2');xlabel('time');ylabel('a.u.')
%@title('time series: nonlinear threshold AR coupling')
%@xlim([0 200])
%@
%@% plot sampling
%@subplot(2,2,2)
%@hist(diff(tx2));hold on;
%@stem(1,Tmax,'r')
%@xlabel('sampling time interval')
%@ylabel('frequency')
%@title('sampling of time series')
%@
%@% plot scatterplot
%@subplot(2,2,3)
%@nest_scatter(tx2,x2,ty2,y2,1,3,1)
%@title('Weighted scatterplot')
%@
%@% plot cross-similarity
%@subplot(2,2,4)
%@lag=-maxlag:dt:maxlag;
%@clear cxy
%@[cxy(:,1)]=similarity(tx2,x2,ty2,y2,lag,'gMI');
%@[cxy(:,2)]=similarity(tx2,x2,ty2,y2,lag,'ESF');
%@[cxy(:,3)]=similarity(tx2,x2,ty2,y2,lag,'gXCF');
%@plot(lag,cxy);legend('gMI','ESF','gXCF')
%@xlabel('lag');ylabel('\rho_{xy}')
%@title('Cross-Similarity')
%@
%@
%@
%@
%@%% Check the similarity as it changes over time with SIMGRAM
%@% first we create and sample two processes that are coupled, for the first
%@% part linearly, then nonlinearly (to see the difference between
%@% linear/nonlinear methods).
%@clear TX X TY Y TX1 TX2 TY1 TY2 C1 C2 L1 L2
%@
%@windw=5000; % length of the ts
%@part_windw=500; %
%@TX1=generate_t(1,1,windw,'gamma',2,1);
%@TY1=generate_t(1,1,windw,'gamma',1,2);
%@[X1 Y1]=car(TX1,TY1,phi,coupl,coupllag,'linear');
%@
%@TX2=generate_t(1,1,windw,'gamma',2,1)+windw;
%@TY2=generate_t(1,1,windw,'gamma',1,2)+windw;
%@[X2 Y2]=tar1(TX2,TY2,phi,coupl*[1 -1],coupllag,'linear',0);
%@TX=[TX1;TX2];TY=[TY1;TY2];X=[X1;X2];Y=[Y1;Y2];
%@
%@% Then check, how MI reflects the coupling (the strength of which is
%@% constant throughout the whole time period:
%@lag=-10:1:10;
%@figure('Name','Similarity estimates over time')
%@subplot(2,1,1)
%@%gMI
%@[C1,L1,T1] = simgram(TX,X,TY,Y,'overlap',0,'windw',part_windw,'lag',lag,'method','gMI','verbose',1);
%@imagesc(T1,L1,C1);title('Sliding window gMI');xlabel('time'); ylabel('lag');
%@
%@hcb=colorbar;colormap(jet)
%@set(hcb,'YTick',[-1 0 1]);set(hcb,'YTicklabel',{'-1','0','1'});
%@
%@% and how gXCF reflects it
%@subplot(2,1,2)
%@[C2,L2,T2] = simgram(TX,X,TY,Y,'overlap',0,'windw',part_windw,'lag',lag,'method','gXCF');
%@imagesc(T2,L2,C2);title('Sliding window gXCF');xlabel('time'); ylabel('lag');
%@
%@hcb=colorbar;colormap(jet)
%@set(hcb,'YTick',[0 1]);set(hcb,'YTicklabel',{'0','-1'});
%@
%@% The first part of the time series follows a symmetric linear coupling
%@% which is picked up both by gMI and gXCF. In the second half, only gMI
%@% returns visible similarity at the lag of coupling.
%@
%@%% And now to check the significance of the similarity over time,
%@% we do this by creating autocorrelated, but mutually uncorrelated,
%@% surrogates for the time series {TX,X}/{TY,Y} with the function AR1Sur.
%@% Then we estimate their similarity for the same lags/time vectors.
%@
%@NoSur=100; % number of surrogates for MC analysis. For most reliable results, set to 2000.
%@XSur=ar1sur(TX,X,NoSur); % create AR1 Surrogates for X/Y
%@YSur=ar1sur(TY,Y,NoSur);
%@
%@% compute similarities in sverliding windows
%@h=waitbar(0,'Surrogate tests');
%@for k=1:NoSur
%@    waitbar(k/NoSur)
%@    %gMI
%@    [temp1] = simgram(TX,XSur(:,k),TY,YSur(:,k),'overlap',0,'windw',part_windw,'lag',lag,'method','gMI');
%@    %gXCF
%@    [temp2] = simgram(TX,XSur(:,k),TY,YSur(:,k),'overlap',0,'windw',part_windw,'lag',lag,'method','gXCF');
%@    Csur1(k,:,:)=temp1;
%@    Csur2(k,:,:)=temp2;
%@    clear temp1 temp2
%@end
%@delete(h)
%@
%@%% then I perform a two-sided test for zero correlation using the critical
%@% values from the surrogates. Everything that is above the critical value
%@% is considered significant similarity
%@
%@% - quantiles for each point from Csur to obtain critical values
%@
%@Q1=quantile(Csur1,[.05 .95],1);
%@
%@
%@% % lower threshold for each time window/lag (only useful for XCF!)
%@% imagesc(squeeze(Q1(1,:,:)))
%@% % upper threshold for each lag/time window
%@% imagesc(squeeze(Q1(2,:,:)))
%@% % - threshold the original estimates by cvs ... for MI, only with upper
%@% % limit
%@CMbinary=[1 1 1;0 0 0]; % binary black/white colormap
%@colormap(CMbinary)
%@
%@figure()
%@subplot(2,1,1)
%@imagesc(T1,lag,C1>squeeze(Q1(2,:,:)))
%@title('significant gMI');xlabel('time'); ylabel('lag');
%@hcb=colorbar;colormap(CMbinary)
%@set(hcb,'YTick',[0 1]);set(hcb,'YTicklabel',{'Yes','No'});
%@
%@subplot(2,1,2)
%@Q2=quantile(Csur2,[.05,.95],1);
%@% correlation ('gXCF') can have positive or negative values -> significant
%@% correlation or anticorrelation is considered a similarity. Therefore, a
%@% two-sided tests for zero correlation is employed, where both quantiles
%@% are used:
%@imagesc(T1,lag,(C2>squeeze(Q2(2,:,:)))|(C2<squeeze(Q2(1,:,:))))
%@title('Significant gXCF'); xlabel('time');ylabel('lag')
%@colormap(CMbinary)
%@hcb=colorbar;
%@set(hcb,'YTick',[0 1],'YTicklabel',{'Yes','No'})
%@
%@
%@%% Another example: Suppose we have a long time series containing a trend.
%@clear tx x 
%@
%@%The temporal sampling is irregular, so it is not possible to use index for cutting
%@%the time series:
%@tx=generate_t(1,1,1000,'gamma',1,2);
%@% Original signal is random with a trend:
%@x = randn(size(tx))+linspace(0,5,length(tx))';
%@
%@figure()
%@h0=plot(tx,x);xlabel('time');xlim([0 1000]);
%@
%@% we cam use gausssmooth to obtain the trend via a low-pass filter of the time series:
%@bandwidth=50; % bandwidth of the Gaussian kernel smoother
%@x_trend=gausssmooth(tx,x,bandwidth);
%@
%@hold on; h1=plot(tx,x_trend,'r');
%@
%@% We can subtract the trend component to extract the high-frequency
%@% variability (high-pass filter)
%@x_fluct=x-x_trend;
%@h2=plot(tx,x_fluct,'g');
%@legend([h0 h1 h2],'orig','trend','high-frequency')
%@
%@%% If we wanted to compute a statistic in one single time window,
%@% say, between t=500 and t=700, we can use the function extractwindow:
%@links=500;%left border
%@rechts=700; %right border, in time units of tx, not in index units!
%@[txw,xw]=extractwindow(tx,x,links,rechts,0);
%@[txw_fluct,xw_fluct]=extractwindow(tx,x_fluct,links,rechts,0);
%@
%@disp(['the mean in this time window is ' sprintf('%g, ', mean(xw)) 'whereas that of the detrended time series is ' sprintf('%g.',mean(xw_fluct))])
%@
%@% if we wnated sliding, possibly overlapping, windows we can use
%@% extractwindow in combination with slidingwindow 
%@% Here: for 100-timeunits wide windows, overlapping 50%:
%@[center,leftborders,rightborders]=slidingwindow(min(tx),max(tx),100,0.5);
%@
%@% now we can compute the mean in each segment
%@for k=1:length(center)
%@    [txw,xw]=extractwindow(tx,x,leftborders(k),rightborders(k),0);
%@    mean_wind(k)=mean(xw);
%@    std_wind(k)=std(xw);
%@    clear txw xw
%@end
%@% and plot it over the trend series
%@hold on;
%@hout=errorbar(center,mean_wind,std_wind);
%@set(hout,'linewidth',2,'color','k','linestyle','none')
%@
%@%% Create age uncertain realizations of layer counted data
%@clear tx x NoEns
%@tx=1:500;tx=tx'; % create sample time vector
%@x=randn(size(tx));
%@% example for percentwise error (increasing with depth)
%@NoEns=100;
%@[LCrec1] = PANlayerchronensemble(x,tx,NoEns,'%',1);
%@% example for absolute counting error error (i.e. +/- 5 years)
%@[LCrec] = PANlayerchronensemble(x,tx,NoEns,'abs',1);
%@xmat=repmat(x,1,NoEns); % create
%@subplot(2,1,1)
%@plot(LCrec(:,2:10),xmat(:,2:10),'k');hold on;
%@plot(LCrec1(:,2:10),xmat(:,2:10),'g');xlim([0 500])
%@ylabel('pseudoproxy measurement value');
%@xlabel('time')
%@subplot(2,1,2)
%@plot(1:length(tx),var(LCrec(:,2:end),0,2),'k'); hold on;
%@plot(1:length(tx),var(LCrec1(:,2:end),0,2),'g');
%@xlabel('time');ylabel('variance (time)');xlim([0 500])
%@legend('percentwise age error increase','absolute layer counting error')
%@
%@%% If you got this far without an error - great! Not? Let's try to fix it.
%@% Feedback to krehfeld@awi.de is highly welcome :)
%<-- ASCII ends here -->