From 39ff2924f57e3a9e3f9e4d1803093aa4f69b977d Mon Sep 17 00:00:00 2001
From: schochastics <david@schochastics.net>
Date: Tue, 16 Jan 2024 10:58:09 +0100
Subject: [PATCH] added metromap layout

---
 DESCRIPTION               |   2 +-
 NAMESPACE                 |   1 +
 NEWS.md                   |   1 +
 R/RcppExports.R           |  24 ++++
 R/data-examples.R         |   8 ++
 R/metro_multicriteria.R   |  77 +++++++++++++
 data/metro_berlin.rda     | Bin 0 -> 11382 bytes
 data/multilvl_ex.rda      | Bin 9067 -> 9105 bytes
 man/layout_as_metromap.Rd |  52 +++++++++
 man/metro_berlin.Rd       |  19 ++++
 src/RcppExports.cpp       |  84 ++++++++++++++
 src/metroLayout.cpp       | 185 ++++++++++++++++++++++++++++++
 src/stress.cpp            | 233 ++++++++++++++++++--------------------
 13 files changed, 563 insertions(+), 123 deletions(-)
 create mode 100644 R/metro_multicriteria.R
 create mode 100644 data/metro_berlin.rda
 create mode 100644 man/layout_as_metromap.Rd
 create mode 100644 man/metro_berlin.Rd
 create mode 100644 src/metroLayout.cpp

diff --git a/DESCRIPTION b/DESCRIPTION
index 86d6601..0a87c9b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -26,6 +26,6 @@ Suggests:
 LinkingTo: 
     Rcpp,
     RcppArmadillo
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.0
 Roxygen: list(markdown = TRUE)
 VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index a572513..f3871ec 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,7 @@ export(annotate_circle)
 export(draw_circle)
 export(layout_as_backbone)
 export(layout_as_dynamic)
+export(layout_as_metromap)
 export(layout_as_multilevel)
 export(layout_igraph_backbone)
 export(layout_igraph_centrality)
diff --git a/NEWS.md b/NEWS.md
index 0536860..e895f58 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,6 +3,7 @@
 * `layout_with_constrained_stress()` and `layout_with_constrained_stress3D()`
   work for disconnected graphs
 * internal code refactoring
+* added `layout_as_metromap()`
 
 # graphlayouts 1.0.2
 
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 9d1b7a0..7f4a637 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -17,6 +17,30 @@ constrained_stress_major3D <- function(y, dim, W, D, iter, tol) {
     .Call(`_graphlayouts_constrained_stress_major3D`, y, dim, W, D, iter, tol)
 }
 
+criterion_angular_resolution <- function(adj, xy) {
+    .Call(`_graphlayouts_criterion_angular_resolution`, adj, xy)
+}
+
+criterion_edge_length <- function(el, xy, lg) {
+    .Call(`_graphlayouts_criterion_edge_length`, el, xy, lg)
+}
+
+criterion_balanced_edge_length <- function(adj_deg2, xy) {
+    .Call(`_graphlayouts_criterion_balanced_edge_length`, adj_deg2, xy)
+}
+
+criterion_line_straightness <- function() {
+    .Call(`_graphlayouts_criterion_line_straightness`)
+}
+
+criterion_octilinearity <- function(el, xy) {
+    .Call(`_graphlayouts_criterion_octilinearity`, el, xy)
+}
+
+layout_as_metro_iter <- function(adj, el, adj_deg2, xy, bbox, l, gr, w, bsize) {
+    .Call(`_graphlayouts_layout_as_metro_iter`, adj, el, adj_deg2, xy, bbox, l, gr, w, bsize)
+}
+
 reweighting <- function(el, N_ranks) {
     .Call(`_graphlayouts_reweighting`, el, N_ranks)
 }
diff --git a/R/data-examples.R b/R/data-examples.R
index ff7feb9..92dad25 100644
--- a/R/data-examples.R
+++ b/R/data-examples.R
@@ -3,3 +3,11 @@
 
 #' @format igraph object
 "multilvl_ex"
+
+#' Subway network of Berlin
+#'
+#' A dataset containing the subway network of Berlin
+#' @format igraph object
+#' @references
+#' Kujala, Rainer, et al. "A collection of public transport network data sets for 25 cities." Scientific data 5 (2018): 180089.
+"metro_berlin"
diff --git a/R/metro_multicriteria.R b/R/metro_multicriteria.R
new file mode 100644
index 0000000..7d59c84
--- /dev/null
+++ b/R/metro_multicriteria.R
@@ -0,0 +1,77 @@
+#' @title Metro Map Layout
+#' @description Metro map layout based on multicriteria optimization
+#' @param object original graph
+#' @param xy initial layout of the original graph
+#' @param l desired multiple of grid point spacing. (l*gr determines desired edge length)
+#' @param gr grid spacing. (l*gr determines desired edge length)
+#' @param w weight vector for criteria (see details)
+#' @param bsize number of grid points a station can move away rom its original position
+#' @details The function optimizes the following five criteria using a hill climbing algorithm:
+#' - *Angular Resolution Criterion*: The angles of incident edges at each station should be maximized, because if there is only a small angle between any two adjacent edges, then it can become difficult to distinguish between them
+#' - *Edge Length Criterion*: The edge lengths across the whole map should be approximately equal to ensure regular spacing between stations. It is based on the preferred multiple, l, of the grid spacing, g. The purpose of the criterion is to penalize edges that are longer than or shorter than lg.
+#' - *Balanced Edge Length Criterion*: The length of edges incident to a particular station should be similar
+#' - *Line Straightness Criterion*: (not yet implemented) Edges that form part of a line should, where possible, be co-linear either side of each station that the line passes through
+#' - *Octiinearity Criterion*: Each edge should be drawn horizontally, vertically, or diagonally at 45 degree, so we penalize edges that are not at a desired angle
+#' @return new coordinates for stations
+#' @references
+#' Stott, Jonathan, et al. "Automatic metro map layout using multicriteria optimization." IEEE Transactions on Visualization and Computer Graphics 17.1 (2010): 101-114.
+#' @author David Schoch
+#' @examples
+#' # the algorithm has problems with parallel edges
+#' library(igraph)
+#' g <- simplify(metro_berlin)
+#' xy <- cbind(V(g)$lon, V(g)$lat) * 100
+#'
+#' # the algorithm is not very stable. try playing with the parameters
+#' xy_new <- layout_as_metromap(g, xy, l = 2, gr = 0.5, w = c(100, 100, 1, 1, 100), bsize = 35)
+#' @export
+layout_as_metromap <- function(object, xy, l = 2, gr = 0.0025, w = rep(1, 5), bsize = 5) {
+    adj <- as_adj_list1(object)
+    adj <- lapply(adj, function(x) x - 1)
+    adj_deg2 <- adj[unlist(lapply(adj, length)) == 2]
+    el <- igraph::get.edgelist(object, FALSE) - 1
+
+    xy <- snap_to_grid(xy, gr)
+
+    bbox <- station_bbox(xy, bsize, gr)
+
+    xy_new <- layout_as_metro_iter(adj, el, adj_deg2, xy, bbox, l, gr, w, bsize)
+    xy_new
+}
+
+# helper ----
+snap_to_grid <- function(xy, gr) {
+    xmin <- min(xy[, 1])
+    xmax <- max(xy[, 1])
+    ymin <- min(xy[, 2])
+    ymax <- max(xy[, 2])
+
+    deltax <- seq(xmin - 4 * gr, xmax + 4 * gr, by = gr)
+    deltay <- seq(ymin - 4 * gr, ymax + 4 * gr, by = gr)
+
+    xdiff <- outer(xy[, 1], deltax, function(x, y) abs(x - y))
+    ydiff <- outer(xy[, 2], deltay, function(x, y) abs(x - y))
+
+    xy_new <- cbind(deltax[apply(xdiff, 1, which.min)], deltay[apply(ydiff, 1, which.min)])
+    dups <- duplicated(xy_new)
+    while (any(dups)) {
+        xy_new[which(dups), ] <- xy_new[which(dups), ] + c(sample(c(1, -1), 1) * gr, sample(c(1, -1), 1) * gr)
+        dups <- duplicated(xy_new)
+    }
+    xy_new
+}
+
+station_bbox <- function(xy, bsize, gr) {
+    cbind(xy - bsize * gr, xy + bsize * gr)
+}
+
+as_adj_list1 <- function(g) {
+    n <- igraph::vcount(g)
+    lapply(1:n, function(i) {
+        x <- g[[i]][[1]]
+        attr(x, "env") <- NULL
+        attr(x, "graph") <- NULL
+        class(x) <- NULL
+        x
+    })
+}
diff --git a/data/metro_berlin.rda b/data/metro_berlin.rda
new file mode 100644
index 0000000000000000000000000000000000000000..8168056904fcf88222a40b320e7729005e6944ad
GIT binary patch
literal 11382
zcmbVxi$Bxv`+w&`(iU<qY??Xd7%4S#*lQ(QmKl*#4nrYDs1({TY_>Us+F03~njCTn
zZy~3g&&SHSq$1={eb+zm`8^)L*JFF!_wHT0?`zlnysrCsT{nW@Wpwhqy~<^K=#SL&
z9sY&)|Ng(8561U*{C__}+B+`n*i-jEe!nMmhsus!*ZzmZt_wSk?AWROKO}al>=1ge
zBlUks>`4Bn#fASRA<S0X<#puZksTLz?66ST({OE%ck+>o+7~E#|DDw9y>R4$O1xt7
zh2!-*ymqU^JlJh^Za3m-wne1LqjS45!Nvnif0sED#7p<SGT|FOuoR{0Imr}tLd9EN
zsdp$Mk5=Jg9xIv4#;9T0fI4OytZ19PE3TDUh&c>kT+PWxX?Q*onUAyO#j0DeQ~@YF
zH(S_>4aW!pWGEa=XJc)6yc{kPfSQ8^efAX&f^FU?jysQJAxJbVLs(J*Pip{BP&HvN
zoyNvMa-(wC(L&+~ES(Nu*;qO<mx|0~SwJP6@U~oJ4nr8u5W>MBY+D!ugs{ax=;A1I
zFdc}6VUR4JMoFkDIBqYPp38ztf>}@i#%@F(Rf?08EFj=q8wxP4I1F^Jog{?s64=Ax
zh~wGy{gOEqY9ivmCp1uq$C!)IN>y#RE<XLOg#>mZ55~B{1M*VTZAyaza6Wl$EvUfy
z<-S}-Y`!MUhVMF)7bgj4fFvQMYH{M`dA5%9LJaJ6DK0MVN-m?3ms<g=fb)29aTuSO
zqaqR*q%Av&;7DL#38mgHJ_6-&ff_ap3494(GDTb+YEI&-VVy9p94^X-r>-<(>!gGR
zb6)^6q!RY>3>nFD($d6p&07oA@m5W2QcXo*0}w|OjO4W{ZQvTQwk-2pQZ)k%19%*7
z+y@?_9)R1tjE;p0t7D+-Pj*W67;iNjBq#7VJ62MurzcjzQPVum4veQ0@D+hJc)B^3
zCk!i*1oHzm@N6ytnC2X=y9VRU@Ny&L8Q<w_5qx~~v{Iv#xez{rV<!o3vaNTqd7UjC
z#SP~s*|ex@Ne$Wc;{xNMS2TgBTt@z0EH9sNlnEsxHbNye;5F_{seVmZy^}K@7g!%v
z#D`a?hSfJ-5ILH*4|PRXH6k$H{5dL`qN<Z7<o*)QkH8(3gf-dn8gg2lF^FD8#K3Op
z2Q_G6tSewQM2!qml+-$@>8|e1iEnDg;eMJ=?TUW{4aG?o;vVppH{)mw3@%3EsMJAP
zKJnlqAu~mv#(XuA)T1yrAbNFXrcFv%{A3ywkZNOG)0Rr5!w7g(IFpi8lXs+D=Zc2*
zD|uP{@;9{j=WfvGUr*GwDw-(IFsrGlc>B)WK^X(0Q5~$Fsb7R`@s=(+*In}Ym5e3O
zUE;g$RJTe&3(u+`W~p{xom@`p3`y$cwU=?ST})csj#xFt>R!>-Z;Us$gZ=fs{z_O6
z<#P8*T2fER#o(lYQn+rKh6tV*_kOn(bN87bBXv=;T`L!evAQE1>G-C|EBp8y_w<K}
zLg}SMM>o$xCnC-jL)@*A%Drsu1McunmDy(o<>7mDrN09C9&2cz#Vb}<<~A#m-s~ea
z{MERLL$qiZOuudrrXNUla@UcWLbl?&(v^>QbArzGNsF`)+gGQc6gDnv@aQm{Uy_$T
zz{+UoEtX36j32bs@YHK(Qe0;_ZdYA=2r<-*Ufy7;WV#9`g5%_a$NAWQNtci@sCn6|
z&DZd`qNncXfE<te>VWLhv&iq&k%>|HP~UD+1rXyO))zEzRrtv2JTkghmm(Aq_z29C
zKu0(^1CDCxH#I>G=}eXV2Rbkrp4tT(&{t9!I&bwp4y!`dIs)Jva*(2>Nf)f5QBTOt
z4IcoUto~>oN0EhQy?2FOe_fstKz1yZG0$tkyeo8nC4EptN}N!~ZY;bSS#@enXBp1d
zJYSsF?!(lfF{Grwv=+MSq&HorIk{-Rl<VQ~vT30+a`ikqK$(eb&E+i*H1Xmnk(!tn
zfzcpHg<G1t0><r0;L={(-1ivvtyW2>Yb>;#eqEc#kknyn1hq4;)w5$GI2gTjMbSJ?
zW=fF)S2v(VaX_$SeN`u%Es_k5Vq?il*yLRV6(H~ZuFYBfRUA44o`<NRdk?6i_+_Pe
zOTk1ix($!%c%$ozf0^w?;Gq|Qm6J3)BTw7~M#dmHNDidX9UADc3+GKIFf+5AJ(%!1
zF1?yPXyUv7`AwC@v9yDtRGfr47YZgrT^q{Yoy!w$U=MsMt;kAcW|}KFvdK8LG$CJe
z(cmtL>xG(*OduLes?f*`;k}^RluAlqT%cqO7{HPtbd;)~0n<=QSTcY@=FrIij9Qg0
z0d~UYjtvyt5N$~x8&PSf9;c1svBMJ1u+Tl_D)GaF=VpkCaA&{S`7vuMrkFRMkVMBQ
zJDOir=GrVz?Vx<sWPO&`QIZljew3?jNQ<_aurW@yni!e7f1}L-^URWK@)_Hw1^eIt
zrUYn+@XU~RS+xq?g-_lNX2zA@?!bh)4Q53yrZ7+AW{$!ct8E?LRn+-k@1&*=K9)>o
z;M&Rt8BTr|%nQ$SRw~?<R=s2#1v)%WdkEV7?CN>Y0dpgL9+M*?l(^;1s^}a+%+p*z
zhogDofQo?;D3TrVxB%nk^|$wQbt+R?#r3d%4@g=Fc|O4fbdZUE`}>Ubo@*(4wZEw7
zJXd<(eJ}Ujfha`~nGCYxo;*xK=zfG%we}n?M!r%>{#5P(R|1LjJXDkc9g>FcK_Zlj
z6+s7DytT__^iVOGqW|r(->-9%dr1C3NlZNHwPng72_b|{Qq98+Ns)QrD9QYzFiYus
zHDX#_>4US<^6_O-1t=1P-2`Eka=ErjB=85BGK(GcdikGn2gQ}osO5|^Y>+vLOuEvn
z{!S4v2?fcC<H*1)5?mDz#83{O0}e#pJuw>e&H5!N=nw;G3b-5?8*V**uY`RZq5RFX
z?MP+_MQyfLmlTpn=&Co6^8mZoVoB#3l0XaKNRS9P3M3mFa^k`2rJfnt<U39=g%{F{
z+^L(FKFu{c&me0a)p^vc3}`8UT@E~WcTz`r(d5pHE3$Eccgy`}2C{D(zx^-Vu=gI$
z&ayUQAHSddALxLO?CC25Z^HV|e#t-GJ6#bh_oNHL=E+x>$I{WDqxuPNs%v)7L4Fry
z{L=aq@s59g1vL*=1Re4_xeJS=fmF|{BA_rZ9ZVJiK&ns(2MnP8wNUNHox#>dYYt)V
zrl))J{R^xmQ@5;7UKWX4lDOlkH||t1Js$cqA$U_t9Q~yhGkR8;<>%M9oHDF6@He)p
zqo=Z`YN|6=Y4F1K#ei%5Y2?sGarjTo{y&UhP5$KBxA#IW-J86+GE}8JoY`$>efN$J
zx;BF5u;mH;Pfl3+hp2Mj8C|y!wBsX^scxT1GNS)%ugBeTm6ay5&WtTT|Ad`2Z-R;g
zZsgR@jogWumsjhSDJkFk+idL3NlYEYw2FFr=?C4S`s(V$`A_8o`>P#E^ON5@gSA3#
zH3XwEInKvdA1iI&pwCEt7ZrJSGjHeco4ln{DV$qkb=Ce_+8Z6+*Am|Ukk`F`qr-ia
z4is7Vl1wVk3LBlUXd)$qcls9oUK^KPFH$-tSGQIWFmwu2S9J$mw>Y+o^k}O3)xC(z
z1E$-zy{aF5FnsqtJnHWG<w?v7{hV;m(OX4pO{B2v`IL~*6Ds8k`kFt8+i&hpn=x<1
zIj{CTmA{$T_w!rSw%D!FOkL-DFOtKLHoD6{v{|ucF~7>YF^CSA-7K%&(EPq!#QrKm
zhcH@#xG*wd)1MWAKtZ(<aImG8pxIs0VhQ7238ywaj_rB2(AE|OBOCV__XR%(Fm~C;
z%xAM!8>^v?h{imL(Be6|;_AqkWS^cM69^EKRtV!Y=Fd=eHBw||u#8+@ES%59cSC*p
z=q^q;*F?u$UVmIFmKU$a%;5(K$8uDy*ckJ65nG&VW87DGuWbc6aA2SQbiN}Feo_EJ
z4zfs_SKU<Y_>zf;J~{x$sqe{cz%#(SY2I4CFb7_UK{7DW>fWxX=9;a0x0>CBl(5=G
zsSJefm0ZAvY{j$U!}L{=5-3SXY|g)rxWhCFD?b<wMuP@EjN8@6f>Sei4IHSQi0!!%
z+NozPw>p7O3cCq64?Z_zxtB>c9E+Tv#1(;(>dp+8Dy*{QzkD`~DJ`zXB~(!MrBQM&
ziSs|!xD8Z)(!6nqQG6<9Urq1n6O;EEb8Z(pCK!d<NzSF(={)ToF|q`Im4iRu9DnlV
zbJh&9e9CBlNU+$szlpIHLr8@DrqT7f<*4&w%%q{Gtm={?iEpVxjxzc)prwpnXV0Xe
z#ED;?!-SI0Jqs;&AW2b%c8Foz?GxQf%DW%kPk9+4^51CnzHLh>BnJk_CzJj1BYVrA
z7c`lMG1oF`dwn~MG|$XBg=#9!9r6eztrVV#4Dt)OG}4T(YigPvfYDt5dl4s@IDBe~
zY;TvYyy=ej{@?0|(}7#P`k(`wFFYelOym4R*ZxdS;(m`lNkeayTfzsW6YwjVj)U(}
zJ~KtW8z)RWCA_J}q$>`k9T@Av?IOi2=2|kt5K0PZL)<X;KmM7A3chD4<J&~$A7!37
zQg%b@xw2lA6>3iuRY(y6Mgm~NVkMns1y-5q?nJ)}Nj=Xz4mTj6hixXv0pXP~fFuX}
ztU~x15_KSnL=`%qCxkfs5&k%EFpKeE$uWD~U@rMOD8EffXOFJM=3fzYr07l<39ur2
zBZUoA&_YL*`S+DRiGDuk;~_h9D<sO&x9jQEb;6y}uop~Xdyh7YWo7|AKvQl$DTg+F
z<)|<mex|DCrl0p5-sCn_toYl(Zkb>1iZq#B@(86U%gvkL3OoG0Z(b}d%Oo!~lSuVZ
zQLN)6Hzl#Gn)&g~MB6S<di()h3^?|<mlTQQH66(w4zzqY0qWEjeKB!1PJh@q=QRV|
zfz&1Jf{Nw*0;tk1p!JECWo)y@lbE&86WN5hF@#QYn>S6`%G-t|fh<z!?n!v?gA|j)
z5(n$;B0=a#2u;a`3MN^dbCCfn9rXH`_C9j*WcD3XSIQ1_4<D3nQ3>%?q?tu187N+}
zGH*uYup3|kdxKz5utgI{Mo%3AA<@9+14&jWa1sgxK?xzjLY9FI_f`tL9HzVo?Miqj
zkJc)~>m~k_OtYag|3!nYM6bk}*pb+V`?-PsEp9u(l;6{oaSZWyq44Xb&ELY=**5@L
zqZ9w#YgDLE0fQgkOb$eEagI@8f2?;dJ+8xCn|bQdUsdkAAl>>Z&Z&QK{B)I*xpl$1
z&1nC<z=oF`91gN1tl<Pi14_9qQn6r<vtD$5xghDa<!bqFzU_?7b@EbZR0&vQK+UmC
zeUSqxXWM~kFD$E6o=ePxt7EcCLBpIAgyc7t?~tC#gHTY9LN|C(=&n1`;w2u>%_B+{
z#m7hzVxfhZk*ep3FO{~q3SFBWl&|C<O^|Ei1g@NI@O)@h7?v0^>fEx-R0qJ!@1<BS
zDlk?E1*UVrw152Tyw(T4%@{bEtJIfW7oa&aT#Ih49QgjjDSBEw^u?($XHG-49cg^T
zZvO5&o2jS6<27wLK^+?#4pqZ%k~4_Xm-<GpD#wIeav#5%JyH9mrzX3<%p$r-;XJ+H
z;PGG6jY=cWw0n{twzY1r9bz@X^4bjSx`|Hqt79Wu1E)PZLT>#0AoDK1W>W4l)5qz*
z#q@Xg7My1W^g2JRET$a(Rrr3O($L_&$Lo~oJ2#X*jCk%Cmw92a+{yef`iks5exxG&
zM7qC5#87pwOOQrO@2UL{zAGHLd&@d<ZB#kzmi6kkv-fvF%(Hu{@k=3iH4{Bro0k<N
z`sPC|k#)~gX2xMF2vEXC3HY?&h|i(31?RmgR0#v&2myjo5<Q0wAh2YVsuBtc6R=h=
z9R<5n1mfB91wB>32Eam9k!%#|A4`C$LIlML083|`K(LS;u%#_s08T&%3ZbD?ku-t1
zp#qyLVbPSEiSngq+mIig&J39n>jVAG2P#rt2m1#!l^fJ|IN!IM4zobdR-8x}&7!te
zdv*rqb?SXP#-DNat11ZiBKKs|z*2otUK#p*RCl8Wd{t!LpxkymYh&_kRcGz#iHZ|r
zzFK=<nG*AT6Qp03d4m5YoNPFi^2%SVdQ#!*+F094&7+FD7MTriE+hI(!+MN8yY3sV
z_-U<mglLs*7p-j<Y5pzzI$|*RDO_@PaQyUKWem6Btl5tefA6Z-?<C*f)OxTeFZ+y`
z7_ug|xV9ue@Y)Xc#ZfJRm?*lH^<zmcn}eTq`p66Rr$+y9rVg4sc@x}OI6szdMdt_!
z?4$|-urvYBqG7&ZtCJBad>k?iSr~zvaeTt@i>;M$;Ms)#e4~SDIW>Deh=kwngv>Z^
ze;gdTh#;&WegcU@DYt2`<ooDHsQ+@N!K481d-`*GTGXxRserzvQ-GQcUyxgXn^8#4
zzl{N`RU1wAKf|wU6#SjB_teZ5h6*`j1k}hQ7DDI90k%*W&1zU?%QHOXdPM*jr3wX;
z=m17lkd;7TbdG9Bhmko7$%4@6fDJv<ZHZyAn6i(-MF|EVLs`vi!yJ~Y4V`o1_}@ke
z1eRVqwP{-1fCo^5i6Et18v$^jvY#H35y+e*em$HEF71PH!M^QYz3^7_vm3(Vz%thA
zy_YyH`&Mg7oA|%X3P!_H8)0mlC2>UrW;o~H(i{OWX#iG`CvjlOHlq`m%MS=_tN@Ah
z`?k3>Y_gM07xd>7;T|qB+Z8Z|WT!!|(f`4XECdm>Wz*UI`XUmL<iF5}VVIe3_|Z?e
zy|eXfW6QTTrGXEp{+jUeeP>(qQURTG2~pWCSl*c=Tn07t1cds`&8Eww3Oc544xaBB
zrdeR|J^ws}1|dTa)7b*Ea9}jSPGV82V0sRK0h0kN0EGzjN)nh3*z?aP{`GTcSdze|
zNMsI@M*8P7|KvF!uz+5J=}3WL*#a+tkf<<ll7QnNS%R=6SOzAT2?A9WMEEGS5IV$4
z@Q(o-s3HgkU<7X~I#nP>B13lx)S*MD0^hO{tN{Yk&x1*1MOCEW1%f=plA?0R5Hbo(
z6N2tk6?igmoQ46ZqF`XEKm~|EcR{7k3Q4n}k%a_K4oMPx3nBe;biio=Y35|hwzZ)@
zdZ`FQaFJ{qz>1I*3k3-#L#Zl?`XI=Ne>%{g3D&RzYy{x~jM~eE*5)`yg_a&GSjeK<
z-mS+fITR=SwmoxXr?@g1x&n@UmU7(c2h9(42EAJtOPoh4!lF1xYIqzbH^ApiU_Z!5
z`IE87CG|<il}>f@HG~6zxBA%a{a&FIQ4f8avwox2a^;!@^piU18~7P&-bp5w#|p%f
zZDS01w-RQ2`oKTAPNFU5QP8=bznx1%7b*W~DH){Lj*M=#@ueY{kl5Fy(Z|duX>k~_
zvdz3pKY#nq?y-0#K+ANWHmF*QVf^CWmq8_y$|-{XC0^a@ID1qWq<Av&%f*l1?YfT|
zLqVVeLHUR%{B)+<>#$W;Yq3*fMNQu4lP3E(w8sm+@Z+UPr&5w{0_?LcN#8MBc0E&1
z$z*A*47Rv7b>J9it{;gXnSd|nP+eqq#Vr0(+XMRD6I;704y6Tys@gnStc^o(6}-)A
z?>BANkyb?-{NKCTgm@^#FN)639@=<W&Fa!$=Y;7e>~-y*nTh927ZnwEi`CFft^~@s
z^H_o%2BTm!n`e(%D6vYqr43yHb%=a!%qgG-9wO$%lG(&$qxzuJV8h)6%)WY>>00Ez
zK#Nue=~N#T`MAX#)MQ{lg>A&b#<A<bX}gE~wo|L?ck{@Y6*+7o?^k<|sjS61%u_KN
zA_(w;fF^*U*L?(n<Jg;smrbU?j3gs@*RBb0LTAh5RPCR`<O2$skn{bJ{*prTV-YX#
z{iOxK@^YWgm0T7VU(JjA4Cng@Dq2^yc)SAFyn|}(A`D_aQKYT?Ai9-baJLbUZAR}t
zW1*C<`6kj096d|O`*NYAs7fD?bJxzluI`POO!MNW^VKh?r3}M01$EH(_-iii9Ash0
zL8c0%*j^;nUZWtCxWP~JRlhJGNf@wKYcD;jo5DaVs>eu309;2V`i7q9ww|lrQe+0#
z2L&j<y2(2*BBWdfttvF%%><4`l)*?FUoUU-EX(d4?Y;dzGbt*E&f4s%77yeptS*is
zat-#pn*8`3%f6(IyA~pPdCK-{&tnE)BxDqoH2HM=mGQem!sC$SW$59M&iX$*Y7*T%
z7d;Jq4v797p7^lYc3|%5kM*aau0mF8V%mSvwAULKysD8mwVIy!j0qu5UlBykb*}cF
z1QaZR`$5AXI`S*_(fe8*l11|z!NXE7`KLkf9<LtDEI{NK)NS;1aGZG15<es%6Jk~j
zD1sCgApVCwA}uOcTR~YPwr1U@G(KV`tM>i*o^|@8ANcz(m4RXW*ZU@33DtDV?3)*E
z$o-fQnLFbEo@3xAj=pQR41|K~7#y)vIP)J^oz|C?J+pd_+oVLfua}xx`z;*nqHgca
z6V*Xqjl?~Vey!BhgGOd+NTtY!T{js?bj0K${m4kUN{m82&G=}Y%~=^pf)R|=J69rf
zJ^y^8DJ?Rq!<+EM*B^fr;;!F+jQqx6cHjO2%fr=w6Dr80AA&<~lhe-V(9ft#vmjiI
z>N2AH`WV>VW6fPuR%F?<61q4tF?)YH*7U%d@Yxmw%I->U8H;^xsz%7FKArIqmZ%UL
zXW>nI*?Rq?y4fy_PeCKSUr{s-;s_jkXw?VZ|Jy&_uQa|S-fwpnKi}^Qenc|!ieJHo
z2Dp!}NhCz`n6gh6ikt7*KARjP_QK7<$dF|VhX{P`tCkYja2nE<N6HgiyfozNEnvht
zIt7gLV*l>}?Wdaq4iUHJ+mH(9A-Blqi)ohA;qQurzcmpoEC*~TJK03(0)oTS-64r}
z3${wXJ4Cqv9tV}7ml8Z2)J_(og@0r*9E(5Yc0npa7;A`n6mXHSoNrZA8hR`sNIl^r
z3>nyR{#64tl9G!_VpB5V+PV^WDvN*#W!8^SL)R!F_3!ab*?Y4CgH7VJDQXlK+j!CO
zYJM>tBLkw-4H&a&;oIsLKh8RRow1L9!N9<wHAT{~#i=Owl2~nhQ8VDvLD8R`_&g5}
zeu+oH@GBjM&4Ks|H(C8iikdA(+{vaB!@T}d#<p~Tn$_;)<mT#zOK<V}){NtL4IgCg
zlQh6KuBkM&IHRoy)%9&a6DKr_pH!E3yxA2K89^~`EpYNa3hyiMVp1(8_*XYThhv`A
zWS=Cb^T}jdtK=)0nOOwt1`k#|um`;}A`fSFg+j^5Gpa4bzg<XE`efg3E%f4=MV;a%
zEhf`7hW`w}@{zRcdC0I4k-Btj`zJh&a>VJHzJ^aDVRJdo;US@^l3+poZfpAGjZNLX
zWfCQFrcY;}G`=`GO8?UQ4DNhrXyCF<w)LurWPz_mWN4_l)BJWj%!?M{ZEV<WE_Uu4
zc2FUznh(!EBYm^)-2ql})&;gCs$7Y(>s>vBFDep<pG#MvhDvJsY7|6Q6-v<z6(cD=
z4YzV?vTd`W>hPx?G&ca>Q;}57w#A59Fdsx^<5~L`?m`bo+f!H(mjm&v=_hYtV-D|E
z9Ivb6BE)ZKEJV14c@wUiZCPwkA|lqFOv%{8p^O70#DJu0V^oo0W1A$HdWGQ{Pw>cF
ze9j?4D+RY8abKR?r4Y1r>_wz#pTG^1qPhrYd{+sa_dV-%|3M+-EzudXZRj!R1dMd&
zf`=HgCznJT7COAnJuj>F`|L=Lx&8B{l+Wg*$BgSRQfa=EkE{#V1s@Ax1CqS#lX-rv
zg#a24=0_iHN}j5PlaRzY?QM>g<r@^GC^eQf1!?oW@TQS}YMLydH3D`aAVgRI$-b~p
zye!#OF>A%oO)=DGnZzQiiLT0gn!b)iA|yE)9Ogj2047jy3!IIY>&~^!14ZXLd+H2}
zoo#dR#=2lic^=s!8pvALqRLqrtV~nWIy_P6L?B(+SBFL~@zzO;CqC5Fr5rgKxbk1#
zZd*Qy9~o-x?0&l6sx)MoBCA1v`PYbK8Nc=UK(IR-wDp8WdbHcHhewK$gcj_?P9qPq
z&8~Fiei`hNI%Y&;Nx6*{>{D^;?WWP8$|lX8JO+}EAmfT*Bho(VfP;?tbfL4`b*2=S
zxu<MTD%P2hy$kc1#t)5FW5TNW(LNk3%az9fw%j5-^1j%?Sn4sZhpK63i_LSB*1r5S
zc<UXyZ>YMy4}wm<JNV4EL;kHHcq{XUs9L+WNp^cN!Ror+ArR~4SK$pg-A<jpAg*Wy
z%Uw5Z*1>8ezL9X}@C4mUYIJ9it!0K)Hr?i@sKaC8q2^|F&nK!yPqJ&!)tGra12ZA5
zoB=3|i-KW~7chJ|;kRs*EGV_hu(!n^g{pyZor31-`2kxZUV<8J9)2j*)}oZwuB{IT
z^+HG#CCuUPVA+Qk8OjAA`(_h)=A;lI<Vk8db$?VhoTqF}iXvD7zZQj4X>oH@pK+^R
z>BWoly;^C=;@B$N(hJeC<)xpb&_Vzd{855v*G-GRNLUdGz_=+!B5^!VUK8D3z@l#W
zz;Lk85{8%Vv3|COQeSGah7nR*9d+W|>_qn~Lk~M@kPyf=xDh1pit+DevU9ItlA>rG
zjwK1Kwp<H0>uP2}iH?*_19Pz*UtLkE1*>M;tF=m{GNEv`t-eTHHwROs&54Mh6B?rk
z)y#c}CjA*ZSjhmqx&`BtQ|OE)`Vob*cO?W|xr<N5_Y-ZM1&+eY7LI>+h1V2SUEk{L
zb@;IoAh@Ks=IY{Jr6Lh8d^9Pcdrek(T-~G?6Oru?hcJEcwmvc{sTASm6$?I-DM>cZ
z6;L!En+{JFC6HpgWzwyCa36qz(4Wz1QO~GK0uB$R*{7zC+}3-t8^*?{(~;sRScN)g
zMv8i@QJLdoYo-Vz{W_2=41abl#IIoaijulNI)$oW&(*l)+aDK4*5f*Qr5h%h^fc|$
z-3^Rf%<+-=T;LKY)kC4w7x!P&=SdAnx6fjgQ=sx{pJ2*BG8J_Uqbzf2gvu}5Url(K
zkjb$oo^e^i8uvXL>+vO9Rsza|U)G?@K)jm0AD8P<DrTIGtqV8&Vu3Xoo`5s#o{Vgo
z#ChE*x{@93cW~uTjdD`x-q`|ncTvVLkr*MF|Lc;{Lxpi3X4`q8NKzfiQTGuCY~9R{
zbvI{@yln`N*;D)Q9V-|v-T5Ird#!3vYPHZWqOJP&(I<q>M8as?6kb4HxwVRyQ@uxJ
zo^zH6F`3_n#=;r`Ep|`p$cAOzL_2)C_`5otL%7A!(#E{y^vLN%I}8(hd5&x>^D@kx
zjlwTnJxOUa_9$pcujp6a<!Ed-gfhuKSJw`!;pL^pEnsN%*@~3NXx9W*hYpq>dLFUS
zlu3oLEd{kFzg;9AR<N}cW#O6jn7ejRJ}saXgqT^S$PfZ;Ve}m4>kc3>mZcRQGJ|H|
zLd~c|FgrRH9IGrckg>~LTt`M<dO<pmHkc7a=rpOU2Llb+He8(T&G$Ms@1c#kzM`3j
zG~}e_x=hHhl><@q&`7z8aK95=b-b?y*X>I3x#cBTD`Z$kB?wW?{{AUuB<1kw7sor%
z3Gg9NMj2<`D%$c(sQ!Vw<;g*bw_ZNpM*YlUn4LVa<ewX6$af<Ib5FmkwIQs08_KQo
z^}mQ|fh6**R=tn?H)#53xyGkk_DTEMo)@r>lU7`@T<V4xmX$$0_FnnIl2|-<>nUgn
z$f&fIUlms_N69Nf?tG2@5u0%JGhL_0qAt18>Ni(`wJIOhfKQxAF0k@L$>jl`N82V$
zL(L@|EqhS=KcMzCCmSi+MNMBll4O}p_!-v_fJ|C+lzF_48dV&8J<z=D+t`rvMOrg@
z(BO{AynT_%7~Zlc?8pW@Vy{EwALjmT!-fzsP(!)?+bql5iGKpVMtb+=#55;o`+Wwz
zg(RvZf6NjjqfU&uL$*IXPX&2u2hKh7PxKm>=MSqBvhSR3Y!%mnHka&AA|}{>>&Ks7
zn~xPWhVaKtB&<5ZXGDF_Pt(r7fl!E}uRS;MA%|E)<%Tjn#wutFrf%oSW1-77tfo^j
zk=1A25aXGf8)LO^v1^CEZcGMee80-AEKSHMS?>I_YE69y8xsofAKn(*%XSSowEFGo
z*ZMbZrQt7rg$|b7e{MEcb;6(j!#>9CrJq>+T9I#vCS^1cR#B&yKRATnj5XX3UOfCG
zBRTWz+d5aBM9^O+PN%^;@K*HtS)$yNAvoTn1U4;OHEQy9LQ5{_Hl~l&ZG3y;=l)$0
zbF}LU)V7m)q7jy_{J(|#I^<Udi+!`X`CkRp+lwd5J~U4}6mVrZZ0dtT`1;J@<!0G{
zESq1#B|?~Ma|zd)&{#CoBCLn{C8Z-`P!h(cwfYb|U*0K*8grJ;BR)4{gd}1=$R*S7
zZ;Y0g-yFU8<C9_ePosCPVJY)x-fq8|2<<g>@br=ost<V|8o4dzpWS6_*BYuan?3j1
zGWF(#GReC6iT5yGQ>))M&}0?kJ=R%j%;CuGdr|%!()DHg_V+}nG|gY#xp?UPdM3HE
z3UNjKZqQq|;P>w<C2r)D&yCzDV7M{DnL|9&Ba^11$Uk{#Ytdo&y8nM=Z_K9NxXvBl
zy2#iV)qN0hp==`Ja*c0Z=UCHHchghV$RC#6yU}xUf-332_|t{KH}zY?XLo-%?fIQR
z3I00%b~7I~pT)5>;Qx3kZSiMuvfpmUP=<7u+{rHaoh7N1Uy}#lDlSjT|5V`CzP3!5
z>HeLu`)9SC_N<iCPWrlLL4Vck?9U~?9<3T;O&G&b*3I*c{h#jVj{G{4Z0jWf1CPM3
z@btb(FFhH#Sn$1rvM*?}TxI&3BX*_f6qC4VB3^ZSh~L0zJNszf5bPYeCabn^Hk;?-
zgZ)kq+ut5qX2=S-WV?RDx?v^3V|60*Q_ZWH%6GTwhU8_fJ`I@kZq&l3e9mW)EsZxc
zq8lrkpb5X{)2!vm0kh8EKTU3o8rq#&a4HzC2}%v|{1Wcja<|O2587Xoyk1KjtiFO+
zi3e6*-P`{n(5GbX3HpagYrkW!2QYIXr^+k=#@Ne|k8NN8A62zDl3Hhp$Zu%}981-&
zqlM$Zu0R8N?O;w?^g;83jVK%>wc8aa46ur)3DQ1}T-}{_&Ei%Ys9XMIOcBk?uFr9G
zQP}8C_$lnwoy&1pb)?b@r-xI|A;0zlaDNW@r0$^%`YY>9+3rwH`{l&ISxWaWp|HTl
zj-S5|Me~hIW!-#l9$Z)d-6rd@q|hDFT9q<*(^}|dez?W<+NFkRh2U4>!EH~|LcdC1
zulbwU?00=6ap%m)mbb;v8{8XT9MVE5)uzh534Ige_M53s$&s^iPh4@t&Z_=P%<nX|
zJk1oF_$+ygQ-Hqem4!^$JAS$x!56!};k*CZD*WtdxM5!)8QrMU_Ozm{_U-+<5`|K?
z5*@>u>L0h=?o_XtO1~3BEd2VXgLCxBUFVijeV~H%?u>Fz%@g;VYO?6+z~}Yr{vY2Q
zOY7GcRVJ(y4g3`SSY2}C)WU>|-H%{Ra{cSMo4*eUfmg2YFFZrcFk!tkdTY&g{Kz(p
zdc7)+nK(8>DX!WVHdPyx=#B2wo0DI0yixEwq<VZ(Xz2Bet%!l1UndUL)mJY__-a+M
z&j??7{kJHI;h5Z(>)Z5U)tdXQ&9TZaprgX=;~C=Lvw<sdL<gZr&YcU^_8RE8vr|cC
zf6u;cx^M7l%5E}rTwaYK2s0gLzr?%3RzGtxW@%q1Aqkn=>2E6|-dPi^7fN6aUnW`W
zR(Vr~!Ix|ovqFcB=4_Up+D~Ry73o~6+0@#eln-5zU%7Wb#b$7f+x7a+poxmZli(8_
z5z&bc5`?}ww@evNKaKv<El2!V{<zTHG~oN>c#n*y|5AF2<HGRNmFc};d@fZ3^uzk1
z%I91DC_CV*kmF3&;-lw)ee)v|dN@`&7ee<snkv;cYY<89JNUxpvCXf-owH8omxOSm
z`GoC8{9uRq=drB?_>cDvzJ8rPjAN>=PV2r35WBWGP&}EM{iEW<w}^LXgdwMMf5^Jd
z+;-MR_e_TD)XOT#KRq>W!Qpp4S|*3C;6G-K{T5qnxwT;YVlBcl++VmRF8yZQQe(#%
z^pDCH;7v`Bc9RoNLp{5VlOKfX-SGN3cNa4t_j$|7vlejnZYsMY_gw+`msfK-H*CEo
zFgHSSaoMu!{?R`{i6KAzTV5K`3{zCTkHJlX9z^V4)|K74F~Q<<Myl*ce>9?hWxi$J
zTbC*Asy2J>`>r7(=+sx|{yhH|JLW0&U6tUgJBj7@*4xmJF<ngf-0swyHv%Wez=JjS
zo=^7FuGHSvQu-`!wrp_p)%}fk<#t)Od<O%U^Rj|Vjcx^Sz+sPzZWo>@J>!P|RkYj|
z;`zg`gVI@T_D6ZGVtpgj?!REona%XALf6ppb1R8U;g3zmG2brsBTMJvx3xf3ra|M1
zf8@-H!^_u+OW!9m+R{TUw{B|H{N$aUtG6Imn!Jq<f4mw#c>UX<?&UEN{8HAhwTo^g
zrc-tPng`cGAO$~duMDND*}(0Xw91OI`h`12I<P-#j;Ei6mBar2s-mgZM@4GIH+20U
D*5Mhm

literal 0
HcmV?d00001

diff --git a/data/multilvl_ex.rda b/data/multilvl_ex.rda
index a95a62b1ce6e3d2577220110e3ec2a6066628d86..5e2c864c0ae520a1515914bbbeaca736a7ea4c95 100644
GIT binary patch
delta 8314
zcmY*-30P8F)OMPw*O?TYHK81UEa!5(4yb9G14u}eA(FG@TxMJ~%-PH#aj39xK%vx}
zP)ogHYD#8`WaKT&A!jo)%dQ`FzwdvZ|K01Xz4ku)S!X|IzkBU<&Z1+U6Ufdt@VJ33
zM%x>+4P4LM^R{4D0Jl7l9%;`WXrH*_o;^Ym*Z1t%v#(a5F3PFzp&r=dAu#`g!=8$C
zfh8gS*B%MsKZQ`<&f1i{+9JJuH<FE2j8B#h(*JLy&^>`Ag!0O&s?6UQ^|(TjUbr@{
zyQkQAkFZ~0Dj+>cSxLN0IOt|zm(V@=Kk|}MB;=-w*^hDH0a3@l_Wp_Avu9+F;U1xz
zuN?oBqA36J+Al0^Xed>jZQNDjm!b_3QZ^3kRZ6wj&MMZ*OQ$9)8yZU+Nbl9EJSlG6
z-*1RFekB}~ELEH(LU?&!4zDa#V(6FZsGTZK=pQleHI_=H9Np*bDIzkmhf?8oYScEJ
z^x(5$XcbMRYTM-0sCQpU(D%yxp9ibRR<}$nz1gb9#oEr&LV1;L=~hlPKj?{O-?=tp
zKiLybTkwb9v~E=ZlHs{O^L{RY92~{I-(B=m0;eUN23JUp(9?d*d#zTm7PdW+yLd$e
zK!UB{;7Ne(qHdq5N>*sc3eVtEcCYKJWmIRdp!K!{oRAz8;=2f)@7{swrfD4XtJ;O-
zj|!QRE`pyIkVo1@C0I~>B%%=Hr39{P?`J)ydh7-(-S-x5M7+FJ<bdLyt_CY&JgU>`
z&=|=g#7zWz_fzH}o>oS$gCUL@bF(|32;(W9Y%|)hz|w6nu#FkV7GX+`q98wqZhHgs
zkPSys2=Jf-0s|6n=trW#i~5gtuM9TK<lhiJKw8sE)@1>d^Z=(;Fp;HLMSOS6Yj`)o
zuZyKmmh1MnTLAfG8hGp`ubM9^5PB>6x!uUq#8e4QUM6!iLsiMpg9Wx}i}7Ofm{gEQ
z!}9`6019!Y-VXu$QRJ7_M)V^iaYWeteFY5XG0(rMS!}Ys0wGK30F8|(3JS*$kV*}F
z0*UODeW;@jdNEl@q-U~$U*o<4q^FV<I|T$Yt|tGj-7Fe<(*RQwz%vR$p>T0q%wZHt
zBGtKDqG5J?GI?2&nTn7syfw(wE((?C4j5yDaLI5KNb8mX0-k1>HqxD}K+DBYOXWey
z1=lfgyP5UGftV7IACv-!CBc;R8|6W|!+z|TRB=Z{s(7syabF<;QG`>5>}h}wi8p}R
zH{17d@r|HA>l6u@%uNcsyVUgL&}n$X;?XNZyhh?MWp{DkM4%CooneF@NYOq?IZ9{*
zat#St=Md=``a54wy$=iway)tZD6w%L*RNq%s#r_eGIWBNrEEmPvlTHV@PTLozEKBi
zgl-)!)xM6QQWRs24UGs{m|9ca-tNUl{8ot+=4ha?ITy~(1X2W3_F83TQreY`LgQKt
zFB;oSHVk3iV<!31S#M8H#io9A+UxCV%H<xh2Zc$8uwJ*S!Ar~Bvao8hh50%&@Q=Y^
zm-HvQ>+|hNQT`zzrY#&?zUgh~a=0n%7*<u;%*?pC`B-t?&NKNam*?tXIi|BkePe+^
zp%<H4j!_JQ2*YhXj+M<3qv6&m{z)|iLU33(jzNTDu~^DuL*hVNn?_}2QDrElxsKnU
zCRI|34#Z?8^|lc!!@@;prL%m|em`R?XWfRgkC1vAy|0f9mx8!6rKQ;0ME(%_!d@}S
z$R$<ZQuR8@V^gzZXJiyYFH)NIM@HY{dUj#t6f&)MfSQ68ZQ{dJ)ebWxTlD!nc71(p
zXbIN02tPJ9#%?N2s!2)uu?%Yr#Fe^+*!zdrr`InU+hkMPdPm+N-eA$8VV92aWfbK}
zSg++_Vz_9;e4yVbkH;Tve2d^qhn7hX`2-Ut1ztB_w_??#!x-+Mz82v*dDFwLHlH1f
zKK#r+e3+>ttl=b7uUcGh^7z>ECJ4FfeK28w)^<VObZ#$JN2u%j(i!QY60B;8QmdLY
z<M8nCNflp;#ob9RZ**jIWIFVc*)fu8Y1!c%P{puw|Jc|7jnT#BHvYAnp<1fy#>i%5
zXNTp6QGe!$9uAie@#TQH-0?R)m>SLF^B1iy$Y;~OU`r;?S@n(O+_uMN_g$(RVSl*T
zc_AS?Gm5hKxapWWQ;b^AS&TJlJA1>sp1aXv+>R+ZlXVx98B*W3dFO__Ya6M)XF))<
z_&Yx+rLMiR7hBJ!N*|GD)J!enSIu#5D_a)zw9?Lbr;x(V-j#?ey@XuVdcj-jj-%2#
z2R{VGc1q!vD0RHqc)oi}Ixn-ir5&gLb<3AC`wW)U$;nkU!zw1Emzuie%o%Z3acyI~
z)!=9|4M+Eugx+Ppbyhiv!JJL4vS{k$%{XpFf%-Pudb<1PBipe2;bHf=DDNab?}yBI
zDbD>_NX;aVJ$=yvmopI}CRsn(-qbn7Te{TB4|6w@Z>pbMw@_P&EegND3l0y9Y(9Fa
z*#fI!+LEAS(Kfki+=h;AE8*+2dCl#b=lglfeoKL|X6<J<AX513VSOi~!pV=ZkSJ_h
z4S#BbRkj?qLi`%u64gX)YU-@1Pw}R3N4vSasj-p8%aiKvRS5~$UhdS`X1Ubc;V}J5
z8ON3otT$~gx-qUXw!MVUbaXduDGw;cJ&S6Z#*7N)%lT~Ki(7V|^KXiaDx40Ff8A_q
z6xR&tEbSz9%{qp?fP{qDcXq`!vu?*7S05f;(V@g8(E3*{nway$mN`yBeC~QXpUC~{
zmc7zc$`9x2i-i;pHZmsHAr#)|$2jH>@sDxsd3{sik;nOS!OKqGoz*j^XJXr4QhD=n
zM$z}Rw^TXH5piWQa~7^mO_RJ>^U}BBTL}qm)1>ftibZL8gnUv;@4GnlXL}j^ji#2z
z?Jw&&+>tpVaW2Xzx*Zf>=Duab85mo_a7R{&qf@IFAGc1#b<B4&db#tM$(6W{@Gbj8
zWiQN4+sip~b(5Rr`mrI$WG!f&y$hG8OsibAD_}Wd_WM5+UTW?AbWz{H92@F3)2HLj
zB2C4#m2oaW3QN-mvvbU%YWlZ1=<LDuwzEzg;`VU3<DtXEg9H5Zvh=d-EfusyHisY8
z;U1@2{vHQs7NK3CKFGQn%+=>Hs93!szkD}wJGtvx;&w{LTIa>jP$hpiO<3zDJqpxV
zs>z!@B?V)ZG!s2*<XB6)>O^b{sofpzpHSa99dWU#{j5$@_NIEy2S$DycOfv&{4J;1
zZKZh&o4^=ah=?&`mA(zjUa8!2O{hspsF{v2S7+XiB9ATwT~R+})Hymex^D6IaCl+)
z%55z5LIQ0nXUny7bTzV_&scFvD4hBjdmQ)Mn2Bq{anSyST^|#&NgPKhb6_zNuKkF7
zojZTpuE%A+o%4x5v<@D8)Oy1B_8--5cER!<CD+F%g6I8F+@<}_CuH*^357SG%H0C>
z`}i_FZtb^wdGH`v+=EUgz;sb=C==u>R76NbNNiIZZV0QQI{9WIf-of@>{4)qpj<RD
zYx;kM)N~v>`Q}&QHMKN7Zuc!`xO1yQT$4J#XnepneL4Ffq;172ReH{AgVDT*i+hvv
zd4aePFt#cj1JyVrfoxG|L$)A0av_PXh5fW%T6Rp|7!KJnbJqWC)DeSbwS~NM!TLuk
zx)%Iw(?7b5N3tgCT*fPXQdClnl@&2^r%5$<DLTc2=rfd3pT|qFFe(?lf&XA9sr<Iu
zxZT&omMS7=E0AB6f~q@G#!2-=;SA#5`~?IRbX`R*Fz)rp1>17bXysp(n7K@}y14Q2
zkw9B2%v@UB_&`^Ry}8u9(LGLzp-iuOnuT2Kfu*PeO9);eg*mtw1NR}(z|faIT!b(5
zr9UDMn0+1N84y`99$;H2qU<5i(B8WA)30E*z7ivXNlR8AZSy#9OBOw-bQ5ozE~=nZ
zeENUed3g7p@<Ga984s}H%_^$mbc7c)@2ssCv=|7ksqV=54C9QUz%@YU-#5H63Xl~v
zARJ2a{AIQQyfO&LieH^npl%Jq1?uV5HNxT}?4W9C)d>3_bG?`D83gZ&a)tR;Pmhe$
z3aA=f0^q?F3;=j!Tt^-#-zTnXYbPP+F)h$a$OBZ7ni2QpXoUmiBJ(;Czcs?5W|RW#
z^){$N)MZ2qd_BpCoh1En^JFqX{s4t$__$j8ve#o;gfZ);t)2)cnxW^*(DQ1~x$WHq
zIykIX%FrwEHm(eNTp4t7ltLRBEmXw1_!u*0Z(pJrxu8Q_2ftWWd3Sw3A_ki?vo0;1
ziL(rIzfN2?KW#-kRmc7q34fQ6?XKhPtuu;Qu1_@E%jlfHWK9y2^f&!l2``FlZoqA-
z*FH%BHMjVuQ>nC7%O&4)DDL;w%bA}H{oW^vJ^B^bHcD08UO@ZT%cNvXZ)AV;y9l0)
zf346;d6VwDeE6AtI*&7R`Rxiw{9SCv(j*$w_NAoQs>s!{E=%K08-v=@H8_cQ7fL5y
z2x+5~!*VtTNh=ngu&ct0J!7$L#V%9Y6(5We>!?%y;T^gS$FX&iD~m~&T3PdvDp~H^
zJ@4ZwMZCA!sej|RG@8+V@-yLSpVT)&Sa~u(^l@{UzV_eery?w6s#KcayuM7)Zx_Vl
zN70=$p9j~lSsw9ur+7|z$T?XaA}uk=)5mE=KfVV>;$$gWw<1O`(|Q^OF($0nPi>ct
zHKS9jmaR**8(pU`v`<0Aff{@hyen(u$6@JVxXSR6mq3DkCcS@{L6nw#GEgddoxmEV
zR~uYqNb!s?wB7|gAJ6tJ8K<SJHv0AM1mcvN>Wv|LLDIS-?56L_LzpHjH)*%(M(_QU
z3esS$a$OoXvvo>af8Xl{Ks7G;d6U8E1~u_vwX94-zhOoa0bjo<H=qa;<`nYbO2$oQ
z`(9$e&ktmav-BI5RFxlvIjv|>r~x)HF|Uu6sHpC9Y<3<<@v9E@YtKSLHid1HGJCtf
zS5Pxd>e@G49L1S2|I|Iu&)~Lij@dOElEQ=WKa(pcmN}v_G=M97x7>LbA{s(BslT?L
zihJg@2vmNAk$WG8xge>NWm_|LhS!Z(W;VL6OS$WSXD;*LmTnP~7DaxzGlBG7aA9de
zKdbf3L!9GQ<ll8Fr64pX0b0f6zQ>)7HnuM9BcE9a$DJzMDD*c;eT)N^>THu-5(06W
z&DR##g^`-A_JfZYJsoL`h&wyUttvj6tot=ULQohAebZ6*(dYXj=v5dVNCou$D~1)7
zX3ixSYjtHzeNP(An4CG|1DB1jz4txs;hp!GuAl1!nVBL(r!{f1O?uBag;3kirga0G
zWt=2W;F&Z`X8%Cb=>wx#${>;-X4iu|d_Gw%@S%u)5J`T9q`Lkn(nfW?AUu%Rbzihd
zLXj{^e+Y)B!9wi6dYCjzhn?D}$Mg@4JY;@{sks}BHNLjapNCO@h8wbgJ?}AXJGB}=
z^0qV~U;mjL-dUY@JoO2;fcZD&`)v~D!l&AwI<*n3b4ir0?YXkcH{O$ccU;c$>%R@7
zKia@*AXA6YA?o9wcUJH0bK8FJYw7xo^TaSBxaQ+!cdc0)$!`o2WXOi{Q+6To{tr$5
zuIS*8ng!soGq4rW^lZ1-COrDZfEUfB(9%-TDzLOSMrrf@&xDCO;mtKDo!(@Vmv8pf
zrhWULmG^dFQjiNYcZ*Ic@F1~)D7ylVU<AG9il-zMg~C@KBBfrij#%RNS4Ptvz{tx0
zyV^A{xOUsq=UT;J0|4YGq`-nMbTxj@G*E<s>%m~aN#`_#gI|Ld2#=Gyv8JpSdGmpt
zmK>RI-2v&~t$?(Hg6t1DXaRtt&S^+OE%-TnAOI<zMkXOeYyoni_D-8?El)pPJ0)c0
z<M$zO^d}6{A0^xQ_UBJ?8{6zv6(Gt5h0@a1g)#HCuN+1pVdBStf41y<`nMFMJOt%2
zc~L7kn2b`Wc%XptVvTz$I3WK_11LB^!FdAH!L~D?P`o+=47KUT=l;NsbxWktw19FC
z5J)>xg?KVJ4*)~~P(YwuW}T<%c(iUB4Iq9UoR$m*(<{JGfH;CHVt*MGt?S^RC+7eV
zu{%(mk(LSfRAuXv9Z(J^FK|H`{Q&|HlSDw}?SScd5&+3oVLgS62etrY1`=@-PEW2l
z@TXS=m<j~QrKM${P-)WQ*O35#2$BS*Do6;u0|53*82<75s1;H8Y;Zk&|3T+GupRjS
z!@)oPbP&K!{6MlmM1I*x009mNj3?>Bg0Z%aw!bRypJ@Wz{AC@?I_69<Bn;#WzKVr^
zb=W?nsOvJun0hYW!h!`Soj}s{5H9E(cAkp|Y8hXE9!oftIA=Lmwoz<I;Lwwhub?g{
zZ@%BRe{LrFFIqb|^#>D1x<6~Ie^!&dxh)XjXHBiq=yIJCHJq79D+&514+&<=>g%(Z
zQti~g&p+b%4dk0YTn#~QwHIB%^4*TqT`63!TwQ;PkUsLt#+kJ7<oo?YcR$~KaK-7-
zm-|^?+Enkn<bKT!N(65wCMI4I=&Mi8>ECvxy$Wcak7AX(J%myEKPDhF*qXMol^D$v
z)Gq!<_{ZZi_8`ZGnz|KSTQ81}>zE3QFsm(|H7m{Q_+(Z(-tp<b8M1q#>6vxi$&E%v
zNAUy2mYSrXA1|h^fmb_r^7qr3U2G(a?V0=>&PKq=_{_Y30;o;E$d0RK7%`$}YzJl|
zxn$(GbNI+>eEVACML2sl!FkqiL&Kl4qQ9ZB(DB01dhZZ%LsMwYvVS@#;aSJS_^rZK
z*pzjOAs6p&WRO@q)Kok~G*Zs$^JjeuxU%vLLL#s3OC;#BQ|!;zW2IB2Fv<lFuK^L|
zH_r27)pP!KeYAu=-yf@fL-Ns*g$z5o9Qw)JSG6dr_3a*bccv}%-iRMAA7Y~L!$!9I
z>&3lmyXb7r_tci=9#+rTdF*+tS2Oi-&w0hkI+<z5`$F*@J=N8pENh=*>4m_<m*`OF
zNm8Dtr}t~?TXzPQN@R36#ObIGb$#>*3UxKN`Rm=VuUuSA@zJXfn|E4RBU_5~6{0k;
znZHOz+e!vup}3fdu+T7p1ciiPBrqVmB19emlfYy>hd2C|`hb#HdQB5TN-#Q%$XaAh
z{>JFei^ir{Zf2c;q@-U=O?$C?`CFkA0I6&kc<ES~E0eXNBOj99H7#)bi{%c!&jrji
z-aK9&HoI}T^R%+zab%S)%BQ*c8LY5N;CCjn?F<tQbH(~XW*ho**BbhNdCd$*sb+X|
zr)ba4&P?E9XMdxpxuK+)UK?DkbQC{1(3C#|JV+7~K^OLZMRj~VMW~*KUOI?;ToV%H
zHI6{I$Hjzs#upTPt{n`F>X5LMuv8D9Na=njc;s6Al&n}vH|=Xo!8<5(GMiC{J$=-#
zJ>)mM=j|slm_knSjQ)9ZnUs6of;G}jd^QUuM{A+B)9B2{);~7x7O<YW<!{~XBkdY}
zE0^UbdVd~a?}W*wPvhHl)}N|!{&iZe9#$;`Ak%oG?H&xbLHYm?u!A6RGZ-!y5!+}2
zkKwEgfp}}PQFuQ-bkRt^9oREJIVSiyBkD?ix$Vqgn$6WatrZ}+A(%nW%shLtW;Duy
zIN9g{UWC@#8HC`F27WuKq(0eVsR;~<8bdki0CJM_#uf|^_n93lDb-Y6Gl6gCc6EOP
zT`9V@cg;o^vSv-Vc6mcHW)m>Beb*2&YvcCi?0x^iO<4bA<LuQJlyioiPBkAc8eX|W
z6KwQ2t%by*s+t82+Lc_-1$52VmkTCVPJHY4Te7WzmM;=lkC$_*OdrihJ$aA%_W=40
zn66opmhH?^xUul)42p{pr#nb$xq|6(0Ex5;El(ji?~Ruy0e17UL)Ri3pDMl(CXHNb
z<3yFs8RCc3)`VBg3Ra^u?E4>T<v)1ubaOKobGCj_L-4O9C*|I}xzt%FSb0zTy^2gG
z`@v|C@gqVb#6+eg?SWmMgo7SHY(byq#ky;?tn89S%hcG{0`bIuXEVq%WO6eiH1)M>
z5Ty8RNc#ek;pG+2r{fTZ1Rx#KyQL?9Nm+y{?R|XyVGvmY1-hl>oJMa@Nb}NiCWHU0
z^r=7OJTVQyn68UzVnUhPCwr>HFq2b`8)9iov6jwDXoGs5V2l4<&->_PI})wNNr7Pf
zUDJE##a{iJDIr=JcDV9Y(BtmRgG32U_}<|kXS1w1TBp)TiP5O5;y+{yKR-Q$@4{R{
z6ghZ3r!N*m4HK&-9jhjB`OEs;cMjIC(v_Z!e!P0b=)G}^m8|_p4!Wb8I_*^U?VRL$
zOHM(o<?8isr;^qHYt{;c4cLOl1*qjS|DlHlOHXR2(8Yof_s-(vQjOGX)N%T%?#V8`
z*P{E(r9&~6Bqr!QQq$~R1KFAuYV}-@L*&*1ofg*reL8bgxXe8+?NfX%U>Y?<ff%f?
zeOl@0MKJl8g}LWNdJk}dgcwJCAMrlXTOZ`wQT^nc;xN8~_9b%>`#f{z^Y+!Rxq<7Y
zQ$~I9<FS2n<{UBGrLq0c+_<kWQs*-F0OPRsOlUaCOf|ECcvTkMgSfJLj(hp2?zf){
zyR!NFP{A3EuWJ>80~@l#Z@9j>#JG2Q{F0e;D&Y@RxZ4-0*;glM%;shWmgYMiK>UGW
z*45Qf<Lek0OkLeh${VYe$GM?jpLgVt->8K<ew~n?Kqu>+eF;R`9K4=bD&f}bRlflb
zkT8BsYULQEtA%>@`VFL*pf2C@{)g<j`OwD`)&G&u@Iv{6hG6<r#)x>nOiR`hPQxT#
z&}YtRaWUsKESk2xQIknw33qrRtxJF|>|3U)DlhAJxP0vZuiUw7886-^P7MBg`)7{h
zsV4(lnwP&v6lOf7q!>0=WV=8%0R)8i{^ul+sR9cYgzEC%iI~8k5k&S*PwRMG!y}f7
z4b@fa)C>yY1l9XCE`pw=gD?b=dEIcF>b`J28n9&e;U_8dY0|wqE(UjM9)-zIiu_1m
z2fR7=7%dC~rMgngMa9*md@e%#fhFRpRr{AuD%I^AIvRXUM^$WDuLMX{#Iy~Pe5921
z=}Uf#hPKIK(<lcMiP3=zfrNE<R3X~l^VD|v3r_UW_<zzW6dYS1t|LDyn?n>ZA9UA`
zP2To{fyPVJWQZr#(%492czxcv^C{E^^vWtpY<~u>Kin9S?WUWv1;CiLorbuEOa~@9
zEkZsdWOlN7Pf_0=axrST*+kOq$!i&qFSy)-&#jt?E64d&^ESsnC8RNRXqz4bpRBZR
z(B_@qUO9W4>qLut@TGP^L+8PlJNctmP@e>++K(ImHJk)@Mib{Ad?~t`8#I@&qWCKQ
z&D$3cYW`@Vf^e3a6QHg4lCeIaVTr8}&x7>E)po3!wU19-Qj8;&&s}NfYzZ!MwAxBa
zS)-ZJJ#%qel-;msioJ0loSlcHuH30#QJvy_j|f@@TU9+y6J$3{Gl2=(vN%Ell)pJN
z)lOdQjJdLkGcYh5>wT7OU^sEQxVWI8pimDIFf<c+T*FK7bs21;l70B34<-+RGB&P{
zE_vQIgq4XQMdb*h-|=4|@c&dXJ3T$4{;1=Tq+YVkp??opRds;$G0EL00Uq?KS5B`8
z?6{uB_g8hhOF^>a*()I2JZ=sC#gsW1;OB)6z&j5%4`Q9)1q3+f6}<NMH!K@^sN#qs
zlbddU5#;N&p#O!W#B)hWNp+Hx+l;=MQDF}+qVra<n#8+fGJ8#k1Cucp?th=>`3>`8
zEr2}c-@me#zl&#&=lD~G!mbUVaTH2PqX9nW_9JG%-4jc)=g*y3vLfd#X>>1Xo;Y`+
z=1Y-1dC5A8=akg9rXefhR5*i<?5HaEd<S*;l-`VkrR?6WHA{6t$R&+_z|X{h92+rp
zTs&L3?uHA?J-?=I!93@F^@8~g&#@SFC&30CoJAW9O^0JjjRXK4v0!}FRNXiyxe8h8
z9Oqu<dTx>P9Q%@sAzvi$8CgDd&zpg#O@Q0;bj%}q<Dpcxw<x(VWuu_M5MgW(ilgSQ
z(2K#3f`PHpFp|QLBT1jD12J1JHizv+L+-}ddkraGo<!D~)_p?;C&px2BaR{h&+HUc
z>YKibRD*oD=B2-_a{FpT9-#g+?4#rD{`T?<*&F`==3@oPg--9r{UzA6l_dOI;ku^$
zw6Ut|6UelQHAIkbufwL%z29=P*J~N;CFhja@BR|ff0aPbCdM(xzVUM=<ja$X1(|n!
zrUFJi;Ozwrm*IpNTXcx2ty?tbcnj3{P?Xv;gTC3w3PVFe2m>2pa!V)jScWv>k5Qe&
zRXXl@oRh_yyeA2BEs%Y8#8j#8;x)ulM6Q15S`!$5t^I!a$IjaSS^;;`GgWem{RzQw
zYGB|bPaW*88HrXGg{w=9yhGay&Rn3NDP|yeU|$uaaHs)dWI`n3*_-EZ3QWWu;*rad
zsc?;$5n4u;5c;ToyGX)nScV`gxc?94!V}yoJU0=rs<Ew-8@+{EseNnpjFq-^D!RWb
z527l;MSRI3I%JVDgZ?OLDu<`q*Pl1+>FE<@*9vI)8?y2NeczeHxz>iI$!y1prPD0H
zC~njor{hW+4fkm|(f2!>cFz~kGTR+(ls)+2Qavu;<#EYLB?HVz<AhlUlW_Sg%BhfG
zN&j47V)+o-WciDN-V^Z}knW)zcDy?CRctrzh<w)?I_cL-5-3>L*lg_aGWQE9+k<q?
z<mZn~N``~mf%uP@q4^k9_xq+|-jPOV&uT7LKv7x@OhD2c|Kcc>H)Y{E!b8dw%54Mu
oE)Z?!y)!5{OMklTiZ@g%+5Y*EuM6D$)_zN`=lYHQhgqlp2QkCU2mk;8

delta 8261
zcmY*-2Ut@{)HVVFOA(`lj!0N3NkHixq$||`0Rl!5LhoHb%tAyVp=?B?D*-|Y2~rHw
z1q3!gKtu=_14>7#QrG@P_xqpc``>$>nRm*ax%bSubIx46^@f$GnMZ)0nkhoT8L@D1
z;xfl<Ci>^kWIMd5GJ2eYr|#rmJscbtYA$kcaNK1xX>Q(&mn=9$l-T9p;JygP{>f&S
zKpY%O$2pF}0X$U~hf>tU&-`J;X>?j)x7MD=++LF1<Ml<!07*@ri<%rL4lv2(_MO+?
zs)-`uQ;y0xFTpjM9FmD=B{?+N9W@32D0_X8g9FXxd2l0h`yU6Em)JLP@Yo;AwdXp^
zA#n^}b`1YFkNHKOA`ULAJL0FUiX;Q86cp|j1$fDGWS&*`QYgA`l8588JVz0V>nx94
zm6_ll74|}!uTP&gx7WM?_{-{~#D&xL;(S~jf8hVJS3mYQ*U8ss!A0q3`8e$OgKpm!
zQj-uzKXSfo<k<h{zGMH3Dd&C1*Jj0XjtczR&h{DVTo+LRMvCi&wd-j9N6p`ZW!KT1
z8x;r8C!6NyEabTnvR~wX2#Q)*@qT%{;Ux}APjLs83yzc0zg=?LfyYjL0@eQLSH*&W
zSf%uGn%j?H)iy^MKa!FX)oA*&NzBs45(#t{Eo)N!D;)@+#OB3#{eVBPI8o_U{bM6{
z=9o4*%rakDfDk7LjB&YIg9Ae=FIXbtQ#Vx8%)i%oXW8+Oidoqp0B#8HvUH@cSvjZ?
ziBd&~nu;|4ybI0ux(flE9EwY0m>++L@U;N({j?b&t(DUTj!(nG5q<%^E;&CCs#OqB
zEy|08wS16{ZGr-P|E{UmVjyq)oEUCe&duUHi{6z>RcZ%_D+3%So0hd(R0mQwnh>e9
z-jo<5zg+5I+9b4B)SK`f*R{-08p^Q5Q&MXhMg>JoSh)HLtdzL=m3EMc(>f8?h%|!|
zeDjr>hDq2q?;fn%IwBRcT$l=k#TL|OiaJr{VOE^w8=^o#*!>G$srf)r1hEUJfDRcz
z5Q;?wfiAr#i3~Lu%!+Ryl`aZ{S$G|~3c8ag%EbhWW+a*}4)_uCy+(;}c<j56$RZL6
zD2Pu5`i|C7k2A%Pf><%gxxA|?&>3LrNG%rH!R=0hn*ZsIMVw545is1!wr48LwG~v3
zEAut&9MumLpJN0d{vx3)XV7G?cDkK=g%$?L5SNn>6$g>i3Pj;}QA7+Xlap3=0Yt-z
znr(_9uyh1x1M<FV$0VwmC*+NWnpe{Zxp!T|d$d^F|F8iWhtUpkuh51TsH28*d2H2?
zJ_?rx_GSvocUp(OQG|q)hHyNB4EAKA2tqEmydp;HO$CqL47$;h;fGU4N{Ew(Xqi*{
zZ<<l#I86=EcFSF>5Ij!HZ$P}j$Tj3TtZX2c2V-9AO?PiZ+}}{sQvFiXR!qs&{wFC(
zVddUeR|GDQ-r*VTqhqVnIX;)pnw%J=QgLGg#pWK{RiR@%_S>R^Dp@ZxauvIqvude*
z9l<<<_Cj0eEDD+VPA$%kc@;-(Eg-)S2|>-xO)zKXy_>gMJ(;fV>miXkPg}Gr%x!&A
zC|!L+xwZ<S+nicj&9_DOR<(7spqW|rlnp9*02dsbq<z=!jn1=*PgG5+25OY*cknlB
zsFoU+R`5*M4LUhWY3suHh~~|hK%A2l4hQqMq3Fr?ni_rs)Goh6#=)Cs&1#4?OiEXr
z5XC-w(2gl3#S$9h4vCOxQ5~OnkL>xr%EilNo0e<N(;^hyDypq3<slZ1>}tgMcyGSZ
zY@tF2YHMq8!P+ZTLN6Bsk;P`4OivFfDGzCt92nu-=IFNQ-6bTI8aFsCkkQsN9~^{4
zcacZiioX<zYqn_XJS88#id$*1AKi2fLNNy#8+SA!*l*ZE)s>Qbo3D##3_|=&)qM7N
zGcYhT;+foXn6`Y9s6xv<KwM>nj;>^(dSE!BZ?vn8k{j8!qZSs9nwXgzA82c%61%!;
zF}BW%(LV=U1|Qyicu&A5t*fo8t7iuNJ_tI=)}URH$J-(%Bh2Kxd6G{cwUOG`NNL=B
zYQka-P<ldx5j{*Mvm<CkOMb=ZaJ6PWQnNcL$<ubFwFp1H8XT#ceXdBnQ(71{x8Oau
zoJdiz9~)flo>__2XR>mVCpKb@gtk-%{qm(0r@WaPkw(_t^MUc&x-zOQ=(~>cvr)}~
z3DWhvT<^RlcJ+%&x2sh52oqcJZKa{m*_~vWZ0q1O)>h&jR1U8vrHh!llT2hTMQV1O
z`6`yA5Uj0t?Sz;!ZNu^{vi=|QF6KIdHM0~J*M@&TIyo>m?>Cz3(>A<sFdZJ#V!$lf
zTFuL?o!d1KH<VEwoQ^g}I~UGu#kK*M+uk!<21Ow&y`-*0X`jm$M9U$+F!9zr_TJzr
zs*-ulYSo7_v+4CdwnHq=;3_5i+3Ae_$%TjxDVZF33VkKeZ#g-x>Tp`#hxfEkW8-_q
z4r(RdKxQQ$wEYjnGlkJ1I_<g9`c@`WCiKL!;Kt$ie%sy?i&1T8nQR>o9co8W>DV-3
zVkJ1P4N0DLSoJXv4yM*NVx03PCweDmwt@oKDm;;I+r?z;YIS%#swgdPqYx*02N|2e
zK`}-r?uH$P{S)np9~|mu>_eM_Tess3ar%In&BQZqKDJ%8tc^I;?I^>n9iGsv%-OB<
z9lxEhs0!_)JT#R`m}AAaw7bRY&W>(R3ygJijW5RQXz7<YFpi`mS&WbIty%R9`|N3_
zvFX!|xadaO5@K>Q0WiD~7}+tM9240wOIwI5W`+5y&2+#fc7pUo^<=V90kLhf^jW?k
z=0c)O_H;rUYb7vFyF=>4Gh6E9PJ9RQrF;@?h&6{^Pmn4w%o<y4Wv#~>357B86Seh<
zn9_7QYs@w5p2ti1G`2xW=p%#mi)1#L$oeBL+l=_6@k~jUN?2@Z_G(h^+$yR$Z_R2t
zx{Y-pb*+LxAN5-tUyEx+o;UBU9Ulr6lM3zcTa3_UZmL?n%}S#&ZiKYXT&;iMOJ8Q|
z|7w_3m>5V~w{$FDYDMqN32O{~`PN`~0pr~bmsd9&2;HUZiy5vE@vKB8UZ8zGTo?Hd
zo^3#L)=8APjrjK2ow9Mp+{eUi-e4+gJKo^_1ER9cRVSw|<sc$t;#WaZg(R8b^$_y+
zK)pHl1-2k8*r|G>`Kd!bFZA6}Auacoz{t=#?qUUpRc9iBtqA=>1Ej;WGj(z~rUR)u
z%~sjW7Mq=zw&e{bF2%IXX5WiW8Q%y@+RiVXTZ?<znSEsLN$^i1Os*p4*5WI=#&=>2
zSlLJN)WMx716J5G`L5rDUDHO}h|<xWd!3Yp$+Fe(Lg%er%UcT|Cc2T^P~^P(O|CX0
zOq`rwxtLtkZ>b!peE6#JRUo_VrLCzSsx{Ts))Y!#A-UQN`63|aqTIj_QO&Aweu9gf
z3pmZ4)CE2C!vf90l6-Lf3ow2;%V9YSrL>dfB4)U8640#1pOoft2v&n(Gd&JyQRq=$
z_qMpEp<um-(Z0j0(OM-V^{<g{kPb5r?S!{iYHizzPF3yIoK=0*$?rnJ5js;kaD6e)
zkzf*Qo7aFa_OWeNAR8F)M5T2ae?96+S>H}ofUfi?kEp_*o0s%O+E2#Rw_kCmcnb#2
zjHA3!rjs9H&kob(O$#R%{Y-rz6f{q~r#_#?*;GEZ{K`NWcuz|AoWEKpU$eAbO*4p&
zQ5PJ04Hm3XFB``E&gM<~Fd01l{DL$M1z*!b9+cR)pRZ{ikBh_=_5_$(UZ)e4GLHNw
zDv_59CGEoe!Gq^6!@Mms%5kTlm)Zd_<&BY|n(i#Vr{d#2Xj2|dc4pV8LDV1^oD#gG
zcD5fwI6f&Z&3EB!4~Ab*=|Wi_hEotK`2WfcHC|WVJ9JxmK+@|qR6!JAK%M1S>QnfW
zI!lkYrDFQ3@gCq;uGpvWa&{BR(g}=r!B!o=^otd*ku3GW_+NES*h)C47I+Cwy6%Sz
zE++z2_Fc+}Bm$7HhTR2l1$)Ocz1Gr(3HHRT!<=4O+JIqB*at3UgyTxerZ)UnoMzdn
zyGR+4d-&v))H@>V20o>b->Ic}9<5YNQhAmi5K#VWz)%o@oEmewtdx6A(chuJkb8&I
zFJ|;%JeW`6EV{I>uvEjeRNb!_KfGa?W}qMULO(8M^8=Iq;mzyjYX3-$DogrU;do)w
z*lT}>u|me!p3zakhntU$7M5-e9yKf_l)mJt*4Ido5j<1Db}zfFjIbno2Qjv{8%Jz5
zCJBvI2fvHxZcb?|rq<%lm1L3glNG&pN(Lwt#(Y2kIxs26Tr__&)@ZeLU|-@wgmE+@
zPM$dTIV(%ys1S8iDtK$fKjGc=(DJGHQgAVu5;Lk0*OAe_?i*or;)shoEY3)f>4xaq
zmn<;GoqObGCG5*$&{@}gg37n3z$W#lk3_BIv)<g444d$K-O|3T(mp53b-KM}5f!TG
z(xG}6{h}2|?3uZKElw&-yL}_0+w(?T*5G10?de*qs;XwmaN+BQ+Q=4GU}4dytAgmF
zea+ZZn2wv6%)^nE<#08v$Cj(8_MONsTJQM5ln4G(YZ@<J<Qr|5RMPty?=5@P1&;qF
z@aPzD%}U_WN6T@<lCo@yj7nHlc;^9AxZT_tz29}$-Y>tk+4iGmGitKY<Dr*TK=Q-8
zhYyAP^4|ebsyGsDl;S5O+&@&za|hKvN~%%26)HHRfgtvM@SgJ?@F<!hVm@hj_5Zk|
zn#r6g0;J%09t3zZzfm?3R1+z4XLfwiUlq_q8Vm+$m2WHJaZ#IZ_$F0*eVGMhgb<oG
zDjkAT_Zp?7pu8LR{~i*9aE_(11jRL}I%%&Eppmm#eC?`Dn^F=Dq1GSe3JJa@X|Zo~
zJtU=$TeVyo!g$pLd3B__$bXZQr41~3eM==^wa6XYBAI}zBi(Sh7+31MJPU}A!eP~U
zx{I@7UK@sJ;WT2RkN--RsIOrRpueV>b!6u;z#3jnay`tZF&@WbJk^(@eU7L5wU3@3
zSJS;3oVVuHu&ed_GOL+$eZkUgyB5)L4L;Sglt?l_j+`h93PD^~+lM>CDLw67d;v+q
z+~8H)#PwVMl+{U)waVAXZ%l4hH>JJ$y84lfaFhGSZ8*X|2N}jPK0m2)MX%)D`&jeA
z5+Ugchc5%LP(xZT@P}52*;A12W7@h`FO~JlnOg_cas84$!1F2&JCL#sjFG(z+`B5E
zW0HyBYV*y+e@m5f+|`vuSoO;WQvDoB%hA|$(xMstnyWpjm*;^Ot(o)XisT#e1@jEx
zLfYv0U1OcS=LAgFci}a((%krHij@#lM0JP=zvhBC>B_kG;UZ&3i)u1F;4h_N*0uW5
zeLnI;-eZd|!X@?bU!*Og>r#4}lYXRC)XpV@%dK>|7VmU%QxxBO-qNa`Hv&FeP%-9r
z983tUmfuy^o{L;GQT}}YcC2H%`+gR0)$My=x1r0dn>!N+r=BCeqEFp=aY%U{Eijmv
z=7=Dnc|RGx_<W}LbLEeFm1yU$vaEf_kBQafp!0gyoLUwm6Iwn<S*ZDIJ9c%lk2wmC
z&*%|_SH9*14mCZ`KV(ryM}B}wB&rF0sQAEy@^fw<4={n*fO9tCN@YJ2!K5_h2XHo=
z_RGc0NTh))Hb2mt^v89<rwW~FPvHUUpfopHI`#>BIX;WR|4whZz(;)L=6d|`2WbD*
z8f<@hx}_RO4hTr4-9co;n90&n#9=c$-yJjHeE<Lku!F()Y7l~XC+}FeNPy&E_sKYD
z^SQh$AQ12hfZM{9A|&<5`cv0;>Q9$vTW@e(6*krYG|l*@zfSd`yT|4SxD6aGr=$d_
z!hc!~A`#qeD5NrUJn+e#D9R117eSKsaF`8bh@W`e0}jTn!!2L{80T3F*vT5~FACuG
z7y$!dA{H*_r~%nS1D9U@bfO&aWTmBxSvg-S2oC_l0dOGjB#sP|A+?KGIHlh)1E4KT
z0RT7*0OY#^+TfI9!jb#{0KWx*+w^2Q7*q}#5wQ*+j@Mvoup|p7JkTY*T(tp;Z89l<
zJIR^J#eKwDf@)7%z(IKclmNJ_CjF#HO?r)`g#`!%%7c4znaKg!FV6M>W+@j9y+I&T
zD^t)fPkR<4Cog&l53)dKocfd9_y6MJ18xg8{&UG^f6+?-HcoL{SlGz@Dqm8f{r}Y#
z7Wk|`PAk|I{yO{q*A%>aE{3Qx4j(3wi09x+L?HhiOJ&omCU1kBu!q+zw+&ZI_8vv^
z*r5waDK4*1@mIST7mM|Lv7R)pSg^Qd)nYFw3pj90pZhl_{okC7>wr%(ikBy1ib}?;
zr8Y=)gp1(zcO$-u&D!=tcVFCS=*?9K+6s0%>~Os;7Mw;raAdk|orzz(&8_)-(FCxT
z2$z6<^<OaIJ9wNXZ9O0PtSh>9<5BxAd|?;8yO#z2lyI~2nrgl+@imOR5UP(<d7_aw
zxR|(nrs=Vsuve2Y#ymPHvpI87myVDeu4>;#cJ(1yNM>9yvO1{E!TIT2_Up1O12{8A
z{{*INa>e_-KG2fk)r%-BEKozy&zw*D=D2YC{!aIIpHrmzo&lHk0a)q?bifjt=8eYT
zeII~LBFDZ%WJCSJ$0xrFn3!)gdwmx=&+_ir*9b%4#tvQ$_pTbQLYT-^BWl8IAR@cj
zp5dpVKS3?MYYwuzYf*R#rHs1#cL|Ml4W<4k4W%w2c-LHRCQ7Y3wv{Cu0a+(lK_#Vx
zurnbSylHe|QUeLPp+DA4-iR$iBmF(v$27<rx-$L<Mr>W#nw#nP_*`>g|LQS6K@TAz
z9t`JZp+rsfo3FB-o9}to27j)7j5fwoHXBqMR2oz+sDM43ro1i)j2=HMf>ZZOR=r$Z
zqt5W^+ouH$PlQ7t5QjVj!lP;3(%1I|V)%8icFE{q-}>c0b7@I0&!sG(JM)Os2psR*
z+b_Ne2NJf<e`3An^2U7Qke&Dm;ZYf_t8d>LrHrygB{W#h7DRKRH9>LbW2X_m?7v>W
z1lsH{1=p@2V0ct0<Y{E=ThCwcxE}qU*%55X_+`Jp;Q!^vPyU*DFh;@tgqT2D>@>1(
zo_&ai>w0)FU;VOg(I#4YdX}-Z0c9?}bTzR-9TlbQaFIj_4b8TekJ_m3z|)qZ^rW*?
zLIZN%HqmmnfAvpK$A!bChOhKu6b=s$rjc=nzma<c)IjsBJt>+`;@>e!{d{Zb{4=K~
zmDewmy6=C!)jc04nm^M2Pr)999uRG4!@iIX|KJ7|+BfPy-88&*UJPgZ{)+*BIJhCq
zYGp-*E`9;o_8YZ(^LyTP5fvFrQNK)W?E;E6Y9rh9l|;C@22Y-ZyWxcgzvmxKyKSkX
z(WhS~ynOmLf$aa_!K;Lg9GT8no`Hr(Z%vQRc$?cxO#xieXC^wXQ0&)9LqNc{Wi%a2
zu|aT~#-L^>mn7unofQV)-d=Oy4Gq-}VBZ>Za)_;B+q75T_g$w?e^j&H6R~6jBWN&O
zmU)Wo9k_)bv-t{W16*gO7L0UJ^EyNqiyL3jK>E6-F;KlH5z<2TVdm!wMCRs$LT)(@
zPc|=`KzNtg`}}R^-REHW-ENoLT+1c`ALNN_YCVij+k5x-%$^CTlTv5wHvAdBTeoCv
zjfqyIeFpAAmg41hjYpi9Zef<(Xt(asp5PB(?MgnDJId)yUfhusHkeJy)^m#{Zs@UI
zcEWGBLbmX5;nwu5OYI`}R=ry5$W83G7d3lm>2ekX1X!4rbDpi<EH411jmv&AUT>Zi
zd&!9&i`2{ZDOgeWrb~}=ZYvh-L_KmIexjUP_TKuzfiUbEVf_rD|KLP)$*LaxnwxxK
zPe(~Pu703$B2cpd)hw@dqWbSjGd#bAGJtmlX+w<Z4*HZy9&~Sf>k&>NUvs~|SbF0|
z;6(oiEo@QkwbM5$x4<s<%PTtHyDWo$0S$J@gRxOQ7hjUS*j_4kmgv;P9?lZ}zr&vV
zZ}eE<SI2?H3dCY|)EbFXu<NS{Ma*o;KkZ%8pA84;PYs8}w*UUrTVIW_pA}L<W6Vs<
zN3`{BqLFD_&e{)O2R>zb4cyl_E}VUA_1oMJS^tyOJw{7tcrq$Ee(>LS<VG)OYxxLB
z*>@_I4wNuf?9Xz{h+|93lKDbJg;;-^?2XlAjBTP)d0qo;{A&0)e?!O9{jtwZJ$o}J
z-|y!qzhT^`o-nnaV}L&sO<8p-`EtJI2*}>hO6akL?h}P}RnDewoJy1HYw3iZS-a(%
zGQtsXUi8psD8Z5F31POANHt&ZF1N$Ze?lkxB{Cmvc+Z@cAPlNU4N2sLfvE#d&-Im@
z<FA{k0^QBppX|)SuRr1stOsTQ4ZtG<^EY2XAlGa<bA1Op*<V)9)q|H@)b}^f@OStm
z75c5}?O(Y>GJ4L6`}%0qn}1>=+;wm-p|eB5AI(T`Jq4LM>axz~%(=&Uvi;Q^T;~FR
zY(<9EmJe&S%?b)FR|-f9G1?;WoURoDfn3U@G5<gn6*V<I@%!FtQud;`fb@Q+RW<Q~
zOe@@n$99K~z8e*9jPCD1O5^@HoyBC%Cz3u{LZ)j-gfCNkU$sSwLx?`B?!zL1fsxYI
zXnx&-iaWizJVB(rp6U+%iuE~id0p1vknbbjzdjlY^QMA2rAi&SC^0u7Js%BERMd`V
zSJbqR*82^v>dhXzc8PSR4`?%RUt3D@6?q1ktKjRDvwr*BWBS*ftRL~;!?Ya-hmIbn
zeRa&jJ||+-TgtO+<d-e+Ks#Hg`;5~X1x<8ac-->y)NngQ`}_xBe_c1WZ*B4i6`E#3
zCt<TKoPWfc0*`<K&?b{wSHF&o)MpPSu%mv{pI-nOvio^OBM{>EpgQ%it8Ef$h`8yU
zyfz@t{2aYxkz0uF8fG;J(F6!N0GS19h_A#8LAUoZFh`T4{Np(C+!BOGFdvSNZHSdP
z+sk#72NN=T+Ju;qFvLj@z{e+@cRL=(j`hC(zK{-+WQCqj{<;;yP-^a}pe+O;{~~+3
z282a5c`aIj2?XCv?3YdHt#fcKUCB8Yr)`Y{g0xMbMSQEx5S7*|>3N;Y{DX<UvDniN
zCK9)2qOnKL*6C}uUAX<H^#VG!TZ7K}79L1nGU3|pIv)t^+wU%yFjA}>-|yN_J|ccd
zc1d(gbhPQlFC{u&{$Qj?yp})sHN5k`LM6MUOL3lE@Mk9On9oKS-NR5>yKQfb6klFc
z>g|gtW+@7ot=Ov4Vb&P#n=<P=24mBS9b>1q_d0jn66SW$ZV6-c(S>8tqJGQK_w@L{
zHHfnD0Kcr`1x?@wM$MB(71Yg~8yj6Q*7XMy#B1jbF6RS_{$Sj%yUcP3td4iCzN{a!
zwDdty(O_t3BKnV^p^S`-+*4HPdTXhy6)*`SS^=w>m>lpe6oEUN+lrZzI|>FGR@^SE
z{NDBSveE3l?>}niUQkeXn%n!!n27S1m(yI+Nh#jT1yGS*AU`RI?@&Swc_^BEZ^lY+
zrXl%W=?=N$9rEB|(pFt}ZB{ZYFe?or`C1|n5}K6-cGZB5uc<u})k%Oz_VrD?NCx-y
z`Qilr!|}VIzP`SHP7C-tC(K9apv9ydZ&&dwhv`UVP8BQOnAp@Gy@$Z)AQ&5N?1+UO
zyy%ZPKUmgta8ss!v1#9=)29k?A~8n3QRFYe?uo-+gb^Fi9QT79_b>ckgdEG7q=gZC
zXPOKh6w^Q5EH`|q-jTOr_}c|1<fzLn+<c^%v9NwqU*?2@s`X`5uok*#)DH?PqQyYU
z*7QBbaw|XRJLFczl(NIns8-l;Jv^$UUi6@_SwllZ&3j+}G4f;g;-(3|@yJ`GKfzb_
z4nia(KsnXq82eCyr6J<_>UheZ`vPsXfv%{A8Isp{Q$z9;G!t(x=qVO}`nq1lN0v-O
zVr&MCM>dSLganMypRG5QyhfvUYJw1{JLWcdep8QcSXFm6`ZJxmUDIANg1i?(bqGDi
zVn_IRW)OS&33;C_$d<hAmg)-P{pTlt9nFS|Ij4Oz;`_4xuW+*OW6o}__qTX3f$=@J
zc&dAnby{!n{|o$R%r$=tynag|zwegf{}ps!!rb0e&@ELRpEP=}A3MmN8Nfob%eO0P
z)?F)VKE#6b^Cn|utmW8=4SN}!A1p+*0VE{UHSc=_JawuGe`aGF<Z|wc$CWDt$emKX
z^~PJ%qv2Kp>;*9B|77SW%J=zR#sZrF{?mPgVRgpB&30s%A{Woy+{=~d32*OgF#|lt
z`Im!r*P9}=UEO9oFgHE68>?0xjw&h%vj)>L*~?yt6Prw;5hKZajY+r%YF_)QX;2q0
z{<f6Z8Tr{HsB6-b*fZ-RMiV=csHB>lh@(rBN&YQQLa)8qGJ0y=!sds7Gq2ycwH&<K
zI>R)g{(S()4IiRkh8$hL!uHYLpS^AW*F*(QY`{7jhvyC-B8y=)lSrJHzP>;^uX)*8
zjK0u1V@)?3LqPJHjz_T{u$HzQsVD8A0BH(z&Ap7_f#@+3tApvw)`crB2@lQbPt@g|
z>3a}odP0ucijiB!lHt{w4Y!P$O+)zjuAw6^?6*a-6KtowgPUhw-}24pY)|nS+F_0$
zGwK>$RkhfE11>-Ppy3tstU@cMC@VCKZKQ5C?QAn$*0WZoboI}NAut*m>MPLe_4P3(
irR<aUb(P!w2p7S#Wp4itA4mUCUBX5X<ObHI@Ba_8sk}u1

diff --git a/man/layout_as_metromap.Rd b/man/layout_as_metromap.Rd
new file mode 100644
index 0000000..5e89e46
--- /dev/null
+++ b/man/layout_as_metromap.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/metro_multicriteria.R
+\name{layout_as_metromap}
+\alias{layout_as_metromap}
+\title{Metro Map Layout}
+\usage{
+layout_as_metromap(object, xy, l = 2, gr = 0.0025, w = rep(1, 5), bsize = 5)
+}
+\arguments{
+\item{object}{original graph}
+
+\item{xy}{initial layout of the original graph}
+
+\item{l}{desired multiple of grid point spacing. (l*gr determines desired edge length)}
+
+\item{gr}{grid spacing. (l*gr determines desired edge length)}
+
+\item{w}{weight vector for criteria (see details)}
+
+\item{bsize}{number of grid points a station can move away rom its original position}
+}
+\value{
+new coordinates for stations
+}
+\description{
+Metro map layout based on multicriteria optimization
+}
+\details{
+The function optimizes the following five criteria using a hill climbing algorithm:
+\itemize{
+\item \emph{Angular Resolution Criterion}: The angles of incident edges at each station should be maximized, because if there is only a small angle between any two adjacent edges, then it can become difficult to distinguish between them
+\item \emph{Edge Length Criterion}: The edge lengths across the whole map should be approximately equal to ensure regular spacing between stations. It is based on the preferred multiple, l, of the grid spacing, g. The purpose of the criterion is to penalize edges that are longer than or shorter than lg.
+\item \emph{Balanced Edge Length Criterion}: The length of edges incident to a particular station should be similar
+\item \emph{Line Straightness Criterion}: (not yet implemented) Edges that form part of a line should, where possible, be co-linear either side of each station that the line passes through
+\item \emph{Octiinearity Criterion}: Each edge should be drawn horizontally, vertically, or diagonally at 45 degree, so we penalize edges that are not at a desired angle
+}
+}
+\examples{
+# the algorithm has problems with parallel edges
+library(igraph)
+g <- simplify(metro_berlin)
+xy <- cbind(V(g)$lon, V(g)$lat) * 100
+
+# the algorithm is not very stable. try playing with the parameters
+xy_new <- layout_as_metromap(g, xy, l = 2, gr = 0.5, w = c(100, 100, 1, 1, 100), bsize = 35)
+}
+\references{
+Stott, Jonathan, et al. "Automatic metro map layout using multicriteria optimization." IEEE Transactions on Visualization and Computer Graphics 17.1 (2010): 101-114.
+}
+\author{
+David Schoch
+}
diff --git a/man/metro_berlin.Rd b/man/metro_berlin.Rd
new file mode 100644
index 0000000..86e8fde
--- /dev/null
+++ b/man/metro_berlin.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data-examples.R
+\docType{data}
+\name{metro_berlin}
+\alias{metro_berlin}
+\title{Subway network of Berlin}
+\format{
+igraph object
+}
+\usage{
+metro_berlin
+}
+\description{
+A dataset containing the subway network of Berlin
+}
+\references{
+Kujala, Rainer, et al. "A collection of public transport network data sets for 25 cities." Scientific data 5 (2018): 180089.
+}
+\keyword{datasets}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index f6e912a..8193d3d 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -69,6 +69,84 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// criterion_angular_resolution
+double criterion_angular_resolution(List adj, NumericMatrix xy);
+RcppExport SEXP _graphlayouts_criterion_angular_resolution(SEXP adjSEXP, SEXP xySEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< List >::type adj(adjSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type xy(xySEXP);
+    rcpp_result_gen = Rcpp::wrap(criterion_angular_resolution(adj, xy));
+    return rcpp_result_gen;
+END_RCPP
+}
+// criterion_edge_length
+double criterion_edge_length(IntegerMatrix el, NumericMatrix xy, double lg);
+RcppExport SEXP _graphlayouts_criterion_edge_length(SEXP elSEXP, SEXP xySEXP, SEXP lgSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerMatrix >::type el(elSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type xy(xySEXP);
+    Rcpp::traits::input_parameter< double >::type lg(lgSEXP);
+    rcpp_result_gen = Rcpp::wrap(criterion_edge_length(el, xy, lg));
+    return rcpp_result_gen;
+END_RCPP
+}
+// criterion_balanced_edge_length
+double criterion_balanced_edge_length(List adj_deg2, NumericMatrix xy);
+RcppExport SEXP _graphlayouts_criterion_balanced_edge_length(SEXP adj_deg2SEXP, SEXP xySEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< List >::type adj_deg2(adj_deg2SEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type xy(xySEXP);
+    rcpp_result_gen = Rcpp::wrap(criterion_balanced_edge_length(adj_deg2, xy));
+    return rcpp_result_gen;
+END_RCPP
+}
+// criterion_line_straightness
+double criterion_line_straightness();
+RcppExport SEXP _graphlayouts_criterion_line_straightness() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(criterion_line_straightness());
+    return rcpp_result_gen;
+END_RCPP
+}
+// criterion_octilinearity
+double criterion_octilinearity(IntegerMatrix el, NumericMatrix xy);
+RcppExport SEXP _graphlayouts_criterion_octilinearity(SEXP elSEXP, SEXP xySEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< IntegerMatrix >::type el(elSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type xy(xySEXP);
+    rcpp_result_gen = Rcpp::wrap(criterion_octilinearity(el, xy));
+    return rcpp_result_gen;
+END_RCPP
+}
+// layout_as_metro_iter
+NumericMatrix layout_as_metro_iter(List adj, IntegerMatrix el, List adj_deg2, NumericMatrix xy, NumericMatrix bbox, double l, double gr, NumericVector w, double bsize);
+RcppExport SEXP _graphlayouts_layout_as_metro_iter(SEXP adjSEXP, SEXP elSEXP, SEXP adj_deg2SEXP, SEXP xySEXP, SEXP bboxSEXP, SEXP lSEXP, SEXP grSEXP, SEXP wSEXP, SEXP bsizeSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< List >::type adj(adjSEXP);
+    Rcpp::traits::input_parameter< IntegerMatrix >::type el(elSEXP);
+    Rcpp::traits::input_parameter< List >::type adj_deg2(adj_deg2SEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type xy(xySEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type bbox(bboxSEXP);
+    Rcpp::traits::input_parameter< double >::type l(lSEXP);
+    Rcpp::traits::input_parameter< double >::type gr(grSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type w(wSEXP);
+    Rcpp::traits::input_parameter< double >::type bsize(bsizeSEXP);
+    rcpp_result_gen = Rcpp::wrap(layout_as_metro_iter(adj, el, adj_deg2, xy, bbox, l, gr, w, bsize));
+    return rcpp_result_gen;
+END_RCPP
+}
 // reweighting
 NumericVector reweighting(IntegerMatrix el, List N_ranks);
 RcppExport SEXP _graphlayouts_reweighting(SEXP elSEXP, SEXP N_ranksSEXP) {
@@ -191,6 +269,12 @@ static const R_CallMethodDef CallEntries[] = {
     {"_graphlayouts_constrained_stress_major", (DL_FUNC) &_graphlayouts_constrained_stress_major, 6},
     {"_graphlayouts_constrained_stress3D", (DL_FUNC) &_graphlayouts_constrained_stress3D, 3},
     {"_graphlayouts_constrained_stress_major3D", (DL_FUNC) &_graphlayouts_constrained_stress_major3D, 6},
+    {"_graphlayouts_criterion_angular_resolution", (DL_FUNC) &_graphlayouts_criterion_angular_resolution, 2},
+    {"_graphlayouts_criterion_edge_length", (DL_FUNC) &_graphlayouts_criterion_edge_length, 3},
+    {"_graphlayouts_criterion_balanced_edge_length", (DL_FUNC) &_graphlayouts_criterion_balanced_edge_length, 2},
+    {"_graphlayouts_criterion_line_straightness", (DL_FUNC) &_graphlayouts_criterion_line_straightness, 0},
+    {"_graphlayouts_criterion_octilinearity", (DL_FUNC) &_graphlayouts_criterion_octilinearity, 2},
+    {"_graphlayouts_layout_as_metro_iter", (DL_FUNC) &_graphlayouts_layout_as_metro_iter, 9},
     {"_graphlayouts_reweighting", (DL_FUNC) &_graphlayouts_reweighting, 2},
     {"_graphlayouts_sparseStress", (DL_FUNC) &_graphlayouts_sparseStress, 6},
     {"_graphlayouts_stress", (DL_FUNC) &_graphlayouts_stress, 3},
diff --git a/src/metroLayout.cpp b/src/metroLayout.cpp
new file mode 100644
index 0000000..6814acf
--- /dev/null
+++ b/src/metroLayout.cpp
@@ -0,0 +1,185 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+double angle_between_edges(NumericVector pvec,NumericVector qvec) {
+  if(pvec[0]==qvec[0] && pvec[1]==qvec[1]){
+    return 0.0;
+  }
+  double dot_pq = pvec[0] * qvec[0] + pvec[1] * qvec[1];
+  double mag_pq = sqrt(pvec[0]*pvec[0]+pvec[1]*pvec[1]) * sqrt(qvec[0]*qvec[0]+qvec[1]*qvec[1]);
+
+  if(dot_pq/mag_pq < -0.99){
+    return M_PI;
+  }
+  if(dot_pq/mag_pq > 0.99){
+    return 0.0;
+  }
+
+  return acos(dot_pq/mag_pq);
+}
+
+// [[Rcpp::export]]
+double criterion_angular_resolution(List adj,NumericMatrix xy){
+  int n = adj.length();
+  double crit=0;
+  double elen;
+  double angle;
+
+  for(int i=0;i<n;++i){
+    IntegerVector Ni = adj[i];
+    if(Ni.length()!=1){
+      NumericMatrix edges_xy(Ni.length(),4);
+      for(int j=0;j<Ni.length();++j){
+        edges_xy(j,0) = xy(i,0);
+        edges_xy(j,1) = xy(i,1);
+        edges_xy(j,2) = xy(Ni[j],0);
+        edges_xy(j,3) = xy(Ni[j],1);
+
+        edges_xy(j,2) -= edges_xy(j,0);
+        edges_xy(j,3) -= edges_xy(j,1);
+
+        elen = sqrt(edges_xy(j,2) * edges_xy(j,2) + edges_xy(j,3) * edges_xy(j,3));
+        edges_xy(j,2) /= elen;
+        edges_xy(j,3) /= elen;
+      }
+      for(int p=0;p<(Ni.length()-1);++p){
+        for(int q=p+1;q<Ni.length();++q){
+          NumericVector pvec = {edges_xy(p,2),edges_xy(p,3)};
+          NumericVector qvec = {edges_xy(q,2),edges_xy(q,3)};
+          angle = angle_between_edges(pvec,qvec);
+
+          crit += std::abs(2.0 * M_PI / Ni.length() - angle);
+        }
+      }
+    }
+  }
+  return crit;
+}
+
+// [[Rcpp::export]]
+double criterion_edge_length(IntegerMatrix el, NumericMatrix xy, double lg){
+  double crit = 0;
+  double elen;
+  NumericVector evec(4);
+  for(int e=0;e<el.nrow();++e){
+    evec = {xy(el(e,0),0),xy(el(e,0),1),
+            xy(el(e,1),0),xy(el(e,1),1)};
+    elen = sqrt((evec[0] - evec[2]) * (evec[0] - evec[2]) +
+                (evec[1] - evec[3]) * (evec[1] - evec[3]));
+    crit += std::abs(elen/lg - 1.0);
+  }
+  return crit;
+}
+
+// [[Rcpp::export]]
+double criterion_balanced_edge_length(List adj_deg2,NumericMatrix xy){
+  int n2 = adj_deg2.length();
+  double crit = 0.0;
+
+  if(n2==0){
+    return 0.0;
+  }
+  NumericMatrix edges_xy(2,4);
+  NumericMatrix elen(2);
+  for(int i=0;i<n2;++i){
+    IntegerVector Ni = adj_deg2[i];
+    edges_xy(0,0) = xy(i,0);
+    edges_xy(0,1) = xy(i,1);
+    edges_xy(0,2) = xy(Ni[0],0);
+    edges_xy(0,3) = xy(Ni[0],1);
+
+    edges_xy(0,2) -= edges_xy(0,0);
+    edges_xy(0,3) -= edges_xy(0,1);
+
+    elen[0] = sqrt(edges_xy(0,2) * edges_xy(0,2) + edges_xy(0,3) * edges_xy(0,3));
+
+
+    edges_xy(1,0) = xy(i,0);
+    edges_xy(1,1) = xy(i,1);
+    edges_xy(1,2) = xy(Ni[1],0);
+    edges_xy(1,3) = xy(Ni[1],1);
+
+    edges_xy(1,2) -= edges_xy(1,0);
+    edges_xy(1,3) -= edges_xy(1,1);
+
+    elen[1] = sqrt(edges_xy(1,2) * edges_xy(1,2) + edges_xy(1,3) * edges_xy(1,3));
+
+    crit += std::abs(elen[1]-elen[0]);
+  }
+  return crit;
+}
+
+// [[Rcpp::export]]
+double criterion_line_straightness(){
+  return 0.0;
+}
+
+// [[Rcpp::export]]
+double criterion_octilinearity(IntegerMatrix el,NumericMatrix xy){
+  double crit = 0;
+  NumericVector edge(4);
+  for(int e=0;e<el.nrow();++e){
+    edge = {xy(el(e,0),0),xy(el(e,0),1),xy(el(e,1),0),xy(el(e,1),1)};
+    crit += std::abs(sin(4.0 * atan(std::abs(edge(1)-edge(3)) / std::abs(edge(0)-edge(2)))));
+  }
+  return crit;
+}
+
+
+double criterion_sum(List adj, IntegerMatrix el, List adj_deg2, NumericMatrix xy,
+                     double lg, NumericVector w){
+  double cn1 = criterion_angular_resolution(adj,xy);
+  double cn2 = criterion_edge_length(el,xy,lg);
+  double cn3 = criterion_balanced_edge_length(adj_deg2,xy);
+  double cn4 = criterion_line_straightness();
+  double cn5 = criterion_octilinearity(el,xy);
+
+  return w[0] * cn1 + w[1] * cn2 + w[2] * cn3 + w[3] * cn4 + w[4] * cn5;
+}
+
+// [[Rcpp::export]]
+NumericMatrix layout_as_metro_iter(List adj, IntegerMatrix el, List adj_deg2,
+                                   NumericMatrix xy, NumericMatrix bbox,
+                                   double l, double gr, NumericVector w, double bsize){
+  double cur_crit, tmp_crit;
+  double lg = l * gr;
+  int n = adj.length();
+  double x, y, xold, yold, xbest, ybest;
+  bool running = true;
+  int k = 0;
+  NumericVector xoptions(8),yoptions(8);
+
+  cur_crit = criterion_sum(adj,el,adj_deg2,xy,lg,w);
+
+  while(running){
+    running = false;
+    for(int v=0;v<n;++v){
+      xold = xy(v,0);
+      yold = xy(v,1);
+      xbest = xold;
+      ybest = yold;
+      xoptions = {xold-gr,xold,xold+gr,xold-gr,xold+gr,xold-gr,xold,xold+gr};
+      yoptions = {yold+gr,yold+gr,yold+gr,yold,yold,yold-gr,yold-gr,yold-gr};
+      for(int i=0;i<xoptions.length();++i){
+        x = xoptions[i];
+        y = yoptions[i];
+        if((x>=bbox(v,0)) && (x<=bbox(v,2)) && (y>=bbox(v,1)) && (y<=bbox(v,3))){
+          xy(v,0) = x;
+          xy(v,1) = y;
+          tmp_crit = criterion_sum(adj,el,adj_deg2,xy,lg,w);
+          if(tmp_crit < cur_crit){
+            cur_crit = tmp_crit;
+            running = true;
+            xbest = xy(v,0);
+            ybest = xy(v,1);
+          }
+        }
+      }
+      xy(v,0) = xbest;
+      xy(v,1) = ybest;
+    }
+    // Rcout << "run: " << k << " min: " << cur_crit << std::endl;
+    k +=1;
+  }
+  return xy;
+}
diff --git a/src/stress.cpp b/src/stress.cpp
index e88c685..83d05bc 100644
--- a/src/stress.cpp
+++ b/src/stress.cpp
@@ -2,129 +2,119 @@
 using namespace Rcpp;
 
 // [[Rcpp::export]]
-double stress(NumericMatrix x, NumericMatrix W, NumericMatrix D){
-  double fct=0;
-  int n=x.nrow();
-  for(int i=0;i<(n-1);++i){
-    for(int j=(i+1);j<n;++j){
-      double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-      fct+=W(i,j)*(denom-D(i,j))*(denom-D(i,j));
+double stress(NumericMatrix x, NumericMatrix W, NumericMatrix D) {
+  double fct = 0;
+  int n = x.nrow();
+  for (int i = 0; i < (n - 1); ++i) {
+    for (int j = (i + 1); j < n; ++j) {
+      double denom = sqrt((x(i, 0) - x(j, 0)) * (x(i, 0) - x(j, 0)) +
+                          (x(i, 1) - x(j, 1)) * (x(i, 1) - x(j, 1)));
+      fct += W(i, j) * (denom - D(i, j)) * (denom - D(i, j));
     }
   }
   return fct;
 }
 
 // [[Rcpp::export]]
-NumericMatrix stress_major(NumericMatrix y,
-                     NumericMatrix W,
-                     NumericMatrix D,
-                     int iter,
-                     double tol) {
+NumericMatrix stress_major(NumericMatrix y, NumericMatrix W, NumericMatrix D,
+                           int iter, double tol) {
   int n = y.nrow();
 
-  NumericMatrix x(n,2);
-  for(int i=0;i<n;++i){
-    for(int d=0;d<2;++d){
-      x(i,d)=y(i,d);
-    }
-  }
-
-  NumericVector wsum(n);
-  for(int i=0;i<n;++i){
-    for(int j=0;j<n;++j){
-      wsum[i]+=W(i,j);
-    }
-  }
-
-  double stress_old = stress(x,W,D);
-
-  for(int k=0; k<iter; ++k){
-    NumericMatrix xnew(n,2); //out or in?
-    for(int i=0;i<n;++i){
-      for(int j=0; j<n;++j){
-        if(i!=j){
-          double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-          if(denom>0.00001){
-            xnew(i,0) += W(i,j)*(x(j,0)+D(i,j)*(x(i,0)-x(j,0))/denom);
-            xnew(i,1) += W(i,j)*(x(j,1)+D(i,j)*(x(i,1)-x(j,1))/denom);
+  NumericMatrix x(clone(y));
+
+  NumericVector wsum = rowSums(W);
+
+  double stress_old = stress(x, W, D);
+  NumericMatrix xnew(n, 2);  // out or in?
+  for (int k = 0; k < iter; ++k) {
+    std::fill(xnew.begin(), xnew.end(), 0);
+    for (int i = 0; i < n; ++i) {
+      for (int j = 0; j < n; ++j) {
+        if (i != j) {
+          double denom = sqrt((x(i, 0) - x(j, 0)) * (x(i, 0) - x(j, 0)) +
+                              (x(i, 1) - x(j, 1)) * (x(i, 1) - x(j, 1)));
+          if (denom > 0.00001) {
+            xnew(i, 0) +=
+                W(i, j) * (x(j, 0) + D(i, j) * (x(i, 0) - x(j, 0)) / denom);
+            xnew(i, 1) +=
+                W(i, j) * (x(j, 1) + D(i, j) * (x(i, 1) - x(j, 1)) / denom);
           }
         }
       }
-      xnew(i,0) = xnew(i,0)/wsum[i];
-      xnew(i,1) = xnew(i,1)/wsum[i];
+      xnew(i, 0) = xnew(i, 0) / wsum[i];
+      xnew(i, 1) = xnew(i, 1) / wsum[i];
     }
-    double stress_new = stress(xnew,W,D);
-    double eps = (stress_old-stress_new)/stress_old;
+    double stress_new = stress(xnew, W, D);
+    double eps = (stress_old - stress_new) / stress_old;
 
-    if(eps<= tol){
+    if (eps <= tol) {
       break;
     }
-    stress_old=stress_new;
-    for(int i=0;i<n;++i){
-      x(i,0)=xnew(i,0);
-      x(i,1)=xnew(i,1);
-    }
+    stress_old = stress_new;
+    x = clone(xnew);
   }
   return x;
 }
 
 // [[Rcpp::export]]
-NumericMatrix stress_radii(NumericMatrix y,
-                           NumericMatrix W,
-                           NumericMatrix D,
-                           NumericVector r,
-                           NumericVector tseq) {
+NumericMatrix stress_radii(NumericMatrix y, NumericMatrix W, NumericMatrix D,
+                           NumericVector r, NumericVector tseq) {
   int n = y.nrow();
   int m = tseq.length();
 
-  NumericMatrix x(n,2);
-  for(int i=0;i<n;++i){
-    for(int d=0;d<2;++d){
-      x(i,d)=y(i,d);
+  NumericMatrix x(n, 2);
+  for (int i = 0; i < n; ++i) {
+    for (int d = 0; d < 2; ++d) {
+      x(i, d) = y(i, d);
     }
   }
 
   NumericVector wsum(n);
-  for(int i=0;i<n;++i){
-    for(int j=0;j<n;++j){
-      wsum[i]+=W(i,j);
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      wsum[i] += W(i, j);
     }
   }
 
   NumericVector rpow(n);
-  for(int i=0;i<n;++i){
-    rpow[i] = 1/(r[i] * r[i]);
+  for (int i = 0; i < n; ++i) {
+    rpow[i] = 1 / (r[i] * r[i]);
   }
 
-  for(int s=0; s<m; ++s){
+  for (int s = 0; s < m; ++s) {
     double t = tseq[s];
-    for(int k=0; k<1; ++k){
-      NumericMatrix xnew(n,2); //out or in?
-      for(int i=0;i<n;++i){
-        for(int j=0; j<n;++j){
-          if(i!=j){
-            double bij = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-            double ai = sqrt(x(i,0)*x(i,0)+x(i,1)*x(i,1));
-            if(ai<0.0001){
-              ai=0;
-            } else{
-              ai=1/ai;
+    for (int k = 0; k < 1; ++k) {
+      NumericMatrix xnew(n, 2);  // out or in?
+      for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < n; ++j) {
+          if (i != j) {
+            double bij = sqrt((x(i, 0) - x(j, 0)) * (x(i, 0) - x(j, 0)) +
+                              (x(i, 1) - x(j, 1)) * (x(i, 1) - x(j, 1)));
+            double ai = sqrt(x(i, 0) * x(i, 0) + x(i, 1) * x(i, 1));
+            if (ai < 0.0001) {
+              ai = 0;
+            } else {
+              ai = 1 / ai;
             }
-            if(bij<0.0001){
-              bij=0;
-            } else{
-              bij=1/bij;
+            if (bij < 0.0001) {
+              bij = 0;
+            } else {
+              bij = 1 / bij;
             }
-            xnew(i,0) += (1-t) * W(i,j)*(x(j,0)+D(i,j)*(x(i,0)-x(j,0))*bij)+t*rpow[i]*(r[i]*x(i,0)*ai);
-            xnew(i,1) += (1-t) * W(i,j)*(x(j,1)+D(i,j)*(x(i,1)-x(j,1))*bij)+t*rpow[i]*(r[i]*x(i,1)*ai);
+            xnew(i, 0) += (1 - t) * W(i, j) *
+                              (x(j, 0) + D(i, j) * (x(i, 0) - x(j, 0)) * bij) +
+                          t * rpow[i] * (r[i] * x(i, 0) * ai);
+            xnew(i, 1) += (1 - t) * W(i, j) *
+                              (x(j, 1) + D(i, j) * (x(i, 1) - x(j, 1)) * bij) +
+                          t * rpow[i] * (r[i] * x(i, 1) * ai);
           }
         }
-        xnew(i,0) = xnew(i,0)/((1-t) * wsum[i] + t * rpow[i]);
-        xnew(i,1) = xnew(i,1)/((1-t) * wsum[i] + t * rpow[i]);
+        xnew(i, 0) = xnew(i, 0) / ((1 - t) * wsum[i] + t * rpow[i]);
+        xnew(i, 1) = xnew(i, 1) / ((1 - t) * wsum[i] + t * rpow[i]);
       }
-      for(int i=0;i<n;++i){
-        x(i,0)=xnew(i,0);
-        x(i,1)=xnew(i,1);
+      for (int i = 0; i < n; ++i) {
+        x(i, 0) = xnew(i, 0);
+        x(i, 1) = xnew(i, 1);
       }
     }
   }
@@ -132,69 +122,68 @@ NumericMatrix stress_radii(NumericMatrix y,
 }
 
 // [[Rcpp::export]]
-NumericMatrix stress_focus(NumericMatrix y,
-                           NumericMatrix W,
-                           NumericMatrix D,
-                           NumericMatrix Z,
-                           NumericVector tseq,
-                           int iter,
+NumericMatrix stress_focus(NumericMatrix y, NumericMatrix W, NumericMatrix D,
+                           NumericMatrix Z, NumericVector tseq, int iter,
                            double tol) {
   int n = y.nrow();
   int m = tseq.length();
 
-  NumericMatrix x(n,2);
-  for(int i=0;i<n;++i){
-    for(int d=0;d<2;++d){
-      x(i,d)=y(i,d);
+  NumericMatrix x(n, 2);
+  for (int i = 0; i < n; ++i) {
+    for (int d = 0; d < 2; ++d) {
+      x(i, d) = y(i, d);
     }
   }
 
   NumericVector wsum(n);
   NumericVector zsum(n);
-  for(int i=0;i<n;++i){
-    for(int j=0;j<n;++j){
-      wsum[i]+=W(i,j);
-      zsum[i]+=Z(i,j);
+  for (int i = 0; i < n; ++i) {
+    for (int j = 0; j < n; ++j) {
+      wsum[i] += W(i, j);
+      zsum[i] += Z(i, j);
     }
   }
 
-  double stress_oldW = stress(x,W,D);
+  double stress_oldW = stress(x, W, D);
   // double stress_oldZ = stress(x,Z,D);
   double stress_old = stress_oldW;
-  for(int s=0; s<m;++s){
+  for (int s = 0; s < m; ++s) {
     double t = tseq[s];
 
-    for(int k=0; k<iter; ++k){
-      NumericMatrix xnew(n,2); //out or in?
-      for(int i=0;i<n;++i){
-        for(int j=0; j<n;++j){
-          if(i!=j){
-            double denom = sqrt((x(i,0)-x(j,0))*(x(i,0)-x(j,0))+(x(i,1)-x(j,1))*(x(i,1)-x(j,1)));
-            if(denom>0.00001){
-              denom = 1/denom;
-            } else{
+    for (int k = 0; k < iter; ++k) {
+      NumericMatrix xnew(n, 2);  // out or in?
+      for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < n; ++j) {
+          if (i != j) {
+            double denom = sqrt((x(i, 0) - x(j, 0)) * (x(i, 0) - x(j, 0)) +
+                                (x(i, 1) - x(j, 1)) * (x(i, 1) - x(j, 1)));
+            if (denom > 0.00001) {
+              denom = 1 / denom;
+            } else {
               denom = 0;
             }
-            xnew(i,0) += ((1-t)*W(i,j)+t*Z(i,j))*(x(j,0)+D(i,j)*(x(i,0)-x(j,0))*denom);
-            xnew(i,1) += ((1-t)*W(i,j)+t*Z(i,j))*(x(j,1)+D(i,j)*(x(i,1)-x(j,1))*denom);
+            xnew(i, 0) += ((1 - t) * W(i, j) + t * Z(i, j)) *
+                          (x(j, 0) + D(i, j) * (x(i, 0) - x(j, 0)) * denom);
+            xnew(i, 1) += ((1 - t) * W(i, j) + t * Z(i, j)) *
+                          (x(j, 1) + D(i, j) * (x(i, 1) - x(j, 1)) * denom);
           }
         }
-        xnew(i,0) = xnew(i,0)/((1-t)*wsum[i] + t*zsum[i]);
-        xnew(i,1) = xnew(i,1)/((1-t)*wsum[i] + t*zsum[i]);
+        xnew(i, 0) = xnew(i, 0) / ((1 - t) * wsum[i] + t * zsum[i]);
+        xnew(i, 1) = xnew(i, 1) / ((1 - t) * wsum[i] + t * zsum[i]);
       }
-      double stress_newW = stress(xnew,W,D);
-      double stress_newZ = stress(xnew,Z,D);
-      double stress_new = (1-t)*stress_newW + t*stress_newZ;
+      double stress_newW = stress(xnew, W, D);
+      double stress_newZ = stress(xnew, Z, D);
+      double stress_new = (1 - t) * stress_newW + t * stress_newZ;
 
-      for(int i=0;i<n;++i){
-        x(i,0)=xnew(i,0);
-        x(i,1)=xnew(i,1);
+      for (int i = 0; i < n; ++i) {
+        x(i, 0) = xnew(i, 0);
+        x(i, 1) = xnew(i, 1);
       }
-      double eps=(stress_old-stress_new)/stress_old;
-      if(eps<= tol){
+      double eps = (stress_old - stress_new) / stress_old;
+      if (eps <= tol) {
         break;
       }
-      stress_old=stress_new;
+      stress_old = stress_new;
     }
   }
   return x;