Векторизация - функция Сумма и Бесселя

Кто-нибудь может помочь векторизовать этот код Matlab? Конкретной проблемой является функция суммы и бессела с векторными входами. Спасибо!

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];
for ii = 1:length(rho_g)
    for jj = 1:length(phi_g)
        % Coordinates
        rho_o = rho_g(ii);
        phi_o = phi_g(jj);
        % factors
        fc = cos(n.*(phi_o-phi_s));
        fs = sin(n.*(phi_o-phi_s));

        Ez_t(ii,jj) = sum(tau.*besselj(n,k(3)*rho_s).*besselh(n,2,k(3)*rho_o).*fc);
    end
end
Источник

3 ответа

Решение

Инициализация -

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

Форма вложенных циклов (скопируйте из своего кода и показана здесь только для сравнения) -

for ii = 1:length(rho_g)
    for jj = 1:length(phi_g)
        % Coordinates
        rho_o = rho_g(ii);
        phi_o = phi_g(jj);
        % factors
        fc = cos(n.*(phi_o-phi_s));
        fs = sin(n.*(phi_o-phi_s));

        Ez_t(ii,jj) = sum(tau.*besselj(n,k(3)*rho_s).*besselh(n,2,k(3)*rho_o).*fc);
    end
end

Векторизованное решение -

%%// Term - 1
term1 = repmat(tau.*besselj(n,k(3)*rho_s),[N*N 1]);

%%// Term - 2
[n1,rho_g1] = meshgrid(n,rho_g);
term2_intm = besselh(n1,2,k(3)*rho_g1);
term2 = transpose(reshape(repmat(transpose(term2_intm),[N 1]),N,N*N));

%%// Term -3
angle1 = repmat(bsxfun(@times,bsxfun(@minus,phi_g,phi_s')',n),[N 1]);
fc = cos(angle1);

%%// Output
Ez_t = sum(term1.*term2.*fc,2);
Ez_t = transpose(reshape(Ez_t,N,N));

Следует отметить, что это векторизация или упрощение кода.

  1. 'fs' не изменяет вывод скрипта, Ez_t, поэтому он может быть удален на данный момент.
  2. Выходными данными, по-видимому, являются 'Ez_t', для которых в коде требуются три основных термина: - tau.* Besselj (n, k (3) * rho_s), besselh (n, 2, k (3) * rho_o) и fc. Они рассчитываются отдельно для векторизации как термины 1,2 и 3 соответственно.
  3. Все эти три термина имеют размеры 1xN. Таким образом, наша цель состоит в том, чтобы вычислить эти три члена без петель. Теперь два цикла выполняются по N раз каждый, что дает нам общее количество циклов NxN. Таким образом, мы должны иметь NxN раз данных в каждом таком термине по сравнению с тем, когда эти термины были внутри вложенных циклов.
  4. Это в основном суть векторизации, выполненной здесь, так как три термина представлены самими "term1", "term2" и "fc".
Источник

Вы можете попытаться векторизовать этот код, что может быть возможно с некоторыми bsxfun или так, но было бы трудно понять код, и вопрос в том, будет ли он работать быстрее, поскольку ваш код уже использует векторную математику во внутреннем цикле (даже если ваши векторы имеют длину только 3). Полученный код станет очень трудным для чтения, поэтому вы или ваш коллега не будете иметь представления о том, что он делает, когда вы посмотрите на него через 2 года.

Прежде чем тратить время на векторизацию, гораздо важнее узнать о движении инвариантного кода цикла, которое легко применить к вашему коду. Некоторые наблюдения:

  • ты не используешь fsтак что убери это.
  • семестр tau.*besselj(n,k(3)*rho_s) не зависит ни от одной из ваших переменных цикла ii а также jj, так что это постоянно. Рассчитайте это один раз перед вашей петлей.
  • Вы, вероятно, должны предварительно выделить матрицу Ez_t,
  • единственные термины, которые изменяются во время цикла fc, который зависит от jj, а также besselh(n,2,k(3)*rho_o), который зависит от ii, Я думаю, что последний стоит гораздо больше времени для расчета, поэтому лучше не рассчитывать это N*N раз во внутреннем цикле, но только N раз во внешнем цикле. Если расчет основан на jj займет больше времени, вы можете поменять циклы ii а также jj, но это, кажется, не имеет место здесь.

Код результата будет выглядеть примерно так (не проверено):

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

% constant part, does not depend on ii and jj, so calculate only once!
temp1 = tau.*besselj(n,k(3)*rho_s); 

Ez_t = nan(length(rho_g), length(phi_g)); % preallocate space
for ii = 1:length(rho_g)
    % calculate stuff that depends on ii only
    rho_o = rho_g(ii);
    temp2 = besselh(n,2,k(3)*rho_o);

    for jj = 1:length(phi_g)
        phi_o = phi_g(jj);
        fc = cos(n.*(phi_o-phi_s));
        Ez_t(ii,jj) = sum(temp1.*temp2.*fc);
    end
end
Источник

Чтобы дать самостоятельный ответ, я скопирую исходную инициализацию

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

и сгенерировать некоторые отсутствующие данные (k(3) и rho_s и phi_s в измерении n)

rho_s = rand(size(n));
phi_s = rand(size(n));
k(3) = rand(1);

тогда вы можете вычислить тот же Ez_t с многомерными массивами:

[RHO_G, PHI_G, N] = meshgrid(rho_g, phi_g, n);
[~, ~, TAU] = meshgrid(rho_g, phi_g, tau);
[~, ~, RHO_S] = meshgrid(rho_g, phi_g, rho_s);
[~, ~, PHI_S] = meshgrid(rho_g, phi_g, phi_s);

FC = cos(N.*(PHI_G - PHI_S));
FS = sin(N.*(PHI_G - PHI_S)); % not used

EZ_T = sum(TAU.*besselj(N, k(3)*RHO_S).*besselh(N, 2, k(3)*RHO_G).*FC, 3).';

После этого вы можете проверить, что обе матрицы одинаковы

norm(Ez_t - EZ_T)
Источник
Другие вопросы по тегам