Mathematica: Матрица в цикле

[Sergey K]

amax := 25
z = Table [FindE[mm, a, 0], {mm, 2, amax}, {a, 1, amax}]

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

pr = z[[1, 25]]
procedure[i_, j_] := 
 Module[{ii = i, jj = j}, 
  While[jj < amax, {If[z[[ii, jj]] == pr, z[[ii, jj]] = 0], jj++}]]
i := 1
While [i < amax, {procedure[i, 1], i++} ];
z

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

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

pr = z[[1, 25]]
i = 1
While [i < 25, {If[z[[i, 25]] == pr, z[[i, 25]] = 0], i++} ]
z

Это просто для фильтрации последнего столбца на битые значения. Однако ошибка Set::noval: Symbol z in part assignment does not have an immediate value и не измененная матрица сопровождаем меня и тут. По справке причину возникновения ошибки не понял. В чем дело?

[Vince Diesel]

Получить нужный результат можно просто подстановкой:

z/. pr -> 0

если pr - максимальное значение.

[Aritaborian]

Sergey K, в Mathematica процедурный подход используют, как правило, от безысходности ;-) Почти всегда можно применить функциональный.

[Sergey K]

Я к сожалению очень долго пользовался процедурным, а функциями - только учусь. подскажите, как начать? вот как это написать функциями?

Цитата:

Получить нужный результат можно просто подстановкой:

большое спасибо, но это часть проблемы. некоторые битые значения надо заменить не на 0, а на соседнюю ячейку. Процедурно представляю, функционально - нет.

[Vince Diesel]

Фигурные скобки в математике служат для образования массивов, а не для выделения блоков.
А без циклов можно сделать так:

f[x_,{i_,j_}]:=If[x<pr,x,<здесь команда, высчисляющая значение z[[i,j]], если x=pr>]
z=MapIndexed[f, z, {2}]

В <> можно использовать индексы {i,j}. Например,

f[x_, {i_, j_}] := If[x < pr, x, If[j > 1, z[[i, j - 1]], 0]];

[Sergey K]

z=MapIndexed[f, z, {2}]

не понятно что делает этот код. точнее именно {2}

[arseniiv]

Нажмите F1, когда курсор стоит на MapIndexed. (Может быть, придётся раскрыть раздел «Details and Options». Ещё есть копия справки онлайн, для этой функции вот) В данном случае это спецификация уровня, говорящая «только на втором сверху уровне» (без скобок было бы «с 1 по 2»).

Уровень показывает, как глубоко внутри выражения находится что-нибудь. Например, в {1, 2, 3} на первом уровне находятся элементы: 1, 2, 3, а на нулевом — List. Других уровней в этом выражении нет. В более сложных может быть: в {1, {2, 3}} на втором уровне 2 и 3.

[Sergey K]

Честно говоря, так и не понял вашей реализации - наверно и потому, что вы не так поняли что мне надо)

amax := 50
i := 1
FunctionFind[mm_, a_] := 
  If[a < mm - 1 , FindE[mm, a, 0], 
   If[a == mm - 1, FindE[mm, a + 1, 0], 0]];
z = ParallelTable [FunctionFind[mm, a], {mm, 2, amax}, {a, 1, amax}]
While[i < amax, z[[i, i + 1]] = z[[i, i]]; i++]
MatrixForm[z]

Вот так я сделал. Считает долго, да. Думал, потому что это из-за неэффективность FindE - она работает перебором значений на диапазоне от А до Б (расчет лишних значений я отсек в условиях функции FunctionFind). Поскольку значения величин в матрице с лева на право возрастают, улучшил и ввел поиск не от А, а от предыдущего значения Б. dE здесь точность, шаг перебора.

amax := 50
i := 1
FunctionFind[mm_, a_] := 
  If[a < mm - 1 , FindE[mm, a, If[a < 1, z[[mm, a - 1]], Ed[a] - dE]],
    If[a == mm - 1, 
    FindE[mm, a + 1, If[a < 1, z[[mm, a - 1]], Ed[a] - dE]], 0]];
z = ParallelTable [FunctionFind[mm, a], {mm, 2, amax}, {a, 1, amax}]
While[i < amax, z[[i, i + 1]] = z[[i, i]]; i++]
MatrixForm[z]

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




Наверное, он не обращается в ходе вычисления z[mm,a] к уже посчитанному z[mm, a-1], а выстраивает всю цепочку заново, весь z ещё не заполнен до конца. Правильно ли я понимаю? Как можно исправить? Наверное Table реализован внутри как обычный цикл, поэтому и медленно.

[Leierkastenmann]

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

А так вот вам пока отвлеченный пример, как по какому-то условию заменять в матрице элемент на соседний без явных циклов. Тут просто в матрице случайных чисел элемент заменяется на соседний справа, если он более 0.5. Соседство в данном случае циклично, для последнего элемента строки соседом является первый. Если все элементы более 0.5, то строка остается неизменной.

n = 1000;
k = 2000;
rt = Table[Random[], {i, n}, {j, k}];
rtnew = FixedPoint[
     MapThread[If[#1 <= 0.5, #1, #2] &, {#, RotateLeft[#]}] &, #, 
     k] & /@ rt;

[Sergey K]

Цитата:

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

я это уже учел вот здесь, проблем ни каких:

FunctionFind[mm_, a_] := If[a < mm - 1 , FindE[mm, a, If[a < 1, z[[mm, a - 1]], Ed[a] - dE]], If[a == mm - 1, FindE[mm, a + 1, If[a < 1, z[[mm, a - 1]], Ed[a] - dE]], 0]];

проблема в том, нельзя ли сделать ещё быстрее:)

z = ParallelTable [FunctionFind[mm, a], {mm, 2, amax}, {a, 1, amax}]

по сути это цикл, как я понимаю. может есть альтернативные способы задания?

[Leierkastenmann]

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

[Sergey K]

Есть определенное выражение (уже ясно что физический неверное, но интересно узнать как пользоваться математикой правильно) которое дает энергии при котором из кластера в mm атомов, вылетает а одноатомных осколков с учетом того, что они часть энергии уносят. Поскольку я не смог решить его аналитический, приходится решать его перебором значений на диапазоне. Это долго, и чтобы не терять время я решил загнать все значения в матрицу z. Поскольку из mm-атомного кластера не может вылететь осколков больше чем mm, то все значения над главное диагональю (такие как z[[mm,a]], где a>m) должны быть заполнены 0 (чтобы не ждать, пока функция поиска весь диапазон проверяет и ничего там не находит). Далее, распад 3 атомного кластера на 2 и на 3 атома по сути одно и тоже - поэтому такие значения должны быть равны.

Все это худо бедно реализует следующий код:

(*Исходные данные*)

Ed[1] := 0
Ed[2] := 4
Ed[3] := 3.7
Ed[4] := 4 (*Эксперимент!*)
Ed[5] := 5.61 (*Эксперимент!*)
Ed[6] := 5.68 (*Эксперимент!*)
Ed[7] := 6.53 (*Эксперимент!*)
Ed[8] := 5.93 (*Эксперимент!*)
Ed[9] := 5.83 
For[n = 10, n < 100, n++, Ed[n] := 5.78]

(*Нахождение перебором*)

Eost[E_, m_] := E - (Ed[m] + (E - Ed[m])/m)

Fun2[En_, m0_] := 
  Module[{m = m0, EE = En}, 
   While[m > 0, 
    If[Eost[EE, m] < 0, Break[]]; {EE = Eost[EE, m]; m--}]; m0 - m];

dE := 0.1
FindE[m0_, n_, En_] := 
  Module[{m = m0, nn = n, EE = En}, 
   While[EE < 5000, If[nn == Fun2[EE, m], Break[]]; {EE = EE + dE}]; 
   EE];

(*Функция ищет энергию, при которой кластер дает n одноатомных осколков*)

(*генерирую матрицу из значений энергии распада для кластеров [mm,a] - масса, осколки*)

amax := 10 (*максимальный размер кластера*)
i := 1

FunctionFind[mm_, a_] := 
  If[a < mm - 1 , FindE[mm, a, If[a < 1, z[[mm, a - 1]], Ed[a] - dE]],
    If[a == mm - 1, 
    FindE[mm, a + 1, If[a < 1, z[[mm, a - 1]], Ed[a] - dE]], 0]];


z = ParallelTable [FunctionFind[mm, a], {mm, 2, amax}, {a, 1, amax}]

While[i < amax, z[[i, i + 1]] = z[[i, i]]; i++]
z
MatrixForm[z]

Собственно, код сделан был для проверки одной идеи в 0 приближении. Она не оправдалась, но интересно а как правильно нужно было записать? императивные программы идут на ура, а вот функционально - мозг пока не работает так.

[Leierkastenmann]

Ну тут самое медленное место - это FindE, и пытаться оптимизировать заполнение таблицы особого смысла не имеет. Один вызов FindE в случае, когда внутренний ее цикл дойдет до 5000, будет стоить гораздо больше, чем выполнение последующих операций с матрицей независимо от того выполняем мы эти операции функционально или процедурно. Так что для начала надо как-то оптимизировать FindE, если конечно постановка задачи это позволит сделать. А функцию FunctionFind в таком виде можно было и не писать, там много лишнего кода, который не выполняется, тот же результат будет в этом случае:

z = Table[
  If[a < mm - 1,
   FindE[mm, a, Ed[a] - dE],
   If[a == mm - 1,
    FindE[mm, a + 1, Ed[a] - dE],
    0]
   ],
  {mm, 2, amax}, {a, 1, amax}]

[Sergey K]

Цитата:

Ну тут самое медленное место - это FindE, и пытаться оптимизировать заполнение таблицы особого смысла не имеет.

к сожалению, даже когда я укоротил диапазон поиска этой функции в 5-10 процентах случаев, введя FindE[mm, a, If[a < 1, z[[mm, a - 1]], Ed[a] - dE]] это улучшило скорость выполнения на несколько секунд! это меня по подтолкнуло к мысли, что тормозит не только она.

но все равно спасибо)