From 993eec07508b7bb1d4c8a903d5c0410ea6ea86f1 Mon Sep 17 00:00:00 2001 From: isaacengel Date: Fri, 10 Dec 2021 15:39:02 +0100 Subject: [PATCH] Tidying up --- core/applyCovConst.m | 8 - core/fromSH.m | 2 +- core/getYinv.m | 63 +++++++ core/toSH.m | 55 +++---- core/toSH_BiMagLS.m | 4 +- core/toSH_FDTA.m | 10 +- core/toSH_MagLS.m | 11 +- core/toSH_MagSun.m | 57 +++++++ core/toSH_SpSub.m | 15 +- core/toSH_SpSubMod.m | 4 +- core/toSH_TA.m | 10 +- thirdparty/AKtools/AKcolormaps.m | 84 ++++++++++ thirdparty/AKtools/AKcolormapsBrewer.m | 52 ++++++ thirdparty/AKtools/AKpCoordinateTransform.m | 79 +++++++++ thirdparty/AKtools/AKpDataInterpolation.m | 92 +++++++++++ thirdparty/AKtools/AKpDrawCoordinates.m | 88 ++++++++++ thirdparty/AKtools/AKpMulti.m | 114 +++++++++++++ .../AKtools/AKpSetColormapAndColorbar.m | 155 ++++++++++++++++++ thirdparty/cbrewer/cbrewer.m | 128 +++++++++++++++ thirdparty/cbrewer/cbrewer_licence.txt | 27 +++ thirdparty/cbrewer/cbrewer_preview.jpg | Bin 0 -> 92995 bytes thirdparty/cbrewer/change_jet.m | 64 ++++++++ thirdparty/cbrewer/colorbrewer.mat | Bin 0 -> 4475 bytes thirdparty/cbrewer/interpolate_cbrewer.m | 36 ++++ thirdparty/cbrewer/plot_brewer_cmap.m | 50 ++++++ utils/autoreg_minphase.m | 23 +-- utils/getRealSHmatrix.m | 41 +++++ utils/getSphDist.m | 7 + utils/hrtf2sofa.m | 3 +- utils/sofa2hrtf.m | 5 +- 30 files changed, 1197 insertions(+), 90 deletions(-) create mode 100644 core/getYinv.m create mode 100644 core/toSH_MagSun.m create mode 100644 thirdparty/AKtools/AKcolormaps.m create mode 100644 thirdparty/AKtools/AKcolormapsBrewer.m create mode 100644 thirdparty/AKtools/AKpCoordinateTransform.m create mode 100644 thirdparty/AKtools/AKpDataInterpolation.m create mode 100644 thirdparty/AKtools/AKpDrawCoordinates.m create mode 100644 thirdparty/AKtools/AKpMulti.m create mode 100644 thirdparty/AKtools/AKpSetColormapAndColorbar.m create mode 100644 thirdparty/cbrewer/cbrewer.m create mode 100644 thirdparty/cbrewer/cbrewer_licence.txt create mode 100644 thirdparty/cbrewer/cbrewer_preview.jpg create mode 100644 thirdparty/cbrewer/change_jet.m create mode 100644 thirdparty/cbrewer/colorbrewer.mat create mode 100644 thirdparty/cbrewer/interpolate_cbrewer.m create mode 100644 thirdparty/cbrewer/plot_brewer_cmap.m create mode 100644 utils/getRealSHmatrix.m create mode 100644 utils/getSphDist.m diff --git a/core/applyCovConst.m b/core/applyCovConst.m index cef0e92..48152d9 100644 --- a/core/applyCovConst.m +++ b/core/applyCovConst.m @@ -111,14 +111,6 @@ end fprintf('\n') -% Hnm_match2(:,:,1) = Hnm(:,:,1).*G1(:,1) + Hnm(:,:,2).*G2(:,1); -% Hnm_match2(:,:,2) = Hnm(:,:,2).*G1(:,2) + Hnm(:,:,1).*G2(:,2); -% -% Hnm_match3 = Hnm.*permute(G1,[1,3,2]); -% Hnm_match3(:,:,1) = Hnm_match3(:,:,1) + Hnm(:,:,2).*G2(:,1); -% Hnm_match3(:,:,2) = Hnm_match3(:,:,2) + Hnm(:,:,1).*G2(:,2); - - Hnm = Hnm_match; %% Plots diff --git a/core/fromSH.m b/core/fromSH.m index 62c6c5a..66fc375 100644 --- a/core/fromSH.m +++ b/core/fromSH.m @@ -39,7 +39,7 @@ kr = k*r; %% Process -Y = AKsh(N, [], az*180/pi, el*180/pi, 'real').'; % SH coefficients +Y = getRealSHmatrix(N,az,el); % SH coefficients Hnm = ffth(hnm,nfft,1); % to frequency domain H = mult3(Hnm,Y); % interpolate if isaligned % re-align if specified diff --git a/core/getYinv.m b/core/getYinv.m new file mode 100644 index 0000000..e3283cf --- /dev/null +++ b/core/getYinv.m @@ -0,0 +1,63 @@ +function Y_inv = getYinv(N,az,el,w,reg_eps) +% Get transformation matrix to obtain the SH coefficients of an HRTF. +% Available options: +% 1. Use provided integration weights (see w) +% 2. Use pseudoinverse method (least squares solution) +% 3. Use Tikhonov regularisation which is useful if there are missing +% directions, e.g. at the bottom of the sphere [1] +% +% SIMPLE USAGE EXAMPLE: +% Y_inv = getYinv(15,az,el); +% +% INPUT: +% N = target SH order +% az = HRIR azimuth (ndirs x 1) in rad +% el = HRIR elevation (ndirs x 1) in rad (0=top, pi/2=front) +% w = quadrature weights (ndirs x 1); if empty, use pseudoinverse +% reg_eps = epsilon used for Tikhonov regularisation according to [1]; if +% empty, set to 0 (don't regularise) +% +% OUTPUT: +% Y_inv = transformation matrix (ndirs x (N+1)^2) +% +% REFERENCES: +% [1] Duraiswaini, R., Dmitry N. Zotkin, and Nail A. Gumerov. +% "Interpolation and range extrapolation of HRTFs." 2004 IEEE +% International Conference on Acoustics, Speech, and Signal +% Processing. Vol. 4. IEEE, 2004. +% +% AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk +% August 2021 + +if ~exist('w','var') + w = []; +end + +if ~exist('reg_eps','var') + reg_eps = 0; % if epsilon not provided, don't regularise +end + +Y = getRealSHmatrix(N,az,el); + +if ~isempty(w) && reg_eps == 0 + + % If integration weights are provided, use them +% Y_inv = 4*pi*w.*Y'; % not compatible with old Matlab versions + Y_inv = mult2(4*pi*w,Y'); + +else + + if reg_eps == 0 + + % If no regularisation requested, just apply the pseudoinverse + Y_inv = pinv(Y); + + else + + % If regualarisation requested, apply it according to [12] + D = (1 + N*(N+1)) * eye((N+1)^2); % reg. matrix + Y_inv = Y'*(Y*Y'+reg_eps*D)^(-1); % regularised weights + + end +end + diff --git a/core/toSH.m b/core/toSH.m index 5c88133..041aeb7 100644 --- a/core/toSH.m +++ b/core/toSH.m @@ -65,6 +65,8 @@ % dualBand = whether to use dual-band tapering (def=false) % Hnm_ref = precomputed high-order Hnm used to calculate diffuse field % EQ. To speed things up if processing many HRTFs +% reg_eps = epsilon used for Tikhonov regularisation according to [12] +% (def=0=don't regularise) % % OUTPUTS: % hnm = HRIRs' SH coefficients (nfftOut x (N+1)^2 x 2 ears) @@ -119,16 +121,14 @@ % [11] Archontis Politis, Microphone array processing for parametric % spatial audio techniques, 2016 Doctoral Dissertation, Department of % Signal Processing and Acoustics, Aalto University, Finland. +% [12] Duraiswaini, R., Dmitry N. Zotkin, and Nail A. Gumerov. +% "Interpolation and range extrapolation of HRTFs." 2004 IEEE +% International Conference on Acoustics, Speech, and Signal +% Processing. Vol. 4. IEEE, 2004. % % AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk % February 2021 -% TODO: implement Tikhonov regularization for irregular grids -% According to Duraiswami 2004, it goes like this: -% reg_eps = 1e-6; % epsilon -% D = (1 + N*(N+1)) * eye((N+1)^2); % Tikhonov regularization matrix -% Y_inv = Y' * (Y * Y' + reg_eps*D)^(-1); % this instead of pseudoinverse - %% Parse inputs p = inputParser; addParameter(p,'mode','Trunc',@ischar) @@ -152,6 +152,7 @@ addParameter(p,'dualBand',false,@isscalar) addParameter(p,'tapering',0,@isscalar) addParameter(p,'Hnm_ref',[]) +addParameter(p,'reg_eps',0,@isscalar) parse(p,varargin{:}) mode = p.Results.mode; @@ -175,6 +176,7 @@ dualBand = p.Results.dualBand; tapering = p.Results.tapering; Hnm_ref = p.Results.Hnm_ref; +reg_eps = p.Results.reg_eps; if isempty(fc) && ~contains(mode,'BiMagLS','IgnoreCase',1) fc = N*c/(2*pi*r); % if fc not provided, use aliasing frequency @@ -226,15 +228,9 @@ alignOption = 1; % 1=phase delay, 2=onset detection, 0=nothing if alignOption == 1 % Option 1: phase delay H = ffth(h,nfft,1); % to frequency domain - Y0 = AKsh(0, [], az*180/pi, el*180/pi, 'real').'; - if ~isempty(w) -% Y0_inv = 4*pi*w.*Y0'; % if integrations weights are provided, use them - Y0_inv = mult2(4*pi*w,Y0'); % use integrations weights if provided - else - Y0_inv = pinv(Y0); % if not, the pseudoinverse will do just fine - end + Y0_inv = getYinv(0,az,el,w,reg_eps); H0 = mult3(H,Y0_inv); % 0th order SH signal -% pd = -unwrap(angle(H0(2:end,1,:)))./(2*pi*f(2:end)); % phase delay in s +% pd = -unwrap(angle(H0(2:end,1,:)))./(2*pi*f(2:end)); % phase delay % pd_mean = mean(pd,'all'); % avg delay % H = H.*exp(1i*2*pi*f*pd_mean); % subtract delay pd = div2(-unwrap(angle(H0(2:end,1,:))),(2*pi*f(2:end))); % phase delay @@ -244,7 +240,6 @@ onsL = AKonsetDetect(h(:,:,1)); % left ear onsets onsR = AKonsetDetect(h(:,:,2)); % right rear onsets ons = min([onsL;onsR]); % ipsilateral ear onsets - % nshift = -floor(min(ons))); % detect earliest onset nshift = -round(median(ons)); % median onset h = circshift(h,nshift); H = ffth(h,nfft,1); @@ -255,30 +250,30 @@ %% Preprocess using the indicated method if strcmpi(mode,'Trunc') % just order truncation - % Order-truncated SH-HRTF (faster than order Nmax + truncation) - Y = AKsh(N,[],az*180/pi,el*180/pi,'real').'; - if ~isempty(w) -% Y_inv = 4*pi*w.*Y'; % use integrations weights if provided - Y_inv = mult2(4*pi*w,Y'); % use integrations weights if provided - else - Y_inv = pinv(Y); % if not, pseudoinverse will do just fine - end + % Order-truncated SH-HRTF (faster than order Nmax + truncation) + Y_inv = getYinv(N,az,el,w,reg_eps); Hnm = mult3(H,Y_inv); elseif strcmpi(mode,'SpSub') % subsampling aka virtual loudspeaker - Hnm = toSH_SpSub(H,N,az,el,w,Nmax); + Hnm = toSH_SpSub(H,N,az,el,w,Nmax,reg_eps); elseif strcmpi(mode,'FDTA') % frequency-dependent time-alignment if fc>=fs/2 warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). FDTA preprocessing has no effect...',fc,fs/2) end - [Hnm,varOut.fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl); + [Hnm,varOut.fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl,reg_eps); elseif strcmpi(mode,'MagLS') % magnitude least squares if fc>=fs/2 warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). MagLS preprocessing has no effect...',fc,fs/2) end - [Hnm,varOut.fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r); + [Hnm,varOut.fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r,reg_eps); + +elseif strcmpi(mode,'MagSun') % magnitude optimisation by Sun + if fc>=fs/2 + warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). MagSun preprocessing has no effect...',fc,fs/2) + end + [Hnm,varOut.fc] = toSH_MagSun(H,N,az,el,fs,w,fc,k,r,reg_eps); elseif strcmpi(mode,'MagLSold') % old smoothing, kept to reproduce old data if fc>=fs/2 @@ -287,13 +282,13 @@ [Hnm,varOut.fc] = toSH_MagLS_old(H,N,az,el,fs,w,fc,k,r); elseif strcmpi(mode,'TA') % time alignment (ear alignment) - Hnm = toSH_TA(H,N,az,el,fs,w,r,earAz,earEl); + Hnm = toSH_TA(H,N,az,el,fs,w,r,earAz,earEl,reg_eps); elseif strcmpi(mode,'BiMagLS') % ear alignment + MagLS if fc>=fs/2 warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). BiMagLS preprocessing has no effect...',fc,fs/2) end - [Hnm,varOut.fc] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl); + [Hnm,varOut.fc] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl,reg_eps); elseif strcmpi(mode,'BiMagLSold') % old smoothing, kept to reproduce old data if fc>=fs/2 @@ -305,7 +300,7 @@ if fc>=fs/2 warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). SpSubMod preprocessing is the same as SpSub...',fc,fs/2) end - [Hnm,varOut.fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax); + [Hnm,varOut.fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax,reg_eps); else error('Unknown method %s. Use one of these: ''Trunc'', ''SpSub'', ''FDTA'', ''MagLS'', ''TA'', ''BiMagLS'', ''SpSubMod''',mode) @@ -364,7 +359,7 @@ if EQ==1 || EQ==2 % these modes require a high-order reference if isempty(Hnm_ref) - Y = AKsh(Nmax,[],az*180/pi,el*180/pi,'real').'; + Y = getRealSHmatrix(Nmax,az,el); if contains(mode,'ear') % time-align reference if required kr = 2*pi*f*r/c; p = earAlign(kr,az,el,earAz,earEl); diff --git a/core/toSH_BiMagLS.m b/core/toSH_BiMagLS.m index 6425e5e..2bc03da 100644 --- a/core/toSH_BiMagLS.m +++ b/core/toSH_BiMagLS.m @@ -1,4 +1,4 @@ -function [Hnm,fc,p] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl) +function [Hnm,fc,p] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl,reg_eps) % Transform HRTF to SH domain at order N by using ear aligning first [1] % and frequency-dependent optimisation above a certain cutoff frequency % (MagLS [2,3]). The idea is to align the HRIRs first to reduce the @@ -83,5 +83,5 @@ H = H.*exp(-1i*p); % apply correction %% Then, apply MagLS -[Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r); +[Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r,reg_eps); diff --git a/core/toSH_FDTA.m b/core/toSH_FDTA.m index f5674ad..02dc49b 100644 --- a/core/toSH_FDTA.m +++ b/core/toSH_FDTA.m @@ -1,4 +1,4 @@ -function [Hnm,fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl) +function [Hnm,fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl,reg_eps) % Transform HRTF to SH domain at order N using frequency-dependent % time-alignment [1]. Alignment is performed with the method in [2]. % @@ -60,12 +60,6 @@ H(ind,:,:) = H(ind,:,:).*exp(-1i*p(ind,:,:)); % apply correction %% Then, get the order-truncated SH-HRTF -Y = AKsh(N,[],az*180/pi,el*180/pi,'real').'; -if exist('w','var') && ~isempty(w) -% Y_inv = 4*pi*w.*Y'; % if integrations weights are provided, use them - Y_inv = mult2(4*pi*w,Y'); % use integrations weights if provided -else - Y_inv = pinv(Y); % if not, the pseudoinverse will do just fine -end +Y_inv = getYinv(N,az,el,w,reg_eps); Hnm = mult3(H,Y_inv); diff --git a/core/toSH_MagLS.m b/core/toSH_MagLS.m index d29f5fa..f960611 100644 --- a/core/toSH_MagLS.m +++ b/core/toSH_MagLS.m @@ -1,4 +1,4 @@ -function [Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r) +function [Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r,reg_eps) % Transform HRTF to SH domain at order N with the MagLS method [1], % following the simplified approach from [2] sec. 4.11.2 and 4.11.3. The % main difference of this implementation vs [1] is that the phase is @@ -56,13 +56,8 @@ f = linspace(0,fs/2,nfreqs).'; % frequency vector %% Get the SH matrices -Y = AKsh(N, [], az*180/pi, el*180/pi, 'real').'; -if exist('w','var') && ~isempty(w) -% Y_inv = 4*pi*w.*Y'; % if integrations weights are provided, use them - Y_inv = mult2(4*pi*w,Y'); % use integrations weights if provided -else - Y_inv = pinv(Y); % if not, the pseudoinverse will do just fine -end +Y = getRealSHmatrix(N,az,el); +Y_inv = getYinv(N,az,el,w,reg_eps); %% Apply simplified MagLS from [2] sec.4.11.2 diff --git a/core/toSH_MagSun.m b/core/toSH_MagSun.m new file mode 100644 index 0000000..75b779c --- /dev/null +++ b/core/toSH_MagSun.m @@ -0,0 +1,57 @@ +function [Hnm,fc] = toSH_MagSun(H,N,az,el,fs,w,fc,k,r,reg_eps) +% Transform HRTF to SH domain at order N with the Magnitude-optimised +% reconstruction method from Sun's thesis [1]. +% +% SIMPLE USAGE EXAMPLE: +% Hnm = toSH_MagSun(H,15,az,el,48000); +% +% INPUT: +% H = HRTF up to Nyquist frequency (nfreqs x ndirs x 2 ears) +% N = target SH order +% az = HRIR azimuth (ndirs x 1) in rad +% el = HRIR elevation (ndirs x 1) in rad (0=top, pi/2=front) +% fs = sampling frequency in Hz +% w = quadrature weights (ndirs x 1); if empty, use pseudoinverse +% fc = cutoff frequency in Hz above which phase is "disregarded" in +% favour of magnitude; if empty (default), use aliasing frequency +% k = length of the transition band in octaves (def=1). E.g. +% if fc=623.89 Hz, smooth between 349.65 and 882.31 Hz. If k==0, +% don't smooth. +% r = head radius in m (def=0.0875) +% +% OUTPUT: +% Hnm = HRTF's SH coefficients (nfreqs x (N+1)^2 x 2 ears) +% fc = see above +% +% REFERENCES: +% [1] Sun, David. "Generation and perception of three-dimensional sound +% fields using Higher Order Ambisonics." (2013). +% +% AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk +% September 2021 + +%% Some parameters +if ~exist('r','var') || isempty(r) + r = 0.0875; % default head radius +end +if ~exist('fc','var') || isempty(fc) + c = 343; % speed of sound (m/s) + fc = N*c/(2*pi*r); % if fc not provided, use aliasing frequency +end +if ~exist('k','var') || isempty(k) + k = 1; % default smoothing = one octave (half to each side) +end +nfreqs = size(H,1); +f = linspace(0,fs/2,nfreqs).'; % frequency vector + +%% Get the SH matrices + +Hmag = abs(H); +Hphase = unwrap(angle(H)); +HphaseAvg = mean(Hphase,2); %HphaseAvg = mult3(Hphase,ones(ndirs))/ndirs; +alpha = 0.5 + k/log(2) * log(f/fc); % 0 for fc1, 1 for fc2 +alpha = min(max(alpha,0),1); % keep between 0 and 1 +HphaseMod = alpha.*HphaseAvg + (1-alpha).*Hphase; +Htarg = Hmag.*exp(1i*HphaseMod); +Y_inv = getYinv(N,az,el,w,reg_eps); +Hnm = mult3(Htarg,Y_inv); diff --git a/core/toSH_SpSub.m b/core/toSH_SpSub.m index ca76b5e..0ec9cec 100644 --- a/core/toSH_SpSub.m +++ b/core/toSH_SpSub.m @@ -1,4 +1,4 @@ -function Hnm = toSH_SpSub(H,N,az,el,w,Nmax) +function Hnm = toSH_SpSub(H,N,az,el,w,Nmax,reg_eps) % Transform HRTF to SH domain at order N using spatial subsampling [1]. % This is equivalent to virtual loudspeaker decoding (first introduced by % [2]) with "mode-matching" (pseudoinverse-based). @@ -35,13 +35,7 @@ %% Process % First, get the SH-HRTF with highest available order -Ymax = AKsh(Nmax, [], az*180/pi, el*180/pi, 'real').'; % SH coeffs -if exist('w','var') && ~isempty(w) -% Ymax_inv = 4*pi*w.*Ymax'; % if integrations weights are provided, use them - Ymax_inv = mult2(4*pi*w,Ymax'); % use integrations weights if provided -else - Ymax_inv = pinv(Ymax); % if not, the pseudoinverse will do just fine -end +Ymax_inv = getYinv(Nmax,az,el,w,reg_eps); Hnm_max = mult3(H,Ymax_inv); % Then, resample to an appropriately sized sparse grid (Gauss) @@ -50,8 +44,7 @@ Hs = mult3(Hnm_max,Ymax_s); % = cat(3,Hnm_max(:,:,1)*Ymax_s,Hnm_max(:,:,2)*Ymax_s); % Finally, convert to SH again at the target low order -Ys = AKsh(N, [], sgrid(:,1)*180/pi, sgrid(:,2)*180/pi, 'real').'; -% Ys_inv = 4*pi*diag(sgrid(:,3))*Ys'; -Ys_inv = mult2(4*pi*sgrid(:,3),Ys'); +Ys_inv = getYinv(N,sgrid(:,1),sgrid(:,2)); Hnm = mult3(Hs,Ys_inv); + \ No newline at end of file diff --git a/core/toSH_SpSubMod.m b/core/toSH_SpSubMod.m index 4454201..a26c18c 100644 --- a/core/toSH_SpSubMod.m +++ b/core/toSH_SpSubMod.m @@ -1,4 +1,4 @@ -function [Hnm,fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax) +function [Hnm,fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax,reg_eps) % Combination of spatial subsampling [1] and frequency-dependent time % alignment [2], as suggested in [3]. % @@ -71,4 +71,4 @@ H(ind,:,:) = H(ind,:,:).*exp(-1i*p(ind,:,:)); % apply correction %% Then, apply spatial subsampling -Hnm = toSH_SpSub(H,N,az,el,w,Nmax); +Hnm = toSH_SpSub(H,N,az,el,w,Nmax,reg_eps); diff --git a/core/toSH_TA.m b/core/toSH_TA.m index d5c509c..c82a87a 100644 --- a/core/toSH_TA.m +++ b/core/toSH_TA.m @@ -1,4 +1,4 @@ -function [Hnm,p] = toSH_TA(H,N,az,el,fs,w,r,earAz,earEl) +function [Hnm,p] = toSH_TA(H,N,az,el,fs,w,r,earAz,earEl,reg_eps) % Transform HRTF to SH domain at order N using time-alignment through phase % correction by ear alignment [1], followed by order truncation. To undo % the alignment after SH interpolation, use the function fromSH with option @@ -50,12 +50,6 @@ H = H.*exp(-1i*p); % apply correction %% Then, get the order-truncated SH-HRTF -Y = AKsh(N,[],az*180/pi,el*180/pi,'real').'; -if exist('w','var') && ~isempty(w) -% Y_inv = 4*pi*w.*Y'; % if integrations weights are provided, use them - Y_inv = mult2(4*pi*w,Y'); % use integrations weights if provided -else - Y_inv = pinv(Y); % if not, the pseudoinverse will do just fine -end +Y_inv = getYinv(N,az,el,w,reg_eps); Hnm = mult3(H,Y_inv); diff --git a/thirdparty/AKtools/AKcolormaps.m b/thirdparty/AKtools/AKcolormaps.m new file mode 100644 index 0000000..c2afea8 --- /dev/null +++ b/thirdparty/AKtools/AKcolormaps.m @@ -0,0 +1,84 @@ +% AKcm = AKcolormaps(AKcm, N) +% creates the N step colormap specified by the string AKcm +% see AKplotDemo.m for examples on how to set colomaps in a plot +% +% I N P U T +% AKcm - string specifiying the desired colormap. +% this can be one of the custom colormaps listed below, or a map +% from the color brewer set - call AKcolormaps() for a list of +% brewer color maps. +% +% Custom color maps +% 'AKbwr': going from blue to white to red (Hagen Wierstorf style) +% 'AKgyr': going from green to yellow to red (The danger scale) +% 'AKwr' : going from white to red +% 'LGBT' : if you like the look of multicolored maps, but 'yet' +% isn't your type you can try 'LGBT', 'LGBTflag', or +% 'LGBTcat' +% +% N - Number of steps / colors in the colormap (default is 128) +% +% O U T P U T +% AKcm - colormap, that is a [N x 3] matrix with rgb values +% +% 06/2016 - fabian.brinkmann@tu-berlin.de, initial dev. +% 10/2016 - helmholz@campus.tu-berlin.de, added brewer colormap support + +% AKtools +% Copyright (C) 2016 Audio Communication Group, Technical University Berlin +% Licensed under the EUPL, Version 1.1 or as soon they will be approved by +% the European Commission - subsequent versions of the EUPL (the "License") +% You may not use this work except in compliance with the License. +% You may obtain a copy of the License at: +% http://joinup.ec.europa.eu/software/page/eupl +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" basis, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. +function AKcm = AKcolormaps(AKcm, N) + +if ~nargin + AKcolormapsBrewer +else + + if ~exist('N', 'var') + N = 128; + end + + switch upper(AKcm) + case 'AKBWR' + AKcm = makeColorMap([.2 .2 1], [1 1 1], [1 .2 .2], N); + case 'AKWR' + AKcm = makeColorMap([1 1 1], [1 .2 .2], N); + case 'AKGYR' + AKcm = makeColorMap([18 159 73]/255, [239 238 36]/255, [238 36 36]/255, N); + case {'LGBT' 'LGBTFLAG' 'LGBTCAT'} + LGBT = [220 38 36 % red + 251 91 19 % orange + 241 221 45 % yellow + 95 185 95 % green + 42 57 187 % blue + 162 42 164]; % ourple + LGBT = flipud(LGBT) / 255; + + x = linspace(1, N, size(LGBT,1)); + + if strfind(upper(AKcm), 'FLAG') %#ok<*STRIFCND> contains is recommended and was introduced in 2016b + AKcm = LGBT; + elseif strfind(upper(AKcm), 'CAT') + AKcm = repmat(LGBT, [ceil(N/size(LGBT,1)) 1]); + AKcm = AKcm(1:N, :); + else + AKcm = interp1(x, LGBT, 1:N, 'linear'); + end + + otherwise + try + brewer_struct = AKcolormapsBrewer(AKcm); + AKcm = cbrewer(brewer_struct, AKcm, N, 'PCHIP'); + catch + error(['AKcolormap ' AKcm ' not defined']) + end + end +end diff --git a/thirdparty/AKtools/AKcolormapsBrewer.m b/thirdparty/AKtools/AKcolormapsBrewer.m new file mode 100644 index 0000000..6fcb66d --- /dev/null +++ b/thirdparty/AKtools/AKcolormapsBrewer.m @@ -0,0 +1,52 @@ +% Only to be called from AKcolormaps. See AKcolomaps for documentation. +% See AKplotDemo.m for examples +% +% 10/2016 - helmholz@campus.tu-berlin.de + +% AKtools +% Copyright (C) 2016 Audio Communication Group, Technical University Berlin +% Licensed under the EUPL, Version 1.1 or as soon they will be approved by +% the European Commission - subsequent versions of the EUPL (the "License") +% You may not use this work except in compliance with the License. +% You may obtain a copy of the License at: +% http://joinup.ec.europa.eu/software/page/eupl +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" basis, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. +function fieldName = AKcolormapsBrewer(color) + +if nargin + + % check if the brewer colormap exists + fieldName = ''; + load('colorbrewer.mat'); + + fields = fieldnames(colorbrewer); + for f = 1:length(fields) + if isfield(colorbrewer.(fields{f}),color) + fieldName = fields{f}; + return; % found name in colorbrewer.(fields{f}) + end + end + + if ~nargout + clear fieldName + end + +else + + % show a plot of all brewer colormaps + if isempty(findall(0,'Type','Figure')) + AKf + end + set(gcf, 'numberTitle', 'off', 'name', 'Brewer colormaps (http://colorbrewer.org/)', 'toolBar', 'none', 'menuBar', 'none') + imagesc(imread(which('cbrewer_preview.jpg'))) + axis off + + if nargout + fieldName = []; + end + +end diff --git a/thirdparty/AKtools/AKpCoordinateTransform.m b/thirdparty/AKtools/AKpCoordinateTransform.m new file mode 100644 index 0000000..c01d490 --- /dev/null +++ b/thirdparty/AKtools/AKpCoordinateTransform.m @@ -0,0 +1,79 @@ +function co = AKpCoordinateTransform(az, el, type) +% function co = hp_coordinate_transformation(az, el, type) +% +% transform different coordinate conventions to matlab default +% if az containes values bigger than 2pi, it is assumed that the unit of az +% is degree. In this case it is converted to radians +% +% See AKplotDemo.m for examples +% +% type +% 1: Matlab default +% 0 <= az < 360 (0=source in front 90 = source to the left of listener) +% 90 <= el <= -90 (90=north pole, 0=front, -90=south pole) +% positive x is at (0|90) (az|el) +% positive y is at (90|90) +% positive z is at (0|0); +% +% 2: Mathematical +% 0 <= az < 360 (0=source in front 90 = source to the left of listener) +% 0 <= el <= 180 (0=north pole, 90=front, 180=south pole) +% positive x is at (0|90) (az|el) +% positive y is at (90|90) +% positive z is at (0|0); + +% AKtools +% Copyright (C) 2016 Audio Communication Group, Technical University Berlin +% Licensed under the EUPL, Version 1.1 or as soon they will be approved by +% the European Commission - subsequent versions of the EUPL (the "License") +% You may not use this work except in compliance with the License. +% You may obtain a copy of the License at: +% http://joinup.ec.europa.eu/software/page/eupl +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" basis, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. + +if max(max(az > 2*pi)) + az = deg2rad(az); + if ~islogical(el) + el = deg2rad(el); + else + el = 0; + end +end + +switch type + case 1 % matlab default, see doc sph2cart + % coordinate transformation + co.az = az; + co.el = el; + % definition of axis normals for balloon plots + co.x = [1 0 0]'; + co.y = [0 1 0]'; + co.z = [0 0 1]'; + % x and y axis for spherical planar plots + co.xtick = 1:45:361; + co.xticklabel = 0:45:360; + co.ytick = 1:45:181; + co.yticklabel = [90:-45:0 -45 -90]; + case 2 % mathematical correct + + % coordinate transformation + co.az = az; + id1 = el > pi/2; + id2 = el <= pi/2; + el(id1) = -(el(id1) - pi/2); + el(id2) = abs(el(id2) - pi/2); + co.el = el; + % definition of axis normals for balloon plots + co.x = [1 0 0]'; + co.y = [0 1 0]'; + co.z = [0 0 1]'; + % x and y axis for spherical planar plots + co.xtick = 1:45:361; + co.xticklabel = 0:45:360; + co.ytick = 1:45:181; + co.yticklabel = 0:45:180; +end diff --git a/thirdparty/AKtools/AKpDataInterpolation.m b/thirdparty/AKtools/AKpDataInterpolation.m new file mode 100644 index 0000000..85e7f46 --- /dev/null +++ b/thirdparty/AKtools/AKpDataInterpolation.m @@ -0,0 +1,92 @@ +function [X, Y, Z, data] = AKpDataInterpolation(co, data) +% interpolates spherical data with an arbitrary grid to a full sphere with +% 1 degree resolution. Works only with the matlab default definition of +% coordinates (see doc sph2cart). Warning: Simple interpolation is +% performed that might distort the data! +% +% Only to be called from AKp.m +% See AKplotDemo.m for examples +% +% fabian.brinkmann@tu-berlin.de, Audio Communication Group TU Berlin, +% DFG research unit 'SEACEN', 7/2012 + +% AKtools +% Copyright (C) 2016 Audio Communication Group, Technical University Berlin +% Licensed under the EUPL, Version 1.1 or as soon they will be approved by +% the European Commission - subsequent versions of the EUPL (the "License") +% You may not use this work except in compliance with the License. +% You may obtain a copy of the License at: +% http://joinup.ec.europa.eu/software/page/eupl +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" basis, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either expressed or implied. +% See the License for the specific language governing permissions and +% limitations under the License. + +% transpose az and el to obtain row vectors +if size(co.az,1) > size(co.az,2) + co.az = co.az'; +end +if size(co.el,1) > size(co.el,2) + co.el = co.el'; +end + +% show original and interpolated data (for debugging only) +do_plot = 0; + +if do_plot + scatter3(co.az, co.el, data) %#ok<*UNRCH> +end + +% the interpolation is done one 3D plane, where az is x axis, el is y axis +% and data is z - axis. For propper interpolation close to the edges of the +% plane, duplicate version of the original data has to be moved to the left +% and right and mirrored above and below. This is done with the full data, +% to make sure that it works in any case. It therefore might be a bit +% slower. This means we need to add 8 surrounding versions of the original +% data. + +% left and right +az_i = [co.az co.az-2*pi co.az+2*pi]; +el_i = [co.el co.el co.el]; +data_i = [data data data]; + +% mirror everthing to the top +id1 = el_i >= 0; +id2 = el_i < 0; + +el1 = el_i; +el1(id1) = el1(id1) + 2*(pi/2-el1(id1)); +el1(id2) = abs(el1(id2)) + pi; + +% mirror everything to the botton +el2 = el_i; +el2(id1) = -el2(id1) - pi; +el2(id2) = el2(id2) -2*(pi/2+el2(id2)); + +% insert mirrored versions +az_i = [az_i az_i az_i]; +el_i = [el_i el1 el2]; +data_i = [data_i data_i data_i]; + +% get interpolator object +warning('off', 'MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId') +data_interpolator = scatteredInterpolant(az_i',el_i',data_i'); +warning('on', 'MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId') + +% get grid for interpolation +az = (0:360)*pi/180; +el = (90:-1:-90)*pi/180; +[az, el] = meshgrid(az, el); + +% get interpolated data +data = data_interpolator(az, el); + +% transform to cartesian coordinates +[X, Y, Z] = sph2cart(az, el, data); + +if do_plot + hold on + mesh(az, el, data) + hold off +end diff --git a/thirdparty/AKtools/AKpDrawCoordinates.m b/thirdparty/AKtools/AKpDrawCoordinates.m new file mode 100644 index 0000000..b6862fd --- /dev/null +++ b/thirdparty/AKtools/AKpDrawCoordinates.m @@ -0,0 +1,88 @@ +function AKpDrawCoordinates(co, type, marker_size, line_width) +% AKpDrawCoordinates(co, type, marker_size, line_width) +% draws the coordinate axis in 3D plots. +% +% See AKplotDemo.m for examples +% +% I N P U T +% co - strct with fields x, y,z (unit vectors pointing in positive +% x, y, and z direction) and L (scalar that gives the axis +% length) +% type - 1: lines are drawn, marker indicates positive x,y and z +% 2: additional circles and lines mark intermediate angles +% marker_size - Size of markers that show positive x, y, and z direction +% line_width - Line width of x, y, and z axis +% +% fabian.brinkmann@tu-berlin.de, Audio Communication Group TU Berlin, +% DFG research unit 'SEACEN', 7/2012 + +% AKtools +% Copyright (C) 2016 Audio Communication Group, Technical University Berlin +% Licensed under the EUPL, Version 1.1 or as soon they will be approved by +% the European Commission - subsequent versions of the EUPL (the "License") +% You may not use this work except in compliance with the License. +% You may obtain a copy of the License at: +% http://joinup.ec.europa.eu/software/page/eupl +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" basis, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. + +if ~exist('marker_size', 'var') + % this is the size for the marker that indicates positive x,y and z + marker_size = 10; +end + +colors = AKcolors; + +switch type + case 1 + hold on + % draw x axis and point for markin positive x + plot3([-co.x(1)*co.L co.x(1)*co.L], ... + [-co.x(2)*co.L co.x(2)*co.L], ... + [-co.x(3)*co.L co.x(3)*co.L], 'color', colors.rgb(4,:), 'LineWidth', line_width) + plot3(co.x(1)*co.L, co.x(2)*co.L, co.x(3)*co.L, '.', 'color', colors.rgb(4,:), 'Markersize', marker_size) + % draw y axis and point for markin positive y + plot3([-co.y(1)*co.L co.y(1)*co.L], ... + [-co.y(2)*co.L co.y(2)*co.L], ... + [-co.y(3)*co.L co.y(3)*co.L], 'color', colors.rgb(7,:), 'LineWidth', line_width) + plot3(co.y(1)*co.L, co.y(2)*co.L, co.y(3)*co.L, '.', 'color', colors.rgb(7,:), 'Markersize', marker_size) + % draw z axis and point for markin positive z + plot3([-co.z(1)*co.L co.z(1)*co.L], ... + [-co.z(2)*co.L co.z(2)*co.L], ... + [-co.z(3)*co.L co.z(3)*co.L], 'color', colors.rgb(3,:), 'LineWidth', line_width) + plot3(co.z(1)*co.L, co.z(2)*co.L, co.z(3)*co.L, '.', 'color', colors.rgb(3,:), 'Markersize', marker_size) + hold off + case 2 + hold on + % draw circles + x = linspace(0, 2*pi, 256); + y = zeros(size(x)); + plot3(cos(x)*co.L, sin(x)*co.L, y, 'color', [.7 .7 .7], 'LineWidth', line_width) + plot3(cos(x)*co.L, y, sin(x)*co.L, 'color', [.7 .7 .7], 'LineWidth', line_width) + plot3(y, cos(x)*co.L, sin(x)*co.L, 'color', [.7 .7 .7], 'LineWidth', line_width) + % draw angles + for aa = [22.5 45 67.5 112.5 135 157.5] + plot3(cosd(aa)*[co.L -co.L], sind(aa)*[co.L -co.L], [0 0], 'color', [.7 .7 .7], 'LineWidth', line_width) + plot3(cosd(aa)*[co.L -co.L], [0 0], sind(aa)*[co.L -co.L], 'color', [.7 .7 .7], 'LineWidth', line_width) + plot3([0 0], cosd(aa)*[co.L -co.L], sind(aa)*[co.L -co.L], 'color', [.7 .7 .7], 'LineWidth', line_width) + end + % draw x axis and point for markin positive x + plot3([-co.x(1)*co.L co.x(1)*co.L], ... + [-co.x(2)*co.L co.x(2)*co.L], ... + [-co.x(3)*co.L co.x(3)*co.L], 'color', colors.rgb(4,:), 'LineWidth', line_width) + plot3(co.x(1)*co.L, co.x(2)*co.L, co.x(3)*co.L, '.', 'color', colors.rgb(4,:), 'Markersize', marker_size) + % draw y axis and point for markin positive y + plot3([-co.y(1)*co.L co.y(1)*co.L], ... + [-co.y(2)*co.L co.y(2)*co.L], ... + [-co.y(3)*co.L co.y(3)*co.L], 'color', colors.rgb(7,:), 'LineWidth', line_width) + plot3(co.y(1)*co.L, co.y(2)*co.L, co.y(3)*co.L, '.', 'color', colors.rgb(7,:), 'Markersize', marker_size) + % draw z axis and point for markin positive z + plot3([-co.z(1)*co.L co.z(1)*co.L], ... + [-co.z(2)*co.L co.z(2)*co.L], ... + [-co.z(3)*co.L co.z(3)*co.L], 'color', colors.rgb(3,:), 'LineWidth', line_width) + plot3(co.z(1)*co.L, co.z(2)*co.L, co.z(3)*co.L, '.', 'color', colors.rgb(3,:), 'Markersize', marker_size) + hold off +end diff --git a/thirdparty/AKtools/AKpMulti.m b/thirdparty/AKtools/AKpMulti.m new file mode 100644 index 0000000..12918e7 --- /dev/null +++ b/thirdparty/AKtools/AKpMulti.m @@ -0,0 +1,114 @@ +% AKpMulti(data, plot_type, varargin) +% can be used for showing multiple plots of the data in a single figure. +% AKpMulti parses the plot_type and passes the input to AKp, for example +% +% AKpMulti(AKnoise, {'t2d' 'm2d'}) +% plots the time signal and spectrum of a pink noise. AKpMulti can also be +% called using AKp, e.g. +% AKp(AKnoise, {'t2d' 'm2d'}) +% +% See AKplotDemo.m for examples +% +% +% I N P U T +% data - see help AKp +% plot_type - cell array with different plot_type as documented in AKp, or +% string (see predefined plot types for more information) +% The string in each cell determine the current plot type, and +% the size of the cell array determines the subplot layout, +% e.g. +% {'t2d'; 'm2d'} +% plot the impulse response and spectrum in a figure with 2 +% vertically arranges subplots +% {'t2d', 'm2d'} +% plot the impulse response and spectrum in a figure with 2 +% horizontally arranges subplots +% {'t2d', 'm2d'; 'p2d' []} +% plots a figure with a 2x2 subplot matrix, where the last +% subplot is empty +% +% varargin - is passed to AKp. See AKp for documentation +% +% P R E D E F I N E D P L O T T Y P E S +% '1a' : {'t2d' 'et2d'; 'm2d' 'p2d'} +% '1b' : {'t3d' 'et3d'; 'm3d' 'p3d'} +% '2a' : {'t2d'; 'm2d'} +% '2b' : {'t3d'; 'm3d'} +% +% O U T P U T +% dataOut - struct holding the output data from AKp. See AKp for +% documentation +% +% 10/2016 - fabian.brinkmann@tu-berlin.de + +% AKtools +% Copyright (C) 2016 Audio Communication Group, Technical University Berlin +% Licensed under the EUPL, Version 1.1 or as soon they will be approved by +% the European Commission - subsequent versions of the EUPL (the "License") +% You may not use this work except in compliance with the License. +% You may obtain a copy of the License at: +% http://joinup.ec.europa.eu/software/page/eupl +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" basis, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. +function dataOut = AKpMulti(data, plot_type, varargin) + +if ~exist('plot_type', 'var') + plot_type = '1a'; +end + +if iscell(plot_type) % <- user defined plots + + % make figure if none exists + if isempty(findall(0,'Type','Figure')) + AKf + end + + % get dimensions + [N, M] = size(plot_type); + + % run the plots + for nn = 1:N + for mm = 1:M + + % check if current subplot is empty + if ~isempty(plot_type{nn,mm}) + subplot(N, M, (nn-1)*M+mm) + dOut = AKp(data, plot_type{nn,mm}, varargin{:}); + + % copy output data of single AKp calls to output data of + % AKpMulti + f = fieldnames(dOut); + for ff = 1:numel(f) + dataOut((nn-1)*M+mm, 1).(f{ff}) = dOut.(f{ff}); + end + end + + end + end + +else % <- pre defined plots + + % re-define plot_type + switch lower(plot_type) + case '1a' + plot_type = {'t2d' 'et2d'; 'm2d' 'p2d'}; + case '1b' + plot_type = {'t3d' 'et3d'; 'm3d' 'p3d'}; + case '2a' + plot_type = {'t2d'; 'm2d'}; + case '2b' + plot_type = {'t3d'; 'm3d'}; + end + + % run AKpMulti again + dataOut = AKpMulti(data, plot_type, varargin{:}); + +end + +% clear output data if not wanted +if ~nargout + clear dataOut +end diff --git a/thirdparty/AKtools/AKpSetColormapAndColorbar.m b/thirdparty/AKtools/AKpSetColormapAndColorbar.m new file mode 100644 index 0000000..e05db8f --- /dev/null +++ b/thirdparty/AKtools/AKpSetColormapAndColorbar.m @@ -0,0 +1,155 @@ +% [dr,c_obj] = AKpSetColormapAndColorbar(cm, cb, cr, dr, cl) +% +% sets properties of the colormap and colorbar +% +% AKpSetColormapAndColorbar('RdBu', 'South', 1, [0 5], 'dB') +% will create a colormap from Red to Blue with a colorbar at the bottom of +% the plot. The range will be restricted to 0-5, with a resolution of 1, +% i.e 5 colors will be used. The colorbar is labeld with 'dB'. +% +% See AKplotDemo.m for examples +% +% I N P U T +% cm : (a) string specifying the colormap (cf. doc colormap) +% (RdBu_flip) append '_flip' to invert the colormap +% (e.g. 'gray_flip' wil have white assigned to the +% smallest and black assigned to the largest value) +% (b) string specifying AKcolormaps (see AKcolormaps.m) +% append '_flip' to invert the colormap +% (c) colormap as [N x 3] matrix with values between zero +% and one, where each row represents a color and N the +% number of steps. cr will be obsolete in this case +% (default is RdBu_flip, from the color brewer set) +% cb : location of colorbar. 0 for not showing the colorbar. +% (default is 'EastOutside', see doc colorbar) +% cr : scalar value specifying the colormap resolution. The +% minimum value defined by 'dr' will be adjusted if cr does +% not perfectly fit the range. +% (By default, a colormap with 128 steps is created) +% dr : Dynamic range of the colormap specified by a two element +% vector [cLow cHigh]. +% (By default the current range is used) +% cl : string giving the label of the colorbar (default = false) +% +% O U T P U T +% dr : dynamic range of the colormap (might have been changed) +% c_obj : handle to the colorbar +% +% fabian.brinkmann@tu-berlin.de + + +% AKtools +% Copyright (C) 2016 Audio Communication Group, Technical University Berlin +% Licensed under the EUPL, Version 1.1 or as soon they will be approved by +% the European Commission - subsequent versions of the EUPL (the "License") +% You may not use this work except in compliance with the License. +% You may obtain a copy of the License at: +% http://joinup.ec.europa.eu/software/page/eupl +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" basis, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. +function [dr,c_obj] = AKpSetColormapAndColorbar(cm, cb, cr, dr, cl) + +if ~exist('cm', 'var') + cm = 'RdBu_flip'; +end +if ~exist('cb', 'var') + cb = 'EastOutside'; +end +if ~exist('cr', 'var') + cr = false; +end +if ~exist('dr', 'var') + dr = get(gca, 'clim'); +end +if ~exist('cl', 'var') + cl = false; +end + +% check if we flip the colormap +if ischar(cm) + if strfind(cm, '_flip') + cm = cm(1:end-5); + cm_inv = true; + else + cm_inv = false; + end +else + cm_inv = false; +end + +% set colormap +if strcmpi(cm(1:2), 'AK') || strcmpi(cm(1:4), 'LGBT') || ~isempty(AKcolormapsBrewer(cm)) + if islogical(cr) + + % fixed resolution of 128 steps + AKcm = AKcolormaps(cm); + c_obj = colormap(AKcm); + + else + % resolution according to cr + num_steps = dr(2):-(abs(cr)):dr(1); + % adjust zmin, to make sure that cr fits into the data range + if num_steps(end) > dr(1) + dr(1) = num_steps(end) - abs(cr); + end + clear num_steps + % still rounding is needed to ensure that an integer number is + % passed to colormap (numerical inaccuracies if zmin and zmax + % are not integer values) + AKcm = AKcolormaps(cm, round(abs(dr(2)-dr(1))/cr)); + c_obj = colormap(AKcm); + end + +elseif ischar(cm) + try + if islogical(cr) + % fixed resolution of 128 steps + eval(['c_obj = colormap(' cm '(128));']) + else + % resolution according to cr + num_steps = dr(2):-(abs(cr)):dr(1); + % adjust zmin, to make sure that cr fits into the data range + if num_steps(end) > dr(1) + dr(1) = num_steps(end) - abs(cr); + end + clear num_steps + % still rounding is needed to ensure that an integer number is + % passed to colormap (numerical inaccuracies if zmin and zmax + % are not integer values) + eval(['obj = colormap(' cm '(round(abs(dr(2)-dr(1))/cr)));']) + end + catch + error(['AKcolormap ' cm ' not defined']) + end +else + c_obj = colormap(cm); +end + +% flip the colormap +if cm_inv + c_obj = colormap(flipud(colormap)); +end + +% make a colorbar +if cb ~= 0 + c_obj = colorbar('location', cb); + + % label it + if iscell(cl) + ylabel(c_obj, cl) + elseif cl ~= 0 + ylabel(c_obj, cl) + end + +end + +% limit the colorrange +set(gca, 'CLim', [dr(1) dr(2)]); + +% clear output if not required +if ~nargout + clear dr c_obj +end diff --git a/thirdparty/cbrewer/cbrewer.m b/thirdparty/cbrewer/cbrewer.m new file mode 100644 index 0000000..26be891 --- /dev/null +++ b/thirdparty/cbrewer/cbrewer.m @@ -0,0 +1,128 @@ +function [colormap]=cbrewer(ctype, cname, ncol, interp_method) +% +% CBREWER - This function produces a colorbrewer table (rgb data) for a +% given type, name and number of colors of the colorbrewer tables. +% For more information on 'colorbrewer', please visit +% http://colorbrewer2.org/ +% +% The tables were generated from an MS-Excel file provided on the website +% http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer_updates.html +% +% +% [colormap]=cbrewer(ctype, cname, ncol, interp_method) +% +% INPUT: +% - ctype: type of color table 'seq' (sequential), 'div' (diverging), 'qual' (qualitative) +% - cname: name of colortable. It changes depending on ctype. +% - ncol: number of color in the table. It changes according to ctype and +% cname +% - interp_method: interpolation method (see interp1.m). Default is "cubic" ) +% +% A note on the number of colors: Based on the original data, there is +% only a certain number of colors available for each type and name of +% colortable. When 'ncol' is larger then the maximum number of colors +% originally given, an interpolation routine is called (interp1) to produce +% the "extended" colormaps. +% +% Example: To produce a colortable CT of ncol X 3 entries (RGB) of +% sequential type and named 'Blues' with 8 colors: +% CT=cbrewer('seq', 'Blues', 8); +% To use this colortable as colormap, simply call: +% colormap(CT) +% +% To see the various colormaps available according to their types and +% names, simply call: cbrewer() +% +% This product includes color specifications and designs developed by +% Cynthia Brewer (http://colorbrewer.org/). +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 06.12.2011 +% ------------------------------ +% 18.09.2015 Minor fixes, fixed a bug where the 'spectral' color table did not appear in the preview + + +% load colorbrewer data +load('colorbrewer.mat') +% initialise the colormap is there are any problems +colormap=[]; +if (~exist('interp_method', 'var')) + interp_method='cubic'; +end + +% If no arguments +if (~exist('ctype', 'var') | ~exist('cname', 'var') | ~exist('ncol', 'var')) + disp(' ') + disp('[colormap] = cbrewer(ctype, cname, ncol [, interp_method])') + disp(' ') + disp('INPUT:') + disp(' - ctype: type of color table *seq* (sequential), *div* (divergent), *qual* (qualitative)') + disp(' - cname: name of colortable. It changes depending on ctype.') + disp(' - ncol: number of color in the table. It changes according to ctype and cname') + disp(' - interp_method: interpolation method (see interp1.m). Default is "cubic" )') + + disp(' ') + disp('Sequential tables:') + z={'Blues','BuGn','BuPu','GnBu','Greens','Greys','Oranges','OrRd','PuBu','PuBuGn','PuRd',... + 'Purples','RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd', 'Spectral'}; + disp(z') + + disp('Divergent tables:') + z={'BrBG', 'PiYG', 'PRGn', 'PuOr', 'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn'}; + disp(z') + + disp(' ') + disp('Qualitative tables:') + %getfield(colorbrewer, 'qual') + z={'Accent', 'Dark2', 'Paired', 'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3'}; + disp(z') + + plot_brewer_cmap + return +end + +% Verify that the input is appropriate +ctype_names={'div', 'seq', 'qual'}; +if (~ismember(ctype,ctype_names)) + disp('ctype must be either: *div*, *seq* or *qual*') + colormap=[]; + return +end + +if (~isfield(colorbrewer.(ctype),cname)) + disp(['The name of the colortable of type *' ctype '* must be one of the following:']) + getfield(colorbrewer, ctype) + colormap=[]; + return +end + +if (ncol>length(colorbrewer.(ctype).(cname))) +% disp(' ') +% disp('----------------------------------------------------------------------') +% disp(['The maximum number of colors for table *' cname '* is ' num2str(length(colorbrewer.(ctype).(cname)))]) +% disp(['The new colormap will be extrapolated from these ' num2str(length(colorbrewer.(ctype).(cname))) ' values']) +% disp('----------------------------------------------------------------------') +% disp(' ') + cbrew_init=colorbrewer.(ctype).(cname){length(colorbrewer.(ctype).(cname))}; + colormap=interpolate_cbrewer(cbrew_init, interp_method, ncol); + colormap=colormap./255; + return +end + +if (isempty(colorbrewer.(ctype).(cname){ncol})) + + while(isempty(colorbrewer.(ctype).(cname){ncol})) + ncol=ncol+1; + end + disp(' ') + disp('----------------------------------------------------------------------') + disp(['The minimum number of colors for table *' cname '* is ' num2str(ncol)]) + disp('This minimum value shall be defined as ncol instead') + disp('----------------------------------------------------------------------') + disp(' ') +end + +colormap=(colorbrewer.(ctype).(cname){ncol})./255; + +end \ No newline at end of file diff --git a/thirdparty/cbrewer/cbrewer_licence.txt b/thirdparty/cbrewer/cbrewer_licence.txt new file mode 100644 index 0000000..627279a --- /dev/null +++ b/thirdparty/cbrewer/cbrewer_licence.txt @@ -0,0 +1,27 @@ +Copyright +I have adopted an 'Apache license' for ColorBrewer and its color schemes, on the initial advice of Frank Hardisty. Alan Isaac recommended updating to Version 2.0 for compatiblity with GPL licenses. Thanks Frank and Alan. + +Here is the license: + +Apache-Style Software License for ColorBrewer software and ColorBrewer Color Schemes + +Copyright (c) 2002 Cynthia Brewer, Mark Harrower, and The Pennsylvania State University. + +Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. +You may obtain a copy of the License at + +http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software distributed +under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +CONDITIONS OF ANY KIND, either express or implied. See the License for the +specific language governing permissions and limitations under the License. + +This text from my earlier Apache License Version 1.1 also remains in place for guidance on attribution and permissions: +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +1. Redistributions as source code must retain the above copyright notice, this list of conditions and the following disclaimer. +2. The end-user documentation included with the redistribution, if any, must include the following acknowledgment: +"This product includes color specifications and designs developed by Cynthia Brewer (http://colorbrewer.org/)." +Alternately, this acknowledgment may appear in the software itself, if and wherever such third-party acknowledgments normally appear. +4. The name "ColorBrewer" must not be used to endorse or promote products derived from this software without prior written permission. For written permission, please contact Cynthia Brewer at cbrewer@psu.edu. +5. Products derived from this software may not be called "ColorBrewer", nor may "ColorBrewer" appear in their name, without prior written permission of Cynthia Brewer. \ No newline at end of file diff --git a/thirdparty/cbrewer/cbrewer_preview.jpg b/thirdparty/cbrewer/cbrewer_preview.jpg new file mode 100644 index 0000000000000000000000000000000000000000..bd2830a54021deef189886b2ac901fd3ad644b7c GIT binary patch literal 92995 zcmeFY2UHYMwk})<2$B(yOrsBoPsj45H*9IU~?YR0I(a5F`i5AUV^LK?Ee{ zoHI>q==7I&=FXiv_pLYc-umyq*1y(mS7CLZQ&nf5v(Mh&{`M~1ByJI)dZ?nV0^s2R z03P@cz|8>n0phDyiLMe76A=+zyGBexMoUgcN=n9foraQ@l?lSe%EZFL&MC;lev6-j zh2^&79saw*Vq#(t9%*?g5jjCoF_Awv!Mk?t8W|}W135W^$W4};BLBx9+z)`77(bNY zJ3ihm;0iS!J~bY$1Au_zB*go-5AfFy-W7ZT!mC8Y*GNdg4%Jk^6+C?WD+Kt2gaia& z?*OnJAfP6^ep6WKDvkC_qFb)CB5z|di8=0-x6(ZwI^+~JcMG~ka)X|Mk%^0&=k^_5 zF>wh=DQTJe50q6@)zlyAJk!BjHtutL zLgJUCR!Tv!mYEZ5#1O)g5M1SPMyW$NV_|yc1H-)cWSJEbW z=}L1;zL-G)eRaEdK)A^wQ=Z<76Mf(89wN%l9v{zI;5fD9iG zTs(Yg01BL6-p+bM@;}Fa(;SqgK)TbZPvH>_(V=2mRIQFUfQB^RG30FM)q}`-pf;iZ zT!(cNfAf;B?erWuv&1S;-B;qSIVK=rBH!6u-5`9_&H4T4)9wmEK@X7}pth?(@2?%c zOy1nz+?nHLl$*?<&h(Sf`srO$yo)#Ftt@-9aN&NykBe>%2TKm`t1)kc}t-~i+)4$uiX zZaF`YIXFMrQ)J6jrhS6K-w&gcYML_IM?@)4VizTt{9QVqUD)q(>;bg=hfjVQ@W_{C-;eB2-I?rpZA$Vy1p zjg|si`>@nI(a6WpRi3k&Blsu|D46_$19Xotdl+8~-}j6w5qujsVA=s|<&1g_o`~U{ z5Yi1IX|U=WVenl#2ON-<&_Dv*&WiZwM=9m9O-CODuzN|R>Rw`UQ=%s1kVs49ksB4; zoIxZ-NEv;$7NV%j&&{F6s3tyZYVJj=gBW0S_wF{6G_!-7%KAn=vswnFOI8!K4aOv{ z0il%)Q@ClQlsCvWsW48{H+7rUYPEEt=jH^72!%vSS>4}FqdQg>2gLd=o1q9zVVmE| zAebcBx*ZOv)xwl0LKL;1n-$}LHc+`1ZAmzCR8=;f-sX;D8{1EgWzyI)vX72b5!? z_nzDA%PojHO))B*U3IUOmz_{jB&5~nK~Q#6A@SQaiZm>^c6wbB5cetw2;W)vNl;uJ zb9!_9c7YwgIIA895FpdBG!SBg9UNeRf=o9J&ib}k&sB&7?2>Ib`FWQ@4Z&d--~iHl z&|UZiZQz6>94ikSU6MO!Mh!G#FtSgG^0k#|f&b1!{}&rM4v5PHhx@%W7YAhYn`uq- zu`$CdN1%4beSN4p0Rf9KPu62C_9*&WFGRZE#wwz1@wfdh+?PiAvIB=Iwvv5g=lg}6 z_a10)yy*BKIQlB(J0b0H+B3GS=CTePa6^t~X^z4H>G_)V>7xJlX+u-aXpXcudp!1t zKDM;gH~#I5`eEXlHl4MBa^viCq`qCG;TL4JA`+})`au3wW%-Ipd|t0N7TBSMqbDag z07a;Z1G25k?qf1N2uGWUZhS9-r01IO(ObXpFCbRY#ShZ(RZaWp@){zC$%e0HZsCA4 zl4u;D(r2c=J^e~%ASxhtR_p`eUQ@!6nB`TH?PoEcR9phc0pEom&eSkRf0>9bL+?4e zr`452UtfP2>}RGP(e zafh6dz}iFe{3$ooSJZ7oiXal55tCY&agsqvhM0)hm-9qCLq6Rx!uB5aS5H;fT^|q_ z3b7HXEax)(aPwN@+as52uh#kt?BLg0xF`FUU{t+Dwj!KvbKT@^KN;0MES&}7Sx22b z2ksjB>^mSEACAaLBgXgSG#$C-KY02%<<70&Bhy*?Q1Gadcq)1h)g>a;i$ApfTj5G; zZj|!+$Rk3!Gwxj3zOXK+N2!_Ik;JNd%8TEGgVVBD_4RZd@X`PWoII=U@{wE`e23-m zmm6McnXni>Vdjb3y@6<5PnzsmTXSuC5oZ1&?Z+b7S|juyuKewP9o*sfGhs8`*w5;U zEU<3)JEm{J0@2}O4h6g4(uAr_X4||J`a5oA2gk25)BsrOlpAzu3xsQj@bsnrZ#+3lRx+#<$%kT(nUVfrqNUAF^=WyF_#!qq1voXQhMc3+RuR6OwZ7w)=x zA2INFd?jKFt&K=;Zy0Rv>ysDeX&U!rm0yl+W5@y3WsZi?FeQVmJgOsRf;CJH=>9v) zW56czXVZYL>U-%|aMQFK244_A@&zxqg2)r+Js*-b{V8Lw@Uz-dC%H9FQ|t za2dk$VUxUC7-l*M{mZ7qysjBD8coLyqdtM^(>j9uxL)ECM)SMOcLqe?Y`@QsOFI_I`p-SVtKFtnKv@_i%OtN__O`hp#p-Ka22Iga(RMg>i7iv~DmcP%j;@ z`E)O?Egi1c=3a~kpZD3ga+6#o`eJ3;S5?as#;C_4Y+GASxGrRq_=mb~JqB=)H3zMEwd^Sb0 zdilGy`0j0iV!xXl^2FCac&}d=KR9go9DC#CCq}Jh5s{qWAQy*?;AL85E#ii)zxpy{ zH=R_u#o$H4mSa%Xn#DA;8{>NygA@du2q2lnd@c6WLaKQ>oPKo8VbxNMdvW!u>wCQ3PQa6m zc;e=32&0Ug0r?o z3};QUO-$vOV|63%QirxzED+Eb0n3=U#y8Q-qM-W{MM2&@jdrT^=NLLGt^6m zw{(*E&05%e(vFGE9!K_bo_8Lp?$$eujtbq)O0t})-t8QM=P`%`z8AwTd5w;v z-zS<)LaTVrn)jNP$C(@XTC1fK#MI7d2k_PyJiC+_m#Y-S{MARAY;tDo8VKBH?P?Mo z1p1F>DsF!2@~?jSSVigU&CfKIk#B7Ahq1KwPc|2>@QZ7N#$;5~HrB;h0Oik;_1G=n zrci1>S6Z49XyNlo)P^J9F4vu)q-Z>wu3iYV$$q-h7JDW0N?!j*_wTXl+|ABQVM&-; z6m)t^!8`73M`OlW(rNGOTD-Rfhtz;;%=4Rv;|r|?%glQ`%Y@PKjjvq#T={=r(D3cq z&TI|T(X$6*E;Y^ob{mzkiuG=!c!({=9- zBpg7}r?Bx2VVDH(MmaEvlrEgNh0-&(OBm(N zc)ZVTcoXlw+E=c@-;)^;rttK5@Tc0K-vZ_vd|eMj(x7`XQEi{WkUm0gk*_t+-bkkP z_D;7+_9e6Z{$-7NEj{MbGymD+9aH$8I1Ui0#WaI^gb6sKyaKwW!wF1K3GC$Qgc69K z-$GZ@d>pW?dhp9T4jEt!a{HeilwTA$aX<&^?0xmHr_jP=$z0T@5Bu1)Gf_c(qBxF} z4Y|`r{gKn-jc#Ys1>1S9?})P~mI&Py$)|dd>|U8K*Jr6*?k(OX=}s*g`s(Pn&>3~dF_6SQ**eLQC&td%Ddxe~lpO6o z$mo);I_DWdGSH4mfvu+5L(b^pKun-)LKF;J-E@NO@)f`~MRf;5AhiC{i_wJkdPa?= zU*R>FtlKZbU6Ty8@8+>zX}l);ZZXmlVE`i&EE9s;k{b7`EbN|-%+2??QnaO+MA7CZ z+O58bz3WShcVZUh(YI%;5=Lt3r%5p+c*9ZU=#E_djVZdMAn4)sR*@+HaF2YS?F6G7NEiAFi%8^u=X_MvnoGF%<$;Z&CM67i>$^(pD$7YQd zrWOJRbub_INuxR1H3{Y)y+zmG7+F{8R$*i$iZOxi_Td0gEldM!&mQ=JNJOWS*J46F z7BE+tK=?7UDuAAeLMd~1bfY2Y6uJxPlQ+kzM)WB7>;w}kzUlW@`%$x z-jT}VhuNvc1}_2;U+4xD2?AL_1rA|eaV(8rDse{LS4)2qG~1xr!1rzqFVHcn zdP;5_IGL^gu?22sR?zV6?m;`*KuG3vY%zjvZtpun@F#VNHp(AIdX-5vC@_sENE@VH z@!|L_;Awm>&i9walM=|6nc1Vv;MA!uX(EO)GfV~!Sbllf)ppc};U2Z88n6!EP-0eP z*PoS|*0@cSIV*3A4_!=bkAk2e_YmFF&hRP_yFxvi<ua^I(atXg`WNF8f2!g_ztz$W!Ialr4F85GfH zaldJ2Hs?8_EGSA}*_ij&za?A&JbrZ4%R#<)?ro-1`oDB8^ye*~aNc)4V^-0)VA2sS z2%#E4GS`aa`&&Szc1E{BY^~@_>qBWcv~0-ayVs}q?q}}WJvQ+^7d)#c?jNtMdX^=~ z370LH`aC}VaigEjbMEobOHrX-YJqU`_H}}a)~PVMI--1{S~T_dmeCPo#Qb*-hkm!k zcaP^fZ>g|{n0E#_UEE`Naj#iLcs;!14#Ihv|2HzO13IN}6+L8!MXKXQE2M5dlVL5Eq%=y@ZP4FE4sNVvZ@A&_xR%>60y{Ce1>s)hjVLfGdCl& zZ{7I9-Hm7b($n$X;3hM`f_F0K%bND~(D>oFN5L-lh}jc|$)~@Jwx7;QrZ?V+*MCL+ z!>!%<(P^2kwX1n{i6yz!KW|tYdC7>SAM7ldMA?vmV%ucy z)TMqlavt?slzb#}Q(bN4-fU+yzq(wkuXv_10LE?CtBnJM+8#NhWfdB~nLNA+528{J zBQNfk#7BAsQRov9xt*3#)guPq7JJBRJVe5lm^Z5f?1gsEUT0V;7_X^m&CA}925$f0 zA=t1b+Pb`r{}X!^$%JZq&eYg4W*$3_1D<)fJ8gW)QuK0S_*pNw8ks_#u^7Zx+>G>I zGQHnPx@=s+_iJq2V4WUf`}ySx3Dq8XZB)@rio2;YEnf7p7OCyL;7IoO34ZK^151{b zmiS8h8%61Rlq%3vgCQ|Z+4f{PrpIy4XRjF~eO;YxX6iq%8N~W>9y#2L0M^iwhX-vf z)T5U0-zQ}irs~w8>3i3EKXYEcyQ``TiGPVU9QTNl)F`h(KtlIVQ?Gh^H`LV*yqc1D zH=uIYFx|?NjkWQ1H##^w?t?f>Bb)PN987USjvRp&3)aE z>urmgWHeW;Me`hbb_e;BlZT$upI@&Y8z!PCZv?Ic5n1UDntk=U#_&f8|0Vu3D*`AB zcxh8G)oTN@1L5J91v12Y>hDU7I<2h$macQ~&|>MdzoQZ!Zb8YB_>%QkPexoCoe3e^ z4?tGmR~gVR|6V!5zf$?VDrk6IZa0j$!?vC?<@TawEQ|@%owsQ5ZJzWkPr-c=q2St>^e*STWQC z4AA>g%>o^Ry-qlCGv3DFrp*Ttf=SjSVr$1df^yFDeoxvb zR|DT?+0C;IEe?Bpg;h?rHuimZtuXh?;HV?6mnY5udlTuo1fhFQ+Hdh7=ym3}rOo?= zHwxDGNMGE>_i~+~8~wBE!&_K=RNHF9j#Y=)gEmi+rq(&cX6BG3ETI6-YT zqqi4c*h$u(-p_t_uk7r-McWlM1{M4w$Up;V6rjdKVh7Zl$6GSh@M(UVnRVsnIInRs z>8Hf2_^pLl=EIt(B~}>?a*Dy-go=mOmdDQ@)dMt%oG02BdY!a*{m8^03iOD84tl2j zIM&g_n`*1qh0Ef3TvZTT0lO4F41BKalm^X^yS4cY*0(|V1MrX zVrOr-f#fzmq&@J)P!mU<`pe8IOzgP&qhDzby@~8s+*Z>fpuF8JNW;Zc4U;&-#^3A- zm+B9Og<};zt`G?@IEHgtuKG_&4`_PPJ|Th>qxlaX!}FW5*%^xGTCoRYvG(H&PftpNK+S$>;rz{ReEf#% zg;z~E=iP71URiHwEstD)9ld`s)yI)W>{PN=TNfdn5l^Dd62WMy)=I!e%n_+iA8-i! zArrcqsN>^C5QB&zsJIslP*5}6-kNluDDehm2%kaXFVnYc^g@0l9y%~RIywlDIku-# zdvUE-`^Qv*%F}K}6Ot23eESx5gnsB+Tql+_Sj*1)!3&r6!Fy_szw5Z@85pc(>!pCm zGw^mi-_Rn5&*O?~O01SpY6ql=dX59-HH%XZ%yo&`U*t1v-sRPRh6mD(pr6chRS!2c zo94z^xSqX#{9LhhV0gxNcCK#lkFGHUavcV~T_N>hpQxHVxK1QQE!fke*4loDh$u30 zanX-t5K7i`V;-VO-UPLMMLjk2_(0h39sh&+`D@b|sk5Q@8J(V=pAbQFeR$d%c*4s( z$Yr&TJo*kq_o#e_B`}-RE#{qJ&+uAleJsEgI`p1~e;50cvnaAjMQw4;)yf8NQZV$b20d|g~;A_T(LjfV3|;tqnb1Bb9SY5SK_6O zm8KPGhzTN=qf%ukxSc@(eSkDxT%d5EvY&n&b@qBw($@`B2%6}JX^cH`*QM#W;fEuw z*xO^wI^$vn8%>tP!^Oi7OentA+$C6vOmhJ?3IB+~Q{oTKY|CLw<`U}OVa5YnNky+v zUEspX8z{`cPr-W|vB^GM>{w-6n@SJd?CUaKx8F@90fOJPh_c3A{aG zrE^GfTggM_zLNe^6@Srv>#-4T^H@acD zb$Q)5DzoY$+DD-)ms;%?B-?rVsngx!k8Ezi?`(VM$4p2U!g^H9z?R)i)ftV1v%9{D z=P>ZW@F6o9UFL<;cFgn+b`X31M~E59sI!sPwc1~o))5z8jSC_0BXLvM#0MO(X=gav z>Za89)tRI2%;Y+SnB2ZT+` zLIQc(3w$Vb#?eYFUPfV&N8Gu#Z&ZH%NPVllz+Ov?Pn}3c^bk0C$k&h7Uv{S5f)NR! z)$asrj?MHY6Pt&`Ik3cb5xDuj+5f~x)omXc3#@Q=o_BnJlKKk|iDR6ctG4y0{Y8|W zzUVh3W$V2H?iR(xPW}L0A_4yI+hsy66AJvc4J^s)+24>)jg1XL`^QqQKdOA!MHayj z`!;AJ<@-w1EdHeJyu=~4e#ryD+LB8k<9x+X(9)nOyoZ^yx2;Bv)whTJlM*{vvXGX% zAI-H^lHynLVa~LJrzX^is%hwxnup0+B9QV?N_@pnvIG9!=sc(zN@rl;^RoD4!!!AY z`uJ07-Ko#c(~5e7^F+aK9^(P$5BVrtv0UL@Y&72|GHPSz1c%;-_HZ#8WRrc&c1vLs z5&p4d0_lu?*Zt`56+M<5K^K{*k*`X2V8)YkCIL#HL9cu#tuJV9dhjet)?7HZUyKV0$Lg<%J=5{3FF5-G-Kc_{|IXgK3mde>0SB2Tr*zS0AmI_u9H8qh zIej)p!L&}NaD=CRhJs$8g0WaiMCpBS4A0^FV9ARvFPp0vme!f2MZoOK{^Lx`tD0+d zkke6Elq5?3F%C#anW50|6%+o2ZIH9^zyZ8JGi_*KR?H(mUveQt!8BpD#{sw*t8oCl z4=dI^UX?!p(W%gxxatWRZyG2+bw*yAof934yY^ep)oTxFaO6oF_NR4XnzKRv#Hbqw zT%keRV&_3-Yg51iI(ZnAT=eeQAW|I-#Q{qFqF`rmw=zMOcn)y@ zhDdA>O&X1MSa0BimG+gFU6yNL1EKX0b;*b(nJ;r?b`}b?F*qQK2hG*Oi~};m)xkXr z&S_W6OE?ZlfJL7adJd87Yd631msFQ^hyS=N!d_nhOA~&#=iq?FSFjb3;tA<)1DSpi zaAuLL^f-XdjBni0fpK=dEJi2KL329*EZevi9SRQK5|rRjMhMIMtqLI+oNvlil2H3? zc2>blnS5Bb=IUcWTJ8#r1c_}l3rRoHse+*#LHqW+bc4&JkSgy^fjo`P$ceVHXvaS- zR^ufdxBOmb#bxlIgUUszwN!NPJpK3_ia|qHH2J>hl@^F4%NFv*4wnh$jAp9VCw<=U zvEgJR_h|fz0%-<&WT=Ws3H`mmak-vyW%C{Pe!5sOw+6^vvLfoqamtwk=*B(_g>>+j zJ^Q_3#3WC#9C6%;d5Do~uC3)uIf)O?Yu3SzjlwptJ6Iebi5wSqD0Fvmuo!M%t1Q!e z_OyENNmf&QCywhuozhD8aTm4Yh`JvTeQfOj?Ce(GsgP@0$cA#tooM#vZzKlnJS z3EQXL;-l`+7v$|75AQ0*v>amt!ORLgtf>Ks9LnFSe$6_1*R4qNnU;HPR$X;VV>dbg z%$32&Ax@2;h8BQYzf=sbzKJegDiBn!JCl9Lmmaa~PC_FV5)bLH1KHTG*w^6V2QZy- z9B>ZghC4oA-=iC3PC9}G@L*1CQIBxI@w66d^Y@Ffq~Al&H8JJjJ%6hlfD+asN6%zx zvMxg_YBD{3T0K!pddI;rU+cl7U;q2CpKtb*myc+V1~zC1J;;8nzRnMQ4c^fq=-bVE zf_jrHkr&A-Yjn{Y8aP1S40|bt%>DH1P0#8LNQcAVN7rlWg8Buco8qHmSrL4bcmL_i z|F2@3yaM}arzoddrf->c7Wh)ivD;SfHD07#=^$IAiMG(^>%fZG8O6gRJ+$ZsGkcdQ zDYnjbc}B)nW>|t>(|ZDMzD%mO67@|_5?JC-ip?2ojyq1TXp9O7W$51?Hu_bO{|Z0o zSxgF*>pt%q=$I+Q%?KUFy_3sAUa^m?27H<~%HbpZHoQiv67DSfC5&S#x}jOvD4opqc%`w3SSEIwSJN3x58}jo2W|rxjs?XYh|P+{fB)4 zj##t_=^g%TZ)Jy%IYhiJa`1U?{QcGByiEUR;~WT;&~M1FgSJG8%Pw=sWTKE*LflI) z{yH18q^~wajJe28hQ|8QU_ZVR^o*ruMNH5_W5{zGd+YKZA$_ zn_b77(mT>`G#M^M3+^8xC4BlPBOJMhPuH4jSrXKVW$z0PhrZ$PlDbu;AcO;Y?;}B@ z`BDcgPjcJ@kJEF^ON)J%F=7S2ZfZlJdR@Q~iyN^q@wdGPCZ}M4} z7g;ylK9$?XT&Z+;SZQB)B5g_c0gCol$6jDPs|~=Rc6X#7q;D9#thEv3uv9hh3rEA$ zk+7S_@nE$`_epyA5#81jmZQG|O=!O%u604&2R53rx+mwML#)uyPh~cEXw@yOD(RiE zzksLaXU?z*%5#Tt8*4MPd9^qE8f?`b%v+R!HvRW-9q&OPSP|4S23qsy;H+lmx-Y{* zIp*XOB5s$FDgApD=yS?t2$`KDO3lvb)py?prxL2s(%rA`EK@mB_-Gl3KF@mFa4kra zDZlQtnDsGMX+$XJy|<}cD4KroGH9jWwVe2c%5SgqmN`{eF5rd(m|?4*cHn0RI3SmL zA|2T;$S=ornAovC;Cs@~#AoRIoQk1~WF~$wX2Jr{V*J5rX)`O{O!FAQ7OfChWQ@t` z3TUx3o|1HY(H_tKQBjU5mYlb1M&*O5hV4K)YhAUmuY#vA;mr5_fu-7y?rA4@X1WN9 zuQ@*UD;+emU6=k6pk2KFVEaK)Wd1!hGGlfjxmbCJA7cHys7IItwBq^ru%eK&c|+|& zt50OYHiT7?>#^wQi;u)SV~a@P!^`RAJ0tDoiBawOvFkQZl?HhS?o*l*%kLDApJ>W< z_{m*9E^$8lz<72M<7_!s9Zc~oNE$e_#g9YUEV0r*J`4*kEcxuslymlERvh#5JiU5s(m5qHmd4<}G_SKsYWqUgF` z@TqBNR7+Bxv4LUt+VphCr?#)}ry~IExZWT5B{oRswrR-_k1c*NS!jh(@13N$eqWJQ z;IT>Ekm{Nun<|1&ySxsQUe)UlapYtoJEvCV)F8U7d+Q%Q;Nrjb#26|Hqn_)h{HVZk z7#6c=O0sgVe8h3%b*Qngz{CBtg`b?|?CaVwjPq`RF9(0Q@A$#m;Ursrtl0MsxfYM$ zt##G0%ZAkNuGYA8_F6|-Hf)Dzb=PYP9j1+!AZTtJ1fwEuoW7k=gl z6KY!65gRc64=l)6@>Q`qwOuUk%2lbCSUL!rLl1U54nu9iLamRA@4psm+6t42hzHR)m97gEA!kpZJ8ahE=aJa6Zw?d+77 znS3?am7lg{+F<3XGYBOH`_D0GEi895QtuLqMlZk?I>2sD!Jdt!qYMY}^f72E;pN~X=tfQ-hAr$+kT#;CiT{QN4@Uf%4tFQ# zUKf)2X+Xdtb8T2#+~ERf7VQ_#&!AXIQ?bpT!sXapW<9K-8#rA*3fX(;+OaBot``}2_F0)rt zX60R=?&{oV-hYc$_H#iGF{O(9_O7emul&cKv&y`Z=CzgSR0sb{c^X|S+60k$9wbNW zeA2TbnJ7tWM_bktQc>wV6Df>Gfw#g*^XaY@(G;>Zt^HXOV&LmCtAvMDT?DDr810_4 zt#$N4%&pK>K*>;`h(AC@2_^Q7_?#zCan0r$VG0wlce~WAJy7sNR|`3q2s~Cg zJGkAwW*N*mlk_%!cFb@cl(`>4^g{u{ghFCRwAGwA^eXCymx2cgFk^hN z-St3ll`w(%8~NCXyqtN>{LOL6+8eQkgRagj6~^kRzbkKZQ7B3itYmEc*8K^T!EgZk zPTA!$Xx}GDO~-HmI=B?R14hdnCITi}vdoBEXP|=IG=MFdU0y3Kt=EKL?|la&1^5Y# zsj07YI6^<~B2|iV32t^I0!BGVgx2LYG^W(6K>&u8)4?`GuR5fbX$WG8er+YJv8rMX zvC!fCvAw^ZOi+j-k1hgLN7vx{FMf2nmP0CfA*9yW4VGPJ_@Q85laG1h!KF6V@#pNC!pbGrQbP z|D~S^4NXhI-thKVMJ<1QtkuG|;RN4h?ganyq$tjxG?YxRPh1W(CydZ)rYe8xQ%d4q* z5ladb_AO!~{-pUr;c}98#YzIZZ&}N3+^GGc%%)NFLWW8QPYudv85waYSdNxla%66k z@kn4|O`j65{QRqTy(!5J;1N{c%s$GXz0Gx!7L{Ru`*X(B?2yM)lx3h3g1NG?r1x zSxzFsZ&e@Sctf50MCv+W_%A8uGbzHBEh~B2RX;0QymhQ#H9|kCh1NlW25Ezm)K4~u zB^9del>G96^Rbyv+sCL@>_J^=P`m0F2hz|n>j4fhADeifp$wB`d$bI^Ie4@}KIcAi z>rZlNU0s33r)}eh%hnU5tY#&=st)yFBKEn;xqj#SC;3VTg`N(e@0Wwfqbaer%V2%U)VKVc%G`$`BK{XCcH4a(PIs;|D`6Q}}Bnr#! zCn)^~{50168z}w8?B%vB@xb$M@U#AJnDluXWLa)j>xW*c`!gKyDH0BX)1m($^v+x$ z9f9eqE)}WTU=f=MbgpPj0+3XiY)Po ze7~U2IqR2gUZ{f5#dbQ#8V!EH`sHCx?z~p>rqEMYzBi$_=`ul$zrP%}R|7vjro^CW zn~R}I0V}JC=6<~R0k(MD+W-;l?lahQG7Pj}#;}Dokdcuufv%hRj7hgd)+#^4j*fuN z?{Wr=o`sMj#s(u)OVxbh3Sl526z!mOkvPxx6ZW&ON*A)%%!3W`GDA|Ouc)5Q$AJ`& ztPq-53z0oNg0(u2o_dDWKc~+x9z?2DK`+U+7Lc@HM>MAq1asBn08OHW$lH54M^sUO zF54H|EUrOFy!kR$vQ?Cp)EEgfaz=A0uz_{*XNe%G6a&&gY{T-MyBkxJ>=mNJbyd8QXYzIX$gMffe$H*u8KD=JTLEh$ZU7i14 z1apqpE3y>8&4*xQ-LyTyG)m@nR?^LXj~V{@|5ZrD#^Lg2S6`>dwCcBQ$}@|znki<9 zaTih#oe)k56=DYIPM_Q_!)Qu`!~JTccPY%ya7X0bDXRhx>XEB5yz~^ioyfsiLBRE1 zpB;rc4Nc3l;Y!8i(0f%__=|RNs{d2{o_DiNnebX=6| z9?60RCJn8fV~b@6Cn1*Ep$aB5Em}dvc8!&u>_=f_J+F(Tgo{#{?dXXB#P|qaT~!3q zY&;5P(|g0FNnNgWEH?Og8AUfFZunWHczNCUWMtHx`l*#6ik>qrl!Yu-YekVI{ViHO ztE^p9!gEZ3OH)X-_Vp53mR>C>u|TsvCnUbE%SU{jO|%l}yVx%vr@E-=G(6Fv7hN$d z%eU5eoxZK*1SSfouD?KM^VI(+a~}~CfIK*cRd^IJQ_8a0ObHANAO}RN$X0mD3JP6( zgxYFqsZETL4k}?jDROSAFRxeJl_sF7Ct!0|q!;5W>N}vcOH7V}PaWh{MHs*HqCHx* zB>ZbtBOUlYAl`EvfK{3>JM9Eh6s0#Bi*UdN5y-sLV^!*m>e9Rj)M0W<0FRVxP&gD$6oL@S5Mo(E{d&^~z>Pe{*%?DGUuG-c9= zz)zWIbR`;=t6enx*OO8f8^qcXI=Y|_HfM2+o?Oe#F){+}%RT}G78fi#m`j*YFVhXk zHZ$g6lXpxPL@`em4q!@pMXSe<^McDXEIwwMEqG#ocvxsIw^Vck0qIIFLihW?Mpdz! z*k_lPe*+Z(p3t{`Y)uNZDt^fibw-Du_zusd!Sq1>_u?8j*a=Ild68}^6njq(emVvH zQg#nUDld)Z-{0yz@03VM_tzg-h6Y1bFW9(*J}Tr?Az>5U*h8$x&%c3+i$b%G=t6n# z@^WQb{h2hM|Ce7l!!trXE+1Vq@FeAXT^C;)Qup5S$s~sq`*9FV2J|flK1AicQf*n)F2yF_YSqqo4y2-z&yZ)PnJ4}!*>WH>864pLKq z9y}0DY`^5aYcfZ74z$fS;l>j}sDs)BARA~r%E zY#z|^QeC}|1LU%CK#MZe_Jp=dXEpkz>H8`_66ZL2^syFdRcA9hRz5_u_XE~94|8-U zY+y82vx+kfMX0GvLT$`m-J{S6>YKuE9I#OS$Wf$qveIBiSpL^TSgui4vikKE`UYN@ zmEFSe64)zHAI;12yUN>qFxQ}xEoUYZ`m^eGp=L>FFI$&4n&VOkbA=yz$Ovsqo@e>} zWv$9oD*g4Jxf#hymurkt;d>SdN2*joF0a_Zk=G!|XJ^SONqIl5Dn?2Ka`>v6?($xj zSz-5B4Ak|CV9B}s4Ei*+zlAR#E7y_EfDCTEzU7x&5<(A_sOrgUhLDXIN43Y_QdD@M z#%OXCpJ)9_^Oefeex|_=8P@%nz^{8w%eLvj_4>`ePl_*Gl)4G<6&d_CBG$2^LseQ7 z%o-*4b?>_gnu)YJ+ShU05G%cynocRKgfs~UhZiIB46BgGO9$jizLjc?BlC89=p5Nu z{#DivCy1*b+kz1o#@QG%`=Uh)9Nbmcw;jFDgsRm#Q%;uZyhtUd-k)5_JU$P;8(_qEJ(es zPnr4$lz|20HAZF~QdJJGH8m+GPg%uhm02Pbtv;){5&eQo0E~=1v-*&wUr|d3k%5v! zo(7iHjbr=dr_o=Y@guHry($>|^{_XEY{g=OP%V39HL7^ZR)K2dX6%7cb2GiLamTYy z5g+TF+%gh5X}YL;IUE7Mj^vuJ6a^HCf|=l(#l0~WUtY3Eni9KyjrzNq)XUaZf;jM$zW&!=i0S_wtLQlXEG)rtp zusT@0`M8Lz3+irhzbaON9@C2Mw;Q(x%>{%MyNI=Zx9{2RmP)IwXYW3$t?WUTSpe%3 z>b}%vYX0|=-P}WtgtCsd@*VT?ox*6HXzJg56-XT%z$K&x7GW)1a_vNw_AP0lKR0t= zXJDI&rG+r2Qz?+`p?Ds73FbA+o*@29LG3@hi9sx+o7T4U?1lvm0QB-)tsDIQ95x8T zt@B-w*YSR)^SFyHpSMI}K;vO$L-fdA#cYg4ls#P+|dL2iO1`@G;hn5BY?4>tLc% z_Nr-4rUQX|z8$y8M6S@nuH|9UoG(^*Fh>qv5vLscIUjT1*%@`qIp+0Zex;xHX77%G zY#1Z-AkLJWmO;|XSMp2&^h1g}aMa<#lFsCoOoaW7U(*4@`SI%9HvqC|doWq$ZCdTX zN_l;{IucUZ4>}b~(5YxjvwyjZ83k`y@{=#<8|!01$=Q`iX(MZBoDlj5hFTq2j!+Vz zVSd8^yH;YvH~3MDxL3g@b9ATnHyHJwk^Qr!|Npw33`UMLRyq8b$`uxFbEaX)B`V@u z#%=QET3@W+ZFGJeSTX6n{6XvbpbcEV>Rmfwck+rk)sLu*J!n7!_|n9ucJwACjxI*A^6W;pT1U$yldiHDshvHc$#A5d?x?kB)cOG ztg1B4Q+-)#h$$VHe(*ERp*tLS;}*ge36lW_A}x;()`X*l3I+4BH_Kan^0lq(az1Cs zdmWhLbu+Ro#{5nmtpx?LTmAsQ-fB~+-}vsne0TDOACN(wvj=)GiUlL^Hh*%96jfi9 z1!z-rW(lV6gmGEKYN{)1kFvdwaThvwm9H6CI^)?yg7a>W$g>{wbUR{d`4J?KKXhKSu(16CP%^ASa)mu zR^gCZFO1SiMkI-+l>ya`<}E|tVuGR1>OqYIUCVp0PS9)$>Ab3}O8;C|Jq1L-|EW2j z`%o}RHq7n-@2`~b$dKw;AoB&a2J1UB$CNxY-yRP(PzZKV1+u^3s6qokV=~xzfH_M? zYkGRjV_(4)tP=UZN)i9B(#HRtzww13 z3zo@s+5;JCnueWP23^jJXrKvY@wZna?{XAI96^vq^WtD+K7oDM4&uFM^_+eNCZ#pp z=~cmc$M6e?23YfC2d6<>xVrQTme=2-iLSgS8-2#7kG;VY0pGj>z7~prA;d!vcY`Y* zZ5dwW!sKHP(9vF5O%L0<=Xui3q*h1IZw{!&FzmvnO5jM1(}fF=i|7F_>=<|J9x;WS z9&;PA-4YzoM#FkQoA}ESZConiykHK!N6nyE)F|S588v3j2 zo30U)1RQcC%fuHGZZyi~Ggc}>3Ir;jUh37(ltMT1urhb}2jVutGO}aojo<=uW|dB` z6lW*;Vhaa+j2M}pug-BjS*w6e{Oa;Z@#f0P0)I(7L#6+vohv$zbp+dDt@Vd>Km`n< zPe3k9*-Ep-o>VMC@Gd@;0{En7Fxz3!VNxF7($SvX$uE^fx&Lu!(P&p#TSb|ic5$c2 zX5ts^NXk3%3h~)gZ&rKA?Zs7iC9?#r(}gHTW1Phgx*|91luuyc7wJmcX40^>--txaRp^ls1beH|ER0QL0)~`?9rqbV+UX#o^TAoopS%iO^EX_2YJg>g#-;x8nj5^Qr_K)QP_Lu6rpD zVCq{LIBnSRtFoE~mx22ka42zxh7x_P_B}xlm#x@qI3WH>+*JB8t|lc9&mE z3~|ngVvMTJIj8EoGwo8iPZFHs$#Bl9`jWfkTYKUC+`fDx?o1Kqh@o%Vd}VTlHk|$? zJ*?5WfJs*6Nx!%1IC^0>RS%?_db^X(G(r-s-6QvS7+ZHX9eM{lU~FdF3gXz_n#0BN z!;VVdbO&{GFUNgu;p*?k)NG9^YG^<{D`FK5SyA&G3{POYGtM@msjeNedra{D|6%Vv zqnhg0w$UIe3L=UiohYCnRRyFI6#)Se=}m}A?;Ytu5KyE@5tJ6Ch;)%&qx2#mz4zWj z4%Q*mz9xh5wRtD!Q?=eE*!v=?SFh)1Xnq^8 zpf>jCR(dbW`<`e~q-D30qtih%mG{WJY>soUe=GHR!|{7B$TA*laPy#dC_x>v&B%B( zTgg>hsJfesup15?(0G5O+Fykew|SD1iw%5F zj9WnLFg=FcquR_?Tui&AY*3JtYaxmVV67k)fJ6OI>;Tgryo<@8^T*_lKCN`%B-u{h zrR*FfNLu?4E!&-9HePiSSOw`Pw^Rj0O-Tp&?ZR=?zvxCr0phV01ojf5q4Iq2FkT$- zU10U);R}Z$>ToJ>_@{u{|5^%|+;E(tBm%IQqhOx>5~}Tg49;5_KwTT&S#*_-wxh+l zT&qBK^#fU2`0}B)tOBqQR{;z0i&#WnwS;eWY(Q@m<618gzApq`Y!S)dF(X_WWBiMF z4VP&~xST`_eiJVTPK(3fGWYd9XNkW9U50LX)_AxPs>WgxdB9L;N1%H36BeAYLQt6` z$}RF)yR==Oo4dCXNKraEOv(vkP9i}Bc$pY9ep}z?juRo40GAAEt#dVU*KF~Wnl=Lm znZaDS=j;zh!~d9{%J6C0xzuWI?B!UKY)d_7#y88uBSitNx4&m&`Q0X7(@$TR7DnEV z*TZ5OPY|wHP4OuicG+fE;){ANiBRjF$(GYcT;x6X`KaQ1Ki49d=Oxmlter^U4;G+Z8m9?N{M2E?n&SFJZ|rLa zQ)9VOVAw6coJKNRIS1|b;bisYyOFwP2C@sYHi}EFWlTMx7 z-09^CPEq9RhT}>wS7~+&?Kt_b%6gQFkS9@5h0UsY_H7fRX%gGi)te{jCAmL6fH9kK z_1Tj0}{B^$+7UpX&Jg@V`uv=XC0;7Kxb9Cki&l7mSUe||M zbzH#IVjhR9yI~i8rIT6*9m9L|CJ8BQx|`IOWFcA`J|%Ge7IwU5cdz8Hi9_20} z&dyCqQz!mY>hR>B8rXiyc5~YYdGa{$A=yL6iT>lj^he&%!7v=F56LrO@MIV$B&=9J zu~cRH&K?kO7fFZwqplj}cCox9W^9UU&QIeCha;DX3;WSCZWqF?8N#vJH$fv}|EnpN z$i;1n@Xw21Mv8a+*m?1OyT|e)Z^GU4LmQP#Vezo#@AhTIq8j&?vHKE znZ=IBSp&B`Ee9no-&^H4Tc)CEeJ+rmU6K(t1M%Ye)v)jY8y97M3WfBT9lk-^h$Sl{dwLBrkUTX6B4|7;ooVK zQvYM>2ipXsmyYftx-c);dWG*oVNtlE@~@f&+rS^~kMoiulM?_`F{H zJ=bz3Z8nRy!gaRfO~DVCnNOJ|S&gg`_)_%TlGp1#>#p>!O`Yj*m@6Ce3d!AzrDx)0? zUq5VUYvv>Ue8@W6zN3by?980`8u@1V3a&UUh+HL=a&NGASVLih7m!#G5o&~q@5pvV z-pmlA$ZSXE#g;1&H&)h0x6CB_=yA#uE}duvgX7ms_b42fIAfjeW@;*U55FLZr%0fz zo4)KTp$Jlkvj7^Tva>~*XLEN)C;sSW&@}xC)9Za|{%L5Gp_F+y%~8X%r7(wCJxT(@?N^rY`IXAT#Weg9wML*lgocT8XveWMjQm46 zE?1!xetHPVgJu%WT&tzRe!@)B=Kq8xmB9D9KtAQ#z4s(abm=j-qH;L5ka+Z}(B0=>X~YYh9rT`$=sbJ_d7~^q z)mgwP4u>Yl!m-)T$o&HhV0QpesM8uZOcq$FlsDN+IxFjmBydv!)X*^ps2Y?JI6fEn z4!1z5e4o$*=&QpN64}3y6x6enMD+>PiU33%_rVc6e^wT>?egdYqb?nwK?!~Z(S!$L z{-QZbd;4HZJmW{glPF2*&smG%I5lYi8_=uz5Yz~u$x=Fj3r-HeiVOfvza&X5x#A%? zPl(as30b8$DM0n913Xa-O3rMGJ4oCEX#g?a2ZGfcl?;|gm#j(5b5=OE8e^3`FlOoa zaM*63pp`+UOfkXNciU*JZ0U(b3*W_tgs6MlAV;n2e*ceU>%aTo_)F@U$}W~q!Mb6B zzGA}QdY_aI2YM~^{z))Ib1gBOWYqw)xg)mE85-QGeO1>U0?HwD5VbcPLLoNKfs003 zcOWK8_xlZ%O9^NE^(vyod3=8#1TH25 zw`%6#b-&9Ko6_nKZLL-+a!H-o`A{ch{jBokNCN8nLXE6el-h#x_)rXLYxpOO1u-eV zn}#Orf_JOnw8dl{QQ10qv*9B@Vc*vd5~Wub5XHQRc?={v4qjLvMA-Xv8lITs>Q~E7 zDb z`MtMIM@^W`p2ee&+E$MO@4XVVS<-k%4>j`g(s#(c{nZ1c;fW5cbb4@)wJGOv4rmVAYc5=h~61#nFPGgOtu;}{Le6m${Wkdf(i z4NLW+kW0EZ%gv*!Z;eO2>78Whps#A{4WrvIc6NIP4)U@rf(3ST+^^^^>*l)r#<@Wi zssY23&x`jYbZ@U~e!W=E>?hi!{#wcrxD^yM6f9{Y3+w`Hg)MC7Dr42hx(o`0!qQid zo#bmi=K9+ES{P3rF8}C_YH^&dogIz+xAW;=M5iBQ245=qZFXVU{JZzmJ-v}?pE?`( z2(qCM_!Y5?d+Ue7bg=&XAxswz*^{mXgTloy`$s2OWL80~;~n&*c=_C4(P5s%*^Q*{ zuKX6rP3l)0kX*Uq+;-mnNm(3#PzIMAdPW*@XZkQ3`&2|{lDWn&C>}{pqi$i{)9aO)2TC}KS4W%a2`ku8rGq2V3kN=&xx70J%{U*>I>$|QAzpF zBp}G%U-_y6Ezspz79AhwvdE<)oeEW`N`$x!8%g96Sy@^K>x6i4ct_AL^$*t;7|TiE zBYTsJO^=-9`lPouc$#tQg7=7;cX2o5^8thmzcBMX+D<9jhvRxMg>qk#o_EIt|JooX zFXS0?eN8_(6J-hXzdvDCAacYKEW&dv?%E$D#-U}z!+>yr_%G=Xcu8T9)9v6^fQ@2rffzIv4+InnkPxay={<;Fc3?gD3YcC-{-6%l zxis3yF)r+Zk!!_Aeh<^NYeWmgBoB6;c#2UQIxY_h8Bh$#l;z19qiVD#I}YfD@(8EG z2^MJVTbvOTD8>F0Mhc?aR7QA`R9hww9dOP)h5FH$ikfgn5;T=9rj#ijlq>1;St*y5 zdaVGV&(lAIJ`fF=kY`-9s_d~pE7yPcL!RxByivS4@=%HS>Yad@NICLP?~Nl>U?lBU z2xZU~mOUa|9hi+IrVOL}xcYm(O?g zY8pjP!C#2yy@W{Gruf+0w(&Syc7N5Y3T%>h#QkSJe``JeV&?+S=VmrvsW15}Ub>T2 z_ctXZsLohuKUFg`zG~RWUQEJCe)>L)6Pt}VQk%Jqq2znCp<7t4gy)dGs&IMHCW#ha zWh%CHU*Ca@dE+6nzJtn5+1V3Cxsb$~&&u6nB9mU5K9uUnMY&b+rh$Ek*=xis!7DAA zbcidPNZz;ESUY5rl;W`aYJra}c*ob+sVoDuys5w5LQ+C3Ie*QAJAoyq6Wh z_Gt2s*Q#pwnV^5^ba8}!;u?_C9Mkb%gXd^xo$;>S$t_2pFu^qLUk91lXX)o9R0nD4 z-4}sTP2mB~3yH5IBImq5a^kvoP;dVTwYn}DbO+Hur>h+sic=koAvnT;+$F(d){-l^ zdSHe@;2;+KvCqar2?9ZIR5;!MhX_b04gH~+BLdxyK?7~g_C-z0J}?I^Dyaa~?!lSd z+JtTH{{!^>Rf6WEKAdjR!Q$X~$yr(EnnyN+nkC0f)~M)>9tIHrN<9U_^lr{?)H?05 z=&WWT_%Y^s71w5W)ZMKtfKZu>UY8p#rACCu{lcgKJd>ku0}xuJ1VNmja$+z_dpd+; ztut>oMgZ~=d{KfkNY=mqC9Uthu2V?U+r{uT&2`+(!n}PPfKS!^hc%*)!9jtkdoPlV zf<4<3>U@rMO-DRbaes3mq#QsL-`;hc99T@|!OtFmRpf)oQ8lV4<93~v&9Qq;DyIii zvR-nwyq42J1?WR0m2L1siJS znLQ)P{B=!K!9{5MS8h@vOGT-$PHmVbu+oI8;sjN(OiU2(=-~IL*|7Wcr~s=@rKLht z^!+PqfPE7-1GxkEsSd(v>Q>UFEA!m#HaOzVg2KfuEPVXj6VrLI>v^I19BT4-H)6hO zoVW`*Ak<26CdniByQ+F|rVS2kMjsJn28pzm(dV0j0#f*x@53Ql(;$En^}-K(AK0=U zda+`O-=B?bsUJP%zzcE)4@z3avj!)1#3@!%+EJ5AZcY%&1^K&TrEf@&^GwC5s6|3G zQ@aYr&ui|i1-M40e+lzZH87!B2@Ye7Ht2ByaphgB| z_KJQd8CE0Es=&~S^(huCos_#@4%XgH;XRc(^*YDDCtbv+p0R-=zOURpO0B;(`P$D( z>*9u1Q3gS5gDg0PejCTwGBthoWwNvfo7}K*tgA-Mg{u#)(lJ-NO+r@}&hu>ueA|2z zp>9~J#G=SG+^YW~(%hcO_9Wse11Ihr{jlMO>AM$;a28Wco^C0DCkj2I!$f76&}k)j%c5fXC2v^Nx3c$xBu;O)}f{9*K>*a5O(e zN&Fg&wA;-n0hZY)sa=6HsyjfwDin3`W<%W=yU(B>^s-UIVK7jB?gYn~yGK2Ia=cAa zlwGEJx}8gkT@?7)E7~gGXZm32$RE%>snAadOTQ3LHDz#I;+4|VzU|(mowKV(R;iTY zWi>Pmt;Z#g2qyJT9i)f42(bFJLzjDn@A6$t;u*;4wNZVe_>ym!)9G6K!@ccn6*$gF zsZ08p8&b`WN>HsTrb}VnY=7=i-}gI@JJ^KDGe3!v$a_VSY*&Feb;3j4qA~5p=d98n z8r0rf>OJ6}hi{;>p`UUQJtC8YCB5KjqKpr`z4Y6PYI5)u`gZ`9{Gme814->4jpLSm zjl9rhc?7|;6`F)JgY-cb(^!dFgU{-?IEZE0`h5X*!TY~@zjJiO(+OUAhv=u%*0C3- zD5NT@`rWBeRo2q{K@Vm&cACf?YWR-_Hk2V_le)q=y27c-pLWPsr2*z{z7G)MY6Ia$j6)UeSMZ1 z;wLSS;CfWA()R5^Z3gY!UefuEr{R*`@2$!ko*7NH8c^%Za*ST;+0*gH%WyK?dH}Pw z8LP8xC;JS>y@2}!i+=5xSii89%5e)r{%+xy$~E=wl1jqa??y&iDsIO-rGKD6UdGP6 zO~*>r8a&>8j*&KygKmwf9>kTVIQqD|~dt)V$;S`$CMjdPfFr+uHh+;wn!qGi>uDC*LKli`C9 z&5@n5w9mFLm7;q2q~9F|(f(SHnvHQT0YlB%YluTlYE2Ov$QIFS!L2!WF+J({P5@SOr&ndkp^vdI6FTO*K!XM?`c{vsBB&|g=uwN_Z>wJfq_;FaEz zFnO0O+OTXvf)@3K(mocHE%b{GZlqB0w; zOvv0+xk+{_I7$<%66AuR(luucL(f~pNFxU)wo(n^&mW(kF3$I4}$l&u|R9979xw3l%z3Yq^LxC_19|Un6&@a~@TE}4NJf(eEW7F2;;d_(s zW{;k-m&d6iCgI?|CPksfWe7Z4!9F$UClA>U7^o(P<>h2JaF6{1;gg!Ezcv>-3#qPH z?4C;@D;HbG$pEv8RRE#A1%U2XGgp!{?crU}cK8pCNw2EPpiI61mBRYD*MMPiN<|%Z z8y+f3;MiP19PI;6s=E6Sm{bQ#MN6F9(L=P@qp%r;@O^v_W_=2^L4&IK_ydqTfa&-} z{sA*Rgc)n?K@ogb@&M&O1_vx&+dZ4bmxt6Jw*XFUj2R%LQSDzeU7dD{s>+*--df?C zNlpOI=+~mR@L)oY4;*rS0D6)0)WSh>l?8safWeJ=c+;!|9wG|KT8KtF5HbG~ru_%K z*8@5Ty(#{?4#MS%T^hFm0c z26gZPlnSVzmhvU8YZV`ut@4!0+h<8w%WQPY1pzT+@A5C*HS)N5}sjc0k801k{zn`!k5FIX{b~QPC z$OlTHA03FTMzGsY?w8;aJv&#L_%wcV^LBmc%=iK7_i5qSz@IR1VerZ(Z^Yd3?=!Vi4F z*G3ZGDA`P+3Pxo_hDF?M;jCJM{%E}g^DK`4j_$q)A3bO%g?}s17lmQA5uZ>mfu%kiM~o7jGa$s%F~3*zNi|?n{CvtI?)66 zlw~kr6kQB1lSS`1P9F@pBLE;^fMNL?4+mH4VLUq%0UZotN%unDjUrtBZhUkQ{+#W~ z)^7H0#~zOobe@iQ`zLG@%)tClZsj~r_g~J6)>k5EAu~H!2M4#Gu(o&X(D#Np3N+#) zYCp3Q<;U2?n&c0 zvb%Y^2)&`iyFHjI326S=U~bG1W`Jn%UYS#o@bWoZQDF}nDCxZ zB4(2lR9^~mE6sMTf?lmK6;GojSPF5Ze|ieZgW+v;eSIou z(AE&X;yvupG2}3e7|GPJ%rh3q?ppJpr3F<$Y%yYq0_-+|BnCXF!gcS8Jx5`~*IU*) zhOaF{jYVqT`cY;ZvT{{o5|=%yXud`EUSM*C@LUTnA%fRBMR^ zp`&pBT`l)Cll&e=k@I=Zi82L?O+E#|O3+67_WgfR7gP$!;bJ#<46$12Q&|*Lr-YLd z)(UPni?`h7uDBkGze_%d|BAIg`Bgf90{GC@PE(y6u==6jY#VD{S2IH$lC2sjuUdAW zih+K@`8{O|kHEiYUKBEx%cWn=qO2^7Z?L+0VTwjL^fyAL9T*zqRa9>>(rCvML^tji z(q9R!-~FVAB9lVCTsFGE z$b8_wqvO?U?Xob-!@SBlb)-6FD-~60$*4S#nyXUgEkA$PB!rpnflm>oI-2ODq3Fsl zwAKa}N9-?v?z~isN!6G)upM3={`Vjw4%Y^Xahc4nvQ+cSU1bHkt;qfTAP7(6{iPrB zc?<$my`zv8Qr9~S{i(66xA=o`5T(0q`j3{BgD#ZC0IFlBH#o?NFo2QE!YiB%Nbr zJly%j@?Iuk{w}Z4H$0TrSRIoNoCWe4w(LapKjbw6d9C@C*U24S=gmsN$dCaYHAX6c z#skC%@N&A?iCDM>hRTE@D%6dV2Ey<+j4#a!%$rrAYSbs!52!3qI0iTl;Qx~EjSjxjGd5Qdl$()V?cQp<8oNH6KZKc5i@x*=wZ;dvH$AQh9sPXWo}_u= z98SBTXv7P|!pO7~)G&&8_LAPGVX;h;&gaZ!#qvGdAnkGg*Gz@vUm6Dg-+gX^i&cx4 z!{+8?>SxGAB zgUa`MIh{vmE0~@&2>OEU)(TdH<0aOq^sVBd4kM{#L%Y(8pc+9s+iriAN^-)ZKM^`@@T0k6?eFO_sx<&q1p-}d-hNDes-}g=l z&feBoWkqhpG&fTedN0A2xbj4B4pWI{h8wYAE5Z^T()~(3M}hSY+rGlwy=!k^NoPIg z(;Qrp=Y1}WyHTlbJl9I1XY>yES`u?3Y`Ts2o)NdmG`!L<)fw@EtXfsa)#pmLeoeVn zRi=TSg}vo@YzxM+OheP+oX^)b-}{0)-=+!$gl$ArN{dRWdieO#E4*Vm+peeF8zen( zvX%KWJ!?vzFF#ln{^#jt4!)Lv72mFQ8pC%UxqVmYw~1Wk4UpqT;gLQAO4v$0s{ zcDcoxW5<&;C-FSOa}aBM0ksLrHh+EQymyginL-bIh&qtM8STf7X7gD=h;vX`_?*>8 zrua@=$k2AUA9_>F5H$cw*W1P}b32RhV~80?uuOjT1a53iOV?1}B8FKy+(3w#r>yx? zQTj)7dX%KpDF}08e%9XGam$xVuv_xjl#+o{ zO*w({>ySYpn)i`P>N}E(K!cZlYdLd*SB6p(B8iSpYHpkg-CnXAo{g(hiPBwr4+1B( zx|kL;B5WC=uQ$k@82IdkYd$+{#NS+trgdbvk^Va z_V)4@G(%B`-36+dH`24s^E=#{&Zg<6z9?l8v3_*2kwJ2TIFpabm?U6*=Ux;xm^RCP_o zB6MG|Q;R%H$U35=UA4eU^+OUcA<{2DHX84(u};XX7QNt=qlu)UY^(7V+ocTxSq|ze zVh8~$ozESqB-*x$_rJJ%79_ZP7Ta`Sr)@Nh`vmwbZQW6F`8hP*Vjesnpc=qwCq)4k znd_eqd8o3!FCEOXhX8ZWNz^FxF;m-4B#<}VzK!qe1*1S3=m|q4k;8%vK3ce)f4oE_CWPa21DoZY)$oAVw=MwEIKzHeq zV3r}`(*I$lLo*))bP*ch03dd5(3bROo+DQHuR+t_)-e7XnLF!|UUb-*w{2@#jj4yP zpy|M=`7nV~{3%E2W3`~~5U&AB1>n~6Si5T+n@Ruy88HMfN;rtey(+i|J|K_!iv%V0 zih`|Uf1CYQQ9M38k*5>~QngD;a({|e5dY;JhoTidvVa}DjcLo4dgRo&PkoY33rmZ} z%u%uw)}3Vkdox*KDcp|1OYg=3C6zz*kBs+!!vD;mX8^38Np7*wH*Ms>m8|&}tHQ%7 z&!y>miYV1Q*N+u{fY7PK!_VEO>8KK#O8K*8yyqS;L@`f_~~v)odcJ_bqV- z8jr}$Jd!yTCU;AUN4leC=mvVy72BJGj!hL(yKO`07gM01V|qE-KDXZB1}so;CBk7J z?66BG>Z~S~TL+RH**72QyaajIOpR1qEmz^+vcYMXk@4D;865ve;++G487wuj&G<|((D{@T#C9>6Po~!;9C#qZ? z?0Ll8{RQaGyv3(lrn+y1U0KOJ_MQ6VEbo^isnkj)zFa5NUyDsgy*`Bsi?bH=l8n;T zie!v`f+0~)7`@5(F(meP^Xdf_;tRO_~OrxTVBgwV5QXgoL|x{v}{EoN~qA+ ziKJI|kWjrze%)OZ=JKA3^Vfi*ntDnMkcUD`^CYo!K{g+8s#&;?^BNKNBRn}eE%qud zb{)h=m^fqptK}x%rjB>PDF!(K~up8c1d5wdFE-+ z<_0U|R6~}L%La0C?cK5Oj)uCH=6lx!E??yMY&|+aO_Si#Ky_PiDXuZVna@*3nNF?y zLcDUgS?m#i`d4-opL59{Ic!e*$={YW872r}y)!Koj=y3yQY^K~k}9}MuBfY1!7}I+ ztx}alb2{MxnJxgRWWU{+x&L%$Cdu=`O#+t}%|BgUNOd}#R<0xWG?&^a|0(Ou7ee9i zsDnAw+nn;Pd^p2O=~5+RFUtiOBbwdPSO6eP6xvVwddu7HI#1 zQEdBY8wp;mhmM0w-6PXkU)Xb>s)LX5+rLFy0<5o1zZ{tt)5}_=g5P{!!fc+4he!dI z>HX(pm+q10>VH8euAhtYa0HcuImEl%x3+fjQ-W4q(&VV7S-@%jq3($9PxYb!d``>O zAL@=zAfGrcZ~YW(_{^1;Ao8E(9jQ$_G_XQ7V~=1Njz<&a{%htGr<32!HG;b1InZr! z=ohR3{DSFN-X3rZ6aj*y!ws{ZFFynrq4`SW-awkOw65!=ZtcBR^RDiI zd_Xdd@h4aWaa5e=%4~e_1s3IQ3tPIE1qs7<8vyNgi%A6!VWw~Pid9ANcWPjYg*pyS z455CtU>GfzOPWz-4rDC*+xs7M7r=kJvT$XO>HdD6iCSfG#Tu_^1{FBIQB0MNQJ7*q8##MjnL4XT9xCV@f&%rUHRZSXL&a>$hV+;C!|4z4$6uNzm?KyCwFiJ{EB76m z|6s6gwp=2pi}O@3FY%6&vEYp;CaFEj!OSWiyaC7QvVvmC`w~bA-2b}%La=kQtdtp4 zX$CYvrq$ND%nxT7cqHJuwWXHA8Ol9=F;!}2nqcQ<_ImU-PzN=;VQ_IJRvTpX)OVDY z-bsKs%X-)rSFec|*x>{QAUxKY_ubwwEM@PVzEZ3DA-1f`?^Ss(2YSK0J zBM9hnB0OADvFp8f>``Bo101m858I>#cA6n9nYVsEN8sUnp*l=~oPEt#S^ zIDwj>2Uk80F)6e`Mm&nxJcS}CKk7uDTeU<9;_Fv4!M|<#7-n1;U%W80nHaGS{2d7E zp5yDrx6grTb@mVhg#+dle1G_&#`2;guEEW`MR6V3hSrR!SD|1~Fj(F%Nd4EhfZF3p zQ;?e=F#9JFHITfC1feY?{SPi-m#6@Sd^9z+6m($#M4X?zlGkg2q8IS!1tZx=*9e!i>BkWNhkC6Pc#D~DU~hSp8m3D zUOjny%9BAUb^+c^z2JgHfeVes5CYMY8sqgIKO%gqCcyc5qX^tZX3bC$D6Nc0Q2Dd67bPaHJ13J*&&H^zLiOf=Y)bD3#TWo*G2pY;Y3Fz9 zDp8%3DEMIJR@Cd=7^&;U-?y6ErCr|LlI!#mE_IFA(g!Xg1(;{(R|xznm}1*{JR*Ln zdzW?FudMZtxNl1z>t#<85geSeRh3V}J6m3D*G3?t$9l@y+{aY>qgWt~bqNNdtv7Q5w!|TI&BTL?Kd93jfh0ml$huA z3=OJP_X>*jI_+A{>PqbV~bVSLYkzFZhieq8x1vxeP!f-3xE^R@gL-~@vK(bX0uZ>$2bv2Cd4qiC_5G!i!cbCU7;@`Sw8}e@; z%~4E+qCOWa9j!hq*%jQWz;njZJmpvBN~&Le zjvlg&JuV^qE{&3R;_=I4Cu$wj5hB*)ux$*9mt4Vy_OtL8!U%RHyb#8;(mn0uiPy*W z-@L`9drYw25}jv|lI^QV8Cd*bDH4|xS`~u#ihq}TmOaXTq(Ucw+4pXvs-f8iX-%A+ zzcg&=;6@C`Rrb17jtZ}4hFl{t$M;7Z0?1fsPWLu>Y4%!cPoSS|J2}C3{f7rfDI8ke zm8@GH9#DRVbr?D%FvGm*zzTB&F)B<8kRjMCf%j8@8sNu$;4QhmR{$&oDwRy#9I2)y zp738Sr!mO@Ha{2-+RE@`IXg;Yo*dk?#083Gx$w2aJq|M4)(oM8Me* zRAb?xQsp7hnGU|ODgf+Mv=Mq=*K4h7{@JaS|4+^CLdph<{m2KP-P6r)0@8z?W1hb_ zWu0<}sv)Z{Vqhx^=-UF!v_j1^Q#$Il$ycE>VAWUa#H?bV7m=mGe}7&8_V@J=3{L+_ zq{{eoD4}Tn+n5Kk3^5Qs+Z(@DGe+rDhbg>tpsYOkEl%wuMZ#l1m8jf2O(IH$s7>PA zkV{9Iu@+3|SFyFPo(x+sOvN)jw&w_eb^4Joum=vM8gG1nJL*0bk?YntSrSr+rgQ^=4gd6R-=!P7bRV{a?B*Vw=i7o; z$NQiE%5KGyc9*PB%cmq~5p0S!-&i+KF^lRI9218DO)Kzpr@smy_ui(;qXSELLCQ_+b?erm4N`JYZ~0d*n_A;PuFS{XuFEu;VI!2Jd?M3 z3dC;oI=RUoAS>h7P}uZW1K9uBkwgHD4?W-1vp`Dts~Pxybr#0~ZS{%ueQ&_1qE>oa|N3Nu*5|n0TU3S8?0yOuZag9T zQFX%O`}8h|13Wm7j4)b{ah(t%jwX7y0=-OJM_x%9LG*yVKfERsO$j!pQ8}1Gz0Lht zX7`sT)j=eLb8gnOXwHLvU=!OH@F%=bH)?D1IdKr0RgNG8t{|4@5X2Lj+!ZO{vf0)` z)c}b92rzaF%I*XuvBR?1|D+!+#7jxO6t{}|VYLqn@CMnAFW*S^ho2^a&7P2=hg|Yy zVMe$p;41#AaA8Hb6;U+Q$n{RW%?&$o=JW&q&hUBNgFRpZjzh1Q>2_5QUP?8*r!G~O zteg}w>5s_!uELR&I{yjG?kWnPv%$81I*18PZ8US;of=H!a6DyX$E0)?9id=toDiwN zDI)YASOcKzDJs4gDaoLXEkpCtp6!L>*8KJ&=R)FEfAir-0aF@JRwW*MWZIR=t!(%K`Z8BjGD zU>;L`M0I##@T)*DyZGs2A;0L*#?ZaL%VmLKxI@YaFbw}kPhNkgMrz-uwD7vT;&vX0 zHV|NRgkI1utt{{fXA}VDFYsUHFRSL6y9H~ARc7t2f`8S?$`h{Gokn%6fwlZsXC6=I zr_9|QiJ;vup8dhOE-z|q&)c@$oZp;gxfy|L{M4Z&IuE=Y3~nHaj#k|X#(bA)0O%-6 zvIFd@6Ns9}YXk7v2{N8?e~5-N05ks+29}dwEqFWPSBR#0TW5Kdb?w+&5tFyw^dZy^ z8c(DHjY>Jl(_U;pBn9$@tJy=sqzi>f{elPbRSmmr{me`r4wyYZ< z3w>hxuc;^)sI|q=rS6oXq_u;VyC&yuxjeRi97hMs;I8e}jhG?l z9727=+;by?VRHobi3KmX9bZ*-LH5!aak+@bPopQbj@`U@*7rTvw{P&3aK+gf9N)*b zOqoojsc|~#59HRtt3|ZelecQBzdzd&f$d2bEk-{$=1GD4yx7-=GEbUrddaqx3-^k& zf!VsQH8OCPWdUW=|G$L)e{BiZ1(m-S+h~e`w6*-FyH-!5@8r!^%reW=EAE#sfvgZqi#z9j!{u0RiC?r5?whb7I zW?)Ded|$z~RYz}KlWb2gH@Qs5hr>;AN;8;G->Kx*uDd+c3exAv8`!;Hv<8sORpuI%;I(Zf@2G=PLMO8{ukFI{W) zE^1Vr75fPrG$*wU24DAiM=jXW=qFHeyf!{8C8ds^t-NB#Bg>yv)c*-D2~s(e_oqJw z&wwv$>LvpL!Q4Kd06?&s^SGjkS7s+~loY*amu|i`ja57}?Pr4M=s=xzK>|n`q7ExQ zS(AWm!kN{JtQ)0RA}!0p*abNw*dCL>623g^4Vu5na_O&{LfYAH;CE4 zc9NH`YeL~u>L9loo1AL-T2vSr=~fTUeS>_uuGnK1^CIU$x3JW4b=LPy6E{-?nl6Ri z*te$QaCq(h6}C;zgTOzPhZ>=Ns5c+_4${o=Le8kYVj^J`KjN}nPp7vQZj^rO(-@HR z>MUC&+E7&Fjq&KyYVO_lR}t2aqnOAm{6Y1g#B z)jS=cF|&~#@c6DFl6s*j>N4(K-V{l{!^6C`rjhLH5`Gd7j?W&7q!V+Fa@4M#vwg-xiEctb=MVRmo z4Ba@FkPvdpDdHw$^da$R?9s&kksi-*24^^>$DuKAmtRIu{+q_ovQ>Z|z zgR0{U^SA%NrjcF>7MpTQ3pEby$w~}4>Rb9JVcu%dpT;Pv?oFO zeCHcbX=``}Tp*u)tNKyaDk?(t-eyY!Kj!&xTHmA zpIq7J5t{Oi*Yf-1X{w`ezRwy-2@@Oo(T8~-$t;=`6uD&q>+m@mb-=B+-&Y-s+XEFv zzF6XO0IECH`r8ZNt(*LWecO*M%)=en?NM#f0%k^=;4K1Jo8SRp6>)ZLviius*~G*0Tp>R}?A5V7fA-nJz{00jr9({0(aek)*NBmtt z^eZDeK?wCj*^@cx^(y^cP%??w#6#aAQ40fEEmpl4Fhlm?0+L{slrFLWM_=cpAG@#K|A4eK4 z?rGDN^j~1O9MxpK6|&ks2aCB-;pk)#fC(tHQlT-x#ZXAaRgK(cSd;9R_;Czo!)cCx z^~3d`FBL1}0>FHoYTMCIYIkQIv@FS;xPZOI{KQwG#co>sdjt zYc?$++>taP6K0~c8JRe+OZj+Eb|Iy_%2r`3Twn$3-o@YcLe=Ndw=kNBW+9TGTsL`H z?07AjS~wuu6uYozb>qYJpz=~aCPAUdX{YnHD%R657;Gz?0_Rji2QP?8$u!Ve%=^Sw z)ve?@>ou@O8{cBE-*WE7XJ6mpATsxkBHWjTX?3h^*<&;KI~SMq5gBq@Ql~;|8e{DB zyUT`8C&v95k5-@R>aMS539(Y;cqu&`qb47>QW2t36xPSQ>F{hAM*3nSwZ9u`W2N%G zbC;9u)2nK!pcak9C)1?iyqBpQ{M(&C9({hIrwH+CjJo@do^?#%Xk4eTi%DVRM1ifA z$+su2H8+g?6#cF?O6TXgA#QH4hT5bo#aTNp?zZPqR@!#*nrLKBKc*+q331Mll$kW3 zsCRwR|FHW)RjALTy-Z`|$)V{}WkZ5~f+4J2q(NA+=aRQ3!f1KIVTg_YuUl%^%z;PbqhwQB>sNZM@EN z&G?BNR`AUu1bT`&@5Sl|(WXbyR* zLZ49(u4Y)D#?RYq8E@{HiEi#$52Gg=MjS#MZnBt0WyoYN^oggHH+?o7ZggO5f<4Op z;48{uM%L5Y!u88w8v>@bup60eP+@Oda*8?a+el@huNkTj!7fU*`@NvB$C;&W*e-HA zwV3KdW;rQ!a39p_|Em$0>mLue%YK9z0X5VgeT_U>@wfl62^u?r8-bJO1tA6$BY(9R zcWMUT&Tj#U2kYif7;xTx+@FF7%uZxDRWKQN8ngm~i1M5VWF1L+sJ)@XI;TUWt zW=KFGBKJ@kc_1|W{ROASNV-1Sw>h|-pr7mq<(RB6wvz$j0?Ac^5bYQUwg|fMy&=>F zFhHp2Ld;L$A<#x}Id55`-l&c%?+FM&Dj!fCV_>aBo%_AdnKX!}2#?>5!bL!sW^3$b z3MdU;417_ycfB7@R*7&b@t*s5W1c<+0kORYhqy9{*(dg*ei_#07=228XerBH!*irM z=3oAniK_9(mJmP<6R`(2J`>Sqt_c2h!G}0& zwZ-qp5;MHw)ZJZC<~ow!iK+0!Qu?Dw9@8w#BRt_fUKG2)c#%;h6jA7fS;5N{TK z1ae^QQ_tkBUVqN5Lbof4so(#8K|#cwGslj|ZfTpJBbnn8gT1bdiU{ilentb;mrym2 zlP%9!jGMn|>EVSM=lJ^c`G-qFL+pbtc+WPE#*pz3Pa3ADF2$u2*df7+E+GTOiMAVw zc4iFM0`WDR?QEhP$2E__%0>$xkR`bBTP+|=!|A0HI0WPSs&drHAmO9;waj{jRARR&<`AEO-USO{ftL8 z-BT*wB$LmwCdL|E8U;tLV!}MMq|h~0Dz(~lN)s4{(S+N{jSfq>HrBI|(&F zlJB>ppMK@ud(Qbi_dMshf2b3j874D(_I}^>u6M1O4ZoHV(N)&g5uxTc5GF9Xhnl^i z*pM)EYf0Qm)A}xZz9V6&kyJtFv=kLY6 zbFtqB=1ND@koJW)np?{o>9KCXUNt$l*C5B;#n44*hvaRS(J7X*)HCN)ZX(46qwoDCNB$x?c$Q7)w&=`Ra$ zpE;ddQg!db+Y}3$yWzrx5YwoK%1{Tz>q|=9AE^2M9odv@={BSGXr7=}IHAPAs2b;# zZVy2sbw_n7qPY*~`oGT&Z0#498m;i3=s5G)VTyo4CK zitluy>t04;_mSHj2;m!ykxLWc7(6I*?odSldBjE&O_SUM7>JliwT=6&u&pOrhqk=_ zjxYQlx6g(Y$3-91Z+V>!ctIH_HL7|?>ZYK>1X4%1s-91E8 z(3M%U3Z9Px*DykMt2aOh`LFJopsN{*D?$Ch@anIQn)my}V8;5TK}xTaZPd%MF-np0 z#qd#!ulguVm<4i+u@3}F*OO`W{#*eCXrnu=9Fd*Aqf%0CVY6pD^={R4ahdfCON#>V)alm6`*()kEjAt%D z*D2wmu=V_TrG?uWX?ulWH-w~g196||cNbTKitO#*A5LUHyPs`O!!Ze&j)lbjt9jw@ zf5T^LQQWsQ_)-m3k*{ju=c*WIa!1xG74{<=GmLdF?87B?^TA8}Qn(Se#cvJ6kpXP* zEd2ggrDhwdCmHo&>@&H-wM4Fv*&v*wu6n_nI7SM4Vy6`tle|-c`w%N6j)DuP+fS?G zkdm_YtgkEUmguZ}r@OAd^i8(Z6Ey*uS;61cJXge_e@-U(C`U2lddErEz`9WmcEO(clnbUkF{`=>{n+*0Ou z)awoPO!)injHp7MDqaO`te!E5c`$!WTq@24MwI^=)o4E z?(>+M+yFP?yq$(_Uu^C6x18 zUZ8ls3YY(S19}zQF8(R#K4j^~?a~2d{PnZ#uEEd%5f?T+Mngt5M`EX!y)48N4EbZa zf#BuVA4Fgd1NF*k!AV3wNi|nfbi|T@%tx~mmR&)dokdNlRV`gYU#KS8u4hkJfRz4c zaqsuP<**w8?W2!fv} zDgEv_(u|#WQzbZUOf82ENk~>U&q0OR*{3iHg-2;Ko)J|z5-oj-`FuF#P>OuEWsy~c zrM(q>XvVF$*<$^7Z9&a^F)dJ`4SD-zBa;No_zs*_59izH&5eXt^nEUT`K{TG>W*&| zx6p}4nN*^0)!oF;&28zUO9*M6ruO!HKwK^)pCkN=w(RjIiX0dS%bd^khuu1T&L@Jh zJ@{|+&Uwgxb7T?bE#?ouGRczJ7UsL)#oM@L_Zo4gWc`ShH%~zUre1%sYxx>qL{7@Z z>>}=oVxL@O6*UDUp3HGFc)KQ@q|3MC;1*S8@IvLpD!aqA=p(As?B_#D<`UTou!eY= zq5LfyD}}A%>PIzw7t5GuRKnKZDja%pKW}!@Mu%gYoNN(6fs2Bkb5*T3c`>f+R4I3f z;ru(Zhi`-=gWI0EwltHaK1I;Wu(Jx+7Qk_04LU-`0v}m~t5e>NuF8dI9HMy2czGY{ zUXP2+7#{WR9(hzlT^fG!sK=}1a86%>(T$qh^WriN8T`l?)T_0*fM)0^*V_sh0YCeX zdW)r-)pfa5QK>F!cC%;n-a2x(yjhokT>-5>U||4H|7a^QXTZIXYrjZ)@#Fu>ph+Tk zHrXX!=fK0>XVd-OGJIDi#GMq_KiY-*ofBDsOlp7oN6%=sPN^h&K`W}imAFXMjx;Dj zWJVDB@Yy^ui}!2j>B6!=pwJ6aBUntXd&>Y~dJ!=S?dWkIR+}Qxa^%w<0{TpokbmmW zghb0~=hv3C!is}7Zdk|6a?ngxHQ*H22VFZ%08uFlps+Tef^#Rb;om<9t9?bJt%lBc zsmRi5Lxuib0sZ&q+s+_ua|Dr)C{epZSMBegG|Vq?InmC530YOGiWfN9y@6p(mrz1V zB(;%6|1v8ZiAqiD;W{7>i7-}Ur^D+KoG8O&)~D6`U+*HAMhuTKPqWI>YOy_kM!n68 z0~k3BbAk z^wqXger$e^2h6<$u#5prYKbYfn+OB#a=y& zXNX57;hMr+u^VPBkoeDBNuh0}+#gQC)o?CoW5L^~ul!CqmeDvd`cFp67q3Da$gYlz zQpY^?)|S|l6Kbvv4QRK;P@ifp95ae=y!s^L{DZO?>i08bzDGEl;NJ`EOm=2EqYOAK8#(EPs@8S$|#`TQ5EuEoZpDZ}8R?NB9)}NBQ$&Ssx zGC0w|Mvduy0Sr)GD~*Tp*$>9AcvQ|l*{^Yz*}tK&~C+V2#w zU!gq|G96w%a{JMlD{eEVrSqE;IZ)A4Gb3!ddCRe)!PXY}rl$hK`S#!BWxv!A{u~r_ zPpwV*Q`H%R6x^kqAWUCl#idj;7yIia4RZ%7AzWv{o226N)TJ;)_RIMfQjhLY*r$5x z_KF03cPjmgzj5s#F0`mVTP9I9xe&pxik~^^$=z~}!<=Or+6bj{7NqvnD#eT%(2dO{ z+SIp|?gVr8uc&{&;yEcq#mx~y*}7~&O%tHk3gvt&-CdEMy`Z&kOOSvEr}o*r6Y zoc9;{T}q1{>lkDg^Lf5DlmT8 zrVIPvc>@zwJ*ACx32`p(UHjJ3TkKq#sF1?|^LlkxIvp}!Mv-yPhi(3&u?gJO<8D>I z4YAwndWy#A?Sk}Z_Vav|YV71;zAd?w`47D8G8YGjegr}496dF?V(bvD$qk`n%R1bL zne%a*NaE+Rjl72y9)E2|NL-cr@LN)|8~y0PACMrZ_Fgzx+~=WpMJhEdJ4fP7niz1r z`%7-;@~2SCAGAmdWk`xzC2p{H_{zk&CGMmTTPF$-negVl{zw++k}j6$%;^>>BevvL zA|pzR?Tt}m<&COo=?z%{3)s9QCCDP(UsA;p?jqNav-;ZEV*4?$F<9HDY*8? znrQ|uGC&6SUjHBhdjP^_V}`xKG;_oMq`>-XQuv?y*VQB#D>$W7J(;W4s6)?o&icB< z*kiwN=XWo)pP4c5yrnDJeLgvP*2hcuCU!hD8z-oW6E;ns)m%enj!Y{RW*Y4AZM4Fam-1!HGBut~o$ok!fkhuD$8lFO2 zTq+J$|JI^$eWj*{qhG`5=ur(o!LSJpq2}6vY_;)r6PaBn^ssWtQz~3WC53?pcPi;u z`Cr?!#0aDxYL!0M+T}(ier~>?WXQ>5RnkkYTf4Bl(AOow@noH}{VH>uBiVDGn-UA$ zWVZtA2Rx0F7G=(sV}-|?9aHP#s+Q%>9zV7ETr0Ddv&N*HBJ z5z`ctKANOxeJWYHx9#gEUWze>a}wg@msyXk-MB6kSl_p8z0A*Koll<@ms(wvrD)GO zZZmsCbJwR_!>5vNWOT+^f_xZpX5s3nfniNqTF!1ps5P_kdCP&W(;vF@u1?a#?|e83 zB03++yP@j0z53L9uH`9)EnTJx^L?>I{W(pWIdpp!hv2a^@=Oazfa146x~{=EMe3&O zFUY9~oF|Y}5i{>FiOT#9Kp&b0Rc`(K!$oZnoy6= zK2;n-r=00$r0F}n9lGizOnf1VB#E-bj?9GzoXpa{yd`!Fmg}bS*b0Wlpw_78tfpTC zdAwJ?@=kljLw<({;^&*z)IS>T4iIwSHz)HHcxbRrj=7@jaIZwp2#QB05r7HV^{6d2kQrcvVAwPOQ=>!o&uxX`PlW=kJU@n) zzUXOTwRD{(s=~_}^mzsN!xl_?VHimmkte+vd&v*s6gB8RnBe$y?;EX{!rCySIAjKM z9M3x=xk~cgZB^-o{m_Th5%U+&2-g|iw+-aZlzQ>+9Qa?KkJ+N8yFg>=MBC@UC0p9+ z;*4azmK$ui9Oz1W_a69qq14!vE(5Gh_Qg~nf|w2a7L6oHiLpVal`71J@85;tFyJ1% zPH$`d$jDP@<9z3ucK*t8{@(15Lz53Amx*qO>3Y1i@3-zEnShbo6_-IWVPqq^0$?XW zyelG@YsY2+j|FpA0B$JTzc>8u29{~ShK}dX0V-IH5a$B$3IrBSytdVs=QH=?2&<8T zd%Z4Z?mXQ)UPJheABTP=AVEGOjRJd$h{?N1dxYO4rXry88k$KZF%KbW)cu;oPQhK6 zx{1>?7uU!6Zhavh*F{y(=`l)>6~8GUv%bU|TOv`M6XR5DJL9C4cn3K~kY)Sv_k*7Z z;a;>uFdb}oo;bh6eZcYu1RUrt{b@zzDh9kj zgWdtY^d*MC4j$Ojubtt?0b6gm zP}9T-VnzNY-S8kEI=_Om1U!Lq#E-Ob(K_}YpWLN{I~vwrl<*3a@Ro*Zb~ZRswt%cX zufm*8^|{e%H{MHycz}>x$WHhSPLhe5u$wo&Rwlu_@+v!W5H`d66Y8fM(B;iRwCLFV zAT;EKHas1~DT(+3Vooz21+0^*^?;K-q8(`E8^7C?n9|N+wfmFRVC`PujQW%9kd=*z z>jlCt+tkmY8`7Y6XIKQ>rAOMc+e^^Huo>5CPrXL91YxPI&nu({2PtuUDesqa-;hp669xi5z|q@=IHg8 zOl)EIOYP_P2*5 z9Moi?k&5{^M+K~E)QD|~DeKtBrw8XEdFnsvXiz35pAR|C3Vgu0aD`dJW1hM=kF=YP zEYVR0vDZ74TVE;iF^fGT_v2!<_w9MpJUzyHV7K!cGnlY*uuLtbWVaP}yecySYH1f) zU>|?~;H=4OJ4|2gPG%|x6Q6xB4>~Vrr}O;f%VrN(T+F)EWp-TeN!n|;MxZDBJdHD& z5qXwy>xz&Il&BHT&|}BXKR9R&xkWgi?yLIc*+}vL^{0%@U3iy@k>Ubg-ZTjt=dmbWFNaX~$5Js~R~U<7;@8 zCP4aCKTp<&$s?>AYM+}Yp(pCm5>-*MrDfqyS{Jl)efm=8zbQCO(0_h5dyRxxpfNczk$GNA?m1sV*>?Z z#tnwc>RaAA0DVMeggD&_VA;C26S$j?)0F`nk`KDXyUL(Kt_~oB`8wnlC5QxGBXC*s z`~6bFaR+QDXMn|pZO%?$yhRLZ0ulb>XFP3t9YgosqCwaJzI5Fu41#bNw(;1GC{U1~ zec70)nKvgg57@lQxa1h@2NCI_eZj0l1Jn?O5m;vmJP(2#(E|dGOW@$CC6)-W)aU`f zwV~||nMZkjXc+z-aJC>>r00XD;HNejx_=lk(?5YXHtzLi!x`T0NB0+x3AR)QxI}Y} za6U0pN@RSspk@`9aJo)aNTt|mOI=oR;b%GQS6S@;fqz|*&H6buRoWP$*F8BW+`5^^ zM-Gqv2DZa{My&O=uC<*WW^Skk#0hK3-8WUo5o!8+I~67$9s`3phBh^!v+-BSF<`mm zCX2v$GRdK^UzQWe0;S;8!}{d?v7DmVZuV@y;|vsx#}s&yubv$`Qt{f;ekz{^*R0o5 zb=EZ5jY)XV#ib%8eASJbneOyaa-lZiO|mixmdDrC>PB|`Ww>Rhkd)uLJ<^t9m6mS?cX4I?lZQpZMS}gm zj3#6~bvdOlIZaMhgiQ~$lVQa3VYc>JN+qMZi}Z3fFIddIxzV8Tk;QrV*vuO?g0*I| zH<;Z7r}8N*pWQQ1JX+#zvu5zm!n?}e$Z~En2j4}}Yi|v6g88{Dr z4%2RUsM4g-=2vivs#-7aRzz*x*C|kLIA{dT)ZQqslv-;`(+s}Q^hQUPU~{Hc8cVjj zp;wXCtQp6&Wr-IpfL8~=nK!ya5gQl3A%>M;1Ui1tWw@-T2Y3bP+lh7HC0)q0o;SF= z`w0+^@$DnONt`)eze6Bh5%{a!rxS-rHQ2uboivaoT+l;I!uOBo<$dN^qw~J#AqZY= z5;&fXZ+s@wtdD3Wpa{!QU{v5Y1E4s+U<`(?-xEm&~(?hlcCfpm``&_$)uO`#t%k^(oV75-CH|3kcn zEn6{AiUEH&C5nji`)*x3@|qhI&8c?*y%?0KO!~_FaJ0sWXHV3@?ZQM_pnj-wAn`nmp5`BJ)jmSivcr~0Z7P-rmr&*x)F8uUn))UMjwQUH-PL0}As=nhT`4 zcKB#8ufI+v(#yXDgN+{$@4*6u6@lw!d_~QegRL z&_tlkaZ*HI=pMfqcvY2$wTZgB)XTJ}7_SB3vMo=fNBE)^u7cC`A+U55c>a7cQnCyL z$S#5FzAULtk44xkm53Ebk3TEcr2|bz)?YyM3SRqGMbLCqnc{f|lUJE_J8iKd9qA&v zE&wZP`&<_#{pooR+qQJL@zJcC3$c=c+({_~|Mz*%zw@B~_RnPPGXF~|bp9ulp}J-2t4o=ekk#|4vf^gpw~?`v;H*ESydIEgAMw} z{B%S(^U=pY8s1R>De4v+&3x zot=R#XrypU0MU($u)%B5)-0Xe@$hfcP36Qz=-v!yXUOREi#wx|LY}C2H58NkKwti2 zu+k}I$LIf;dy4*mn14fbC!h%XKyDBz+neIBo~B>Kmij>h5ot%I5mX)GHkrz4oe2G% zioxX&4*+v)i?R?9w4A_=HFdSS^^9*sAil#RP7Bss8uEa4M_Au8dPL^Q-x4;!zre{k zQSToDz|e6hP7_*Z!W3Dih+gIFx*GNA4aoV6=ikdIji%OLIiRjm^5BETh}No1LdvjyM|(l5-IK6FBi}a}ERRlI&rVFc ztSTTXc_F|)+~t}@Qs9*2Kj55l5WxXiA9QyBQ)zC>s8`}XHCKD|!yK*@2 z7UaWl-~WJc1Q>&r0zu#o@e%hjZ1k4FnEQJQPxDTkxtVKO5gKlLZIREkjJ+eMgOc^A zD#R}C&1n?xVOftO&4EsU^={y5<1`*`{E$ z-DxX>kfO})BazjesnxTK7W^E%%)V!QjzKaGDYBEjeab7?xVz(SzRHWXZa*mJ=Qr&| zdiJ%j%eH(b<8x)?pltrftaY|LydGtkczBQ>#NQ)a4Hc=e zs!?T38t?C7&wcChde7@8y?UQaz&)GnxzU_&I!<7@8BQaB*GI6v+=-13*E7nUzA3FT z*P%Z}>99umk?drb_&hn;Af9(xEEFx5FqO~Hz9V;)%wFuI#l;<`$9y7bV(;;I%lwzHkO$sBVrq zmi095R-v_3p>=H8?W}-%c|&MH{JT}?k7OZ_rxN}WwV_#zH~=O7B$Llpp42_dsn~qY zMR2iI^+S{-)cQaAb2&nBrk10;tC2s{h_V?ZQvt*T#1CQk)K+;%J@o{Yhrt!RXNhN^ z`_3T6LK3Wi@;2}(fa70-LTWnW<;pwS++Fc*X*vJeO&fc3Ui&~m$esvhnOp-&2RXK0 zl3SkCf$lFgdsqo#F_Z|FguxuI6psLD6(-)R+g{~_CPyM5Hygwez^njbOl1I=Hj>q?nc9o;v|=zCWx=SPPmfP(W+Ieis7xA`@~ z_*nURZ1?iD%HZBuU5gS~;f|f;j-B+?#odb?UUUhUvOmeGaY-pMj)Cg|IT{H3`l#^; zaP8Ww4@Nz6`jU#K2c(;tBSK)+V2?~0aEW(8jIdN*W9B5Hs`k+#Fde`g`{p2w(@L5@ z&zj)17N1>6fTm@S> zdw2%4VfJfV0cnwB73n+^MUH~kF~dN{!YEu#A{RIjp@#VYbq`fk(@HVyna{p5uD5}& zviVk)HfkU8J_ihrJ9sfd0wDoV-H4GpP_ga1d)ctvh(oc zlbrqv=Rp0D>kkN>VZY{15!~|i*RfT0sa`f zpS*MB_ZUV)RESKjH?frpU9iPDI8h)DtnBr*6I&D@5{*IODwZ=+SHT)Ce9*k@4tgkO z9AKCnY4E#mAn-9@-=GW_S_20q`XKkjKp@OuoJ*Tv5b6f2N|lQp;awBY$7$7oQ+?nf zs&s;j9~ff?;G+}22CK44FFvT;#m&JhNvBCG6a7Bq{Ahjm$LB8Csi{K{?SI#%W z&zKW#v@RTV6DJGjCM*6&l>8ZW#C1poI~OdPnOf^YMPnpV$=}OUct7hKv%de8X;d!# zI1H1z`-;Dqm*1gx{odzWSb3|vSvROujc?v~nlky;-TjS^(>pM_ZhC;yju>>i6&2|~ zH#80q?BHgCgWar=-vdJ@cRi?|--Xl^5k2_SQByqidriPY{=1fs9B8B?Qv#GOmLv%e zEp5K)ES~j`_-oI4ai_JjvURe9t{%W(JKjQc1(5>N~s0qS@eeKv>J7NKWr}LRQPKCWdJLMo~E33p(W{SWo zXKwh2UqWV|AS)V=^zu#$>oWYuo%+A7sr+~45>rHW*KJRl)bE_RDbXW?Npq}{U)A4B zVU9m>xWR|wo_!554Tnf6@wTrQ!voz_11@DkgD3m*b040qd4J7yM^%^VD{scQ!yz)t z{<-xNwIj`BmNTd193V2jNuOa9xL#IA^lU(&b3k#QhSmcq>Dr5O`c!VC zV`Nt*&o)g8xg&*fg;-6+0X#DW2!JmoN-dNe!%q}lEA2hk;jr(=_xVKpBeJsVvJj6( zK$1Ts&UkHXdHFV>UR|qA=Q3~jI^4jLl)iWO^zr0OEUTDFhVU$u*NPe)+=WjNPO=-{ z$_~dtZ);hB@Tb2FjbG+^6|Q?t>)|)V!L1>pJD8M!`#Z4*tip*X+#+d4NGjN0Rc#qs z>RIcI9(eK=omTrJ+VK-Rcvw#a@4z%L@9r8y%CdKM`vl7uIT8?V8^k=}OGf z8@+!hVnYn41;vCfLN|8-2r_}^jRipyz=Q$h!c(8&=BND?0T)gS(O0r8oj?=iL2ih z4!RKd$%tSD-=G4iKNJ&V0p0v9ZDHNUoxQ)U3XX_HJVYO$>hYmqHTQdshVN4%W~xx6 z0}7?VjD1t>CpM6OEGaZ!Yb7Lf8*vifdoZ-eHY?UJLGgWc)7h)C0+D zHu9dPhVp%IqNw&n8kf}x?(s0eNHox|LONmjmwF|5&9okY4P#CTbwdM{OUA*c6I>dgRL&}P8bJqY8hx_~>>=msu6MNz${KMuO&MB@)Wn~SZ_dU` zVTJHo)!1jgh-5jz+(9zVExwI!tU=6*XP%eD^E^6pLn=n0`O}glNfQW(4n|<8Xu0c{afKP<` zT4+h*PEfUmhbHXX^n~9=zquo8_+iShO+of>YCqIzTQEXttz}p%hG1pMFE}cS`Ue{f z*6;GPq1dG3`(OX7MN%%IcTH7P5Vp@58ora~9 zPfx$S!R}c96o)Qhyo6O5^e@e%QOMoOn3XXNKecNq!F=m2WVvdd^q6<-!CT>$ff&|N zBhb-l_Fqd7lLK^pUc1!Aov16KG3>?93UZVT&01P1Azcv09lBx9OBkKhNf_mrrF+Uw zgc?8Sku<0j^5YKj{j4?5?cuH~=rNi9x@=s!uj&CGk%r07?{it8hzF88^9)tN_jn?< zf~Lp&JX@=*Zoz^1j&s>gXxPsBWtp!#o_pSpkFLhd+d9743@p5H8WfFxuh-7IDC{)2 zRHL{(-3xj;V{B0r=m!dxmM3)>aL?ZZ`R6%W9BPQ`pe{x&h~R*z;{6(3*A3|f`CBu0 zh{*t*vY=I9MagclVa$!Zu}fI-qW-)$b25{n{EzU~SLF2(8kp z8?ZbE%m^%25VXhjJ9GTC*;9Ir6THHD-jR>94FRzgrc`A7&ikyw3LRrLL0BH382^af||~ zKjlU~PQT^~i1AQF<=Z69+$2wcg@pm77gXOL5ZP7}RudC}lxdCDqW(I~mziTJT_f#N z@JwG*Kr^^#S-PpvP1HC zD%@2FD+3P}fKJ3gSdXD=kzx1}hxKT-IJui*A53WhR$4;62M%av!1MY65HWNZ_JjYh zw6ri~UA!)~Tf|d*nJYziG??2tXeU<0a^nQVTXuV2r|dfy5#YwWao|z?>40P1TQ``= zvCg#{e=hK;ivC4JsTu)bzuK#-%lg{sCTUD6<13Rn;Uxifb_Q(+fRD2cey==EC*i@6 z6J6Ic{lEQ)W1pB7eI{!Zt=&S_4u_{&>WVpn`hp{=zEJv7_rB~pPQe-QvB52~Cerz4 z;M(CYhzirFk)Su!mxPBAMLNsT10WB${3q|C_yY5ah6n)0L2Dcd7)Rh$9AncPPz;7k zLXGtKB@E}Fk6-Oj9276jIOp%3=mvr?;L2s4dQ*>X2VP5ya8IDd^9wQ}ddvZk{Qx8% zyn-6I%IP={KZTOcta~h8;bf!0BTvscC>wmIksr;PmXg`Qt;zAlA@)H)_ZQU0rK3+!pER=C8Ci>SU9_94w}@8ajC{4HX)cfwey9y! z6@K)j+kH$^CK~ZA)FeRHcV(9$+TJ+ln`+n6egCsYm7%=RJ(rPGE9#8+7T3 zOkY1Ha#6E{4d2OmbLu9!z)gzy!PITP#B_7c1VReKvrnd;!c81^W7(0xwNocwF!?86#AWY%YE-% zD1L9&!Jv?R_UtS+=ICv*u4j;=3L@J^r*B)^(N$Ne*OM!p5ud)76?hZ^c?w~9n7vq` z8kcNXQWKp%qW7iCIJv~{p~f&72RE5G#C-T3#cRCWP+t=i^2FKrFm>T@y$hwuw?>}l zwts@H3}>kUpJ?VRn*k}OMmjMLiRqf74SN0UGJa_)gYJ4OHgfc+$jGPm-3zJKup zvWm(kzKbIrAN)9MYR?NpJZxo(4}At$3CWQs5IFK2jVtuxF=Kw0@G6TG%fg3i-u00` zanMC(=A^5%>N7B=bnJ&9w)t@g$-aYY@KJ&ayFp$1ix{Y;zO-~~VNxz_2PLoLfy@wz zSJ#ZUAa<94xEB5+$**71mP&@_)E7wvhwMqX5u_2$ISeguiOE6CZ!4%q66u88SX zqKF`wW(BcEuCS?3Y0#YWw|bSy+Y!-3kn+#a_->Ywhz`M)&`l6Ahe7dHY;f&u_)yDY za!l_F%XAzfW`5wWsMQ^bXE%_??;~d90a6g6cHF|DbiS~)2{wqr^?cgTTtnrDRXGM8 zAMeJ1?7P(>Ul*VTVmAP4K#dRP0?R(|K&_y_;}A2JfbW3H=qDYJzB33d!oxa9d(%k~J$NVS`j+8|`cGInQ~3sI*h#7`Ae<_8si}O{mMi z&n+^SZv6pKe2-Y!2BJ#HBZwY`Uo3{fQhPrnVYf0mR^VY$g^}PWm zEjFyHmkF#I1vZHJlnD?8P{RJ&%%XXkF#K@2E>}Qr&GFsO-~cQG1^mDE)$YmF_C0>7 zMPdn7S23%>ueVYCF7TR8{>~%Vcm(JGlK-_~;_@lXq)YA^$p6P|c#bxir~zG6F!&}8 zHi+PI@*-vf?=Gb6nLd7~o%NO2Tt>vBNY@f^APQYq+zZ9zJu#Nho}FEA(#$?LKBsL0 z*4^23G{Fk7GK{F^nuFD)`uB>9+3-0oSqDF`x}>f2_1zP{q?F<73&Yo!h`GYC3-j|} z{d$d>H=3sf7PV8+>+PhmHd?Om! z&nW#z`1F6}|29ye6{V@kuJldQ%KZW><~FL=fGnQH`~1kuwLCYQZiw&DEky#fxYz|B z0VKHU3Q2qdi^GZOY$cO|fhM$lTS(ev$x(x%JzaP2^ryD5;$<0@Pp=l+jNZO)+V&(| zjB+bTkZBP_n*E+oS_ad51!W9pT3Tw9D{?vm!^pCE*HHw&0~;)bMpjE=kNXAEC|1At1eTK3C)4q zcBhwPITISv_$BM*t4}A@y?oIDdU9;Ogwx_v2)!Tei44P|uA%C_)CVK<6qu7M#dVV(eQP(N>%J)u59nnvzRpIi^N73y!efW?{)!$q`+G1`(me zuLH%yIopQ8%Bq#3FWSb(9>p%!K5e4a;f%cVq#>dF+}e#iM19l{^0bTqSLkYV3uZ0) ztp4IgD~)c=p0U0U*P4G=*P0~Hv{Y#`#{Zb%6uhC*KTufG%IwUoWQS}OiTl*|j%#zc zi~x$z5T8KJSh}H!SE2_(h9kPQTqRt-oUCN1h-xyVDF;jkZWB%LJz9R)G6A+_6p0>K zgE1mz?tv&TA53pw9drY9!`=`tFN0e0OQSx-H2Q~lMylE3U-4(`1mj{aPT9Z(k7xu@ zBD>*$Z#tqGPx2qCHqx$Pw}>!6G_OFnzakQVzQ*6SokSmtl--E}QT-;v!SvqX>H`qU zH+R6-LxKO{5?Y4dZ5YUBBs|tha33}y*$?Tk+Ul^n({??k>jaTh^XSBS*=An&=Sn2* zAdUbSCRJkOvPL!efZL*onuP5i1Cpd(aN@0qpTc3MM^iNN`iq|H9T3E$tB+W&AoL3pFB&+8p`=!U_XZWpMtFKv@l^u~R^N1T&ZED3Mi< z&?|UdT7?tu9<{>&NZ!_RViOB&&;ZhO;I6zX2vF$k1Hm<*fjXVHn3Wpb4S9qIKotX? zV+4jLB4)rkN~uO8YbUb%F-OF_7&_W%#mhxJV*@23U>hJ2@ss5WD&~%}p13~^TLNp+ z>(uR7Xjjb^ENckyfY9+hsQ`f%?ty{16PSmCn4#e-$U+jIgDeEdtQ`faOQE^>fCWLc=i41c1(OkpIA%kT}yH&N`q|Gc>2?AA6@JeRh+%h|}S-8?aAA zE|7>R*Sc}~Cg5Zw7a#zSGz{=;CHBJIInma+lf*gTg?<^t*t~VHABI1cKoA=sQ%GOf zg78mdigv)t&Qh=n6kY9ro^1ii1OaE@O$E-^O;iE^eS>xN`CNOs|a_3C)V2e>&#>Io9F-MlgiU1shG2J)s72BsBs9r>nD< z#FQP?H`i>HkfwdZS5qMmu3Wit^6jxTYUtHGph7hifSBLa0(#qp;BM2+PhHu}l(9(k zXFf~$bT(|-zz*;ELq)JZU&j|s9Wd#z=$08SCqa_MN(bXb|0Ccq4JQ97t(8X^0Q#NU&4RdZ$aW^!KCo=1bCdtpza?k=5v{m4SYmZlUxrL*(F>T(G$1i`8etE*hNpH+LUOL92Y(W=SWvhP=K zvgV%kl68%r#I0Qko=&B~i#d`Do+aBNx(5Q(>#d}KqrAB6%41na7Z`R|>|hnYTuIVK zBgF|kZxR_kh08heH|d3wx|DOQWDkpv%x(h-B9%YWN%v;_#VFY3;5-qKPqG-G)WvZz zBj)ARRKqR0cb%Gdt@Vapn99E!V$VginEVtyoB*(?Az+_uZi1PKBwO}6XF1=uhOM&P zZKi5~OjL+o5BsHgm)FWQ%iiVj+i*3_FWC$rVh#nfJjoAsU!F*t?hopDyV`#~+KOC9 z40GL(289RV{^_e5oglPj&gFcQG7#*QqZYILuNkirg)^U^MD*MSb1UjAa1&aPX2n$A z!Kmg7q^6Boo$2H|Iu=o1$ns8y_?ggORmrba!Pj+Ip zM|Gh#W##4Z9J=bCUUD+}ke~F~ZXK07B=_a$Y(P-S#8KSb_qrNk6mg(zpp;cuH?ZUO zLEWwQ#ReBPw9i>MGrsT%8?ixCox_+T(rC{89f4*&+3p&Y(7n9+^*4n5_ zEw0(JCvb6Jv`(E6UA*I^b!JI{$U1qC;sZ2DhQgRpbZ2aWIn(UKXs84DBJN`Gxb#r#(9*t0$c(#dA z;p;mraP#qO-ym#xGQxNnYnOP}!0WG;UNT$ynmuPfQN-by&Ovob_|2QG#0yv*KeX~y z=RB*Y`f?LbB{*klA*(l|p4c*e5${r?n?ka9%A-%hzf{oG_<5y`szy@|Vn_VOj-cM% z)PQu{zI3%<5}zFU4wN1n%h#rG)o9+h;K$dYYSzyTtq)#lO!0d7J?0zNWFgzUB^G~Z zBOV_Qac}-lp2(;5n#(vYl_8F=dE2;)H6iL$fu?VhB)-%uFx;;dmw*YJviCKb@NsIl zwtFQRlHq;Pz%}At|C86{WyO9p3U8XS+1+egC|}SXJ|qD^{N?uyar+BZeWC zp}zwk&!mZp;-ce1aRzKDer9|(Z#{e4?$-$M5|h>JwY${sWX=3NNRe+II9=-Vr*tjo zzdI9fp5=vokJAFrEJTEe&ag5+@c)=I|Ji={M+M9vs|b-R@@_}Of0{??rv*@WGFbgc z7CDSZ`9XaCstpfdZnSns?^HTPKemszrC|KBxWJ>Zl_1KhVWYt>^Ht<{#G^86%GbuO zXe1;1tjY0}>pZixhmh7EEp{ZG_MHyH-5 z<1aoTCRb%9er#s-yp1K~zU#QFG|h7KAiaA+e3Z-1TjP<|U`+S9FlKU^NNJYS)My>< z4^#x66emb%d;mgZ_c@@6av1c!LVV;-iap%d%t)%WuVmnmnxZ zJ&`KjbiEAMhbY)u-`fST%6Gi=tw7`hMMC5*ZMHVtk}zlAQ#=Uzh+bUr2HaYvARO!#m3nqj~CwgI%IW(Jn`5H>$k zv*QTI__&%9pCk9!;6$EK+zsig0Icf9hYA8T&X}Lp%OGGw@H6%D`_R2T!Rq)bX0W2KfQhrhotDvODC?>gs?^Ig{r0@-0_@BNhf zUiZ4!s&qV8bmM%Qt$yNjKp$yRUk0H_MlH`jm5DyuoXCu*IDh8(ud=5UnQC{6mSJxY z;^O>reMLEVBW|}uhwlk;yT4BtG8c#6EZbjTJ|`7Lp4tQ`;@WFyr@ z=J%5{3@=9}%Lorf%LK%i7=P;<%5+Mvk77+)_;{N0Le$|3$*4G0@PMs~Kn(mAw>*lK zWY&0)MiNcUyr7xN{wTCOa0PRGZaX5sha7TkiNAk!@_PfOc*ZgsJT!k0lfAD7f%j8E zSfrLfFr-g29zST`vsuJfTOyh<7rv9;)Pd#f&pA1599ou~kIS`zrdP&c7!bEwDpM6U zq~EbKjeV@nH3S1l^3)E~sJCGEdW@AmXMQJkOb)?iTf1 z>V%JxLg)Kxa#Z|zJA9&}JGM0gS8Zr5+#!ittQYY;&F>J@?>lQFaFFNVIi0ob4vYj@PjH z?Cp1tx{ZtC%Pf9;>{CBJcBUW$5p57t8Ev3BFA*G8I&ovBTTO1y`03ykzBJeK3c#t> zyJ`~v zkMi5V&j9CbFOJyd=bJTu%+>#9prOvzxJN-mSQeY#zG_*P37zy)#}|#i+;Mzm^JZh4 zC%zah{Es8?pGNBc^y^trb=-vlYQblcK=Tlb&ME8&FOM))*>7sGSgi6 zcxceZyoc_!nQJ_zAu12P`<5-A`1W;TPkYS}b6MAdu%)5CxDy@8sWAVkvFDp!`X|}_ z&$;$5X2^bHvlpUj(=ad?K3+nz|D*6Uf&at1qU1GQ*UvpQ%i1@L1fM=be zz|I5yF<+J{HK>~AD#?~`dl5W&;t6(mDSY!TDFtzoQOC1yD0Y9LReWi+6JJ2I@XbvAQkyJJ%khq0_^C7fv1j2cBhx2@ zse1kZd$d?@Yx%fsYt1DN&qqJuUQK{EXO+4yo!y0V(mMX2QopCe^d0{T(eeEQ(a~XI zq{=>*{ot6D8{Q66Fg&UX547W^0dUJ^gU|X&X`6@olaleDk5RWkYxbk$sp1Dg& zxi=|oLs>VbVd=-eA?$Dw4!KmbV3J18K&v|d4Z~1chpa|t^e#mE9b}xD|bX-IRHOBUYfmIJz(WQ}dE8j-d4Y26Q(>Fq8R*K~bc|0hCO^UmoL|2Hos> z&{tg2qZXPD3_jnJ4xF}`{2DAT0ckjVu{S%^QUQWLT_oAM%p}X=i3|?x4E@0E-498A z1!C55rw{R-ns}Rl^#I{m?mHy;HBw1IS%Z*uA0PoeIZP!LO)M{^9*H7HEYiV4sI|Bo zh;8^^WGwas5~NSwCO#WSGP!ous1TmJ+naS`4`m3Mx723Q)gMt=tbaQB^xdUtF|p}H zq8)*QWHXvs6X2;Ld>q&NIj^*{TRzI-14HiN6U)>V)BG*=j0ksmrn(vAQ8IEYM8mNW zsUPZ5eQo}j9H;!(`F!^$YT{p~5phpC+?}(hJ#Tg9Vu&}&+9&icV0xfll1yWzQOPP{ zb=ui^8rFw$uAcH`svGcB0Uq-fSeZPE$Nb&9jV-mnWB&EokV1#oP;ne%+`A4Scwls3 zCfX@tBO&t}y=beVbkU-7*tqQtTqKokB6S_YG_>RvjfTAO49we-w|HnhQ&h{x>WVb{ zXyN8V{*o{Fp7T$e8$%CB9}t-WR3Dhv)p_-s#3ULp((~{5H>~vaCcX&e4i=3%hp;|K zN}RNJV-coW0MaHt)6w%^%h(GFB9&93Zhmp&3=6uP{X&^DPMP~%*Cav~WD)GlD8LF0 z%MvJ7StjDa;+&;Q4zY7;mwRmg&_w|i2u&CiQzY$YEbrOe3me01F;iL?klU=?G2NNI zL}cptSeeL-;}8uPza7KtOT*hoyS)q;o`U+Ci0F<-KKwjyShr!?XxXOq)Eo8DcPSyL z$r`v1TF9%{)g4#gvilA31w-Puj~JJ3x~WKDX(j5yfbTMbB#Ag`H~sv>pd8VR^e%j- z8q>ULx^LWOiyjQWd-I^{M#RCX&9xtv1MH4&ur<%98Vo+r8Z#XjMLk8zb(Wk0bc)O< zG6QH0v>3=h(LlB7AhMwbgS{`GuNhO>{H6+G_r|V@s{qXg3|u6-`q9Gj&u=Zif+OxB zX(dkD5ZC>8zLC@P1N&j=2Y?X>@kZh_z=xof30jd6z}5)rxRzSDTEYsk954(cvwYT4S6rBmC^rU z;y%s!@!t?&6)HdU{#)wKv&Ay2`u;6%m#E;!;4w@X$|ys5eKtmtIV3L@=WKwYay!Z$EE;tA zaMCh47F@p5C~!3NM(sENS@1>23nBtF*DBf;JPu*1)$;~9kv)>4rcEIx1tcov?LX^2e>-U7`lPWOR@?{^ACTCXWF@_TubGXL8=fYx*G zNI(m>plu65{fP@^uIwvQ$5fvaW&Kv7#A?Z#HhkGB<6@;+qW?0tC#(>tFJIDJ2hKueJMqasFOoCe{|AW zB=$?3*G)a2*t^CLSZ8qof9&1hU#e+{wM5Ws4CCJ>mt4cy?W}vyHIk4#^@hj9ge}y$ z&{{}sfGS*9>Jq`Y-1}_HC@w}fyKRY~K_qFhO@UA_w|`7jQ0n6Rf-H4GhR5M(yGgM8lErCkQeRb>+z_iRK_2*ryqd=A$r{Md(|BOwf6b@`L=f zC$7@ypkm!oe!}Y6o(z876`*oviYslEB7PUl9-7m;jB}H^^zfqcID#H=sb{wzu?KkA zA$xUDy+B`L!*RuXW4SN{=Mfj{mi|)qP-cnF20@Yj$Z4hVkZWl?ClFL$5q5!yt?+Ms zR#~O^W)l1Km~}*aPcz@iOaTv`Bb(`d^p8w+J$>bhVQNeiQ!n{iPRvbww~rK*Y*?2n z;|D&EdGQpBP3FKkmcdZ`i3z4N=;eM?73Ovs9^O?*` zr)4SE@iHSRfJd#3#kK&OxNyO0tx=+{ya^@H)y3(@IH1;XYfo^T8sd~C3O<551k}cM z|0_o&c@;OV!tEP8ukD}IW|mty&| z)n;~hbs6Am`NS5QBOs7s&w&bn09W6-bBklmwM?s*3I|qCk}YuLv-(If5iFGe`e#L~ z($;mN#wQ^lC*1uREVo=AK%|I<-}{Ei zpFkQ4-KfK+nG$tJiet;HLyB$G*HZoAW(}x39*8i3O5ut%1PeToxtnECRtJ_6d7x1M zR)Yx;yORyI!7sQ#sOJmTZAz&3nrdaAJIy!2ExtbHof2qPDgxXGd}0}h5lg{RkdGwW z=L)O)mL0-RW$({SoEbJPd*%`)*q7{8{5o~?-Lr2TVIY|Lm%12I##QoL3$X&Yyjv}a z#G0x<6%keDXwQ1CB2kj#!?{oC!P6(GKAF5;vSJ56%6xp|29F%M#l6qp?KkQ>x^x%I2%0=Z)fL-qmQROH zGF{>Q%=2E7&SE7rtj;C~ReI7@e*ZZuPW6z7P~><0Ij+ z7hjw?#F5(=VgJ@qcOZf0!2ScgybO%4=w17%*rdK;bBg^`T9w&Tiv3hnqt%=GV=wF{ z5dCOw*XH7Xi;3!Zdeh-E&I`$pX0jZ8)XSUs{Hny#A-^rRlM!kRwhR^s^jFNvZccPR z-Drm~ujRJwf3={#dU6^NG|_aJ?TO*qZf3)zO|q*22pP}83I#>8A98^9X9BVcwIW~$ znz&JGB3MzZ5<#6D6b0cKb`DuK&UclQ89E3ABi1zu8?f2dTG(uF+Uy*(EX$(#3g!!I z3|tCInb-7}(ZO;%n_%J8`Q^j3}!A)o9Af=a9t2rQ!ui)N*)0Au0}G(QNf zoy)TR1fc&VTg+GVhqF~EL&-IM~rC5*i8-@@D(%Q!+a z`9=3{hz;-&6!^U@Pql#ga#ZFcp{NZHY}@;Ecf94mB?LBR3W6)DZA^l?Cs5lbvT8k* zHhDeiHvukF_GHWUzf(}a^5F>B5NsK*CD25kU`HN&U>&tbRsLL^X4p}okZ#9!wqi(3 zn(nhNu^Dqvu7BW}mi6U@yGiY1&IKe^gC%u=0mi;*YR+dE8+XK<)MbL`?amUev00H* zLUcyI5G~G!RGP0|T}nT5>tZF|_;g53vg|Z9M*Riu$XPrWHuPM(l(=2I@OxEy8GE(V zsW)|7Z$38A?-4b&6|1)EP#y_F?en}qaMUnV_ww>c88o;zZ1Ch!wECHxrpx(~Z}#Z< zEY>x;?RAb`OesgwTj=|xN^kn~gZ)Xp&1w;-mInDZJu+a{Ia9XELthtadzGe`7&LY+ zYIGRUN6gYQb+vi;Vqr0$qITunngYie*S+CQ7wy8I*pC(wOWplO0Wxr4c4n!H-+D4m z`JrJARkV;s=tX;J_ruR(h*|t&Ie41j+#acB26JY+n{}24H9v{J80URMrF&#?D{$oV zVY0izhcOQqJlMWptSOI__({uH7MqImOdbiK7vl)?5#!ivYqqC+J~ah+))LNn z(Gq31vE?VAQATjwM>0M0_OcS=^ON%B^m~+N5=xB-n|{JLdgActYY9(eS(J8T;{Z-- zQ+@iXN&~mmjsnHi;vFFm+Bb`B&K z=-M9ye#^+e{lxM{>A$pX%;07Dlvz}nSycYmfwn1=wJFOz_OSa?Jy9fKmC~L8GN#>bG z+fJkHcbMY(rr!PB^A`>t((2YO4tL%|$>mt;zE0LP9Yr$d_Ma@QK%*087eCz*{G`;p z9R*M_vyGQa-1^ZZm1R+$3lW{^w31^Fj=0hb%}i_YE`iGloL}e2p^mAE<5{CEBZF@t^;2F2(%5R z;6x%gRK{e$;%6-%8?bBrXIegMhKFuHYZZ24)Q}e{|A)tTKoB=xJ`ER&u8351yU3f< zzAL4-(42NDoQwl?mN5He7=)1@GTll+mE5VV#zU+>iXhHA3AYZ%;D1B-xaxw9>L={T zM&9jBAEckYC}>mJ6?v?;CUQEI*bW;&iz#q}L6X-=c{?evty@ywYobNf%q?J30UsPI zV*(0QI3f+)XvhR0zmkgs*?G+l`vSK!547ejdN~gcY|U2>udw*TrC^$<$5!xbZIJCe2E!Wy?qG#}pYQ zo(;ty-jX)!ehly5L87(}-A9BHg5Biscklk%6obuBXhIFZBwwKIOO~1OW#%T%F4JVx zNFbBjqa-VIOu9-1ZsoGM(fSr;vmKtvOn8!m69rrJ-la$rN3bdag&|P9SYx>G-zj16 zu@9H`iKGA-7nOZ+jrW$$&~LI^#?I_{awR<1t7=Bl%kp6$0t*=f&Gqgp%wb!R`5ID} z&y9(4>^bqsUg+TLoel5ftzZ(;<@;s+N-jMU|LT~{RAIf6UcrFcl8dPK<*HvbuAbON z`Q}>D$ax(tZ?=9zs?cxbYkC8f6uiu@cdK8UcxP`JgBhu*n2&j33LUOyR}?6hd>baq zEyY_vEByrFoCSO>rk2Unqh2Z%I^L&>+KeSsWBh%3IbF*?(NjD;w_-TMUnz@Abv7E8 z_gyLdIP0<|5y{n(cg@K$$l+Z)FprK18y0+_cRKp7eZpA+`iec)XCVB_Xm&4C|4O4E2!E2DQUIQ4Q)#qHN&+4cy`cSQ zR(rPJ;e$?!RTH-&Se$-nIpgRpj8GQIbl z^5UK^7g_qH5%A)pS@rBj%@BoLMdqLPGVE>eU?70UHCILLYz&iC9}bgwx~;K~-U_UZ zWCrgKTWw)j(h51BlKuSCI`n$gy$amN@kP}98~AZ4Xj`KRKuxcLwHNmpbWX=Yc*V+m zt@NagR|o^owz@fJHrDclvb4_$_f#@>xNZvhAXw^#HsNpB-w+9B6JSR5KnnjBJU|?! zC{wgh#!43H@&si25Jio`PmwJAtcWiwl2^uNhAR_-&u1DlJ0wT2q$an{IcmsTxXb>) zcK13o1GEUR`)lZroB?pYPZ|#kXZTlyL-b1NF?~XFWwWSG%es=(rO}BPvZ zR}OtgPV~o;KaBXKIE(5gS{9jZ4EyIvb2qJc51RAm$K*RV&geP|T@IdRVVZ;vlFgHe zHrJuOuCe1S>sF^G4sFav3codMnA96Fot_O3^V5_SL##0zo-N~R7ogW44zLm_R#3DN z?x8au_e}0Qc_QJQX!y$mViTcnr^lL~b}Ojf$(yN)84&q=$;Ytq?8nDcSC&riH~4Z< zlR!&}liOXd}cb{hP<=jPQDN=lvr^HkF95v%&cS!9W8M@^~ai91T z;1-XrJW-m+l?aZ%_-ps6TczX(%;hpnCCA;X^Y=055cpZi_Q{!`2^vmgba((+>jMCt$p?ztb5X|tzT0!h(yUwSYya8MDztob(9>?~b;LV|r%hDSRO!rR zGWfFCGz)jDsZ2Izcf*3h62*|Op_g-LoChw+JVc}^b3LIJncUF#ihy%QRd3h?MO)3Z zXPsg7ib#5{v%(Y><{8GO#v+IKC~NfE^fY~E$mJW^9ls$g)PIp)&vl65r|Q+EwG=H8GBt2N3HFhg>VJ5`0TmXCz*`yr6#mdher0oR3WOsGl_*zG&-p$Yj`4 z3$GTFrJP_mkC*pii)Uf1Idk3NsozV1Z$7s6-Ue*r`@t-Jj$F9~b3;A1p_%@`(f<(i zd(V(Oljg2E(H9zG>|tll^=|Vxp@A#+?I4|y&#ejSmueD+WA_=n-)mE3ArV$K`*BY1 zjkY<|#hRlHRz*DvlWZXt4Y1ZwypTKcn3}i!$f(%I8)p`70TL_6F&Wa1`>y7KPP*vu zd^cdYS>&1Uf9EP>15agmwjZkw!f1vcH&GprZr+GLjD@yN~lM`82`m z6hW!*T)yriKIh3L0k|mB7eN%;@|eOdp;!*z=7syFO7aq7EQjV==NZAhqrA`Y9X4AXWXgw^yrh?T25~fdueV}|JyOQgVb|X zuRZGJ?<Y;LSUJ@-r8;^sgdJ#Wb6NAH08oqgrR zHbx%|k^h*&7TjjV{yO*89`LoDp3*s*DHEt8`w(>&VAIKAnID{Z>t1>_s>mO|_L&bd z1zU47a7_RNGi2`uZWW{v?5J@av!n+0@=5f=*U_rMA%!D#)iNXq^oKHo(T;)0KoT3< zs>taYl4_ak36vxhA`3odH3)HjaMJQ-k`_jL)ZQVlKZ43Hfsxyc*~O!i z@ZQfZ?|)UpNe^W~FmFCjLT|;VSXP?_Yyn6mp7T`KkU?27vwb|YWUi)473F`^=@vFEW%1?2D;v7BFzJilRcJ(Uw{&Fqj z?`;mDm{}n(d2crpb=UzI^d1G@Yy@a_x1l^+_5yo()C#tc4Nc}_pn+OI4e}f)+l*4?INBNg-6qa0X4}6vaT>}q z-Z`Y@)5-atq(6K9AwBwk<8{6!vhldIIr)yw2xfncSzkxjvUKBO(>*@UKIaJ@P8sL0 z0^=#QlMsUag}kjT(P4}$JjR!>ll*=rShH(cK2$8pDupWg%HaetarRgSxfTodi!H*~ zU*tX*&1v+{9lM*?rg-+E(RsN@dt%S~0|P6HHxL8yg}6hmq;goWAa3|tXF^fjUIxy` zq9I?xgP+VDWckWz%EB_oSN9)^Mc!iYj4yy_iH+k#jTG1Ps+FCv;Hv|w8_p-6hCi~& zmPLv24A6V_B4~R$0p(`!SVfmUobb_2!eldQ!;Sp*r1zQtjz_bC_uvFAVuq499sF9C zI8>I48WdqtaXuojzMd>hNd0E*?H8jBMw@=8fa<~eDqL>W4uEIdL^^ndjL;v_gzIpf zB7$bM(nAk2Dn6mF7Gu_E+IeAI%={BHh#FPip$NtU5P?t) z!tm)VAmoSMv?$fheXLojEO|<{>jL9Z(-C`b1@LMReDshY?2TJF3HdsJ$&m_esy^1~ zC|yqtxoM(QdQ=iaR1(H*ZG#DE1v9b5@wCQUVcnjAx^EgFA@Sy5f6&IewRtbS%lK~k zM*GHwQTO*X6j24}+k?xfRfid&@;mX+?{ddjq4*-16!l^uxFfUI8#U!F2uwmwg8B7> ze3Sqm!yN(_uL7FPlmV#G9n`A7a}DFNK*z;17y)cnrdw5hVag;9aZJMp?Z`*ABe020ZXHTp$LFmvwKY@c(= zv|(fMBgizu;fdZ#J&R4dw?gSeo2jTFTT z-=CD8(aLGpKxMHc%6M&!`u6Zq$F{RHun;%ofV}~T3kEzez5GcYtodV>`Z}74oK3=P z#ieX`2`q1Kox+m2Wk9g`9{oEn4@`}p?{arY>pk2bgPI|=qF6eZIA3x9s1rNQ&P57V z2!Mdm`$Qwy2eTz8NZ|lV3T~$WC9^Ik_Y)Ux&Q}#Wif1!MTzMtk3?FCk1tk|Sl_*pO zf(UR&RI&E0Gyu6N`%$W?}r03y|kDjKeiV?a0R29=3W`ibf+5y}*VGx!+9cwKOx}zzoZp7EElK(ff2hb! zZeR!`njK)uXQ-`sR%~u0_1jEudsgb>?~JOi)WE1JPa)kz0Tny@Pf$%IMCqi~7C6t= z{^7)enO8D$8TP}JdRYo(b>uz}`?#Q5>Sc~G3nC4FNY=hlpEzJ;d~V?$>AKO-&Y`P) znXfMC#O!55z6L{2_D`6SwiRf+Tlf8D0 zQo7<~5@VymhJwjx&R#ac=ahiJOco+XQda8BQzxq>7y@kA+f{I+>$5CjQFV^xKKt$K zX*Bt!^w-p%p;p8U!T$6NjP2U_w79dHuWlCTzD^96t`FqNs>t{m$>ggdD}w1xJ3MYl z#A)Tn8PXZ#-W3-ZSL5fT+JAtBzw%7+QI;=$JhB^sp{{3&@*{`Fc&2=ApXapaF>&J9 z`-|)8&`Vd@)d{9!j{3U3R;GC9#aed<>-teyOJUJNlfT;Ar#SB)DH(6?${|}O>CGGE zM*03K!e_RWiM%<4udB~i&y|>PIr-YwRr;_Z)7&qvg5_WGj~Lt~S}ORtzLT<7d~ozL z&|m+ha4)-q{!w#@U;7dz5{tYN=6*@8&FW(kieD`+6QcP(Ks~VCLK8vsfYtn~FdMXJ zaJTZZF3Yx)nVrV50a4x;aG~WMs%pv?1at`RLX8Of()pK z+muUifk48t=u2ReO3Jh8(19QqS+-EII|m1!&p&@?SM7}^?Za%iY`FBECNV@zc*QiZ z8%2c_Us>x9z0}Tx`H>X?-CP`{v8RI4*z*`P_87elDc%@AnO*HZ_T&oXAuE!S<5L4_ zza{pRMerM3*9MksnYtT-NtfP-gx*+P7ZKopQ9b1)WV5cSZ}d25x5-CO?|uqWOQrX- zIE%l$m6uOH>;{0mYV@bOP&10`+Tn5j*h263K_HxrTHuE{8LGXdqW;g29*TR z;$u-7BNWN5IxsQ$8$$Cp#8UM4EaJ=t{eg?Nu9S>P9$W$pN6~&4Ov2SwItK za!ofGEJ)lS+5(N>5Fr6RB7N)+QC>$cd(2~#ov6DMmmvH$Umt38FMfE38d{M8tRN)b z$JOQ!LQvlJj!g6Rmj8yK^*o4~IKJf>ZBb&c8%@{`CPa@>bHbFlM|S(pSq>d>g`N zmT-YI2t$4gt|c&&gYI_R-;V-Xy4pY+@(FHBL=O}2(v|!y3gjH2R($+c12dWdkYpOR+DXO(g#T2;w z8oqe+ix4fvGDy`}>y+9pJLS31fCtwK<6HE<{Y z5*waFyzVSTA97KpT9^Iw=gzb3H@8*GtMer+vwfT%j`J0~QZVMcLa`s6pl$_qmVc1T z|4wEiFxL!8w%|Fy5!_*OYc>qS>-gi9W zVhBumLVIQ&zhwo`4CBqq4Owox9$DWTbL84MyOuI0w)_=cR%`gzaqjQ1fU8HjWkNpy zj`8aO1f2xlQ!=MyIsvQ8CUzj}5_`+|)VS36ez^zY76%V3E5!)!K`h{>Tig)w}dS$ zX#h5QMTCG4z^n+H|KQJbiZxT5 zasJzKiK`yi;6CZHgp4~?VE!E$f*z@G>gMR*Nerv{^y$J~s2<0`_RjW@<8g1PwsrM= zm}9fYH~W$b`dHiQkwTLOFKJJ#ZbsR1$G!8B-Mhv8EYFD1yi9`ls_X=Qz`b&1&YESS zL`i4xdCmgY-F@FU@drqy3YPZ(zm_$s0@=H&mmj3@ieGPlt1|z%*5~Jsqd6-d~b7PrYk^-yJ!Od@~qG~1Kjb(*xnJIo# z&I9Em*E6#8{G{IbalFbH%4;MuTP~kK6IIsdw*;1iJK*~4r`>!9j#S_K^^0=E3Y9g( z5S!1QC+^%qQsy6A`PA+H!F zxrDA(mx;7{(j>^GS-doC+?!8loo*5M{Y3f|wJD9@Gt0AAPjP`ro!yjB17A-^kq)H5 zL9a^ol!9FiO@5&RZ5-K^4B6+TK~ztv#|?I;3z{1d)^?WIr&lcO9qd3kK<6xhsIGNwM8DxpRA9tcq;e65S4??t@N{Ao6@#md;N7!ZixxO>cA_RY1V_ofZtGVr##Kp7=b^R^s8GB$ z4>|N<$j`aRsFxYzg|ENx>d0gQq+W%B)PLba>VV`gX&0$CWb{5u%QS=84^7yD+yP4i zC_G`^vEA{WAZ86KUxq0GsDh&58Wa#>?rFBk%M)KJ*ZT9u@_+w+R=L2vWlt;k+Ro}4 z5twryZG1ILSzdnQ4TRffwH&%zofLhOYQenMw^f=+odk$?9?-5Zb?0Xw)Hb{;Of7D7 zf6B1owT@@>Q0z%7kKiAyl>D{kEfDzte1>t0s|WB_%r)+WjGf}Yon7>mntKz#HZ)|e zXl^nLpi+=7(|R{~`F}5txveAkw_k^u=m$BAIxH&Z+GY?g+2Yqw_KkzU_5`h6X1`3#zR>`SAWPJkN{` z9}NLY0D(1{y;3e$j`2J_bq5{6W;#vxLXZ9_T99pmbv=-Wx)G*zy?Y0ber9t)qya_o z=98+AvrvA0shX<(s6;9`ot-r5S;k2$91B`>Q5k%1!Rjt=<@`B$?ac**5PhW1N#v`G zwz9QiGT(rP80N_+VcOsG)JHSc$Pl(}caxF?gQ|sh1-sPvKNvVqez!MCi=tjnYwVv6c}uWT)8}k|AudAb9`-*5!>D; zE<;Ri4_^+fV1}1bid$@jMx=Elp=R}Zu7Q2Q5DPy}e)nT;vme6k02oIGbiL2!1rZO zIDDY1f8}=m2D%TIF95esx}B6l$0_v3*0^kcH>FM~HFNE+K;qvJuaT~jOnczYNw7-P zuU=UUGg_K3m9*hPU}4X{WrW_b>C)K_aTTZ@WrUVgWkRUU{KC)?ncRS4MLti)(9ly1&cyGE|$ zqu}5O8dCkqL%}|WZFRY&h6?fNm2nsHzho1*(G$4#v4T5&C)NTUJiN-`uQRgvn9WDL z=`dgC@034=%fxXI;b~)bo`ZW0#Fwdzu6r7{+vKSq=#x$9H88BtpP7<$jaHzWtiAk% zfrVog`wM^xB8O#g0BNjIFRW`19+MWK%2C^!XAwCT9d@6A;k|+?c-s1A6Sv17t%4U zCE{7~=~q|G`T|c6Mq273JX89LnioP;;HTZt<|MBCbBwMSLHqPnMpxyPzd)e-eX|mQ zgXN7Wg9uJdW#5A7SwrrgCERqM4QZsyD$c~MD(_@gb-vWZn@35#cS^l}PI-!hyW|=otJDG=X3(3^kD_TSToBTg0VnCol!UO@5sCXRK;(jfuF7j zxvSG=Ez|5p6F%9HZDc(grcQrA>ZmwVk1mD`;-FV@;QIP;@N31ZY=$X;C3eZgXJ5~Y z26A-y2fRAJ6gPYWuh^N3rY7o^=fj8&4jV>`^VZR!u{M{A?qkxp)=dzENkR!}3LTrV@~e4-aY+)3R}ix^&YRPk5gdGu$4 zf%wY@)j0Y{j7JX&@tmeN<$hu-(WqNJbgETK+5|f%;UAuCcaAuug$OcUJ!a8 zcQp$&7D+L^9#P9ooCVRSY4)?am9e4be5d+LlP|bmsh=wu2QCb2AgX)jeWdk*Dao)i zSI6E^3Mu70akuEegk5>9+5)!$$FWc$s?IP%%c19i{X^UDK2^=m<=mitrN9^5PscO|sZ36r_D8P@5_MR+$vfu2 z3KGCX#6S><6gC8uYi6}$zX5@IHM@)avR3pa=qjZQxhRl$-;=_AAKUlSHrI;A?{G zDx4v=HDn(#ZhlD$`!7cT>fq0R0DC*W)IfNM^;MEO*o+U`>ZJ3wJ8>8PA8!j)_hW7U zBHx8+huu*ix4oAwZ*@UcnNJAc^I|m3t(wu-unFs(YC-4B98__#pnz;OD@IVT(b(n^iSy&pA|CT`WwDI7SBmt}hhvc=J3 z_+biYPUhIw?eQE&Eet_UD_Z~2I7Hc-psbbs>C7lxcq|ckZ-wVEm?7dm_)?U33 zc*W_QB?+T|th~s=;)G@ZjG~-fvZOTprm^xArVH_5QOzuFpX-*-tjd^{|?fNzepXfLPn&F;c}rkIW=+kHSge%CJU+0V$d z6~!*DU`yAS`RCs%azB$o@mYHS;jG=hRT~Y(T#U0wat)Yxz!eVq(-jUBI-eseaBCn* ztzCf#+AX^W%9svJp+*gX_`_3Z>yjYoYX!p`M+*HoPsR;?P@TaJmeSlu(0OpvJha%B#RGM1bU#ndhb9Uv_kO9Ul z;>+X@!~409NzL3HYKcQ3pW*{~mHzLL!asOk{wJE^|ErF{aubqE$5Lf3&XgT#A_VvGp;YF=bId0ms%X;&-#<<6X4e}<*4 z(6>)o6CAW5)k>@4oJR~FM~vGXp7&-i?Ot1yYihsX*_^9BQCocFwDZ;tKQT=PL<^0q zyxn$lWeD9(qBndoFJ|$J;OjiuH_>CNRac&j9I1Wi%E_CJ82j`;D}qGpf?^~eGj%hy zw8MU5>kg9_mQ4t!;zE3I?hjZwfD%8k_t+CNYVhoN5WBz;2t=;NAOUuC^W))z)3eIA zl;Kb9Hzo>&^?t|?9pS07mhnsLaDY-RvTm!yWZNu8f@Y=8`P-d5us1i?P-?yg8wnuwp`+kosqNXp*;;B{ih6xR|5a-fUnB@# zCVvlJDw2#@#?urZ;jnLKaQL74pLoTkj6PiNzel10*cr%y$*OMarWVNGxwMZKvL3Ip zrpLo&rr<-An*uNyfLLG#3+<9l{`RxTM6S-g*07xv-Ops``#SR9_dVf@c%^d#;nwtC zP^TW?B8HFh*iBjVG^?*sG@?z(7k)zkf}b%6ktiD@-&cSU?$)+GsA|8AGub`G_GS2A zJ6qa1ldi}~A zsopBj^kE2a1plqG!~crcb)!8&@EXS3440n=zPuAysw|6wP-Xj`K1R3BC9iOK&lgEh zO*^1RzsYo;ihEGM*Z_j0uM!;3rgW#s45>WWx1g?hpbGqw!|>GORTVFE4NRp%wX+O* zua!@2LGKbuQns5cs`X@^Y^&njmkiQ(16>3^KG~_w{ zPCNJpbKNo1_pe!I&?iQWsGeC%*mx#2#dtz=pZ!R%=Qoelw$icB>w<^OBpS;HETYfF zcuBmMe>K<1b-$1s@rj#-b9|tMCWb6bU~4DKc9&dA7*_sBM5*Q@AH;qOUs#M_*K7#> zY?4S<@I{BK&-S1c_Y*1gq*>&i5{=|$)Z5g14K8fQZ+Wg>v^3d&?O7PB?_*ObL`oUe zwxBBDB(n}?u9eEPRqwZ7Je9$BH1aFUnX4>csb$L3Q~YvIED0^hObh|2sIHmug>K8-@3Llq?r@>#g85p`G< zd`%VlATdbJ54O5jqf&sn&=r-he6o7{-55_|hk8BIiXGEZFDMntDjGB z9EqQJWt;Hyl``2At0ejD+3<-2?nhRs{dyQ8&3s)A2? z#+YZej-$=`if;$g@3~qddzF$?sVN%~+T@Gtwy>d38rs_4`h6=v@sk zo_%*~7Psd=h2Qx6HeX6ytY!!`eiN>02geVC_1uP#;=uCY@!(RZ1Vr47+NmI?A1=L= zt1A9FHNTdFX4-}T%J3(^UHz(arP_eK5hz04#^UD=r=K;(4UrM7t{~}%RaV`8keb5`s5GJN2 z4E}GtZIJ6N`{HJpL6~%u^4Hm`DxK z^TIx*pWt9NxtbB4M6z0jl z``>5faej6|EP|iy*Ae!&)FQQD8#MX|HJ`L0P2z}Yb9w4pdxrphi3V&us&1L~W86G> zkL@0&-Qc&+0T=NccMDX8t_*KB0HXW^!HFcc|HnJ0&GCSj!W7FJ?sPV0F4$ZrCR%k} z!iac75Sos_?W_!sE0VdGajW&i@RkiB;|5>fhTVeuhN8U97+;}_btqM^qZTY8x-r>J zUiDu=jqB+c)eyDQ@OF7Av@5}ek(jWVTM25n+s)eH0kEpZ*}{rsi=7a1{>3P#t8&Iq zy!JWdWrJL|CLJn+US^#iv*b$jy%!=@KFb>C;Z7^HvBt#2gzlOL>`JjhjjygTW9*sn zM8A3UFZrL%;(K#zL-i)99sEu;J_2{7yBnHrR|Aj%*n6wG)X~s{QwIu2oCWDm%`abc zJm1Wqrm1Vqw<4vVmj*=L*G-+kR%`@lUl>(Mcn1flCAUPY20bq!joV#OtM%%zKtVl7 z7`)1oW8E3YBEr_CxE10A@|3y_>bP*J;C`16byMcUO>jz0h~NnQ{b#4I8U}FvDRdEB z^A*HSiS7{j!j28>1#$)wC8ak+;&avkSqIF(<$g+_U!`Hhv<~jl_vSa-u?&-nX2Xmk zpAvVX)JLFZkjjtEA&L4$%j9kI;Y~2hiA590CDFFD?)0e6Sf#qgffWbdTPZyMUwdaB z4rSZ-aZ;8_lI+V=lyIZ$`=lC6wi_{sD6*!p)=*qZ!XVj_JtEm7`!;2d5EYSa5ZOn> zm>JV^j=Q?+exBobj`w)qVW3xL6gxMpgFo-v;hrmm-Crv^0lutLMu(!wVmK|_UN71r#@k@bJ!ArLXBJT zIDK@BYUL;Eb~TOB(Ahuy(L6Az$tB@FE5EC!S=?S{aZ z04B0nCq+d=jHdDSNXvQ@OF_)0rljDbO%m*nb(hwzlaAg(u$H&cZw}d2x zqK}y%=!BxtG=Z1^c+8{7MSd-z@rDW%^TE`Tqc;h$o(B7Xw7@;_xw=+0jiCIq{(1l0W*t+n-HCuXoi zfS_0BQ&uln9p7Nt3lx1p&Ixjamv@&F^qG~~Z?alxDs`$YwSxi&V|NY+lE+( zw|I3{$PLVV-kEcJU)E2Tk1l_qHU;Wex7-BYP6RdG{~GDdoX~3QmREMBu>Vo;tmoNJ z>@gPaVyv_#ie8KT6tB&={Pd(UnJXfS-NT}v90V4G01#gtsvevy105UwL=Qpt2#Wht?c#|eIj#HogD{1PVjn8Zty~e%(4(iWGE+Kzp{q(?SoJ|(_-0{>?;$yC8*c=CI zxxfz;1)zw7ShN+8C??}Cc(zuhBOK&v9_QYd2e?>j-ye|kSCaqU_s2PTK_GEmDi_! zYWY(El$V3R?+hkqK(3GlIsJw$BE_^i!V`x7rJM943WsnI^ut#!vl(Rl5%Od~>1@Q0 zc&#Eugqw;O8zYG2(r-3-fMMVfSeeuRWKzJ8dEi_wEX+Vsq^^3*Bpl}mFbaXN4RJ7U zvVz60;R7DNNe`BFBT$nOd@^@La*-!bZ{hHIN&$N_^vkQ%M$#%`g|pN~(mB(^oxe7c z7Cz&GOPPaJ0;eRObU#FPMLYMU@)L4!Sk@Vez%qScW;G(~d|?6IU=z~e zez={3&5Pf~*cE_}bdQHNKfEa{u3t?OZT-TxVDcW?e`>V53Zr}{#Inj^fOqB2j`Z-O zEgsVhLN;^*M~ftv?ZYnNvADEQ)QzLh27H!J-(vF=%+&7^bTSsNviopg61HKV?EVU3 zq3{#V&ovSmFzL~Yt{!LauDR+-reykr4vWg|>17(_-6kFJJY_(Jb~YB-#mPFGs!WI* ze+K-ZWJ8xMtFn9!TQd@$T%h}j=f_{#@4a>Po^)kYTT6XjHX&YXGf~!0vYBH9@h-U- zab>sSvDx%8?gHlywANHPMaWZX4!-O}Z9hS$TcTqeISESfkU znr%%!g=ZObC3i3E!WLWIm{(~#U3>wp0u}gcuypp{3}A7~$a`=0VJ6YLEGXIEsvuu` zw@78mlZQLFRg!6x-eG8RfVb$#I;9d$x#iO7()*W?%zeBTP~-A=lu+BABQ(e0x{s`w1CLO=mx(x3&{cOviqpeuIpgOT$uiq<}yI9 zSoF(VXO6*#ouD9pExZ`C?|i=Kqs=v9VucB8sp9Rsi*lL+ZknWKYs5wb7~Zj)e1qU& zWnywW;0m6}GfmryH=f0jfsat-!$jo6U)%_$m;Mq=1V{h(g9$UrZ^e~w$-IST=B-m! zR!-}C83xy;h9~$y8kfy* z83xk0+P$un9zS@@$L2cYs{>Jyzn$RQ<|8{rtO!!^n}F{oIJ~H+=M6->@mdnHbe6K# zbTf$Oyw7!rFZ|Z1q~c-!OEih=RJRQoH5TDrYv?**uwz&qRq%o{&fDqOF8)=Uee5B8 z#m}uD+zhxawSW5xv{z;+1kPIOqRLKY6^IV?jVXJZSSu|Wke{3z;Jz{)o~l4!SMtX` zVa4Ds?DJvsuQcxw)ZTq#vP+|U@SiJNT)LV7es3Wi&sqAzj%UvEqim%}uA=Fh5kV)> zntb~r8Hy^rTh%E^IeFY4-qx>a`DSQNVY)Ml0$_$ z7GR+O$uT@^b35|fFIr21U{#~>#UM|mBJ%C2t%dxcmF&;YEVxQi?49AucuQXMgdZ*W9(wF1J-l}bf+LpnsU;_g7o+xB||C1YxVB0zV9hi}SzBNme5JKd^ z9DSk+iu(Ud@D%aB1+ovvEX1XPz;#G+wx&Dzx<-%qM4P>k$_T-UN0055K6p8%1ukh< zlqRbu-H=$cz&h@za^Y}>)7&MtmMevw)f@Fzm0UI*F{KsUUsfOPNox=gUxtPF)1!mU z-%UlK#wzMK`lDJWOEJZlEZLv8UF?XXTlBrv=+a2CNZeVEZVU)ho1kb{4tY=yBhOY zRSoH{G*hp^HiQh;XojcG^>X;dP~4n2WcXGkMT*oxQlaQZV0}5QV?{2+ZU4N(Y?%d9 zD!J1&@erNd3>_jYAZ2Atq- z{wj{eTZ(*cC2TcJ{ftc$;j0kGL`~TE`t8CivJY1MR0{@Nixv;+gcK~YKX~S&>&DW_ zIO4|8?)lk`c20n2`OeLWT|U#Vy^k*KV_-VU*w&`i*8Y@%A(_u+jL6ebPskoik++P| zsOcZW5EV~LiQGLIGEK*Pz$$S2*d)19v4z>wpZ1M-LH&VQrF255)41L#!;&N8Bfa|( zhopJ(BRQQQp=jRJy;z-M6ghT<;tO@fW$k|dQV!H?3hO)SKjCq|NQY!yI#4^Ur_J5kK9 z4j1cVjuz;?0tulCq4XT6L$Z}y9!}CQg4V-MF>W@@x{&i4#H;R6w-@j;=vMBRD>uxXqb4^~YV@uD6 z%$|IBDaxRpN(D7dlWoTqBwKl4Yx)J#`C_7-eD~D)BU&^x)ERCy-t#4ng^;z!3gnwD ztl|o@RYwQsQ-sRjKJ>=Sig0p!U$Wfxo}1P^c`xlwO>G2?r!T{;c`xPbgvWJn%i^)> z>`m2mt%K&{w8fkUnA4P#ZBwQdC%A~nh6q8;x{B+rwoMEM<&~tiKh4Ly54@2F09|5dxobwM5zDaTkjOu*#tRD%b9DBIa1sn(zhL;VMSXy)T46_ z6Fjqq&|ttf%F7cAn&dy{_amPTZNy7e#~z5c(cC6;K83SCaUQf#hg%ce9}yy7R9gov zNuZ)z^*+}o9Q4{JAI}|eZC|DUQy0xtsR|NYjuwuq%hnrsON;g+Wp$;H#Pw_2`5ui; zC~!|aRXrj25OiN!)e|GJX$I<({JsF0t7h{E1tsUu!LgA&LUos|SWGluxw{)!{?k~9J-XSAlu}a}C=_R8w_-ml&BCZqO+*0|jK5CNv^qr@0Nbsy!Yf$q=TJ4+5AGyKjhRw}>dZ_uNef?O zujIlfN1BGN}q~EAndp(3v_KuR1Q!%rr*36(PVq#2( zfk9si!U%KjpxpaWaaFiQ@Vp$5XQ{J~>_v=SY3fJl=?C!spr-Ae{p-^Ht5?fAj=xj} zop>}b*X*j@u4meVjTSr${7DhzG-^556D@tBG9l}Jqi2Qjv=)5Lc_bsO#4kd44LRgF ziVsA4Qf`iJL8Ws$sMw8myBO`E=^TNplvj$?5~Xc#2`+K>mSHcpE~mYcPTxTNX2p4FvimMZo~ zAT`bQXeKUu!Q9C9Y90~GEd5<)sT;Yi9P0O_|33lS#lWPNy(i}J2PX%eK!egSgQIP( z``FDe<=2=4Y#D?F`6he*?W^y=!~>&j!NwgR2Z)Lq4+g05PJkN6O7mV@Xpy%K|I+Me zZDu}u)x&&rQ2mEK|)i{fHI{^O^_*ZT^*!r~Vmoec8MEkYAu=2k;@`ot4mZ~7 zbCVitjrs!o9afmOyv-Dl15?3n`8ymepl?tZi;Ogxnqbj+LqCBAjgA*F;eKnt(}7w$ zk^*WLfm#W$Fph?rz~B+B8b0+$_p$@wuoP?yRi8em0Yq3h83SaDbNh`=j+9^4j}=(f zfRhgm%9BB^cMDvOF<@|%mYaV1vIXA%?S9xSFqneSb)dzEwL2K64Q#m49|#yYjtRfC zNeA})GN5MBYro;vUDG@-0*R9D+EUm@f!x6PEqQPjYN!F#2*8^B`65UXRQrmMe6g=6 zYg?d~%8x<>7aVsnR=Z0H%IpAlA}C|bBl82}C`3+30yLMxdW>2z64*T&w}>-*!DP{` z#4+qU83y)iasZHF_S7M3(?(>^Bo2x#D*Tsn2Zb3J{rr);x`kSWHxc_yFa5KYWO(A| z-|fCHl*oxukM$O9R8^ld$KoSq;43pA_=S;OOfm6Xbv;tIFI7unalE}iJ_m9ioX{J& z$4S&fToEiIkrde+#7H->*)A~^6Mr&>B(}(4gwlB5EdsA0_Iwmxhb;WS#HR^5)P>d@Ly@^3I>(osKX2w! zvCr0XqgTR>hL4`^A^W+%q=Vk{SRKz($l`jL89uBae>DmeoYdhfetg9b{Q1zj$thyr z{XMxENqL4#!s|)D=K+sgm^-RGEvg-is0w@<6>Lu(6_EtfEerqv+v|WqfHXlEsMMsu z2Nk+ENRm|kUP2etume{_qz41m2TsU^Pzk;S?`I6*tkKp#97tBO?3ow1ZiqJNXHG9GCwG&@mUe;AUGm!r=bl5+p*?zM zc^}9u#GITvB%gM$k?x+`xxQ*1(UAWiCZPXorl!9?&)?(r?=kTA82Eb({BMi__}77d E0K>dkf&c&j literal 0 HcmV?d00001 diff --git a/thirdparty/cbrewer/change_jet.m b/thirdparty/cbrewer/change_jet.m new file mode 100644 index 0000000..b8d4ecb --- /dev/null +++ b/thirdparty/cbrewer/change_jet.m @@ -0,0 +1,64 @@ +% This script help produce a new 'jet'-like colormap based on other RGB reference colors + +% ------- I WAS ASKED --------------- +% "is there a chance that you could add a diverging map going from blue to green to red as in jet, +% but using the red and blue from your RdBu map and the third darkest green from your RdYlGn map?"" +% +% ANSWER: +% You should construct the new colormap based on the existing RGB values of 'jet' +% but projecting these RGB values on your new RGB basis. +% ----------------------------------- + +% load colormaps +jet=colormap('jet'); +RdBu=cbrewer('div', 'RdBu', 11); +RdYlGn=cbrewer('div', 'RdYlGn', 11); + +% Define the new R, G, B references (p stands for prime) +Rp=RdBu(1,:); +Bp=RdBu(end, :); +Gp=RdYlGn(end-2, :); +RGBp=[Rp;Gp;Bp]; + +% construct the new colormap based on the existing RGB values of jet +% Project the RGB values on your new basis +newjet = jet*RGBp; + +% store data in a strcuture, easier to handle +cmap.jet=jet; +cmap.newjet=newjet; +cnames={'jet', 'newjet'}; + +% plot the RGB values +fh=figure(); +colors={'r', 'g', 'b'}; +for iname=1:length(cnames) + subplot(length(cnames),1,iname) + dat=cmap.(cnames{end-iname+1}); + for icol=1:size(dat,2) + plot(dat(:,icol), 'color', colors{icol}, 'linewidth', 2);hold on; + end % icol + title([' "' cnames{end-iname+1} '" in RGB plot']) +end + +% plot the colormaps +fh=figure(); +for iname=1:length(cnames) + F=cmap.(cnames{iname}); + ncol=length(F); + fg=1./ncol; % geometrical factor + X=fg.*[0 0 1 1]; + Y=0.1.*[1 0 0 1]+(2*iname-1)*0.1; + + for icol=1:ncol + X2=X+fg.*(icol-1); + fill(X2,Y,F(icol, :), 'linestyle', 'none') + hold all + end % icol + text(-0.1, mean(Y), cnames{iname}, 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 10, 'FontName' , 'AvantGarde') + xlim([-0.4, 1]) + axis off + set(gcf, 'color', [1 1 1]) + ylim([0.1 1.05.*max(Y)]); + end % iname + diff --git a/thirdparty/cbrewer/colorbrewer.mat b/thirdparty/cbrewer/colorbrewer.mat new file mode 100644 index 0000000000000000000000000000000000000000..ec59ef47624450f36bf71566afd6fa552df80ce5 GIT binary patch literal 4475 zcmb7`XE>X0|A!-_QPlQBtyGH`6(VMBV%M&!YP7c2UNst{s$$e$r8c!%d(?_em5N!a zM#UaQ%o5|z?>YW2o|n&k+%N9qdU0LX=RCik!3W@uWp6p zF79s4dy0zWq;#R{pLhkfQ@=5|%8#ja638+u8tXo1zuy#-Z$b&v+A?-qXZXhBtU#p{g_bS5`T8^AXfvlq4ei5B`s zBS|eVxViAu2e*de0?ifR_kgHWu5D1KnAp&Osz-%8?Vn1?8up?=7C zgj6$owv}2DYeDCTG-QxH7|Umf8P3w7E}GD3rx3tGe-qr>WW`9o>)ers{$!H^(|t?8 zDeYPO)OW4Y#g6?QA;ma7%(a1`oZ7Yc)Ou1kH}O&oB6(+Z zXBRXmt5^#;6T-Lgx3tlLX0~; z(obiFNyFYt3HlG(XR4FRrL7JVlYOIyW}G?7MwU2_SI9fP*mhpuA_lZxF7FT?WPy@j znLvM>rb%^PC|u~BQmI9BI}3C^S8#ZsalhV5&f$K}xGnrqzs4DPxpwuE*<2W_LtPZ` z6y2DC3S$cOY&yh#?gCEOFzUlPje|y$94ViG1qTP}@_!mS*5C)8qAlLr>zKMenQET9 zy}C0Px0IQiTB~5DJ$~Gj!A-r*PQrJKDqGu5alL^X9W6wI;6`jII~>a!`US#Q$o=?5 z$xfRQOx?kQy$xB8y85nm!QttNiE`=)JyKlrw=O~CCGD>hpKb?O)j0Tzs9&0YQ9Q%m zB<|JZY`6Z)_!1b*70(ouHW-q~0Lv)nsVB|pHb`8)7B@mn=wtTlt0Z&$+^<~8eSit- z{XQWCdMR17I`c8#6PsEXyrL%}DsQ%ZkL5Vc@zbNtL36IATdS^nOP|Vr-aLo<+WGx1 z!E}CSmKq1xs0Dj2_>Pr*h^e~c4))~vo|D{O- zxh=y}u!n3}H=>F%O!Xl+wQOfEmM;CBYWYaeO}Xskgw{l_+b^-L6Y$Jg3pX{);x{`D zY0(}aod?qqe-%d^L<>`<)!7mhMnS5`ez)B z1s*`n>cUpyACv}qO=P@((GZ(ZM4?244^O=@iLD7$|J~S1I%KF#6t@oz9(_zX#11x` zz$u2uO!&I^5IT5Za%yDT0a;)?0lls7ql9WEuldVLMbU3+p4oQC2Q0ZWyN^ef?OE>( z5aUx>iU-g-)ypLoZsYNHw`mAfXqS*9rr%5G;1!XK4u!<%ACc9p%VHOAD{vZmJih}oKwc9T&*jh zsUx=zmb@FJ9@}+?EzE5hj7sHWW34``55*bZ0f9x1LfM1f-c&p(Acxw>Cu7?lFMcCA z6%xDNRX=7%5i^-WCj7WkN&%y}jIDTBp3*kvqe(R+;6W@bUt9U}2yN2@<#4y?z;v2d zCWtMnjY=(I->m&ZW>XwnvvF|IdgBCW|Aixe!=_jVSxLZ)Es}(jjk-KL8JM#P z)|PlTL0rb)g}*>WHsGBs+pA$e8OJ4DPacAr@K=0y!2O1j3{$+8*kp`F=NX`OJP0vG znF!|^B~bcShl~Sw--PGRuc9boe9ojI+XHzk|CVzk3{1pKHYLMCs#57bU!AeZvEuYG zcd79HRfK}zjnL@xENX%g+6zHl3|stwQ_?q(jZZ5RIjTiWkqi}G&p$3yo;Hy&BL9mk*9m>jVJT!h6QuP%<}<9_flCCK`cDo}R8%0fM6 zGMX>dx;22YcT!BWw4^l*Han^xp|27I{UbUu(SFYJaHz4gB+49#(FWrhBgf@dbq=R^V zYC@cc|9W1auK5ZLdOnSBfNioWB()`)KAO%S_>+ZYuUoE3|M}~xwPg-54VEwVdi^~S z)D5td&NtZzrue$<;}K>>b1@C(s^H71-|XknktF#NSPI2{j?CMUSl1?VGDXaF-6<}= zxcr^Ixac0Y(0Kfi9lVonmQkI@W<~n0>IXCEN&EG5nuPc%CUC|c?v;o z3$vQ&l^a<)=(Jog`I8*^A}j)She>(M^l3!$;Nl^LRc)%CLl&bq+o$@vQY8}fZR zni^WWty>UOA;TShEw97k?HJ>%bd=ftNx+3&+{&?!&A7D<@gnx9;lT`}b)cqXZG-sv z{r2KdW?l(r?;0yJ2^P26F8USjJvIoy_IRtEu>_V$a{CCMAD?niJkSpN$aZnGslC|u6rNo_~1c%jdKR%`LvqtjBfJO zcvX<=Y}ZCu!HQV-$J4lXe0fQtOl7t0^zs%W+ev`gUt>nlB7K9L#S7nIl4etUA59hx z(c<#$!u$C`UxtOCRW#zu59f1J<>qEc=4=CBW`|;EIK{2nhL8qPlIfQ_0@g|eh=0Vz z;zRG}IUCR(DFP(p8Yse;garjYFHx(*B<8iagl1nq4~a&>mc4sC=zBPC)WIxt;P}-8 zOJ^`SnJQ2x{Ls>Gv?+jx4~xqa@)n<-N-Nb}+za0a`YS4BbL$({Fprkrn>`@(eE3`5 zEBh@^q4DP1)Eb5b0ME6#9-~8yWUbYDb%mzQf&@ zQhsivH~&DIvS+jNCL&LY1NdheMyBL`r9bN$-awu)(voAyQtAI7XWsR0_55Wdo;e*W zN4{7nZAym$fJ&2&_@whR%5e%j%d_oP#%qp&~NzaW%8MO=9X++UlB=j zUMz0=ewAAwAW0m+D7S^>C*iI%MdXQWKS!TIa%EweiF}p@tOCqPCII9?`itTTC&_Qq zDI(K)<@!l2c+^`S$RF-+JuS6|VZc;myueE{NWbB)uQ`unkx99aCKK>2_G3UJmM>KQ z@llsua-%G%{gaY*kNr&x_&Wg}%1_|V>plPYnltbr^cmbrm>-|tRS=BD!eGyVz_ zDjqgT+4Zc;=0M;A+c6v|H|C&T^#CP$g?#wlj`MVElYD!fwC$z$?p$U%ayrePM0?mj zN3Svm863QEVhs{M{nom_LI3rx!l-JT$j?6&+5My(C+2#tI>oG#X}7ge+t=%qf^r;; zi8y@ig9*lj!~caoZ>=)L<&KZyA4Umn9{=E2Q~$$hlOsMcH9IQTx5BUnI-~{BG1+KE99oia zMIHoy+t2%}1Eu+B=0ONvZCL+sRWw>sItCl%=N^&!csQk*@Grsh&Hn)};GvfQ literal 0 HcmV?d00001 diff --git a/thirdparty/cbrewer/interpolate_cbrewer.m b/thirdparty/cbrewer/interpolate_cbrewer.m new file mode 100644 index 0000000..4a8db91 --- /dev/null +++ b/thirdparty/cbrewer/interpolate_cbrewer.m @@ -0,0 +1,36 @@ +function [interp_cmap]=interpolate_cbrewer(cbrew_init, interp_method, ncolors) +% +% INTERPOLATE_CBREWER - interpolate a colorbrewer map to ncolors levels +% +% INPUT: +% - cbrew_init: the initial colormap with format N*3 +% - interp_method: interpolation method, which can be the following: +% 'nearest' - nearest neighbor interpolation +% 'linear' - bilinear interpolation +% 'spline' - spline interpolation +% 'cubic' - bicubic interpolation as long as the data is +% uniformly spaced, otherwise the same as 'spline' +% - ncolors=desired number of colors +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 14.10.2011 + + +% just to make sure, in case someone puts in a decimal +ncolors=round(ncolors); + +% How many data points of the colormap available +nmax=size(cbrew_init,1); + +% create the associated X axis (using round to get rid of decimals) +a=(ncolors-1)./(nmax-1); +X=round([0 a:a:(ncolors-1)]); +X2=0:ncolors-1; + +z=interp1(X,cbrew_init(:,1),X2,lower(interp_method)); +z2=interp1(X,cbrew_init(:,2),X2,lower(interp_method)); +z3=interp1(X,cbrew_init(:,3),X2, lower(interp_method)); +interp_cmap=round([z' z2' z3']); + +end \ No newline at end of file diff --git a/thirdparty/cbrewer/plot_brewer_cmap.m b/thirdparty/cbrewer/plot_brewer_cmap.m new file mode 100644 index 0000000..a5cab9e --- /dev/null +++ b/thirdparty/cbrewer/plot_brewer_cmap.m @@ -0,0 +1,50 @@ +% Plots and identifies the various colorbrewer tables available. +% Is called by cbrewer.m when no arguments are given. +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 14.10.2011 + + + +load('colorbrewer.mat') + +ctypes={'div', 'seq', 'qual'}; +ctypes_title={'Diverging', 'Sequential', 'Qualitative'}; +cnames{1,:}={'BrBG', 'PiYG', 'PRGn', 'PuOr', 'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral'}; +cnames{2,:}={'Blues','BuGn','BuPu','GnBu','Greens','Greys','Oranges','OrRd','PuBu','PuBuGn','PuRd',... + 'Purples','RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd'}; +cnames{3,:}={'Accent', 'Dark2', 'Paired', 'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3'}; + +figure('position', [314 327 807 420]) +for itype=1:3 + + %fh(itype)=figure(); + subplot(1,3,itype) + + for iname=1:length(cnames{itype,:}) + + ncol=length(colorbrewer.(ctypes{itype}).(cnames{itype}{iname})); + fg=1./ncol; % geometrical factor + + X=fg.*[0 0 1 1]; + Y=0.1.*[1 0 0 1]+(2*iname-1)*0.1; + F=cbrewer(ctypes{itype}, cnames{itype}{iname}, ncol); + + for icol=1:ncol + X2=X+fg.*(icol-1); + fill(X2,Y,F(icol, :), 'linestyle', 'none') + text(-0.1, mean(Y), cnames{itype}{iname}, 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 10, 'FontName' , 'AvantGarde') + xlim([-0.4, 1]) + hold all + end % icol + %set(gca, 'box', 'off') + title(ctypes_title{itype}, 'FontWeight', 'bold', 'FontSize', 16, 'FontName' , 'AvantGarde') + axis off + set(gcf, 'color', [1 1 1]) + end % iname + ylim([0.1 1.05.*max(Y)]); +end %itype + +set(gcf, 'MenuBar', 'none') +set(gcf, 'Name', 'ColorBrewer Color maps') \ No newline at end of file diff --git a/utils/autoreg_minphase.m b/utils/autoreg_minphase.m index e387363..1e4e3ce 100644 --- a/utils/autoreg_minphase.m +++ b/utils/autoreg_minphase.m @@ -1,4 +1,4 @@ -function EQ = autoreg_minphase(H,T,fs,maxAmp,frac,flims) +function [EQmp, EQ] = autoreg_minphase(H,T,fs,maxAmp,frac,flims) % Regularized inverse of an impulse response [REFERENCE PENDING]. % % IMPORTANT: this function does not deal with phase and the output is @@ -17,7 +17,8 @@ % flims = inversion limits in Hz (def=[100 18000]) % % OUTPUT: -% EQ = minimum phase EQ filter (nfreqs x nchannels) +% EQmp = minimum phase EQ filter (nfreqs x nchannels) +% EQ = linear phase EQ filter (nfrqes x nchannels) % % REFERENCES: % TODO @@ -37,19 +38,20 @@ end %% Calculate the inverse filter -H = abs(H); % disregard phase -T = abs(T); % disregard phase X = H./T; % this is what we want to invert +Xmag = abs(X); nfreqs = size(H,1); f = linspace(0,fs/2,nfreqs).'; alpha = 1/(2*db2mag(maxAmp))^2; -Xhat = AKfractOctSmooth(X,'amp',fs,frac,false); +Xhat = AKfractOctSmooth(Xmag,'amp',fs,frac,false); sigma = zeros(nfreqs,1); -sigma(Xhat>=X) = Xhat(Xhat>=X)-X(Xhat>=X); +sigma(Xhat>=Xmag) = Xhat(Xhat>=Xmag)-Xmag(Xhat>=Xmag); -EQ = X ./ ( X.^2 + alpha + sigma.^2 ); +EQ = conj(X) ./ ( abs(X).^2 + alpha + sigma.^2 ); +EQmag = abs(EQ); +EQph = angle(EQ); %% Fade in/out to 1 outside the inversion limits % Low-frequency boundary @@ -65,7 +67,7 @@ fcvec = f(fc1_ind:fc2_ind); % f vector (lin scale) w = interp1(fcvec_log,w_log,fcvec,'spline'); % weights (lin scale) w = [zeros(fc1_ind-1,1);w(:);ones(nfreqs-fc2_ind,1)]; % for all freqs - EQ = (1-w)*1 + w.*EQ; % weighted sum + EQmag = (1-w)*1 + w.*EQmag; % weighted sum end % High-frequency boundary @@ -81,10 +83,11 @@ fcvec = f(fc1_ind:fc2_ind); % f vector (lin scale) w = interp1(fcvec_log,w_log,fcvec,'spline'); % weights (lin scale) w = [ones(fc1_ind-1,1);w(:);zeros(nfreqs-fc2_ind,1)]; % for all freqs - EQ = (1-w)*1 + w.*EQ; % weighted sum + EQmag = (1-w)*1 + w.*EQmag; % weighted sum end +EQ = EQmag.*exp(1i*EQph); % preserve original phase %% Make minimum phase -EQ = makeMinPhase(EQ); +EQmp = makeMinPhase(abs(EQ)); end diff --git a/utils/getRealSHmatrix.m b/utils/getRealSHmatrix.m new file mode 100644 index 0000000..ca9a429 --- /dev/null +++ b/utils/getRealSHmatrix.m @@ -0,0 +1,41 @@ +function Y = getRealSHmatrix(N,az,el) +% Get matrix with real spherical harmonics. +% INPUT: +% N = SH order +% az = azimuth (ndirs x 1) in rad +% el = elevation (ndirs x 1) in rad (0=top, pi/2=front) +% +% OUTPUT: +% Y = real SH matrix ((N+1)^2 x ndirs) +% +% AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk +% August 2021 + +az = az(:).'; % force row +el = el(:).'; + +nsh = (N+1)^2; +ndirs = size(az,2); + +Y = zeros(nsh,ndirs); + +for n=0:N + for m=-n:n + ind = n^2 + n + m + 1; % ACN ordering + % Elevation term + Pn = legendre(n,cos(el.')); + % Azimuth term + if m<0 + ym = sqrt(2)*sin(abs(m)*az); + elseif m>0 + ym = sqrt(2)*cos(m*az); + else + ym = ones(1,ndirs); + end + % Normalisation term (N3D) + Norm = sqrt( (2*n+1)/(4*pi) * factorial(n-abs(m))/factorial(n+abs(m))); + % Final formula + Y(ind,:) = (-1)^m * Norm .* Pn(abs(m)+1,:) .* ym; + end +end + diff --git a/utils/getSphDist.m b/utils/getSphDist.m new file mode 100644 index 0000000..fe83123 --- /dev/null +++ b/utils/getSphDist.m @@ -0,0 +1,7 @@ +function dist = getSphDist(azA,elA,azB,elB) +% Great circle distance between two points in spherical coordinates (rad) +elA = pi/2 - elA; % elevation is expected as (pi/2=front) but the formula below assumes (0=front) +elB = pi/2 - elB; +dist = acos(sin(elA).*sin(elB)+cos(elA).*cos(elB).*cos(azA-azB)); + + diff --git a/utils/hrtf2sofa.m b/utils/hrtf2sofa.m index b0bb134..4615c88 100644 --- a/utils/hrtf2sofa.m +++ b/utils/hrtf2sofa.m @@ -44,4 +44,5 @@ SOFA_obj.API.Dimensions.Data.IR = 'MRN'; % time-measurements-receivers SOFA_obj.GLOBAL_SOFAConventions = 'SimpleFreeFieldHRIR'; SOFA_obj.GLOBAL_RoomType = 'free field'; -SOFA_obj.GLOBAL_Comment = 'Generated with ''hrtf2sofa()'' from an HRIR'; \ No newline at end of file +SOFA_obj.GLOBAL_Comment = 'Generated with ''hrtf2sofa()'' from an HRIR'; +SOFA_obj.GLOBAL_Title = 'HRTF'; \ No newline at end of file diff --git a/utils/sofa2hrtf.m b/utils/sofa2hrtf.m index e709dc7..5bf163c 100644 --- a/utils/sofa2hrtf.m +++ b/utils/sofa2hrtf.m @@ -11,7 +11,7 @@ % fs = sampling frequency in Hz % az = azimuth in rad (dirs x 1) % el = colatitude in rad (dirs x 1) -% r = head radius in m (def=0.085) +% r = head radius in m (def=0.0875) % % AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk % February 2021 @@ -33,6 +33,9 @@ az = SOFA_obj.SourcePosition(:,1)*pi/180; % azimuth in rad el = pi/2-SOFA_obj.SourcePosition(:,2)*pi/180; % colatitude in rad r = 0.0875; % using this as default for now +% For negative collatitudes, reverse them and add pi to the azimuth +az(el<0) = mod(az(el<0)+pi,2*pi); +el(el<0) = -el(el<0); %% Rearrange IR dimensions dimIn = SOFA_obj.API.Dimensions.Data.IR; % dimensions, e.g. 'MRN'