From f594da109f767347cccbe7b7835d7092be2adfa1 Mon Sep 17 00:00:00 2001 From: Cornelius <60546655+loyalLogician@users.noreply.github.com> Date: Wed, 21 Jan 2026 15:29:06 +0100 Subject: [PATCH] first try --- notebooks/03_analysis.ipynb | 33 +++++-- src/CA_model.py | 121 +++++++++++++++++++++++ src/__pycache__/CA_model.cpython-311.pyc | Bin 0 -> 10646 bytes src/__pycache__/__init__.cpython-311.pyc | Bin 0 -> 294 bytes src/__pycache__/analysis.cpython-311.pyc | Bin 0 -> 4240 bytes 5 files changed, 147 insertions(+), 7 deletions(-) create mode 100644 src/__pycache__/CA_model.cpython-311.pyc create mode 100644 src/__pycache__/__init__.cpython-311.pyc create mode 100644 src/__pycache__/analysis.cpython-311.pyc diff --git a/notebooks/03_analysis.ipynb b/notebooks/03_analysis.ipynb index a03ce92..2ee26e3 100644 --- a/notebooks/03_analysis.ipynb +++ b/notebooks/03_analysis.ipynb @@ -10,7 +10,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 1, "id": "237bf30d", "metadata": {}, "outputs": [], @@ -24,7 +24,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 2, "id": "14500544", "metadata": {}, "outputs": [], @@ -43,12 +43,12 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 5, "id": "db6528cc", "metadata": {}, "outputs": [], "source": [ - "size = 100 # width and height of the grid\n", + "size = 500 # width and height of the grid\n", "p = 0.5 # starting fraction of vegetation\n", "update_rule = CA.update_Scanlon2007 # function containing update rule\n", "true_frac=0.1 # 'natural' (equilibrium) fraction of vegetation\n", @@ -68,7 +68,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 4, "id": "285dcba6", "metadata": {}, "outputs": [], @@ -85,6 +85,25 @@ ")" ] }, + { + "cell_type": "code", + "execution_count": 5, + "id": "764891c9", + "metadata": {}, + "outputs": [], + "source": [ + "grids = CA.evolve_CA_fast(\n", + " size=size,\n", + " p=p,\n", + " update_rule=CA.update_Scanlon2007_fast,\n", + " true_frac=true_frac,\n", + " k=k,\n", + " M=M,\n", + " N_steps=N_steps,\n", + " skip=skip,\n", + ")" + ] + }, { "cell_type": "markdown", "id": "f831e03a", @@ -155,7 +174,7 @@ ], "metadata": { "kernelspec": { - "display_name": "CLS", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -169,7 +188,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.7" + "version": "3.11.4" } }, "nbformat": 4, diff --git a/src/CA_model.py b/src/CA_model.py index 2f7e221..b925053 100644 --- a/src/CA_model.py +++ b/src/CA_model.py @@ -12,6 +12,7 @@ # Import necessary modules import numpy as np +from scipy.ndimage import convolve import matplotlib.pyplot as plt from numba import jit @@ -50,6 +51,126 @@ def update_basic(grid, cells_to_update): return grid +@jit(nopython=False) +def compute_weight_matrix(M, k, dmin=1): + """ + Precompute the distance-based weights for the neighborhood. + Returns a (2M+1) x (2M+1) matrix with weights (dmin/dist)**k, + and zero at the center (a cell doesn't count itself). + """ + size = 2 * M + 1 + weights = np.zeros((size, size)) + + for dx in range(-M, M + 1): + for dy in range(-M, M + 1): + dist = np.sqrt(dx**2 + dy**2) + if 0 < dist <= M: + weights[dx + M, dy + M] = (dmin / dist) ** k + + return weights + +@jit(nopython=False) +def compute_local_density(grid, weights): + """ + Compute the local weighted tree density for every cell simultaneously. + + For each cell, this calculates: + rho = sum(weights * neighbor_states) / sum(weights) + + Using convolution, this becomes a single operation over the whole grid. + """ + # Numerator: weighted sum of occupied neighbors + weighted_neighbors = convolve( + grid.astype(float), + weights, + mode='wrap' # periodic boundary conditions + ) + + # Denominator: sum of all weights (same for every cell with periodic BC) + total_weight = np.sum(weights) + + # Local density at each cell + rho = weighted_neighbors / total_weight + + return rho + +@jit(nopython=False) +def update_Scanlon2007_fast(grid, cells_to_update, true_frac, weights): + """ + Optimized Scanlon update rule using precomputed weights and convolution. + """ + size = grid.shape[0] + + # Compute local density for ALL cells at once (the expensive part, now fast) + rho = compute_local_density(grid, weights) + + # Global vegetation fraction + frac_occ = np.mean(grid) + + # Avoid division by zero + if frac_occ == 0 or frac_occ == 1: + return grid + + # Update only the selected cells + for i, j in cells_to_update: + rho_local = rho[i, j] + random_nr = np.random.random() + + if grid[i, j] == 0: # currently empty + prob_flip = rho_local + (true_frac - frac_occ) / (1 - frac_occ) + prob_flip = np.clip(prob_flip, 0, 1) # clamp to valid range + if random_nr < prob_flip: + grid[i, j] = 1 + + else: # currently occupied + prob_flip = (1 - rho_local) + (frac_occ - true_frac) / frac_occ + prob_flip = np.clip(prob_flip, 0, 1) + if random_nr < prob_flip: + grid[i, j] = 0 + + return grid + + +def evolve_CA_fast( + size=500, + p=0.5, + update_rule=update_Scanlon2007_fast, + true_frac=0.1, + k=3.0, + M=25, + f_update=0.2, + N_steps=200, + skip=100, + seed=0, +): + """ + Optimized evolution function. + """ + np.random.seed(seed) + grid = initialize_CA(p, size) + grids = [] + + # Precompute weights once (not every iteration!) + weights = compute_weight_matrix(M, k) + + # Fixed: use size**2 for total number of cells + N_update = int(f_update * size**2) + + for n in range(N_steps): + # Randomly select cells to update + cells_to_update = np.column_stack([ + np.random.randint(0, size, N_update), + np.random.randint(0, size, N_update) + ]) + + grid = update_rule(grid, cells_to_update, true_frac, weights) + + if n >= skip: + grids.append(grid.copy()) + + return grids + + @jit(nopython=False) def update_Scanlon2007( grid, cells_to_update, true_frac=0.3, k=3.0, M=25, dmin=1, BC="periodic" diff --git a/src/__pycache__/CA_model.cpython-311.pyc b/src/__pycache__/CA_model.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..ee35e26f3db7a97585b37114f929c80c94b23440 GIT binary patch literal 10646 zcmb_iYit`=cAnw;O}#}+c05*PTeKzFa-7(fo%QNr#}Cdft_** z?XU{}v$;4KVY3ieH@havsc9XjO~9#HFr8b1}CAd^l;1!W0shBADb9_SJVp1d;cAbnTV`5xP zDDtptfQyi*z(o^LCCbO5w|FHgC7@ImuZx6_)$tg1ac~X%gd62$V52{r6;lef?30yf zEC!soB#|OSl;vQx$@ws!h{1$Jh2vv`wLD=YDIs!GDd;I)mtrX#o|C4yr~<>|1I|1g zmWY6hSERbmk5MMh2HA=uEG3}#brJ5(SEI^rHDOsf`uHjAUf5*F4#PO$-VC4b{tn1T zj8ao|)eyRTZ4aY`V$?ItRpdc06!nOyUT9TTWv5IVv^S~DEl?SiY3%j6`7_}9z~G!V zsVD-LTYr7B`g%l}8npA z){3McC&LEI{qzM*@-no8ZRhdLP4bh`SX7zCmq2^ax4^5MG`WnOpA1wL4v#V z7&j$os#xUcEwRETgO-XZk*wIXxm9f8X(<{OD^{Kmel}<(Xe@{mNX3GZC*x{Cl76p_ z4PUt^i$uP1bw*0c`>ueKx#Ij^LIsCcFC3@?LJaI+#v`B3o6n1phClRG~P#c}R@ zG!F8QB2%py&}s(OY6e~>(p536T#-rmN{tp(yoLgZp|KIt43lbuk9-=)9P>LHIyF9*LmS2*{^*B1*HE`D{YIQC{~ z?9GC^yXbzazyzKVRo&FBkF~?D06H;h^ku&Ez`gj~h;^tu9j zoi?S-DwD5~!fGwkU?S5Nl^ti^3DaKIv`w|BR)M8bXS!tr*|SZVHW{7McD)`-)MwQ@ zuwJ$49L&g}uD{c$M-J=jU8-Gm+`>91r7omtx9S#5*BRmvU{%O}%5((E=dTEC-Ls%F zR15xzO6&1dB210Mjp=B3niIruQ~<5GG0iL7G@ne0D4`O_vJ3_zCQ`ddI}Rs#!pAvT z0X-S!#CQ^A9Ly!B5K-j%hBRxXw@QQ(Vl*;6DG_;?i%vl+w68FL%)N-c_H#ykKiDd6 zQVfF~mNmBhoCG>TZbW5KXHO?U(DVdNkInG)*aC#=52%h+hKY+6OBiG#=plZns`zm{ zIiyG-&5DCk-Gt9<#S=~u0=7GZ+EZa~R@m8!DV3~P1nEYiVu^`UO2tZG?n=d7vkKV* zT+w1zFVcj=cr@-IumRC9>4lg9hwmPq zn^y!O;=^0_jmRU6Z^Id?AKyl^o; zxiGP~E%#c<-=DFR{oT2dJewc6XI*SrP^+a{+x)qE?TbD4Jh@P*X&_@QcLZ|#??m&x z#f}||?ZuA4jC;l1k=s#lbMRl@)Rh^}HfKk&ZMRQlPOUUNqT-PjSP^{zVch=eWrrJl z{_X{;a42}UKnVpC6;N}X0XNi6HLK>S27{lrD5$*pr7~_%UX9CVptuhVcLGZJcG~)b z+m5Dfb#pN0PSqLYi7F+0<9zE~P}&Y^We=Nfnr7&%h)rc?=o&~z!o)K3dsLHZzvTuV z>@DPkDo**R&=UFS>NqxkS+!6F{52Kh*9ppcG6l#2J%JDfjY@>Y0nLrvK($e13RAGU zswUHtBn;C2RW|~DY~PvPLqYCl^#)K2iQeRHM3rf-%BN3=M-zK+kl?OeGn&IfBy&q7 zpyUdjLKxHm{6e3uq=F>MiJb};mQo1?lv5UCQ$cO)$M_Daza)u`g%*g-nm!5^|fZtWk+*tZe+n)^u18> zy|Btytgld#F+FPO%r)mss6`9MiY+~*mY!9{?1X9{nbAkB?RO94_T{cEyp-RUzjp7X zg;$EL!BT55GhX&LW%p;&ne+$n$qgaZ{>r=eEcM<${BWqa>qyoPyw5Gz@+QzmTh3Q( z+WGye9r?Q$PwU)Ss#DMQU-Pg)z*dIJfOu`-w-hvl`c*H$Dl4#`n;H?3f-0+DSB<_Z zTi+M_qzU{Gyj<{vD7$+1~YoG8EJ*p3HRgciFW1Trdfz%`6uLaR1e@qNA(i=qfw>WpC@M-O^=y z!T@=Sq+sdzosVhjxchQ$=bfW>kIqeii*awssYTcBxv>?uZ|;;Pnud4)e}F8SrT-1c zK|Qn**cvEgQ4AZ4LaUnV?RUNHUJC-$M(b-s%REEvzKr)K?VbSmKMxic{J`H*aAr?C z!0rM(ss~iHYOyBf`Pvo1g0!nX;!V3?)(!lsGa>s_Cn$yufuNv&qmpQFtkx@CHM|*s zGzf4S)X=TERQC)Wp4fqtJ|cy0ad>u8sA&T*iGhyo4=@HoLb@tC#xE%(G~arID#vR{ zB^rn5L*S}`W=-VFrKsv84X-5_z6JeAjY?_SK^=nz)oi2_mPqAP(+5jj z#qEb4@?Z9S5i9ndE(T6#PCg2BePYi$7l#)4VqkA6us3tE><{D)E&F>5{@(JI7qXV@ zxg49laL4sM%%rP3_x8e;#e+-H;_jng9xU$q@myePItq+@ynxpiNA)CunAXmUN7#UXtq9X-I)Gsg*mWDPS z1x<9d9v~Y;XMjWyzBc$mui~w?17)<=*McFsllIkm>Z^u_-C&J~ch+6h2zcww#Q(Q{ zV}2wJ$u%%vF4YeZ-18Y2!v}^z1heN;9W&Hyk)5hLU$c8_lT}^$+IMPH)SJ-y9RjHQ zq3T8qacgITxmn9tP3cjoZbNneXANPxspZgZhAfH#%`H$USTRQv3h9HU5Y&NNk0%fY zgl8&tGzw4$J{jX#fdrsUut~EhIRA=UHwBp5BS&$t-AHy}onxvRGTP4qLoq1>qMeB* z$pki`!KIaEUNo2PI9qMbCa58kF)*xxmg=|W$bGor>zEBgn{{{Km6o z*|R0b=iBl(i=MqD&)yA+SFPsVt|yGS+x-+t#t%lasXOy0x&8UAOTmY47X8Oc{$m-- zirbU5d~hgpsO)L}#l!~_xwcX|GFMVlyEd;r0gpQ4xTL)Vl2r07cJAXq#g){9E z*fcoqw5QI}SRgF5Y!H?{N_&5B+uwl2sDVjq07dl*W*SWTwC`txrG5ZI799+|3jw!N z^XgMM&q}O*8j2U1}o}EC<2iI zT1v2rkc0BnmjcMpXH*@t7MP)RrTskfe+KKY>g#xA!?E|$u@`>y*bbVd(^n0-Gh2Pe zt!5toHveC+1~;xN-BkBhs4c{eY?K@oGH!GSuurpMNEl?bVI1ws8+X+b{BDrc7PT3m zC6xqtYj9wp-OHv^%zG@|NmuRe88CE<+5~pIS%;zY+kl9zYKz*czkdx{I#u9;urJHb zlV@}HpCTqEVO5vo9p8$cms2!T{v4}(<5fQIHeZZ6kJ$ivu zxhP}a0nbhh+^HXM^dUCWQw}E+5N&ca_b5V^0@8`#1?X7^HKdxucn!`#ctLj4K^Jc(ast8Da^az0giE8gFUHpf+4g60xjLL z1KwZ7Y0Q;ehUcK2dN$2NQmBp2kOYw&K!Qj+*rvswh|wBc9qEU%ifMF=Y{p8owqyj! zD3&|n_#lJ{8i3XbI7xu>5s8pNtlxuVFOp#-KSeT*=<+}`zC2TYC$wy&AZ_mfyEtFt9Y4}oc`=w&wQqGAH0OWs%m)o}&+PANCLMk|$&Za+vPksog z_E*pFnPcGmS{5$lm4(+9`K6xwQx7foLxs&pAX9y~>j`7pfn-w`4!ILIb1wryNwDbq zql_7oX)ovci>*6LtveU@FE!skP;7lM<9J3C1k-H_Xi2qvVDa3qw&dv_6ng?%_*8Ol z$u}>%dkXHJa-cgmyLhm;^~Hzm!xLY4iZ2{1299M;u51tv{qKr*Q*Il)``(gwYkp|i z+gtGVLION@exWn}){+xLt?6R#tHr>p&tqalLv=9@B6}L+EQ2BZac84{fEed_F-(25 zvp$B&0(yZtSFbfcwm}S1i_)r{Fx}&UC<^>$eXND0knjiHeqi{@v`;lDh?Mn9M*?<; zLd=*-k660rTo`ECx>vXYXkHHSI8wbBMfXoegH+ta+H=+U@0AO-vnjR*gOR4cXhc(YNz$J}d&=ni# zXSf`uTj&vC^$44ugI^aTBA)#T`r)_1$H?0tKjf3(;h-Db@H6?)gLAGgg6Y)+ta@^m zBYZS5#m8bngXK)^)C_NqNCEl34$D`h0qUEgliE)nz;mNlhT};%+>5~aVeVRX8UjJ^ zU35XRF~FJH0SGTdQ!*TdPQeZ!57<%cUwigKz5jC$w|3#$df`O1FDc1#bTT#@ ztjYk02?Yr#DFY}t#`N^SQ=23>-${fCwLj`GFaq5W3`dQ=7^@nqTk8v`(2b7VLOCh$hYYbp2z!=(1}k+pY%!MDK|ll$~lr{c+3PX9fHmW$)`G4teR2^WFClpo2~Aa2(+50>$j=G) zokqA%1nBz&KJxcKAkyFeLnHl1{+8#5^R)#E7u@~j?p=!~m##hBR_HpCIh{FO_69zkSoZc6 zynRqr=-LfM{w|F83;x0UZXnwKkB#|*CW3(#!T2c1YKdYic^ivfMREcOI#Lbdwp8q) zkRXLap^6m*D>`nw%X3Y#!$$qILq)49A4Yn52O0!VkOfu42=_&B+*?>Z4%t z#Dcu3{Pg185E48}A;*zmN{0Ll3BI!XjCSyij_UW0BM`Yj!ssU?4v*+9<$naaYGzq> z)#PF=kX`}O#jvdfrbYWNGd+d+|1z_+Q2$?M+6(pnt2PVU1c#J$#K)MdbFCS%z_u0H zwi4U6YU;!WK%STxQ@c^Je|^pgY^SUX8mgov1IT1hTCv_AZrb= z_;bZ1wI?YMa@s~H}DgTVkM)bc}jVtd&lR43}B1QIkl3N*SWN~4T zJ-J3hU=o>RzK;vnGA~svpCP|5==XsT+9@GE$!hdlzs;`rBSJvFtewY8e(O|cp^eC^ z*l7d&^iE<-jy=f0E($(#);I9Nsy>Y^tBma$8I5Zqcn>TOX$J}q?@hF1!V?S%0jv(|eGH5k!6M5nQavN!A8flEin#?pyzg(K1 z78pX~%+M71@@6`+O4sdQL?u<9G%N$RL8I~@|=ZbC8L;yqLG=Gool*`bvtJ|cJvWZ^iq=XnXxYz zUoe$f-GmWm=hE+W~(e&wAsS&C8vGz5>Rg=hiOe9R6_S) zbiYl}nCm>5#Y0ZZJ}}bQAQuDcG)`lF$_}+fEh0x~570t;wZ=bhzvs0ndDnSoqSs4 zLG&K%e!p#Fu!pb~>Ri#nwODXp%1Q3ynWjl4rgN8Hp*zo5z4F5nr^UB4UK5`Y0|c#W z`<{WDY1``_@ZN6=0Pp3%C%cP>-!(1ArXFbLvRTv6HQ6+5N4Dl<>RWO~&pWbX$-1Im zlPz#nJn@LcS@G-?xUG`5a+wb9wh#k^T5?*+6b{H%#w^HL*vM&l`<#JckAtsHTzW*j z@ewiZBjQCPpDSLcCcb(bS-5n1uN(H{8jU<3`7L4rj@ME ziSL#FN9(wSp3S8N^I zGYi+PtbJ(atfl7Cdd9J5E;~8R&}T;9QEZSU7jy;72a_{rtaR4Y^YY~akkj&IBb_tp z!CEup&6)A7nelhAwV_?C$%jiHZisQWYswR^)Wj>Uc%|O6qpDVMRlOB03tx99Dy6y@`PpPSw{om<>~n7I z+^5mI(e?K}?Rv<#dk%ed*q!=8gJXL_j}d^U1l+u^bf(<*O$c@M{bI6efAZR`*H%)M zR0HwCBn4%r9`C94R=!_m>tgid$VZXt$SQMdbZwv<@x;-ZIO>X{WG4PmysCQQu9~>Z z6?YMe^pWK5IJh>p{@&ejK%jf#p_+KerM=FF%gplmfABm0&ci7@epij(P?1S~C2oW>F$Dou!l1dXHzno7?DavMR&q7oZJMF>QatU43W zLOZFJPz;s0Vhn9dwqyl_9CInj`fMiab9z3Tb%q1`N{9mnBwXiOGU1A(wY# z#_g2}vNeLNy|5^BX8e*6ltHUPXff=7v0JSqltK>i^37HXt#$0vqJlPnZskLbE5bM^ z{Q-Xuta6{&hqgtrB!W(R7TFT#^tPx_FScdwJe%*RlpqNA!v%$kaNoh(!>)51v(A12a{0BCZftEY-A9xsl8bYNz9mpfge(uRh&yJd7mRyKI&j2ooX_N2_MOp7~&22JSqIXpkS;Ne8I^C@2v}7Y!MUvR>F@tc94|j!H~=MLtXAp zR7X}8Dho?*L9r(ld7OeWzmeFv%7608tyfmUWwv}~V_;}?c#Zw-zK4lF4E%23%fa6d zes%iG7rn6)wXqZKz)5f5WNqMNIacrA(?CoFqyQ+NtxGT7xxCi@%d5Y>>Pq{e*QEXH z2R-RPbL3H4g9F}w^8Kq}kKY4o5q%fkks}XnZ|CvaPN-Y`8&X&0?W%I8&)q%l4NcUB zCLZ>BLx((RvL;O~or8lUUf2-(>tehv4%TBYxwJPTp#(fe06Zn&=EbG6<>5vT5~T); z3w;}6tRlL-d)Fq{=UwrrCmyYdM**Sqy41TN4X*BRrIVj0?hd+RZ#28N_e5>)37|H6 zUFxYP-uSEVN9|9szr?)6_Zkr)5`BySJSCt(z&8?#NR?>yh$jx!#35H4B3nE=KH-Tk z*Tk1y@nt&i_6>XDNKG7Z#Suzakt$MIqUxvz5xE0Ei|L4H?24%3`}Dd-iI_GLd|b8MuaSJd_QihMz%1K(Tub*uXhAo zS91b9*GPlwV-jLMPyRB{_hmEMh;x2{sBYAezZG9kTUyT4Pv9AtB5x(T4ghYSVHzyU z@C`(Oh?qVXC7QjC_PD`bM=-nD|JUqp7X^A94ZFeK2%+TI`s7-^fqIx21yyYe@4)sk Y%|oq(SFzhS(+TmS$7 literal 0 HcmV?d00001