text:
function [K] = ku0u0(x, y, xp, yp, hyp, ubar, vbar, ubarp, vbarp, dt, i)

logsigma = hyp(1);
logthetax = hyp(2);
logthetay = hyp(3);

a1 = hyp(4);
a2 = hyp(5);

logsigmap = hyp(6);
logthetaxp = hyp(7);
logthetayp = hyp(8);

n_x = size(x,1);
n_y = size(y,1);
n_xp = size(xp,1);
n_yp = size(yp,1);

x = repmat(x,1,n_xp);
y = repmat(y,1,n_yp);
xp = repmat(xp',n_x,1);
yp = repmat(yp',n_y,1);

ubar = repmat(ubar,1,n_xp);
vbar = repmat(vbar,1,n_yp);
ubarp = repmat(ubarp',n_x,1);
vbarp = repmat(vbarp',n_y,1);

switch i

case 0

K=dt.^2.*exp(1).^(logsigmap+(-2).*logthetaxp+(-1/2).*exp(1).^((-1).* ...
  logthetaxp).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetayp).*(y+(-1) ...
  .*yp).^2).*(exp(1).^logthetaxp+(-1).*(x+(-1).*xp).^2)+a1.*dt.*exp(1).^( ...
  logsigma+(-5).*logthetay+(-1/2).*exp(1).^((-1).*logthetax).*(x+(-1).*xp) ...
  .^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2).*(exp(1).^((-2) ...
  .*logthetax+logthetay).*ubar.*((-1).*exp(1).^logthetax.*(x+(-1).*xp).*( ...
  exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.* ...
  dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^(2.*logthetay).* ...
  ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).* ...
  yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).* ...
  exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).* ...
  logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+ ...
  3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1) ...
  .^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+a1.*dt.*exp(1).^( ...
  logthetax+2.*logthetay).*ubarp.*(exp(1).^logthetay+(-1).*(y+(-1).*yp) ...
  .^2)+(-2).*dt.*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2))+vbar.*(exp(1).^logthetay.*(a1.*dt.* ...
  exp(1).^((-1).*logthetax+logthetay).*(3.*exp(1).^(logthetax+logthetay).* ...
  vbarp+(-2).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).* ...
  exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^2)+(-2).*exp(1).^(2.*logthetay) ...
  .*(y+(-1).*yp)+2.*dt.*exp(1).^(a2+(-2).*logthetax).*((-6).*exp(1).^(2.* ...
  logthetax+logthetay)+(-1).*exp(1).^(logthetax+2.*logthetay)+exp(1).^(2.* ...
  logthetay).*(x+(-1).*xp).^2+2.*exp(1).^(2.*logthetax).*(y+(-1).*yp).^2) ...
  .*(y+(-1).*yp))+(-1).*(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).* ...
  (y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^( ...
  2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*(y+ ...
  (-1).*yp)))+exp(1).^(logsigma+(-4).*logthetay+(-1/2).*exp(1).^((-1).* ...
  logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).* ...
  yp).^2).*(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp) ...
  .^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^(2.* ...
  logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+(-1) ...
  .*dt.*exp(1).^(a2+logsigma+(-6).*logthetay+(-1/2).*exp(1).^((-1).* ...
  logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).* ...
  yp).^2).*((-2).*exp(1).^((-2).*logthetax+2.*logthetay).*(exp(1).^(2.*( ...
  logthetax+logthetay))+dt.*(6.*exp(1).^(a2+2.*logthetax+logthetay)+exp(1) ...
  .^(a2+logthetax+2.*logthetay)+a1.*exp(1).^(logthetax+2.*logthetay).* ...
  ubarp.*(x+(-1).*xp)+(-1).*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).^2+3.* ...
  a1.*exp(1).^(2.*logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-6).*exp(1) ...
  .^(a2+2.*logthetax).*(y+(-1).*yp).^2))+exp(1).^((-2).*logthetax+2.* ...
  logthetay).*((-1).*exp(1).^logthetax+(x+(-1).*xp).^2).*(exp(1).^(2.* ...
  logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^(( ...
  -1).*logthetax+logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+ ...
  3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1) ...
  .^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*( ...
  3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+( ...
  -1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.* ...
  logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+(-2).*dt.*exp(1) ...
  .^(a2+(-2).*logthetax+4.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).* ...
  yp).^2)+(-2).*dt.*exp(1).^((-3).*logthetax+4.*logthetay).*(a1.*exp(1) ...
  .^logthetax.*ubarp+(-2).*exp(1).^a2.*(x+(-1).*xp)).*(x+(-1).*xp).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2)+(exp(1).^(2.*logthetay).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4)).*((-1).*exp(1).^logthetay+(y+(-1).*yp).^2) ...
  +(-2).*exp(1).^logthetay.*(a1.*dt.*exp(1).^((-1).*logthetax+logthetay).* ...
  (3.*exp(1).^(logthetax+logthetay).*vbarp+(-2).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).*exp(1).^logthetax.*vbarp.*(y+( ...
  -1).*yp).^2)+(-2).*exp(1).^(2.*logthetay).*(y+(-1).*yp)+2.*dt.*exp(1).^( ...
  a2+(-2).*logthetax).*((-6).*exp(1).^(2.*logthetax+logthetay)+(-1).*exp( ...
  1).^(logthetax+2.*logthetay)+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+2.* ...
  exp(1).^(2.*logthetax).*(y+(-1).*yp).^2).*(y+(-1).*yp)).*(y+(-1).*yp)); ...

case 1 % logsigma

K=exp(1).^(logsigma+(-6).*logthetay+(-1/2).*exp(1).^((-1).*logthetax).*(x+ ...
  (-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2).*(a1.* ...
  dt.*exp(1).^logthetay.*(exp(1).^((-2).*logthetax+logthetay).*ubar.*((-1) ...
  .*exp(1).^logthetax.*(x+(-1).*xp).*(exp(1).^(2.*logthetay).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4))+a1.*dt.*exp(1).^(logthetax+2.*logthetay).* ...
  ubarp.*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+(-2).*dt.*exp(1).^(a2+ ...
  2.*logthetay).*(x+(-1).*xp).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2))+ ...
  vbar.*(exp(1).^logthetay.*(a1.*dt.*exp(1).^((-1).*logthetax+logthetay).* ...
  (3.*exp(1).^(logthetax+logthetay).*vbarp+(-2).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).*exp(1).^logthetax.*vbarp.*(y+( ...
  -1).*yp).^2)+(-2).*exp(1).^(2.*logthetay).*(y+(-1).*yp)+2.*dt.*exp(1).^( ...
  a2+(-2).*logthetax).*((-6).*exp(1).^(2.*logthetax+logthetay)+(-1).*exp( ...
  1).^(logthetax+2.*logthetay)+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+2.* ...
  exp(1).^(2.*logthetax).*(y+(-1).*yp).^2).*(y+(-1).*yp))+(-1).*(exp(1).^( ...
  2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1) ...
  .^((-1).*logthetax+logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).* ...
  xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1) ...
  .^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*( ...
  3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+( ...
  -1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.* ...
  logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*(y+(-1).*yp)))+ ...
  exp(1).^(2.*logthetay).*(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1) ...
  .*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1) ...
  .^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+(-1) ...
  .*dt.*exp(1).^a2.*((-2).*exp(1).^((-2).*logthetax+2.*logthetay).*(exp(1) ...
  .^(2.*(logthetax+logthetay))+dt.*(6.*exp(1).^(a2+2.*logthetax+logthetay) ...
  +exp(1).^(a2+logthetax+2.*logthetay)+a1.*exp(1).^(logthetax+2.* ...
  logthetay).*ubarp.*(x+(-1).*xp)+(-1).*exp(1).^(a2+2.*logthetay).*(x+(-1) ...
  .*xp).^2+3.*a1.*exp(1).^(2.*logthetax+logthetay).*vbarp.*(y+(-1).*yp)+( ...
  -6).*exp(1).^(a2+2.*logthetax).*(y+(-1).*yp).^2))+exp(1).^((-2).* ...
  logthetax+2.*logthetay).*((-1).*exp(1).^logthetax+(x+(-1).*xp).^2).*( ...
  exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.* ...
  dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^(2.*logthetay).* ...
  ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).* ...
  yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).* ...
  exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).* ...
  logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+ ...
  3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1) ...
  .^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+(-2).*dt.*exp(1) ...
  .^(a2+(-2).*logthetax+4.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).* ...
  yp).^2)+(-2).*dt.*exp(1).^((-3).*logthetax+4.*logthetay).*(a1.*exp(1) ...
  .^logthetax.*ubarp+(-2).*exp(1).^a2.*(x+(-1).*xp)).*(x+(-1).*xp).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2)+(exp(1).^(2.*logthetay).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4)).*((-1).*exp(1).^logthetay+(y+(-1).*yp).^2) ...
  +(-2).*exp(1).^logthetay.*(a1.*dt.*exp(1).^((-1).*logthetax+logthetay).* ...
  (3.*exp(1).^(logthetax+logthetay).*vbarp+(-2).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).*exp(1).^logthetax.*vbarp.*(y+( ...
  -1).*yp).^2)+(-2).*exp(1).^(2.*logthetay).*(y+(-1).*yp)+2.*dt.*exp(1).^( ...
  a2+(-2).*logthetax).*((-6).*exp(1).^(2.*logthetax+logthetay)+(-1).*exp( ...
  1).^(logthetax+2.*logthetay)+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+2.* ...
  exp(1).^(2.*logthetax).*(y+(-1).*yp).^2).*(y+(-1).*yp)).*(y+(-1).*yp))); ...

case 2 % logthetax

K=(1/2).*exp(1).^(logsigma+(-6).*logthetay+(-1/2).*exp(1).^((-1).* ...
  logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).* ...
  yp).^2).*(a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(x+(-1).*xp).^2.* ...
  (exp(1).^((-2).*logthetax+logthetay).*ubar.*((-1).*exp(1).^logthetax.*( ...
  x+(-1).*xp).*(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).* ...
  yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^(2.* ...
  logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+a1.* ...
  dt.*exp(1).^(logthetax+2.*logthetay).*ubarp.*(exp(1).^logthetay+(-1).*( ...
  y+(-1).*yp).^2)+(-2).*dt.*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2))+vbar.*(exp(1).^logthetay.*(a1.* ...
  dt.*exp(1).^((-1).*logthetax+logthetay).*(3.*exp(1).^(logthetax+ ...
  logthetay).*vbarp+(-2).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1) ...
  .*yp)+(-3).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^2)+(-2).*exp(1).^( ...
  2.*logthetay).*(y+(-1).*yp)+2.*dt.*exp(1).^(a2+(-2).*logthetax).*((-6).* ...
  exp(1).^(2.*logthetax+logthetay)+(-1).*exp(1).^(logthetax+2.*logthetay)+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+2.*exp(1).^(2.*logthetax).*(y+( ...
  -1).*yp).^2).*(y+(-1).*yp))+(-1).*(exp(1).^(2.*logthetay).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4)).*(y+(-1).*yp)))+exp(1).^((-1).*logthetax+ ...
  2.*logthetay).*(x+(-1).*xp).^2.*(exp(1).^(2.*logthetay).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4))+(-2).*a1.*dt.*exp(1).^((-3).*logthetax+2.* ...
  logthetay).*((-1).*exp(1).^(2.*(logthetax+logthetay)).*ubar.*(x+(-1).* ...
  xp).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^( ...
  logthetax+logthetay).*(exp(1).^(logthetax+2.*logthetay).*ubar.*ubarp+( ...
  -2).*exp(1).^(2.*logthetay).*ubar.*ubarp.*(x+(-1).*xp).^2+(-1).*exp(1) ...
  .^(logthetax+logthetay).*(3.*ubar.*vbarp.*(x+(-1).*xp)+ubarp.*(3.*vbar.* ...
  (x+(-1).*xp)+ubar.*(y+(-1).*yp))).*(y+(-1).*yp)+2.*exp(1).^logthetay.* ...
  ubar.*ubarp.*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^logthetax.*( ...
  ubarp.*vbar+ubar.*vbarp).*(x+(-1).*xp).*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^a2.*((-6).*exp(1).^(logthetax+3.*logthetay).*ubar.*(x+(-1).*xp)+3.* ...
  exp(1).^(3.*logthetay).*ubar.*(x+(-1).*xp).^3+(-3).*exp(1).^(2.*( ...
  logthetax+logthetay)).*(ubar.*(x+(-1).*xp)+vbar.*(y+(-1).*yp))+6.*exp(1) ...
  .^(logthetax+2.*logthetay).*(x+(-1).*xp).*(vbar.*(x+(-1).*xp)+ubar.*(y+( ...
  -1).*yp)).*(y+(-1).*yp)+(-3).*exp(1).^(2.*logthetay).*ubar.*(x+(-1).*xp) ...
  .^3.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax+logthetay).*(6.*ubar.*(x+(-1) ...
  .*xp)+vbar.*(y+(-1).*yp)).*(y+(-1).*yp).^2+(-2).*exp(1).^(logthetax+ ...
  logthetay).*vbar.*(x+(-1).*xp).^2.*(y+(-1).*yp).^3+(-1).*exp(1).^(2.* ...
  logthetax).*ubar.*(x+(-1).*xp).*(y+(-1).*yp).^4))+(-2).*dt.*exp(1).^(a2+ ...
  (-4).*logthetax+2.*logthetay).*(exp(1).^(2.*(logthetax+logthetay)).*( ...
  exp(1).^logthetax+(-2).*(x+(-1).*xp).^2).*(exp(1).^logthetay+(-1).*(y+( ...
  -1).*yp).^2)+dt.*(6.*exp(1).^(a2+3.*logthetax+2.*logthetay)+6.*exp(1).^( ...
  a2+2.*logthetax+3.*logthetay)+6.*a1.*exp(1).^(2.*logthetax+3.*logthetay) ...
  .*ubarp.*(x+(-1).*xp)+(-18).*exp(1).^(a2+logthetax+3.*logthetay).*(x+( ...
  -1).*xp).^2+(-3).*a1.*exp(1).^(logthetax+3.*logthetay).*ubarp.*(x+(-1).* ...
  xp).^3+4.*exp(1).^(a2+3.*logthetay).*(x+(-1).*xp).^4+3.*a1.*exp(1).^(3.* ...
  logthetax+2.*logthetay).*(ubarp.*(x+(-1).*xp)+vbarp.*(y+(-1).*yp))+(-6) ...
  .*exp(1).^(a2+2.*(logthetax+logthetay)).*(2.*x.^2+(-4).*x.*xp+2.*xp.^2+( ...
  y+(-1).*yp).^2)+(-6).*a1.*exp(1).^(2.*(logthetax+logthetay)).*(x+(-1).* ...
  xp).*(vbarp.*(x+(-1).*xp)+ubarp.*(y+(-1).*yp)).*(y+(-1).*yp)+(-12).*exp( ...
  1).^(a2+3.*logthetax+logthetay).*(y+(-1).*yp).^2+24.*exp(1).^(a2+2.* ...
  logthetax+logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+18.*exp(1).^(a2+ ...
  logthetax+2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+3.*a1.*exp(1) ...
  .^(logthetax+2.*logthetay).*ubarp.*(x+(-1).*xp).^3.*(y+(-1).*yp).^2+(-4) ...
  .*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).^4.*(y+(-1).*yp).^2+(-1).*a1.* ...
  exp(1).^(3.*logthetax+logthetay).*(6.*ubarp.*(x+(-1).*xp)+vbarp.*(y+(-1) ...
  .*yp)).*(y+(-1).*yp).^2+2.*a1.*exp(1).^(2.*logthetax+logthetay).*vbarp.* ...
  (x+(-1).*xp).^2.*(y+(-1).*yp).^3+2.*exp(1).^(a2+3.*logthetax).*(y+(-1).* ...
  yp).^4+a1.*exp(1).^(3.*logthetax).*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^4+ ...
  (-4).*exp(1).^(a2+2.*logthetax).*(x+(-1).*xp).^2.*(y+(-1).*yp).^4))+(-1) ...
  .*dt.*exp(1).^(a2+(-1).*logthetax).*(x+(-1).*xp).^2.*((-2).*exp(1).^(( ...
  -2).*logthetax+2.*logthetay).*(exp(1).^(2.*(logthetax+logthetay))+dt.*( ...
  6.*exp(1).^(a2+2.*logthetax+logthetay)+exp(1).^(a2+logthetax+2.* ...
  logthetay)+a1.*exp(1).^(logthetax+2.*logthetay).*ubarp.*(x+(-1).*xp)+( ...
  -1).*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).^2+3.*a1.*exp(1).^(2.* ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-6).*exp(1).^(a2+2.* ...
  logthetax).*(y+(-1).*yp).^2))+exp(1).^((-2).*logthetax+2.*logthetay).*(( ...
  -1).*exp(1).^logthetax+(x+(-1).*xp).^2).*(exp(1).^(2.*logthetay).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4))+(-2).*dt.*exp(1).^(a2+(-2).*logthetax+4.* ...
  logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+(-2).*dt.*exp(1) ...
  .^((-3).*logthetax+4.*logthetay).*(a1.*exp(1).^logthetax.*ubarp+(-2).* ...
  exp(1).^a2.*(x+(-1).*xp)).*(x+(-1).*xp).*(exp(1).^logthetay+(-1).*(y+( ...
  -1).*yp).^2)+(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).* ...
  yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^(2.* ...
  logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*(( ...
  -1).*exp(1).^logthetay+(y+(-1).*yp).^2)+(-2).*exp(1).^logthetay.*(a1.* ...
  dt.*exp(1).^((-1).*logthetax+logthetay).*(3.*exp(1).^(logthetax+ ...
  logthetay).*vbarp+(-2).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1) ...
  .*yp)+(-3).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^2)+(-2).*exp(1).^( ...
  2.*logthetay).*(y+(-1).*yp)+2.*dt.*exp(1).^(a2+(-2).*logthetax).*((-6).* ...
  exp(1).^(2.*logthetax+logthetay)+(-1).*exp(1).^(logthetax+2.*logthetay)+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+2.*exp(1).^(2.*logthetax).*(y+( ...
  -1).*yp).^2).*(y+(-1).*yp)).*(y+(-1).*yp))+(-2).*dt.*exp(1).^((-2).* ...
  logthetax+4.*logthetay).*(exp(1).^(a2+logthetax)+a1.*exp(1).^logthetax.* ...
  ubarp.*(x+(-1).*xp)+(-2).*exp(1).^a2.*(x+(-1).*xp).^2).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2));

case 3 % logthetay

K=exp(1).^(logsigma+(-6).*logthetay+(-1/2).*exp(1).^((-1).*logthetax).*(x+ ...
  (-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2).*(exp( ...
  1).^((-2).*logthetax+3.*logthetay).*(exp(1).^(2.*logthetax+logthetay).*( ...
  3.*exp(1).^logthetay+(-2).*(y+(-1).*yp).^2)+dt.*(6.*exp(1).^(a2+2.* ...
  logthetax+logthetay)+3.*exp(1).^(a2+logthetax+2.*logthetay)+3.*a1.*exp( ...
  1).^(logthetax+2.*logthetay).*ubarp.*(x+(-1).*xp)+(-3).*exp(1).^(a2+2.* ...
  logthetay).*(x+(-1).*xp).^2+6.*a1.*exp(1).^(2.*logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-6).*exp(1).^(a2+2.*logthetax).*(y+(-1).*yp).^2+( ...
  -2).*exp(1).^(a2+logthetax+logthetay).*(y+(-1).*yp).^2+(-2).*a1.*exp(1) ...
  .^(logthetax+logthetay).*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+2.*exp(1) ...
  .^(a2+logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+(-1).*a1.*exp(1).^( ...
  2.*logthetax).*vbarp.*(y+(-1).*yp).^3))+a1.*dt.*exp(1).^((-3).* ...
  logthetax+2.*logthetay).*(exp(1).^(2.*logthetax+logthetay).*((-4).*exp( ...
  1).^(2.*logthetay).*ubar.*(x+(-1).*xp)+(-9).*exp(1).^(logthetax+ ...
  logthetay).*vbar.*(y+(-1).*yp)+3.*exp(1).^logthetay.*ubar.*(x+(-1).*xp) ...
  .*(y+(-1).*yp).^2+2.*exp(1).^logthetax.*vbar.*(y+(-1).*yp).^3)+a1.*dt.* ...
  exp(1).^logthetax.*(4.*exp(1).^(logthetax+3.*logthetay).*ubar.*ubarp+9.* ...
  exp(1).^(2.*(logthetax+logthetay)).*vbar.*vbarp+(-4).*exp(1).^(3.* ...
  logthetay).*ubar.*ubarp.*(x+(-1).*xp).^2+(-3).*exp(1).^(logthetax+2.* ...
  logthetay).*(3.*ubar.*vbarp.*(x+(-1).*xp)+ubarp.*(3.*vbar.*(x+(-1).*xp)+ ...
  ubar.*(y+(-1).*yp))).*(y+(-1).*yp)+(-12).*exp(1).^(2.*logthetax+ ...
  logthetay).*vbar.*vbarp.*(y+(-1).*yp).^2+3.*exp(1).^(2.*logthetay).* ...
  ubar.*ubarp.*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+2.*exp(1).^(logthetax+ ...
  logthetay).*(ubarp.*vbar+ubar.*vbarp).*(x+(-1).*xp).*(y+(-1).*yp).^3+ ...
  exp(1).^(2.*logthetax).*vbar.*vbarp.*(y+(-1).*yp).^4)+dt.*exp(1).^a2.*(( ...
  -12).*exp(1).^(logthetax+3.*logthetay).*ubar.*(x+(-1).*xp)+4.*exp(1).^( ...
  3.*logthetay).*ubar.*(x+(-1).*xp).^3+(-9).*exp(1).^(2.*(logthetax+ ...
  logthetay)).*(ubar.*(x+(-1).*xp)+vbar.*(y+(-1).*yp))+(-30).*exp(1).^(3.* ...
  logthetax+logthetay).*vbar.*(y+(-1).*yp)+9.*exp(1).^(logthetax+2.* ...
  logthetay).*(x+(-1).*xp).*(vbar.*(x+(-1).*xp)+ubar.*(y+(-1).*yp)).*(y+( ...
  -1).*yp)+(-3).*exp(1).^(2.*logthetay).*ubar.*(x+(-1).*xp).^3.*(y+(-1).* ...
  yp).^2+2.*exp(1).^(2.*logthetax+logthetay).*(6.*ubar.*(x+(-1).*xp)+ ...
  vbar.*(y+(-1).*yp)).*(y+(-1).*yp).^2+10.*exp(1).^(3.*logthetax).*vbar.*( ...
  y+(-1).*yp).^3+(-2).*exp(1).^(logthetax+logthetay).*vbar.*(x+(-1).*xp) ...
  .^2.*(y+(-1).*yp).^3+(-1).*exp(1).^(2.*logthetax).*ubar.*(x+(-1).*xp).*( ...
  y+(-1).*yp).^4))+dt.*exp(1).^(a2+(-4).*logthetax+logthetay).*(exp(1).^( ...
  2.*logthetax+logthetay).*(12.*exp(1).^(2.*(logthetax+logthetay))+5.*exp( ...
  1).^(logthetax+3.*logthetay)+(-5).*exp(1).^(3.*logthetay).*(x+(-1).*xp) ...
  .^2+(-18).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-4).*exp( ...
  1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+4.*exp(1).^(2.*logthetay) ...
  .*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+2.*exp(1).^(2.*logthetax).*(y+(-1).* ...
  yp).^4)+dt.*(45.*exp(1).^(a2+4.*logthetax+2.*logthetay)+15.*exp(1).^(a2+ ...
  2.*logthetax+4.*logthetay)+24.*exp(1).^(a2+3.*(logthetax+logthetay))+ ...
  15.*a1.*exp(1).^(2.*logthetax+4.*logthetay).*ubarp.*(x+(-1).*xp)+(-30).* ...
  exp(1).^(a2+logthetax+4.*logthetay).*(x+(-1).*xp).^2+(-5).*a1.*exp(1).^( ...
  logthetax+4.*logthetay).*ubarp.*(x+(-1).*xp).^3+5.*exp(1).^(a2+4.* ...
  logthetay).*(x+(-1).*xp).^4+12.*a1.*exp(1).^(3.*(logthetax+logthetay)).* ...
  (ubarp.*(x+(-1).*xp)+vbarp.*(y+(-1).*yp))+(-12).*exp(1).^(a2+2.* ...
  logthetax+3.*logthetay).*(2.*x.^2+(-4).*x.*xp+2.*xp.^2+(y+(-1).*yp).^2)+ ...
  45.*a1.*exp(1).^(4.*logthetax+2.*logthetay).*vbarp.*(y+(-1).*yp)+(-12).* ...
  a1.*exp(1).^(2.*logthetax+3.*logthetay).*(x+(-1).*xp).*(vbarp.*(x+(-1).* ...
  xp)+ubarp.*(y+(-1).*yp)).*(y+(-1).*yp)+(-90).*exp(1).^(a2+4.*logthetax+ ...
  logthetay).*(y+(-1).*yp).^2+(-36).*exp(1).^(a2+3.*logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+24.*exp(1).^(a2+logthetax+3.*logthetay).*(x+ ...
  (-1).*xp).^2.*(y+(-1).*yp).^2+36.*exp(1).^(a2+2.*(logthetax+logthetay)) ...
  .*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+4.*a1.*exp(1).^(logthetax+3.* ...
  logthetay).*ubarp.*(x+(-1).*xp).^3.*(y+(-1).*yp).^2+(-4).*exp(1).^(a2+ ...
  3.*logthetay).*(x+(-1).*xp).^4.*(y+(-1).*yp).^2+(-3).*a1.*exp(1).^(3.* ...
  logthetax+2.*logthetay).*(6.*ubarp.*(x+(-1).*xp)+vbarp.*(y+(-1).*yp)).*( ...
  y+(-1).*yp).^2+(-20).*a1.*exp(1).^(4.*logthetax+logthetay).*vbarp.*(y+( ...
  -1).*yp).^3+3.*a1.*exp(1).^(2.*(logthetax+logthetay)).*vbarp.*(x+(-1).* ...
  xp).^2.*(y+(-1).*yp).^3+15.*exp(1).^(a2+4.*logthetax).*(y+(-1).*yp).^4+ ...
  4.*exp(1).^(a2+3.*logthetax+logthetay).*(y+(-1).*yp).^4+2.*a1.*exp(1).^( ...
  3.*logthetax+logthetay).*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^4+(-4).*exp( ...
  1).^(a2+2.*logthetax+logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^4+a1.* ...
  exp(1).^(4.*logthetax).*vbarp.*(y+(-1).*yp).^5))+(-1).*dt.*exp(1).^a2.*( ...
  (-2).*exp(1).^((-2).*logthetax+2.*logthetay).*(exp(1).^(2.*(logthetax+ ...
  logthetay))+dt.*(6.*exp(1).^(a2+2.*logthetax+logthetay)+exp(1).^(a2+ ...
  logthetax+2.*logthetay)+a1.*exp(1).^(logthetax+2.*logthetay).*ubarp.*(x+ ...
  (-1).*xp)+(-1).*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).^2+3.*a1.*exp(1) ...
  .^(2.*logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-6).*exp(1).^(a2+2.* ...
  logthetax).*(y+(-1).*yp).^2))+exp(1).^((-2).*logthetax+2.*logthetay).*(( ...
  -1).*exp(1).^logthetax+(x+(-1).*xp).^2).*(exp(1).^(2.*logthetay).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4))+(-2).*dt.*exp(1).^(a2+(-2).*logthetax+4.* ...
  logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+(-2).*dt.*exp(1) ...
  .^((-3).*logthetax+4.*logthetay).*(a1.*exp(1).^logthetax.*ubarp+(-2).* ...
  exp(1).^a2.*(x+(-1).*xp)).*(x+(-1).*xp).*(exp(1).^logthetay+(-1).*(y+( ...
  -1).*yp).^2)+(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).* ...
  yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^(2.* ...
  logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*(( ...
  -1).*exp(1).^logthetay+(y+(-1).*yp).^2)+(-2).*exp(1).^logthetay.*(a1.* ...
  dt.*exp(1).^((-1).*logthetax+logthetay).*(3.*exp(1).^(logthetax+ ...
  logthetay).*vbarp+(-2).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1) ...
  .*yp)+(-3).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^2)+(-2).*exp(1).^( ...
  2.*logthetay).*(y+(-1).*yp)+2.*dt.*exp(1).^(a2+(-2).*logthetax).*((-6).* ...
  exp(1).^(2.*logthetax+logthetay)+(-1).*exp(1).^(logthetax+2.*logthetay)+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+2.*exp(1).^(2.*logthetax).*(y+( ...
  -1).*yp).^2).*(y+(-1).*yp)).*(y+(-1).*yp)).*((-6)+(1/2).*exp(1).^((-1).* ...
  logthetay).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^logthetay.*(exp(1).^((-2).* ...
  logthetax+logthetay).*ubar.*((-1).*exp(1).^logthetax.*(x+(-1).*xp).*( ...
  exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.* ...
  dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^(2.*logthetay).* ...
  ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).* ...
  yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).* ...
  exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).* ...
  logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+ ...
  3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1) ...
  .^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+a1.*dt.*exp(1).^( ...
  logthetax+2.*logthetay).*ubarp.*(exp(1).^logthetay+(-1).*(y+(-1).*yp) ...
  .^2)+(-2).*dt.*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2))+vbar.*(exp(1).^logthetay.*(a1.*dt.* ...
  exp(1).^((-1).*logthetax+logthetay).*(3.*exp(1).^(logthetax+logthetay).* ...
  vbarp+(-2).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).* ...
  exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^2)+(-2).*exp(1).^(2.*logthetay) ...
  .*(y+(-1).*yp)+2.*dt.*exp(1).^(a2+(-2).*logthetax).*((-6).*exp(1).^(2.* ...
  logthetax+logthetay)+(-1).*exp(1).^(logthetax+2.*logthetay)+exp(1).^(2.* ...
  logthetay).*(x+(-1).*xp).^2+2.*exp(1).^(2.*logthetax).*(y+(-1).*yp).^2) ...
  .*(y+(-1).*yp))+(-1).*(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1).* ...
  (y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1).^( ...
  2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*(y+ ...
  (-1).*yp))).*((-5)+(1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2)+ ...
  exp(1).^(2.*logthetay).*(exp(1).^(2.*logthetay).*(exp(1).^logthetay+(-1) ...
  .*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+logthetay).*(exp(1) ...
  .^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1) ...
  .^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*(( ...
  -4)+(1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2));

case 4 % a1

K=dt.*exp(1).^(logsigma+(-6).*logthetay+(-1/2).*exp(1).^((-1).*logthetax) ...
  .*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2).*( ...
  exp(1).^logthetay.*(exp(1).^((-2).*logthetax+logthetay).*ubar.*((-1).* ...
  exp(1).^logthetax.*(x+(-1).*xp).*(exp(1).^(2.*logthetay).*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).*logthetax+ ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*(3.*exp(1).^(2.*( ...
  logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^( ...
  3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*( ...
  y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+ ...
  exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4))+a1.*dt.*exp(1).^(logthetax+2.*logthetay).* ...
  ubarp.*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+(-2).*dt.*exp(1).^(a2+ ...
  2.*logthetay).*(x+(-1).*xp).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2))+ ...
  vbar.*(exp(1).^logthetay.*(a1.*dt.*exp(1).^((-1).*logthetax+logthetay).* ...
  (3.*exp(1).^(logthetax+logthetay).*vbarp+(-2).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).*exp(1).^logthetax.*vbarp.*(y+( ...
  -1).*yp).^2)+(-2).*exp(1).^(2.*logthetay).*(y+(-1).*yp)+2.*dt.*exp(1).^( ...
  a2+(-2).*logthetax).*((-6).*exp(1).^(2.*logthetax+logthetay)+(-1).*exp( ...
  1).^(logthetax+2.*logthetay)+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+2.* ...
  exp(1).^(2.*logthetax).*(y+(-1).*yp).^2).*(y+(-1).*yp))+(-1).*(exp(1).^( ...
  2.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1) ...
  .^((-1).*logthetax+logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).* ...
  xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1) ...
  .^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*( ...
  3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+( ...
  -1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.* ...
  logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*(y+(-1).*yp)))+ ...
  a1.*dt.*exp(1).^((-2).*logthetax+2.*logthetay).*(exp(1).^logthetay.* ...
  ubar.*(exp(1).^(logthetax+logthetay).*ubarp.*(exp(1).^logthetay+(-1).*( ...
  y+(-1).*yp).^2)+(-1).*(x+(-1).*xp).*(exp(1).^(2.*logthetay).*ubarp.*(x+( ...
  -1).*xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).* ...
  exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^3))+exp(1).^logthetax.*vbar.*(3.*exp( ...
  1).^(logthetax+2.*logthetay).*vbarp+(-3).*exp(1).^(2.*logthetay).* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-6).*exp(1).^(logthetax+logthetay).* ...
  vbarp.*(y+(-1).*yp).^2+exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).* ...
  yp).^3+exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^4))+(-1).*dt.*exp(1).^( ...
  a2+(-3).*logthetax+logthetay).*((-2).*exp(1).^(2.*(logthetax+logthetay)) ...
  .*(exp(1).^logthetay.*ubarp.*(x+(-1).*xp)+3.*exp(1).^logthetax.*vbarp.*( ...
  y+(-1).*yp))+(-2).*exp(1).^(logthetax+3.*logthetay).*ubarp.*(x+(-1).*xp) ...
  .*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+exp(1).^(2.*logthetay).*(( ...
  -1).*exp(1).^logthetax+(x+(-1).*xp).^2).*(exp(1).^(2.*logthetay).* ...
  ubarp.*(x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).* ...
  yp)+(-1).*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).* ...
  exp(1).^logthetax.*vbarp.*(y+(-1).*yp).^3)+exp(1).^(2.*logthetax).*((-1) ...
  .*exp(1).^logthetay+(y+(-1).*yp).^2).*(exp(1).^(2.*logthetay).*ubarp.*( ...
  x+(-1).*xp)+3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).* ...
  exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^3)+(-2).*exp(1).^(2.*logthetax+ ...
  logthetay).*(3.*exp(1).^(logthetax+logthetay).*vbarp+(-2).*exp(1) ...
  .^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).*exp(1).^logthetax.* ...
  vbarp.*(y+(-1).*yp).^2).*(y+(-1).*yp))+exp(1).^((-1).*logthetax+3.* ...
  logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.*exp(1).^( ...
  logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1).^logthetay.* ...
  ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1).^logthetax.*vbarp.*(y+ ...
  (-1).*yp).^3));

case 5 % a2

K=dt.*exp(1).^(logsigma+(-6).*logthetay+(-1/2).*exp(1).^((-1).*logthetax) ...
  .*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1).*yp).^2).*( ...
  exp(1).^a2.*(2.*exp(1).^((-2).*logthetax+2.*logthetay).*(exp(1).^(2.*( ...
  logthetax+logthetay))+dt.*(6.*exp(1).^(a2+2.*logthetax+logthetay)+exp(1) ...
  .^(a2+logthetax+2.*logthetay)+a1.*exp(1).^(logthetax+2.*logthetay).* ...
  ubarp.*(x+(-1).*xp)+(-1).*exp(1).^(a2+2.*logthetay).*(x+(-1).*xp).^2+3.* ...
  a1.*exp(1).^(2.*logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-6).*exp(1) ...
  .^(a2+2.*logthetax).*(y+(-1).*yp).^2))+(-1).*exp(1).^((-2).*logthetax+ ...
  2.*logthetay).*((-1).*exp(1).^logthetax+(x+(-1).*xp).^2).*(exp(1).^(2.* ...
  logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^(( ...
  -1).*logthetax+logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+ ...
  3.*exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1) ...
  .^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*( ...
  3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+( ...
  -1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.* ...
  logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+2.*dt.*exp(1).^( ...
  a2+(-2).*logthetax+4.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp) ...
  .^2)+2.*dt.*exp(1).^((-3).*logthetax+4.*logthetay).*(a1.*exp(1) ...
  .^logthetax.*ubarp+(-2).*exp(1).^a2.*(x+(-1).*xp)).*(x+(-1).*xp).*(exp( ...
  1).^logthetay+(-1).*(y+(-1).*yp).^2)+(-1).*(exp(1).^(2.*logthetay).*( ...
  exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+a1.*dt.*exp(1).^((-1).* ...
  logthetax+logthetay).*(exp(1).^(2.*logthetay).*ubarp.*(x+(-1).*xp)+3.* ...
  exp(1).^(logthetax+logthetay).*vbarp.*(y+(-1).*yp)+(-1).*exp(1) ...
  .^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp).^2+(-1).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^3)+dt.*exp(1).^(a2+(-2).*logthetax).*( ...
  3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+3.*logthetay)+( ...
  -1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.* ...
  logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)).*((-1).*exp(1) ...
  .^logthetay+(y+(-1).*yp).^2)+2.*exp(1).^logthetay.*(a1.*dt.*exp(1).^(( ...
  -1).*logthetax+logthetay).*(3.*exp(1).^(logthetax+logthetay).*vbarp+(-2) ...
  .*exp(1).^logthetay.*ubarp.*(x+(-1).*xp).*(y+(-1).*yp)+(-3).*exp(1) ...
  .^logthetax.*vbarp.*(y+(-1).*yp).^2)+(-2).*exp(1).^(2.*logthetay).*(y+( ...
  -1).*yp)+2.*dt.*exp(1).^(a2+(-2).*logthetax).*((-6).*exp(1).^(2.* ...
  logthetax+logthetay)+(-1).*exp(1).^(logthetax+2.*logthetay)+exp(1).^(2.* ...
  logthetay).*(x+(-1).*xp).^2+2.*exp(1).^(2.*logthetax).*(y+(-1).*yp).^2) ...
  .*(y+(-1).*yp)).*(y+(-1).*yp))+a1.*dt.*exp(1).^(a2+(-2).*logthetax+ ...
  logthetay).*(exp(1).^logthetay.*ubar.*(x+(-1).*xp).*((-2).*exp(1).^(2.* ...
  logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+(-1).*exp(1).^(( ...
  -1).*logthetax).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^( ...
  logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+( ...
  -6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^( ...
  logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1) ...
  .*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4))+ ...
  vbar.*((-15).*exp(1).^(2.*(logthetax+logthetay))+(-3).*exp(1).^( ...
  logthetax+3.*logthetay)+3.*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+10.* ...
  exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(2.*logthetay).*(x+(-1).*xp) ...
  .^2.*(y+(-1).*yp).^2+(-1).*exp(1).^(2.*logthetax).*(y+(-1).*yp).^4).*(y+ ...
  (-1).*yp))+(-1).*dt.*exp(1).^(a2+(-4).*logthetax).*((-2).*exp(1).^(a2+ ...
  2.*logthetax+4.*logthetay).*(exp(1).^logthetay+(-1).*(y+(-1).*yp).^2)+ ...
  4.*exp(1).^(a2+logthetax+4.*logthetay).*(x+(-1).*xp).^2.*(exp(1) ...
  .^logthetay+(-1).*(y+(-1).*yp).^2)+(-2).*exp(1).^(a2+2.*(logthetax+ ...
  logthetay)).*(6.*exp(1).^(2.*logthetax+logthetay)+exp(1).^(logthetax+2.* ...
  logthetay)+(-1).*exp(1).^(2.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^( ...
  2.*logthetax).*(y+(-1).*yp).^2)+exp(1).^(a2+2.*logthetay).*((-1).*exp(1) ...
  .^logthetax+(x+(-1).*xp).^2).*(3.*exp(1).^(2.*(logthetax+logthetay))+ ...
  exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).* ...
  xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).* ...
  exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay) ...
  .*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp) ...
  .^4)+exp(1).^(a2+2.*logthetax).*((-1).*exp(1).^logthetay+(y+(-1).*yp) ...
  .^2).*(3.*exp(1).^(2.*(logthetax+logthetay))+exp(1).^(logthetax+3.* ...
  logthetay)+(-1).*exp(1).^(3.*logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^( ...
  2.*logthetax+logthetay).*(y+(-1).*yp).^2+(-1).*exp(1).^(logthetax+2.* ...
  logthetay).*(y+(-1).*yp).^2+exp(1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+ ...
  (-1).*yp).^2+exp(1).^(2.*logthetax).*(y+(-1).*yp).^4)+(-4).*exp(1).^(a2+ ...
  2.*logthetax+logthetay).*((-6).*exp(1).^(2.*logthetax+logthetay)+(-1).* ...
  exp(1).^(logthetax+2.*logthetay)+exp(1).^(2.*logthetay).*(x+(-1).*xp) ...
  .^2+2.*exp(1).^(2.*logthetax).*(y+(-1).*yp).^2).*(y+(-1).*yp).^2)+exp(1) ...
  .^(a2+(-2).*logthetax+2.*logthetay).*(3.*exp(1).^(2.*(logthetax+ ...
  logthetay))+exp(1).^(logthetax+3.*logthetay)+(-1).*exp(1).^(3.* ...
  logthetay).*(x+(-1).*xp).^2+(-6).*exp(1).^(2.*logthetax+logthetay).*(y+( ...
  -1).*yp).^2+(-1).*exp(1).^(logthetax+2.*logthetay).*(y+(-1).*yp).^2+exp( ...
  1).^(2.*logthetay).*(x+(-1).*xp).^2.*(y+(-1).*yp).^2+exp(1).^(2.* ...
  logthetax).*(y+(-1).*yp).^4));

case 6 % logsigmap

K=dt.^2.*exp(1).^(logsigmap+(-2).*logthetaxp+(-1/2).*exp(1).^((-1).* ...
  logthetaxp).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetayp).*(y+(-1) ...
  .*yp).^2).*(exp(1).^logthetaxp+(-1).*(x+(-1).*xp).^2);

case 7 % logthetaxp

K=(-1/2).*dt.^2.*exp(1).^(logsigmap+(-3).*logthetaxp+(-1/2).*exp(1).^((-1) ...
  .*logthetaxp).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetayp).*(y+( ...
  -1).*yp).^2).*(2.*exp(1).^(2.*logthetaxp)+(-5).*exp(1).^logthetaxp.*(x+( ...
  -1).*xp).^2+(x+(-1).*xp).^4);

case 8 % logthetayp

K=(1/2).*dt.^2.*exp(1).^(logsigmap+(-2).*logthetaxp+(-1).*logthetayp+( ...
  -1/2).*exp(1).^((-1).*logthetaxp).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1) ...
  .*logthetayp).*(y+(-1).*yp).^2).*(exp(1).^logthetaxp+(-1).*(x+(-1).*xp) ...
  .^2).*(y+(-1).*yp).^2;

otherwise

        K = zeros(n_x, n_xp);
end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function handleFileName(obj)

%--IF EMPTY FILENAME, CHECK FOR PLOT TITLES--%
try
    if isempty(obj.PlotOptions.FileName)
        for t = 1:obj.State.Figure.NumTexts
            if obj.State.Text(t).Title
                str = get(obj.State.Text(t).Handle,'String');
                interp = get(obj.State.Text(t).Handle,'Interpreter');
                obj.PlotOptions.FileName = parseString(str,interp);

                % untitle.html if \text exist (special chars)
                if ~isempty(strfind(obj.PlotOptions.FileName, '\text'))
                    obj.PlotOptions.FileName = 'untitled';
                end
            end
        end
    end
catch
    obj.PlotOptions.FileName = 'untitled';
end

%--IF FILENAME IS STILL EMPTY SET TO UNTITLED--%
if isempty(obj.PlotOptions.FileName)
    obj.PlotOptions.FileName = 'untitled';
end

end

author: plotly


text:
function [x, infos] = kl_bmd_nmf(V, rank, in_options)
% Block mirror descent method for KL-based non-negative matrix factorization (KL-BMD-NMF).
%
% The problem of interest is defined as
%
%           min f(V, W, H),
%           where
%           {V, W, H} >= 0.
%
% Given a non-negative matrix V, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       V           : (m x n) non-negative matrix to factorize
%       rank        : rank
%       in_options
%
% Output:
%       x           : non-negative matrix solution, i.e., x.W: (m x rank), x.H: (rank x n)
%       infos       : log information
%           epoch   : iteration nuber
%           cost    : objective function value
%           optgap  : optimality gap
%           time    : elapsed time
%           grad_calc_count : number of sampled data elements (gradient calculations)
%
%
% This file is part of NMFLibrary
%
% This file has been ported from
% BMD.m at https://github.com/LeThiKhanhHien/KLNMF written originally by LTK Hien.
%
% Ported by H.Kasai on June 28, 2022
%
% Change log:
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.metric.type   = 'kl-div';
    local_options.myeps         = 1e-16;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = 'KL-BMD';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % initialize for this algorithm
    lambdaH = (1./(sum(V))); % the row which is the sum of columns of X
    lambdaH = repmat(lambdaH, rank, 1);
    lambdaW = (1./(sum(V,2)))';
    lambdaW = repmat(lambdaW, rank , 1);

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('KL-BMD: Epoch = 0000, cost = %.16e, optgap = %.4e\n', f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update H
        rj = sum(W)';
        rjc = repmat(rj, 1, n);
        bAv = V./(W*H+eps);
        cj = W' * bAv;
        H = H ./ (1 + (lambdaH .* H) .* (rjc - cj));
        H = H + (H<options.myeps) .* options.myeps;

        % update W
        rj = sum(H, 2);
        rjc = repmat(rj, 1, m);
        bAv = V' ./ (H' * W' + eps);
        cj = H * bAv;
        Wt = W' ./ (1 + (lambdaW .* W') .* (rjc - cj));
        W = Wt';
        W = W + (W<options.myeps) .* options.myeps;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
%InterfaceUtil
% Checks if a class (or class instance) properly implements a specificed
% interface.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 9/12/2010
classdef InterfaceUtil_Test < TestCase
    methods
        function this = InterfaceUtil_Test(name)
            this = this@TestCase(name);
        end
    end

    %tests
    methods
        function testAssertInterface(~)
            import edu.stanford.covert.util.InterfaceUtil;

            %ex 1-4
            InterfaceUtil.assertInterface('double', 'double');
            InterfaceUtil.assertInterface('double', metaclass(zeros(0,0,'double')));
            InterfaceUtil.assertInterface('double', meta.class.fromName('double'));
            InterfaceUtil.assertInterface(zeros(1,1,'double'), 'double');

            %ex 5
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            InterfaceUtil.assertInterface('double', metadata);

            %ex 6
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            metadata.Methods{1}.Name2 = 'plus';
            assertExceptionThrown(@() InterfaceUtil.assertInterface('double', metadata), 'Interface:invalidInput');

            %ex 7
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Properties{1}.Name = 'data';
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            assertExceptionThrown(@() InterfaceUtil.assertInterface('double', metadata), 'Interface:incompleteInterface');

            %ex 8
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            metadata.Methods{1}.Access = 'protected';
            assertExceptionThrown(@() InterfaceUtil.assertInterface('double', metadata), 'Interface:incompatibleInterface');

            %ex 9
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            metadata.Methods{1}.Access = 'protected';
            assertExceptionThrown(@() InterfaceUtil.assertInterface('double', metadata), 'Interface:incompatibleInterface');

            %ex 10
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            metadata.Methods{1}.Access = {'protected','public'};
            InterfaceUtil.assertInterface('double', metadata);

            %ex 11
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            metadata.Methods{1}.Abstract = 1;
            assertExceptionThrown(@() InterfaceUtil.assertInterface('double', metadata), 'Interface:incompatibleInterface');

            %ex 12
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            metadata.Methods{1}.Abstract = [0 1];
            InterfaceUtil.assertInterface('double', metadata);

            %ex 13
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'plus';
            metadata.Methods{1}.Abstract = [];
            InterfaceUtil.assertInterface('double', metadata);

            %ex 14
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'ind2sub';
            metadata.Methods{1}.InputNames = 3;
            InterfaceUtil.assertInterface('edu.stanford.covert.util.SparseMat', metadata);

            %ex 15
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'ind2sub';
            metadata.Methods{1}.InputNames = 4;
            assertExceptionThrown(@() InterfaceUtil.assertInterface('edu.stanford.covert.util.SparseMat', metadata), 'Interface:incompatibleInterface');

            %ex 16
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Methods = cell(0,1);
            metadata.Methods{1} = struct;
            metadata.Methods{1}.Name = 'ind2sub';
            metadata.Methods{1}.InputNames = {'this','inds','siz'};
            InterfaceUtil.assertInterface('edu.stanford.covert.util.SparseMat', metadata);

            %ex 17
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Properties{1}.Name = 'inds';
            metadata.Properties{1}.GetMethod = 1;
            metadata.Methods = cell(0,1);
            InterfaceUtil.assertInterface('edu.stanford.covert.util.SparseMat', metadata);

            %ex 18
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Properties{1}.Name = 'inds';
            metadata.Properties{1}.GetMethod = 0;
            metadata.Methods = cell(0,1);
            assertExceptionThrown(@() InterfaceUtil.assertInterface('edu.stanford.covert.util.SparseMat', metadata), 'Interface:incompatibleInterface');

            %ex 19
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Properties{1}.Name = 'inds';
            metadata.Properties{1}.GetMethod = [0 1];
            metadata.Methods = cell(0,1);
            InterfaceUtil.assertInterface('edu.stanford.covert.util.SparseMat', metadata);

            %ex 20
            metadata = struct;
            metadata.Properties = cell(0,1);
            metadata.Properties{1}.Name = 'inds';
            metadata.Properties{1}.GetMethod = [];
            metadata.Methods = cell(0,1);
            InterfaceUtil.assertInterface('edu.stanford.covert.util.SparseMat', metadata);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function cvx_optval = huber_pos( x, M, t ) %#ok

%HUBER   Internal cvx version.

%
% Check arguments
%

error( nargchk( 1, 3, nargin ) );
if ~cvx_isconvex( x ),
    error( 'Disciplined convex programming error:\n    HUBER_POS is convex and nondecreasing in X, so X must be convex.', 1 ); %#ok
end
if nargin < 2,
    M = 1;
elseif isa( M, 'cvx' ),
    error( 'Second argument must be numeric.' );
elseif ~isreal( M ) || any( M( : ) <= 0 ),
    error( 'Second argument must be real and positive.' );
end
if nargin < 3,
    t = 1;
elseif ~isreal( t ),
    error( 'Third argument must be real.' );
elseif cvx_isconstant( t ) && nnz( cvx_constant( t ) <= 0 ),
    error( 'Third argument must be real and positive.' );
elseif ~cvx_isconcave( t ),
    error( 'Disciplined convex programming error:\n    HUBER_POS is convex and nonincreasing in T, so T must be concave.', 1 ); %#ok
end
sz = cvx_size_check( x, M, t );
if isempty( sz ),
    error( 'Sizes are incompatible.' );
end

%
% Compute result
%

cvx_begin separable
    variables v( sz ) w( sz )
    minimize( quad_over_lin( w, t, 0 ) + 2 .* M .* v )
    x <= w + v; %#ok
    w <= M * t; %#ok
    v >= 0; %#ok
cvx_end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: SSD versus translations, linearInter, HNSP level=4
%
%==============================================================================

clear, close all, help(mfilename)

setup2DHNSPData;
level = 4; m = ML{level}.m;
imgModel('set','imgModel','linearInter');
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
xc    = getCellCenteredGrid(omega,m);
Rc    = imgModel(R,omega,xc);

trafo('set','trafo','translation2D');
figure(1); clf;
[w1,w2] = ndgrid(0.2*linspace(-1,1,21),0.2*linspace(-1,1,21));
dc  = zeros(size(w1));
for j=1:numel(dc),
  yc = trafo([w1(j);w2(j)],xc);
  Tc = imgModel(T,omega,yc);
  dc(j) = SSD(Tc,Rc,omega,m);
  viewImage(Tc,omega,m); FAIRpause(1/100)
end;
figure(1); clf; surf(w1,w2,dc); hold on; grid off; contour(w1,w2,dc)
title(sprintf('translation, m=[%d,%d]',m)); view(-135,33);
%==============================================================================


text:
function u = acosh(a)
%ACOSH        slope inverse hyperbolic cosine acosh(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  u = a;

  u.r = acosh(a.r);
  u.s = slopeconvexconcave('acosh','1./sqrt(sqr(%)-1)',a,0);

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
%MDS Trainale mapping for multidimensional scaling, a variant of Sammon mapping
%
%   [W,J,STRESS] = MDS(DT,Y,OPTIONS)
%   [W,J,STRESS] = MDS(DT,K,OPTIONS)
%   [W,J,STRESS] = DT*MDS([],K,OPTIONS)
%   [W,J,STRESS] = DT*MDS(K,OPTIONS)
%              X = DS*W
%
%
% INPUT
%   DT      Square (M x M) dissimilarity matrix used for training
%   DS      N x M dissimilarity matrix between testset and trainset
%   Y       M x K matrix containing starting target configuration, or
%   K       Desired output dimensionality (default 2)
%   OPTIONS Various parameters of the  minimization procedure put into
%           a structure consisting of the following fields: 'q', 'optim',
%           'init','etol','maxiter', 'isratio', 'itmap', 'st' and 'inspect'
%           (default:
%
%              OPTIONS.q       = 0
%              OPTIONS.optim   = 'pn'
%              OPTIONS.init    = 'cs'
%              OPTIONS.etol    = 1e-6 (the precise value depends on q)
%              OPTIONS.maxiter = 100
%              OPTIONS.isratio = 0
%              OPTIONS.itmap   = 'yes'
%              OPTIONS.st      = 1
%              OPTIONS.inspect = 2).
%
% OUTPUT
%   W       Multidimensional scaling mapping
%   J       Index of points removed before optimization
%   STRESS  Vector of stress values
%   X       N x K dataset with output configuration
%
% DESCRIPTION
% Finds a nonlinear MDS map (a variant of the Sammon map) of objects
% represented by a symmetric distance matrix D with zero diagonal, given
% either the dimensionality K or the initial configuration Y. This is done
% in an iterative manner by minimizing the Sammon stress between the given
% dissimilarities D (DT or DS) and the distances DY in the K-dimensional
% target space:
%
%   E = 1/(sum_{i<j} D_{ij}^(Q+2)) sum_{i<j} (D_{ij} - DY_{ij})^2 * D_{ij}^Q
%
% If D(i,j) = 0 for any different points i and j, then one of them is
% superfluous. The indices of these points are returned in J.
%
% There is a simplified interface to MDS, called SAMMONM. The main
% differences are that SAMMONM operates on feature based datasets, while MDS
% expects dissimilarity matrices; MDS maps new objects by a second
% optimization procedures minimizing the stress for the test objects, while
% SAMMONM uses a linear mapping between dissimilarities and the target
% space. See also PREX_MDS for examples.
%
% OPTIONS is an optional variable, using which the parameters for the mapping
% can be controlled. It may contain the following fields:
%
%   Q        Stress measure to use (see above): -2,-1,0,1 or 2.
%   INIT     Initialisation method for Y: 'randp', 'randnp', 'maxv', 'cs'
%            or 'kl'. See MDS_INIT for an explanation.
%   OPTIM    Minimization procedure to use: 'pn' for Pseudo-Newton or
%            'scg' for Scaled Conjugate Gradients.
%   ETOL     Tolerance of the minimization procedure. Usually, it should be
%   MAXITER  in the order of 1e-6. If MAXITER is given (see below), then the
%            optimization is stopped either when the error drops below ETOL or
%            MAXITER iterations are reached.
%   ISRATIO  Indicates whether a ratio MDS should be performed (1) or not (0).
%            If ISRATIO is 1, then instead of fitting the dissimilarities
%            D_{ij}, A*D_{ij} is fitted in the stress function. The value A
%            is estimated analytically in each iteration.
%   ITMAP    Determines the way new points are mapped, either in an iterative
%            manner ('yes') by minimizing the stress; or by a linear projection
%            ('no').
%   ST       Status, determines whether after each iteration the stress should
%   INSPECT  be printed on the screen (1) or not (0). When INSPECT > 0,
%            ST = 1 and the mapping is onto 2D or larger, then the progress
%            is plotted during the minimization every INSPECT iterations.
%
% Important:
% - It is assumed that D either is or approximates a Euclidean distance
%     matrix, i.e:
%            D_{ij} = sqrt (sum_k(x_i - x_j)^2).
% - Missing values can be handled; they should be marked by 'NaN' in D.
%
% EXAMPLES:
% opt.optim = 'scg';
% opt.init  = 'cs';
% D  = sqrt(distm(a)); % Compute the Euclidean distance dataset of A
% w1 = mds(D,2,opt);   % An MDS map onto 2D initialized by Classical Scaling,
%                      % optimized by a Scaled Conjugate Gradients algorithm
% n  = size(D,1);
% y  = rand(n,2);
% w2 = mds(D,y,opt);   % An MDS map onto 2D initialized by random vectors
%
% z = rand(n,n);       % Set around 40% of the random distances to NaN, i.e.
% z = (z+z')/2;        % not used in the MDS mapping
% z = find(z <= 0.6);
% D(z) = NaN;
% D(1:n+1:n^2) = 0;    % Set the diagonal to zero
% opt.optim = 'pn';
% opt.init  = 'randnp';
% opt.etol  = 1e-8;    % Should be high, as only some distances are used
% w3 = mds(D,2,opt);   % An MDS map onto 2D initialized by a random projection
%
% REFERENCES
% 1. M.F. Moler, A Scaled Conjugate Gradient Algorithm for Fast Supervised
%    Learning', Neural Networks, vol. 6, 525-533, 1993.
% 2. W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
%    Numerical Recipes in C, Cambridge University Press, Cambridge, 1992.
% 3. I. Borg and P. Groenen, Modern Multidimensional Scaling, Springer
%    Verlag, Berlin, 1997.
% 4. T.F. Cox and M.A.A. Cox, Multidimensional Scaling, Chapman and Hall,
%    London, 1994.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, PREX_MDS, MDS_CS, SAMMONM, TSNEM

%
% Copyright: E. Pekalska, R.P.W. Duin, e.pekalska@37steps.com
%

function [w,J,err,opt,y] = mds(varargin)

  mapname = 'MDS';
	argin = shiftargin(varargin,{'scalar'});
  argin = setdefaults(argin,[],2,[]);

  if mapping_task(argin,'definition')
    w = define_mapping(argin,'untrained',mapname);

  else % Train, evaluate or extend a mapping.

    [D,y,options] = deal(argin{:});
    opt = mds_setopt(options);

  % YREP contains representative objects in the projected MDS space, i.e.
  % for which the mapping exists. YREP is empty for the original MDS, since
  %  no projection is available yet.

    yrep = [];

    if (isdataset(D) | isa(D,'double'))
      [m,mm] = size(D);

      % Convert D to double, but retain labels in LAB.
      if (isdataset(D)), lab = getlab(D); D = +D; else, lab = ones(m,1); end;

      if (ismapping(y))

        % The MDS mapping exists; this means that YREP has already been stored.
        pars = getdata(y); [k,c] = size(y);

        y     = [];        % Empty, should now be found.
        yrep  = pars{1};  % There exists an MDS map, hence YREP is stored.
        opt   = pars{2};  % Options used for the mapping.
        II    = pars{3};  % Index of non-repeated points in YREP.
        winit = pars{4};  % The Classical Scaling map, if INIT = 'cs'.
        v     = pars{5};  % Weights used for adding new points if ITMAP = 'no'.
        n = c;            % Number of dimensions of the projected space.

        % Initialization by 'cs' is not possible when there is no winit
        % (i.e. the CS map) and new points should be added.
        if (strcmp(opt.init,'cs')) & (isempty(winit))
          prwarning(2,'OPTIONS.init = cs is not possible when adding points; using kl.');
          opt.init = 'kl';
        end

        % If YREP is a scalar, we have an empty mapping.
        if (max(size(yrep)) == 1)
          y    = yrep;
          yrep = [];
        end

        % Check whether D is a matrix with the zero diagonal for the existing map.
        if (m == mm) & (length(intersect(find(D(:)<eps),1:m+1:(m*mm))) >= m)
          w = yrep;         % D is the same matrix as in the projection process;
          return            % YREP is then the solution
        end

        if (length(pars) < 6) | (isempty(pars{6}))
          yinit = [];
        else
          yinit = pars{6};   % Possible initial configuration for points to
                            % be added to an existing map
          if (size(yinit,1) ~= size(D,1))
            prwarning(2,'the size of the initial configuration does not match that of the dissimilarity matrix, using random initialization.')
            yinit =[];
          end
        end

      else

        % No MDS mapping available yet; perform the checks.

        if (~issym(D,1e-12))
          prwarning(2,'D is not a symmetric matrix; will average.');
          D = (D+D')/2;
        end

        % Check the number of zeros on the diagonal

        if (any(abs(diag(D)) > 0))
          error('D should have a zero diagonal');
        end
      end

    else    % D is neither a dataset nor a matrix of doubles
      error('D should be a dataset or a matrix of doubles.');
    end

    if (~isempty(y))

      % Y is the initial configuration or N, no MDS map exists yet;
      % D is a square matrix.

      % Remove identical points, i.e. points for which D(i,j) = 0 for i ~= j.
      % I contains the indices of points left for the MDS mapping, J those
      % of the removed points and P those of the points left in I which were
      % identical to those removed.

      [I,J,P] = mds_reppoints(D);
      D = D(I,I);
      [ni,nc] = size(D);

      % NANID is an extra field in the OPTIONS structure, containing the indices
      % of NaN values (= missing values) in distance matrix D.

      opt.nanid = find(isnan(D(:)) == 1);

      % Initialise Y.

      [m2,n] = size(y);
      if (max(m2,n) == 1)    % Y is a scalar, hence really is a dimensionality N.
        n = y;
        [y,winit] = mds_init(D,n,opt.init);
      else
        if (mm ~= m2)
          error('The matrix D and the starting configuration Y should have the same number of columns.');
        end
        winit = [];
        y = +y(I,:);
      end

      % The final number of true distances is:
      no_dist = (ni*(ni-1) - length(opt.nanid))/2;

    else

      % This happens only if we add extra points to an existing MDS map.
      % Remove identical points, i.e. points for which D(i,j) = 0 for i ~= j.
      % I contains the indices of points left for the MDS mapping, J those
      % of the removed points and P those of the points left in I which were
      % identical to those removed.

      [I,J,P] = mds_reppoints(D(:,II));
      D = D(I,II);
      [ni,nc] = size(D);
      yrep = yrep(II,:);
      n = size(yrep,2);

      % NANID is an extra field in the OPTIONS structure, containing the indices
      % of NaN values (= missing values) in distance matrix D.

      opt.nanid = find(isnan(D(:)));

      % Initialise Y. if the new points should be added in an iterative manner.

      [m2,n] = size(yrep);
      if (~isempty(yinit))             % An initial configuration exists.
        y = yinit;
      elseif (strcmp(opt.init, 'cs')) & (~isempty(winit))
        y = D*winit;
      else
        y = mds_init(D,n,opt.init);
      end

      if (~isempty(opt.nanid))         % Rescale.
        scale = (max(yrep)-min(yrep))./(max(y)-min(y));
        y = y .* repmat(scale,ni,1);
      end

      % The final number of true distances is:
      no_dist = (ni*nc - length(opt.nanid));
    end

    % Check whether there is enough data left.

    if (~isempty(opt.nanid))
      if (n*ni+2 > no_dist),
        error('There are too many missing distances: it is not possible to determine the MDS map.');
      end
      if (strcmp (opt.itmap,'no'))
        opt.itmap = 'yes';
        prwarning(1,'Due to the missing values, the projection can only be iterative. OPTIONS are changed appropriately.')
      end
    end

    if (opt.inspect > 0)
      opt.plotlab = lab(I,:);  % Labels to be used for plotting in MDS_SAMMAP.
    else
      opt.plotlab = [];
    end

    if (isempty(yrep)) | (~isempty(yrep) & strcmp(opt.itmap, 'yes'))

      % Either no MDS exists yet OR new points should be mapped in an iterative manner.

      printinfo(opt);
      [yy,err] = mds_sammap(D,y,yrep,opt);

      % Define the linear projection of distances.

      v = [];
      if (isempty(yrep)) & (isempty(opt.nanid))
        if (rank(D) < m)
          v = prpinv(D)*yy;
        else
          v = D \ yy;
        end
      end
    else
      % New points should be added by a linear projection of dissimilarity data.
      yy = D*v;
    end

    % Establish the projected configuration including the removed points.

    y = zeros(m,n); y(I,:) = +yy;
    if (~isempty(J))
      if (~isempty(yrep))
        y(J,:) = +yrep(II(P),:);
      else
        for k=length(J):-1:1        % J: indices of removed points.
          y(J(k),:) = y(P(k),:);    % P: indices of points equal to points in J.
        end
      end
    end

    % In the definition step: shift the obtained configuration such that the
    % mean lies at the origin.

    if (isempty(yrep))
      y = y - ones(length(y),1)*mean(y);
      y = prdataset(y,lab);
    else
      w = prdataset(y,lab);
      return;
    end

    % In the definition step: the mapping should be stored.

    opt = rmfield(opt,'nanid');   % These fields are to be used internally only;
    opt = rmfield(opt,'plotlab'); % not to be set up from outside
    w   = prmapping(mfilename,'trained',{y,opt,I,winit,v,[]},[],m,n);
    w   = setname(w,mapname);

  end

return

% **********************************************************************************
%                             Extra functions
% **********************************************************************************

% PRINTINFO(OPT)
%
% Prints progress information to the file handle given by OPT.ST.

function printinfo(opt)

  if opt.st < 1
    return
  end
  %fprintf(opt.st,'Sammon mapping, error function with the parameter q=%d\n',opt.q);

  switch (opt.optim)
    case 'pn',
      %fprintf(opt.st,'Minimization by Pseudo-Newton algorithm\n');
    case 'scg',
      %fprintf(opt.st,'Minimization by Scaled Conjugate Gradients algorithm\n');
    otherwise
      error(strcat('Possible initialization methods: pn (Pseudo-Newton), ',...
                   'or scg (Scaled Conjugate Gradients).'));
  end

return

% **********************************************************************************

% [I,J,P] = MDS_REPPOINTS(D)
%
% Finds the indices of repeated/left points. J contains the indices of
% repeated points in D. This means that for each j in J, there is a point
% k ~= j such that D(J(j),k) = 0. I contains the indices of the remaining
% points, and P those of the points in I that were identical to those in J.
% Directly used in MDS routine.

function [I,J,P] = mds_reppoints(D)

  epsilon = 1e-20;      % Differences smaller than this are assumed to be zero.

  [m,mm] = size(D);

  I = 1:m; J = []; P = [];

  if (m == mm) & (all(abs(D(1:m+1:end)) <= epsilon))
    K = intersect (find (triu(ones(m),1)), find(D < epsilon));
    if (~isempty(K))
      P = mod(K,m);
      J = fix(K./m) + 1;
      I(J) = [];
    end
  else
    [J,P] = find(D<=epsilon);
    I(J)  = [];
  end

return;

% **********************************************************************************

% MDS_SETOPT  Sets parameters for the MDS mapping
%
%   OPT = MDS_SETOPT(OPT_GIVEN)
%
% INPUT
%   OPT_GIVEN Parameters for the MDS mapping, described below (default: [])
%
% OUTPUT
%   OPT       Structure of chosen options; if OPT_GIVEN is empty, then
%               OPT.q       = 0
%               OPT.optim   = 'pn'
%               OPT.init    = 'cs'
%               OPT.etol    = 1e-6 (the precise value depends on q)
%               OPT.maxiter = 100
%               OPT.itmap   = 'yes'
%               OPT.isratio = 0
%               OPT.st      = 1
%               OPT.inspect = 2
%
% DESCRIPTION
% Parameters for the MDS mapping can be set or changed. OPT_GIVEN consists
% of the following fields: 'q', 'init', 'optim', 'etol', 'maxiter','itmap',
% 'isratio', 'st' and 'inspect'. OPTIONS can include all the fields or some
% of them only. The fields of OPT have some default values, which can be
% changed by the OPT_GIVEN field values. If OPT_GIVEN is empty, then OPT
% contains all default values. For a description of the fields, see MDS.

%
% Copyright: Elzbieta Pekalska, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function opt = mds_setopt (opt_given)

  opt.q       = 0;
  opt.init    = 'cs';
  opt.optim   = 'pn';
  opt.st      = 1;
  opt.itmap   = 'yes';
  opt.maxiter = 100;
  opt.inspect = 2;
  opt.etol    = inf;
  opt.isratio = 0;
  opt.nanid   = [];    % Here are some extra values; set up in the MDS routine.
  opt.plotlab = [];    % Not to be changed by the user.

  if (~isempty(opt_given))
    if (~isstruct(opt_given))
      error('OPTIONS should be a structure with at least one of the following fields: q, init, etol, optim, maxiter, itmap, isratio, st or inspect.');
    end
    fn = fieldnames(opt_given);
    if (~all(ismember(fn,fieldnames(opt))))
      error('Wrong field names; valid field names are: q, init, optim, etol, maxiter, itmap, isratio, st or inspect.')
    end
    for i = 1:length(fn)
      opt = setfield(opt,fn{i},getfield(opt_given,fn{i}));
    end
  end

  if (isempty(intersect(opt.q,-2:1:2)))
    error ('OPTIONS.q should be -2, -1, 0, 1 or 2.');
  end

  if (opt.maxiter < 2)
    error ('OPTIONS.iter should be at least 1.');
  end

  if (isinf(opt.etol))
    switch (opt.q)                 % Different defaults for different stresses.
      case -2,
        opt.etol = 1e-6;
      case {-1,0}
        opt.etol = 10*sqrt(eps);
      case {1,2}
        opt.etol = sqrt(eps);
      end
  elseif (opt.etol <= 0) | (opt.etol >= 0.1)
    error ('OPTIONS.etol should be positive and smaller than 0.1.');
  end

  if (~ismember(opt.optim, {'pn','scg'}))
    error('OPTIONS.optim should be pn or scg.');
  end

  if (~ismember(opt.itmap, {'yes','no'}))
    error('OPTIONS.itmap should be yes or no.');
  end

return

% MDS_SAMMAP Sammon iterative nonlinear mapping for MDS
%
%   [YMAP,ERR] = MDS_SAMMAP(D,Y,YREP,OPTIONS)
%
% INPUT
%   D       Square (M x M) dissimilarity matrix
%   Y       M x N matrix containing starting configuration, or
%   YREP    Configuration of the representation points
%   OPTIONS Various parameters of the  minimization procedure put into
%           a structure consisting of the following fields: 'q', 'optim',
%           'init','etol','maxiter', 'isratio', 'itmap', 'st' and 'inspect'
%           (default:
%             OPTIONS.q       = 0
%             OPTIONS.optim   = 'pn'
%             OPTIONS.init    = 'cs'
%             OPTIONS.etol    = 1e-6 (the precise value depends on q)
%             OPTIONS.maxiter = 100
%             OPTIONS.isratio = 0
%             OPTIONS.itmap   = 'yes'
%             OPTIONS.st      = 1
%             OPTIONS.inspect = 2).
%
% OUTPUT
%   YMAP     Mapped configuration
%   ERR     Sammon stress
%
% DESCRIPTION
% Maps the objects given by a symmetric distance matrix D (with a zero
% diagonal) onto, say, an N-dimensional configuration YMAP by an iterative
% minimization of a variant of the Sammon stress. The minimization starts
% from the initial configuration Y; see MDS_INIT.
%
% YREP is the Sammon configuration of the representation set. It is used
% when new points have to be projected. In other words, if D is an M x M
% symmetric distance matrix, then YREP is empty; if D is an M x N matrix,
% then YMAP is sought such that D can approximate the distances between YMAP
% and YREP.
%
% Missing values can be handled by marking them by NaN in D.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, MDS, MDS_CS, MDS_INIT, MDS_SETOPT

%
% Copyright: Elzbieta Pekalska, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [y,err] = mds_sammap(Ds,y,yrep,opt)

  if (nargin < 4)
    opt = [];                 % Will be filled by MDS_SETOPT, below.
  end
  if isempty(opt), opt = mds_setopt(opt); end
  if (nargin < 3)
    yrep = [];
  end

  % Extract labels and calculate distance matrix.

  [m,n] = size(y);
  if (~isempty(yrep))
    replab = getlab(yrep);
    yrep = +yrep;
    D = sqrt(distm(y,yrep));
  else
    D = sqrt(distm(y));
    yrep = [];
    replab = [];
  end

%   if (isempty(opt.plotlab))
%     opt.plotlab = ones(m,1);
%   end

  it   = 0;           % Iteration number.
  eold = inf;         % Previous stress (error).

  % Calculate initial stress.

  [e,a]= mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio); err = e;

  %if opt.st > 0, fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e); end

  tt = sprintf('Processing %i iterations: ',opt.maxiter);
  prwaitbar(opt.maxiter,tt)

  switch (opt.optim)
    case 'pn',        % Pseudo-Newton minimization.

      epr   = e;      % Previous values of E, EOLD for line search algorithm.
      eepr  = e;
      add   = 1e-4;   % Avoid G2 equal 0; see below.

      BETA   = 1e-3;   % Parameter for line search algorithm.
      lam1   = 0;
      lam2   = 1;      % LAM (the step factor) lies in (LAM1, LAM2).
      LMIN   = 1e-10;  % Minimum acceptable value for LAM in line search.
      lambest = 1e-4; % LAM = LAMBEST, if LAM was not found.

      % Loop until the error change falls below the tolerance or until
      % the maximum number of iterations is reached.

      while (abs(e-eold) >= opt.etol*(1 + abs(e))) & (it < opt.maxiter)

        % Plot progress if requested.

%        if (opt.st > 0) & (opt.inspect > 0) & (mod(it,opt.inspect) == 0)
%          if (it == 0), figure(1); clf; end
%          mds_plot(y,yrep,opt.plotlab,replab,e);
%        end

        yold = y; eold = e;

        % Calculate the stress.

        [e,a]      = mds_samstress (opt.q,yold,yrep,Ds,[],opt.nanid,opt.isratio);

        % Calculate gradients and pseudo-Newton update P.

        [g1,g2,cc] = mds_gradients (opt.q,yold,yrep,a*Ds,[],opt.nanid);
        p        = (g1/cc)./(abs(g2/cc) + add);
        slope    = g1(:)' * p(:);

        lam = 1;

        % Search for a suitable step LAM using a line search.

        while (1)

          % Take a step and calculate the delta error DE.

          y      = yold + lam .* p;
          [e,a]  = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);
          de     = e - eold;

          % Stop if the condition for a suitable lam is fulfilled.

          if (de < BETA * lam * slope)
            break;
          end

          % Try to find a suitable step LAM.

          if (lam ~= 1)
            r1 = de - lam*slope;
            r2 = epr - eepr - lam2*slope;
            aa = (2*de + r1/lam^2 - r2/lam2^2)/(lam2-lam1)^3;
            bb = ((-lam2*r1)/lam^2 +(lam*r2)/lam2^2)/(lam2-lam1)^2;

            if (abs(aa) <= eps)
              lamtmp = -0.5 * slope/bb;
            else
              lamtmp = (-bb + sqrt(max(0,(bb^2 - 3*aa*slope))))/(3*aa);
            end

            % Prevent LAM from becoming too large.

            if (lamtmp > 0.5 * lam), lamtmp = 0.5 * lam; end

          else
            lamtmp = -0.5 * slope/(e - eold - slope);
          end

          % Prevent LAM from becoming too small.

          lam2 = lam;
          lam  = max(lamtmp, 0.1*lam);
          epr  = e;
          eepr = eold;

          if (lam < LMIN)
            y = yold + lambest .* p;
            [e,a] = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);
            break;
          end
        end

        it  = it + 1;
        err = [err; e];

        prwaitbar(opt.maxiter,it,[tt num2str(it) ', error: ' num2str(e)]);
        %if opt.st > 0, fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e); end

      end

    case 'scg',              % Scaled Conjugate Gradient minimization.

      sigma0  = 1e-8;
      lambda  = 1e-8;       % Regularization parameter.
      lambda1 = 0;

      % Calculate initial stress and direction of decrease P.

      [e,a]   = mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio);
      g1      = mds_gradients(opt.q,y,yrep,a*Ds,D,opt.nanid);   % Gradient
      p       = -g1;                                            % Direction of decrease
      gnorm2  = g1(:)' * g1(:);
      success = 1;

      % Sanity check.

      if (gnorm2 < 1e-15)
        prwarning(2,['Gradient is nearly zero: ' gnorm2 ' (unlikely); initial configuration is the right one.']);
        return
      end

      % Loop until the error change falls below the tolerance or until
      % the maximum number of iterations is reached.

      while (abs(eold - e) > opt.etol * (1 + e)) & (it < opt.maxiter)

        g0     = g1;            % Previous gradient
        pnorm2 = p(:)' * p(:);
        eold   = e;

        % Plot progress if requested.

%        if (opt.inspect > 0) & (mod(it,opt.inspect) == 0)
%          if (it == 0), figure(1); clf; end
%          mds_plot(y,yrep,opt.plotlab,replab,e);
%        end

        if (success)
          sigma = sigma0/sqrt(pnorm2);      % SIGMA: a small step from y to yy
          yy    = y + sigma .*p;
          [e,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio);
          g2    = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid); % Gradient for yy
          s     = (g2-g1)/sigma;        % Approximation of Hessian*P directly, instead of computing
          delta = p(:)' * s(:);         % the Hessian, since only DELTA = P'*Hessian*P is in fact needed
        end

        % Regularize the Hessian indirectly to make it positive definite.
        % DELTA is now computed as P'* (Hessian + regularization * Identity) * P
        delta = delta + (lambda1 - lambda) * pnorm2;

        % Indicate if the Hessian is negative definite; the regularization above was too small.
        if (delta < 0)
          lambda1 = 2 * (lambda - delta/pnorm2);  % Now the Hessian will be positive definite
          delta   = -delta + lambda * pnorm2;     % This is obtained after plugging lambda1
          lambda  = lambda1;                      % into the formulation a few lines above.
        end

        mi = - p(:)' * g1(:);
        yy = y + (mi/delta) .*p;                  % mi/delta is a step size
        [ee,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio);

        % This minimization procedure is based on the second order approximation of the
        % stress function by using the gradient and the Hessian approximation. The Hessian
        % is regularized, but maybe not sufficiently. The ratio Dc (<=1) below indicates
        % a proper approximation, if Dc is close to 1.

        Dc = 2 * delta/mi^2 * (e - ee);
        e  = ee;

        % If Dc > 0, then the stress can be successfully decreased.
        success = (Dc >= 0);
        if (success)
          y       = yy;
          [ee,a]  = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio);
          g1      = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid);
          gnorm2  = g1(:)' * g1(:);
          lambda1 = 0;

          beta = max(0,(g1(:)'*(g1(:)-g0(:)))/mi);
          p   = -g1 + beta .* p;                 % P is a new conjugate direction

          if (g1(:)'*p(:) >= 0 | mod(it-1,n*m) == 0),
            p = -g1;                            % No much improvement, restart
          end

          if (Dc >= 0.75)
            lambda = 0.25 * lambda;             % Make the regularization smaller
          end

          it  = it + 1;
          err = [err; e];
          prwaitbar(opt.maxiter,it,[tt num2str(it) ', error: ' num2str(e)]);
          %fprintf (opt.st,'iteration: %4i   stress: %3.8d\n',it,e);

        else        % Dc < 0
          % Note that for Dc < 0, the iteration number IT is not increased,
          % so the Hessian is further regularized until SUCCESS equals 1.
          lambda1 = lambda;
        end

      % The approximation of the Hessian was poor or the stress was not
      % decreased (Dc < 0), hence the regularization lambda is enlarged.
      if (Dc < 0.25)
        lambda = lambda + delta * (1 - Dc)/pnorm2;
      end
    end
  end
  prwaitbar(0);

return

% **********************************************************************************

%MDS_STRESS - Calculate Sammon stress during optimization
%
%   E = MDS_SAMSTRESS(Q,Y,YREP,Ds,D,NANINDEX)
%
% INPUT
%   Q         Indicator of the Sammon stress; Q = -2,-1,0,1,2
%   Y         Current lower-dimensional configuration
%   YREP      Configuration of the representation objects; it should be
%             empty when no representation set is considered
%   Ds        Original distance matrix
%   D         Approximate distance matrix (optional; otherwise computed from Y)
%   NANINDEX  Index of the missing values; marked in Ds by NaN (optional; to
%             be found in Ds)
%
% OUTPUT
%   E         Sammon stress
%
% DESCRIPTION
% Computes the Sammon stress between the original distance matrix Ds and the
% approximated distance matrix D between the mapped configuration Y and the
% configuration of the representation set YREP, expressed as follows:
%
% E = 1/(sum_{i<j} Ds_{ij}^(q+2)) sum_{i<j} (Ds_{ij} - D_{ij})^2 * Ds_{ij}^q
%
% It is directly used in the MDS_SAMMAP routine.

%
% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [e,alpha] = mds_samstress (q,y,yrep,Ds,D,nanindex,isratio)

  % If D is not given, calculate it from Y, the current mapped points.

  if (nargin < 5) | (isempty(D))
    if (~isempty(yrep))
      D = sqrt(distm(y,yrep));
    else
      D = sqrt(distm(y));
    end
  end

  if (nargin < 6)
    nanindex = [];    % Not given, so calculate below.
  end

  if (nargin < 7)
    isratio = 0;      % Assume this is meant.
  end

  todefine = isempty(yrep);

  [m,n]  = size(y); [mm,k] = size(Ds);
  if (m ~= mm)
    error('The sizes of Y and Ds do not match.');
  end
  if (any(size(D) ~= size(Ds)))
    error ('The sizes of D and Ds do not match.');
  end
  m2 = m*k;

  % Convert to double.
  D  = +D; Ds = +Ds;

  % I is the index of non-NaN, non-zero (> eps) values to be included
  % for the computation of the stress.

  I = 1:m2;
  if (~isempty(nanindex))
    I(nanindex) = [];
  end

  O = [];
  if (todefine), O = 1:m+1:m2; end
  I = setdiff(I,O);

  % If OPTIONS.isratio is set, calculate optimal ALPHA to scale with.

  if (isratio)
    alpha = sum((Ds(I).^q).*D(I).^2)/sum((Ds(I).^(q+1)).*D(I));
    Ds    = alpha*Ds;
  else
    alpha = 1;
  end

  % C is the normalization factor.
  c = sum(Ds(I).^(q+2));

  % If Q == 0, prevent unnecessary calculation (X^0 == 1).
  if (q ~= 0)
    e = sum(Ds(I).^q .*((Ds(I)-D(I)).^2))/c;
  else
    e = sum(((Ds(I)-D(I)).^2))/c;
  end

return

% **********************************************************************************

%MDS_GRADIENTS - Gradients for variants of the Sammon stress
%
%   [G1,G2,CC] = MDS_GRADIENTS(Q,Y,YREP,Ds,D,NANINDEX)
%
% INPUT
%   Q          Indicator of the Sammon stress; Q = -2,-1,0,1,2
%   Y         Current lower-dimensional configuration
%   YREP      Configuration of the representation objects; it should be
%             empty when no representation set is considered
%   Ds        Original distance matrix
%   D         Approximate distance matrix (optional; otherwise computed from Y)
%   nanindex  Index of missing values; marked in Ds by NaN (optional;
%             otherwise found from Ds)
%
% OUTPUT
%   G1        Gradient direction
%   G2        Approximation of the Hessian by its diagonal
%
% DESCRIPTION
% This is a routine used directly in the MDS_SAMMAP routine.
%
% SEE ALSO (PRTools Guide)
% MDS, MDS_INIT, MDS_SAMMAP, MDS_SAMSTRESS

%
% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [g1,g2,c] = mds_gradients(q,y,yrep,Ds,D,nanindex)

  % If D is not given, calculate it from Y, the current mapped points.

	y = +y;
	yrep = +yrep;
  if (nargin < 5) | (isempty(D))
    if (~isempty(yrep))
      D = sqrt(distm(y,yrep));
    else
      D = sqrt(distm(y));
    end
  end

  % If NANINDEX is not given, find it from Ds.

  if (nargin < 6)
    nanindex = find(isnan(Ds(:))==1);
  end

  % YREP is empty if no representation set is defined yet.
  % This happens when the mapping should be defined.

  todefine = (isempty(yrep));
  if (todefine)
    yrep  = y;
  end

  [m,n]  = size(y); [mm,k] = size(Ds);
  if (m ~= mm)
    error('The sizes of Y and Ds do not match.');
  end
  if (any(size(D) ~= size(Ds)))
    error ('The sizes of D and Ds do not match.');
  end
  m2 = m*k;

  % Convert to doubles.

  D  = +D; Ds = +Ds;

  % I is the index of non-NaN, non-zero (> eps) values to be included
  % for the computation of the gradient and the Hessians diagonal.

  I = (1:m2)';
  if (~isempty(nanindex))
    I(nanindex) = [];
  end
  K    = find(Ds(I) <= eps | D(I) <= eps);
  I(K) = [];

  % C is a normalization factor.

  c = -2/sum(Ds(I).^(q+2));

  % Compute G1, the gradient.

  h1 = zeros(m2,1);
  if (q == 0)                    % Prevent unnecessary computation when Q = 0.
    h1(I) = (Ds(I)-D(I)) ./ D(I);
  else
    h1(I) = (Ds(I)-D(I)) ./ (D(I).*(Ds(I).^(-q)));
  end
  h1 = reshape (h1',m,k);
  g2 = h1 * ones(k,n);          % Here G2 is assigned only temporarily,
  g1 = c * (g2.*y - h1*yrep);    % for the computation of G1.

  % Compute G2, the diagonal of the Hessian, if requested.

  if (nargout > 1)
    h2 = zeros(m2,1);
    switch (q)
      case -2,
        h2(I) = -1./(Ds(I).*D(I).^3);
      case -1,
        h2(I) = -1./(D(I).^3);
      case 0,
        h2(I) = - Ds(I)./(D(I).^3);
      case 1,
        h2(I) = - Ds(I).^2./(D(I).^3);
      case 2,
        h2(I) = -(Ds(I)./D(I)).^3;
      end
    h2 = reshape (h2',m,k);
    g2 = c * (g2 + (h2*ones(k,n)).*y.^2 + h2*yrep.^2 - 2*(h2*yrep).*y);
  end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function w = flattopwin (L, sym)
% Author: Paul Kienzle <pkienzle@users.sf.net> (2004)
% This program is granted to the public domain.
%
% flattopwin(L, [periodic|symmetric])
%
% Return the window f(w):
%
%   f(w) = 1 - 1.93 cos(2 pi w) + 1.29 cos(4 pi w)
%            - 0.388 cos(6 pi w) + 0.0322cos(8 pi w)
%
% where w = i/(L-1) for i=0:L-1 for a symmetric window, or
% w = i/L for i=0:L-1 for a periodic window.  The default
% is symmetric.  The returned window is normalized to a peak
% of 1 at w = 0.5.
%
% This window has low pass-band ripple, but high bandwidth.
%
% According to [1]:
%
%    The main use for the Flat Top window is for calibration, due
%    to its negligible amplitude errors.
%
% [1] Gade, S; Herlufsen, H; (1987) 'Use of weighting functions in DFT/FFT
% analysis (Part I)', Bruel & Kjaer Technical Review No.3.

  if nargin == 0 || nargin > 2
    help(mfilename);
  end % if

  divisor = L-1;
  if nargin > 1
    if strcmp(sym, 'periodic')
      divisor = L;
    elseif strcmp(sym, 'symmetric')
      divisor = L-1;
    else
      error('second argument must be ''periodic'' or ''symmetric''');
    end
  end % if

  x = 2*pi*(0:(L-1))'/divisor;
  if L==1
    w = 1;
  else
    % these coefficients come from wikipedia, and match better with what
    % matlab outputs (JM)
    a0 = 0.21557895;
    a1 = 0.41663158;
    a2 = 0.277263158;
    a3 = 0.083578947;
    a4 = 0.006947368;
    w = a0 - a1*cos(x) + a2*cos(2*x) - a3*cos(3*x) + ...
            a4*cos(4*x);
    %w = (1-1.93*cos(x)+1.29*cos(2*x)-0.388*cos(3*x)+0.0322*cos(4*x))/4.6402;
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x=choosenk(n,k)
%CHOOSENK All choices of K elements taken from 1:N [X]=(N,K)
% The output X is a matrix of size (N!/(K!*(N-K)!),K) where each row
% contains a choice of K elements taken from 1:N without duplications.
% The rows of X are in lexically sorted order.
%
% To choose from the elements of an arbitrary vector V use
% V(CHOOSENK(LENGTH(V),K)).

% CHOOSENK(N,K) is the same as the MATLAB5 function NCHOOSEK(1:N,K) but is
% much faster for large N and most values of K.

%   Copyright (c) 1998 Mike Brookes,  mike.brookes@ic.ac.uk
%      Version: $Id: choosenk.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

kk=min(k,n-k);
if kk<2
   if kk<1
      if k==n
         x=1:n;
      else
         x=[];
      end
   else
      if k==1
         x=(1:n)';
      else
         x=1:n;
         x=reshape(x(ones(n-1,1),:),n,n-1);
      end
   end
else
   n1=n+1;
   m=prod(n1-kk:n)/prod(1:kk);
   x=zeros(m,k);
   f=n1-k;
   x(1:f,k)=(k:n)';
   for a=k-1:-1:1
      d=f;
      h=f;
      x(1:f,a)=a;
      for b=a+1:a+n-k
         d=d*(n1+a-b-k)/(n1-b);
         e=f+1;
         f=e+d-1;
         x(e:f,a)=b;
         x(e:f,a+1:k)=x(h-d+1:h,a+1:k);
      end
   end
end

author: decouples