Rajoitettu vaihteluväliset kävelyt

Tutkitaan $n$ :n pituisia jonoja symboleista $-1$ ja $1$ . Kun tällainen mielletään kävelyksi, määräytyy sille kunkin askeleen jälkeen korkeus, jolla kävely on (lähdetään nollasta). Kävelylle $w$ lasketaan vaihteluväli $h$ erotuksena $h(w) = \max_j h_j(w) - \min_j h_j(w)$ , missä $h_j$ on kävelyn korkeus $j$ askeleen jälkeen.

Kaksi esimerkkikävelyä ja niiden vaihteluvälin lasku. Tässä $n=7$ .

Kaksi esimerkkikävelyä ja niiden vaihteluvälin lasku. Tässä $n=7$ .

Kysymys: kuinka monta kävelyistä on vaihteluväliltään korkeintaan $r$ , missä $r \geq 0$ on annettu parametri.

Ratkaisu

Tehdään ensin kokeilua varten algoritmi, joka käy kaikki kävelyt (s.o. -1, 1 -jonot) läpi. Kielenä SageMath

def brw(r, n):
    count = 0
    for w in cartesian_product([(-1,1)]*n):
        ss = [sum(w[:j]) for j in range(n+1)]
        if max(ss)-min(ss) <= r: count += 1
    return count

for r in range(7): print ([brw(r, n) for n in range(10)])

Huomataan, että arvoille r = 2, r = 3, r = 4 ja r = 5 jono löytyy OEIS:sta. Sieltä löytyy generoivat funktiot ja tapaukselle $r=2$ kaava

brw(2, n) = $2^{\lfloor \frac{n+2}{2} \rfloor} + 2^{\lfloor \frac{n+1}{2} \rfloor} - 2$

Myös tapauksille $r=3,4$ on OEIS:ssa lähes yhtä siistit kaavat (niissä esiintyy Fibonacci-luvut).

Koska emme kuitenkaan keksi tästä yleistä kaikille $r$ toimivaa kaavaa, ryhdytään kehittämään tehokkaampaa algoritmia. Kun jonoa $w$ luetaan alusta lähtien ja tarkastellaan korkeutta, jolla ollaan askeleen $j$ jälkeen, niin mieleen tulee tietenkin ~~Markovin ketju~~ suunnattu verkko, jossa ~~tilana~~ solmuina on $h_j$ . Tämä ei kuitenkaan itsessään riitä, vaan tilaan tarvitaan myös tieto siitä missä maksimi- ja minimikorkeuksissa on käyty siihen mennessä. Siirtymät tietenkin tehdään vain sallituissa rajoissa. Seuraava Sage-koodi tekee tämän ja laskee kävelyjen määrän siirtymämatriisin $n$ . potenssista (summataan lähtötilaa vastaava rivi).

def brw(r, n):
    bds = [(i,j) for i in range(-r, 1) for j in range(r+1) if j-i<=r]
    states = [] #triples (i, j, k) where (i,j) are the min, max -heights so far and k is current height
    for (i,j) in bds:
        states += [(i,j,k) for k in range(i, j+1)]
    matN = len(states)
    startState = states.index((0,0,0))
    arr = [[0]*matN for _ in range(matN)]
    for ind1, (i,j,k) in enumerate(states):
        i2 = i-1 if k-1<i else i
        if j-i2<=r:
            s2 = (i2, j, k-1)
            arr[ind1][states.index(s2)] += 1
        j2 = j+1 if k+1>j else j
        if j2-i<=r:
            s2 = (i, j2, k+1)
            arr[ind1][states.index(s2)] += 1
    return sum((Matrix(arr)**n)[startState])

Tilojen määrä kuitenkin kasvaa nopeasti, kun r kasvaa ( $O(r^2)$ ), joten tutkitaan verkon rakennetta ja huomataan siitä laskemista helpottavia rekursioita.

Piirretään verkko koodilla

def getBrwMat(r):
    bds = [(i,j) for i in range(-r, 1) for j in range(r+1) if j-i<=r]
    states = []
    for (i,j) in bds:
        states += [(i,j,k) for k in range(i, j+1)]
    matN = len(states)
    ret = [[0]*matN for _ in range(matN)]
    for ind1, (i,j,k) in enumerate(states):
        i2 = i-1 if k-1<i else i
        if j-i2<=r:
            s2 = (i2, j, k-1)
            ret[ind1][states.index(s2)] += 1
        j2 = j+1 if k+1>j else j
        if j2-i<=r:
            s2 = (i, j2, k+1)
            ret[ind1][states.index(s2)] += 1
    return Matrix(ret), states


myR = 3
a, states = getBrwMat(myR)
startState = states.index((0,0,0))
g = DiGraph(a)
vColors = {}
heights = {}
for r,color in enumerate(['blue', 'purple', 'pink', 'orange', 'red']):
    vsOfH = [j for j,s in enumerate(states) if s[1]-s[0]==r]
    vColors[color] = vsOfH
    heights[r] = vsOfH
g.show(figsize=8, vertex_size=70, vertex_labels=None,
       vertex_colors=vColors, iterations=700,
      layout='ranked', heights=heights)

Tutkittava verkko, kun $r=3$ . Lähtötila on alhaalla.

Tutkittava verkko, kun $r=3$ . Lähtötila on alhaalla.

Verkon rakenteellahan on selvä merkitys: ”putket” ovat sen hetkisiä (min, max) -rajoja. Putkessa voi ”seilata” edes takaisin ja päästä siirtyä aina seuraavaan putkeen eli siirtää (min, max) -rajoja yhdellä. Kunnes max-min = r, jolloin vaihteluväli ei saa enää kasvaa, ja jäädään seilaamaan siihen $r+1$ -putkeen, johon päädyttiin. Huom. ”putki” on itse asiassa verkkona (suuntaamaton) polku.

Keskitytään aluksi yhteen $m$ -putkeen ja sen kävelyihin eli yllä mainittuun ”seilaamiseen”. Merkitään $p(m, j, n)$ :lla ensimmäisestä päästä ( $j=0$ ) lähtevien kävelyjen määrää $m$ pituisessa putkessa, jotka päätyvät putken solmuun $j$ . Huomaa, että tämä on tietysti symmetrinen, jos putki käännetään.

Tutkitaan sitten koko verkkoa ja tehdään hieman yleisempi määritelmä: brwRec(m, r, n) = kävelyiden määrä verkolla, jossa ensimmäinen putki on $m$ -pituinen ja lähdetään tämän putken (jommasta kummasta) päästä. Huomaa, että haluttu lukumäärä on tällöin erikoistapaus brw(r, n) = brwRec(1, r, n). Nyt kävely voi seilata ensimmäisessä putkessa jonkun aikaa ja siirtyä sitten, päätyessään jompaan kumpaan päätyyn, seuraavalle tasolle, jolle jää loput askeleet (-1 sillä myös tasolle siirtymä kuluttaa yhden askeleen). Kävely voi myös olla siirtymättä seuraavalle tasolle ja päätyä lopuksi mihin tahansa putken somuun. Saadaan siis rekursioyhtälö

$\text{brwRec}(r, n) \\= \sum_{k=0}^{n-1} (p(m, 0, k) + p(m, m-1, k)) \text{brwRec}(m+1, r-1, n-1-k) \\+ \sum_{j=0}^m p(m, j, n)$

Koodataan nämä rekursiot ja käytetään memoisatiota:

pMemo = {}
def p(m, j, n):
    if n==0: return 1if j==0 else 0
    memoKey = (m,j,n)
    if memoKey in pMemo:
        return pMemo[memoKey]
    ret = 0
    if j>0: ret += p(m, j-1, n-1)
    if j<m-1: ret += p(m, j+1, n-1)
    pMemo[memoKey] = ret
    return ret

brwMemo = {}
def brwRec(m, r, n):
    if n==0: return 1
    memoKey = (m,r,n)
    if memoKey in brwMemo:
        return brwMemo[memoKey]
    ret = 0
    if r>0:
        for j in range(n):
            ret += (p(m, 0, j) + p(m, m-1, j)) * brwRec(m+1, r-1, n-1-j)
    ret += sum(p(m, j, n) for j in range(m)) #stay on the first path
    brwMemo[memoKey] = ret
    return ret

brw = lambda r,n: brwRec(1, r, n)

Tämä algoritmi toimii jo kohtuullisessa ajassa suurehkoillekin r:n arvoille, mutta tässä ongelmaksi tulee suuret n:n arvot. Matriisimenetelmässä taas suuret n:n arvot eivät tuottaneet ongelmaa, koska matriisin potenssiinkorotus on nopeaa ( $O(log n)$ ). Koitetaan saada molempien maailmojen paremmat puolet ja käytetään keskikokoisille $r$ :n ja $n$ :n arvoille rekursiivista metodia ja pienille $r$ :n arvoille matriisimenetelmää. Arvolle $r=2$ voidaan käyttää myös valmista kaavaa. Koska ratkaisujen lukuarvot kasvavat nopeasti hyvin suuriksi, annetaan testin vuoksi vastaukset jossain suuressa modulossa, esim 7347249713. Tehdään laskenta siis renkaassa $\mathbb{Z} / 7347249713\mathbb{Z}$ . Tämä onnistuu Sage-koodilla

MOD = 7347249713
MY_RING = Integers(MOD)
ZERO = MY_RING(0)
ONE = MY_RING(1)
TWO = MY_RING(2)

pMemo = {}
def p(m, j, n):
    if n==0: return ONE if j==0 else ZERO
    memoKey = (m,j,n)
    if memoKey in pMemo:
        return pMemo[memoKey]
    ret = ZERO
    if j>0: ret += p(m, j-1, n-1)
    if j<m-1: ret += p(m, j+1, n-1)
    pMemo[memoKey] = ret
    return ret

brwMemo = {}
def brwRec(m, r, n):
    if n==0: return ONE
    memoKey = (m,r,n)
    if memoKey in brwMemo:
        return brwMemo[memoKey]
    ret = ZERO
    if r>0:
        for j in range(n):
            ret += (p(m, 0, j) + p(m, m-1, j)) * brwRec(m+1, r-1, n-1-j)
    ret += sum(p(m, j, n) for j in range(m)) #stay on the first path
    brwMemo[memoKey] = ret
    return ret

def brwMatrixMethod(r, n):
    bds = [(i,j) for i in range(-r, 1) for j in range(r+1) if j-i<=r]
    states = []
    for (i,j) in bds:
        states += [(i,j,k) for k in range(i, j+1)]
    matN = len(states)
    startState = states.index((0,0,0))
    arr = [[0]*matN for _ in range(matN)]
    for ind1, (i,j,k) in enumerate(states):
        i2 = i-1 if k-1<i else i
        if j-i2<=r:
            s2 = (i2, j, k-1)
            arr[ind1][states.index(s2)] += 1
        j2 = j+1 if k+1>j else j
        if j2-i<=r:
            s2 = (i, j2, k+1)
            arr[ind1][states.index(s2)] += 1
    return sum( (Matrix(MY_RING, arr)**n)[startState])


def brw(r, n):
    if r==2: return (TWO**((n+3)/2) if n%2 else 3*TWO**(n/2)) - 2
    #if r==3: TODO
    if r<5: return brwMatrixMethod(r, n)
    return brwRec(1, r, n)

testCases = [
    (1, 3),
    (2, 3),
    (2, 4),
    (2, 7),
    (3, 5),
    (3, 12),
    (4, 15),
    (5, 16),
    (2, 20),
    (2, 109),
    (3, 67),
    (3, 174),
    (11, 23),
    (12, 43),
    (15, 38),
    (21, 45),
    (3, 923),
    (2, 242398),
    (4, 923847),
    (2, 12345678910),
]

for (r, n) in testCases:
    print (r, n, brw(r,n))

Tämä tulostaa testitapauksiin vastaukset

"""
1 3 2
2 3 6
2 4 10
2 7 30
3 5 26
3 12 1028
4 15 12328
5 16 37144
2 20 3070
2 109 868658622
3 67 1447052816
3 174 1595158214
11 23 8099490
12 43 1135120293
15 38 2365381958
21 45 5576619597
3 923 1918313630
2 242398 3965505254
4 923847 2261069692
2 12345678910 4322229703
"""

Tehtäviä

Yllä, arvolle $r=2$ (kakkosen potensseja sisältävä) kaava toimi myös kun laskenta tehtiin modulossa. Pohdi kuinka se onnistuisi muille kaavoille (ts. jos potenssiin korotettava luku ei ole kokonaisluku).
Perustele tapauksen $r=2$ kaava $2^{\lfloor \frac{n+2}{2} \rfloor} + 2^{\lfloor \frac{n+1}{2} \rfloor} - 2$ .
Vinkki: Käytä edellä ollutta rekursiota. Huomaa, että 3-putken kävelyt ovat muotoa askel keskelle ja takaisin jompaan kumpaan päätyyn jne. Käytä geometrisen summan kaavaa.
Laske verkon spektri. On tunnettua, että polun spektri on $\{2\cos \frac{j\pi}{n+1}, j=1,\dots, n\}$ . Onko näillä joku yhteys (onhan verkon osat polkuja)?
Keksi kaava yleiselle $r$ .
Tutki vaihteluvälin jakaumaa $n$ :n pituisille ${-1, +1}$ -kävelyille. Kuinka sen keskiarvo käyttäytyy?

Sanan osajonojen lukumäärä

Sanan osajono on sanasta järjestyksessä poimituista merkeistä muodostettu uusi sana. Esimerkiksi sanalla ’ABBA’ on osajonot (11 kpl) [”, ’A’, ’AA’, ’AB’, ’ABA’, ’ABB’, ’ABBA’, ’B’, ’BA’, ’BB’, ’BBA’]. Huomaa, että esim. ’AB’ saadaan kahdella eri tavalla, mutta lasketaan vain kerran.

Kuinka laskea sanan $w$ osajonojen lukumäärä $f(w)$ ? Oletetaan että $w$ :ssä esiintyvien merkkien lukumäärä on jokin ennalta kiinnitetty pienehkö luku. Halutaan löytää lineaariaikainen algoritmi.

Ratkaisu

Olkoon $w$ sana, jonka osajonot halutaan laskea. Olkoon $w$ :n pituus $n$ ja siinä esiintyvät merkit $\{c_1, c_2, \dots, c_m\}$ . Kun samaistetaan jokainen $w$ :n osajono $v$ siihen indeksien osajonoon, josta $v$ löytyy kaikkein aikaisimmin, vastaa jokainen sanan osajono solmusta $-1$ lähtevää kävelyä suunnatulla graafilla $G$ , jonka solmut ovat $V(G) = \{-1, 0,1,\dots, n-1 \}$ ja solmusta $j$ lähtee kaaret $(j, k_c)$ missä $k_c$ on indeksi, jossa merkki $c$ esiintyy $w$ :ssä seuraavan kerran lukien indeksistä $j+1$ lähtien (ja jos $c$ ei enää esiinny, niin kaarta ei ole) ja tämä vastaavuus on bijektiivinen.

Esimerkiksi sanalle w = ”ABBA” graafi näyttää seuraavalta

Sanan ”ABBA” graafi. theFolls on dicti, jonka ao Python-koodi antaa, koodissa indeksit ovat siirretty yhdellä ja itse asiassa theFolls antaa seuraavan etsimisen alkaa itse siitä indeksistä missä ollaan.

Olkoon $A$ edellä kuvatun suunnatun graafin vierekkäisyysmatriisi. Osajonojen lukumäärä on tällöin ensimmäisen rivin summa seuraavasta matriisista

$\begin{aligned} \sum_{j=0}^\infty A^j = \left(I - A \right)^{-1} \end{aligned}$

Huomaa, että edellä summa on itse asiassa vain äärellinen, sillä kaaret vievät aina eteen päin, joten graafin kävely voi olla korkeintaan $n$ :n pituinen ja $A^{n+1}=0$ , joten summan suppenemisen eikä $I-A$ :n kääntyvyyden kanssa tule ongelmia.

Merkitään $B=I-A$ . Mehän tarvitsemme vain sen ensimmäisen rivin summan, joten koko matriisia $A$ ei tarvitse muodostaa, eikä summaa $\sum_{j=0}^{n}A^j$ laskea, vaan riittää ratkaista

$B^{T}x = e_1$ ,

missä $e_1$ on sarakevektori, jonka ensimmäinen komponentti on 1 ja loput nollia. Kun tämä kirjoitetaan Gauss-Jordanin menetelmän vaatimaan muotoon (jossa ei ratkaista koko käänteismatriisia, vain sen ensimmäinen rivi), niin huomataan, että ratkaiseminen on erittäin helppoa, sillä $B$ on alakolmiomatriisi, jossa on ykkösiä diagonaalilla ja jokaisessa sarakkeessa korkeintaan $m$ :ssä paikassa $-1$ (ja muuten nollia). Näin ollen saadaan seuraava lineaariaikainen algoritmi (Python):

def theFolls(w):
    """map c -> list a where a[j] is the next index of c in w
    when start searching from the index j"""
    alph = list(set(c for c in w))
    n = len(w)
    ret = {c: [n]*n for c in alph}
    prevs = {c: -1 for c in alph}
    for i, c in enumerate(w):
        for j in range(prevs[c]+1, i+1):
            ret[c][j] = i
        prevs[c] = i
    return ret

def countSubsequences(w):
    """How many words can be obtained from w as subsequences"""
    n = len(w)
    folls = theFolls(w)
    edges = [
        [folls[c][i]+1 for c in folls if folls[c][i]<n]
        for i in range(n)
    ]
    ret = 0
    a = [0 for _ in range(n+1)]
    a[0] = 1
    for i in range(n):
        for j in edges[i]:
            a[j] += a[i]
    return sum(a)

Junarata-tehtävä 2

Tiedosto fenwickTree.py

class FenwickTree:

    def __init__(self, n):
        self.n = n+1
        self.arr = [0]*self.n

    def incr(self, ind, val):
        """Increment the value at index ind by val"""
        if ind<0: raise Exception("Index out of bounds")
        curr = ind+1
        while curr<self.n:
            self.arr[curr] += val
            curr += curr&(-curr)

    def getSum(self, upto):
        """Get the prefix sum up to index upto"""
        s = 0
        ind = upto+1
        while ind>0:
            s += self.arr[ind]
            ind -= ind&(-ind)
        return s

Tiedosto junarata6.py

from fenwickTree import FenwickTree

class IntvalSet:
    """Set where the included natural numbers
        are given as intervals [a, b]"""
    def __init__(self):
        self.intvals = []

    def has(self, x):
        for (a,b) in self.intvals:
            if a<=x and x<=b: return True
        return False

    def add(self, a, b):
        self.intvals.append((a, b))

    def getAll(self, lowerBdd=0):
        """generates all numbers in the set.
            If lowerBdd is given,
            generates only numbers >= lowerBdd"""
        for (a,b) in self.intvals:
            for x in range(max(a, lowerBdd), b+1):
                yield x

    def __str__(self):
        return "[]" if len(self.intvals)==0 else " U ".join(["[%d, %d]" %x for x in self.intvals])


def getData(num):
    ret = []
    with open("junar"+str(num)+".in") as f:
        ret = f.read().split("\n")
    return [int(x) for x in ret[1:] if len(x)>0]

def getPoss(p):
    """Get list of positions of elements in the
        permutation p of [1..n]
        e.g. getPoss([2,5,3,1,4]) = [None, 3, 0, 2, 4, 1]
    """
    ret = [None]*(len(p)+1)
    for i,x in enumerate(p):
        ret[x] = i
    return ret

"""
Let's denote types of orderings of k-1, k and k+1
(or k and k+1; k-1 and k in cases k=1; k=n  like this:

123: 0
132: 1
213: 2
231: 3
312: 4
321: 5
_12: 6
_21: 7
12_: 8
21_: 9
"""
def getTyyppi(a, b, c):
    if a<b:
        if b<c: return 0
        if a<c: return 1
        return 3
    else: #b<a
        if b>c: return 5
        if a<c: return 2
        return 4


def getNumOfRounds(poss):
    """How many rounds does it take to gather the numbers.
        @param poss: index-positions of the permutation
    """
    curr = 1
    prevInd = poss[1]
    r = 1
    for curr in range(1, len(poss)):
        currInd = poss[curr]
        if currInd<prevInd: r += 1
        prevInd = currInd
    return r


def countScore2Pairs(arr):
    n = len(arr)
    poss = getPoss(arr)
    toConsider = []
    t = None
    tempSet = None
    for i, x in enumerate(arr):
        if x==1: t = 6 if poss[x+1] > i else 7
        elif x==n: t = 8 if poss[x-1] < i else 9
        else: t = getTyyppi(poss[x-1], i, poss[x+1])

        tempSet = IntvalSet()
        
        if t==2:
            tempSet.add(poss[x-1], poss[x+1]-1)
        elif t==5:
            tempSet.add(0, poss[x+1])
            tempSet.add(poss[x-1], n-1)
        elif t==1:
            tempSet.add(poss[x-1]+1, poss[x+1])
        elif t==7:
            tempSet.add(0, poss[x+1])
        elif t==9:
            tempSet.add(poss[x-1], n-1)
        
        toConsider.append(tempSet)


    """
    mat = [[1 if toConsider[i].has(j) else 0 for j in range(n)]
           for i in range(n)]

    print ("To Consider Matrix:")
    #print (mat)
    for row in mat: print row
    """

    tree = FenwickTree(n)

    events = [[] for _ in range(n+1)] #grouped by time
    #print ("Events len is "+str(len(events)))
    for i, c in enumerate(toConsider):
        for (j1, j2) in c.intvals:
            #print ("Pushing events for interval [%d, %d]" %(j1, j2))
            events[j1].append([j1, i, 1])
            events[j2+1].append([j2+1, i, -1])
    
    numPairs = 0
    for k, evs in enumerate(events[:-1]):
        for (j, i, d) in evs:
            tree.incr(i, d)
        for (j1, j2) in toConsider[k].intvals:
            numPairs += tree.getSum(j2) - tree.getSum(j1-1)
        #subtract the adjacent ones
        x = arr[k]
        if (x>1 and toConsider[k].has(poss[x-1])
            and toConsider[poss[x-1]].has(k) ): numPairs -= 1
        if (x<n and toConsider[k].has(poss[x+1])
            and toConsider[poss[x+1]].has(k) ): numPairs -= 1
    
    return numPairs/2


def junarata(arr):
    n = len(arr)
    poss = getPoss(arr)
    origRounds = getNumOfRounds(poss)
    numPairs = countScore2Pairs(arr)
    return "%d %d %d" %(n, origRounds-2, numPairs)



import random

#a = range(1, 10); random.shuffle(a); print(junarata(a));

#Solve putkaposti
for i in range(1, 10):
    print(junarata(getData(i)))



#Test execution time
"""
import time
startT = time.time()
print (junarata(getData(9)))
print ("Time: "+str(time.time()-startT))
"""