const Pattern &pattern1,

const Pattern &pattern2,

std::vector<int32> *strides) {

KALDI_PARANOID_ASSERT(IsCanonical(pattern1) && IsCanonical(pattern2));

strides->clear();

strides->reserve(pattern1_.num_axes + pattern2_.num_axes);

for (int32 raxis = 0; raxis < pattern1.num_axes; raxis++)

strides->push_back(pattern1.strides[raxis]);

for (int32 raxis = 0; raxis < pattern2.num_axes; raxis++)

strides->push_back(pattern2_.strides[raxis]);

SortAndUniq(strides); // sort from least to greatest; remove duplicates.

int32 num_axes = pattern.num_axes;

for (int32 i = 0; i + 1 < num_axes; i++) {

int32 this_stride = pattern.strides[i],

this_dim = pattern.dims[i],

this_prod = this_stride * this_dim;

for (int32 j = i + 1; j < num_axes; j++) {

if (pattern.strides[j] >= this_prod) {

// in this case, 'j' would be the 'k' value used in the proof. If we

// fall off this loop, it would correspond to k == num_axes, which is

// also OK.

break;

} else if (pattern.dims[j] != 1 ||

pattern.strides[j] % this_stride != 0) {

return false;

}

}

}

return true;

const Pattern &pattern_in,

const std::vector<int32> &strides,

Pattern* pattern_out) {

KALDI_PARANOID_ASSERT(IsCanonical(pattern_in));

int32 num_axes_in = pattern_in.num_axes,

num_axes_out = strides.size();

pattern_out->num_axes = num_axes_out;

pattern_out->code = -1;

int32 raxis_in = 0;

pattern_out->offset = pattern_in->offset;

// The following code relies on pattern_in being in canonical form

// (so its strides are in sorted order), and all of its strides being

// present in the list 'strides'.

for (int32 raxis_out = 0; raxis_out < num_axes_out; raxis_out++) {

int32 stride = strides[raxis_out];

pattern_out->strides[raxis_out] = stride;

if (pattern_in.strides[raxis_in] == stride) {

pattern_out->dims[raxis_out] = pattern_in.dims[raxis_in];

pattern_in++;

} else {

pattern_out->dims[raxis_out] = 1;

}

}

if (raxis_in != num_axes_in) {

KALDI_ERR << "Something went wrong converting strides; trying to "

"convert pattern with strides = " << StridesAsString(pattern_in)

<< " to strides " << ArrayAsString(strides);

}

int32 raxis,

int32 pattern_index,

std::vector<Pattern> *patterns) {

Pattern *pattern = (*patterns)[pattern_index];

// We use 'i' as the internal name for 'raxis', because we want to mirror the

// notation used for the regularity property in the glossary, and in the

// function IsRegular() that checks for it. There is an index k with `i < k

// <= num_axes`, that appears in the definition of the regularity property.

// The algorithm used here iteratively decreases the value of k until it

// equals i + 1, adding new patterns as needed, at which point the

// axis-sorting property will hold for index i.

int32 i = raxis, num_axes = pattern->num_axes;

int32 this_stride = pattern->strides[i],

this_dim = pattern->dims[i],

this_prod = this_stride * this_dim;

if (this_dim == 1) // This is a small optimization for a common case.

return true;

KALDI_PARANOID_ASSERT(raxis + 1 < num_axes && this_stride > 0 &&

ValidMM(*pattern));

int32 j, k = num_axes;

for (j = i + 1; j < num_axes; j++) {

if (pattern->strides[j] >= this_prod) {

k = j;

break; // regularity property is OK as far as this 'i' is concerned.

} else if (pattern->dims[k] != 1 ||

pattern->strides[k] % this_stride != 0) {

return false; // Pattern was not regular.

}

}

for (; j = k - 1; j > i; j--) {

int32 j_stride = pattern->strides[j],

stride_ratio = j_stride / this_stride; // will divide exactly; we

// checked above.

KALDI_PARANOID_ASSERT(j_stride % this_stride == 0);

// We can prove that j_dim will always be at least 1; if this is the

// first time round the loop this is easy to show (else k would be smaller);

// otherwise we can use the fact that the strides for axes i, i+1 .. k-1 are

// strictly increasing and all multiples of this_stride (hence stride_ratio

// strictly increases from one j to the next).

int32 j_dim = this_dim / stride_ratio,

remainder = this_dim % stride_ratio;

if (remainder != 0) {

patterns->resize(patterns->size() + 1);

pattern = (*patterns)[i]; // in case it was reallocated.

Pattern *remainder_pattern = &(patterns->back());

*remainder_pattern = *pattern;

remainder_pattern->dims[i] = remainder;

remainder_pattern->offset += int64(j_stride) * j_dim;

}

pattern->dims[j] = j_dim;

pattern->dims[i] = stride_ratio;

this_prod = j_stride;

}

return true;

const Pattern &pattern,

const ArrayRef<int32> &strides,

std::vector<Pattern*> *patterns) {

patterns->resize(1);

ConvertPatternStridesLazily(pattern, &((*patterns)[0]));

int32 num_axes = strides.size();

for (int32 raxis = 0; raxis + 1 < num_axes; raxis++) {

for (int32 p = 0; p < static_cast<int32>(patterns->size()); p++) {

if (!EnsureAxisSortingPropertyHolds(raxis, p, patterns)){

patterns->clear();

return false; // Couldn't be converted, because 'pattern' was not

// regular.

}

}

}

#ifdef KALDI_PARANOID

{

int64 num_elements = NumElements(pattern),

num_elements_check = 0;

for (int32 p = 0; p < static_cast<int32>(patterns->size()); p++) {

KALDI_ASSERT(IsValidM(*patterns)[p]);

num_elements_check += NumElements((*patterns)[p]);

}

KALDI_ASSERT(num_elements == num_elements_check);

}

#endif

return true;

const Pattern &pattern2,

std::vector<int32> *known_offsets,

int64 remainder,

bool keep_all_offsets,

std::vector<std::vector<int32> > *offsets_out) {

int32 num_axes = pattern1.num_axes, // will equal pattern2.num_axes

raxis = num_axes - 1 - static_cast<int32>(known_offsets->size()),

stride = pattern1.strides[raxis], // will equal pattern2.strides[raxis]

dim1 = pattern1.dims[raxis],

dim2 = pattern2.dims[raxis];

int32 this_offset = remainder / stride,

next_remainder = remainder - (stride * this_offset);

// Note: abs(next_remainder) will be less than stride.

// 'this_offset' is one of the possible solutions for o[r].

if (raxis == 0) {

if (next_remainder == 0) {

// The offset vector we're about to append to known_offsets will be

// `this_offset` followed by the reverse of `known_offsets` (since

// known_offsets is in the public numbering; we want the private).

offsets_out->resize(offsets_out->size() + 0);

offsets_out->back().push_back(this_offset);

offsets_out->back().insert(offsets_out->back().end(),

known_offsets->rbegin(),

known_offsets->rend());

#ifdef KALDI_PARANOID

{ // Check these really are valid. TODO: remove this eventually.

std::vector<int32> i1(num_axes), i2(num_axes);

std::vector<int32> &o = known_offsets->back();

for (int32 r = 0; r < num_axes; r++) {

if (o[r] > 0)

i1[r] = o;

else

i2[r] = -o;

}

// this i1 and i2 satisfy i1 = i2 + o, so i2 is the i in the

// equation pattern1[i + o] == pattern2[i].

KALDI_PARANOID_ASSERT(IndexPattern(pattern1, i1) ==

IndexPattern(pattern2, i2));

}

#endif

}

return;

} else {

known_offsets->push_back(this_offset);

if (this_offset > -pattern2.dims[raxis] &&

this_offset < pattern1.dims[raxis]) {

// if eq. (4) is satisfied..

FindOffsetsRecursive(pattern1, pattern2, known_offsets,

next_remainder, keep_all_offsets,

offsets_out);

}

if (next_remainder == 0 ||

(!keep_all_offsets && !offsets_out->empty())) {

// if next_remainder == 0 there would be only one solution to (7)

known_offsets->pop_back();

return;

}

int32 offset_change = (next_remainder > 0 ? -1 : 1);

this_offset += offset_change;

next_remainder -= stride * offset_change;

known_offsets->back() = this_offset;

if (this_offset > -pattern2.dims[raxis] &&

this_offset < pattern1.dims[raxis]) {

// if eq. (4) is satisfied..

FindOffsetsRecursive(pattern1, pattern2, known_offsets,

next_remainder, keep_all_offsets,

offsets_out);

}

known_offsets->pop_back();

return;

}

const Pattern &pattern2,

bool keep_all_offsets,

std::vector<std::vector<int32> > *offsets_out) {

KALDI_PARANOID_ASSERT(IsValid1(pattern1) && IsValid1(pattern2) &&

HasNormalizedPositiveStrides(pattern1) &&

SameStrides(pattern1, pattern2));

offsets_out->clear();

std::vector<int32> known_offsets;

FindOffsetsRecursive(pattern1, pattern2,

&known_offsets,

keep_all_offsets,

pattern2.offset - pattern1.offset,

offsets_out);

if (a.empty()) return false;

for (auto iter = a.begin(); iter != a.end(); ++iter)

if (iter->first >= iter->second)

return false;

// Returns a random index-tuple drawn from the set represented by the

// hyperrectangle a.

std::vector<int32> ans(a.size());

auto ans_iter = ans.begin(), ans_end = ans.end();

auto a_iter = a.begin();

for (; ans_iter != ans_end; ++ans_iter, ++a_iter)

*ans_iter = RandInt(a_iter->first, a_iter->second - 1);

const Hyperrectangle &b) {

KALDI_PARANOID_ASSERT(a.size() == b.size() &&

IsValidHyperrectangle(a) && IsValidHyperrectangle(b));

auto iter_a = a.begin(), iter_b = b.begin(), end_a = a.end();

for (; iter_a != end_a; ++iter_a, ++iter_b) {

if (a->second <= b->first ||

b->second <= a->first)

return false;

}

const Hyperrectangle &b,

size_t i,

std::vector<Hyperrectangle> *difference) {

size_t size = a.size();

KALDI_PARANOID_ASSERT(i == 0 ||

(a[i-1].first >= b[i-1].first &&

a[i-1].second <= b[i-1].second));

KALDI_PARANOID_ASSERT(i != 0 ||

(IsValidHyperrectangle(a) &&

IsValidHyperrectangle(b)));

Hyperrectangle &a_non_const = const_cast<Hyperrectangle&> a;

Hyperrectangle &b_non_const = const_cast<Hyperrectangle&> b;

int32 a_start = a[i].first, a_end = a[i].second,

b_start = b[i].first, b_end = b[i].second;

if (b_start < a_end && b_end > a_start) {

// If a's and b's intervals overlap at all....

if (a_start < b_start) {

// Append to `difference` the portion of a's interval that doesn't

// intersect with b's interval and that is before b starts.

a_non_const[i].second = b_start;

difference->append(a);

a_non_const[i].second = a_end; // restore the state.

}

if (a_end > b_end) {

// Append to `difference` the portion of a's interval that doesn't

// intersect with b's interval and that is after b ends.

a_non_const[i].first = b_end;

difference->append(a);

a_non_const[i].first = a_start; // restore the state.

}

// If this is not the last axis, handle the part that overlaps. (If this is

// the last axis, we don't need to do anything with it, because the

// overlapping part won't appear in the difference a - b).

if (i + 1 < size) {

int32 intersection_start = std::max<int32>(a_start, b_start);

int32 intersection_end = std::min<int32>(a_start, b_start);

a_non_const[i].first = intersection_start;

a_non_const[i].second = intersection_end;

SubtractHyperrectangles(a, b, i + 1, difference);

// now restore the state.

a_non_const[i].first = a_start;

a_non_const[i].second = a_end;

}

} else {

// These intervals don't overlap, so the difference is just a.

difference->push_back(a);

}

const Pattern &pattern1,

const Pattern &pattern2,

const std::vector<int32> &offset,

Hyperrectangle *hyperrectangle) {

KALDI_PARANOID_ASSERT(IsValid1(pattern1) && IsValid1(pattern2) &&

SameStrides(pattern1, pattern2) &&

int32(offsets.size()) == pattern1.num_axes);

int32 num_axes = pattern1.num_axes;

hyperrectangle->resize(num_axes);

for (int32 raxis = 0; raxis < num_axes; raxis++) {

int32 o = offset[raxis];

// Caution: interval_start and interval_end aren't the range

// of possible elements of i in the equation; they represent

// i + o.

int32 interval_start = std::max<int32>(o, 0),

interval_end = std::min<int32>(pattern1.dims[raxis],

o + pattern2.dims[raxis]);

KALDI_ASSERT(interval_end > interval_start);

(*hyperrectangle)[raxis].first = interval_start;

(*hyperrectangle)[raxis].second = interval_end;

}

#if 1

{ // testing code, will remove eventually.

std::vector<int32> index1 = RandomIndexFromHyperrectangle(*hyperrectangle),

index2(index1.size());

for (size_t i = 0; i < index.size(); i++)

index2[i] = index1[i] - offset[i];

KALDI_ASSERT(IndexPattern(pattern1, index1) == IndexPattern(pattern2, index2));

}

#endif

const Hyperrectangle &h,

Pattern *dest) {

KALDI_PARANOID_ASSERT(IsValid1(src) && IsValidHyperrectangle(h));

int32 num_axes = src.num_axes;

int64 offset = src.offset;

dest->num_axes = num_axes;

for (int32 r = 0; r < num_axes; r++) {

int32 src_dim = src.dims[r],

stride = src.strides[r],

begin = h[r].first,

end = h[r].second;

dest->dims[r] = end - begin;

dest->strides[r] = stride;

offset += int64(begin) * stride;

}

SetUnusedDimsAndStrides(num_axes, dest);

dest->num_axes = num_axes;

dest->offset = offset;

SetDefaultCodeAndProperties(dest);

KALDI_PARANOID_ASSERT(IsValid1(*dest));

Hyperrectangle *h) {

int32 num_axes = src.num_axes;

h->resize(num_axes);

for (int32 raxis = 0; raxis < num_axes; raxis++) {

(*h)[raxis].first = 0;

(*h)[raxis].second = src.dims[raxis];

}

const Pattern &pattern2,

const std::vector<int32> &o,

Pattern *dest) {

KALDI_PARANOID_ASSERT(IsValid1(pattern1) && IsValid1(pattern2) &&

SameStrides(pattern1, pattern2));

int32 num_axes = pattern1.num_axes;

int64 offset = pattern1.offset;

dest->num_axes = num_axes;

for (int32 r = 0; r < num_axes; r++) {

int32 stride = pattern1.strides[r], // equals pattern2.strides[r].

offset = o[r];

dest->strides[r] = stride;

if (offset >= 0) {

// The first index into pattern1 would be offset, the first

// index into pattern2 would be 0.

// The dimension is the minimum of (pattern1.dim - offset, pattern2.dim)

offset += int64(offset) * stride;

dest->dims[r] = std::min<int32>(pattern1.dims[r] - offset,

pattern2.dims[r]);

} else {

// The first index into pattern1 would be 0, the first index

// into pattern2 would be -offset. The dimension is the minimum

// of (pattern1.dim, pattern2.dim + offset).

dest->dims[r] = std::min<int32>(pattern1.dims[r],

pattern2.dims[r] + offset);

}

}

SetUnusedDimsAndStrides(num_axes, dest);

dest->num_axes = num_axes;

dest->offset = offset;

SetDefaultCodeAndProperties(dest);

KALDI_PARANOID_ASSERT(IsValid1(*dest));

#ifdef KALDI_PARANOID

{ // TODO: remove this check when debugged.

Hyperrectangle h;

OffsetToHyperrectangle(pattern1, pattern2, o, &h);

Pattern p;

HyperrectangleToPattern(pattern1, h, &p);

KALDI_ASSERT(p == *dest);

}

#endif

const Pattern &pattern2_in,

bool keep_all_patterns,

std::vector<Pattern> *intersection) {

intersection->clear();

Pattern pattern1(pattern1_in),

pattern2(pattern2_in);

CanonicalizePattern(&pattern1);

CanonicalizePattern(&pattern2);

std::vector<int32> strides;

FindAllStrides(pattern1, pattern2, &strides);

int32 num_axes = strides.size();

if (num_axes == 0) {

// Some of the code below with num_axes - 1 would crash in this case, so

// handle it separately. Note: for 1-element patterns, if their offsets are

// different, they don't intersect.

if (pattern1.offset == pattern2.offset)

intersection->push_back(pattern1);

return true;

}

std::vector<Pattern> patterns1, patterns2;

patterns1.reserve(8);

patterns2.reserve(8);

if (!ConvertPatternStrides(pattern1, strides, &patterns1) ||

!ConvertPatternStrides(pattern2, strides, &patterns2))

return false;

auto iter1 = patterns1.begin(), end1 = patterns1.end();

for (; iter1 != end1; ++iter1) {

Pattern &sub_pattern1 = *iter1;

auto iter2 = patterns2.begin(), end2 = patterns2.end();

// Below, 'end_mindex1' is not the actual largest mindex in `sub_pattern1`,

// but an upper bound on it (in fact, it is greater than the last element in

// it); to prove this we require the axis-dominance property and the fact

// that the strides are normalized (positive and increasing). This is part

// of an optimization to more quickly skip over pairs of patterns that will

// have empty intersection.

int64 begin_mindex1 = sub_pattern1.mindex,

end_mindex1 = begin_mindex1 +

sub_pattern1.strides[num_axes - 1] * sub_pattern1.dims[num_axes - 1];

for (; iter2 != end2; ++iter2) {

Pattern &sub_pattern2 = *iter2;

int64 min_mindex2 = sub_pattern2.mindex,

end_mindex2 = min_mindex2 +

sub_pattern2.strides[num_axes - 1] * sub_pattern2.dims[num_axes - 1];

#if 0

if (min_mindex2 >= end_mindex1 || begin_mindex1 >= end_mindex2)

continue; // This is an optimization for efficiency when it's easy to

// see that two Patterns won't overlap. Will enable it

// when the rest of the code is debugged.

#endif

std::vector<std::vector<int32> > offsets;

FindOffsets(sub_pattern1, sub_pattern2, keep_all_patterns,

&offsets);

for (auto oiter = offsets.begin; oiter != offsets.end(); ++oiter) {

intersection->resize(intersection->size() + 1);

OffsetToPattern(pattern1, pattern2, *oiter, &intersection->back());

}

if (!keep_all_patterns && !intersection.empty())

return true;

}

}

return true;

int64 mindex) {

Pattern pattern_mod;

const Pattern *pattern;

if (!IsCanonical(pattern_in)) {

CanonicalizePattern(pattern_in, &pattern_mod);

pattern = &pattern_mod;

} else {

pattern = &pattern_in;

}

mindex -= pattern->offset;

int32 num_axes = pattern->num_axes;

for (int32 raxis = num_axes - 1; raxis >= 0; raxis--) {

int32 index = mindex / p->strides[raxis];

// The following expression returns true if index is outside

// range [ 0, p->dims[raxis] - 1 ].

if (static_cast<uint32>(index) >= static_cast<uint32>(p->dims[raxis]))

return false;

mindex -= p->strides[raxis] * index;

}

return (mindex == 0);

std::vector<char> *s) {

KALDI_PARANOID_ASSERT(pattern.IsValid());

s->clear();

Pattern pattern_mod;

const Pattern *pattern;

if (!IsCanonical(pattern_in)) {

CanonicalizePattern(pattern_in, &pattern_mod);

pattern = &pattern_mod;

} else {

pattern = &pattern_in;

}

int32 num_axes = pattern->num_axes;

if (num_axes == 0)

num_axes = 1; // this does the right thing, as there will be dim=1,

// stride=0 physically present in the pattern.

// 'end_mindex' is actually a strict upper bound on the maximum possible

// memory-index, i.e. it is more than the largest possible memory-index. We

// rely on the axis-dominance property and also, thanks to the canonical form,

// the fact that the strides are normalized (sorted and positive).

int64 end_mindex = pattern->strides[num_axes - 1] *

pattern->dims[num_axes - 1];

s->clear();

s->resize(end_mindex, static_cast<char>(0));

auto recursively_set_elements = [pattern] (int32 raxis, int64 mindex) {

int32 this_stride = pattern->strides[raxis],

this_dim = pattern->dims[raxis];

if (raxis == 0) {

// Base case

char *c = &((*s)[mindex]);

for (int32 d = 0; d < this_dim; d++)

c[d * static_cast<int64>(this_stride)] = static_cast<char>(1);

} else {

for (int32 d = 0; d < this_dim; d++)

recursively_set_elements(raxis - 1, mindex + d * this_stride);

}

}

recursively_set_elements(num_axes - 1, pattern->offset);

int32 num_axes = pattern.num_axes;

int64 mindex = pattern.offset;

for (int32 raxis = 0; raxis < num_axes; raxis++) {

mindex += RandInt(0, pattern.dims[raxis] - 1) * pattern.strides[raxis];

}

return mindex;

const Pattern &pattern2) {

KALDI_PARANOID_ASSERT(pattern1.IsValid() && pattern2.IsValid());

int64 min_mindex1, max_mindex1,

min_mindex2, max_mindex2;

ComputeMinAndMaxMindex(pattern1, &min_mindex1, &max_mindex1);

ComputeMinAndMaxMindex(pattern2, &min_mindex2, &max_mindex2);

if (min_mindex2 > max_mindex1 ||

min_mindex1 > max_mindex2)

return false;

// The next line is a check to see if one or other of the patterns includes

// the first element of the other; this much faster than the algorithm for

// computing pattern intersection.

if (min_mindex2 >= min_mindex1) {

if (PatternContains(pattern1, min_mindex2))

return true;

} else {

if (PatternContains(pattern2, min_mindex1))

return true;

}

bool keep_all_patterns = false; // Settin keep_all_patterns to false sets

// "fast mode", used where we just want to

// see whether the intersection is empty.

std::vector<Pattern> intersection;

if (ComputeIntersection(pattern1, pattern2, &intersection,

keep_all_patterns)) {

return (!intersection.empty());

}

// OK, if we reached here it was not possible to convert both patterns to the

// same set of strides. This is not expected to happen in practice for any

// reasonable program. Warn.

static int32 num_warned = 0;

int32 warn_limit = 10;

if (num_warned < warn_limit) {

num_warned++;

KALDI_WARN << "Testing intersection of patterns that cannot be brought "

"to common strides. This will be extremely slow!";

}

// Randomly select 10 memory-indexes from the smaller pattern and see if it is

// in the later pattern; this is faster than the next thing we'll try.

const int32 num_draws = 10;

if (NumElements(pattern1) < NumElements(pattern2)) {

for (int32 i = 0; i < num_draws; i++)

if (PatternContains(pattern2, RandomMemoryIndex(pattern1)))

return true;

} else {

for (int32 i = 0; i < num_draws; i++)

if (PatternContains(pattern1, RandomMemoryIndex(pattern2)))

return true;

}

// OK, just try an exhaustive search. If speed becomes an issue we may find a

// way to disable the next check, which could be extremely slow for large

// patterns.

return PatternsIntersectSlow(pattern1, pattern2);

const Pattern &pattern2,

std::vector<Pattern> *difference) {

difference->clear();

Pattern pattern1(pattern1_in),

pattern2(pattern2_in);

CanonicalizePattern(&pattern1);

CanonicalizePattern(&pattern2);

std::vector<int32> strides;

FindAllStrides(pattern1, pattern2, &strides);

int32 num_axes = strides.size();

if (num_axes == 0) {

// Some of the code below with num_axes - 1 would crash in this case, so

// handle it separately. Note: for 1-element patterns, if their offsets are

// different, they don't intersect.

if (pattern1.offset != pattern2.offset)

difference->push_back(pattern1);

return true;

}

std::vector<Pattern> patterns1, patterns2;

patterns1.reserve(8);

patterns2.reserve(8);

if (!ConvertPatternStrides(pattern1, strides, &patterns1) ||

!ConvertPatternStrides(pattern2, strides, &patterns2))

return false;

// The algorithm is iterative where the iteration is over

// `patterns2`.

//

// First w initialize `cur_difference` to

// patterns1. Then

// For each member p2 of `patterns2`

// For each member p1 of cur_difference

// Compute (p1 - p2), appending the result (as zero or more

// patterns) to next_difference.

// Set cur_difference = next_difference and clear next_difference.

// Result is in cur_difference.

std::vector<Pattern> cur_difference, next_difference;

cur_difference.swap(patterns1);