MATLAB: переписать код в алгебраических терминах

stream_24

Подскажите, п-та, можно ли переписать следующий участок кода без циклов
 

N = length(x); M = length(y);
up = zeros(N,M); um = up;
for n=1:N
for m=1:M
up(n,m) = x(n) + y(m);
um(n,m) = x(n) - y(m);
end
end

Спасибо!

Vlad128

тут поможет функция meshgrid

stm7518204

N = length(x); M = length(y);
up = zeros(N,M); um = up;
    for m=1:M
     up(:,m) = x(:) + y(m);
     um(:,m) = x(:) - y(m);
    end
Так работает?

stream_24

а от второго цикла как избавиться?;-)

stm7518204

а от второго цикла как избавиться?;-)
просто через оператор : у меня никогда не удавалось избавиться от второго цикла... Однако 5 лет назад был случай, когда удалось одному парню сделать это... но больше у него этот результат не повторялся..

Vlad128

В качестве шутки:
в случае с суммой можно взять логарифмы, перемножить матрицы и взять (поэлементную) экспоненту :grin:
На разность тоже адаптируется :D

Vlad128

Типа up = log( exp( x ) * exp( y' ) );
Типа um = log( exp( x ) * exp( -y' ) );
это если size( x ) = [ N 1 ], size( y ) = [ M 1 ]
Возможно, у log основание другое, надо проверить.

stm7518204

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

Vlad128

Перепутал, экспоненту и логарифм местами надо поменять. В посте исправлено.

narkom

что-то типа должно сработать:
X = repmat(x(1.M);
Y = transpose(repmat(y(1,N;
up = X+Y;
um = X-Y;
transpose если не комплексные можно заменить на '

stm7518204

Типа up = log( exp( x ) * exp( y' ) );
Типа um = log( exp( x ) * exp( -y' ) );
это если size( x ) = [ N 1 ], size( y ) = [ M 1 ]
Возможно, у log основание другое, надо проверить.
работает!

narkom

ох не устойчивые у тебя шутки :)

Vlad128

никто не спорит :grin:

Vlad128

Кстати, нет матлаба под рукой, но что-то мне кажется, что тут кошерно заюзать bsxfun, не помню, правда, развернет ли он так, как надо, но практически уверен.

stm7518204

что-то типа должно сработать:
X = repmat(x(1.M);
Y = transpose(repmat(y(1,N;
up = X+Y;
um = X-Y;
transpose если не комплексные можно заменить на '
Тоже работатет

Vlad128

А, я понял, это намек. Попробуй, пожалуйста
up = bsxfun( @plus, x, y' );

stm7518204

Вопрос отцам матлаба...
Есть функция find которой можно искать число (например 1 )в массиве. типа [u,v] = find(A==1);
А как устроить тоже самое, но просматривать надо значение филдов у structure.
типа
s(1).a =1;
s(2).a=2;
....
Как без цикла записать в стиле find(s.a ==1)...... а то такая команда не работатет

lenmas

Может быть, что-то типа reshape(x,length(x1)*ones(1,length(y)+ones(length(x1)*reshape(y,1,length(y?

Vlad128

ты думаешь, умножение на ones будет лучше, чем repmat?
Я, кстати, не удивлюсь, если в этом гавене через repmat быстрее работает, чем через bsxfun.

stm7518204

А, я понял, это намек. Попробуй, пожалуйста
up = bsxfun( @plus, x, y' );
up = bsxfun( @plus, x, y' ); - работат!

um = bsxfun( @minus, x', y ); работает только вот так вот надо видоизменить..

narkom

reshape(x,length(x1)*ones(1,length(y)+ones(length(x1)*reshape(y,1,length(y
чем это репмата отличается? ну кроме более высокой цены? :)
ПС reshape(x,length(x1) == x(:)

Vlad128

В общем, bsxfun — это правильное решение.

narkom

Я, кстати, не удивлюсь, если в этом гавене через repmat быстрее работает, чем через bsxfun.
вполне возможно. Зато читается легче. У меня голова кругом идет после решейпов и репматов одних за другим.
Кстати, долго ломал голову как сделать без циклов. Не допер:
как преобразовать [X_1 X_2 .. X_n] в [X_1; X_2; .. X_n]
X_i квадратные блоки заданной величины.

stm7518204

Это круто!
Но кто подскажет про find в филдах ?

Vlad128

Почитай где-нибудь здесь :grin:
http://www.mathworks.com/access/helpdesk/help/techdoc/matlab...

stm7518204


Спасибо зачитался

Vlad128

Честно говоря тоже сразу не приходит в голову, но здесь работает отмазка: число итераций цикла значительно меньше числа действий для выполнения этого преобразования, поэтому можно циклом :grin:
Да даже скорее так и надо, если нету стандартной функции для этого дела.

lenmas

ПС reshape(x,length(x1) == x(:)
Мне нужен именно столбец был. А так ничем, задача-то ставилась без циклов написать. Или тут быстродействие было прописано в условиях задачи? :grin:

Vlad128

transpose если не комплексные можно заменить на '
есть ".'"

lenmas

как преобразовать [X_1 X_2 .. X_n] в [X_1; X_2; .. X_n]
X_i квадратные блоки заданной величины.
Такое ощущение, что здесь собрание неудачников, которые книжку не смогли прочитать по матлабу :grin:
Если что, то я про себя тоже :)

narkom

вообще логично :) итераций мало я подумал слекга и забил. Все равно не так много времени берет. Учитывая, что потом квадратичная программа вычисляется с этими матрицами то не особо критично :)

narkom

Мне нужен именно столбец был.
ну дык x(:) и дает столбец! .... или строку... :grin:

Vlad128

А так ничем, задача-то ставилась без циклов написать. Или тут быстродействие было прописано в условиях задачи? :grin:
математик :smirk:

Vlad128

x(:) — столбец.
А что, x(:).' не работает?

narkom

А что, x(:).' не работает?
да черт знает, надо шуми спросить :) не пробовал

Boris

Как без цикла записать в стиле find(s.a ==1)
 find([s.a]==1) 

нет?

Vlad128

да, кстати, comma-separated list можно загнать в array, ты должно быть прав.

stream_24

Спасибо!
Решение с meshgrid наиболее универсальное, только плохо работает с M=1.

Vlad128

не-не, юзай bsxfun, оно наиболее универсальное, что с ним не так?

lenmas

математик :smirk:
Служу Советскому Союзу!
А так предпочитаю прозрачность и краткость кода быстродействию. С современными компами можно особо не париться над производительностью :)

narkom

А так предпочитаю прозрачность и краткость кода быстродействию. С современными компами можно особо не париться над производительностью
брехня. Ты б ещё суперкомьютеры как довод приводил. Чипы в самолетах или машинах в бортовых компьютерах слабее чем айпод. У меня программы неделями гоняются, и при этом они уже оптимизированы. Тут конечно ещё матлаб играет роль, но оптимизировать код очень полезно.

lenmas

брехня. Ты б ещё суперкомьютеры как довод приводил. Чипы в самолетах или машинах в бортовых компьютерах слабее чем айпод. У меня программы неделями гоняются, и при этом они уже оптимизированы. Тут конечно ещё матлаб играет роль, но оптимизировать код очень полезно.
Ну, ты лучше знаешь, что тебе нужно :)
У меня скрипты на матлабе максимум полчаса выполнялись.
А что, в самолетах или в бортовых компах работают программы на матлабе? :shocked:

narkom

нет, ты прав :) Но это просто пример к оптимизации.
К тому же сейчас тренд идет на увеличение количества процессоров в пк и тут векторная арифметика очень полезна. Циклы могут заметно замедлить программу. Вообще от применения зависит, конечно же. Тебе может и нужно, другим очень нужно. Так что спорить не вижу смысла. Вполне допускаю, что многим быстродействие не особо интересно. Мне в частности нужно, у меня размерности матриц на близки к десяткам тысяч.

MammonoK

пиши на фортране если тебе критично быстродействие :)

Vlad128

repmat более читаем, чем умножение на ones, bsxfun более читаем, чем repmat.

narkom

к матлабу очень много всяких тулбоксов и полностью переписывать тулбоксы нет ни малейшего желания :)

stream_24

для FORTRANа "тулбоксов" как минимум не меньше:)

stream_24

Это для данного примера более читаемо bsxfun, а если надо не прибавить/вычесть, а записать формулу посложнее (например, плотность двумерного нормального распределния)?
ПО тестам производительности и наглядности кода пока что выигрывает meshgrid. Ещё раз спасибо!

Vlad128

Это для данного примера более читаемо bsxfun, а если надо не прибавить/вычесть, а записать формулу посложнее (например, плотность двумерного нормального распределния)?
up = bsxfun( @(x,y) <тут любая формула>, x, y' ) ?

narkom

но ведь не в количестве дело :) нее, согласен, что если уж совсем туго, то надо переписывать. Хотя бы тяжелые части. Ну хотя бы в си++ и включать в матлаб. Хотя порою и это не помогает.

Vlad128

не, матлаб для разработок промышленного кода вообще не подходит, имхо. Он помогает прикинуть алгоритм и т.д.
По крайней мере сейчас он слишком сырой, оптимизирует код ужасно, с такой-то информацией.

lenmas

для FORTRANа "тулбоксов" как минимум не меньше:)
На С больше, поэтому часто приходится самому лабать.

lenmas

В матлабе удобно мне кажется решать только диффуры, то-есть интегрировать разные системы, моделировать системы управления, в механике той же очень удобно задавать системы уравнений Эйлера сразу в векторном виде и так далее. Для перемножения матриц как-то нерентабельно в системах реального времени его гонять.

Vlad128

fortran и C хорошо совмещаются

lenmas

fortran и C хорошо совмещаются
Ну вот люди, которые занимаются реальными квантомеханическими расчетами, говорят, что библиотек для фортрана меньше, дело касается распараллеливания.

stm7518204

Как без цикла записать в стиле find(s.a ==1)
code: find([s.a]==1)
нет?
Работает! Круто! Спасибо!
to all.
Не понял про x(:).' - это про что?

Vlad128

Не понял про x(:).' - это про что?
Это более короткий аналог reshape(x,1,length(x

stm7518204

Это более короткий аналог reshape(x,1,length(x
если x-столбец?

Vlad128

если х — что угодно, не обязательно столбец.
Если известно, что х — столбец, то можно просто х.'
А если известно, что строка, то можно просто x

lenmas

если х — что угодно, не обязательно столбец.
Если известно, что х — столбец, то можно просто х.'
А если известно, что строка, то можно просто x
Забавно, не знал, что двоеточие всегда обрабатывается как строка :)
Просто часто в матрицах использовал диапазоны, а там где строка, а где и столбец.

Vlad128

Забавно, не знал, что двоеточие всегда обрабатывается как строка
как столбец.
Вообще, если индекс один, то он по внутреннему представлению: там все матрицы вытянуты в порядке размерностей. Для двумерной матрицы это все столбцы по порядку, для трехмерной — все слои, каждый слой — как двумерная матрица и т.д.

lenmas

ля двумерной матрицы это все столбцы по порядку
Это-то я хорошо усвоил, что в памяти массивы упорядочены по столбцам --- сначала первый столбец (когда первый индекс меняется) записывается, потом второй после него (когда изменился на единицу второй индекс, и опять пошел первый индекс меняться) и т/д. Про одномерный массив вроде в этом свете стало понятно.

stream_24

Для быстрого получения (в том числе и не программистами) результатов MATLAB - то, что нужно. Кроме того, MATLAB можно успешно использовать только для графичечкого вывода.

Vlad128

MATLAB можно успешно использовать только для графичечкого вывода.
тут не понял немного. Почему только для графического?

stream_24

Все ресурсоёмкие вычисления производить на C, а для графики/видео использовать MATLAB.

Vlad128

ну тогда можно и гнуплот или еще чего использовать.
Я матлаб рассматриваю исключительно как коллекцию реализованных алгоритмов, где можно затестить применение этих самых алгоритмов. А потом уже писать на C/C++, да.

lenmas

А потом уже писать на C/C++, да.
Что, и Рунге-Кутта будешь переписывать на С/С++/Fortran? :grin:

Vlad128

А он еще не написан?
Тем более нашел алгоритм, он же отстойный.

lenmas

А он еще не написан?
Ну, ты ж пишешь, что сам будешь его переписывать :grin:
В матлабе откуда ты знаешь, как реализован алгоритм, и совпадает ли он на самом деле с реализованным в библиотеке.

lenmas

Тем более нашел алгоритм, он же отстойный.
А какой ты считаешь не отстойным?

Boris

А чего, большая проблема Рунге-Кутта написать? :)

Vlad128

А какой ты считаешь не отстойным?
ладно, я имел в виду стандартный рунге-кутта второго порядка, он простой и беспонтовый.

Vlad128

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

lenmas

А чего, большая проблема Рунге-Кутта написать? :)
Проблема не в том, чтобы написать, а в том, чтобы не заниматься лишней деятельностью :)
Если на такие мелкие подзадачки отвлекаться постоянно, то на основную не останется сил.
Твой вопрос хорош, если Рунге-Кутта изучать как самоцель, в книжках копаться, разбираться, вспоминать.
Я предпочитаю воспользоваться готовым матлабовским. :)

lenmas

ладно, я имел в виду стандартный рунге-кутта второго порядка, он простой и беспонтовый.
Мне в матлабе хватает даже 4-5-ого порядка (ode45 что ли называется).

lenmas

не писал я такого.
"А потом уже писать на C/C++, да." --- Я это понял так. :)

Vlad128

Что значит "даже"? :grin:

Vlad128

Неправильно понял. Повторное использование кода никто не отменял.

lenmas

Что значит "даже"? :grin:
Там вроде еще есть 6-7-ого порядка. Еще какие-то специальные для жестких систем, но это уже для специальных задач.

Vlad128

А, понял.
Ну так и как ты предлагаешь писать код, который должен работать в реальном времени, я так и не понял.

lenmas

Ну так и как ты предлагаешь писать код, который должен работать в реальном времени, я так и не понял.
Мы уже от реального времени отошли вроде в обсуждении.
А в реальном времени используют максимум ломаные Эйлера, какие тебе там высокие порядки! :grin: Там ошибки данных такие, что гладкие методы просто убьются насмерть :)

Vlad128

Реальное время реальному времени рознь. И не только ведь в диффурах дело :grin:
Не, ну т.е. ты предлагаешь на матлабе писать, я правильно понял?
Ну там, обработку изображений, распознавание.

lenmas

Не, ну т.е. ты предлагаешь на матлабе писать, я правильно понял?
Ну там, обработку изображений, распознавание.
Да, спокойно использовать в моделировании. Как показывает опыт (да, признаюсь, мой опыт :) ошибки гуляют в мозгах авторов моделей, но никак не в алгоритмах. Если хочется проверить правильность вычислений при отлаженном алгоритме, то попробовать другим методом, с большими точностями.

narkom

можно так. Приблизить дифур (или интегро-дифференциальное ур-ие) системой линейных уравнений, а она уже реализуется на уровне железа. Даже считать ничего не надо :o

Vlad128

да не в правильности дело, а в тормознутости и тяжелости.

Boris

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

Vlad128

ode45 не всегда приятен в использовании, например если неизвестная функция мыслится не как вектор, а как матрица.
там оберточка на это дело пишется в 6 строк. Могу подкинуть :grin:

Boris

Подкинь :o

Vlad128

function Q = matrixOde45(A, tspan, X0, options)
    N = size( X0, 1 );
    [t,Q] = ode45( @(t,y) matrixPull( A, t, y, N tspan, X0( options );
    Q = reshape( Q, length( tspan N, N );
end

function X = matrixPull(A, t, y, N)
    X = A( t, reshape( y, N, N ) );
    X = X(:);
end

A — матричная функция от t.чушь какую-то написал. A — функция от t и матрицы.
Например, A = @(t,X) [ t 0; 0 t ] * X.
reshape в конце сделан, чтобы была группировка по времени (элементы с близким временем расположены близко в памяти что логично и в той задаче, для которой была написана эта функция, было удобно.
Фишка в том, что reshape при вычислении matrixPull вообще должен оптимизироваться в ничто. Но это же, блин, матлааааб.

Boris

Спасибо :)
* Задумался, откуда у меня взялась аллергия на такие обертки.

Vlad128

* Задумался, откуда у меня взялась аллергия на такие обертки.
доложи о результатах, интересно.
Оставить комментарий
Имя или ник:
Комментарий: