[Insight-users] Speed optimizations in mutual information registration

Hannu Helminen Hannu . Helminen . 0001 at dosetek . varian . com
Tue, 28 May 2002 08:57:27 +0300


--=====================_53708578==_
Content-Type: text/plain; charset="us-ascii"; format=flowed


[Note: This is a repost with files compressed in ZIP format. My first 
message was rejected because there is a size limit of 40K for messages. 
This is of course good for a mailing list but it also means that this list 
is not very effective for sending code snippets!]

Hi Luis!

Thanks for the explanation. It must be frustrating when just as you have 
managed to reduce the number of template parameters down, a newbie (me) 
comes in and suggests adding new templates!

>As far as I can see, your analysis seems to indicate
>that we may be doing unnecesary constructions of
>itkArray<>'s (which produces the hit in performance).

Yes, this seems to be the case. For example, this innocent-looking statement
       derivative += ( derivB - (*aditer) ) * weight;
involves copying arrays multiple times. E.g. the binary operator- creates 
the result in one copy, and without return value optimization, the compiler 
will generate a second one for the return value. So simply hoisting the 
variable declarations to upper scope is not enough, but one needs to 
carefully eliminate all temporary copies.

A simple approach which would speed up all calculations that use itk::Array 
would be to create a custom memory management system for this class. Some 
C++ experts (e.g Andrei Alexandrescu in his book Modern C++ Design) claim 
that great speed-ups could be achieved with this method. But only 
measurements would tell the truth (I have not done any).

Yet another alternative would be to use templates inside the mutual 
information function. This way we could hard-code a few fast special cases, 
and let a default generic handler do the rest. Since this is a bit hard to 
explain, I will post some example code.

---
Attached please find the code snippets that I have used in my experiments. 
(I have attached the source files in full - about 60K - hopefully this is 
not too much traffic in this mailing list. Diffs would have been much 
smaller, but I do not know if everyone can handle them.)

itkMultiResolutionImageRegistrationMethod-timing.txx shows the timing code 
I have used.

itkMutualInformationImageToImageMetric-simple.txx contains minimal changes 
to the original code that should show the speed-up I have claimed. In my 
test case, the speed spent in the metric function went from approx. 200 
seconds to approx. 40 seconds. However, this case is hard-coded for 7 
dimensions only.

itkMutualInformationImageToImageMetric-generic.txx and .h: This is a 
example of how we could have fast and generic code co-exist in the same 
file. I made a small helper template that does the work of 
GetValueAndDerivative, then instantiated it over itk::Array<double> and 
itk::FixedArray<double, 7>. More explicit instantiations could be added by 
changing the switch statement on line 560.

Of course, these code snippets are intended only to illustrate my point. 
But if any of the ideas proves useful, my changes are in the public domain.

Thanks!

Hannu Helminen

--=====================_53708578==_
Content-Type: application/zip; name="Snippets.zip";
 x-mac-type="705A4950"; x-mac-creator="705A4950"
Content-Transfer-Encoding: base64
Content-Disposition: attachment; filename="Snippets.zip"

UEsDBBQAAAAIAC2IuyyjUSUCJBMAAJBUAAAyAAAAaXRrTXV0dWFsSW5mb3JtYXRpb25JbWFnZVRv
SW1hZ2VNZXRyaWMtZ2VuZXJpYy50eHjtHGtzG7fxu2f8H7ZuR0NKFCll4npKSuxQsh0rtWyPJCdN
Mx3NkQeKZx0P7D1IMa3/e3cXwAH3IEXKctLMlHmIBBaLfWF3AexdZ/f4sT5Pnzx9AvAhljexN+0C
wFmUBDeTFC7FzVREqZcGMoIduBA3QZLG6ueVlOFtkNLIc+lnoaCB8KeL08txQD+C9PY8SzMvPIvG
Mp7yoLOpdyOuJP85F2kcjNrp3V1rDn8iNG+96CbDHkR0urdHLS+9VKFFxOr7NwcH33QOvu0c/gUO
X3Sf/7n7/Lka/IOIE5yhyzSIeaB+HLa/PaB+gjiVs2XMbDVGTUaU83kqo0TGaZBN2zAIQ2CwBGKR
iHgu/DYNvxQCzq7+lmNBylOQMUzSdNbtdBaLRRtZbsv4pvPm6vxtxwJO0imgCMAXqReESVuRg5+r
SZBAIsfpwosF4HefxBsMs1T48OPZ1Zv3H69g8O4n+HFwcTF4d/VTDxZBOpFZCmIuItBo0gkOns7C
AEchptiL0iXIMZy/ujh9g8MGJ2dvz65+ImJfn129e3V5Ca/fX8AAPgwurs5OP74dXBhUHz5efHh/
+aqt2CXE3lDOBYxy2UUyDUYiYYamksi26mXO7je3DT+7nadP/hiMI1+M4Xoza7pGa8JBOCKIxHaD
aK5oFGa+gGcbmu7kWWnQqZx7cYDi/0GMUhlXAXj0hRf5ckoml56lAleTjH9ErZ4ho3elIfMo7OB/
10jCpIrtOy9LEpzubyKORPg6i0asBIJ7+iTypiKZeSNBC/Hpk39TW2cX1bzL1p7GGZFIv1HMqUDz
wfUFRzAKvSSBq9fBnfCZ3JZpOpfzILrhNug/fbKBhI5cNO54HN7tboCg0dSUA0yv32XToYjfjy9n
CO2Flx7SjIZ4DAc9AkhxMe33L0VaD9eA5wfQ7BlkRZkBYvGXKDK0hpGXpEfF7t1+Q62Peol3u+/E
otFsfyfSDzKIUKeNpjOVlcFliqr3Yv8l+iflRJH49rc9BefIZzWgBg0i9NZDbxiEAa516jw4ODjU
/Z2O+j+8lmEoF4gUF2lAsgh+UdgC8m2kcB9SqYFHXjjKsI3BWcm+iIM5NsxFoqZ9mTecamB0Asdg
m41ErpYzoaXCNH2m/z3AzMjO5jLwH8HaPsSomUsRjhtJ6ne7EleA8KY7IJMW0NKLUhQS/Wmis8P1
wYaHbjCbiZgpc1HQIA3NypMJHB3pFvr2rN4Gu/DMAV9p0djHNIrID+vRrzGp8hzrrO/+idbZZHmm
tfZ7/1SlNVVCXvJxhaVWwq2NTTs7dAgcxyKWNYXFREkbEiXuh7rAR7TNNU7r6ZOMcpQIVyqyS1xA
U5tmMIYGNxwfr7alJi70NIsj7RqUk8R8LRhjtsCrE8z6D3EIugOYorOYIloU1eE9nrfRoPn7cNiE
vzIpmHHlng9xYgIV/KLyiAR9BYltzqFROxSFadBWcI01bPTcASebDSibwscoIBWFS616pmtMejEO
T069IGrDK280MTDkDDArS1AcXcYC+5Vxcy/MhNs5kjFSOJORT950yuouw/7WNsfsWYiXzDvGuYIg
MVXAnDUS8Y5ZLuD6RzYp9PWUojGWdbnNkZ2M4kMfFLCBYw0XmwBTWZ9+NAqOrKXNGJ2AbaS9iYyc
wGvGtr+TV/IEeyNl7rbDXXZmvd1jUIZhyq6gXlLdbmDIpy+9jYawSK/zgejJcHVpibfxV8NwNcRd
F3hh+D5LkwDzwGPATE7oTjSGBk/qDB5qzlX7H44Jdw/29vhnk0YBaE3ykr1Bf6nG+uyd71SPFTTr
khSounEqK2mRcm9ObQVj3arBT2OXyGm27Sw/UGcJdxVt7cqCYtK1379CJAktFEW6/DBZJgGmOxxB
GoqNltkB1dDCcIogS4GznrifJTL1ZrgY+DfSbjBOr3MKHGL07BtOF4zJNs8o4s0kJ177/bOE9q++
OMnGY14jzuxNo1vWbpEzh3Qj5jLqVyRM9EwlpD2DqmCBYy9MRN71WX0R2LYNBQcuglzJnzIVk2L2
CzCTSUDOTvXv7RnryMX02SxSkpdDpDV0jRiDJ3aryKRXuOKUZqNWqUemk1guwItA3I3EzE6Ou6pX
puXcG8WycXT0bFBA6SO9JDaL2CDF0Epgrvk+y6X7uRq5vtNJDLp6jEznwkuy+IuCyJdv3hQNZPWP
gA35YztoqNgCH7wYvSUqNqEJdmCW/66EH1L0H9z1Dv/5DzW49qyb8qVXMgeVIQE6ftzIjBAqSJMj
cDiEfrf7DxFLV0FsRVPvFpWNUKydNMc/8RKVB2SYB6AFWPppYMEfYBCy3DYKnNo8aoQbOcybTE6S
oDkMqEvveetDuc2v7kF1sjmqkwIqvQlU3E/ZCNzjIALzZTbEafy38uYymzJitV3tVXuVVZCEV0F8
z+sp72a6twyt3oODsqeiss1aTVzeHtXwwUQMi0ScFJIDjv9DnQBYEJ0CQE/3YRIw1FnAsJQG0B8j
8lxfjKx45NCrgFrlbQSuNVkPq6CZG6/EzcDhxjPceJobz+XGCTt6Xk4QmKVeXY/iIJ8eHHja6mDw
GtZnKfvY5dV3NaFjsN1zFNSrnbW88XUjswPYLFOttVEguxJ0HborfSXC1+3te/Uz30u6hnRptza3
d1yrroKhGZhccRZIm1cBC8bRej1/JmdG64oOB5RbHEAo5cyAFLzX/jH23TQcWk3grrgxF9Th1oHV
dLqA3+eploKsoe7EoU5bcERdwOeY3HDvTiYfiQmOSCYyJHPbP2g/RzEpZLuKpsppo2KB4m6R4b6D
CgNuQWjFvlxTrhBckGq+5qUQYkBOMbSG40LaFkQoDz/w6WYC5AhDLl+SBFE+lqAxxv4iIuzBPHKB
f/x0YgXKSAYKvjavMzYPfn6gFWCElxKSKaaRpdwNCrnDVH0nzbgS2StJb78gDUZnRnaOlUYKjXva
ZJSu8k2vSmU0UM/NJDvwkUJzmkW4BMMlMdrtMi2DOPaWaA4oB5LhLEaKMEml/FjOhLruS8BLGAkP
syPacCmxjcRBYz3UREzpRSwwKR0Juj1s87g3Aok+Q0EnKEs0R501jMUCxtpDUM6Eeqa7sEgIH4U0
XBIQjzcp4iDy7XkzhgEzuA0nGc3Pd1J8/kJJjgwzlYkghoXM0LTEXUqLieZ2mDdYNPnhsm0T66M8
Rr8KxZR02uKzOB5+GfxCKWxJmEcW0kJpYeL2qqGy3I0G7YjDFmwF/01TrZuNwNFg0CmyZdGNXoMY
C3g7hn+OLCT+3NtrOhkzjvo5+CdCikP6uw/iG/zLiD47lpjjZzP8aiLd3Vqk2vtNszAN6Oo0/i2k
tuvM/4ii20YWRoJ7x42vbI9bSEpJB11cvVGJw19LKp1tpGLS22AeJNLYU8XOkPnD9gF0DNy2VqQt
53id6VgBsefjOgP814OJCJEzfSpAk9a61bYeJSBL+LDCcfojL4KUtr0ezEW8hNCLcc89o+oJdJ8I
K+7EiJ0upMFUtBiVjxllKjWsN5UZcoqg+soVwxV67SUFHDlS0UZR8NGZ3WoAd91LSGYYISCbqW22
ubrESTFkeBgZUzXxMCONo1j+laFj1wGGswe6l9DHMPmeG/SFeQK3kVxgjFeBYySns4AyJeSnXRFp
jLFWRt6Q7wqmGDynQSK6sFDCK9M+pjt3HJ/IqXA3+2pmoiqBhi9GgQ1/hCdutgBTEDpqC2HojW4Z
SSpdzWAOQ9A3IqJjDBRKklN7LtHyMe5nMY5f8l5f2wDmWxQVh6TPcRxQaFQdcqxGblAHUUospksj
sdgmFaTbhSLanFfxPNiN3KJhILteoknArOJc8HVdDmxzBC9WCYEXJnygxhnIQrBhCkz8cJS4m8lE
JTAztmnKwghfYoREpxqKAIUrgRku02DUoiIZHI5fg5R4oeIXxmNVpQ8I3aXUVkut6o0KZ3G2tXj+
pQS+gaTf8JQnqFhyBopiupQCvVotVyQQyqN8SbnlQsa3ZHlcx0RU0H3Ntod2RarhHOc6pal69hyG
0Zr2brf+SGPFSYdGMw9iPkaiW64VOZ/y4848xdNC2LG6aimXee/HQeam7ju8X3wAEkutwmMrKZr6
qNsVWW4XDob8hJAROJ9CxwaISqTYT7FjFabCdN3uWTTLUnvfoT7WSPKenloQACVlwe70mqx0pC0H
NrCX3fzsZ9sBJz3DE/NSrtrRFQSVkwoexQf3K7bQDkCkAawN9tyoX3cOlvdZudWfqVgk609fPvce
6nzcZg4ipLvtPRJ0tfvczomtcy9F/7Yx2vwo1BiyZQv4RiDywrwlN1J2NqakSdiFkWhfU2PgsNOq
IoQdfTvRc/B+dSf2Rc5rM6fF/Kw0s6rl5C3qMqUarc6K9WiYiFHd6rLt5jN0BsDRlpIxDrnizsNN
vY9pWETlvS2NTCV5xTxL5zU8NKD0A9M0f8kBMKF65pFSTk6H4GE/UPkmnZ45ZrOjdhKUmjXXLTZX
Coy7JIQjB4Ic6VYTq+SfMmlJ6dQEF0q7vAVaNa/atihn0oc1cxcBd1gjFRqKQKgR8ogNIx1QeqQN
ruqx+5DH9E8b+4PVLqZVxLhROarxMSSJr8DVxpU9W/O2rW8bO4hrvFzgNmmXp6xjtTux3qRRFzVN
gVeeHOz3a4l2KSs6qPpNdG2EpvKXfD9dYEYdzljEzunDV1L45rFtC31/zZjz0HCzUaS54DMvrZu5
Lg0pX8zXE+De0Xc6FiUcl1Ldhi0fe1exDV1ABjVmb1H2cke4rXu3OFq1pqnmtlDt10EYNraTgHNF
XxI02W4DqrOuZJkGPJxZGr2OTerXDB5AmXNUm55dR77fR7XBV6gz2O9vUWjwFUoM7PxaE7NYWGWY
KqqAHJH7yAC7Y3OdauiglJxrtVUZ8FHF5PrOYnXPBqCuXaMfOAGCQavNmxYMO9Kqmc9Vr5/Lt37Z
8PAvWzuMYvUC0iGvpkZiv6+LJFqaTjiuk8nQKYqFFdUULfri11SJdHZ3d+H12d/PX3X5LnRC93se
LDD5w5Td49CmcvY09vj5MVptSSpmf4Vd+uj6AuWLa2N9XVmFLY5sWTHn1684gGltau9hleCUBq7w
7Cdf6NVPNtBUpTYn19SWtTkvRSSnQUSm6BRcrC26cYZ8rwtUV5beVEqAti7TcRj7PZTpGFxo17+v
Mp0i4b9WmU7V/FbV4lSsrghoanKMuvNxKwuBDMTXL9qpMnlP8U6FVyuzjqpXrwuUHCf5RuZEM7KJ
Y6zYuOsYlT9ypn/EMLF5oHj4yi50LQQ9d7yuz1bMVLoe2U38v5pvGzfhag6Ky75Ts7x6xXEmSK3y
FiUcxgiKSFhmBTr2C+ibpUnpXv1BInZ9l13aezS9Xo37TnbSRE4qBnqvR1tbCUggpVpA+H8x4BcW
A3Y6X1YOqB4e4Q1PnEVUP3At4hid8b0oDIY8a2XDZ5VuWUWoBhZqCFVTfQWhY73FQaUO81qC6+Rf
/HjeOp9Sg1wd9NSczbyA1Yd5L4qVMGVs9x/elR6nGUpUvFraVFxhSYFHeMTmkZ7TrD/KM9VMqx+U
aeHc9mxoR+kcG4vy3nEPUQsPvtOJy0SMbulJqYBufbnygtAQW6ufu3HF8MAncVzT/3UeswEoP2Vd
PiE7LqJbcXbYy0lOtIk5pwj5KxMUcXU9TCffsbMAS+uq8OzN/94TQYXzHS9L5fUsjY8e7Ua4Dwpu
gJg/pOropdCiK3I0NQn6/dEEGjU7cmVtVB8FL7rK1nQxz3ERo37rB0AkFo9zTs897hVaC170yVvQ
LMbnD2Ph3WqfO/ayMP2ViSxe8a0kTq9PRdV+37m7odMXNJ5esduY3jHsFCs6yiAnLshJCaT04pZq
ZlqCX/nugrXFHXZ8Kd9Zccxi4ddsCnjwfZsGZ+Z1EfV4kzTe4qr1aCtqV8yo+sjjRBgdU0r3cOue
XDUPL9gRT/8nIqzD371PoLZKKQu4MRSKQdS9nZnrbYK5zFDH5iM5HQYUd+aC319mHW+9+GGt/KGq
gPwAoXA2rytfA5Yl3YBi+EPrplyZXtoxw3jAo3Wtow25OL2anUop1AzqLWaq3gKdAxWYDqnyMBKm
gJUfgI5LI7h4lvQuYy9e5o9qAOXJoToQGQsqWeWCXHq/GT9NInxGkddxMmbzJjWNQkZIAxU9qksI
L1paDhRbVAacZDMqWyYT/d4byWHgRfS+qd/cHmuPfEzCVzrw0cZJX9fmdtXnpp1ZV72/YO17C2bO
/nklPhJj8W0IBa/FRm7exoAS4tczmHkazeqk1ZcRONjsiykUGP/u8bx8/lVI8T6pXcUnVLCD4qUp
LcKevb1Scupg/fkT7TRUIZF6YVgoo5t+w+6K4shvWPYJvpmHUB04S2+LMwlBHTl9tVIdg+hptqDx
UHbKGBVPT3UysSo/Nac+A/2mD1cqLZy/xKPJCfNbtVI9qVl4rDb3B6NRJu8278An/asmI89XsWOc
hOZRc3sl8wLKWyX8WxR+XbHJbUX0zrr8+Zakbu6fH6pS96i1iBs3+UZiqGBq2a3TfM8M12/v+AzE
rY4k+bkTaTCx7xSkf/6I7cGYvv0XUEsDBBQAAAAIAGpfuyxXYwDexQsAAJMoAAAwAAAAaXRrTXV0
dWFsSW5mb3JtYXRpb25JbWFnZVRvSW1hZ2VNZXRyaWMtZ2VuZXJpYy5o7Vrrb9s4Ev8eIP8D0Vvk
7K5rp2l7i3W6Bdyk3Xq3eVzstldggYC26JhbSRQkyo5bdP/2mxlSEinLjwB7384f8hCHv3kPZyj3
Hv/yd30ODw4PGLtO1V3Koz5jbBhn8m6u2UjcRSLWXEsVsyN2I+5kplPz71ip8IvUuPNCBXkocCP7
4eZsNJP4j9RfLnKd83AYz1Qa0aZhxO/EWNGvC6FTOe3OOwv2A4K85/FdDs8B5uzHH/HJOdcGFGDN
3yfHxye94+e9kxfs5Fn/+b/6J8/N5o8izQC/TxKIhTT/PO0+O8Z1pDhTySolpVrTNgGVWp6pOFOp
lnnUZYMwZESWsVRkIl2IoIvbR0Kw4fj3EqWr7zVTKZtrnfR7veVy2QWFuyq9670bX7zvVYRzHTEw
AAuE5jLMukYc+IznMmOZmuklTwWDvwM0rpzkWgTs03D87urDmA0uP7NPg5ubweX48ylbSj1XuWZi
IWJmYfQcNkdJKGEXIKU81iumZuzizc3ZO9g2eD18Pxx/RmHfDseXb0Yj9vbqhg3Y9eBmPDz78H5w
U0Bdf7i5vhq96Rp1EZhP1EKwaWm7WGk5FRkpFCkUu3IuabY72Pb8PO4dHvxDzuJAzNjt7X7BdDuH
PbBBxuJBe5BTPA3zQLBHsKspSh/VaM7UgqcSbP1RTLVK1wmulYw1Pa6jg0r36/S/izQW4ds8npIt
1xlCHqY8PJezmUhFPBUknkd/eBDzSGQJnwrMvsODb/is9/gx+2Ma8ixje5gDIgHIJ6kEs5+pKIFg
zCgSItrr+ptNhF4KCES9VBCBAAGUCp5C5mCdADEDhCPIPVhDmG3lRzgFT85m8h4inscBBOJCxndG
hDUJuqUIlG/GEPCHFpAzHFMNIjwllm8RcWhQVokgbHx+QfjVggspNsiRMUzqTGi2kJxFAvI2yCCt
dMWk1cZdhAPPHSatdtcIGxm7RPwL4OWZwLSGXynLEjGVM8z4MeR7hiYqsYYx6J0o0A3CkiSsaCTB
BGikiCcswRiFZE5VRIoaVQo7Ehq5wrVvoCIuYwPssnKxxYKHOdjW1ibYRljATEDJhfLEdV2PO6Gs
tlYoGa/x7pY297T2xdhg9nKHb3V3K5i9dOwfUEljZOz6ASInjyg8Qf41w8w5+DxMBQ9WEIHChOtU
xTEUCGMWsoYrK5TyTAbkVl0GZ5d54WV0YL8K/RGMikFTpglt25gpyzkcw/Xtgzg4F6lcANnCxTK5
pfR8U66jnSW4JSh3Z3Qa4SkJTtRWPxMzpXcSDg2FAI0zP2emPJzmISGbRJlwDB1lfW5kBvU+SrAU
8f4kwjCzekFaE12SqgmfyBAjKrCRVZ6hJbTINDgIXZBn4C3CuObpVygiSxkHaunI9hqBZjwPNRSY
X3meZVDi2RcqzWWA28gsOFp8cwACxhUsgmsTMltRwc94Pp2vSJMMjupYPCGNod8wiTAFNiYZVmXa
WbYJxARFGRY9Sg6g58w/L4xupQcgsv11cLWxaqXrxbqTQcOatai2ZxxKJWVV1gcfUY0tPChKcZEc
nqSk5J8oJ0uC2T8z+r++i3bgiQadRVkboPsMu64PXhxb3llRFdChhRMEn85RKNfsWHiAGs05qeoA
MgOTXObRRKRXs1EC6vJwZKAbTPNvMIw0XVSDp/FBImKoKzZep3MFTZGlJgDjOtDdNAlT0WVXCW4O
C+KlhEbT4BQw5tww2gw1OqNQY66Wsak6PFW5PZWUBSxYdJzUtSatPFuKLrGvDSmBwxWQZKQb/Ifu
yQAPJYRO3JYEa4FKjWEMVStl4j6BQoCzQdZhSyx+gDCzsqGctLvYhjDH3edsqdIvUDYglamBtAdl
jFKG8qspkoSEpZvHhQRfRapM6mj4ydMADAduJcVg+Wn3uDoZSpa1EGg8gUcW77yAa7UJCHn5R3ID
ZRUuQ4wjMycZibTXaRS1rU/EVNI67LpblTU2HA477FOXtZ7+/PNPBpg9GoTyLkZUNlmxC34vI/m1
1Nkmr9NPPTKiYHWI6QEs/wauwt+q7OVS4I+zUYedPG+dtPtPn/305OmL56Uqf2T1ylI+Lsrh+rIh
AVNBcCbekGgaOwolaOeLlou9NLYZV57o2Cduq/Xq8MA8hbnr9s1/rq9uxnu1kWhnluSTEP5eX37p
sV3jCC2z2do3cxp2z4XzrUuxBYQKBQc1KsWK//eSDcIqnJ26ux4oIRvlkHkkiAcziniqr83p8BKZ
vIJxzvy7mQ66k0wzS41TsK62FMpfmMMY83UKzQ0FoJ6Dp+9Mr6Amf+LxP+M4Bq0Km8D4cSmWF3ya
qhbitx3Emzx+AtXINtlulWphSlB5goSxGdt2IMewwWDuYetOg2ldMRAM20w4NiTyoy64sm7du/gb
xyuHpN8vu0oEY87HWzjdH8fafx2nyZfboX7jUzWBfC1Ea1zYjed2yJ6W9YXdUBeCZ3kq6rZizFnY
jVL1rzUgf2E30HXZnNaA/IXdQFW21oD8hT0MVGW6j1RbeIhIblp7ItXyfX/RfERvIYsbKghdeNBp
R2slnyxPEpXqjYnmWw8jEXDqwVMRlcvNyjQQWkgajAxujahcaYaseaVZQodoh4RezUAwOLDJYi5i
JWG5tBfaVa4b4BzhXLjCc2/iPCqnTO8SAFs/7CMK5wmk/ObinRck7JcGQ/kU3x2eMKnaprsaNAsR
uIaOPzLFouC8UJLGW2estdeZ5nTzk/nIGUk7hqxWT46cB6xtQBzxSvlo0iiEcCtbOanvFGAnPCWN
awg8hXHOCi1BYKYA6J3TrG6Qhnl/t0ms6Rx9jhgBdWp22mGmkVUjppELHZiZoasY6OwNl8xqZBRf
JMTjcvYrrpX86dGOmO4dgDP7dy3GeJsQzLTpUWLwJzR4yFhGJuyfFiDndh41ToFNL449a2+eLVmO
7xhiGlg1CsLaDbG+xUpl7wO0VHhN/8Oa+XV8fm3fI72Cn5fLxa2cXDXMV9ZN5k49s+ZwLicmsGEp
Az33rka2u8aiOHdA3XUL07g4l5DvTUNjJXQplD9FmgGS5s4dkyNjFuBjwboWEeGSrzJ2h70v3RMD
LA6kjmvAtGe4x7pm29jYYYGC8UIUiQZuhDSaQrGWOntp1l71+5cqTpSZyy8kzJqdTYQRv2+1yc9N
QbKPJBuDxL0M/n+M/M0xsuUOogyR/1WE7C/Cxtiwnp0Vb5/KYr6lVFsTbizYpX/tmdB83+Cb9IqG
T6uQT9ip35C6JnjovsMDEFvTRT7dLuwxfrYMv4VM6abmr722sG/fT8tT5ToFq+Hw3Mp00O8rsJXg
0RFTUOaxlYT0kfTLPX4TOpTF3lJWVwBH7VPW6yU5zHQZ3g3GStN7ZbrZojd5RiwFQwHOfL88ZK9z
q2KiaVA76kgFCIey2TPlx7xtMre/nbWyRJXAxhVj9FIGflhibu9JpyrFdxQqDvDM844+r4czlzve
eYqPv+EPR3jmk7Tazqxl+r7j7nG745Te6qnZD5j1Np5Ium9lGLagqh1DzLHvhvivGrdvRa/cOAow
79PE5NTstfm9+VNTat9tdbVpn9vee+rAWaXBvVgyV1tnQkqABb1pf+khvGIbEE/dVCiYO/wUYOL3
F6p3Gwxi8klZwHgCGX8vo6KEOWXLBlzxVqAogqyVBLPytirKNUdTbdCXPtGteTw43VvE11tELK/+
We3y3wvzhwj22grmNZUbPtFtc0tKEbAzbghgW7vyAJgtB1s5KtOI7Bf6ft+5fItu/cWHaCHj6ypY
HNd+iCX6BMoj1Ei8My2mkMb375Mcv+jhjxlEX6l3TsWxtcGVRyX+2nhmY+TMO6SLdmnlTZz4LsdU
VLOp51VyEquEqSbCrBg1G2rUEVufJSsJ/aRdu9Vwr6TOlEqhpgP1jUjwy1r2i2qmEJqCtY3m1OWw
/es1L+tXGJ098NkrR88CyblhKZVqpnLiMbqtSAprK3PFZmrrdzh9Gb5KrH0BqPoKFb5DuRhcfhi8
vx1ejsaDy/FwMB5eXda+ZbTPN/b0/T1+OQnYyRmxKP/6L1BLAwQUAAAACAD9VrYsM3vTe8gHAADm
IAAANAAAAGl0a011bHRpUmVzb2x1dGlvbkltYWdlUmVnaXN0cmF0aW9uTWV0aG9kLXRpbWluZy50
eHjVWW1v28gR/m7A/2HOdwioOJHs9NoCkh1Ap7Mb4WxLkOQe8smgyJXEmuKy3KVkN8h/78wuX5ak
Xsi7oEgdICZXM8/M7sw8Myt33l5/q5/Tk9MTgHHEl5G97gLAMBDeciVhypZrFkhbejyANzBhS0/I
SL/OOPefPUma99yNfUaK8NNkMF149OLJ5/vYl96ECe7HpDJc20tmYtwzueJuW768vNvAT4R0ZwfL
GKUQa3B+Tiu/2lIjI7Z+/nBx8aFz8XPn8u/42P3rz93Lv2jlf7JIIGpXucE2nn65bP+NPiaBAQ9f
I7Uxy2kpnGynAx4IHkkvXreh7/ugxARETLBow9w2qU8Zg+HstwwFHZfAI1hJGXY7ne1228ZNt3m0
7Hya3d91csGVXMMCJV0mbc8Xbe0O/sxWngDBF3JrRwzw2aXD8eaxZC78Ppx9Gj3OoP/wGX7vTyb9
h9nnHmw9PLRYAtuwABIYuULldeh7qIVIkR3IV+ALuL+ZDD6hWv+X4d1w9pmcvR3OHm6mU7gdTaAP
4/5kNhw83vUnKdT4cTIeTW/aersEbM/5hoGTnV3ApecwoTa05uR2gI9rFVG1s+MJV/Pnbef05Edv
EbhsAU+18+kJ8wn1UMkLWGM9shg4fuwyOKufwquzkt6EOTFm44aVAMavWGKeq3BuPV+ySOka2lfS
W7P26iMtBvaaidB2GFXT6ckXWuu8xUi9VQkro9iRPKJ3PCnJMAOwQuAK5GvISBdmt94L09beGav3
fOMFS7UMaKjeJq9MMBMCEbrdehhWK9kFwPoph8O8u4aLHnQ6sLIFSA5zBmGEJlzM6PmrysMYS7Gt
NU3/m2nOsDYEpSs0tTkMMFghxxPGvG/mLcOSdjSHNdMchZgM3n9Y1EiTdFHuUTA4mn1AYKiK1WKj
bKLqqZMNtbRol8OVwJBLlcUZJlm3+8C2VqsHlWilmtdQXTQ10xR5iNdzFo0Wd8h2vkC1y57+YBBH
EXYmta5OJtWYSh7iwsL2BcsWh4EnPdvPgj+20SLD7RNk/kIeWJet3jGl0eKBvWS29+nf2ULWt3jM
Zhvj5Vtw0b5YQCMPdynuca0senrylf7LOCexhwlJGSeYlBhC3XyoOBweBMzBrjlncsuwOzl8HfIA
wyTa8KdJasM999tQVb4Nq4XeR3wL1s2Lw0ICGM3/hXvIaQqrYWqjwis4K+Y8C1rzFhb8UOCvlu6e
X/QvJOsM7952Im5dXZ0Z0tjnsYli+eJ0EcgzHRWAr9pigm5uvwa8KX4MH7QBTUt1sLVkCRb24Oak
VQM6Fz6KrsBz/q4BngvXOvECw9eAL8jXslCl0LqJU6TkxhnUwJzWamgPI1/mZCTlkrn6lKrlD2vs
dKEeI0516f+QGjLS6v3HfzCZdp1cG6VrHB3hwtoTOAs7q4wFPe0ThCka2IELMrV4lvTJZCv0SJzD
ZBwqal1Uw08yO5Lp/cdp7rvumFa1h6Y9YLf6MAhjaRWHswMaj6GLZH5nR0sm5JgL4c19xcc8yLp4
YTfrHdml0avZ2mw/e/SzDRX49IBO8y0pdkwA1bMybGBbe0xhqo1iSf4Vi6eVO5jj5adfDM8fQ8sS
3oICqVYlTZbTJWay5CE/k0Nr6O37j7oaUOCXeLFgEXPT06+YMzv5rtDwtLtotazZKEcHXMjbOHBk
4mTaD3u7hBNLmA9eKl+PzHbNUWpGJf8KX6iMI47XahowvpdBifw0dSwgbNDDkTFs4z2UZu3yJu1I
iv+PbaKne/YJB+4VLo+RGoDu68PAmI2vaYLuVSSShNKnkMqQFH2RUu2fdNsrrV1V2K+XNjBT7vzc
aFb6AatiGGz4M8Peg72LJbGg75BkOxNR30bZvs+3+Aseh3TNpI4vUUt9n0TBzIf6tMNl+lEWB/BV
o9WfSISlWiX7N2TRgmHqgn5vQVq+CmdAkza1dLQXMXXhnVMntUFQFCL27xi5WYsnjV+Fp5UNDl/S
h3nE7Ode+vY1NSKj14oo3XvzINa51aSajs+d5ycJ8hLDpl6sVmYz3X2BqUp6H3bpldPq/BoslHyP
ZlrQSXLLGtyNBr9Nn8Y3k6fpzaBlblWZoFHEgtL15g2wKBtxjSNoemk9olW9T2oFPOnQFgJY6hUl
moOJp4la7V5dy9DL2rFzuUn5KrlqxafI9EQFZqE2i1ahxP8H8VIb7wf0KX0PjXXCHYeYAAfGoHIc
7YKannOUUDaZilziUCoUjgz7dEI+WXM0Y61sjY9E+2C8UX+JFESOIsEg3ydefhsPEV1k80JlVNff
cgfYzTWjaZ1d950qN2PkL3fWWM1bUI1YVA9LSLfbdejvA1dXcFbhkDNarqzimtJjgev39uMUsttA
KqzvwPr6fbpWHFfGEbL9lPmL0+9nMMl8spTfOK9iN1u/AS7eYTt3MfcwY+lXC5tSQB1RZdk0DrFt
+VhzJoTS0uLJjMsFHUmyRGepJ+CuPsB0IG5j+Yy5aoXYpyvxqIBkhZfhZCtNobKMz6CylaZQ5q0l
QzMXmwLmccvg8qWmYEbg89PP1/64b8l1C0ENyR33sT/h8B4TVYmDNnZaKfJp2UKJbesAmoxdhiuy
uQkGsAduP5GXwQ/8AaJO7tboF/Utmk3muPE9nadsbl8rrkQlYV31SM0XP4DSnznp34+47iEX/xdQ
SwMEFAAAAAgA5Xi7LH7A+J/xDgAA0j8AADEAAABpdGtNdXR1YWxJbmZvcm1hdGlvbkltYWdlVG9J
bWFnZU1ldHJpYy1zaW1wbGUudHh47RtrUxtH8rur+A8dknKtkJDAFZ/rJCAlMLblGJsCOTlf6opa
aUdozWpHtw9hfPF/v+6emd3Zh4QgJE6qQh7SznT39Gu6e3pWna39h/rbeLTxCOA0kpeRO+sCwCCM
/ctpAuficibCxE18GcJjOBOXfpxE6nEoZXDlJ4R5Ir00EIQI350dnU98evCTq5M0Sd1gEE5kNGOk
wcy9FEPJHyciifxxO/n0qbWA74jMGze8THEGCR01mzTy3E0UWSSsvj/Z2XnS2fm+s/tP2H3WffqP
7tOnCvknEcW4Qpd5EAtfPey2v9+heYI4kvObiMVyxg0mlMl5JMNYRomfztrQDwJgsBgiEYtoIbw2
oZ8LAYPhjxkV5DwBGcE0SebdTuf6+rqNIrdldNl5NTx508kBp8kMUAXgicT1g7it2MG/4dSPIZaT
5NqNBOB3j9Trj9JEePDzYPjq3fsh9N9+gJ/7Z2f9t8MPPbj2k6lMExALEYImk0wReTYPfMRCSpEb
JjcgJ3ByfHb0CtH6h4M3g+EHYvbFYPj2+PwcXrw7gz6c9s+Gg6P3b/pnhtTp+7PTd+fHbSUuEXZH
ciFgnOkulIk/FjELNJPEdm5elux2d1vzb6uz8ehbfxJ6YgIX63nTBXoTIiGGH4q7IdFa4ThIPQGb
a7rudLOEdCQXbuSj+n8S40RGVQDGPnNDT87I5ZJBInA3yehntOoABf1UQlmEQQf/u0AWplVqL900
jnG5H0UUiuBFGo7ZCAS38Sh0ZyKeu2NBG3Hj0f9orLOFZt5ib0+ilFikZ1RzItB9cH/BHowDN45h
+ML/JDxmt2WGTuTCDy95DA42Hq2hoT2bjI2P6N3uGgSchuYcYHbxNp2NRPRucj5HaDc4d5FndMR9
2OkRQIKbafvgXCT1cA483YFGzxAr6gyQineDKkNvGLtxslec3jpw1P6o13i3+1ZcO432S5GcSj9E
mzoNa6lcB+cJmt6NvOcYn1QQRebb3/cUnKWf5YAa1A8xWo/ckR/4uNdpcmdnZ1fPdzrq//BCBoG8
RqK4SX3Shf9ZUfMptpHBPUikBh67wTjFMQZnI3si8hc4sBCxWvZ5NnCkgTEI7EM+bDQyvJkLrRXm
6Qv97x5uRn62kL73AN52GqFlzkUwceLE63Yl7gDhzh6DjFtAWy9MUEn00cBgh/uDHQ/DYDoXEXNm
kyAkDc3GkzHs7ekR+rZZ74Nd2LTAl3o0zjGPIvSCevIrXKq8xirvu32hVT5ZXmml/96+VGlPlYiX
Ylxhq5Voa2fTwQ4DAuexkHVNaTFW2oZYqfu+IfABfXNF0Np4lFKNEuJORXFJCmho1/Qn4PDA/v5y
X2rgRk/SKNShQQVJrNf8CVYLvDvB7P8AUTAcwAyDxQzJoqp2b4m8jkPrH8BuA35gVrDiyiIf0sQC
yv+s6ogYYwWpbcGpUQcURanfVnDOCjF6NsLheghlV3gf+mSi4EabnvmakF1MwJMz1w/bcOyOpwaG
ggFWZTGqo8tUYLuCt3CDVNiTYxkhh3MZehRNZ2zuMuzX9jkWL4d4zrJjnisoEksFrFlDET022wXs
+MguhbGeSjSmsqq22csXo/xwAArYwLGFi0OApaxHD04hkLW0G2MQyAfpbCJDK/Ea3PZLOZSHOBsq
d88n7G1n9tstDmUEpuoK6jXV7fqGffrSWwuFVXqRIWIkw92lNd7GJ8dINcJTF7hB8C5NYh/rwH3A
Sk7oSXQGhxe1kEdacjX+zT7R7kGzyY8NwgLQluQte4nxUuF6HJ0/qZlc0WxLMqCaxqVyTYuEZzNu
KxTrdg3+OVvETqOdr/ITTZZoV8nW7iwoFl3bB0MkEtNGUazL0+lN7GO5wxnEUWK0zAmohheGUwzl
HFj7iedZIzN3jpuBn5F3Q3F2kXFgMaNXX3M5f0K+OaCMN5dceG0fDGI6v3riMJ1MeI9YqzeMbdm6
Rcks1o2ay6SPSZkYmUpEe4ZUwQMnbhCLbOqL+iJw7C4c7NgEMiN/TFVOijguwFzGPgU7Nd9sGu/I
1PTFbFLSl8Vk7uiaMCZPnFaZSe9wJSmtRqNSYybTSF6DG4L4NBbzfHE8VR2bkRN3HElnb2+zXyDp
Ib+ktpywIYqplcBs993MtPulmrle6iIGQz1mphPhxmn0m5LIbz+8KR7I6x+AGsrHfuCo3AKnboTR
Eg0b0wKPYZ49V9IPGfobe7/Dr7/SgO3PeijbeiV3UBUSYODHg8wYofwk3gNLQjjodv8tImkbiL1o
5l6hsRGKrZNk9KdurOqAFOsA9ICcf0IsxANMQrm0TkHSvI4a40EO6yZTk8ToDn2a0mfe+lSe11e3
kDpcn9RhgZQ+BCrpZ+wEdjuIwDyZjnAZ7428PE9nTFgdV3vVWeUVpOFlEK95P2XTzPcdU6t776Ts
qqycV60mL9+d1OjeTIyKTBwWigPO/yNdAOQgugSAnp7DImCkq4BRqQygD6PyzF5MrNhy6FVAc+Ot
Ba4tWQ+roFkatyRN35LGNdK4WhrXlsZKO3pdLhBYpF7djJIgWx4seDrqYPIa1Vcp2zjl1k81oGOo
3dIK6tWuWj742pnZAmyUudbWKLBdSboW35W5EuOrzva9+pVvZV1D2rznPtfcrzVXwdEMTGa4HEi7
V4EK5tF6O3+hYEb7ipoDKiz2IZBybkAK0Wt7H+cuHYtXk7grYcwGtaS1YDWfNuDrrNRSkDXcHVrc
aQ8OaQq4j8kDt55kMkwscEQ8lQG52/ZO+ymqSRHbUjxVuo1KBMq7RYEPLFKYcAtKK85llrKVYINU
6zU3gQATcoKpNZgUyjY/RH14vkc3EyDHmHL5ksQPM1yCxhz7WYQ4g3XkNX54yTRXKBPpK/jaus74
PHhZQ8vHDC8lxDMsI0u1GxRqh5n6TpaxNdIsaW+7oA0mZzA7+8oihcGmdhllq+zQq0oZDdTLKknU
wntKzEka4gYMbkjMbpc56UeRe4POgFogDc4j5AdLVKqO5Vyoy74Y3JiJMFqO0YZziWOkDMJ10Q4R
FReRwJJ0LOjusM14rwSyPEA1x6hJdEZdM0zENUx0fKCKCa1MN2GhEB6qaHRDQIxvCsR+6OXdZkwC
BrkNhymtzzdS3H2hEkcGqapDkMK1TNGxxKeEthKtbQlvqGj2g5t2XlbvZRn6OBAzsmiLO3GMfu5/
pgK2pMy9HDKH0srEw5Wjaty1kB6L3RbcCf5JQ+2atcDRXTAksl/RfZ5Dgvl8GMOPvRwSH5vNhlUv
I9Yv/n8QUuzS5zaIJ/jJhL5YfpjRZyf83VS6dWeV6tg3S4PEp4vT6Gtobcta/wFVdxddGA02953f
2R/voCmlHQxw9U4ldv8orXTuohVT3PoLP5bGnyp+hsLvtnegY+Du6kXac/ZXuU6uoEoDYSQx76kK
DzOadXP3EE2FB+pM18d6Y8HlrYEWrp3n3ceq2MPBnIgaz681z3iPVToKJZMrG7bg2YGFqvPtQnev
1uwdlJcG+/KUoJ28tW2qN6s70NB1Rk4GDxxB4KyxfOHcPhXjK+qA+eE8TWJOucgMGW95P8U29j07
LMUa6Y9onwCUb89KKuVzikVuid57SjPg1FDA8+ezbHuiUPxSD/7rUimFJSb+i/5FFSXKFXnbY+mp
ZuC1jK64VIVnNt8y5PqDqKneI99tRimWTTNxIaIIj8WbH2SKEQC3QxqLAvpmIw8HnU7RuYB2tVOn
BS3gcscnTIahL9rpdireiNrUgNrAsQ47dS8OKFPWzbBVB+Sb7G5O8fz5p+6LGWJ/ic7YV+6J/d7d
MLTDHMv5zBSm38/tdvvlFs6/5uCvMIkRutTkvaeurPeKm+nAdmrDYCZEGXfF1lpOqGZcs9nPp2IG
rQ6ve0VuKb1mPdtJvMxMNWGFh5bGFkJa3c5rafqwXyfLyLq+hSV9vxZ98Wr6mZ2trS14MfjXyXGX
T+1TOou6cI0HP4zCLlUt7Bd4NoxcftOR9lqciPkPsEV/HRONKRKYMCUs/py6BmB+jdfKNZQ1Cjod
RGFuGzpq8ryeZdchTpxsJldUPsTe3oKMUltfcxe6ETUGY98/HFI2qBrs1lygsG2rLm0537Hj/FyE
cuaH5G5WG3FlK9lCea2vXZc2lCuN7Xs3n/8avWdDC73tr9V7LjL+R/Weq963rMFccboioGk0G3Nn
eEu72wbi9+9EV4W8pSNdkTXXmX4Joy6nckqlgAmHWpB1YmjFx+0YmgetjOU8TOaTeZy0xnSgVCN2
lFxns98lNa2fnO4fIgpT14Leyl81l/eTK1MPHG/+vuu6S7yxLQfF+NGp2ae9Ip5JdsvCTomGcYIi
EdZZgY/tAvlGaVHqO91LxXYQzGNEk5ZXe1Ih63qogZJUHPTW0LjynoxASjdl8PdV2de/KqOv7LRs
jjvejynEwu2YGqq/G7M8r4hUmjA/uLmI/xs5t8SDjHgpFTE1k2hadsPu2EpFbt4JNHmE6K1698tc
/+X4G3+Kjm2eEm9/h6tVPo/YPVko/v7CbmUudCgxLRZ1nB/L2cinDt9C/QKQ5rMmKkefYiO5wMjC
qi308hUDZNVKoWegX+FThQ96pAhj2vHcz6N3v/Hcz9j6dcK8uYnLq9WpyadWUL8DnPINZSRQsBs8
LmHQCYW+A1WvEEYlDG4ykt1l5EY32XUnUKswUNXXRLhJqn5bSL8Q5BtZ4TEJ3McWZfNbRE2CmpDc
oFTNETe8ySVQYtFdbZzOce2EXPS1O5Yj3w3pF1tf3R9r60tzgVCqLrVz0teVdwXVNw+tVZe9Abzy
zd+5lWOX0iM1Ft8nLgQkdnLzPjNqiF9wNus4jeqi1dd5LWr5q90KjJ97vC6XyIVm+kfgu6qPaGCL
xHN/RjsBKyX42GyWkpNF9ZePdA8a00+L9U/uAhleHjh59I0wUObiE3yjmDvoF72F31uaJkUdOwdq
p1oO0dNigXNfccoUlUwbOmMv622byrCv35W3tdLC9Usymm5i1qXLLEqG6nbNxmOz2Q9MRrm8PfwY
Puqnsn/au9hyTiLzYLcouc4LJK+U8q9Q+VXqOF5RvbUvf7kirZuu+H1Nah/HirSxmDAaQwPTyFad
5bNSWb///gVIWp1Jstq09Ktc+udbHPcn9O3/UEsBAhQAFAAAAAgALYi7LKNRJQIkEwAAkFQAADIA
AAAAAAAAAQAgALaBAAAAAGl0a011dHVhbEluZm9ybWF0aW9uSW1hZ2VUb0ltYWdlTWV0cmljLWdl
bmVyaWMudHh4UEsBAhQAFAAAAAgAal+7LFdjAN7FCwAAkygAADAAAAAAAAAAAQAgALaBdBMAAGl0
a011dHVhbEluZm9ybWF0aW9uSW1hZ2VUb0ltYWdlTWV0cmljLWdlbmVyaWMuaFBLAQIUABQAAAAI
AP1Wtiwze9N7yAcAAOYgAAA0AAAAAAAAAAEAIAC2gYcfAABpdGtNdWx0aVJlc29sdXRpb25JbWFn
ZVJlZ2lzdHJhdGlvbk1ldGhvZC10aW1pbmcudHh4UEsBAhQAFAAAAAgA5Xi7LH7A+J/xDgAA0j8A
ADEAAAAAAAAAAQAgALaBoScAAGl0a011dHVhbEluZm9ybWF0aW9uSW1hZ2VUb0ltYWdlTWV0cmlj
LXNpbXBsZS50eHhQSwUGAAAAAAQABAB/AQAA4TYAAAAA
--=====================_53708578==_
Content-Type: text/plain; charset="us-ascii"; format=flowed


--=====================_53708578==_--