Commit fa2630c0 authored by Jean-Luc Shaw's avatar Jean-Luc Shaw

Added a current extrapolation routine currextrap.m which uses the genetic...

Added a current extrapolation routine currextrap.m which uses the genetic optimisation routine genopt.m . Added a wrapper function for gradient.m that calculates horizontal current shear given a vector, matrix of 3d array because I am just that lazy.
parent cea4e2a9
function MOM = genopt(Q,X,dom,NP,NI,pc,pm)
% Synthax: MOM = genopt(@Q,@X,dom,NP,NI,pc,pm)
function [A,HS] = genopt(Q,C,dom,NP,NI,pc,pm)
% Synthax: [A,HS] = genopt(@Q,@X,dom,NP,NI,pc,pm)
%
% Optimisation routine following a genetic algorithm. Returns the best parameter
% set which solves the problem expressed by the anonymous function 'Q', when
......@@ -33,12 +33,11 @@ for kk = 1:NI
PP = PP(RK,:) ;
%TEST
if X(HS(1)) ; break; end
if HS(1) < C ; break; end
%SELECTION
DICE = randi([2 4]) ;
MOM = PP(1, :) ;
DAD = PP(DICE,:) ;
DAD = PP(randi([2 4]),:) ;
%REPRODUCTION
PP(1,:) = MOM ;
......@@ -62,5 +61,9 @@ for kk = 1:NI
end
end
end
disp(['Last iteration: ' num2str(kk)]) ;
% OUTPUT
A = MOM ;
HS = HS(1) ;
end
function [ue,F] = currextrap(u,N2,p,z,varargin)
% Synthax: [ue,F] = currextrap(u,N2,p,z,varargin)
%
% Uses the dynmodes routine and genetic optimisation to fit the first
% 4 vertical modes on current profile 'u', where 'N2' is the buoyancy
% frequency, 'p' is the pressure, and 'z' is the depth (positive) of
% the associated current profile.
%
% Output 'ue' is the current profile fit provided and 'F' is a 5
% element vector representing the weights of the 4 first modes plus
% the speed offset. This is the solution set of the genetic optimising
% routine.
%
% Optional arguments should be supplied as parameter/value pairs.
% Supported options are:
%
% - section,[string] Give higher importance to fitting the
% surface 'top', bottom 'bottom', or the
% entire velocity profile 'all'. In
% effect this applies a weighing function
% which decreases as zc/z, increases as
% z/zc, where zc is a constant,
% or is flat with respect to depth
% to the compute the scoring function.
%
% - zc,[float] Depth where the weighing function becomes
% unity. Defaults to 1.
%
% - itermax,[int] Set the maximum number of iterations
% of the genetic algorithm. Defaults to
% 500.
%
% - popsize,[int] Set the population size for the genetic
% algorithm. Defaults to 20.
%
% - nsol,[int] Number of solution sets computed and
% from the genetic algorithm. Only the
% highest scoring solution is returned
% Defaults to 1.
% Dependecies:
% - dynmodes.m : https://github.com/sea-mat/dynmodes
% - genopt.m : https://github.com/jeanlucshaw/matlab_stuff/blob/master/math/genopt.m
%OPTIONS
%default
zc = 1 ;
NP = 20 ; % Genetic algo population size
NI = 500 ; % Genetic algo max iterations
NS = 1 ; % Genetic algo number of solutions averaged
%user
for ii = 1:2:nargin-5
switch varargin{ii}
case 'section'
switch varargin{ii+1}
case 'top'
W = zc/z;
case 'bottom'
W = z/zc;
case 'all'
W = ones(size(z));
end
case 'itermax'
NI = varargin{ii+1} ;
case 'popsize'
NP = varargin{ii+1} ;
case 'nsol'
NS = varargin{ii+1} ;
case 'zc'
zc = varargin{ii+1} ;
end
end
I = find(isfinite(N2)) ;
[wmodes,pmodes,ce]=dynmodes(N2(I)',p(I)',1) ;
%REGRID
pm = interp1(1:0.5:1+0.5*numel(I),pmodes,z,'linear','extrap') ;
%FIT
I = isfinite(u) ;
Q = @(a) sum( W(I).*sqrt( (a(1)*pm(I,1)+a(2)*pm(I,2)+a(3)*pm(I,3)+a(4)*pm(I,4)+a(5) - u(I)).^2 ) ) ;
A = nan([NS,5]) ;
HS = nan([NS,1]) ;
for ii = 1:NS
[A(ii,:),HS(ii)] = genopt(Q,0.01,[-0.2 0.2;-0.2 0.2;-0.2 0.2;-0.2 0.2;-2 2],NP,NI,0.8,0.2) ;
end
%OUTPUT
[~,I] = min(HS) ;
F = A(I,:) ;
ue = F(1)*pm(:,1)+F(2)*pm(:,2)+F(3)*pm(:,3)+F(4)*pm(:,4)+F(5) ;
end
function S = uvz2S(u,v,z)
% Synthax: S = uvz2S(u,v,z)
%
% Wrapper for the gradient function which takes as input
% horizontal velocities 'u' and 'v' as well as their
% corresponding depths 'z'. Note that 'u','v', and 'z'
% are all assumed to be of same size.
%
% Returns the horizontal shear:
%
% S = sqrt( (du/dz)^2 + (dv/dz)^2 )
%
D = ndims(u) ;
M = nan([1 D]) ;
M = size(u) ;
I = find( M > 1 ) ;
switch numel(I) ;
case 1
[dudz] = gradient(u,z) ;
[dvdz] = gradient(v,z) ;
case 2
[~,dudz] = gradient(u,1,z) ;
[~,dvdz] = gradient(v,1,z) ;
case 3
[~,dudz,~] = gradient(u,1,z,1) ;
[~,dvdz,~] = gradient(v,1,z,1) ;
end
S = sqrt( dudz.^2 + dvdz.^2 ) ;
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment