From 0584e705bff5f8a68c601081e589e18b1353e75a Mon Sep 17 00:00:00 2001
From: Lalit Rajendran <lrajendr@kolmogorov.ecn.purdue.edu>
Date: Mon, 3 May 2021 15:29:13 -0400
Subject: [PATCH] moved common functions to module helper_functions

---
 __pycache__/helper_functions.cpython-38.pyc | Bin 14057 -> 10317 bytes
 helper_functions.py                         | 446 ++++++++++++++++++++
 sample_script.py                            | 395 +----------------
 3 files changed, 453 insertions(+), 388 deletions(-)
 create mode 100755 helper_functions.py

diff --git a/__pycache__/helper_functions.cpython-38.pyc b/__pycache__/helper_functions.cpython-38.pyc
index 24413b3dc8c67738effaa72b616d9ae006432c08..221662d9f3849a243cc88cf880af9024e88ac27e 100755
GIT binary patch
literal 10317
zcmbta+m9R9d7lesIOLLh*Y4_KX)HUo=&rP~>|U~U?AVrsDzP?dtlh{1x06}UaCfxB
z8S2b%wKRqgg<`Z1wOgbJ(x7?J_ANml+K0aM4+#3$hd#|y(G&<&pg`XO6bRhkcg}D)
z+@%Z`CBbvQpYQVfzQgTivo7KD<XeCFeeFF-`Xyz`9~GHz;urq{fh)OkAXzex6-(i<
z+O1fcOm#}2cdJ&Fr&Xj4%iw7_s981CSAu$X#+u<}8qy7`fwb;ce;`>+_hq-{)_)*d
zv+j)BK*^ljbY~IIyL0Y5;stlXT|~U-F1gEym)x(ptL_?lEW59`=iKv1opCR?FQ9J4
zz35&-yy}|ni-_0U%kE2v&${x5QtQ=c*o(B)lCx^;9C}?Rp=1y_u1LOm(+gui8Q6Z9
zczeuA{3x{VMSdJ3?uTuUB@U_v==IXHrr-PE!*a8>PSfR|8N?cX@8cK0EK5=%J(OM9
zRb2HmWvD08Q{;|iSxWV%id(@~K30d-B3D7qC~_KdwIZh@SH}oyTBQ+?s-#jXr^>O)
zbK0?lajQ>p;(2}5HS(02j#8qJI?>0N=u<2Dj8dWx)>fLQHqobE^ckf@AMBvir#{hV
zhWeyBZmhBO3`ZrklxGb@&kRcA$~yjD?0@1tQ>c(A^;rc$rg@M0PLOGAKeCTRG{{us
zKB>k2Ue_^^Y#Eu_3zAI77=7ktGmh)paeucLL`hb0j=VUl_u5eqv91^IX9iF1IxMTl
zUSd0O&ub@{miS5FwPrH;K_-91W-&RNqhOwvAttjiS*Bu^Om&W4r{IQ#;0+4i{F%fq
zq3Qn|Z0tlqA2V$@Ub62oyS;IbMGrAYyz$<5KHMk~-G&o7!65eIjhMAJcBA;JE2#Bq
ziCVAroc4jU=WXnJK@Y=qFaS|;e62Ug8nj1S>^r8F#=nC=(q`qlJS)B`r3~cjU!JP4
z7txaufw$@WM<2R|g34b*hUl7Tv7;P=9+#z|>dL7~RR5_W!jAHml&ZMvit_|_T{rQ$
zoiE!=BD2jrC-F?j-1S3;4NQkIXJ9_|lYKLc!mFQnEHb^o>v~}ln@$jgdnQiO>_mNL
z?neF4b(rtPbw2L9Y%kvADVil@>t@(Cui*Szn`YSW?t09OIwp%AW2DeEargah7*iPy
zIk5n&Y}d;d=N3bD{2;Iooq!LzZbm&id?%pcM{7lG-=mG8yAy*hL(iowpCMu{)-c-6
zHS?y|aryx^h6b-+6Px(1m-Ja!+Q%;MW%GmC<m23Mzv?lGzA+sKefUNgu^hkKXJ6pc
zQZgm%WtAxO;#LJ|hF?1~Mtijy<pprM7J-RjvL$hUcDc<rP!V(C*j$vC6?O&rHj&Qg
zQ$?a6qZ`O15)l!J%g@wdB~{W2NJ&do9_wiZ!c#q#QBz5^RPP9KI)ZrM8!D3%xY&fW
zXMWq`r|N~Bh_yY_^-;7Rpf)F(0j@uEea=1qAL=P+gQU`R^GXN$CRsP*wi7rEgu%Ry
z$2=$(Gv{aX%E8`7_mAIb6+OA|B6%{ce@R4`e)uKBna?4gV&`_??0UiGB%PDM^E(3{
zY{k4i5+&wi2r6&NM2|uIMAd$`7scL~DQMxvB>4fDH*R9@f|qg5q*uihCl7<O3)nH?
zgk5macs#ynvloGj<0&ULZYWOvD^I+=_5Q8z-hcPb*8Pvn|NA01`?O^CDu_3$>_$=0
z5_F#FPMi#SUaQK!hC22NHQ*ZSx;0a}=xnZ7xwx0CIdP8G(ztQa@+9-=q_cM0_b^Ay
z;oib*&5#YXuCQw;Ym@#QeQHP$-*JHQF%raYXDL*8sPeGFLyd>}S*fFctPhnTcyg!>
zb-)#gIMS^yL!^_qPqk&~yvXb40YSjEMruIV8_t=j!rlqm{ea~Ch#3IGT%XJK4vV^A
z)3_Hn?OfdVV<H-ooKTh~X@6o@LFaeEUOy>G$?Z+E=O1~&<PZhLct_BeK4}pruzL=3
zx?bY3coWUJ1|`)4JwOs2oNEToAnGSvOGKETnd2^eNFJP$)O3SK7*Z&1bF{ad<KGk<
zpEbMAUg&o)5A7yXE&FCR!`E*2{G`2~Rb4Mg9Q%`Ob}w+;IDurgT_^A|HMFmD%HxF4
zs%7$aw)&i%W%6Bi9oKBB*e@d?%if^iO$s=oA%s$r<gZ}g;&SZK!eg5*1<YN#lo-8)
zaBsY+zO1m{L%vPGZS+|{f`mUOw7&xtfM|3gwKa-CY#<;Nq@%kaw(~?VGSVQq=hWAr
z9F$alq7ViK*->gnjFgW9LPfqxBO^BQxJEI?Z%njBtJ7`9<GA%D^g)Xm>OE{u^qn#f
zh+pglvlV?CyzguxKZPPo8`v*K#U7d>wMyyPQ5tG(j<?i`mTIY`I?;08ZKiW+)tx;6
zszQ3cl%`SCbnaN5lSCcQ&*%9$^n<;{`^}->yhx`;T1)Hc3~J{)gvb}1A5K%Nkl}|g
zk0zOmP|HbLA0ex&04)G?iO|BmeIQ%ntA6MwaSn7Q(b@wL|Ixbn5k<KYfd*+u#Xw*u
z-z62uS!EK^5VqM~H+RMzA7H2xt1Y?~LP3aFbQFm%?|woNi&YduJ%<a;D<ogAkvL)f
zHh`T+QpSNIC8<{KbW+05mby-SFy7JxKFu*}$u{7vpMwGDG^QV>5QD`wxCC1#L&_eE
zXBLNfLaN1+VvW{2R#~`7o@srRkw6tA+Z#~67?wt)4_{oS99?IxL#%N5mn}}nP%f}C
zWpJHsP>ZpgBLT=rGg?)#w!2nty+P>y3g!HBdvEj4umAh!gIjNNXxyrBp~hZ9WmbtF
zu_QD6xB~=}c&&v@<CL0doKo3uQzLl?9tW~O?d?%whmrA>Rp<)YStQvx3JO!wGCEE>
ziJ1L*k-AZ&-Y7$bqjF|&5(2!MbLOOwwv)&n@f=(h;JH7=?6HeLS~O)7)?dA<%gPrD
z>_HtNeqU57qC_uB6jk|+qSo=#e^b#)?Z23;c};#zZ^#XGUY=JGE9@5LIb}GK1;CX8
zsUQ;S59KGoAO{LtPn3ti6Ho+F2WBH2CTgO<$w8Q34Pg%qRc=D6&ZViSLYR=mQPcTh
zQqqKMx=mOx!8D^Ww-t*vx40F^!9KWq8f((cKZhL!Uvt#Bd1}e98ckTkfRvsS<`DYI
zeFmcvmQ1<zx_R3PVh9uhD3ATv6Zk?vQzes~%m}^5@xB8zs%2{6g_*{4g0ia&4iMC3
zwS+lgOqd>&KL)K8+9JD&yE;j8Gr7%2+xTbHkFck*h`X;S>@LbCuqSa~fj$2T8FG%L
za5{b_50w-*epZ50jNdkLfCYIDH&4I-mZ0K7`3`l&q`+m@Se1~k7Wo>_14_7c(9=&7
z71UF8XGrZ*4;Lz#a%&POs>*APQf;$TYoNz0@^h;aTyU^J$>L0_q$oR8D9Cx|^Jz3l
zUPdm{C%}@)LG3-9uVeOD1e=r)t)PlYqJ*{Wg#e}I6=(`j6y;4SC@8KxaM=xMLxk@~
z?x-UlB-ak*E+rv5Fh!T=P-j^ak3jWFJzT7cwC<*PTH6*4a=BROSkZWWLN}MzM!p;&
z_Dz!raMOgtEI+qgAC{R-kUV^b<gVZc$c1gDan2Mnmqavs7eUq>E1JQLo0)u+$pcHd
z`@T>S+OF@!M;4C;>;vjWI2OFKko%Cw{x+skHF>}|Gcj^`wCZdNQFf_(>{D(ktGG_$
zuy?5QHz_C<$ljwY$w)E%350%;mT^MRw<E?pSX8hieu;kZ8Uks)p~B&(k|{_c*^n>F
z7gdPm9U7TD)xY}ezy9&*uYd9PZ?`Eo`dmh$VAFpAl?(6QvHUR*_ppH+F)IW%FKG^E
z6Cl+mgkK?^5zi-;hZ;wc3K$$%AA)=#sit#J<eSnT{q1lOmOwLI1Ov}f*a6!Qmmu^{
zmMw$rmezo*!M>GrS)^9+&afm6&j78S!MA<@l#3M^>6v4^Tcov5X|xqrL){9uD^}B0
z)YsD0V;NS!8rD(c>sT9}P1n-1Sj9rRn66NYe3Xy=_0b1tPguA?-pt{-iQM^g;pt4O
z!OM9L>z(CZO)Ak*a(+}Y(c=P6qX8j=UAXZ7<2<TxZoPn6Ux1ym@U(eM7RJT&;?vn<
zzPI$k)NwBHaf}n=m>9<_jq_sqVlmF8(KuM^d5n4~y_mkp(Y>*CAH1Tu<e}h1zbEeg
z-`>d+|NTy;wqpjLpbl+Aj<g48(V7!feh2AH<NhJ`2b8YG?LiFFIQFtSy4l?GcErBJ
z+d=B>d0EwgQ#JyIQsaI%)5xD2XQ~Ttyyp7OUd+BtjSWh3SL)?Vc@SpGM`2d)1byGd
znDGJ7Uk~^WzkbhB4tti~?-6HN`XP@qZOe8J_pAn3FR{H|?88x<DO&^lz~m1z+0QCl
zy(G4jt)8Xbb$a$)OKT@~4=)8rfRDbL=>rj0Xbh{?4xkX>bltNIKRhI^ORTw3%6{zm
zd;8cEV1cD<hnBJvTKYDBMY8nwqeHvj%QPNa`cA3%FfTqlsZ+UiW8b^KV;_1!)TURe
zlCZLG+{BZb6qB(n&Cf#75^pcd4S5+3)(eQCdI9Ac@)~hI>K75O$n&Tj#q2R1bfwLL
z-3e$$xRhRBNbJP_QQ~?6U1YbzkFUU+uX26`C@Ym!c$t<c?E4et8ZXx;%5+{<4X;d;
zR;RQyL`%q;8kG6Z(Y{U!8vJeKX{hsIEnG{~hj0nIO~44Fl&-;Hj5o1i9kGrWe9Q9<
z#FTFW093)}bx4+3SmTfn;CDzfls41GL@8w3LyQeMg7+=NkPqPaS$7V+JYUo`^13GK
z=6Kx#+RnhvUN}}DJr{>dNdwYm38j#BXjOt0?!3G3SsjkH<>8s(%5ZfE3FR)P%V<k$
znoF0D@p1^8ela~$ltJp{Wh+J5N>R32ltKDYJDNLPCfS&-a!Ih{{2EuyRc+aQ&S6Zh
zzr>RZM=(H$?E~_~jX@0ejS<{ty@-4CfI_?E5pXyue-(^VhAJHIlsv@*=JRz3XJCN}
z%B2%LGKqsq++)IV7$$|sEk{4s;JX;Vz7SiCr@%WJ#9Etn40m3Xn4{d-hrk=3@&G>I
zGKAb1BvYoC+*jZ?5TJ!n&UW5>CvSclN12lyax|596S?Pho`hAUku05GXuKt!f(V{J
zqxa2c^At2G;D+iEQtTcAOWU^N#I>{?is1LYO*I@HsW2*lX7IQQM`aH#%6^xTI}nWy
zsL@-Lo@XCYT18;hf$Dr9ARrcDh}5}Gumk@99?q^0k7(DS=UsZ-%Z$K>Z32%lB$wyo
z@fW}@Vwpj&w{*lps@3T6N{lcQM0@0Fvco80Z`0)GI4rYfaF3%~W%p@0+q9hO1DhX)
zRsCp`*d8Tz@<elV?4o6!&1<vD5iJ@aR$UqJ^Z;qAj-#+KX>6C%1NMmKW;C}h$i>nH
zxj=*Pi-pP%abq+u?bOP1qK@@xltT(mpv3c2P~wPxZ;=xD7CK2wb*PVJxq{z{f)qd^
z6o>-rU4v?&uf8I$s2Ai5>b$DUhNP6gUo2LrH@*t{9HX@fHjX}Y`*hP>p<F=%uuq1&
z6F=DaDnh&=z;Q}>jhwB3LDdwt)z6eIPE&kGc*rMB+!Ib8?#TIFjKw|Knfz_;$F>Y|
zFOj$U<RPq&o!PmEk~gE5hSNT6GU#IdKB>T={X@($@mfhQmGqLr_e1!a_5)`%-aq&r
zRKQE9p${eY7pN4A*95%(Iu(%D%F^$Ob#hnGi5t5-^5MV+7Y~vX{u{ddAB~rYBmXOm
z*rsKTKE%PqwtR^_B!EO=U_+~Ec#qdkzQ-4}cv>Yh+?8DEQ-!oH3}~Kq@uGiLf~UQj
zD70Oez_>b?!d3ipMZ=fS><l!U0VA0|_~f#LWC9~Y;x^~}y@mh0`X@j6@uP3xW%*o5
zB-@iBIm>CAt`dLkDCY**oZy1u6|AKSH2P;W6<r$)uPGk3tWLxsG!$E>w%nI-l~N?<
z*m(+0T*al+bi<!OMlT{irF$T>B<Xc9O<ldLURKyokU#ZrAbN1)_wjQi5+WM>=LaWT
z;SX*V7Y}YvD)i4FoG$(kbgHJB;2#x5YkHEakP`$Uxj-uP&Vdl`4o|8;Kls-La&&{Z
zA09rhZ4Eb?J2X=5i3XH3%GEm>931s6R*)^*#M;^KQ@}ZpyH_XQl5iVDj+|}yL9}M3
znYP(V$u(DcfMb76ou^!%xT6!-^YU}9r;G{We<<!qJtv=2*X4C(USWTOd7zWH&QiD$
z%)TdzNXxa(Wg5hFZ;-1b7&7#4n!w-97RG7(%ZC4tNEb1qr>@M1+kS6w&5sxf38CZk
z1W%2AlKs|A;WWQ9{=?aJ;p-m!;RGz<q|At`#h-6Ta`I<f;f=mU*}D{QE1tkS7ro-(
z`4tcu@_m-5njfHvip+Sc8@c_!dz)A$rqLnfN#NtlW7W7?Uou`YmW_9f^G2)w_4*y7
Iitm~K1;8Y_J^%m!

literal 14057
zcmeHO+m9nxTCbbS<x9J7p1#cVOe)KUot{lkhY8tLT1ht9orDc!C)vz}#M$7j_9=IJ
z+Ace%s(ac)86L8Wgv6jo@B-pt<|Qa^D<pX04TSgucv#f~V#5O(B%VQGlkoe#DwnI=
z-I+<kMM&#zpL6Qm>Ri9?_g#E_W~Qv*v+?KMA9^E2`Dc2_esnxMkFWn%2!tZkuF_CB
z)*2eedc)wj&?s<hHi{fu4U6MauiU7p)K=?Odozt0PU}e58a1R1QTUP4m=&fdBAyeL
zC?TE~Wl=%AAgW>p@uD~@=EVGu)J9z_h(+X_5_NG3@sc<vmc%kjmc@!#Mb3(NM4S<8
z$XONYRb}12kFhIP*45~1mtD{A1jA;>3)~&)1f9Oue68Q{eZ>1-%awtHq9LmO(s5OP
z^OdXFYOMlwoP0`%>G?bO`fnfzl)EavTBzJpLnGAc%AR^p{ed<rjLftQIVjW9vf{`R
zYN&*T&}?hZV7)^%EQS`Pk*l^f%tv$nMpcxuoqS%xz-c0O(6>9n^@2`&Xgi*;Z@aR;
z(ccej(eVder{(rg=-a`r6WDF1+qJjdZvU3;*vY(Y-;*enuD{#wigNmbm#@5b{o2*7
z^m)U+{<wW(%N}$N-R@N9=aG4Fvn`#TdvQI}!td<#oF>}E3AD}*%ATx!$Mx+@C%Ha6
z*-VGGW%oK>r?=m;F_Zo+XS?h2_}hJ&gPZRjE3@s|tv)tixY8B&fz$0|X87v$YsXJ9
zeI43l2kcK)=Y{LnzWeIyndNx}eFygI&P~@7((&x`T^XmJAM8u9?_PApemt#z!71U-
zx*loXAS&G2#poj4+4k2<k$OE+Z$!mp<;IB{(G0d}vJMS%zbX5-e2QEBF5-ja`GB84
zP~;5W@zLJqjed8(=lYwDi}jJs*5+%nf75LR{^qw|y}CKsf#Ig(Io+Y(@i%?h+T8B@
z8#q;TZ8kF3W@F&A_M9DebJy(-&|Nz@h5p6CFglYPVl%gLRgA`GqO-83uIW{^s>xa8
z<t~$n3<F=jOy5O<E=+up_LNZ}z9M`*w6O3@T*_A|DOD*KwRDQ<;)GhLw+*^HS0r`3
zu1UJAvPfxDT9l?6ENNu2f*@MT4tuiTQL*ERPRsRqUE`oe^_SDblArK-6#7&{(Q2xt
zT5(vC3&?5FJI9}SBvK0U$9RY*9U6NErcl^Z?rEbE=c%EIJWvMCD|22cETg=_=}KrJ
zT}4{M#H*{IK%*Ic(vzne@@H0*(JYpzhIkH(H5*!Yb)nr;XBA$&x$HAa=a~BF&K?$!
zuXdL5wE9xwFxhnf9B*hJbhcyqYWLgmh3p6RzzG6ZdKc~QINjEM*9mah{Xjaysmn!_
ztpkFFdo0}@qHUSG8((HBYfIPNKzr8#6#&8W_qWj*PBdC}d=&S4+(*xK{C&Fi*&X$H
z3R?pQw90;C%Ra<AHu_R@JPd0*=mSg%6KZ!{5JVc)Rxlj6RQ=LryWQ=Bk~$LO&-V4^
zmOcC&eWCiveO;$<Wk>yhfua7%oq_1)7y1F6q5jF89fF6po1np4cGqnOHi#W{chYuc
zwkorZnNIwc({<$I@uFszGpXYoZrQT4v-{=VM7^B&CV~aYl(zdpfWz~JCY9;tL9e3W
z1pUG1eouU8dyc=iWpCra%ArjUICZ#${-6&6<?Rq2h_}PHnX&L)Xvx0sZ*WE2Ke}4C
z9Nz{7$Hncr*+YKZ-t*~CCV0YHN;lY-I2G5yek3;Z2FLbBJh)bH{56!+^YnR(uPD>&
z$k^`pyXyv1?Wof2?{sirnz;3mf$Qx>rsGo=>+_NNMx?$uq3DfLddkwIbJ!@Q+mt4q
zL-d7ZjMT$O9Y*SHNffzJ*>9%AtWi10Jectu8)iJmX@*cwoADBWlwi4^L6`m#0%fhP
z)>W`Sy{6W5@CAKNozvu5l;pU=DP#}_XIA)+nnLWXDAZM@4n$q2w7#mGQAWh%j>`uI
z4b7YCyQ<LkG`SmSH}%jAm7928MPG=GRi&+o!X4r++Fc_ofm@V3BQ*J$*)C9Pc~lW*
zSQ151{IMnsMyKg>SP|L^dI8^{JR`0(eh$l2I>`&IkZ?cu)TU9WlU`R9VPPayjHD2j
z!%7?c2=gk@d`C0k%w6-IBFcNld#bz;R)gYQE3C$I;26E#181u%-&D8{aK*TfyO1X^
zubQZawfJoqGdLiRG0hX&_!(nH>3EEo7uC@w&lq1NtYNLR5%JRe%GWTbqL`&IBc7x9
z%y_rvX}6HOz%wHL`$gN%L2~$_ZRa6je9^WGq8!d+m(4xe?{c^>P6t)7_^zrd!3;yL
z5@80cQ59=9k6JULey4(yv^c7BIh8JQNj)r4y><lv?UeIp`U%p_K)RrO!1hkyFvAB4
z=rUDH04hL3D%XMD2zLV0j|wqwJZ^*cW|{tu?1)(S$ci2RH^R=KbI^1lrlbU(amHrP
z31sKczSRkK?Lh9kB+2XpMG@@E6Z?lD^n_?OZ1)i4zeB<yC@QlS7bwjq5sgU}u@~Sk
z5{-$iByR+>r`69k2EHS~j~X)vkGq$<eU`W70Q6n(7Rf=s>v%C(j6q%-f)C_xd%ThR
zKKKMPFb@(N@3BAhP&eqb>`o7YM@9rYMtqAmcd!dF28U$4%Xp`2^DdLPOj2Sj*dWgl
z9V~ziIoWgV^S$f=pa~p4Zc{YGqD(b`xVpEq+54?uooX_v8J9nqEsuq|q~pm0y3fdA
z1c4vwZD|K5Hcy9;?#q`DTcZj~{>?!rXzkvIDnf}Cup?fljkRW#3N5Fkf)BJvbt21i
zZ~42<z<uco<gU}mmuIL|B^`M43Fh*05!u)7gUKn8-e_KqW^(;TMd5Y>r+GW7?Q|XC
z2Yt_N3a1N^*=s%#6|eTY2RNZ07}4p>GB)4Y$G)e7i>mQkX^hRo+|w|cjmtzwc3Qh#
zH(0Mlde`-$0uMh1A4O7-U@}n&!iER5<-3uAk@%6`1B5WW0kTC@pb#xd3*$@T__9v(
zDhwp3l%zSXpOVXzze2%rBB(W{0?~5zo@Fj+RGwTfxk~+OSBCf=Ml)29P9$2Kto0;w
zHcd9W4zxiSG@bQ08fN$BF#9l!bkvJVX_wQ;=PS5#v3%D|4seqwR#tb`#aR9K5h%+m
z>IyK>IVb`Qy=Ke-72&(~RbvUcwqDoj+8lowDXE&Oeae5j_OYg${51?<sN+k_GS%1S
z5E`FohW2yKF!^f~PAoAAm01{trr>Rg_DR{0FJq)F(p`){bc|D6W`SjWeq0u-G{6z7
ziF95KXA-;vOa{DDW0`%DTZl0Oav{gZb?+L&1g@DA#gKHT){nKY%ot|!3=G98k^=ax
z3A{3Q1g`*H0S%OMIHi){6lf1HlKCSz<-u(hj^Go#0q;T2jElHuJTtZff?Yf^F4RQ%
zRJ{cCqAxMC7-K2mEsU{#r2R4POK}+9nzm=8HadlMHN!>94eJyxE8(d)#tJVTK}^T=
zo}0}yvirQVnQni~`@ZmOmc?AS6q;h5G2C3ZJWdBdri|f$(IHWojL(4LaK4sMYgQ~W
zhFck}ayj9;6)sst-MMxd2(JFKr+$IpHdrG6|BBvHnD7f|MRp5tyZ;?%CC`DY$VVu!
z5ky7e$Iwba>CXm~<T1ucSL82K#jjDYPQe8V9z!rT7JP5aqk*LuV7(;2PIVqfuwIHa
z;Fl=lDGHvZU;<PO{~Z~~Ey`rDa{9Bu3ZoB6>W#>d&<j#X8kH8OTS&|2sIPBP@H_=C
zQ1ERENVg#`BWRcq9(v7RN#aXM{G@!5YP>|juTelIfoafj{=uN5=|e)<`&TsfKSYo~
z4}D2rD%J~7zvv7xK4plZnuyI$4XEpD`J5612LO(^27rYT=a^}Z#^9tJgOeC&#PvVT
zmz=Zy8~Y+#$~0`bhFM0%L2Tl@)?zaR%^V*(WhsF9hc5?UYNJ_(Dz$Jn2Vmx=a$^89
z&j2Q=n*f-Fpa2$Qg{5V%@^{tIB0ocUlRVeKTCG*(vhw!dj!ubEcxum-KZWvv;=dOa
zu%X2m_;_>l1ha89g+en_3koAEssA5IeFjFk`kW(bN&T(3KB=0xMJcQkl!{wG>A~+{
zur<r$ct5S*;QCfv|DE^ar%C;H-e;(>z;JBg2yB7E2(yKvD5Exlx`@?Aiwt#^j=-@8
zw^>G?P#Y%j3BBW)l~9PeSHy{`YP6E8i&{tPO!^V^1csgBv90FX9`l|QwBT)N+a+#$
z8tYO7Y+QkI6TM<(R^t?)m=>NMrvb;#aLqzmQ_s}Yb2Zmc(+bbrQ}33bI$Yk<-qU?D
zty$q(Hc!^X3d8U+*%?P?!<qMVDa3Mk2GD4QVbXHAmPykHD%3BV%TUJ8=#!sU(e8{m
z&FxmvZY@qjMRtbZ6_jFYG!prziL%UO5VZh^K7*EP;k-D@K<nJ-toN0mmVhyWV8_<s
zaw<K?t<M6CEif26Cx0FHe2TBg-=Mr#C?JV@2VZ>t>05A}7)!p*?f$+ep!#XH;W^+(
z23dI^;daEeE07P#5)<jLT|{PVD>;!D$?u_Sd4&QpyuP%v^#@OU_=m$6pJN3<%#W^9
zDKn+n6PZ$^lh{lmw0xBUvfRpVAc)E=VmG^;JvXvCej92kaFLX8Nsu;W{gH>|Uh*~S
z`1=&RPQg_KQOW7<^d&^tUgOk2_6N=mjD%Tn_C(&+C`eezM5EB6C^Bd@GyKlNGc1hY
zc;nQiG?nFg%<!1?5y{S2oi98D*V}PDS32DvKqKS-6wPydkIW9>9Bi0|m|P|UHOFO@
zVhcGnZNYpf&S;n<Y-c&69_O6I88h{#IO7HNf_`2#)Usx3#Mk0c5`QaWRAa%s3a&j)
zLrDK2=GG#;)A&Qyd^#uuF#Jytz~BMDfK_Fr1}co>@EM3r-XkqgZ^FzjbfLq%0Y@9q
zJQ&Ym7%^au&~v4QsZt|X3dfyX+d{6iI8|!qO5x~}Yg^2fmZnOrTq(r<T-y?Ogh7r(
z<Ukab+C_+F6=xY(j?sPIc<&v+_lZ@aINrGLK52TP=-ThXn$Cvw%P-t``6{gSeKMp&
zGwIkDp@)ZQ9F8|o)WTZG8tPc{OX}#c>N2B~zZmP|v%2`i08c9YE!*4g4dA;(7WJ&g
znKZR04U5Vq19^)`7Fxf}Q_Kl3TZv(kOWuI;5-Qy{DPpT*Y6ftFSV8<2IWoO%=R1r|
z0hlUlpmT<K-#)fQ4BNAQPiUPp=ub42<oa}MCoi%cV7$U(t9?8Mw);*-$+OGpb|2my
zDJIBo;g(#Hk0Ola5RxbE6A8Uf<oiB3sN5&_ob`fa(Sb5zip-+JlT=iW4BF5}1smFo
zeJK1$1>fLxj|#l*$7(zcn(H*OUxxL8&L?!e(EbMf<|WZ2E&hi@(1;@_^)hUy=hSnP
zuO@FIKj&0L8dmx;G584*0j!_op9lM6>>#8q0Um`0Lj7K0RN!*B5ZzP#Z-BYuGTj5~
z{)|=z+L`3thkjvYMLUn%A@zWE1+*(B?TTne01?(bY8PMOwvOv-{%i|7!Fz1ule@z3
zf`M@Z{F_}?=q%AC^c~WOk8kqH=CIa2JO{j>nOu2N9CJUgNRWd<*F9iYza*0`JbEVY
z9HtQIh_WX;%O0Ds3hnj<ZaC~i;{%+noU?nlet^Tx&Ix{N7jBDb|B0VMGv{9hPsG@L
z0rZP>^Kf3IIm1AeeS>drJ7iQ!EI*TLH?+@RgC%YAC8z7V6C=&kD<H-Gk9ix8b>y_9
z*_BT>A@ScBQ&jSG%8q|br*C!Gsgt)x*l_E>A~J!VnB&H;YIk<v5!HOyH^$ZBN0NO*
zx+~*OGh=Fg&KG2wV#GAI>6fK%%TYJF1`ki+ZlBMdh>1@hj?1%;`N{LMk6paf{-~BJ
zT=4v4n@msa0~r;Q;eE82<ZNWd{y`zq`T)ytblZ)Ln|-K3O!u(k2YzJroIMvi>aSbT
z>WSAadL%c(OxMX3YCJl=1G$RXlP~8$+vHafz=ts&_0&d3XUF-+Zd0^AF3!IH*s9Z{
z`gDuw&=EQUdQ=4R@JT<UoPyGWR3)0JXdfGes#?}FIiHoM$=^o%7E$r>hkP~3i;}7O
zUn4P1p9uuDiKc6znNW5sp|ln#?G{s7q>fVPsg^NyhX*E8cX-M82iz8^hiFR_9-f-q
zwv4ua!EKQOg$F!?@`5S6+`|-ZZ|0~EolQEFpbh?Da%#ufKq|6}R<cKdM?s&00fP04
zJfQR~1dWyP(aH6W!!N5;y!sF)iv5%1A$2*V;5G$QB8h$S8ED>e#g6+lodl8yl@;3e
zF@QHnSolw)I6e=Vd<E&B8?ej8C*{}@O@`<(2mk3J5LV%{@evYOa}`Pebx)ToA-rRS
z7U-~uDS>gbK;Zul`RA)62nO&^GRV(~!$KQARwF&sgb9@elvpG9^K-h0bRkJYJvB0s
zCd(4)h%#(a#pD_F_}RiU>_BOGDR~~3m%~zE-YkkjXo%_^_!ZS*g&I}Da*|u%T#Iw7
z$g{!<R0%Vpc1Op^XJ9cb-7Md%P;RJ;SuuCVfRFDiu*1CbDBu!P(G)kpfmMKqPdi_*
zyK9X4f$}{mMGW(a@5Cq?-ZxoC4(>HU@>CtvA69d)`aHvDSs~*|lJ&rd-}x52d1LH8
zhU@HhK_UdZk6ebk<AEbP1a<J6h7No+2nY%mNvIC{mBcV7$0w$E1e{inI_2{$qQWa3
zZ~yT2=}o@{2j<r9CXD@_*Y`GCiRb(#ejQ?sCwPTHGyDTOx8V`aihBag4GWs)7H*cm
zbLZdw^6tO?-A|soU8O95ALvdx-nxc1fg^*+fJilnv_tyQi3QhkunTV|GHGs(gu~-$
zbjR)zat*1-*h8K5YGlDvB_N;8$n=MP;PzlEaQRmrjoBVlHvr3H2p1W&C(y^c@X8>#
z4rfU?C+t95g5Lw|(gI>ItqvUbEAdKW6$7_{(a@4!66<q~ncTKYmwVYohpud}kpia(
zzhOAU7w!k%ZTtv`T{RBJiDAP$d;)Bt-7tr7TsY*%!Z42Bh@&?tlF!nlr^w0@omMK1
zU4_P;1!SD?ZPZt4&=Y8L(B>mfGzO$KSW=hNI@}24A5o{YMJ<5n^g~h&Eeg|!3rG-+
zAen#_+3>I&>GVSa`2LY%fJ0LGhD5E6c>B52`=8=7{#gOpQ`QYEgv7VS$EDkZWD`}l
z;bM-TOR$tWJ;TL>aAxb=rRfn3Q?vyYllT+rS(BtOySA?4pI^RyIx4X69Ph7=BM}w*
z_$>x-h<pci<?kbC%=0c^%xto}LxsPGAgXc)7s+j<k5d)L@|#rj^~*7gczNPpcYX3}
zsqg3fYuH%CIye4V2mj0=<{IQeA;}=YdK#9`n4L4-j(3M<msu4)0c5j?tY><CvEOx{
sBkhiV5djod>0gypM4Gi`)vYJ2-?4lYe_H-srD#27J*9kHe%iYDU&80b1^@s6

diff --git a/helper_functions.py b/helper_functions.py
new file mode 100755
index 0000000..87ad450
--- /dev/null
+++ b/helper_functions.py
@@ -0,0 +1,446 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib
+# matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from scipy.io import savemat
+from scipy.io import loadmat
+import timeit
+
+from DensityIntegrationUncertaintyQuantification import Density_integration_Poisson_uncertainty
+from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty
+
+
+def plot_figures(X, Y, rho_x, rho_y, rho, sigma_rho):
+
+    # --------------------------
+    # Plot the results
+    # --------------------------
+    fig = plt.figure(1, figsize=(12,8))
+    plt.figure(1)
+    ax1 = fig.add_subplot(2,2,1)
+    ax2 = fig.add_subplot(2,2,2)
+    ax3 = fig.add_subplot(2,2,3)
+    ax4 = fig.add_subplot(2,2,4)
+
+    # plot x gradient
+    plt.axes(ax1)
+    plt.pcolormesh(X, Y, rho_x)
+    plt.colorbar()
+    ax1.set_aspect('equal')
+    plt.title('rho_x')
+
+    # plot y gradient
+    plt.axes(ax2)
+    plt.pcolormesh(X, Y, rho_y)
+    plt.colorbar()
+    ax2.set_aspect('equal')
+    plt.title('rho_y')
+
+    # plot density
+    plt.axes(ax3)
+    plt.pcolormesh(X, Y, rho)
+    plt.colorbar()
+    ax3.set_aspect('equal')
+    plt.title('rho')
+
+    # plot density uncertainty7
+    plt.axes(ax4)
+    plt.pcolormesh(X, Y, sigma_rho)
+    plt.colorbar()
+    ax4.set_aspect('equal')
+    plt.title('sigma rho')
+
+    # plt.tight_layout()
+
+    return fig
+
+
+def create_border_array(nr, nc, fill_val=1):
+    """
+    Function to create a binary array with non-zero elements along the four boundaries
+
+    Args:
+        nr, nc (int): number of rows and columns of the binary array to be created
+        fill_val (int, optional): non-zero value to be assigned to the border elements. Defaults to 1.
+
+    Returns:
+        border_array: the 2d array with non-zero elements along the borders        
+    """
+
+    # create array with fill val
+    border_array = np.ones((nr, nc)) * fill_val
+    border_array[1:-1, 1:-1] = 0
+
+    return border_array
+
+
+def set_bc(nr, nc, rho_0, sigma_rho_0):
+    """
+    Function to set matrices to enforce dirichlet boundary conditions
+
+    Args:
+        nr, nc (int): number of rows and columns
+        rho_0 (float, scalar): reference density (kg/m^3)
+        sigma_rho_0 (float, scalar): uncertainty in reference density (kg/m^3) 
+
+    Returns:
+        dirichlet_label: binary array specifying regions where reference density will be imposed
+        rho_dirichlet (float 2D): array containing reference density values at the non-zero dirichlet label points
+        sigma_rho_dirichlet (float 2D): array containing reference density UNCERTAINTY values at the non-zero dirichlet label points
+    """
+    
+    # define dirichlet boundary points (minimum one point) - here defined to be all boundaries
+    # THIS NEEDS TO BE MODIFIED FOR EACH EXPERIMENT
+    dirichlet_label = create_border_array(nr, nc, fill_val=1)
+
+    # # sample to set left boundary to be 1s
+    # dirichlet_label = np.zeros((nr, nc))
+    # dirichlet_label[:, 0] = 1
+    # dirichlet_label[10, 0] = 1
+
+    # set density and density uncertainty at these boundary points
+    # HERE SET TO BE THE SAME VALUE FOR ALL BOUNDARY POINTS
+    rho_dirichlet = rho_0 * dirichlet_label
+    sigma_rho_dirichlet = sigma_rho_0 * dirichlet_label
+
+    # convert to boolean array
+    dirichlet_label = dirichlet_label.astype('bool')
+
+    return dirichlet_label, rho_dirichlet, sigma_rho_dirichlet
+
+
+def calculate_gradients_from_displacements(U, experimental_parameters):
+    """
+    Function to calculate density gradients from displacements using the BOS equation
+
+    Inputs:
+        U: pixel displacements
+        experimental_parameters: dictionary of optical layout and other parameters
+    
+    Returns:
+        rho_x: density gradient (kg/m^4)
+    """
+    # extract experimental parameters
+    M, Z_D, l_p, delta_z, gladstone_dale, n_0 = experimental_parameters['magnification'], experimental_parameters['Z_D'], experimental_parameters['pixel_pitch'], experimental_parameters['delta_z'], experimental_parameters['gladstone_dale'], experimental_parameters['n_0']
+
+    # gradient calculation from BOS equation
+    rho_x = U * l_p * 1 / (Z_D * M * delta_z) * 1 / gladstone_dale * n_0
+    
+    return rho_x
+
+
+def convert_displacements_to_physical_units(X_pix, Y_pix, U, V, sigma_U, sigma_V, experimental_parameters, mask):
+    """
+    Function to convert displacements to density gradients and co-ordinates to physical units.
+
+    Args:
+        X_pix, Y_pix (float): co-ordinate grid in pixels        
+        U, V (float): X, Y displacements in pixels
+        sigma_U, sigma_V (float): X, Y displacements uncertainty in pixels
+        experimental_parameters (dict): structure of all the experimental parameters
+        mask (float): binary array of the integration mask
+
+    Returns:
+        X, Y: co-ordinate grid in meters
+        rho_x, rho_y (float): X, Y density gradients (kg/m^4)
+        sigma_rho_x, sigma)rho_y (float): X, Y density gradient uncertainties (kg/m^4)        
+    """
+
+    # calculate co-ordinates in the density gradient plane (m)
+    X = (X_pix - experimental_parameters['x0']) * experimental_parameters['pixel_pitch'] * 1 / experimental_parameters['magnification_grad']
+    Y = (Y_pix - experimental_parameters['y0']) * experimental_parameters['pixel_pitch'] * 1 / experimental_parameters['magnification_grad']
+    
+    # calculate density gradients from displacements (kg/m^4)
+    rho_x = calculate_gradients_from_displacements(U, experimental_parameters)
+    rho_y = calculate_gradients_from_displacements(V, experimental_parameters)
+
+    # calculate density gradient uncertainty from displacement uncertainty (kg/m^4)
+    # sigma_rho_x = calculate_gradients_from_displacements(sigma_U, experimental_parameters)
+    # sigma_rho_y = calculate_gradients_from_displacements(sigma_V, experimental_parameters)
+
+    # calculate density gradient uncertainty from displacement uncertainty (kg/m^4)
+    factor_1 = experimental_parameters['pixel_pitch'] * \
+                  1 / (experimental_parameters['magnification'] * experimental_parameters['gladstone_dale'] *
+                       experimental_parameters['n_0'] * experimental_parameters['Z_D'])
+    
+    factor_2 = experimental_parameters['pixel_pitch'] * \
+             1 / (experimental_parameters['gladstone_dale'] *
+                  experimental_parameters['n_0'] * experimental_parameters['Z_D']) * \
+             1/experimental_parameters['magnification']**2 * experimental_parameters['sigma_M']
+
+    factor_3 = experimental_parameters['pixel_pitch'] * \
+             1 / (experimental_parameters['magnification'] * experimental_parameters['gladstone_dale'] *
+                  experimental_parameters['n_0']) * \
+             1/experimental_parameters['Z_D']**2 * experimental_parameters['sigma_Z_D']
+
+    sigma_rho_x = np.sqrt((sigma_U * factor_1)**2 + (U * factor_2)**2 + (U * factor_3)**2)
+    sigma_rho_y = np.sqrt((sigma_V * factor_1)**2 + (V * factor_2)**2 + (V * factor_3)**2)
+
+    # --------------------------
+    # enforce mask and ensure valid values for uncertainties
+    # --------------------------
+    # set uncertainties in masked regions to be zero (e.g. for image matching)
+    sigma_rho_x[mask == False] = 0.0
+    sigma_rho_y[mask == False] = 0.0
+
+    # set nan values to be zero
+    sigma_rho_x[~np.isfinite(sigma_rho_x)] = 0.0
+    sigma_rho_y[~np.isfinite(sigma_rho_y)] = 0.0
+
+    # set zero values to a small non-zero number (OR WLS will blow up)
+    sigma_rho_x[sigma_rho_x == 0] = 1e-8 
+    sigma_rho_y[sigma_rho_y == 0] = 1e-8
+
+    return X, Y, rho_x, rho_y, sigma_rho_x, sigma_rho_y
+
+
+def create_mask(nr, nc, Eval):
+    """
+    Function to create the integration mask
+
+    Args:
+        nr (int): number of rows
+        nc (int): number of columns
+        Eval (float): binary array of the mask
+
+    Returns:
+        mask (bool): integration mask as a boolean array (True for integration, False otherwise)
+    """
+    
+    # --------------------------
+    # set mask
+    # --------------------------
+    # create binary mask. True for flow, False for blocked region.
+    if len(Eval.shape) == 1:
+        # mask = np.reshape(a=Eval, newshape=(X.shape[1], X.shape[0]))
+        mask = np.reshape(a=Eval, newshape=(nc, nr))
+        mask = np.transpose(mask)
+    else:
+        mask = Eval
+
+    return mask
+
+
+def load_displacements_correlation(filename, displacement_uncertainty_method):
+    """
+    Function to load the displacements from Prana processing
+
+    Args:
+        filename (str): file name
+        displacement_uncertainty_method (str): name of the displacement uncertainty method
+
+    Returns:
+        X_pix, Y_pix: pixel co-ordinate grid
+        U, V: X, Y displacements
+        sigma_U, sigma_V: X, Y, displacement uncertainties
+        Eval: Processing mask (binary array)
+    """
+    
+    # load the data:
+    data = loadmat_functions.loadmat(filename=filename)
+
+    # extract co-ordinates
+    X_pix = data['X'].astype('float')
+    Y_pix = data['Y'].astype('float')
+
+    # extract displacements (sign conversion for the SAMPLE DATASET ONLY)
+    U = data['U']
+    V = data['V']
+
+    # extract displacement uncertainties
+    sigma_U = data['uncertainty2D'][displacement_uncertainty_method + 'x']
+    sigma_V = data['uncertainty2D'][displacement_uncertainty_method + 'y']
+
+    # bias correction for MC
+    if displacement_uncertainty_method == 'MC':
+        sigma_U = np.sqrt(sigma_U ** 2 - data['uncertainty2D']['biasx'] ** 2)
+        sigma_V = np.sqrt(sigma_V ** 2 - data['uncertainty2D']['biasy'] ** 2)
+
+    # extract mask
+    Eval = data['Eval'] + 1
+
+    return X_pix, Y_pix, U, V, sigma_U, sigma_V, Eval
+
+
+def grid_PTV_velocity_uncertainty_2D(Xn, Yn, fluid_mask, xp, yp, up, vp, up_unc, vp_unc, N_avg=3, dist_epsilon=1e-6):
+  # This function maps the velocity and its uncertainty from PTV data (on unstructured particle locations) to
+  # a predefined Cartesian grid using inverse-distance based weighted average.
+  # Inputs:
+  # Xn, Yn: 2d array specifying the Cartesian grid.
+  # fluid_mask: 2d array of a binary mask specifying the flow region. The particle data will be interpolated only onto the grid points with fluid_mask==True
+  # xp, yp: 1d array of particle locations.
+  # up, vp: 1d array of particle velocity values.
+  # up_unc, vp_unc: 1d array of paticle velocity uncertainty (std).
+  # N_avg: the number of particle data points used to determine the velocity for each grid point.
+  # dist_epsilon: the minimum distance allowed for setting up the weights (to avoid signularity).
+  # Outputs: 
+  # Un, Vn: the 3d arrays of velocity on grid.
+  # Un_unc, Vn_unc: 3d arrays of velocity uncertainty on grid.
+  # Cov_u, Cov_v: the sparse matrices of covariance matrices of interpolated velocity.
+
+  Ny,Nx = np.shape(Xn)
+
+  j,i = np.where(fluid_mask==True)
+  Npts = len(j)
+  Np = len(xp)
+
+  # Setup the linear operator to map the p velocity to grid using inverse distance based average
+  Map_M = scysparse.csr_matrix((Npts,Np),dtype=np.float)
+  for ct_pt in range(Npts):
+    x_grid = Xn[j[ct_pt],i[ct_pt]]
+    y_grid = Yn[j[ct_pt],i[ct_pt]]    
+    # find closest particles for linear interpolation.
+    dist = ((xp - x_grid)**2 + (yp - y_grid)**2)**0.5
+    dist[dist<dist_epsilon] = dist_epsilon
+    close_arg = np.argsort(dist)[:N_avg]
+    inv_dist = 1.0 / dist[close_arg]
+    inv_dist_weight = inv_dist / np.sum(inv_dist)
+    for col in range(N_avg):
+      Map_M[ct_pt,close_arg[col]] = inv_dist_weight[col]
+  
+  # Peform the interpolation.
+  Un = np.zeros((Ny,Nx))
+  Vn = np.zeros((Ny,Nx))
+  
+  Un[j,i] = Map_M.dot(up)
+  Vn[j,i] = Map_M.dot(vp)
+  
+  # Perform the ucnertainty propagation.
+  Un_unc = np.zeros((Ny,Nx))
+  Cov_up = scysparse.diags(up_unc**2,format='csr',dtype=np.float)
+  Cov_u = Map_M * Cov_up * Map_M.transpose()
+  Un_unc[j,i] = Cov_u.diagonal()**0.5
+
+  Vn_unc = np.zeros((Ny,Nx))
+  Cov_vp = scysparse.diags(vp_unc**2,format='csr',dtype=np.float)
+  Cov_v = Map_M * Cov_vp * Map_M.transpose()
+  Vn_unc[j,i] = Cov_v.diagonal()**0.5
+
+  return Un, Vn, Un_unc, Vn_unc, Cov_u, Cov_v
+
+
+def load_displacements_tracking(filename, dot_spacing, displacement_uncertainty_method='crlb'):
+    """
+    Function to load displacments from tracking processing and interpolate them to a grid
+
+    Args:
+        filename (str): name of the file containing the displacements
+        displacement_uncertainty_method (str): uncertainty quantification method. defaults to crlb
+        dot_spacing : dot spacing in pixels
+
+    Returns:
+        X_grid, Y_grid: co-ordinate grid (pix)
+        U_grid, V_grid: X, Y displacements (pix)
+        sigma_U_grid, sigma_V_grid: X, Y displacement uncertainties (pix)
+    """
+    # ----------------------------
+    # load and extract data
+    # ----------------------------
+
+    # load the data:
+    data = loadmat_functions.loadmat(filename)
+
+    # extract track info
+    X_track = data['tracks'][:, 0]
+    Y_track = data['tracks'][:, 2]
+    if validation:
+        U_track = data['tracks'][:, 15]
+        V_track = data['tracks'][:, 16]
+    else:
+        U_track = data['tracks'][:, 13]
+        V_track = data['tracks'][:, 14]
+
+    # extract uncertainties
+    if displacement_uncertainty_method == 'crlb':
+        sigma_U_track = data['uncertainty2D']['U_std']
+        sigma_V_track = data['uncertainty2D']['V_std']
+    else:
+        sigma_U_track = np.zeros_like(a=X_track)
+        sigma_V_track = np.zeros_like(a=X_track)
+
+    # ----------------------------
+    # identify extents
+    # ----------------------------
+    xmin = X_track.min()
+    xmax = X_track.max()
+    ymin = Y_track.min()
+    ymax = Y_track.max()
+
+    # ----------------------------
+    # interpolate tracks onto a regular grid
+    # ----------------------------
+    # calculate expected number of dots along x and y
+    num_dots_x = np.round((xmax - xmin + 1) / dot_spacing)
+    num_dots_y = np.round((ymax - ymin + 1) / dot_spacing)
+
+    # generaet 1D co-ordinate arrays
+    x_grid = np.linspace(start=np.ceil(xmin), stop=np.floor(xmax), num=num_dots_x)
+    y_grid = np.linspace(start=np.ceil(ymin), stop=np.floor(ymax), num=num_dots_y)
+
+    # generate 2D co-ordinate grid
+    [X_grid, Y_grid] = np.meshgrid(x_grid, y_grid)
+
+    # interpolate results onto a regular grid
+    U_grid, V_grid, sigma_U_grid, sigma_V_grid, Cov_u, Cov_v = grid_PTV_velocity_uncertainty_2D(X_grid, Y_grid,
+                    np.ones_like(X_grid), X_track, Y_track, U_track, V_track, sigma_U_track, sigma_V_track, N_avg=8)
+
+    # set nan elements to zero
+    U_grid[np.logical_not(np.isfinite(U_grid))] = 0.0
+    V_grid[np.logical_not(np.isfinite(V_grid))] = 0.0
+    sigma_U_grid[np.logical_not(np.isfinite(sigma_U_grid))] = 0.0
+    sigma_V_grid[np.logical_not(np.isfinite(sigma_V_grid))] = 0.0
+
+    return X_grid, Y_grid, U_grid, V_grid, sigma_U_grid, sigma_V_grid
+
+def calculate_density_gradient_from_displacement(disp, pixel_pitch, Z_D, M, delta_z, K, n_0):
+    
+    rho_grad = disp * pixel_pitch * 1 / (Z_D * M * delta_z) * 1 / K * n_0
+
+    return rho_grad
+
+
+def convert_positions_to_physical_grid(X_pix, x0, pixel_pitch, M_grad):
+
+    X = (X_pix - x0) * pixel_pitch * 1 / M_grad
+
+    return X
+
+
+def calculate_density_from_gradient(X, Y, rho_x, rho_y, rho_interp):
+
+    # ----------------
+    # create masks
+    # ----------------
+    # boundary mask
+    dirichlet_label = np.ones_like(X)
+    dirichlet_label[1:-1, 1:-1] = 0
+
+    # integration mask
+    mask = np.ones_like(X)
+
+    # boundary condition
+    rho_dirichlet = rho_interp(X[0, :], Y[:, 0])
+
+    # calculate density 
+    rho, sigma_rho = Density_integration_Poisson_uncertainty(X, Y, mask, rho_x, rho_y, dirichlet_label, rho_dirichlet, uncertainty_quantification=False, sigma_grad_x=0.0, sigma_grad_y=0.0, sigma_dirichlet=1e-12)
+
+    return rho
+
+
+def calculate_density_from_displacements(X_pix, Y_pix, U, V, integration_parameters, rho_interp):
+
+    # convert positions to physical grid
+    X = convert_positions_to_physical_grid(X_pix, integration_parameters['x0'], integration_parameters['pixel_pitch'], integration_parameters['magnification_grad'])
+    Y = convert_positions_to_physical_grid(Y_pix, integration_parameters['y0'], integration_parameters['pixel_pitch'], integration_parameters['magnification_grad'])
+
+    # convert displacements to density gradients
+    rho_x = calculate_density_gradient_from_displacement(U, integration_parameters['pixel_pitch'], integration_parameters['Z_D'], integration_parameters['magnification'], integration_parameters['delta_z'], integration_parameters['gladstone_dale'], integration_parameters['n_0'])
+    rho_y = calculate_density_gradient_from_displacement(V, integration_parameters['pixel_pitch'], integration_parameters['Z_D'], integration_parameters['magnification'], integration_parameters['delta_z'], integration_parameters['gladstone_dale'], integration_parameters['n_0'])
+
+    rho = calculate_density_from_gradient(X, Y, rho_x, rho_y, rho_interp)
+
+    return X, Y, rho_x, rho_y, rho
+
diff --git a/sample_script.py b/sample_script.py
index aa60319..e52c516 100755
--- a/sample_script.py
+++ b/sample_script.py
@@ -15,388 +15,7 @@
 
 import loadmat_functions
 
-def plot_figures(X, Y, rho_x, rho_y, rho, sigma_rho):
-
-    # --------------------------
-    # Plot the results
-    # --------------------------
-    fig = plt.figure(1, figsize=(12,8))
-    plt.figure(1)
-    ax1 = fig.add_subplot(2,2,1)
-    ax2 = fig.add_subplot(2,2,2)
-    ax3 = fig.add_subplot(2,2,3)
-    ax4 = fig.add_subplot(2,2,4)
-
-    # plot x gradient
-    plt.axes(ax1)
-    plt.pcolormesh(X, Y, rho_x)
-    plt.colorbar()
-    ax1.set_aspect('equal')
-    plt.title('rho_x')
-
-    # plot y gradient
-    plt.axes(ax2)
-    plt.pcolormesh(X, Y, rho_y)
-    plt.colorbar()
-    ax2.set_aspect('equal')
-    plt.title('rho_y')
-
-    # plot density
-    plt.axes(ax3)
-    plt.pcolormesh(X, Y, rho)
-    plt.colorbar()
-    ax3.set_aspect('equal')
-    plt.title('rho')
-
-    # plot density uncertainty7
-    plt.axes(ax4)
-    plt.pcolormesh(X, Y, sigma_rho)
-    plt.colorbar()
-    ax4.set_aspect('equal')
-    plt.title('sigma rho')
-
-    # plt.tight_layout()
-
-    return fig
-
-
-def create_border_array(nr, nc, fill_val=1):
-    """
-    Function to create a binary array with non-zero elements along the four boundaries
-
-    Args:
-        nr, nc (int): number of rows and columns of the binary array to be created
-        fill_val (int, optional): non-zero value to be assigned to the border elements. Defaults to 1.
-
-    Returns:
-        border_array: the 2d array with non-zero elements along the borders        
-    """
-
-    # create array with fill val
-    border_array = np.ones((nr, nc)) * fill_val
-    border_array[1:-1, 1:-1] = 0
-
-    return border_array
-
-
-def set_bc(nr, nc, rho_0, sigma_rho_0):
-    """
-    Function to set matrices to enforce dirichlet boundary conditions
-
-    Args:
-        nr, nc (int): number of rows and columns
-        rho_0 (float, scalar): reference density (kg/m^3)
-        sigma_rho_0 (float, scalar): uncertainty in reference density (kg/m^3) 
-
-    Returns:
-        dirichlet_label: binary array specifying regions where reference density will be imposed
-        rho_dirichlet (float 2D): array containing reference density values at the non-zero dirichlet label points
-        sigma_rho_dirichlet (float 2D): array containing reference density UNCERTAINTY values at the non-zero dirichlet label points
-    """
-    
-    # define dirichlet boundary points (minimum one point) - here defined to be all boundaries
-    # THIS NEEDS TO BE MODIFIED FOR EACH EXPERIMENT
-    dirichlet_label = create_border_array(nr, nc, fill_val=1)
-
-    # # sample to set left boundary to be 1s
-    # dirichlet_label = np.zeros((nr, nc))
-    # dirichlet_label[:, 0] = 1
-    # dirichlet_label[10, 0] = 1
-
-    # set density and density uncertainty at these boundary points
-    # HERE SET TO BE THE SAME VALUE FOR ALL BOUNDARY POINTS
-    rho_dirichlet = rho_0 * dirichlet_label
-    sigma_rho_dirichlet = sigma_rho_0 * dirichlet_label
-
-    # convert to boolean array
-    dirichlet_label = dirichlet_label.astype('bool')
-
-    return dirichlet_label, rho_dirichlet, sigma_rho_dirichlet
-
-
-def calculate_gradients_from_displacements(U, experimental_parameters):
-    """
-    Function to calculate density gradients from displacements using the BOS equation
-
-    Inputs:
-        U: pixel displacements
-        experimental_parameters: dictionary of optical layout and other parameters
-    
-    Returns:
-        rho_x: density gradient (kg/m^4)
-    """
-    # extract experimental parameters
-    M, Z_D, l_p, delta_z, gladstone_dale, n_0 = experimental_parameters['magnification'], experimental_parameters['Z_D'], experimental_parameters['pixel_pitch'], experimental_parameters['delta_z'], experimental_parameters['gladstone_dale'], experimental_parameters['n_0']
-
-    # gradient calculation from BOS equation
-    rho_x = U * l_p * 1 / (Z_D * M * delta_z) * 1 / gladstone_dale * n_0
-    
-    return rho_x
-
-
-def convert_displacements_to_physical_units(X_pix, Y_pix, U, V, sigma_U, sigma_V, experimental_parameters, mask):
-    """
-    Function to convert displacements to density gradients and co-ordinates to physical units.
-
-    Args:
-        X_pix, Y_pix (float): co-ordinate grid in pixels        
-        U, V (float): X, Y displacements in pixels
-        sigma_U, sigma_V (float): X, Y displacements uncertainty in pixels
-        experimental_parameters (dict): structure of all the experimental parameters
-        mask (float): binary array of the integration mask
-
-    Returns:
-        X, Y: co-ordinate grid in meters
-        rho_x, rho_y (float): X, Y density gradients (kg/m^4)
-        sigma_rho_x, sigma)rho_y (float): X, Y density gradient uncertainties (kg/m^4)        
-    """
-
-    # calculate co-ordinates in the density gradient plane (m)
-    X = (X_pix - experimental_parameters['x0']) * experimental_parameters['pixel_pitch'] * 1 / experimental_parameters['magnification_grad']
-    Y = (Y_pix - experimental_parameters['y0']) * experimental_parameters['pixel_pitch'] * 1 / experimental_parameters['magnification_grad']
-    
-    # calculate density gradients from displacements (kg/m^4)
-    rho_x = calculate_gradients_from_displacements(U, experimental_parameters)
-    rho_y = calculate_gradients_from_displacements(V, experimental_parameters)
-
-    # calculate density gradient uncertainty from displacement uncertainty (kg/m^4)
-    # sigma_rho_x = calculate_gradients_from_displacements(sigma_U, experimental_parameters)
-    # sigma_rho_y = calculate_gradients_from_displacements(sigma_V, experimental_parameters)
-
-    # calculate density gradient uncertainty from displacement uncertainty (kg/m^4)
-    factor_1 = experimental_parameters['pixel_pitch'] * \
-                  1 / (experimental_parameters['magnification'] * experimental_parameters['gladstone_dale'] *
-                       experimental_parameters['n_0'] * experimental_parameters['Z_D'])
-    
-    factor_2 = experimental_parameters['pixel_pitch'] * \
-             1 / (experimental_parameters['gladstone_dale'] *
-                  experimental_parameters['n_0'] * experimental_parameters['Z_D']) * \
-             1/experimental_parameters['magnification']**2 * experimental_parameters['sigma_M']
-
-    factor_3 = experimental_parameters['pixel_pitch'] * \
-             1 / (experimental_parameters['magnification'] * experimental_parameters['gladstone_dale'] *
-                  experimental_parameters['n_0']) * \
-             1/experimental_parameters['Z_D']**2 * experimental_parameters['sigma_Z_D']
-
-    sigma_rho_x = np.sqrt((sigma_U * factor_1)**2 + (U * factor_2)**2 + (U * factor_3)**2)
-    sigma_rho_y = np.sqrt((sigma_V * factor_1)**2 + (V * factor_2)**2 + (V * factor_3)**2)
-
-    # --------------------------
-    # enforce mask and ensure valid values for uncertainties
-    # --------------------------
-    # set uncertainties in masked regions to be zero (e.g. for image matching)
-    sigma_rho_x[mask == False] = 0.0
-    sigma_rho_y[mask == False] = 0.0
-
-    # set nan values to be zero
-    sigma_rho_x[~np.isfinite(sigma_rho_x)] = 0.0
-    sigma_rho_y[~np.isfinite(sigma_rho_y)] = 0.0
-
-    # set zero values to a small non-zero number (OR WLS will blow up)
-    sigma_rho_x[sigma_rho_x == 0] = 1e-8 
-    sigma_rho_y[sigma_rho_y == 0] = 1e-8
-
-    return X, Y, rho_x, rho_y, sigma_rho_x, sigma_rho_y
-
-
-def create_mask(nr, nc, Eval):
-    """
-    Function to create the integration mask
-
-    Args:
-        nr (int): number of rows
-        nc (int): number of columns
-        Eval (float): binary array of the mask
-
-    Returns:
-        mask (bool): integration mask as a boolean array (True for integration, False otherwise)
-    """
-    
-    # --------------------------
-    # set mask
-    # --------------------------
-    # create binary mask. True for flow, False for blocked region.
-    if len(Eval.shape) == 1:
-        # mask = np.reshape(a=Eval, newshape=(X.shape[1], X.shape[0]))
-        mask = np.reshape(a=Eval, newshape=(nc, nr))
-        mask = np.transpose(mask)
-    else:
-        mask = Eval
-
-    return mask
-
-
-def load_displacements_correlation(filename, displacement_uncertainty_method):
-    """
-    Function to load the displacements from Prana processing
-
-    Args:
-        filename (str): file name
-        displacement_uncertainty_method (str): name of the displacement uncertainty method
-
-    Returns:
-        X_pix, Y_pix: pixel co-ordinate grid
-        U, V: X, Y displacements
-        sigma_U, sigma_V: X, Y, displacement uncertainties
-        Eval: Processing mask (binary array)
-    """
-    
-    # load the data:
-    data = loadmat_functions.loadmat(filename=filename)
-
-    # extract co-ordinates
-    X_pix = data['X'].astype('float')
-    Y_pix = data['Y'].astype('float')
-
-    # extract displacements (sign conversion for the SAMPLE DATASET ONLY)
-    U = data['U']
-    V = data['V']
-
-    # extract displacement uncertainties
-    sigma_U = data['uncertainty2D'][displacement_uncertainty_method + 'x']
-    sigma_V = data['uncertainty2D'][displacement_uncertainty_method + 'y']
-
-    # bias correction for MC
-    if displacement_uncertainty_method == 'MC':
-        sigma_U = np.sqrt(sigma_U ** 2 - data['uncertainty2D']['biasx'] ** 2)
-        sigma_V = np.sqrt(sigma_V ** 2 - data['uncertainty2D']['biasy'] ** 2)
-
-    # extract mask
-    Eval = data['Eval'] + 1
-
-    return X_pix, Y_pix, U, V, sigma_U, sigma_V, Eval
-
-
-def grid_PTV_velocity_uncertainty_2D(Xn, Yn, fluid_mask, xp, yp, up, vp, up_unc, vp_unc, N_avg=3, dist_epsilon=1e-6):
-  # This function maps the velocity and its uncertainty from PTV data (on unstructured particle locations) to
-  # a predefined Cartesian grid using inverse-distance based weighted average.
-  # Inputs:
-  # Xn, Yn: 2d array specifying the Cartesian grid.
-  # fluid_mask: 2d array of a binary mask specifying the flow region. The particle data will be interpolated only onto the grid points with fluid_mask==True
-  # xp, yp: 1d array of particle locations.
-  # up, vp: 1d array of particle velocity values.
-  # up_unc, vp_unc: 1d array of paticle velocity uncertainty (std).
-  # N_avg: the number of particle data points used to determine the velocity for each grid point.
-  # dist_epsilon: the minimum distance allowed for setting up the weights (to avoid signularity).
-  # Outputs: 
-  # Un, Vn: the 3d arrays of velocity on grid.
-  # Un_unc, Vn_unc: 3d arrays of velocity uncertainty on grid.
-  # Cov_u, Cov_v: the sparse matrices of covariance matrices of interpolated velocity.
-
-  Ny,Nx = np.shape(Xn)
-
-  j,i = np.where(fluid_mask==True)
-  Npts = len(j)
-  Np = len(xp)
-
-  # Setup the linear operator to map the p velocity to grid using inverse distance based average
-  Map_M = scysparse.csr_matrix((Npts,Np),dtype=np.float)
-  for ct_pt in range(Npts):
-    x_grid = Xn[j[ct_pt],i[ct_pt]]
-    y_grid = Yn[j[ct_pt],i[ct_pt]]    
-    # find closest particles for linear interpolation.
-    dist = ((xp - x_grid)**2 + (yp - y_grid)**2)**0.5
-    dist[dist<dist_epsilon] = dist_epsilon
-    close_arg = np.argsort(dist)[:N_avg]
-    inv_dist = 1.0 / dist[close_arg]
-    inv_dist_weight = inv_dist / np.sum(inv_dist)
-    for col in range(N_avg):
-      Map_M[ct_pt,close_arg[col]] = inv_dist_weight[col]
-  
-  # Peform the interpolation.
-  Un = np.zeros((Ny,Nx))
-  Vn = np.zeros((Ny,Nx))
-  
-  Un[j,i] = Map_M.dot(up)
-  Vn[j,i] = Map_M.dot(vp)
-  
-  # Perform the ucnertainty propagation.
-  Un_unc = np.zeros((Ny,Nx))
-  Cov_up = scysparse.diags(up_unc**2,format='csr',dtype=np.float)
-  Cov_u = Map_M * Cov_up * Map_M.transpose()
-  Un_unc[j,i] = Cov_u.diagonal()**0.5
-
-  Vn_unc = np.zeros((Ny,Nx))
-  Cov_vp = scysparse.diags(vp_unc**2,format='csr',dtype=np.float)
-  Cov_v = Map_M * Cov_vp * Map_M.transpose()
-  Vn_unc[j,i] = Cov_v.diagonal()**0.5
-
-  return Un, Vn, Un_unc, Vn_unc, Cov_u, Cov_v
-
-
-def load_displacements_tracking(filename, displacement_uncertainty_method='crlb', experimental_parameters):
-    """
-    Function to load displacments from tracking processing and interpolate them to a grid
-
-    Args:
-        filename (str): name of the file containing the displacements
-        displacement_uncertainty_method (str): uncertainty quantification method. defaults to crlb
-        experimental_parameters (dict): dictionary of all experimental parameters
-
-    Returns:
-        X_grid, Y_grid: co-ordinate grid (pix)
-        U_grid, V_grid: X, Y displacements (pix)
-        sigma_U_grid, sigma_V_grid: X, Y displacement uncertainties (pix)
-    """
-    # ----------------------------
-    # load and extract data
-    # ----------------------------
-
-    # load the data:
-    data = loadmat_functions.loadmat(filename)
-
-    # extract track info
-    X_track = data['tracks'][:, 0]
-    Y_track = data['tracks'][:, 2]
-    if validation:
-        U_track = data['tracks'][:, 15]
-        V_track = data['tracks'][:, 16]
-    else:
-        U_track = data['tracks'][:, 13]
-        V_track = data['tracks'][:, 14]
-
-    # extract uncertainties
-    if displacement_uncertainty_method == 'crlb':
-        sigma_U_track = data['uncertainty2D']['U_std']
-        sigma_V_track = data['uncertainty2D']['V_std']
-    else:
-        sigma_U_track = np.zeros_like(a=X_track)
-        sigma_V_track = np.zeros_like(a=X_track)
-
-    # ----------------------------
-    # identify extents
-    # ----------------------------
-    xmin = X_track.min()
-    xmax = X_track.max()
-    ymin = Y_track.min()
-    ymax = Y_track.max()
-
-    # ----------------------------
-    # interpolate tracks onto a regular grid
-    # ----------------------------
-    # calculate expected number of dots along x and y
-    num_dots_x = np.round((xmax - xmin + 1) / experimental_parameters['dot_spacing'])
-    num_dots_y = np.round((ymax - ymin + 1) / experimental_parameters['dot_spacing'])
-
-    # generaet 1D co-ordinate arrays
-    x_grid = np.linspace(start=np.ceil(xmin), stop=np.floor(xmax), num=num_dots_x)
-    y_grid = np.linspace(start=np.ceil(ymin), stop=np.floor(ymax), num=num_dots_y)
-
-    # generate 2D co-ordinate grid
-    [X_grid, Y_grid] = np.meshgrid(x_grid, y_grid)
-
-    # interpolate results onto a regular grid
-    U_grid, V_grid, sigma_U_grid, sigma_V_grid, Cov_u, Cov_v = grid_PTV_velocity_uncertainty_2D(X_grid, Y_grid,
-                    np.ones_like(X_grid), X_track, Y_track, U_track, V_track, sigma_U_track, sigma_V_track, N_avg=8)
-
-    # set nan elements to zero
-    U_grid[np.logical_not(np.isfinite(U_grid))] = 0.0
-    V_grid[np.logical_not(np.isfinite(V_grid))] = 0.0
-    sigma_U_grid[np.logical_not(np.isfinite(sigma_U_grid))] = 0.0
-    sigma_V_grid[np.logical_not(np.isfinite(sigma_V_grid))] = 0.0
-
-    return X_grid, Y_grid, U_grid, V_grid, sigma_U_grid, sigma_V_grid
-
+import helper_functions
 
 def main():
     
@@ -478,10 +97,10 @@ def main():
     # load displacements and uncertainties from file        
     if displacement_estimation_method == 'c':
         # correlation
-        X_pix, Y_pix, U, V, sigma_U, sigma_V, Eval = load_displacements_correlation(filename, displacement_uncertainty_method)    
+        X_pix, Y_pix, U, V, sigma_U, sigma_V, Eval = helper_functions.load_displacements_correlation(filename, displacement_uncertainty_method)    
     elif displacement_estimation_method == 't':
         # tracking
-        X_pix, Y_pix, U, V, sigma_U, sigma_V = load_displacements_tracking(filename, displacement_uncertainty_method, experimental_parameters)    
+        X_pix, Y_pix, U, V, sigma_U, sigma_V = helper_functions.load_displacements_tracking(filename, experimental_parameters['dot_spacing'], displacement_uncertainty_method)    
 
     # account for sign convention
     if dataset_type == 'synthetic':
@@ -490,16 +109,16 @@ def main():
 
     # create mask array (1 for flow, 0 elsewhere) - only implemented for Correlation at the moment
     if displacement_estimation_method == 'c':
-        mask = create_mask(X_pix.shape[0], X_pix.shape[1], Eval)
+        mask = helper_functions.create_mask(X_pix.shape[0], X_pix.shape[1], Eval)
     elif displacement_estimation_method == 't':        
         mask = np.ones_like(a=U)
 
     # convert displacements to density gradients and co-ordinates to physical units 
-    X, Y, rho_x, rho_y, sigma_rho_x, sigma_rho_y = convert_displacements_to_physical_units(X_pix, Y_pix, U, V, sigma_U, sigma_V, experimental_parameters, mask)
+    X, Y, rho_x, rho_y, sigma_rho_x, sigma_rho_y = helper_functions.convert_displacements_to_physical_units(X_pix, Y_pix, U, V, sigma_U, sigma_V, experimental_parameters, mask)
 
     # define dirichlet boundary points (minimum one point) - here defined to be all boundaries
     # This is specific to the current dataset
-    dirichlet_label, rho_dirichlet, sigma_rho_dirichlet = set_bc(X_pix.shape[0], X_pix.shape[1], experimental_parameters['rho_0'], experimental_parameters['sigma_rho_0'])
+    dirichlet_label, rho_dirichlet, sigma_rho_dirichlet = helper_functions.set_bc(X_pix.shape[0], X_pix.shape[1], experimental_parameters['rho_0'], experimental_parameters['sigma_rho_0'])
         
     # calculate density and uncertainty
     if density_integration_method == 'p':
@@ -523,7 +142,7 @@ def main():
                         }, long_field_names=True)
 
     # plot results
-    fig = plot_figures(X, Y, rho_x, rho_y, rho, sigma_rho)
+    fig = helper_functions.plot_figures(X, Y, rho_x, rho_y, rho, sigma_rho)
     
     # save plot to file
     fig.savefig('sample-result.png')