Counting Subtrees

Arrangements, sorting, packing, partitions, critical path analysis, networks, graphs, ...
Post Reply
szymczak
Posts: 15
Joined: Wed Nov 23, 2011 8:43 am

Counting Subtrees

Post by szymczak » Tue Sep 29, 2015 6:23 pm

Consider the following rooted binary tree.
Image
For each node v, consider the subtree rooted at v. For the tree above we get the following subtrees (ignoring duplicates).
Image
Note how the second and third subtree above are considered distinct.
We notice that there are 5 distinct subtrees.
Let T(k) be the number of rooted binary trees with exactly k distinct subtrees.
For example, T(1) = 1, T(2) = 3, T(6) = 14487.

Find T(16) mod 109.
Image

User avatar
jaap
Posts: 528
Joined: Tue Mar 25, 2008 3:57 pm
Contact:

Re: Counting Subtrees

Post by jaap » Wed Sep 30, 2015 6:09 am

Unfortunately that sequence (including the 16th term) is a bit too easy to find in the OEIS, which in turn links to the extensive stackexchange discussion that generated the sequence.

szymczak
Posts: 15
Joined: Wed Nov 23, 2011 8:43 am

Re: Counting Subtrees

Post by szymczak » Wed Sep 30, 2015 2:58 pm

Yea I wrote the long SE post. I was hoping someone here could come up with a more elegant solution, a way to find say T(100). Should have added a twist so it's not on OEIS.
Image

User avatar
mpiotte
Administrator
Posts: 1914
Joined: Tue May 08, 2012 4:40 pm
Location: Montréal, Canada

Re: Counting Subtrees

Post by mpiotte » Sun Oct 04, 2015 8:44 pm

Interesting question. We could have made a PE problem out of this.
Dynamic programming is once again effective for this kind of problem.
I get the following result with O(n6) time and O(n5) memory complexity, keeping things under 1 minute:

Code: Select all

T(  1) =                                                                                        1
T(  2) =                                                                                        3
T(  3) =                                                                                       15
T(  4) =                                                                                      111
T(  5) =                                                                                     1119
T(  6) =                                                                                    14487
T(  7) =                                                                                   230943
T(  8) =                                                                                  4395855
T(  9) =                                                                                 97608831
T( 10) =                                                                               2482988079
T( 11) =                                                                              71321533887
T( 12) =                                                                            2286179073663
T( 13) =                                                                           80984105660415
T( 14) =                                                                         3144251526824991
T( 15) =                                                                       132867034410319359
T( 16) =                                                                      6073991827274809407
T( 17) =                                                                    298815244349875677183
T( 18) =                                                                  15746949613850439270975
T( 19) =                                                                 885279424331353488224511
T( 20) =                                                               52902213099156326247243519
T( 21) =                                                             3349349968648038146756090367
T( 22) =                                                           224008474237110735818328494463
T( 23) =                                                         15784481675552990217004355155455
T( 24) =                                                       1168980983609271828189092132883711
T( 25) =                                                      90789628326447871556044333389045759
T( 26) =                                                    7379671591913867687413635289032264447
T( 27) =                                                  626614739073537117360355723435583202303
T( 28) =                                                55485577115128371124823191707413779511295
T( 29) =                                              5115454359795311832797712178681313437868031
T( 30) =                                            490310519553284322702752182845108549592875519
T( 31) =                                          48791468725761851199449519655205582018214453247
T( 32) =                                        5034369951561638858570924846904585834050350695423
T( 33) =                                      537964973713998388999889514657023448872676106387455
T( 34) =                                    59467875273261002743066517406088534107058019534881791
T( 35) =                                  6793164751908330197524197123822029953827509352828620799
T( 36) =                                801109989079284363555583411034360284129353799268061933567
T( 37) =                              97439774800234157527936865903950535420883900577338857119743
T( 38) =                           12212966842762677445920937938343431481488602466294029830690815
T( 39) =                         1576102568966667715777124386738220705637159432149738884321779711
T( 40) =                       209258483862486797315415155447759265226105288742426622806920744959
T( 41) =                     28562201839019312638597793617092469246748808899696158195441127292927
T( 42) =                   4004990457079402844735981023665557864353516616675783648506722317185023
T( 43) =                 576526026299295223790052290449858497256713186201245773065940447813943295
T( 44) =               85146064699321983042897937056238854596395499358114811035260281811412877311
T( 45) =            12893542166293539662215995925168478556182067594359044316243370721375672729599
T( 46) =          2000727533814254009821605140674456289993982732539619481011084097733785705259007
T( 47) =        317957574456214497639799228877982309585978351841711260500637449128815475171262463
T( 48) =      51723046966659730222763456537636334604630427786738469658935909161758736064153665535
T( 49) =    8608171476830272116837045761769711695845434641798904200741065837473162376652967903231
T( 50) = 1464998441941572191008123746853973707755444234070174247285019683005025893928381236264959
Total time =    56.47 s.
Alternatively, I can also do O(n7) time with only O(n4) memory.
Anyone else willing to confirm?
Image

szymczak
Posts: 15
Joined: Wed Nov 23, 2011 8:43 am

Re: Counting Subtrees

Post by szymczak » Mon Oct 05, 2015 5:51 pm

Confirmed up to T(20). Those are all the values I had though. I was thinking about submitting it as a potential Euler problem but I didn't have a good solution for it. Need to figure out your DP.
Image

User avatar
mpiotte
Administrator
Posts: 1914
Joined: Tue May 08, 2012 4:40 pm
Location: Montréal, Canada

Re: Counting Subtrees

Post by mpiotte » Tue Oct 06, 2015 12:16 am

Results for 51 to 100 in 2h20:

Code: Select all

T(51)=254835085445543722215634944512650148042285342540187050263916802228774246942593912081809407
T(52)=45287850941429082230865074942447214927480513849149509424582203044956495749331824540996009983
T(53)=8218963331111956374791801588528351546948056492993708307325513881997526207328619264182647914495
T(54)=1522596907111625798328686575910961990495813331669959225982902238125888016145066514677839522004991
T(55)=287814090727233343624724823200646668827772592350698110858734950111852996110935883732866646513418239
T(56)=55492147073182888752732748303425634686236489102328973223305697845440886854038884262659866188537724927
T(57)=10908943524318124866584558646577574495553717331120678665106248223683710369781228810542474617915602632703
T(58)=2185801382287953120354455770936078446546062996422979632701095128528097200693145017748413064359853671776255
T(59)=446237646438675493197277447367311468775663138709802489429277567738843688655751379079922408800073713837146111
T(60)=92790843834597575121485103242705040622688846873730368075858749213891410572441916279604111799486761449540812799
T(61)=19646641171352959128943032240587277469876366612142648815606393436583923192190783828436826394378226808238100185087
T(62)=4234297317141392818776911727668104602518488571957721872591943900782418973864516124091250043305321553783485588701183
T(63)=928654755551200768501604583930786589349974337764115289022292847199094382526705850766775504030165111160736477280731135
T(64)=207196035108109817419415207529410552319657099478997256005811238055683316176562878461295520559720967110173930205393190911
T(65)=47015512965942941486348873225844767761817421240038369180684005370273252393794291683498718907840707989451350827899436400639
T(66)=10847152207421585934354391653976470993032677409411645474042644823439928801374144082861597910459472615120528085963582445191167
T(67)=2543847202437147943123431892103657235324535685483601527456696700482985141653880991051275786510918494677879660473516395487821823
T(68)=606256586544372096615517709260195923184537323966906378521746448397643397418750686429305600939763256111873776724831013276956491775
T(69)=146792836641992755140949619041659265194932531118677411078798580252172017123178764277827372524214237407641120120899832528614543327231
T(70)=36102077128783856959458374601567918829170934800282306711028605033992963158443958589798870173789422328512656453175904351672706046361599
T(71)=9016489971634686129975111308337566526313370570256573117977358040749860454759660046995889040008808845073103840385576376269929730558394367
T(72)=2286245119110048953986679359440690960661397501936954984932288490549292397629633943067234258744834776961083676101294254298720870046387666943
T(73)=588427341086117363145869973606117129388791944015177812079253508311440080714882269523329960315250997919393787263379619837499401477929289383935
T(74)=153693373181961024156183699261159229332498828612139388088088745189580538270567166315550720910963090823791657644408653610031518128020391216545791
T(75)=40730500267690831653594440708975666112777665312814993294015508096272887377484130284652377559589198909310096450695441584575421619756130760424161279
T(76)=10949623299336364303604974431085058392596398853862021280565818257282910272468689554353647415429386940744044923521631732561680283705330058130660261887
T(77)=2985439981349112881966444020503233551718324073512607849410831795577100961120063015748039172099907850914318560285188112314434512081124401389970615435263
T(78)=825399915024681556607625427652174413646246831616488253415790565473334494396064627837712297536563191673019723789717703780677649014794676109032959399428095
T(79)=231359196682252971620348967126049803799281177361644710778152431156321965453021650877647845776379878550157036867866741068936823820595923110946597997182451711
T(80)=65735028246138729486111958146802926444641136913057796673618618301813606529291936341059827649359277626195785453717614323378616900798604520795521409802085335039
T(81)=18928585688860145587610148324230761367391260801847422051653321664351830341359231985244452433198930143540282532691478895599229419250437593801642958544640722075647
T(82)=5523013871241715479945855035168627612834284381097964153704930535910572020227555358884568744174883417447026810140729028754874120900524941556411545165911664984326143
T(83)=1632668684755395169078949317322595412246175362026555791784319666400016733947245545905211998058819628605027648176353368659376734697688275978755080183362128727062347775
T(84)=488892143468568249491626529094674672405558366621449853452020286719845023387976786398253014293042170537641627495881457860280747792395466015782797952863675532002524659711
T(85)=148269547620695602327411393167658345704244757721363029851313252963024137475342178368067201246719896412952991907871610438316386838424952207095851398724505631337231134228479
T(86)=45535204619014192923397238751714292535848318131816373557799673250766054116479046447118139105405058046738530346176804327067789322714123896012402837154344368845573041608982527
T(87)=14159019469403608334295290534756351645535814250289612571014265391221839785647232197098984709773899437836762934443180419011521085400487770844585438747551810923772712883043958783
T(88)=4457027461420224550250224177702627598591419828202762445302276523564238570838863418412685981558313220009218939930992608827885311170032638336495432873243051685056072772191535497215
T(89)=1420106373451041197128706760959469423251147066674187382035973037270663937609069489102489346907294161752941819692424327112735873734795689846893046727961083091947353931366386821496831
T(90)=457929393554747857904738028107197246892301742012245327563897290337381949602749294549169912891326457606789580170867916210270702659044552334424647076098693610541935002307733299543408639
T(91)=149423287370483314511628840599093582643271044102235204655747515666474963152324645262910252800644723735598763396498163386314417300879446237150635318720554288176213041717249143768380080127
T(92)=49331187797795147340599872468700165856169475256580009818021273647256623229797357374780555374648544220704635823268646924491220546089555451971400516369455557916642651001718178457215674875903
T(93)=16475973528303040263834696975864047068631891522197163561555432124848322365283985033522566336809649894684667307384426494090663648109002485355048606922876782179092654796792580176661489622450175
T(94)=5566097368536569872461666010363123942115725963564114547045811936105465550753961104488123945346276199649713699812608302950604946037279419282243078297770075960457029858541662151137527449247547391
T(95)=1901804680891364522653195739383310106941797935709273838428473871954079573372741841127578543145019408602982471511376419070301745554594184574503992306504891408045378417435901893569808837260104171519
T(96)=657117158009218553023967905328870809671140757015463607947560260156787681924584431099344818995651643370308471034988220459436791777620450369228415634985667664992297136076947064192295256376246269902847
T(97)=229577337946512462203729298756201222174834619525408637608026400807460449720958705002750282656270250271861262373008163295259337771557693281807273580955420060018984634441004592388056028408828789600026623
T(98)=81091097716583798198367504419675193302839237890226926251564079120408267913597110561334222993977241436483266420378997047202752749995965263392456792524940331094197041643412734736794702961010010484570587135
T(99)=28955096488971376189813894570849386252788973564824795476013139629315218166255385930937462771759852727500957723974812026319151930550590150964866098472781333195419544617904119392156067643028895288111778496511
T(100)=10450457132583269349450617394714442608751206265530678410381635192788309450706495588686041179258450375873075914837789821383155323561722331425003064223479600883976025606894381013869321432750565338950806462791679
Total time =  8371.34 s.
Dynamic Programming approach part I:
Expand
We shall count all possible trees with n distinct sub-trees by building sets of sub-trees bottom-up. A set of sub-trees is counted in T(n) if it contains exactly n distinct sub-trees, and each sub-tree except one (the root) is part of one or more larger sub-trees in the set. Sets of sub-trees are characterized by the following numbers:
- h : the maximum height of any sub-tree in the set.
- c : the number of distinct sub-trees part of one or more larger sub-trees (connected sub-trees).
- u : the number of distinct sub-trees with height<h not part of any larger sub-tree (unconnected sub-trees).
- t : the number of distinct sub-trees with height=h (top sub-trees). Since they have maximum height, none can be part of a larger sub-tree.

Let Ch(c, u, t) be the number of sets of sub-trees with characteristics (h, c, u, t).
Then $T(n) = \sum_h C_h(n - 1, 0, 1)$.

We also have:
$\displaystyle C_{h+1}(c,u,t) = \sum_{{i,j,k} \atop {i + j + k = c + u}} C_h(i, j, k) B(i,j,k,u,t)$
This formula essentially says that each set with characteristics (h, i, j, k), after adding t sub-trees of height h+1, will form a set with characteristics (h+1,c,u,t) in B(i,j,k,u,t) different ways. The constraint i+j+k=c+u comes naturally from the fact that each sub-tree at height &le; h becomes either connected (c) or unconnected (u) sub-trees at height=h+1. B(i,j,k,u,t)=0 for constructions that are impossible.

Note that h,c,u,t all have size O(n). The triplet (i,j,k) has size O(n2) due to constraint i+j+k=c+u. Thus computing T(n), knowing B, takes time O(n6).

I shall elaborate on the evaluation of B in a future post.
Image

User avatar
mpiotte
Administrator
Posts: 1914
Joined: Tue May 08, 2012 4:40 pm
Location: Montréal, Canada

Re: Counting Subtrees

Post by mpiotte » Tue Oct 06, 2015 11:14 pm

Dynamic Programming approach part II:
Expand
Given a set including the empty sub-tree, of i connected sub-trees, j unconnected sub-trees and k sub-trees of top height h, we want to build t sub-trees of height h+1. Each new sub-tree is composed of an ordered pair of smaller sub-trees, at least one of them of height h. The number of such distinct pairs of sub-trees is (i+j+k)2-(i+j)2, i.e. the total number of pairs minus the pairs having no sub-tree at height h. There are ${(i+j+k)^2-(i+j)^2} \choose t$ ways to build our t new distinct sub-trees.

Let A(i,j,k,t) be the number of ways to build t new sub-trees of height h+1 using all of the j unconnected and k top sub-trees. A(i,j,k,t) can be computed as the number of all possible sub-tree constructions minus those that do not use all of the unconnected and top sub-trees.
$\displaystyle A(i,j,k,t) = {{(i+j+k)^2-(i+j)^2} \choose t} - \sum_{{p,q} \atop {(p,q) \ne (j,k)}} {j \choose p}{k \choose q}A(i, p, q,t)$

We can now use A to evaluate B. Let's say we build our new t sub-trees by picking exactly p sub-trees from the j unconnected and the q from the k top sub-trees. This creates a set containing t top sub-trees, u=j+k-p-q unconnected sub-trees and c=i+p+q connected sub-trees.
$\displaystyle B(i,j,k,u,t) = \sum_{{p,q} \atop {p+q+u=j+k}} {j \choose p}{k \choose q}A(i,p,q,t)$

A is memorized and takes O(n4) memory .
If B is memorized, it takes O(n5) memory and evaluating T(n) takes O(n6) time.
If B is not memorized, memory requirement is lowered to O(n4), but the time becomes O(n7).

To reduce memory usage, A and B are computed and stored modulo a 32-bit prime number, and thus each value requires only 4 bytes. T(n) is also computed modulo the same prime. The computation is then repeated for sufficient different prime moduli and the exact result is recovered from the Chinese Remainder Theorem.
I'm sure there are ways to simplify this approach, but I did not investigate further. I'll post some code once I have cleaned it up a little.
Image

szymczak
Posts: 15
Joined: Wed Nov 23, 2011 8:43 am

Re: Counting Subtrees

Post by szymczak » Tue Oct 06, 2015 11:52 pm

Brilliant solution! I think I might be able to simplify B but I need to think on it more. Unfortunately it'll have to wait until tomorrow.. I've barely had enough time to read your posts these last two days :(
Image

szymczak
Posts: 15
Joined: Wed Nov 23, 2011 8:43 am

Re: Counting Subtrees

Post by szymczak » Thu Oct 08, 2015 6:30 pm

I wrote my own (slow) code and confirmed your values up to T(50). It uses your algorithm so how well it "confirms" is debatable. For the sister sequence T' where we don't distinguish between
Image
the only thing that needs to change is the number of new subtrees k(i+j)+k(k+1)/2 in the calculation of A. The code is in Python2.7 and it uses numpy. Calculating T(50) takes about 10 minutes.

Code: Select all

T'(1)=1
T'(2)=2
T'(3)=6
T'(4)=25
T'(5)=137
T'(6)=945
T'(7)=7927
T'(8)=78731
T'(9)=906705
T'(10)=11908357
T'(11)=175978520
T'(12)=2893866042
T'(13)=52467157456
T'(14)=1040596612520
T'(15)=22425725219277
T'(16)=522102436965475
T'(17)=13064892459014192
T'(18)=349829488635512316
T'(19)=9983594851603361440
T'(20)=302595243049622291132
T'(21)=9709654327333904247720
T'(22)=328903501349364758535742
T'(23)=11730834245562877653553046
T'(24)=439499477724515631342816210
T'(25)=17259105843445284845953724977
T'(26)=709001872452143272423199623310
T'(27)=30412531388984531527650718225872
T'(28)=1359882311145651991023001603605331
T'(29)=63287078411450376102052573665416495
T'(30)=3061000444512448121107521821538907542
T'(31)=153659477474364121511372973556664569000
T'(32)=7995672159062043048918465357429752114537
T'(33)=430762266692393567329188251346127405562256
T'(34)=24000860250584249352744670773943113490702210
T'(35)=1381564098165671982201523273255799997524948462
T'(36)=82081551157095967346141807570437481606862885930
T'(37)=5028632681971636699786149067669638356244385933912
T'(38)=317399328521955656168741128730106175577741456307716
T'(39)=20623209084695165932602694167147745052335058810134905
T'(40)=1378357817942999597505867314474520294477794992100237467
T'(41)=94689537403423524528866753957217500202434426019729924875
T'(42)=6681458463034996867636810855811215755925054063022843457462
T'(43)=483926889203663723347220556626654810090447112583817037678646
T'(44)=35954338094829256894524920526678058655784779349619905842274125
T'(45)=2738562812650855914400529909833599741271576709207369964088421487
T'(46)=213718857122989323339364170908419911918841623238860998005357418665
T'(47)=17079369249041898531359197994299681178883933073815373088069470515392
T'(48)=1396948203185223071392625392717685335255334095391664072757423526911131
T'(49)=116882303493124092348455672322477033799790643961967597785062203772221860
T'(50)=9999238552244796357461982943895456551979492719134019786696565696441397152
T'(51)=874245083489161484924827693815494801763032648833194756817087540598837277555
T'(52)=78082571579635930573669814620079895085447362885948654627295419729906106936582
T'(53)=7121047151134377177271918046725429059209297757174655919396212949551250522643656
T'(54)=662863162428457144080702189230716513945533188450535494739460713421524917906238225
T'(55)=62953867306295553469403153262600744872892877082183242322537987285573803277838780172
T'(56)=6097813253866341850267643517792216477211951028774395471171801920580989201460832793335
T'(57)=602170879766887225098497786760958485739987872241127276273402244011625891600226741146657
T'(58)=60604630083061404980235881407169169630862930869653235565773812719485573416414875308629191
T'(59)=6214185799485070374623688284668871769488914746648225437661589893135077776042641048339896045
T'(60)=648951368449440444086127770295302998655584699961021393031326152093639801258957097962989693552
T'(61)=69000394787029810249613325306581568357736077965705583853528531527540370037411244297689599570030
T'(62)=7467395844421242060029940871627887215912045441528843683702244511310868417341296936501534359152382
T'(63)=822310891297469761405329050745738161611261495782082830687788811997496827873119118469126508534507696
T'(64)=92114391217755012834305284336117081650072713636231943752916793614218708970949544927017381695434424865
T'(65)=10493559214257242256526457352556131378333286574545327589425798910700798275838126647311299643028471635697
T'(66)=1215360640648237760639718339151422176513325564738518244436562071102838557397903931301249284466119513861045
T'(67)=143074194363064208644657327225914416800855433373973027852043072470374249429872869715070169583738163954344137
T'(68)=17115180726876180450731922408823710770636530849589657482602321235329081665528720187379096876726373367017386781
T'(69)=2079980551609732309595154934543650430068950077275253213511414410894279837525855761335787978914460331864889113575
T'(70)=256738736796213362308527599608757065176561083782237398354650655974608096220515215875658163796782791914144917865737
T'(71)=32179377277830043835868719196550917494190748987905081552396926458215872684015547511739274087781392766356937542473135
T'(72)=4094686441906531182091415133948943821993414121670322212671922946341533297714845796664465236171508742284040962606519652
T'(73)=528841943194356796809513822339846662472514992024371643639311767919782445952963783453853430358680924615671760419065882781
T'(74)=69310889092971535131136521334524312335340050079207074130818235116138588601709144037724993477120401490590970020109919919765
T'(75)=9216330550066055115593707202968113470938861144433599860381353450786857511821018631580608603714687022852864894832438572716007
T'(76)=1243107587282513240015523325992473472219888428789811444909910259254120332111907597962522383754574765924431842603199251605229965
T'(77)=170046989689877952007767633162728817556301178528193742909205121127242154105496583912331107219510839524692751033549089473418645127
T'(78)=23586116798766560891370993889548421899731997481568867835181572042262811623537310162340970452243924975550807349164743727933966245891
T'(79)=3316584624490571969875069497921439883891376783130125456415144307074796058025586538441924137680152938388360999226943603003205727617985
T'(80)=472710033613557597987526205304465938761095903914050658546993862065860713113526131183830707315189953367887597235923760016100997254848037
Here is the code subtree.py

Code: Select all

from subtree_helper import *


def A(i,j,k,t):
    newt = (i+j+k)**2 - (i+j)**2 if plane else k*(i+j)+k*(k+1)/2
    return  ( binom(newt, t) - sum( binom(j,p) * binom(k,q) * A(i,p,q,t)
                                    for p in xrange(min(j,t)+1)
                                    for q in xrange(min(k,2*t)+1)
                                    if (p,q) != (j,k) ) )


def B(i,j,k,u,t):
    return sum( binom(j,p) * binom(k,q) * A(i,p,q,t)
                for p in xrange(min(j,t)+1)
                for q in [j+k-p-u] )


def C(h,c,u,t):
    if (h,c,u,t) == (1,1,0,1): return 1
    if h <= 1: return 0
    return sum( C(h-1,i,j,k) * B(i,j,k,u,t)
                for i in xrange(1,c+1)
                for j in xrange(c+u-i)
                for k in [c+u-i-j] )


def T(n):
    return sum( C(h,n,0,1) for h in xrange(1,n+1) )




if __name__ == "__main__":

    n = 20
    plane = True                                         # CHILDREN ORDERED?
    toc = time()

    binom = Memo_Array([4*n*n+1,n+1]) (binom)            # DECORATE FUNCTIONS
    B = Memo_Array ([min(n,50)+1]*5, partial=n>50) (B)
    A = Memo_Array ([n+1]*4) (A)
    C = Memo_Array ([n+1]*4) (C)
    T = Memo_Array ([n+1]*1) (T)

    for i in xrange(1,n+1):
        print 'T({}) = {}'.format(i,T(i))

    h,m,s = hms(time()-toc)
    print 'Time elapsed {}.{}.{}'.format(h,m,s)
Here is the helper code subtree_helper.py

Code: Select all

import operator as op
import numpy as np
from functools import wraps
from math import factorial as fact
from time import time


def prod(L): return reduce(op.mul,L,1)
def perm(n,k): return prod(xrange(n,n-k,-1))
def binom(n,k): return 0 if k<0 or k>n else perm(n,k)/fact(k)



def hms(t):
    h = int(t/24/60)
    m = int((t-(h*24*60))/60)
    s = round(t-(h*24*60+m*60),3)
    return h,m,s



def Memo_Array(dims, fndims=[], partial=False, nptype='object', negv=lambda x:0):
    """Memoize into a numpy array
            - partial for when the cache isn't large enough to store all the values
            - nptype is the datatype for the numpy array
            - negv provide a base case for negative arguments
    """
    def factory(func):
        cache = func.cache = np.zeros(dims+fndims,dtype=type)
        hit = func.hit = np.zeros(dims,dtype='bool_')
        @wraps(func)
        def wrapper(*args):
            if min(args) < 0: return negv(args)
            if partial and reduce(lambda x,i: x or args[i]>=dims[i], xrange(len(args)), False):
                return func(*args)
            if not hit[args]:
                cache[args] = func(*args)
                hit[args] = True
            return cache[args]
        return wrapper
    return factory
Image

User avatar
mpiotte
Administrator
Posts: 1914
Joined: Tue May 08, 2012 4:40 pm
Location: Montréal, Canada

Re: Counting Subtrees

Post by mpiotte » Sat Oct 10, 2015 5:20 pm

Here's my own code in Java.
Things to be careful with:
- There is an offset of 1 in the height because the empty sub-tree with zero node is considered the smallest tree.
- The order of the indexes in the tables A, B and C is different from the equations above.

Code: Select all

import java.math.*;

public class CountingSubtrees
{
  public static void main(String[] args)
  {
    try
    {
      new CountingSubtrees().run(args);
    }
    catch (Exception e)
    {
      e.printStackTrace();
    }
  }

  private CountingSubtrees run(String[] args) throws Exception
  {
    long    start = System.currentTimeMillis();
    int     size  = 51; // Size has extra 1 because empty tree is counted in the size.
    if (args.length > 0) size = Integer.parseInt(args[0]) + 1;
    long[]  mods;
    if (size <= 51)
    {
      mods = new long[]{1000000007, 1000000009, 1000000021, 1000000033, 1000000087, 1000000093, 1000000097, 1000000103,
                        1000000123, 1000000181};
    }
    else if (size <= 101)
    {
      mods = new long[]{1000000007, 1000000009, 1000000021, 1000000033, 1000000087, 1000000093, 1000000097, 1000000103,
                        1000000123, 1000000181, 1000000207, 1000000223, 1000000241, 1000000271, 1000000289, 1000000297,
                        1000000321, 1000000349, 1000000363, 1000000403, 1000000409, 1000000411, 1000000427, 1000000433};
    }
    else
    {
      throw new RuntimeException("Build large enough mods table to recover full value using chineese remainder theorem.");
    }
    int[][] result = new int[size][mods.length];
    for (int iMod = 0; iMod < mods.length; ++iMod)
    {
      long mod = mods[iMod];
      System.out.format("Processing #%02d/%02d mod = %10d time = %8.2f s.\n", iMod + 1, mods.length, mod, 0.001*(System.currentTimeMillis() - start));

      // Compute binomial coefficients
      int[][] bin = new int[size*size][];
      for (int i = 0; i < bin.length; ++i)
      {
        bin[i]    = new int[i + 1];
        bin[i][0] = 1;
        bin[i][i] = 1;
        for (int j = 1; j < i; ++j)
        {
          bin[i][j] = bin[i - 1][j] + bin[i - 1][j - 1];
          if (bin[i][j] >= mod) bin[i][j] -= mod;
        }
      }

      // Compute A. Order of indexes is A[i][t][k][j].
      int[][][][] A = new int[size - 1][][][]; // c, s, u, t, r
      for (int i = 0; i <= size - 2; ++i) // At least one in t and one in s.
      {
        A[i] = new int[size - i][][];
        for (int t = 1; i + t <= size - 1; ++t) // At least one in t
        {
          A[i][t] = new int[size + 1 - i - t][];
          for (int k = 1; i + t + k <= size; ++k)
          {
            A[i][t][k] = new int[size + 1 - i - t - k];
            for (int j = 0; i + t + k + j <= size; ++j) // At least one in t
            {
              int np = (i + j + k)*(i + j + k) - (i + j)*(i + j);
              if (t > np) continue;
              long z = 0;
              for (int q = 1; q <= k; ++q) for (int p = 0; p <= j; ++p) z = (z + (long)bin[k][q]*bin[j][p]%mod*A[i][t][q][p])%mod;
              z = bin[np][t] - z;
              if (z < 0) z += mod;
              A[i][t][k][j] = (int)z;
            }
          }
        }
      }

      // Compute B. Order of indexes is B[i][j][k][t][u].
      int[][][][][] B = new int[size - 1][][][][];
      for (int i = 0; i < size - 1; ++i)
      {
        B[i] = new int[size - i - 1][][][];
        for (int j = 0; i + j < size - 1; ++j)
        {
          B[i][j] = new int[size - i - j][][];
          for (int k = 1; i + j + k < size; k++)
          {
            B[i][j][k] = new int[size + 1 - i - j - k][];
            for (int t = 1; i + j + k + t <= size; ++t)
            {
              B[i][j][k][t] = new int[j + k];
              for (int q = 1; q <= k; ++q) for (int p = 0, u = j + k - q; p <= j; ++p, --u)
              {
                B[i][j][k][t][u] = (int)((B[i][j][k][t][u] + bin[k][q]%mod*bin[j][p]%mod*A[i][t][q][p])%mod);
              }
            }
          }
        }
      }

      // Dynamic programming:
      // Iteration: height (h)
      // State: number of connected subtrees (c), number of unconnected subtrees (u), top level subtrees (t)
      // Value: number of configurations
      int[][][] C = new int[size][size][size];
      C[0][0][1]  = 1; // Empty tree
      for (int h = 1; h <= size; ++h)
      {
        int[][][] old = C;
        C             = new int[size][size][size];
        for (int i = 0; i < size - 1; ++i) for (int j = 0; i + j < size - 1; ++j) for (int k = 1; i + j + k < size; k++)
        {
          long y = old[i][j][k];
          for (int t = 1; i + j + k + t <= size; ++t) for (int u = 0, c = i + j + k; u < j + k; ++u, --c)
          {
            C[c][u][t] = (int)((C[c][u][t] + y*B[i][j][k][t][u])%mod);
          }
        }
        for (int i = 1; i < size; ++i) result[i][iMod] = (int)((result[i][iMod] + C[i][0][1])%mod);
      }
    }

    // Generate final results using chinese remainder theorem.
    for (int i = 1; i < result.length; ++i) System.out.format("T(%3d) = %s\n", i, crt(mods, result[i]).toString());
    System.out.format("Total time = %8.2f s.\n", 0.001*(System.currentTimeMillis() - start));
    return this;
  }

  private BigInteger crt(long[] mod, int[] a)
  {
    BigInteger x = BigInteger.ZERO;
    BigInteger n = BigInteger.ONE;
    for (int i = 0; i < mod.length; ++i) n = BigInteger.valueOf(mod[i]).multiply(n);
    for (int i = 0; i < mod.length; ++i)
    {
      BigInteger p = BigInteger.valueOf(mod[i]);
      BigInteger q = n.divide(p);
      x            = BigInteger.valueOf(a[i]).multiply(q.modInverse(p)).mod(p).multiply(q).add(x);
    }
    return x.mod(n);
  }
}

Post Reply