skip to Main Content

The permutations of a mixed radix number can be ordered to achieve Grayness (in the sense of Gray code) with optimal balance and span length. Each of these constraints will be explained in turn. In my examples, I use a mixed radix number consisting of a base 2 digit, a base 3 digit, and a base 4 digit. This set is called [234], and it has 2 × 3 × 4 = 24 permutations. The permutations are listed below, in ascending order. For compactness, the digits are shown as rows, with the top row corresponding to the set’s first digit. The leftmost column is the first permutation 000, the next column is the second permutation 001, then 002, 003, 010, 011, 012, 013, and so on.

2:  0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1
3:  0 0 0 0 1 1 1 1 2 2 2 2 0 0 0 0 1 1 1 1 2 2 2 2
4:  0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3

In the above set, multiple digits may change from one permutation to the next. For example, between the fourth and fifth permutations (003 and 010), two digits change at once. To make a Gray set, we must reorder the permutations so that only one digit changes at a time. This constraint includes the wraparound from the first to the last permutation. Below is [234] reordered to be Gray:

2:  0 1 1 0 0 1 1 0 0 1 1 1 0 0 0 1 1 1 0 0 1 1 0 0 
3:  0 0 1 1 2 2 2 2 0 0 1 1 1 1 1 1 0 2 2 2 2 0 0 0 
4:  0 0 0 0 0 0 1 1 1 1 1 2 2 1 3 3 3 3 3 2 2 2 2 3

The above set is Gray, but not balanced. To be balanced, each of the set’s digits must change the same number of times, or as close as possible. In the above set, the 2’s place changes 10 times, the 3’s place changes 7 times, and the 4’s place also changes 7 times. A set’s imbalance is the absolute value of the difference between its minimum and maximum digit changes, in this case 10 – 7 = 3. Below is [234] reordered to have optimal balance; each digit changes 8 times, so the imbalance is now zero:

2:  0 1 1 0 0 1 1 0 0 1 1 1 0 0 0 1 1 1 1 1 0 0 0 0 
3:  0 0 1 1 2 2 2 2 0 0 1 1 1 1 1 1 2 0 0 2 2 2 0 0 
4:  0 0 0 0 0 0 1 1 1 1 1 2 2 1 3 3 3 3 2 2 2 3 3 2

The above set is Gray and balanced, but digits get stuck for longer than we’d like. For example, the 4’s place stays zero for the first six permutations. This constitutes a span, with a length of six. In the above set, the maximum span length is six. For optimal granularity, the maximum span should be as short as possible. Below is [234] reordered so that the maximum span length is four instead of six:

2:  0 1 1 0 0 1 1 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 
3:  0 0 1 1 1 1 0 0 0 0 1 2 2 2 2 1 1 1 0 2 2 2 2 0 
4:  0 0 0 0 1 1 1 1 2 2 2 2 0 0 2 2 3 3 3 3 1 1 3 3

The above ordering of [234] is Gray, optimally balanced, and minimizes stuck digits. It’s the best we can do for this particular set. But for larger sets, such as [345], optimal solutions are much harder to find, because my code is too slow. Can a MIP solver do better? The solution should be coded in a language that’s supported by one of the solvers available at NEOS, because these are the only high-quality solvers I have access to (for example CPLEX via AMPL, GAMS, LP, MPS, or NL). The application is atonal music theory, hence only sets with ranges of twelve or less are relevant. The complete list of sets I’m trying to optimize is here.

EDIT: Some commenters asked about my code, so I’m enclosing it below. I use Visual Studio 2012, but this code should compile fairly easily in any C++ compiler. I use x64 (64-bit code).

// Copyleft 2023 Chris Korda
// This program is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the Free
// Software Foundation; either version 2 of the License, or any later version.

// BalaGray.cpp : Defines the entry point for the console application.
// This app computes balanced Gray code sequences, for use in music theory.

#include "stdafx.h" // precompiled header
#include "stdint.h" // standard sizes
#include "vector"   // growable array
#include "fstream"  // file I/O
#include "assert.h" // debugging

using namespace std;

#define MORE_PLACES 0   // set non-zero to use more than four places
#define DO_PRUNING 1    // set non-zero to do branch pruning and reduce runtime

class CBalaGray {
public:
// Construction
    CBalaGray();

// Attributes
    int     GetPermCount() const { return static_cast<int>(m_arrPerm.size()); }

// Operations
    void    Reset();
    void    Calc(int nPlaces, const uint8_t *arrRange);

protected:
// Constants
    enum {
#if MORE_PLACES
        MAX_PLACES = 8,
#else
        MAX_PLACES = 4,
#endif
        MAX_RANGE = 255,
        ULONGLONG_BITS = 64,
    };
    enum {  // pruning thresholds may require manual tuning; see notes in set list
        PRUNE_MAXTRANS = 18,
        PRUNE_IMBALANCE = 3,    
    };

// Types
    union PERM {    // permutation; size depends on MAX_PLACES
        uint8_t b[MAX_PLACES];  // array of places
#if MORE_PLACES
        uint64_t    dw; // double word containing all places
#else
        uint32_t    dw; // double word containing all places
#endif
    };
    struct STATE {  // crawler stack element
        uint8_t iPerm;      // permutation index
        uint8_t iGray;      // Gray neighbor index
        PERM    nTrans;     // transition counts, one per place
    };
    typedef vector<PERM> CPermArray;
    typedef vector<STATE> CStateArray;
    typedef vector<uint8_t> CPlaceArray;    // enough for atonal music theory

// Member data
    int     m_nPlaces;  // number of places
    int     m_nGrayPerms;   // number of Gray permutations reachable from a permutation
    int     m_nGrayStrideShift; // stride of Gray permutations array, as a shift in bits
    CPlaceArray m_arrRange; // array of ranges, one for each place
    CPermArray  m_arrPerm;  // array of permutations
    CPlaceArray m_arrGray;  // 2D table of permutations reachable from each permutation
    CStateArray m_arrState; // array of states; crawler stack
    ofstream    m_fOut; // output file

// Helpers
    int     Pack(const PERM& perm) const;
    void    MakePerms(int nPlaces, const uint8_t *arrRange);
    void    MakeGrayTable();
    void    DumpGrayTablePerms() const;
    void    DumpPerm(const PERM& perm) const;
    void    DumpPerms() const;
    void    DumpSet() const;
    void    WriteBalanceToLog(int nImbalance, int nMaxTrans, int nMaxSpan);
    void    WriteSequenceToLog(int iDepth);
    bool    IsGray(PERM p1, PERM p2) const;
    int     ComputeBalance(int iDepth, int& nMaxTrans, PERM& nTransCounts) const;
    int     ComputeMaxSpan(int iDepth) const;
};

CBalaGray::CBalaGray()
{
    m_fOut.open("BalaGrayIter.txt", ios_base::out); // open output file
    assert(m_fOut != NULL);
    Reset();
}

void CBalaGray::Reset()
{
    m_nPlaces = 0;
    m_arrRange.clear();
    m_arrState.clear();
}

int CBalaGray::Pack(const PERM& perm) const
{
    int nPacked = perm.b[m_nPlaces - 1];    // init total to first place
    for (int iPlace = m_nPlaces - 2; iPlace >= 0; iPlace--) {   // for each subsequent place
        nPacked *= m_arrRange[iPlace];  // multiply total by places's range
        nPacked += perm.b[iPlace];  // add place to total
    }
    return nPacked;
}

void CBalaGray::MakePerms(int nPlaces, const uint8_t *arrRange)
{
    m_nPlaces = nPlaces;
    m_arrRange.resize(nPlaces);
    int nPerms = 1;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        assert(arrRange[iPlace] > 1);   // radix must be at least binary
        m_arrRange[iPlace] = arrRange[iPlace];  // store range
        nPerms *= arrRange[iPlace]; // update range
    }
    m_arrPerm.resize(nPerms);
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {
        PERM    perm;
        perm.dw = 0;
        int nVal = iPerm;
        for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
            int nRange = m_arrRange[iPlace];
            perm.b[iPlace] = nVal % nRange;
            nVal /= nRange;
        }
        m_arrPerm[iPerm] = perm;
    }
}

void CBalaGray::MakeGrayTable()
{
    // Build 2D table of permutations reachable from each permutation.
    // One row for each permutation, one column for each Gray neighbor.
    // Each element is a permutation index, and must be dereferenced.
    int nPlaces = m_nPlaces;
    int nGrayPerms = 0;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        nGrayPerms += m_arrRange[iPlace] - 1;   // one less than place's range
    }
    // Compute stride of Gray permutations array; to avoid multiplication,
    // round up stride to nearest power of two and convert it to a shift.
    unsigned long   iFirstBitPos;
    _BitScanReverse(&iFirstBitPos, nGrayPerms - 1);
    int nStrideShift = 1 << iFirstBitPos;
    m_arrGray.resize(m_arrPerm.size() << nStrideShift);
    int nPerms = GetPermCount();
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {  // for each permutation
        int iCol = 0;
        PERM    rowPerm, colPerm;
        rowPerm.dw = m_arrPerm[iPerm].dw;
        for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
            int nRange = m_arrRange[iPlace];    // place's range
            for (int iVal = 0; iVal < nRange; iVal++) { // for of place's values
                if (iVal != rowPerm.b[iPlace]) {    // if value differs from row value
                    colPerm.dw = rowPerm.dw;    // column permutation is same as row
                    colPerm.b[iPlace] = iVal;   // except one place differs (Gray)
                    m_arrGray[(iPerm << nStrideShift) + iCol] = Pack(colPerm);
                    iCol++; // next column
                }
            }
        }
    }
    m_nGrayPerms = nGrayPerms;  // save in member var
    m_nGrayStrideShift = nStrideShift;
}

void CBalaGray::DumpPerm(const PERM& perm) const
{
    printf("[");
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        printf("%d ", perm.b[iPlace]);
    }
    printf("]");
}

void CBalaGray::DumpPerms() const
{
    int nPerms = GetPermCount();
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {
        for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
            printf("%d ", m_arrPerm[iPerm].b[iPlace]);
        }
        printf("n");
    }
}

void CBalaGray::DumpGrayTablePerms() const
{
    int nPerms = GetPermCount();
    for (int iPerm = 0; iPerm < nPerms; iPerm++) {  // for each permutation
        DumpPerm(m_arrPerm[iPerm]);
        printf(": ");
        for (int iGray = 0; iGray < m_nGrayPerms; iGray++) {    // for each Gray neighbor
            int iPerm2 = m_arrGray[(iPerm << m_nGrayStrideShift) + iGray];
            DumpPerm(m_arrPerm[iPerm2]);
        }
        printf("n");
    }
}

void CBalaGray::DumpSet() const
{
    printf("[");
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        printf("%d", m_arrRange[iPlace]);
    }
    printf("]n");
}

void CBalaGray::WriteBalanceToLog(int nImbalance, int nMaxTrans, int nMaxSpan)
{
    m_fOut << "balance = " << nImbalance << ", maxtrans = " << nMaxTrans << ", maxspan = " << nMaxSpan << 'n';
}

void CBalaGray::WriteSequenceToLog(int iDepth)
{
    int nPerms = GetPermCount();
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        for (int iPerm = 0; iPerm < nPerms; iPerm++) {
            m_fOut << int(m_arrPerm[m_arrState[iPerm].iPerm].b[iPlace]) << ' ';
        }
        m_fOut << 'n';
    }
    m_fOut << 'n';
}

__forceinline bool CBalaGray::IsGray(PERM p1, PERM p2) const
{
    // Returns true if the given permutations differ by exactly one place.
    bool    bDiff = false;
    int nPlaces = m_nPlaces;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        if (p1.b[iPlace] != p2.b[iPlace]) { // if places differ
            if (!bDiff) {   // if first difference
                bDiff = true;   // set flag
            } else {    // not first difference
                return false;   // not Gray; early out
            }
        }
    }
    return bDiff;
}

void CBalaGray::Calc(int nPlaces, const uint8_t *arrRange)
{
    assert(nPlaces >= 0 && nPlaces <= MAX_PLACES);
    Reset();
    MakePerms(nPlaces, arrRange);
    MakeGrayTable();
//  DumpPerms();
//  DumpGrayTablePerms();
    int nPermGrays = m_nGrayPerms;
    int nGrayStrideShift = m_nGrayStrideShift;
    DumpSet();
    int nPerms = GetPermCount();
    printf("nPlaces=%dn", nPlaces);
    printf("nPerms=%dn", nPerms);
    int nBestImbalance = INT_MAX;
    int nBestMaxTrans = INT_MAX;
    int nBestMaxSpan = INT_MAX;
    m_arrState.resize(nPerms);
    uint64_t    nPasses = 0;
    uint64_t    nPermUsedMask[2] = {0}; // need 128 bits, as number of permutations may exceed 64
    int iDepth = 2; // first two levels are constant to save time; all sequences start with 0, 1
    m_arrState[1].iPerm = 1;
    m_arrState[1].nTrans.b[0] = 1;
    nPermUsedMask[0] = 0x3;
    int nStartDepth = iDepth;
    while (1) {
        nPasses++;
        int iPrevPerm = m_arrState[iDepth - 1].iPerm;
        int iGray = m_arrState[iDepth].iGray;
        int iPerm = m_arrGray[(iPrevPerm << nGrayStrideShift) + iGray]; // optimized 2D table addressing
        int iUsedMask = iPerm >= ULONGLONG_BITS;    // index selects one of two 64-bit masks
        uint64_t    nPermMask = 1ull << (iPerm & (ULONGLONG_BITS - 1));
        if (!(nPermUsedMask[iUsedMask] & nPermMask)) {  // if this permutation hasn't been used yet on this branch
            m_arrState[iDepth].iPerm = iPerm;   // save permutation index on stack
            int nMaxTrans;
            PERM    nTransCounts;
            int nImbalance = ComputeBalance(iDepth, nMaxTrans, nTransCounts);
            if (iDepth < nPerms - 1) {  // if incomplete sequence
#if DO_PRUNING
                // these constants may require tuning, see notes below
//              if (nMaxTrans > PRUNE_MAXTRANS || nImbalance > PRUNE_IMBALANCE) {   // slightly faster
                if (nImbalance > PRUNE_IMBALANCE) {
                    goto lblPrune;  // abandon this branch
                }
#endif
                // crawl one level deeper
                nPermUsedMask[iUsedMask] |= nPermMask;  // mark this permutation as used
                m_arrState[iDepth].nTrans.dw = nTransCounts.dw; // save current transition counts on stack
                iDepth++;   // increment depth to next permutation
                m_arrState[iDepth].iGray = 0;   // reset index of Gray neighbors
                m_arrState[iDepth].iPerm = 0;   // reset permutation index
                continue;   // equivalent to recursion, but less overhead
            } else {    // reached a leaf: complete sequence, a potential winner
                // if branch doesn't wrap around Gray
                if (!IsGray(m_arrPerm[m_arrState[0].iPerm], m_arrPerm[m_arrState[nPerms - 1].iPerm])) {
                    goto lblPrune;  // abandon this branch
                }
                // if max transition count or imbalance are worse than our current bests
                if (nMaxTrans > nBestMaxTrans || nImbalance > nBestImbalance) {
                    goto lblPrune;  // abandon this branch
                }
                int nMaxSpan = ComputeMaxSpan(iDepth);  // compute maximum span length
                // if max transition count and imbalance equal our current bests
                if (nMaxTrans == nBestMaxTrans && nImbalance == nBestImbalance) {
                    if (nMaxSpan >= nBestMaxSpan) { // if max span didn't improve
                        goto lblPrune;  // abandon this branch
                    }
                }
                // we have a winner, until something better comes along
                nBestMaxTrans = nMaxTrans;  // update best max transition count
                nBestImbalance = nImbalance;    // update best imbalance
                nBestMaxSpan = nMaxSpan;    // update best maximum span length
                printf("balance = %d, maxtrans = %d, maxspan = %dn", nImbalance, nMaxTrans, nMaxSpan);
                WriteBalanceToLog(nImbalance, nMaxTrans, nMaxSpan);
                WriteSequenceToLog(iDepth);
            }
        }
        m_arrState[iDepth].iGray++; // increment Gray neighbor index
        if (m_arrState[iDepth].iGray >= nPermGrays) {   // if no more Gray neighbors for this permutation
lblPrune:
            if (iDepth <= nStartDepth) {    // if we're at same level where we started
                break;  // exit main loop
            } else {    // sufficient levels remain above us
                iDepth--;   // back up a level
                // restore bitmask that keeps track of which permutations we've used on this branch
                int iPerm = m_arrState[iDepth].iPerm;   // number of permutations may exceed 64
                int iUsedMask = iPerm >= ULONGLONG_BITS;    // index selects one of two 64-bit masks
                uint64_t    nPermMask = 1ull << (iPerm & (ULONGLONG_BITS - 1));
                nPermUsedMask[iUsedMask] &= ~nPermMask; // mark this permutation as available again
                m_arrState[iDepth].iGray++; // increment was skipped by continue statement above
                if (m_arrState[iDepth].iGray >= nPermGrays) {   // if no more Gray neighbors
                    goto lblPrune;  // keep backing up
                }
            }
        }
    }
    printf("done!n");
}

__forceinline int CBalaGray::ComputeBalance(int iDepth, int& nMaxTrans, PERM& nTransCounts) const
{
    int nPlaces = m_nPlaces;
    PERM    nTrans;
    nTrans.dw = m_arrState[iDepth - 1].nTrans.dw;   // load latest transition counts from stack
    // compare current state to previous state
    PERM    sPrev, sCur;
    sPrev.dw = m_arrPerm[m_arrState[iDepth - 1].iPerm].dw;
    sCur.dw = m_arrPerm[m_arrState[iDepth].iPerm].dw;
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        if (sCur.b[iPlace] != sPrev.b[iPlace]) {    // if place transitioned
            nTrans.b[iPlace]++; // increment place's transition count
        }
    }
    nTransCounts = nTrans;  // order matters; counts passed back to caller must exclude wraparound
    // account for wraparound; compare current state to initial state, which is assumed to be zero
    for (int iPlace = 0; iPlace < nPlaces; iPlace++) {  // for each place
        if (sCur.b[iPlace]) {   // if place transitioned
            nTrans.b[iPlace]++; // increment place's transition count
        }
    }
    // now that we have latest transition counts, compute their min and max
    int nMin = nTrans.b[0]; // initialize min and max to first transition count
    int nMax = nTrans.b[0];
    for (int iPlace = 1; iPlace < nPlaces; iPlace++) {  // for each transition count, excluding first
        int n = nTrans.b[iPlace];
        if (n < nMin)   // if less than min
            nMin = n;   // update min
        if (n > nMax)   // if greater than max
            nMax = n;   // udpate max
    }
    nMaxTrans = nMax;
    return nMax - nMin; // return difference
}

__forceinline int CBalaGray::ComputeMaxSpan(int iDepth) const
{
    int arrSpan[MAX_PLACES];
    int arrFirstSpan[MAX_PLACES];
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        arrSpan[iPlace] = 1;    // initial span length is one
        arrFirstSpan[iPlace] = 0;   // first span length not set
    }
    int nMaxSpan = 1;
    PERM    sFirst, sPrev;
    sFirst.dw = m_arrPerm[m_arrState[0].iPerm].dw;  // store first state
    sPrev.dw = sFirst.dw;
    for (int iState = 1; iState <= iDepth; iState++) {  // for each state, excluding first
        PERM    s;
        s.dw = m_arrPerm[m_arrState[iState].iPerm].dw;  // compare this state to previous state
        for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
            if (s.b[iPlace] != sPrev.b[iPlace]) {   // if place transitioned
                if (arrSpan[iPlace] > nMaxSpan) // if span length exceeds max
                    nMaxSpan = arrSpan[iPlace]; // update max span length
                if (!arrFirstSpan[iPlace])  // if first span length hasn't been set
                    arrFirstSpan[iPlace] = arrSpan[iPlace]; // save first span length
                arrSpan[iPlace] = 1;    // reset span length
            } else {    // place didn't transition
                arrSpan[iPlace]++;  // increment span length
            }
        }
        sPrev = s;  // update previous state
    }
    // wrap around from last to first state
    for (int iPlace = 0; iPlace < m_nPlaces; iPlace++) {    // for each place
        if (sFirst.b[iPlace] != sPrev.b[iPlace]) {  // if place transitioned
            if (arrSpan[iPlace] > nMaxSpan) // if span length exceeds max
                nMaxSpan = arrSpan[iPlace]; // update max span length
        } else {    // place didn't transition
            arrSpan[iPlace] += arrFirstSpan[iPlace];    // compute wrapped span length
            if (arrSpan[iPlace] > nMaxSpan) // if span length exceeds max
                nMaxSpan = arrSpan[iPlace]; // update max span length
        }
    }
    return nMaxSpan;
}

void test()
{
// All cases want PRUNE_IMBALANCE = 3 unless specified otherwise below.
// Pruning greatly reduces runtime, but the results may not be optimal.
// Proven means exited normally with pruning disabled (DO_PRUNING = 0).
//
//  const uint8_t arrRange[] = {2, 10};     // proven
//  const uint8_t arrRange[] = {3, 9};
//  const uint8_t arrRange[] = {4, 8};
//  const uint8_t arrRange[] = {5, 7};
//  const uint8_t arrRange[] = {6, 6};
//  const uint8_t arrRange[] = {2, 9};      // proven
//  const uint8_t arrRange[] = {3, 8};
//  const uint8_t arrRange[] = {4, 7};
//  const uint8_t arrRange[] = {5, 6};
//  const uint8_t arrRange[] = {2, 8};      // proven
//  const uint8_t arrRange[] = {3, 7};
//  const uint8_t arrRange[] = {4, 6};
//  const uint8_t arrRange[] = {5, 5};
//  const uint8_t arrRange[] = {2, 7};      // proven
//  const uint8_t arrRange[] = {3, 6};      // proven
//  const uint8_t arrRange[] = {4, 5};      // proven
//  const uint8_t arrRange[] = {2, 6};      // proven
//  const uint8_t arrRange[] = {3, 5};      // proven
//  const uint8_t arrRange[] = {4, 4};      // proven
//  const uint8_t arrRange[] = {2, 5};      // proven
//  const uint8_t arrRange[] = {3, 4};      // proven
//  const uint8_t arrRange[] = {2, 4};      // proven
//  const uint8_t arrRange[] = {3, 3};      // proven
//  const uint8_t arrRange[] = {2, 3};      // proven
//  const uint8_t arrRange[] = {2, 2};      // proven
//  const uint8_t arrRange[] = {2, 2, 8};
//  const uint8_t arrRange[] = {2, 3, 7};
//  const uint8_t arrRange[] = {2, 4, 6};
//  const uint8_t arrRange[] = {2, 5, 5};
//  const uint8_t arrRange[] = {3, 3, 6};
//  const uint8_t arrRange[] = {3, 4, 5};
//  const uint8_t arrRange[] = {4, 4, 4};
//  const uint8_t arrRange[] = {2, 2, 7};
//  const uint8_t arrRange[] = {2, 3, 6};
//  const uint8_t arrRange[] = {2, 4, 5};
//  const uint8_t arrRange[] = {3, 3, 5};
//  const uint8_t arrRange[] = {3, 4, 4};
//  const uint8_t arrRange[] = {2, 2, 6};       // proven
//  const uint8_t arrRange[] = {2, 3, 5};
//  const uint8_t arrRange[] = {2, 4, 4};
//  const uint8_t arrRange[] = {3, 3, 4};
//  const uint8_t arrRange[] = {2, 2, 5};       // proven
//  const uint8_t arrRange[] = {2, 3, 4};       // proven
//  const uint8_t arrRange[] = {3, 3, 3};
//  const uint8_t arrRange[] = {2, 2, 4};       // proven
//  const uint8_t arrRange[] = {2, 3, 3};       // proven
//  const uint8_t arrRange[] = {2, 2, 3};       // proven
//  const uint8_t arrRange[] = {2, 2, 2};       // proven
//  const uint8_t arrRange[] = {2, 2, 2, 6};
//  const uint8_t arrRange[] = {2, 2, 3, 5};    // slow
//  const uint8_t arrRange[] = {2, 2, 4, 4};
//  const uint8_t arrRange[] = {2, 3, 3, 4};    // slow; wants PRUNE_IMBALANCE = 4
    const uint8_t arrRange[] = {3, 3, 3, 3};    // slow
//  const uint8_t arrRange[] = {2, 2, 2, 5};
//  const uint8_t arrRange[] = {2, 2, 3, 4};
//  const uint8_t arrRange[] = {2, 3, 3, 3};
//  const uint8_t arrRange[] = {2, 2, 2, 4};
//  const uint8_t arrRange[] = {2, 2, 3, 3};    // slow
//  const uint8_t arrRange[] = {2, 2, 2, 3};    // proven
//  const uint8_t arrRange[] = {2, 2, 2, 2};    // proven
//
// *** following cases require MORE_PLACES to be non-zero ***
//
//  const uint8_t arrRange[] = {2, 2, 2, 2, 4}; // wants PRUNE_IMBALANCE = 2
//  const uint8_t arrRange[] = {2, 2, 2, 3, 3}; // wants PRUNE_IMBALANCE = 4
//  const uint8_t arrRange[] = {2, 2, 2, 2, 3}; // wants PRUNE_IMBALANCE = 2
//  const uint8_t arrRange[] = {2, 2, 2, 2, 2};
//  const uint8_t arrRange[] = {2, 2, 2, 2, 2, 2};
//
    CBalaGray   bg;
    bg.Calc(_countof(arrRange), arrRange);
    fgetc(stdin);
}

int _tmain(int argc, _TCHAR* argv[])
{
    test();
    return 0;
}

2

Answers


  1. Within CPLEX, I would use CPOptimizer.

    For instance, in OPL

    // 2 2 2 2
    
    
    
    using CP;
    
    
    
    int Size=4;
    int r[1..Size]=[2,2,2,2];
    
    int States=prod(i in 1..Size) r[i];
    
    int fig[1..States][1..Size];
    
    execute
    {
      var index=0;
      for(var f1=1;f1<=r[1];f1++)
      for(var f2=1;f2<=r[2];f2++)
      for(var f3=1;f3<=r[3];f3++)
      for(var f4=1;f4<=r[4];f4++)
      {
        index++;
        fig[index][1]=f1;
        fig[index][2]=f2;
        fig[index][3]=f3;
        fig[index][4]=f4;
      }
    }
    
    dvar int x[1..States] in 1..States; // list of States in the right order
    dvar int change[1..States] in 1..Size; // the figure that is different next time
    
    dexpr int nbChanges[i in 1..Size]=count(change,i);
    dexpr int inbalance=max(i in 1..Size) nbChanges[i]-min(i in 1..Size) nbChanges[i];
    
    dvar int+ nochangeForThatManyTimes[1..States][1..Size] in 1..maxint;
    dexpr int maxspan=max(i in 1..States,j in 1..Size)  nochangeForThatManyTimes[i][j];
    
    minimize staticLex(inbalance,maxspan);
    subject to
    {
      x[1]==1;
      
      allDifferent(x);
      
      // Gray
      forall(i in 1..States,j in 1..Size) 
      ((fig[x[i]][j]==fig[x[(i<States)?(i+1):1]][j])==(j!=change[i]));
      
      
      forall(i in 2..States,j in 1..Size)  
      {
        (j==change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1);
         (j!=change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1+nochangeForThatManyTimes[i-1][j]);
      }
    
    
    
    forall(j in 1..Size)  
    (j==change[States]) ==
         (nochangeForThatManyTimes[1][j]==1);
     
       
    }
    
    execute
    {
    for(var i=1;i<=States;i++) 
    {
      for(var j=1;j<=Size;j++) write(fig[x[i]][j]-1);
      writeln();
    }  
    writeln();
    writeln("inbalance = ",inbalance);
    writeln("maxspan = ",maxspan);
    }
    

    gives

    0000
    1000
    1010
    1011
    1001
    1101
    0101
    0001
    0011
    0010
    0110
    0111
    1111
    1110
    1100
    0100
    
    inbalance = 0
    maxspan = 6
    

    and with 3,3,3,3 and a 12000 time limit I got

     OBJECTIVE: 1; 14
    0000
    1000
    2000
    2100
    2110
    2010
    2012
    0012
    0010
    0011
    2011
    1011
    1010
    1012
    1022
    1002
    0002
    0102
    0202
    0212
    0210
    2210
    2200
    2220
    0220
    0222
    1222
    2222
    2202
    2212
    2112
    1112
    0112
    0110
    0120
    0020
    0021
    0022
    0122
    0121
    0111
    0211
    0221
    1221
    1220
    1210
    1211
    1212
    1202
    1102
    1122
    1120
    1110
    1111
    1121
    1021
    1020
    2020
    2120
    2122
    2102
    2002
    2022
    2021
    2121
    2221
    2211
    2111
    2101
    2001
    1001
    0001
    0101
    1101
    1201
    2201
    0201
    0200
    1200
    1100
    0100
    
    inbalance = 1
    maxspan = 14
    

    and after 10 hours

    OBJECTIVE: 1; 10
    0000
    0001
    0002
    0012
    0010
    0011
    0111
    0121
    0221
    1221
    1021
    0021
    0022
    0222
    0122
    0102
    0112
    0212
    0202
    0201
    0101
    1101
    1001
    1011
    1012
    1112
    1110
    0110
    2110
    2010
    2020
    2000
    2002
    1002
    1022
    1020
    1120
    0120
    2120
    2121
    2221
    2222
    2202
    2200
    0200
    1200
    1000
    1010
    1210
    2210
    0210
    0211
    1211
    1201
    1202
    1102
    1100
    0100
    2100
    2101
    2102
    2112
    2012
    2022
    2122
    1122
    1121
    1111
    2111
    2011
    2021
    2001
    2201
    2211
    2212
    1212
    1222
    1220
    2220
    0220
    0020
    
    inbalance = 1
    maxspan = 10
    

    In order to get

    inbalance = 1
    maxspan = 9
    

    with

    0000
    1000
    1001
    1201
    1200
    1210
    0210
    2210
    2211
    2201
    2200
    2000
    2010
    1010
    0010
    0110
    0112
    0102
    0202
    0212
    0012
    1012
    1011
    1111
    0111
    0101
    1101
    1100
    1102
    1002
    1022
    1020
    1120
    1121
    1021
    0021
    0011
    0211
    0221
    0222
    0122
    1122
    2122
    2102
    2202
    1202
    1222
    1212
    1211
    1221
    2221
    2222
    2212
    2012
    2022
    0022
    0002
    2002
    2001
    2101
    2121
    0121
    0120
    0020
    0220
    1220
    2220
    2020
    2021
    2011
    2111
    2112
    1112
    1110
    2110
    2120
    2100
    0100
    0200
    0201
    0001
    

    I slightly improved the model

    execute
    {
      cp.param.timelimit=36000;
     
    }
    
    using CP;
    
    
    
    int Size=4;
    int r[1..Size]=[3,3,3,3];
    int maxr=max(i in 1..Size) r[i];
    
    int States=prod(i in 1..Size) r[i];
    
    int fig[1..States][1..Size];
    
    int which[i1 in 1..r[1]][i2 in 1..r[2]][i3 in 1..r[3]][i4 in 1..r[4]];
    
    execute
    {
      var index=0;
      for(var f1=1;f1<=r[1];f1++)
      for(var f2=1;f2<=r[2];f2++)
      for(var f3=1;f3<=r[3];f3++)
      for(var f4=1;f4<=r[4];f4++)
      {
        index++;
        fig[index][1]=f1;
        fig[index][2]=f2;
        fig[index][3]=f3;
        fig[index][4]=f4;
        which[f1][f2][f3][f4]=index;
      }
    }
    
    
    
    dvar int x[1..States] in 1..States; // list of States in the right order
    
    dvar int y[1..States] in 1..States;
    
    
    dvar int change[1..States] in 1..Size; // the figure that is different next time
    dvar int move[1..States] in 1..maxr;
    
    dexpr int nbChanges[i in 1..Size]=count(change,i);
    dexpr int inbalance=max(i in 1..Size) nbChanges[i]-min(i in 1..Size) nbChanges[i];
    
    dvar int+ nochangeForThatManyTimes[1..States][1..Size] in 1..maxint;
    dexpr int maxspan=max(i in 1..States,j in 1..Size)  nochangeForThatManyTimes[i][j];
    
    minimize staticLex(inbalance,maxspan);
    
    subject to
    {
      
    //   inverse(x,y);
    //  allDifferent(y);
      
      x[1]==1;
      forall(i in 1..States) move[i]<=r[change[i]];
      //inverse(x,y);
      
      change[1]==1;
      
      
      allDifferent(x);
      
      // Gray
    //  forall(i in 1..States,j in 1..Size) 
    //  {
    //  (j!=change[i]) == (fig[x[i]][j]==fig[x[(i<States)?(i+1):1]][j]);
    //  (j==change[i]) == ((fig[x[i]][j]+move[i]-1) mod r[j]+1==fig[x[(i<States)?(i+1):1]][j]);
    //}  
     
     forall(i in 1..States)
       
       x[(i<States)?(i+1):1]
       ==which
       [(fig[x[i]][1]+(1==change[i])*move[i]-1) mod r[1]+1]
       [(fig[x[i]][2]+(2==change[i])*move[i]-1) mod r[2]+1]
       [(fig[x[i]][3]+(3==change[i])*move[i]-1) mod r[3]+1]
       [(fig[x[i]][4]+(4==change[i])*move[i]-1) mod r[4]+1]
       ;
       
       
     
      
      
      inferred(change);
      inferred(move);
      inferred(nochangeForThatManyTimes);
     
     
      inbalance>=States mod 2;
      
      forall(i in 2..States,j in 1..Size)  
      {
        (j==change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1);
         (j!=change[i-1]) =>   (nochangeForThatManyTimes[i][j]==1+nochangeForThatManyTimes[i-1][j]);
      }
    
    
    
    forall(j in 1..Size)  
    (j==change[States]) ==
         (nochangeForThatManyTimes[1][j]==1);
     
       
    }
    
    execute
    {
    for(var i=1;i<=States;i++) 
    {
      for(var j=1;j<=Size;j++) write(fig[x[i]][j]-1);
      writeln();
    }  
    writeln();
    writeln("inbalance = ",inbalance);
    writeln("maxspan = ",maxspan);
    }
    

    NB: You can use CPLEX for free in the cloud with this OPL API

    Login or Signup to reply.
  2. Well, after some tinkering, I have an Integer Program running for this, that I think is producing quality results. Tried a couple approaches…each had differing limitations

    It is a little grotesque in parts as the counting of repeat digits is quite cumbersome.

    It really bogs down for things with ~30 states or more, so it’s not going to make it to the finish line. 🙂 I think it is much more nimble if I remove the repeat counting, and I’ll tinker a bit more. In the interim, here are some results for the cases not marked as proven on your web page. The (4, 6) run (second run) is an improvement, the other 2 are now "proven" as stated, perhaps with a different sequence, I didn’t x-check.

    I’ll update later with any other improvements.

    starting run:  (3, 7)
    WARNING: Initializing ordered Set PR_flat with a fundamentally unordered data
        source (type: set).  This WILL potentially lead to nondeterministic
        behavior in Pyomo
    
    Problem: 
    - Name: unknown
      Lower bound: 1.3
      Upper bound: 1.3
      Number of objectives: 1
      Number of constraints: 3744
      Number of variables: 3539
      Number of binary variables: 3557
      Number of integer variables: 3560
      Number of nonzeros: 3
      Sense: minimize
    Solver: 
    - Status: ok
      User time: -1.0
      System time: 353.01
      Wallclock time: 305.85
      Termination condition: optimal
      Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
      Statistics: 
        Branch and bound: 
          Number of bounded subproblems: 1018
          Number of created subproblems: 1018
        Black box: 
          Number of iterations: 667719
      Error rc: 0
      Time: 306.00031781196594
    Solution: 
    - number of solutions: 0
      number of solutions displayed: 0
    
    11
    01
    02
    00
    10
    20
    21
    22
    12
    13
    23
    03
    05
    25
    26
    24
    14
    04
    06
    16
    15
    max imbalance: 1.0
    max repeats: 3.0
    starting run:  (4, 6)
    WARNING: Initializing ordered Set PR_flat with a fundamentally unordered data
        source (type: set).  This WILL potentially lead to nondeterministic
        behavior in Pyomo
    
    Problem: 
    - Name: unknown
      Lower bound: 0.2
      Upper bound: 0.2
      Number of objectives: 1
      Number of constraints: 4854
      Number of variables: 4619
      Number of binary variables: 4640
      Number of integer variables: 4643
      Number of nonzeros: 3
      Sense: minimize
    Solver: 
    - Status: ok
      User time: -1.0
      System time: 34.21
      Wallclock time: 34.89
      Termination condition: optimal
      Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
      Statistics: 
        Branch and bound: 
          Number of bounded subproblems: 1
          Number of created subproblems: 1
        Black box: 
          Number of iterations: 14167
      Error rc: 0
      Time: 34.923232078552246
    Solution: 
    - number of solutions: 0
      number of solutions displayed: 0
    
    10
    13
    33
    34
    14
    15
    35
    32
    02
    03
    23
    22
    12
    11
    01
    00
    30
    31
    21
    24
    04
    05
    25
    20
    max imbalance: 0.0
    max repeats: 2.0
    starting run:  (5, 5)
    WARNING: Initializing ordered Set PR_flat with a fundamentally unordered data
        source (type: set).  This WILL potentially lead to nondeterministic
        behavior in Pyomo
    
    Problem: 
    - Name: unknown
      Lower bound: 1.3
      Upper bound: 1.3
      Number of objectives: 1
      Number of constraints: 5256
      Number of variables: 5011
      Number of binary variables: 5033
      Number of integer variables: 5036
      Number of nonzeros: 3
      Sense: minimize
    Solver: 
    - Status: ok
      User time: -1.0
      System time: 915.71
      Wallclock time: 634.99
      Termination condition: optimal
      Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
      Statistics: 
        Branch and bound: 
          Number of bounded subproblems: 1764
          Number of created subproblems: 1764
        Black box: 
          Number of iterations: 1855323
      Error rc: 0
      Time: 635.0473001003265
    Solution: 
    - number of solutions: 0
      number of solutions displayed: 0
    
    11
    01
    31
    33
    34
    44
    04
    03
    00
    40
    41
    42
    22
    23
    43
    13
    12
    02
    32
    30
    10
    20
    21
    24
    14
    max imbalance: 1.0
    max repeats: 3.0
    
    Login or Signup to reply.
Please signup or login to give your own answer.
Back To Top
Search