При программировании сколько-то интересных алгоритмов, то и дело приходится встречаться с выкладками, которые в тензорном виде либо формулируются проще, чем в матричном, либо которые в матричном виде сформулировать вообще нельзя без лишних циклов.
Каждый раз, когда я с болью выходил из ситуации, пытаясь перевернуть формулу, чтобы её можно было представить в матричном/векторном виде, я задавался вопросом, а что если вообще отказаться от этой идеи и сразу писать код в тензорном виде? Посидев вечерок, кое-что для себя выработал, хочу обсудить. Основной интерес для меня представляет синтаксис использования, чтобы можно было как можно легче проводить паралли между тензорными формулами и кодом.
Рассмотрим пару примеров того, что получилось, на примере формул из этого пейпера: http://matthias-mueller-fischer.ch/publications/strainBasedDynamics.pdf
Начнём с простого, в формуле (5) вводится матрица из трёх трёхкомпонентных векторов-столбцов.
Tensor< Index2<3, 3>, float> positions; //Index2 — второго ранга, 3, 3 — границы индексов
for (auto index : positions.Index( )) //index — комбинированный индекс, состоящий из index.i(), index.j() { positions.Get( index) = particleSystem->GetParticle( triangles[index.j( )])->pos[index.i( )]; }
Tensor< Index2<3, 2>, float> positions; //итерирование будет автоматически осуществляться по индексам [0..2], [0..1]
обычные матричные операции реализуются несколько длиннее, чем в матричном виде, но, что очень важно, все транспонирования задаются в явном виде и ошибки вроде умножения на матрицу слева или справа допустить гораздо труднее, потому что в тензорном виде не играет роли порядок умножения, важен лишь порядок индексов. К примеру, матричная формула вида
\(S=F^TF\), записывается в тензорном виде так:
\(S_{ij}=F_{ki}F_{kj}\) и прямым образом переводится в код, в явном виде указывая суммирование по повторяющемуся индексу:
Tensor< Index2<2, 2>, float> strain; for (auto index : strain.Index( )) { strain.Get( index) = 0; for ( auto sum : Index1<3>( )) { strain.Get( index) += deformations.Get( sum.i( ), index.i( )) * deformations.Get( sum.i( ), index.j( )); } }
Matrix<2, 2, float> strain = deformations * deformations.Transposed();
Tensor< Index4<2 /*strain i*/, 2/*strain j*/, 3/*gradient*/, 2/*point*/>, float > strainGradient; for (auto index : strainGradient.Index( )) { strainGradient.Get( index) = defPositionsInv.Get( index.l( ), index.i( )) * deformations.Get( index.k( ), index.j( )) + defPositionsInv.Get( index.l( ), index.j( )) * deformations.Get( index.k( ), index.i( )); }
Вот такие дела. Интересно, что думают другие ребята, программировавшие мотан — пробовали что-то подобное? Что в моём решении можно улучшить?
Suslik
> Начнём с простого, в формуле (5) вводится матрица из трёх трёхкомпонентных векторов-столбцов.
Строго говоря, такое тензором называть нельзя, ибо тензор — это набор чисел + закон их трансформации.
Но ладно, не суть.
> комбинированный индекс, состоящий из index.i(), index.j()
Лучше уж сразу index(0), index(1) — нету ограничения на размерность, да и все эти индексы по-любому немые, т. е. назначение буквы произвольно.
Еще хорошо бы индексы сделать типизированными, а не просто числами, чтобы была автоматическая проверка корректности.
Тут, кстати, есть сложность с настоящими тензорами, ибо там допускается суммирование между ко- и контравариантными индексами, а одинаковый тип суммировать нельзя.
Это можно разрулить, если ввести функцию суммирования по повторяющемуся индексу, но могут быть проблемы с оптимизацией.
}:+()___ [Smile]
> Лучше уж сразу index(0), index(1) — нету ограничения на размерность, да и все
> эти индексы по-любому немые, т. е. назначение буквы произвольно.
хоть это и удобнее с точки зрения программинга, но я к коду всегда прилагаю фотку с формулами на листочке, и если там индексы латиницей, то проще сравнивать.
> Еще хорошо бы индексы сделать типизированными, а не просто числами, чтобы была
> автоматическая проверка корректности.
вот это — очень крутая идея, я согласен. но я не придумал, как её реализовать без существенного усложнения кода. мне кажется, если реализовать эту штуку, это решит огромное количество проблем в духе умножения на матрицу слева или справа и с транспонированием квадратных матриц.
> Тут, кстати, есть сложность с настоящими тензорами, ибо там допускается суммирование между ко- и контравариантными индексами, а одинаковый тип суммировать нельзя.
это в принципе можно сделать типом индекса, но я не уверен, что оно того стоит.
> Это можно разрулить, если ввести функцию суммирования по повторяющемуся индексу, но могут быть проблемы с оптимизацией.
теоретически, да — все повторяющиеся индексы можно просуммировать автоматически. но пока я выбрал решение, которое склоняется к более ручному, но в то же время более контролируемому и прозрачному, где суммирование делается руками. если кто-то придумает интересное решение с автоматическим резолвингом, будет прикольно.
}:+()___ [Smile]
в общем, сделал компайлтаймовый type safety для индексов! пример использования. сначала объявляются все размерности:
struct WorldName {}; //просто пустые структуры, чтоб одну к другой привести нельзя было struct RefName {}; struct ParticleName {}; struct EdgeName {}; typedef Dim<WorldName, 3> WorldDim; //мировые координаты typedef Dim<RefName, 2> RefDim; //материальные координаты typedef Dim<ParticleName, 3> ParticleDim; //номер частицы в треугольнике typedef Dim<EdgeName, 2> EdgeDim; //номер ребра в треугольнике, перенесённого в начало координат (первая вершина = 0)
Tensor< Index2<WorldDim, RefDim>, float > deformations; //матрица деформаций. F — якобиан частных производных мировых координат(WorldDim, трёхмерные) по материальным координатам (RefDim, двумерные) deformations = ..; //инициализацию опустим Tensor< Index2<EdgeDim, RefDim>, float > defPositionsInv; //обратная матрица материальных координат, EdgeDim — столько рёбер у треугольника, RefDim — столько у них материальных координат defPositionsInv = ..; //тоже опустим Tensor< Index4<RefDim /*strain i*/, RefDim/*strain j*/, WorldDim/*gradient j*/, EdgeDim/*ref edge k*/>, float > strainGradient; //градиент тензора напряжений for (auto index : strainGradient.Index( )) { strainGradient.Get( index) = defPositionsInv.Get( index.l( ), index.i( )) * deformations.Get( index.k( ), index.j( )) + //каждый индекс i(), j(), k(), l() — типизированный и проверяется на соответствие в компайлтайме! defPositionsInv.Get( index.l( ), index.j( )) * deformations.Get( index.k( ), index.i( )); //если перепутать местами любой индекс, оно просто не скомпилится }
если писать rendering-related преобразования, например, то можно объявить мировые и view-координаты разными индексами:
struct Dim<WordName, 3> WorldDim; //world-координаты в мировой системе координат struct Dim<ViewName, 3> ViewDim; //view-координаты в системе координат камеры
Tensor<Index2<WorldName, ViewName>, float> viewMatrix;
Suslik
> если там индексы латиницей, то проще сравнивать.
А если в формуле ik без j или, вообще, αβ? А это достаточно часто, ибо разные типы индексов обозначают разными алфавитами.
Тогда уж надо как-то так:
Tensor<float, 2, 2, 3, 2> strainGradient; for(auto [i, j, k, l] : strainGradient.dimensions( )) strainGradient( i, j, k, l) = defPositionsInv( l, i) * deformations( k, j) + defPositionsInv( l, j) * deformations( k, i);
> в общем, сделал компайлтаймовый type safety для индексов! пример использования. сначала объявляются все размерности:
Я бы сделал как-то так:
struct SpatialDimension { enum { size = 3 }; static constexpr auto x = Index<SpatialDimension>(0); static constexpr auto y = Index<SpatialDimension>( 1); static constexpr auto z = Index<SpatialDimension>( 2); }; Tensor<float, SpatialDimension, SpatialDimension> matrix;
Еще я подумал, что нужны симметричные/антисимметричные варианты.
Т. е. чтобы можно было объявить симметричный тензор 3х3 с 6-ю компонентами так, чтобы (1,2) и (2,1) ссылались на одну ячейку памяти.
}:+()___ [Smile]
> strainGradient(i, j, k, l)
я использую доступ по комбинированному Get(index) вместо Get(index.i(), index.j(), index.k(), index.l()), потому что индекс внутри представлен просто как linearOffset, а его отдельные компоненты берутся через остатки от деления. поэтому вместо того, чтобы мучить оптимизатор постоянной сборкой взад-вперёд линейного индекса в его составляющие и обратно, я просто обращаюсь по линейному индексу, где это возможно. но согласен, через tuple'ы выглядит красивее. прям иногда хочется махнуть на оптимизатор и написать такое, потому что это прикольно.
}:+()___ [Smile]
> А если в формуле ik без j или, вообще, αβ?
я просто формулу на листочке записываю, используя те же индексы, что в программе :) но вообще да, с tuple'ами, где ты сам выбираешь имена индексов, прикольно. не знаю, как это реализовать без оверхеда. но tuple'ы для индексов выглядит здорово, да.
}:+()___ [Smile]
> Еще я подумал, что нужны симметричные/антисимметричные варианты.
> Т. е. чтобы можно было объявить симметричный тензор 3х3 с 6-ю компонентами так,
> чтобы (1,2) и (2,1) ссылались на одну ячейку памяти.
да я уж тоже думал об этом.. но блин, больно много кода нового получится.
реализовал метод из пейпера выше. удивительно, но он реально заработал после второй компиляции. после первой был баг, я в тензорное представление неправильно переводил, потому что не всё хранится в нём, но сами тензорные формулы работали правильно сразу, как только скомпилились. невиданно!
Suslik
> не знаю, как это реализовать без оверхеда. но прикольно, да.
Я думаю, что в большинстве случаев оптимизатор все эти циклы развернет и все индексы заменит на константные смещения.
Но чтобы наверняка, то надо убирать циклы целиком чтобы было как-то так:
TENSOR_INDEX(SpatialDimension, i); TENSOR_INDEX( SpatialDimension, j); TENSOR_INDEX( SpatialDimension, k); dst( i, j) = src1( i, k) * src2( k, j);
для сравнения с альтернативными подходами. в своё время я писал разрывный метод галёркина на произвольном функциональном базисе на неструктурированных сетках, где тензор шестого ранга тензором шестого ранга погоняет. я писал в plain-c стиле, разворачивая все суммирования руками обычными циклами, по типу:
for(int cellIndex = 0; cellIndex < 4; cellIndex++) { for( int valueIndex = 0; valueIndex < valuesCount; valueIndex++) { for( int dimIndex0 = 0; dimIndex0 < dimsCount; dimIndex0++) { for( int dimIndex1 = 0; dimIndex1 < dimsCount; dimIndex1++) { cellFlux[cellIndex].values[valueIndex] += incidentCells[cellIndex].functions[dimIndex0 + dimIndex1 * dimsCount].values[valueIndex] * orientation[cellIndex].refFlux[valueIndex]; } } } }
по опыту могу сказать, что это работает, но на уровне, что неправильный код часто выглядит странно, а не на уровне, что он не компилится. и ворочить такой код — не очень удобно.
}:+()___ [Smile]
> TENSOR_INDEX
если начинаются макросы, то всё, я даю по тормозам. ужаснота.
}:+()___ [Smile]
> Tensor<float, 2, 2, 3, 2> strainGradient;
> for(auto [i, j, k, l] : strainGradient.dimensions())
strainGradient(i, j, k, l) = defPositionsInv(l, i) * deformations(k, j) + defPositionsInv(l, j) * deformations(k, i);
попытался такой способ реализовать. надо сказать, не слишком успешно. даже объявление tuple'а таким образом:
auto [i, j, k, l] =...
SpaceDim i; SpaceDim j; RefDim k; ParticleDim k; std::tie(i, j, k, l) = strainGradient.dimensions( );
UPD: мдее, эта штука называется structured binding и появится только в C++17:
http://en.cppreference.com/w/cpp/language/structured_binding
Suslik
А от когда это , векторы и матрицы перестали быть тензорами ?
P.S: и да ты умножаешь тензор в ковариятном базисе на транспонированый тензор в коварятном базисе , при таком умножении ты получишь тензор их сумарного ранга S(ijkk), для того чтоб выполнялось правило Энштейна , а именно сумирования по повторяющюмуся индексу , чтоб на выходе был свёрнутый тензор S(ij), у тебя в уравнении с права должен быть тензор в контервариянтном базисе , ну то бишь индексы должны быть сверху, ну или наоборот тензор в коварятном базисе надо умножать на тензор контерварятном базисе , но результат будит уже другим , так как умножения не комотутивно ! Если мы говорим о метрических(геометрических) тензорах .
короче, я сдался на попытке сделать std::tuple итератор. самое близкое, что у меня получилось на pre-17 стандарте — это так:
almost::Tensor4< RefDim /*strain i*/, RefDim/*strain j*/, WorldDim/*gradient*/, EdgeDim/*edge*/, float > strainGradient; for (auto index : strainGradient.GetIndex( )) { std::make_tuple( i, j, k, l) = index.Unpack( ); strainGradient.Get( i, j, k, l) = defPositionsInv.Get( l, i) * deformations.Get( k, j) + defPositionsInv.Get( l, j) * deformations.Get( k, i); }
также заметил, что есть очень заметная разница в производительности в зависимости от того, как хранить комбинированный индекс — как линейный, разбивая его остатками от деления на индивидуальные индексы, либо хранить как индивидуальные индексы и собирать его в один линейный при доступе. так вот второй вариант генерит существенно более быстрый код, хотя я надеялся, что разницы не будет.
werasaimon
> коварятном
коварятные базисы в соседнем треде, извини.
Suslik
> самое близкое, что у меня получилось на pre-17 стандарте
Я бы просто сделал несколько вложенных циклов, по типу
for(auto i : dim<Spatial>( )) for( auto j : dim<Spatial>( )) for( auto k : dim<Spatial>( ))
> так вот второй вариант генерит существенно более быстрый код, хотя я надеялся, что разницы не будет.
Хмм... компилятор не осиливает развернуть циклы? Тогда да, умножение будет быстрее деления.
короче, много читал пейперы по мультигридам. много кода прототипировал. в итоге психанул и вернулся вот к этой затее. реализовал версию на мультииндексе. потом перевёл её на вариадик-темлпейты. посмотрел, какой код из этого генерит оптимизатор. взгрустнул. в итоге пошёл по той идее, которую смаел в последнем посте завещал:
struct Values {}; //is just used as a unique type to allow compile-time type safety using ValuesDim = almost::StaticDimension<Values, 3>; struct NullspaceVectors {}; using NullspaceVectorsDim = almost::StaticDimension<NullspaceVectors, 7>; struct Basis {}; using BasisDim = almost::DynamicDimension<Basis>; ValuesDim valuesDim; NullspaceVectorsDim nullspaceDim; BasisDim basisDim(1); //size specified in runtime auto valueToNullspace = almost::MakeStaticTensor<float, ValuesDim, NullspaceVectorsDim>( ); //on stack, static array auto nullspaceToBasis = almost::MakeDynamicTensor<float>( nullspaceDim, basisDim); //basisDim is dynamic, allocate on heap auto valuesInBasis = almost::MakeDynamicTensor<float>( valuesDim, basisDim); //basisDim is dynamic, allocate on heap for ( auto i : valuesDim) for ( auto k : basisDim) valuesInBasis.Get( i, k) = ( i.value == k.value) ? 1.0f : -42.0f; for ( auto i : valuesDim) for ( auto j : nullspaceDim) for ( auto k : basisDim) valuesInBasis.Get( i, k) += valueToNullspace.Get( i, j) * nullspaceToBasis.Get( j, k);
Suslik
> посмотрел, какой код из этого генерит оптимизатор
А чем компилировал? Не студийным ли компилятором?
MrShoor
> А чем компилировал? Не студийным ли компилятором?
если компилятор не может сгенерить нормальный набор инструкций для кода, то я это считаю проблемой кода, а не компилятора.
Тема в архиве.