Commit cea4e2a9e9c995a948644fb96e8da645257a33df

Authored by Jean-Luc Shaw
1 parent cc34cf0a
Exists in master

Several new additions. rho_t.m calculates the circular correlation coefficient f…

…rom Fisher and Lee (1983). genopt.m is a general optimisation routine using a genetic algorythm. seabirdCTD.m is designed to be a low fuss reader for seabird .cnv files.
math/genopt.m 0 → 100644
... ... @@ -0,0 +1,66 @@
  1 +function MOM = genopt(Q,X,dom,NP,NI,pc,pm)
  2 +% Synthax: MOM = genopt(@Q,@X,dom,NP,NI,pc,pm)
  3 +%
  4 +% Optimisation routine following a genetic algorithm. Returns the best parameter
  5 +% set which solves the problem expressed by the anonymous function 'Q', when
  6 +% either the criteria expressed by the anonymous function 'X' has been met or
  7 +% NI iterations have been completed. 'X' should be formulated in terms of the
  8 +% values returned by 'Q'. The 'dom' variable is a matrix of size [p 2], where
  9 +% p is the number of parameters passed to the function 'Q'. It's rows are the
  10 +% domains from which the initial population will be generated for each parameter.
  11 +% Variable 'NP' is the size of the population and must be an even number. Variables
  12 +% 'pc' and 'pm' are respectively the probability of mixing genetic material and
  13 +% probability of mutation. Recommended values are 0.8 and 0.2 .
  14 +
  15 +%INIT POPULATION
  16 +NA = numel(dom(:,1)) ;
  17 +PP = nan(NP,NA) ;
  18 +SC = nan(NP,1) ;
  19 +RK = nan(NP,1) ;
  20 +for ii = 1:NP
  21 + for jj = 1:NA
  22 + PP(ii,jj) = dom(jj,1) + (dom(jj,2)-dom(jj,1))*rand(1) ;
  23 + end
  24 +end
  25 +
  26 +for kk = 1:NI
  27 +
  28 + %EVALUATION
  29 + for ii = 1:NP ; SC(ii) = Q(PP(ii,:)) ; end
  30 +
  31 + %RANKING
  32 + [HS,RK] = sort(SC) ;
  33 + PP = PP(RK,:) ;
  34 +
  35 + %TEST
  36 + if X(HS(1)) ; break; end
  37 +
  38 + %SELECTION
  39 + DICE = randi([2 4]) ;
  40 + MOM = PP(1, :) ;
  41 + DAD = PP(DICE,:) ;
  42 +
  43 + %REPRODUCTION
  44 + PP(1,:) = MOM ;
  45 + PP(2,:) = DAD ;
  46 + for ii = 3:2:NP
  47 + for jj = 1:NA
  48 + % CROISEMENT
  49 + if rand(1) < pc
  50 + DICE = rand(1) ;
  51 + PP(ii,jj) = MOM(jj)*DICE + DAD(jj)*(1-DICE) ;
  52 + PP(ii+1,jj) = DAD(jj)*DICE + MOM(jj)*(1-DICE) ;
  53 + else
  54 + PP(ii,jj) = MOM(jj) ;
  55 + PP(ii+1,jj) = DAD(jj) ;
  56 + end
  57 +
  58 + % MUTATION
  59 + if rand(1) < pm
  60 + PP(ii,jj) = dom(jj,1) + (dom(jj,2)-dom(jj,1))*rand(1) ;
  61 + end
  62 + end
  63 + end
  64 +end
  65 +
  66 +end
... ...
math/nlinfit.m
... ... @@ -27,6 +27,7 @@ elseif sz(1) == 1 &amp; sz(2) &gt; 1 | sz(2) == 1 &amp; sz(1) &gt; 1
27 27 start = pars ;
28 28 end
29 29  
  30 +
30 31 % Minimize parameters
31 32 A = fminsearch(@(a)sum( (y - rhs(a,x,varargin{:})).^2 ),start) ;
32 33 yfit = rhs(A,x,varargin{:}) ;
... ...
math/rho_t.m 0 → 100644
... ... @@ -0,0 +1,33 @@
  1 +function [rho] = rho_t(x,y,degrad)
  2 +% [rho] = rho_t(x,y,degrad)
  3 +%
  4 +% Explicitly calculates the correlation coefficient rho_t detailed in
  5 +% Fisher and Lee (1983), between the two circular variables 'x' and 'y'.
  6 +% Variable 'degrad' is a string and should be "deg" if 'x' and 'y' are
  7 +% in degrees or 'rad' if they are in radians.
  8 +%
  9 +% Reference: Fisher, N.I., LEE, A.J.,(1983), A correlation coefficient
  10 +% for circular data, Biometrika, 70, 2, pp. 327-32
  11 +
  12 +
  13 +numerator = 0 ;
  14 +denA = 0 ;
  15 +denB = 0 ;
  16 +
  17 +switch degrad
  18 + case 'deg'
  19 + x = x*pi/180 ;
  20 + y = -y*pi/180 ;
  21 + case 'rad'
  22 + % do nothing
  23 +end
  24 +
  25 +for ii = 1:numel(x)-1
  26 + numerator = numerator + sum( sin(x(ii)-x(ii+1:end)).*sin(y(ii)-y(ii+1:end)) ) ;
  27 + denA = denA + sum( sin(x(ii)-x(ii+1:end)).^2 );
  28 + denB = denB + sum( sin(y(ii)-y(ii+1:end)).^2 ) ;
  29 +end
  30 +
  31 +rho = numerator./( sqrt(denA).*sqrt(denB) ) ;
  32 +
  33 +end
... ...
util/seabirdCTD.m
... ... @@ -179,7 +179,7 @@ z = movmean(A{:,Iz},10/sr) ;
179 179 if startsWith(showz,'y'); figure; plot(dnum,z,'b'); hold on; end
180 180  
181 181 dz = gradient(z,dnum*24*3600) ; % compute downward velocity
182   -z( [1:nvalues]' < warmup/sr | dz<=0.001 ) = nan ; % get rid of start and end values
  182 +z( [1:nvalues]' < warmup/sr | dz<=0.005 ) = nan ; % get rid of start and end values
183 183 [~,Im] = min(z) ; % find min depth
184 184 [~,IM] = max(z) ; % find max depth
185 185 I = find(dz > dzthres & [1:nvalues]' > Im & [1:nvalues]' < IM) ; % select down-going between min z and max z
... ...