Home Explore Help
Register Sign In
jarjonam
/
tidal_tails
1
0
Fork 0
Files
Tree: 1e44bcb59d
Branches Tags
tidal_tails
/
docs
/
run_simulated_annealing.html

run_simulated_annealing.html 35 KB
History Raw

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739 
  1. <!doctype html>
    2. 

  2. <html lang="en">
    3. 

  3. <head>
    4. 

  4. <meta charset="utf-8">
    5. 

  5. <meta name="viewport" content="width=device-width, initial-scale=1, minimum-scale=1" />
    6. 

  6. <meta name="generator" content="pdoc 0.5.2" />
    7. 

  7. <title>run_simulated_annealing API documentation</title>
    8. 

  8. <meta name="description" content="Simulated annealing algorithm used to match the simulation of the
    9. 

  9. Antennae to the observations by comparing binarized images." />
    10. 

  10. <link href='https://cdnjs.cloudflare.com/ajax/libs/normalize/8.0.0/normalize.min.css' rel='stylesheet'>
    11. 

  11. <link href='https://cdnjs.cloudflare.com/ajax/libs/10up-sanitize.css/8.0.0/sanitize.min.css' rel='stylesheet'>
    12. 

  12. <link href="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.12.0/styles/github.min.css" rel="stylesheet">
    13. 

  13. <style>.flex{display:flex !important}body{line-height:1.5em}#content{padding:20px}#sidebar{padding:30px;overflow:hidden}.http-server-breadcrumbs{font-size:130%;margin:0 0 15px 0}#footer{font-size:.75em;padding:5px 30px;border-top:1px solid #ddd;text-align:right}#footer p{margin:0 0 0 1em;display:inline-block}#footer p:last-child{margin-right:30px}h1,h2,h3,h4,h5{font-weight:300}h1{font-size:2.5em;line-height:1.1em}h2{font-size:1.75em;margin:1em 0 .50em 0}h3{font-size:1.4em;margin:25px 0 10px 0}h4{margin:0;font-size:105%}a{color:#058;text-decoration:none;transition:color .3s ease-in-out}a:hover{color:#e82}.title code{font-weight:bold}h2[id^="header-"]{margin-top:2em}.ident{color:#900}pre code{background:#f8f8f8;font-size:.8em;line-height:1.4em}code{background:#f2f2f1;padding:1px 4px;overflow-wrap:break-word}h1 code{background:transparent}pre{background:#f8f8f8;border:0;border-top:1px solid #ccc;border-bottom:1px solid #ccc;margin:1em 0;padding:1ex}#http-server-module-list{display:flex;flex-flow:column}#http-server-module-list div{display:flex}#http-server-module-list dt{min-width:10%}#http-server-module-list p{margin-top:0}.toc ul,#index{list-style-type:none;margin:0;padding:0}#index code{background:transparent}#index h3{border-bottom:1px solid #ddd}#index ul{padding:0}#index h4{font-weight:bold}#index h4 + ul{margin-bottom:.6em}#index .two-column{column-count:2}dl{margin-bottom:2em}dl dl:last-child{margin-bottom:4em}dd{margin:0 0 1em 3em}#header-classes + dl > dd{margin-bottom:3em}dd dd{margin-left:2em}dd p{margin:10px 0}.name{background:#eee;font-weight:bold;font-size:.85em;padding:5px 10px;display:inline-block;min-width:40%}.name:hover{background:#e0e0e0}.name > span:first-child{white-space:nowrap}.name.class > span:nth-child(2){margin-left:.4em}.name small{font-weight:normal}.inherited{color:#999;border-left:5px solid #eee;padding-left:1em}.inheritance em{font-style:normal;font-weight:bold}.desc h2{font-weight:400;font-size:1.25em}.desc h3{font-size:1em}.desc dt code{background:inherit}.source summary{color:#666;text-align:right;font-weight:400;font-size:.8em;text-transform:uppercase;cursor:pointer}.source pre{max-height:500px;overflow:auto;margin:0}.source pre code{overflow:visible}.hlist{list-style:none}.hlist li{display:inline}.hlist li:after{content:',\2002'}.hlist li:last-child:after{content:none}.hlist .hlist{display:inline;padding-left:1em}img{max-width:100%}.admonition{padding:.1em .5em}.admonition-title{font-weight:bold}.admonition.note,.admonition.info,.admonition.important{background:#aef}.admonition.todo,.admonition.versionadded,.admonition.tip,.admonition.hint{background:#dfd}.admonition.warning,.admonition.versionchanged,.admonition.deprecated{background:#fd4}.admonition.error,.admonition.danger,.admonition.caution{background:lightpink}</style>
    14. 

  14. <style media="screen and (min-width: 700px)">@media screen and (min-width:700px){#sidebar{width:30%}#content{width:70%;max-width:100ch;padding:3em 4em;border-left:1px solid #ddd}pre code{font-size:1em}.item .name{font-size:1em}main{display:flex;flex-direction:row-reverse;justify-content:flex-end}.toc ul ul,#index ul{padding-left:1.5em}.toc > ul > li{margin-top:.5em}}</style>
    15. 

  15. <style media="print">@media print{#sidebar h1{page-break-before:always}.source{display:none}}@media print{*{background:transparent !important;color:#000 !important;box-shadow:none !important;text-shadow:none !important}a[href]:after{content:" (" attr(href) ")";font-size:90%}a[href][title]:after{content:none}abbr[title]:after{content:" (" attr(title) ")"}.ir a:after,a[href^="javascript:"]:after,a[href^="#"]:after{content:""}pre,blockquote{border:1px solid #999;page-break-inside:avoid}thead{display:table-header-group}tr,img{page-break-inside:avoid}img{max-width:100% !important}@page{margin:0.5cm}p,h2,h3{orphans:3;widows:3}h1,h2,h3,h4,h5,h6{page-break-after:avoid}}</style>
    16. 

  16. </head>
    17. 

  17. <body>
    18. 

  18. <main>
    19. 

  19. <article id="content">
    20. 

  20. <header>
    21. 

  21. <h1 class="title"><code>run_simulated_annealing</code> module</h1>
    22. 

  22. </header>
    23. 

  23. <section id="section-intro">
    24. 

  24. <p>Simulated annealing algorithm used to match the simulation of the
    25. 

  25. Antennae to the observations by comparing binarized images.</p>
    26. 

  26. <details class="source">
    27. 

  27. <summary>Source code</summary>
    28. 

  28. <pre><code class="python">&#34;&#34;&#34;Simulated annealing algorithm used to match the simulation of the
    29. 

  29. Antennae to the observations by comparing binarized images.&#34;&#34;&#34;
    30. 

  30. 

  31. import numpy as np
    32. 

  32. import scipy
    33. 

  33. import pickle
    34. 

  34. from scipy import ndimage
    35. 

  35. from fast_histogram import histogram2d
    36. 

  36. from scipy.signal import convolve2d
    37. 

  37. from numba import jit
    38. 

  38. from matplotlib.image import imread
    39. 

  39. import datetime
    40. 

  40. 

  41. from simulation import Simulation, Galaxy
    42. 

  42. 

  43. T = .25 # Initial tempreature
    44. 

  44. STEPS = 0
    45. 

  45. DECAY = .998 # Exponential decay factor
    46. 

  46. 

  47. def simToBitmap(sim, theta, phi, view, scale, x, y, galaxy):
    48. 

  48.     &#34;&#34;&#34;Obtain a bitmap of one galaxy as viewed from a given direction.
    49. 

  49.        The binning has been chosen so that the scale and the offset (x,y)
    50. 

  50.        are expected to be approximately 1 and (0, 0) respectively.
    51. 

  51. 

  52.         Parameters:
    53. 

  53.             sim (Simulation): Simulation to project.
    54. 

  54.             theta (float): polar angle of viewing direction.
    55. 

  55.             phi (float): azimuthal angle of viewing direction.
    56. 

  56.             view (float): rotation angle along viewing direction.
    57. 

  57.             scale (float): scaling factor.
    58. 

  58.             x (float): x offset
    59. 

  59.             y (float): y offset
    60. 

  60.             galaxy (int): galaxy to plot. Either 0, or 1. They are assumed    
    61. 

  61.                 to have the same number of particles.
    62. 

  62.     &#34;&#34;&#34;
    63. 

  63.     # Obtain components in new x&#39;,y&#39; plane
    64. 

  64.     r_vec = (sim.project(theta, phi, view) - [[x+.12,y+1.3]]) * scale
    65. 

  65.     # Select a single galaxy. We match them separately in the algorithm.
    66. 

  66.     if galaxy==0: r_vec = r_vec[2:len(r_vec)//2-1] #omit central masses
    67. 

  67.     if galaxy==1: r_vec = r_vec[len(r_vec)//2-1:]
    68. 

  68.     # Use a fast histogram, much faster than numpy ()
    69. 

  69.     H = histogram2d(r_vec[:,0], r_vec[:,1], 
    70. 

  70.         range=[[-5,5], [-5,5]], bins=(50, 50))
    71. 

  71.     im = np.zeros((50, 50))
    72. 

  72.     H = convolve2d(H, np.ones((2,2)), mode=&#39;same&#39;) # Smooth the image
    73. 

  73.     im[H&gt;=1] = 1 # Binary the image
    74. 

  74.     
    75. 

  75.     return im
    76. 

  76. 

  77. @jit(nopython=True) # Numba annotation
    78. 

  78. def bitmapScoreAlgo(im1, im2):
    79. 

  79.     &#34;&#34;&#34;Computes the f1 score betwen two binarized images
    80. 

  80. 

  81.     Parameters
    82. 

  82.         im1 (arr): nxm ground truth image
    83. 

  83.         im2 (arr): nxm candidate image
    84. 

  84.     
    85. 

  85.     Returns
    86. 

  86.         f1 score
    87. 

  87.     &#34;&#34;&#34;
    88. 

  88.     TP = np.sum((im1==1.0) &amp; (im2==1.0))
    89. 

  89.     TN = np.sum((im1==0.0) &amp; (im2==0.0))
    90. 

  90.     FP = np.sum((im1==0.0) &amp; (im2==1.0))
    91. 

  91.     FN = np.sum((im1==1.0) &amp; (im2==0.0))
    92. 

  92.     if TP==0: return 0
    93. 

  93.     precision = TP/(TP+FP) 
    94. 

  94.     recall = TP/(TP+FN)
    95. 

  95.     return 2*precision*recall / (precision + recall)
    96. 

  96. 

  97. # The matching algorithm attempts to improve the match by shifting the
    98. 

  98. # image by one pixel in each direction. If none improve the score the
    99. 

  99. # f1-score of said local maximum is retuned. To make this highly efficient
    100. 

  100. # (as this is run millions of time) we explicitely write functions to 
    101. 

  101. # shift an image by 1 pixel in each direction, as these can then be compiled
    102. 

  102. # using Numba (jit annotation) to low-level code.
    103. 

  103. # Performance is crucial here and must sadly be prioritized over conciseness
    104. 

  104. @jit(nopython=True)
    105. 

  105. def shiftBottom(im, im2):
    106. 

  106.     &#34;&#34;&#34;Shifts an image by one pixel downwards.
    107. 

  107.     
    108. 

  108.     Parameters:
    109. 

  109.         im (arr): the nxm image to shift by one pixel.
    110. 

  110.         im2 (arr): an nxm array where the new image will be placed.
    111. 

  111. 

  112.     Returns:
    113. 

  113.         A reference to im2
    114. 

  114.     &#34;&#34;&#34;
    115. 

  115.     im2[1:] = im[:-1]
    116. 

  116.     im2[0] = 0
    117. 

  117.     return im2
    118. 

  118. 

  119. @jit(nopython=True)
    120. 

  120. def shiftTop(im, im2):
    121. 

  121.     &#34;&#34;&#34;Shifts an image by one pixel upwards.&#34;&#34;&#34;
    122. 

  122.     im2[:-1] = im[1:]
    123. 

  123.     im2[-1] = 0
    124. 

  124.     return im2
    125. 

  125. 

  126. @jit(nopython=True)
    127. 

  127. def shiftLeft(im, im2):
    128. 

  128.     &#34;&#34;&#34;Shifts an image by one pixel to the left.&#34;&#34;&#34;
    129. 

  129.     im2[:,1:] = im[:,:-1]
    130. 

  130.     im2[:,0] = 0
    131. 

  131.     return im2
    132. 

  132. 

  133. @jit(nopython=True)
    134. 

  134. def shiftRight(im, im2):
    135. 

  135.     &#34;&#34;&#34;Shifts an image by one pixel to the right.&#34;&#34;&#34;
    136. 

  136.     im2[:,:-1] = im[:,1:]
    137. 

  137.     im2[:,-1] = 0
    138. 

  138.     return im2
    139. 

  139. 

  140. @jit
    141. 

  141. def bitmapScore(im1, im2, _prev=None, _bestscore=None, _zeros=None):
    142. 

  142.     &#34;&#34;&#34;Computes the bitmap score between two images. This is the f1-score
    143. 

  143.        but we allow the algorithm to attempt to improve the score by
    144. 

  144.        shifting the images. The algorithm is implemented recursively.
    145. 

  145. 

  146.     Parameters:
    147. 

  147.         im1 (array): The ground truth nxm image.
    148. 

  148.         im2 (array): The candidate nxm imgae.
    149. 

  149.         _prev: Used internally for recursion.
    150. 

  150.         _bestscore: Used internally for recursion.
    151. 

  151.         _zeros: Used internally for recursion.
    152. 

  152.     
    153. 

  153.     Returns:
    154. 

  154.         f1-score for the two images.
    155. 

  155.     &#34;&#34;&#34;
    156. 

  156.     # When the function is called externally, initialize an array of zeros
    157. 

  157.     # and compute the score for no shifting. The zeros array is used for
    158. 

  158.     # performance to only create a new array once.
    159. 

  159.     if _bestscore is None:
    160. 

  160.         _bestscore = bitmapScoreAlgo(im1, im2)
    161. 

  161.         _zeros = np.zeros_like(im2)
    162. 

  162.     # Attempt to improve the score by shifting the image in a direction
    163. 

  163.     # Keeping track of _prev allows to not &#39;go back&#39; and undo and attempt
    164. 

  164.     # to undo a shift needlessly. 
    165. 

  165.     if _prev is not 0: # try left
    166. 

  166.         shifted = shiftLeft(im2, _zeros)
    167. 

  167.         score = bitmapScoreAlgo(im1, shifted)
    168. 

  168.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    169. 

  169.             _prev=1, _bestscore=score, _zeros=_zeros)
    170. 

  170.     if _prev is not 1: # try right
    171. 

  171.         shifted = shiftRight(im2, _zeros)
    172. 

  172.         score = bitmapScoreAlgo(im1, shifted)
    173. 

  173.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    174. 

  174.             _prev=0, _bestscore=score, _zeros=_zeros)
    175. 

  175.     if _prev is not 2: # try top
    176. 

  176.         shifted = shiftTop(im2, _zeros) 
    177. 

  177.         score = bitmapScoreAlgo(im1, shifted)
    178. 

  178.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    179. 

  179.             _prev=3, _bestscore=score, _zeros=_zeros)
    180. 

  180.     if _prev is not 3: # try bottom
    181. 

  181.         shifted = shiftBottom(im2, _zeros)
    182. 

  182.         score = bitmapScoreAlgo(im1, shifted)
    183. 

  183.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    184. 

  184.             _prev=2, _bestscore=score, _zeros=_zeros)
    185. 

  185.     # Return _bestscore if shifting did not improve (local maximum).
    186. 

  186.     return _bestscore
    187. 

  187. 

  188. 

  189. def attemptSimulation(theta1, phi1, theta2, phi2, rmin=1, e=.5, R=2, 
    190. 

  190.         disk1=.75, disk2=.65, mu=1, plot=False, steps=2000):
    191. 

  191.     &#34;&#34;&#34;Runs a simulation with the given parameters and compares it 
    192. 

  192.     to observations of the antennae to return a score out of 2.
    193. 

  193. 

  194.     Parameters:
    195. 

  195.         theta1 (float): polar angle for the spin of the first galaxy.
    196. 

  196.         phi1 (float): azimuthal angle for the spin of the first galaxy.
    197. 

  197.         theta2 (float): polar angle for the spin of the second galaxy.
    198. 

  198.         phi2 (float): azimuthal angle for the spin of the second galaxy.
    199. 

  199.         rmin (float): closest distance of approach of orbit.
    200. 

  200.         e (float): eccentricity of the orbit.
    201. 

  201.         R (float): initial separation
    202. 

  202.         disk1 (float): radius of the uniform disk of the first galaxy
    203. 

  203.         disk2 (float): radius of the uniform disk of the first galaxy
    204. 

  204.         mu (float): ratio of masses of the two galaxies
    205. 

  205.         plot (float): if true the simulation will be saved to data/annealing/
    206. 

  206.         steps (float): number of times the f1 score will be computed, along
    207. 

  207.             random viewing directions per 100 timesteps.
    208. 

  208. 

  209.     Returns:
    210. 

  210.         f1-score: score obtained by the simulation. 
    211. 

  211.     &#34;&#34;&#34;
    212. 

  212. 

  213.     # Initialize the simulation
    214. 

  214.     sim = Simulation(dt=1E-3, soft=0.1, verbose=False, CONFIG=None, method=&#39;bruteForceNumbaOptimized&#39;)
    215. 

  215.     galaxy1 = Galaxy(orientation = (theta1, phi1), centralMass=2/(1+mu),
    216. 

  216.         sim=sim, bulge={&#39;model&#39;:&#39;plummer&#39;, &#39;particles&#39;:0, &#39;totalMass&#39;:0, &#39;l&#39;:0},
    217. 

  217.         disk={&#39;model&#39;:&#39;uniform&#39;, &#39;particles&#39;:2000, &#39;l&#39;:disk1, &#39;totalMass&#39;:0},
    218. 

  218.         halo={&#39;model&#39;:&#39;NFW&#39;, &#39;particles&#39;:0, &#39;rs&#39;:1, &#39;totalMass&#39;:0})
    219. 

  219.     galaxy2 = Galaxy(orientation = (theta2, phi2), centralMass=2*mu/(1+mu),
    220. 

  220.         sim=sim, bulge={&#39;model&#39;:&#39;plummer&#39;, &#39;particles&#39;:0, &#39;totalMass&#39;:0, &#39;l&#39;:0},
    221. 

  221.         disk={&#39;model&#39;:&#39;uniform&#39;, &#39;particles&#39;:2000, &#39;l&#39;:disk2, &#39;totalMass&#39;:0},
    222. 

  222.         halo={&#39;model&#39;:&#39;NFW&#39;, &#39;particles&#39;:0, &#39;rs&#39;:1, &#39;totalMass&#39;:0})
    223. 

  223.     sim.setOrbit(galaxy1, galaxy2, e=e, R0=R, rmin=rmin) # define the orbit
    224. 

  224. 

  225.     # Run the simulation manually
    226. 

  226.     # Initialize the parameters consistently with the velocities
    227. 

  227.     sim.rprev_vec = sim.r_vec - sim.v_vec * sim.dt
    228. 

  228.     # Keep track of the best score
    229. 

  229.     bestScore = 0
    230. 

  230.     # and its corresponding binarized image and parameters
    231. 

  231.     bestImage, bestParams = 0, 0
    232. 

  232.     hasReachedPericenter = False
    233. 

  233. 

  234.     # Run until tmax = 25
    235. 

  235.     for i in range(25001):
    236. 

  236.         sim.step()
    237. 

  237.         if i%100==0: # Every $\Delta t = 0.1$
    238. 

  238.             # Check if we are close to pericenter
    239. 

  239.             centers = sim.r_vec[sim.type[:,0] == &#39;center&#39;]
    240. 

  240.             if np.linalg.norm(centers[0] - centers[1]) &lt; 1.3*rmin: 
    241. 

  241.                 hasReachedPericenter = True
    242. 

  242.             # Do not evaluate the f1-score until pericenter.
    243. 

  243.             if not hasReachedPericenter: continue
    244. 

  244. 

  245.             # Check multiple (steps) viewing directions at random
    246. 

  246.             localBestScore = 0
    247. 

  247.             localBestImage, localBestParams = 0, 0
    248. 

  248.             for j in range(steps):
    249. 

  249.                 # The viewing directions are isotropically distributed
    250. 

  250.                 theta = np.arccos(np.random.uniform(-1, 1))
    251. 

  251.                 phi = np.random.uniform(0, 2*np.pi)
    252. 

  252.                 # Rotation along line of sight
    253. 

  253.                 view = np.random.uniform(0, 2*np.pi)
    254. 

  254.                 scale = np.random.uniform(0.6, 1.0)
    255. 

  255.                 x = np.random.uniform(-1.0, 1.0) # Offsets
    256. 

  256.                 y = np.random.uniform(-1.0, 1.0)
    257. 

  257.                 # Get images for each galaxy and compute their score separately
    258. 

  258.                 im1 = simToBitmap(sim, theta, phi, view, scale, x, y, galaxy=0)
    259. 

  259.                 im2 = simToBitmap(sim, theta, phi, view, scale, x, y, galaxy=1)
    260. 

  260.                 score = bitmapScore(GT1, im1) + bitmapScore(GT2, im2)
    261. 

  261.                 if score &gt; localBestScore: 
    262. 

  262.                     localBestScore = score
    263. 

  263.                     localBestImage = [im1,im2]
    264. 

  264.                     localBestParams = [i, theta, phi, view, scale, x, y]
    265. 

  265. 

  266.             if bestScore &lt; localBestScore: 
    267. 

  267.                 bestScore = localBestScore
    268. 

  268.                 bestImage = localBestImage
    269. 

  269.                 bestParams = localBestParams
    270. 

  270.             if plot:
    271. 

  271.                 sim.save(&#39;annealing&#39;, type=&#39;2D&#39;)
    272. 

  272.     
    273. 

  273.     print(&#39;Best score for this attempt&#39;, bestScore)
    274. 

  274.     print(&#39;using viewing parameters&#39;, bestParams)
    275. 

  275.     
    276. 

  276.     return bestScore
    277. 

  277. 

  278. 

  279. ################################################################################
    280. 

  280. ################################################################################
    281. 

  281. 

  282. # Generate a (50, 50) bitmap for each galaxy
    283. 

  283. # They are stored globally in GT1 and GT2 (Ground Truth)
    284. 

  284. im = imread(&#39;literature/figs/g1c.tif&#39;)
    285. 

  285. im = np.mean(im, axis=-1)
    286. 

  286. im = scipy.misc.imresize(im, (50,50))
    287. 

  287. GT1 = np.zeros((50, 50))
    288. 

  288. GT1[im &gt; 50] = 1
    289. 

  289. 

  290. im = imread(&#39;literature/figs/g2c.tif&#39;)
    291. 

  291. im = np.mean(im, axis=-1)
    292. 

  292. im = scipy.misc.imresize(im, (50,50))
    293. 

  293. GT2 = np.zeros((50, 50))
    294. 

  294. 

  295. # Define the limits and relate scale of the variations in each parameter
    296. 

  296. # In the same order as attemptSimulation
    297. 

  297. # phi1, theta1, phi2, theta2, rmin (fixed), e, R (fixed), disk1, disk2
    298. 

  298. LIMITS = np.array([[np.pi, 2 * np.pi], [-np.pi, np.pi],
    299. 

  299.                    [0, np.pi], [-np.pi, np.pi],
    300. 

  300.                    [1,1], [.5,1.0], [2.2,2.2], [.5,.8], [.5,.8]])
    301. 

  301. VARIATIONS = np.array([.08, .15, .08, .15, 0, .01, 0, .01, .01])
    302. 

  302. 

  303. # Choose a random starting point and evaluate it
    304. 

  304. bestparams = [np.random.uniform(*l) for l in LIMITS]
    305. 

  305. log = []
    306. 

  306. bestscore = attemptSimulation(*bestparams, steps=500)
    307. 

  307. print(&#39;Starting with score&#39;, bestscore, &#39;with parameters&#39;, bestparams)
    308. 

  308. log.append([bestscore, bestparams, True])
    309. 

  309. 

  310. for i in range(STEPS):
    311. 

  311.     T = T * DECAY #exponential decay
    312. 

  312.     # Perturb the parameters
    313. 

  313.     params = bestparams + np.random.normal(scale=VARIATIONS) * 2 * T / 0.04
    314. 

  314.     # Allow the angles from $-\pi$ to $\pi$ to wrap around
    315. 

  315.     for j in [1,3]:
    316. 

  316.         params[j] = np.mod(params[j] - LIMITS[j][0], LIMITS[j][1] - LIMITS[j][0])
    317. 

  317.         params[j] += LIMITS[j][0]
    318. 

  318.     # Clip parameters outside their allowed range
    319. 

  319.     params = np.clip(params, LIMITS[:, 0], LIMITS[:, 1])
    320. 

  320.     # Evaluate the score for this attempt, use more steps as i progresses
    321. 

  321.     # so as to reduce the noise in the evaluation of the score
    322. 

  322.     score = attemptSimulation(*params, steps=500 + i)
    323. 

  323.     # Perform simulated annealing with a typical exponential rule
    324. 

  324.     if score &gt; bestscore or np.exp(-(bestscore-score)/T) &gt; np.random.rand():
    325. 

  325.         # Flip to this new point
    326. 

  326.         print(&#39;NEW BEST ____&#39;, i, T, score, params)
    327. 

  327.         bestscore = score
    328. 

  328.         bestparams = params
    329. 

  329.         log.append([score, params, True])
    330. 

  330.     else: # Remain in the old point
    331. 

  331.         log.append([score, params, False])
    332. 

  332.     # Save the progress for future plotting
    333. 

  333.     pickle.dump(log, open(&#39;data/logs/logb.pickle&#39;, &#34;wb&#34; ) )   </code></pre>
    334. 

  334. </details>
    335. 

  335. </section>
    336. 

  336. <section>
    337. 

  337. </section>
    338. 

  338. <section>
    339. 

  339. </section>
    340. 

  340. <section>
    341. 

  341. <h2 class="section-title" id="header-functions">Functions</h2>
    342. 

  342. <dl>
    343. 

  343. <dt id="run_simulated_annealing.attemptSimulation"><code class="name flex">
    344. 

  344. <span>def <span class="ident">attemptSimulation</span></span>(<span>theta1, phi1, theta2, phi2, rmin=1, e=0.5, R=2, disk1=0.75, disk2=0.65, mu=1, plot=False, steps=2000)</span>
    345. 

  345. </code></dt>
    346. 

  346. <dd>
    347. 

  347. <section class="desc"><p>Runs a simulation with the given parameters and compares it
    348. 

  348. to observations of the antennae to return a score out of 2.</p>
    349. 

  349. <h2 id="parameters">Parameters</h2>
    350. 

  350. <dl>
    351. 

  351. <dt><strong><code>theta1</code></strong> :&ensp;<code>float</code></dt>
    352. 

  352. <dd>polar angle for the spin of the first galaxy.</dd>
    353. 

  353. <dt><strong><code>phi1</code></strong> :&ensp;<code>float</code></dt>
    354. 

  354. <dd>azimuthal angle for the spin of the first galaxy.</dd>
    355. 

  355. <dt><strong><code>theta2</code></strong> :&ensp;<code>float</code></dt>
    356. 

  356. <dd>polar angle for the spin of the second galaxy.</dd>
    357. 

  357. <dt><strong><code>phi2</code></strong> :&ensp;<code>float</code></dt>
    358. 

  358. <dd>azimuthal angle for the spin of the second galaxy.</dd>
    359. 

  359. <dt><strong><code>rmin</code></strong> :&ensp;<code>float</code></dt>
    360. 

  360. <dd>closest distance of approach of orbit.</dd>
    361. 

  361. <dt><strong><code>e</code></strong> :&ensp;<code>float</code></dt>
    362. 

  362. <dd>eccentricity of the orbit.</dd>
    363. 

  363. <dt><strong><code>R</code></strong> :&ensp;<code>float</code></dt>
    364. 

  364. <dd>initial separation</dd>
    365. 

  365. <dt><strong><code>disk1</code></strong> :&ensp;<code>float</code></dt>
    366. 

  366. <dd>radius of the uniform disk of the first galaxy</dd>
    367. 

  367. <dt><strong><code>disk2</code></strong> :&ensp;<code>float</code></dt>
    368. 

  368. <dd>radius of the uniform disk of the first galaxy</dd>
    369. 

  369. <dt><strong><code>mu</code></strong> :&ensp;<code>float</code></dt>
    370. 

  370. <dd>ratio of masses of the two galaxies</dd>
    371. 

  371. <dt><strong><code>plot</code></strong> :&ensp;<code>float</code></dt>
    372. 

  372. <dd>if true the simulation will be saved to data/annealing/</dd>
    373. 

  373. <dt><strong><code>steps</code></strong> :&ensp;<code>float</code></dt>
    374. 

  374. <dd>number of times the f1 score will be computed, along
    375. 

  375. random viewing directions per 100 timesteps.</dd>
    376. 

  376. </dl>
    377. 

  377. <h2 id="returns">Returns</h2>
    378. 

  378. <p>f1-score: score obtained by the simulation.</p></section>
    379. 

  379. <details class="source">
    380. 

  380. <summary>Source code</summary>
    381. 

  381. <pre><code class="python">def attemptSimulation(theta1, phi1, theta2, phi2, rmin=1, e=.5, R=2, 
    382. 

  382.         disk1=.75, disk2=.65, mu=1, plot=False, steps=2000):
    383. 

  383.     &#34;&#34;&#34;Runs a simulation with the given parameters and compares it 
    384. 

  384.     to observations of the antennae to return a score out of 2.
    385. 

  385. 

  386.     Parameters:
    387. 

  387.         theta1 (float): polar angle for the spin of the first galaxy.
    388. 

  388.         phi1 (float): azimuthal angle for the spin of the first galaxy.
    389. 

  389.         theta2 (float): polar angle for the spin of the second galaxy.
    390. 

  390.         phi2 (float): azimuthal angle for the spin of the second galaxy.
    391. 

  391.         rmin (float): closest distance of approach of orbit.
    392. 

  392.         e (float): eccentricity of the orbit.
    393. 

  393.         R (float): initial separation
    394. 

  394.         disk1 (float): radius of the uniform disk of the first galaxy
    395. 

  395.         disk2 (float): radius of the uniform disk of the first galaxy
    396. 

  396.         mu (float): ratio of masses of the two galaxies
    397. 

  397.         plot (float): if true the simulation will be saved to data/annealing/
    398. 

  398.         steps (float): number of times the f1 score will be computed, along
    399. 

  399.             random viewing directions per 100 timesteps.
    400. 

  400. 

  401.     Returns:
    402. 

  402.         f1-score: score obtained by the simulation. 
    403. 

  403.     &#34;&#34;&#34;
    404. 

  404. 

  405.     # Initialize the simulation
    406. 

  406.     sim = Simulation(dt=1E-3, soft=0.1, verbose=False, CONFIG=None, method=&#39;bruteForceNumbaOptimized&#39;)
    407. 

  407.     galaxy1 = Galaxy(orientation = (theta1, phi1), centralMass=2/(1+mu),
    408. 

  408.         sim=sim, bulge={&#39;model&#39;:&#39;plummer&#39;, &#39;particles&#39;:0, &#39;totalMass&#39;:0, &#39;l&#39;:0},
    409. 

  409.         disk={&#39;model&#39;:&#39;uniform&#39;, &#39;particles&#39;:2000, &#39;l&#39;:disk1, &#39;totalMass&#39;:0},
    410. 

  410.         halo={&#39;model&#39;:&#39;NFW&#39;, &#39;particles&#39;:0, &#39;rs&#39;:1, &#39;totalMass&#39;:0})
    411. 

  411.     galaxy2 = Galaxy(orientation = (theta2, phi2), centralMass=2*mu/(1+mu),
    412. 

  412.         sim=sim, bulge={&#39;model&#39;:&#39;plummer&#39;, &#39;particles&#39;:0, &#39;totalMass&#39;:0, &#39;l&#39;:0},
    413. 

  413.         disk={&#39;model&#39;:&#39;uniform&#39;, &#39;particles&#39;:2000, &#39;l&#39;:disk2, &#39;totalMass&#39;:0},
    414. 

  414.         halo={&#39;model&#39;:&#39;NFW&#39;, &#39;particles&#39;:0, &#39;rs&#39;:1, &#39;totalMass&#39;:0})
    415. 

  415.     sim.setOrbit(galaxy1, galaxy2, e=e, R0=R, rmin=rmin) # define the orbit
    416. 

  416. 

  417.     # Run the simulation manually
    418. 

  418.     # Initialize the parameters consistently with the velocities
    419. 

  419.     sim.rprev_vec = sim.r_vec - sim.v_vec * sim.dt
    420. 

  420.     # Keep track of the best score
    421. 

  421.     bestScore = 0
    422. 

  422.     # and its corresponding binarized image and parameters
    423. 

  423.     bestImage, bestParams = 0, 0
    424. 

  424.     hasReachedPericenter = False
    425. 

  425. 

  426.     # Run until tmax = 25
    427. 

  427.     for i in range(25001):
    428. 

  428.         sim.step()
    429. 

  429.         if i%100==0: # Every $\Delta t = 0.1$
    430. 

  430.             # Check if we are close to pericenter
    431. 

  431.             centers = sim.r_vec[sim.type[:,0] == &#39;center&#39;]
    432. 

  432.             if np.linalg.norm(centers[0] - centers[1]) &lt; 1.3*rmin: 
    433. 

  433.                 hasReachedPericenter = True
    434. 

  434.             # Do not evaluate the f1-score until pericenter.
    435. 

  435.             if not hasReachedPericenter: continue
    436. 

  436. 

  437.             # Check multiple (steps) viewing directions at random
    438. 

  438.             localBestScore = 0
    439. 

  439.             localBestImage, localBestParams = 0, 0
    440. 

  440.             for j in range(steps):
    441. 

  441.                 # The viewing directions are isotropically distributed
    442. 

  442.                 theta = np.arccos(np.random.uniform(-1, 1))
    443. 

  443.                 phi = np.random.uniform(0, 2*np.pi)
    444. 

  444.                 # Rotation along line of sight
    445. 

  445.                 view = np.random.uniform(0, 2*np.pi)
    446. 

  446.                 scale = np.random.uniform(0.6, 1.0)
    447. 

  447.                 x = np.random.uniform(-1.0, 1.0) # Offsets
    448. 

  448.                 y = np.random.uniform(-1.0, 1.0)
    449. 

  449.                 # Get images for each galaxy and compute their score separately
    450. 

  450.                 im1 = simToBitmap(sim, theta, phi, view, scale, x, y, galaxy=0)
    451. 

  451.                 im2 = simToBitmap(sim, theta, phi, view, scale, x, y, galaxy=1)
    452. 

  452.                 score = bitmapScore(GT1, im1) + bitmapScore(GT2, im2)
    453. 

  453.                 if score &gt; localBestScore: 
    454. 

  454.                     localBestScore = score
    455. 

  455.                     localBestImage = [im1,im2]
    456. 

  456.                     localBestParams = [i, theta, phi, view, scale, x, y]
    457. 

  457. 

  458.             if bestScore &lt; localBestScore: 
    459. 

  459.                 bestScore = localBestScore
    460. 

  460.                 bestImage = localBestImage
    461. 

  461.                 bestParams = localBestParams
    462. 

  462.             if plot:
    463. 

  463.                 sim.save(&#39;annealing&#39;, type=&#39;2D&#39;)
    464. 

  464.     
    465. 

  465.     print(&#39;Best score for this attempt&#39;, bestScore)
    466. 

  466.     print(&#39;using viewing parameters&#39;, bestParams)
    467. 

  467.     
    468. 

  468.     return bestScore</code></pre>
    469. 

  469. </details>
    470. 

  470. </dd>
    471. 

  471. <dt id="run_simulated_annealing.bitmapScore"><code class="name flex">
    472. 

  472. <span>def <span class="ident">bitmapScore</span></span>(<span>im1, im2)</span>
    473. 

  473. </code></dt>
    474. 

  474. <dd>
    475. 

  475. <section class="desc"><p>Computes the bitmap score between two images. This is the f1-score
    476. 

  476. but we allow the algorithm to attempt to improve the score by
    477. 

  477. shifting the images. The algorithm is implemented recursively.</p>
    478. 

  478. <h2 id="parameters">Parameters</h2>
    479. 

  479. <dl>
    480. 

  480. <dt><strong><code>im1</code></strong> :&ensp;<code>array</code></dt>
    481. 

  481. <dd>The ground truth nxm image.</dd>
    482. 

  482. <dt><strong><code>im2</code></strong> :&ensp;<code>array</code></dt>
    483. 

  483. <dd>The candidate nxm imgae.</dd>
    484. 

  484. <dt><strong><code>_prev</code></strong></dt>
    485. 

  485. <dd>Used internally for recursion.</dd>
    486. 

  486. <dt><strong><code>_bestscore</code></strong></dt>
    487. 

  487. <dd>Used internally for recursion.</dd>
    488. 

  488. <dt><strong><code>_zeros</code></strong></dt>
    489. 

  489. <dd>Used internally for recursion.</dd>
    490. 

  490. </dl>
    491. 

  491. <h2 id="returns">Returns</h2>
    492. 

  492. <p>f1-score for the two images.</p></section>
    493. 

  493. <details class="source">
    494. 

  494. <summary>Source code</summary>
    495. 

  495. <pre><code class="python">@jit
    496. 

  496. def bitmapScore(im1, im2, _prev=None, _bestscore=None, _zeros=None):
    497. 

  497.     &#34;&#34;&#34;Computes the bitmap score between two images. This is the f1-score
    498. 

  498.        but we allow the algorithm to attempt to improve the score by
    499. 

  499.        shifting the images. The algorithm is implemented recursively.
    500. 

  500. 

  501.     Parameters:
    502. 

  502.         im1 (array): The ground truth nxm image.
    503. 

  503.         im2 (array): The candidate nxm imgae.
    504. 

  504.         _prev: Used internally for recursion.
    505. 

  505.         _bestscore: Used internally for recursion.
    506. 

  506.         _zeros: Used internally for recursion.
    507. 

  507.     
    508. 

  508.     Returns:
    509. 

  509.         f1-score for the two images.
    510. 

  510.     &#34;&#34;&#34;
    511. 

  511.     # When the function is called externally, initialize an array of zeros
    512. 

  512.     # and compute the score for no shifting. The zeros array is used for
    513. 

  513.     # performance to only create a new array once.
    514. 

  514.     if _bestscore is None:
    515. 

  515.         _bestscore = bitmapScoreAlgo(im1, im2)
    516. 

  516.         _zeros = np.zeros_like(im2)
    517. 

  517.     # Attempt to improve the score by shifting the image in a direction
    518. 

  518.     # Keeping track of _prev allows to not &#39;go back&#39; and undo and attempt
    519. 

  519.     # to undo a shift needlessly. 
    520. 

  520.     if _prev is not 0: # try left
    521. 

  521.         shifted = shiftLeft(im2, _zeros)
    522. 

  522.         score = bitmapScoreAlgo(im1, shifted)
    523. 

  523.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    524. 

  524.             _prev=1, _bestscore=score, _zeros=_zeros)
    525. 

  525.     if _prev is not 1: # try right
    526. 

  526.         shifted = shiftRight(im2, _zeros)
    527. 

  527.         score = bitmapScoreAlgo(im1, shifted)
    528. 

  528.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    529. 

  529.             _prev=0, _bestscore=score, _zeros=_zeros)
    530. 

  530.     if _prev is not 2: # try top
    531. 

  531.         shifted = shiftTop(im2, _zeros) 
    532. 

  532.         score = bitmapScoreAlgo(im1, shifted)
    533. 

  533.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    534. 

  534.             _prev=3, _bestscore=score, _zeros=_zeros)
    535. 

  535.     if _prev is not 3: # try bottom
    536. 

  536.         shifted = shiftBottom(im2, _zeros)
    537. 

  537.         score = bitmapScoreAlgo(im1, shifted)
    538. 

  538.         if score &gt; _bestscore: return bitmapScore(im1, shifted, 
    539. 

  539.             _prev=2, _bestscore=score, _zeros=_zeros)
    540. 

  540.     # Return _bestscore if shifting did not improve (local maximum).
    541. 

  541.     return _bestscore</code></pre>
    542. 

  542. </details>
    543. 

  543. </dd>
    544. 

  544. <dt id="run_simulated_annealing.bitmapScoreAlgo"><code class="name flex">
    545. 

  545. <span>def <span class="ident">bitmapScoreAlgo</span></span>(<span>im1, im2)</span>
    546. 

  546. </code></dt>
    547. 

  547. <dd>
    548. 

  548. <section class="desc"><p>Computes the f1 score betwen two binarized images</p>
    549. 

  549. <dl>
    550. 

  550. <dt><strong><code>Parameters</code></strong></dt>
    551. 

  551. <dd>im1 (arr): nxm ground truth image
    552. 

  552. im2 (arr): nxm candidate image</dd>
    553. 

  553. <dt><strong><code>Returns</code></strong></dt>
    554. 

  554. <dd>f1 score</dd>
    555. 

  555. </dl></section>
    556. 

  556. <details class="source">
    557. 

  557. <summary>Source code</summary>
    558. 

  558. <pre><code class="python">@jit(nopython=True) # Numba annotation
    559. 

  559. def bitmapScoreAlgo(im1, im2):
    560. 

  560.     &#34;&#34;&#34;Computes the f1 score betwen two binarized images
    561. 

  561. 

  562.     Parameters
    563. 

  563.         im1 (arr): nxm ground truth image
    564. 

  564.         im2 (arr): nxm candidate image
    565. 

  565.     
    566. 

  566.     Returns
    567. 

  567.         f1 score
    568. 

  568.     &#34;&#34;&#34;
    569. 

  569.     TP = np.sum((im1==1.0) &amp; (im2==1.0))
    570. 

  570.     TN = np.sum((im1==0.0) &amp; (im2==0.0))
    571. 

  571.     FP = np.sum((im1==0.0) &amp; (im2==1.0))
    572. 

  572.     FN = np.sum((im1==1.0) &amp; (im2==0.0))
    573. 

  573.     if TP==0: return 0
    574. 

  574.     precision = TP/(TP+FP) 
    575. 

  575.     recall = TP/(TP+FN)
    576. 

  576.     return 2*precision*recall / (precision + recall)</code></pre>
    577. 

  577. </details>
    578. 

  578. </dd>
    579. 

  579. <dt id="run_simulated_annealing.shiftBottom"><code class="name flex">
    580. 

  580. <span>def <span class="ident">shiftBottom</span></span>(<span>im, im2)</span>
    581. 

  581. </code></dt>
    582. 

  582. <dd>
    583. 

  583. <section class="desc"><p>Shifts an image by one pixel downwards.</p>
    584. 

  584. <h2 id="parameters">Parameters</h2>
    585. 

  585. <dl>
    586. 

  586. <dt><strong><code>im</code></strong> :&ensp;<code>arr</code></dt>
    587. 

  587. <dd>the nxm image to shift by one pixel.</dd>
    588. 

  588. <dt><strong><code>im2</code></strong> :&ensp;<code>arr</code></dt>
    589. 

  589. <dd>an nxm array where the new image will be placed.</dd>
    590. 

  590. </dl>
    591. 

  591. <h2 id="returns">Returns</h2>
    592. 

  592. <p>A reference to im2</p></section>
    593. 

  593. <details class="source">
    594. 

  594. <summary>Source code</summary>
    595. 

  595. <pre><code class="python">@jit(nopython=True)
    596. 

  596. def shiftBottom(im, im2):
    597. 

  597.     &#34;&#34;&#34;Shifts an image by one pixel downwards.
    598. 

  598.     
    599. 

  599.     Parameters:
    600. 

  600.         im (arr): the nxm image to shift by one pixel.
    601. 

  601.         im2 (arr): an nxm array where the new image will be placed.
    602. 

  602. 

  603.     Returns:
    604. 

  604.         A reference to im2
    605. 

  605.     &#34;&#34;&#34;
    606. 

  606.     im2[1:] = im[:-1]
    607. 

  607.     im2[0] = 0
    608. 

  608.     return im2</code></pre>
    609. 

  609. </details>
    610. 

  610. </dd>
    611. 

  611. <dt id="run_simulated_annealing.shiftLeft"><code class="name flex">
    612. 

  612. <span>def <span class="ident">shiftLeft</span></span>(<span>im, im2)</span>
    613. 

  613. </code></dt>
    614. 

  614. <dd>
    615. 

  615. <section class="desc"><p>Shifts an image by one pixel to the left.</p></section>
    616. 

  616. <details class="source">
    617. 

  617. <summary>Source code</summary>
    618. 

  618. <pre><code class="python">@jit(nopython=True)
    619. 

  619. def shiftLeft(im, im2):
    620. 

  620.     &#34;&#34;&#34;Shifts an image by one pixel to the left.&#34;&#34;&#34;
    621. 

  621.     im2[:,1:] = im[:,:-1]
    622. 

  622.     im2[:,0] = 0
    623. 

  623.     return im2</code></pre>
    624. 

  624. </details>
    625. 

  625. </dd>
    626. 

  626. <dt id="run_simulated_annealing.shiftRight"><code class="name flex">
    627. 

  627. <span>def <span class="ident">shiftRight</span></span>(<span>im, im2)</span>
    628. 

  628. </code></dt>
    629. 

  629. <dd>
    630. 

  630. <section class="desc"><p>Shifts an image by one pixel to the right.</p></section>
    631. 

  631. <details class="source">
    632. 

  632. <summary>Source code</summary>
    633. 

  633. <pre><code class="python">@jit(nopython=True)
    634. 

  634. def shiftRight(im, im2):
    635. 

  635.     &#34;&#34;&#34;Shifts an image by one pixel to the right.&#34;&#34;&#34;
    636. 

  636.     im2[:,:-1] = im[:,1:]
    637. 

  637.     im2[:,-1] = 0
    638. 

  638.     return im2</code></pre>
    639. 

  639. </details>
    640. 

  640. </dd>
    641. 

  641. <dt id="run_simulated_annealing.shiftTop"><code class="name flex">
    642. 

  642. <span>def <span class="ident">shiftTop</span></span>(<span>im, im2)</span>
    643. 

  643. </code></dt>
    644. 

  644. <dd>
    645. 

  645. <section class="desc"><p>Shifts an image by one pixel upwards.</p></section>
    646. 

  646. <details class="source">
    647. 

  647. <summary>Source code</summary>
    648. 

  648. <pre><code class="python">@jit(nopython=True)
    649. 

  649. def shiftTop(im, im2):
    650. 

  650.     &#34;&#34;&#34;Shifts an image by one pixel upwards.&#34;&#34;&#34;
    651. 

  651.     im2[:-1] = im[1:]
    652. 

  652.     im2[-1] = 0
    653. 

  653.     return im2</code></pre>
    654. 

  654. </details>
    655. 

  655. </dd>
    656. 

  656. <dt id="run_simulated_annealing.simToBitmap"><code class="name flex">
    657. 

  657. <span>def <span class="ident">simToBitmap</span></span>(<span>sim, theta, phi, view, scale, x, y, galaxy)</span>
    658. 

  658. </code></dt>
    659. 

  659. <dd>
    660. 

  660. <section class="desc"><p>Obtain a bitmap of one galaxy as viewed from a given direction.
    661. 

  661. The binning has been chosen so that the scale and the offset (x,y)
    662. 

  662. are expected to be approximately 1 and (0, 0) respectively.</p>
    663. 

  663. <p>Parameters:
    664. 

  664. sim (Simulation): Simulation to project.
    665. 

  665. theta (float): polar angle of viewing direction.
    666. 

  666. phi (float): azimuthal angle of viewing direction.
    667. 

  667. view (float): rotation angle along viewing direction.
    668. 

  668. scale (float): scaling factor.
    669. 

  669. x (float): x offset
    670. 

  670. y (float): y offset
    671. 

  671. galaxy (int): galaxy to plot. Either 0, or 1. They are assumed
    672. 

  672. <br>
    673. 

  673. to have the same number of particles.</p></section>
    674. 

  674. <details class="source">
    675. 

  675. <summary>Source code</summary>
    676. 

  676. <pre><code class="python">def simToBitmap(sim, theta, phi, view, scale, x, y, galaxy):
    677. 

  677.     &#34;&#34;&#34;Obtain a bitmap of one galaxy as viewed from a given direction.
    678. 

  678.        The binning has been chosen so that the scale and the offset (x,y)
    679. 

  679.        are expected to be approximately 1 and (0, 0) respectively.
    680. 

  680. 

  681.         Parameters:
    682. 

  682.             sim (Simulation): Simulation to project.
    683. 

  683.             theta (float): polar angle of viewing direction.
    684. 

  684.             phi (float): azimuthal angle of viewing direction.
    685. 

  685.             view (float): rotation angle along viewing direction.
    686. 

  686.             scale (float): scaling factor.
    687. 

  687.             x (float): x offset
    688. 

  688.             y (float): y offset
    689. 

  689.             galaxy (int): galaxy to plot. Either 0, or 1. They are assumed    
    690. 

  690.                 to have the same number of particles.
    691. 

  691.     &#34;&#34;&#34;
    692. 

  692.     # Obtain components in new x&#39;,y&#39; plane
    693. 

  693.     r_vec = (sim.project(theta, phi, view) - [[x+.12,y+1.3]]) * scale
    694. 

  694.     # Select a single galaxy. We match them separately in the algorithm.
    695. 

  695.     if galaxy==0: r_vec = r_vec[2:len(r_vec)//2-1] #omit central masses
    696. 

  696.     if galaxy==1: r_vec = r_vec[len(r_vec)//2-1:]
    697. 

  697.     # Use a fast histogram, much faster than numpy ()
    698. 

  698.     H = histogram2d(r_vec[:,0], r_vec[:,1], 
    699. 

  699.         range=[[-5,5], [-5,5]], bins=(50, 50))
    700. 

  700.     im = np.zeros((50, 50))
    701. 

  701.     H = convolve2d(H, np.ones((2,2)), mode=&#39;same&#39;) # Smooth the image
    702. 

  702.     im[H&gt;=1] = 1 # Binary the image
    703. 

  703.     
    704. 

  704.     return im</code></pre>
    705. 

  705. </details>
    706. 

  706. </dd>
    707. 

  707. </dl>
    708. 

  708. </section>
    709. 

  709. <section>
    710. 

  710. </section>
    711. 

  711. </article>
    712. 

  712. <nav id="sidebar">
    713. 

  713. <h1>Index</h1>
    714. 

  714. <div class="toc">
    715. 

  715. <ul></ul>
    716. 

  716. </div>
    717. 

  717. <ul id="index">
    718. 

  718. <li><h3><a href="#header-functions">Functions</a></h3>
    719. 

  719. <ul class="two-column">
    720. 

  720. <li><code><a title="run_simulated_annealing.attemptSimulation" href="#run_simulated_annealing.attemptSimulation">attemptSimulation</a></code></li>
    721. 

  721. <li><code><a title="run_simulated_annealing.bitmapScore" href="#run_simulated_annealing.bitmapScore">bitmapScore</a></code></li>
    722. 

  722. <li><code><a title="run_simulated_annealing.bitmapScoreAlgo" href="#run_simulated_annealing.bitmapScoreAlgo">bitmapScoreAlgo</a></code></li>
    723. 

  723. <li><code><a title="run_simulated_annealing.shiftBottom" href="#run_simulated_annealing.shiftBottom">shiftBottom</a></code></li>
    724. 

  724. <li><code><a title="run_simulated_annealing.shiftLeft" href="#run_simulated_annealing.shiftLeft">shiftLeft</a></code></li>
    725. 

  725. <li><code><a title="run_simulated_annealing.shiftRight" href="#run_simulated_annealing.shiftRight">shiftRight</a></code></li>
    726. 

  726. <li><code><a title="run_simulated_annealing.shiftTop" href="#run_simulated_annealing.shiftTop">shiftTop</a></code></li>
    727. 

  727. <li><code><a title="run_simulated_annealing.simToBitmap" href="#run_simulated_annealing.simToBitmap">simToBitmap</a></code></li>
    728. 

  728. </ul>
    729. 

  729. </li>
    730. 

  730. </ul>
    731. 

  731. </nav>
    732. 

  732. </main>
    733. 

  733. <footer id="footer">
    734. 

  734. <p>Generated by <a href="https://pdoc3.github.io/pdoc"><cite>pdoc</cite> 0.5.2</a>.</p>
    735. 

  735. </footer>
    736. 

  736. <script src="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.12.0/highlight.min.js"></script>
    737. 

  737. <script>hljs.initHighlightingOnLoad()</script>
    738. 

  738. </body>
    739. 

  739. </html>
    740.