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

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

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

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

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

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

stm7518204

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

Vlad128

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

lenmas

ПС reshape(x,length(x1) == x(:)

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

Vlad128

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

есть ".'"

lenmas

как преобразовать [X_1 X_2 .. X_n] в [X_1; X_2; .. X_n]
X_i квадратные блоки заданной величины.

Такое ощущение, что здесь собрание неудачников, которые книжку не смогли прочитать по матлабу
Если что, то я про себя тоже

narkom

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

narkom

Мне нужен именно столбец был.

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

Vlad128

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

математик

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

математик

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

narkom

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

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

lenmas

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

Ну, ты лучше знаешь, что тебе нужно
У меня скрипты на матлабе максимум полчаса выполнялись.
А что, в самолетах или в бортовых компах работают программы на матлабе?

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?

Vlad128

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

lenmas

А он еще не написан?

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

lenmas

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

А какой ты считаешь не отстойным?

Boris

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

Vlad128

А какой ты считаешь не отстойным?

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

Vlad128

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

не писал я такого.

lenmas

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

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

lenmas

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

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

lenmas

не писал я такого.

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

Vlad128

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

Vlad128

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

lenmas

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

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

Vlad128

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

lenmas

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

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

Vlad128

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

lenmas

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

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

narkom

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

Vlad128

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

Boris

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

ode45 не всегда приятен в использовании, например если неизвестная функция мыслится не как вектор, а как матрица.

Vlad128

ode45 не всегда приятен в использовании, например если неизвестная функция мыслится не как вектор, а как матрица.

там оберточка на это дело пишется в 6 строк. Могу подкинуть

Boris

Подкинь

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

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

доложи о результатах, интересно.