jarjonam
/
tidal_tails


			
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394
							<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1, minimum-scale=1" />
<meta name="generator" content="pdoc 0.5.2" />
<title>analysis.segmentation API documentation</title>
<meta name="description" content="Segmentation algorithm used to identify the different structures
that are formed in the encounter. This file can be called from the
command line to …" />
<link href='https://cdnjs.cloudflare.com/ajax/libs/normalize/8.0.0/normalize.min.css' rel='stylesheet'>
<link href='https://cdnjs.cloudflare.com/ajax/libs/10up-sanitize.css/8.0.0/sanitize.min.css' rel='stylesheet'>
<link href="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.12.0/styles/github.min.css" rel="stylesheet">
<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>
<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>
<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>
</head>
<body>
<main>
<article id="content">
<header>
<h1 class="title"><code>analysis.segmentation</code> module</h1>
</header>
<section id="section-intro">
<p>Segmentation algorithm used to identify the different structures
that are formed in the encounter. This file can be called from the
command line to make an illustrative plot of the algorithm.</p>
<details class="source">
<summary>Source code</summary>
<pre><code class="python">&#34;&#34;&#34;Segmentation algorithm used to identify the different structures
that are formed in the encounter. This file can be called from the
command line to make an illustrative plot of the algorithm.
&#34;&#34;&#34;

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches

import utils


def segmentEncounter(data, plot=False, mode=&#39;all&#39;):
    &#34;&#34;&#34;Segment the encounter into tail, bridge, orbitting and
    stolen particles as described in the report.

    Parameters:
        data: A data instance as saved by the simulation to a pickle file
        plot: If true the segmentation will be plotted and shown. Highly
            useful for debugging.
        mode (string): If mode is &#39;all&#39; all parts of the encounter will be
            identified. If mode is &#39;bridge&#39; only the bridge will be
            identified. This is useful when there may be no tail.

    Returns:
        masks (tupple): tupple of array corresponding to the masks of the 
            (bridge, stolen, orbitting, tail) particles. One can then use
            e.g. data[&#39;r_vec&#39;][bridgeMask].
        shape (tupple): tupple of (distances, angles) as measured from the
            center of mass and with respect to the x axis. They define the
            shape of the tail
        length (float): total length of the tail.
    &#34;&#34;&#34;
    nRings = 100 # number of rings to use when segmenting the data

    # Localize the central masses
    r_vec = data[&#39;r_vec&#39;]
    centers = r_vec[data[&#39;type&#39;][:,0]==&#39;center&#39;]
    rCenters_vec = centers[1] - centers[0]
    rCenters = np.linalg.norm(rCenters_vec)
    rCenters_unit = rCenters_vec/np.linalg.norm(rCenters_vec)
    # Take particles to be on the tail a priori and
    # remove them as they are found in other structures
    particlesLeft = np.arange(0, len(r_vec))


    if plot:
        colour = &#39;#c40f4c&#39;
        f, axs = plt.subplots(1, 3,  figsize=(9, 4), sharey=False)
        f.subplots_adjust(hspace=0, wspace=0)
        axs[0].scatter(r_vec[:,0], r_vec[:,1], c=colour, alpha=0.1, s=0.1)
        axs[0].axis(&#39;equal&#39;)
        utils.plotCenterMasses(axs[0], data)
        axs[0].axis(&#39;off&#39;)

    # Step 1: project points to see if they are part of the bridge
    parallelProjection = np.dot(r_vec - centers[0], rCenters_unit)
    perpendicularProjection = np.linalg.norm(r_vec - centers[0][np.newaxis] 
        - parallelProjection[:,np.newaxis] * rCenters_unit[np.newaxis], axis=-1)
    bridgeMask = np.logical_and(np.logical_and(0.3*rCenters &lt; parallelProjection,
        parallelProjection &lt; .7*rCenters), perpendicularProjection &lt; 2)

    # Remove the bridge
    notInBridge = np.logical_not(bridgeMask)
    r_vec = r_vec[notInBridge]
    particlesLeft = particlesLeft[notInBridge]

    if mode == &#39;bridge&#39;:
        return (bridgeMask, None, None, None), None, None

    # Step 2: select stolen particles by checking distance to centers
    stolenMask = (np.linalg.norm(r_vec - centers[0][np.newaxis], axis=-1) 
        &gt; np.linalg.norm(r_vec - centers[1][np.newaxis], axis=-1))
    # Remove the stolen part
    notStolen = np.logical_not(stolenMask)
    r_vec = r_vec[notStolen]
    particlesLeft, stolenMask = particlesLeft[notStolen], particlesLeft[stolenMask]

    # Step 3: segment data into concentric rings (spherical shells really)
    r_vec = r_vec - centers[0]
    r = np.linalg.norm(r_vec, axis=-1)
    edges = np.linspace(0, 30, nRings) # nRings concentric spheres 
    indices = np.digitize(r, edges) # Classify particles into shells

    if plot:
        axs[1].scatter(r_vec[:,0], r_vec[:,1], c=colour, alpha=.1, s=.1)
        axs[1].axis(&#39;equal&#39;)
        axs[1].scatter(0, 0, s=100, marker=&#34;*&#34;, c=&#39;black&#39;, alpha=.7)
        axs[1].axis(&#39;off&#39;)
    
    # Step 4: find start of tail
    start = None
    for i in range(1, nRings+1):
        rMean = np.mean(r[indices==i])
        rMean_vec = np.mean(r_vec[indices==i], axis=0)
        parameter = np.linalg.norm(rMean_vec)/rMean

        if plot:
            circ = patches.Circle((0,0), edges[i-1], linewidth=0.5,edgecolor=&#39;black&#39;,facecolor=&#39;none&#39;, alpha=.7)
            axs[1].add_patch(circ)
            txtxy = edges[i-1] * np.array([np.sin(i/13), np.cos(i/13)])
            axs[1].annotate(&#34;{:.2f}&#34;.format(parameter), xy=txtxy, backgroundcolor=&#39;#ffffff55&#39;)
        
        if start is None and parameter&gt;.8 : 
            start = i #Here starts the tail
            startDirection = rMean_vec/np.linalg.norm(rMean_vec)
            if not plot: break;

    if start is None: #abort if nothing found
        raise Exception(&#39;Could not identify tail&#39;)

    # Step 5: remove all circles before start
    inInnerRings = indices &lt; start
    # Remove particles on the opposite direction to startDirection. 
    # in the now innermost 5 rings. Likely traces of the bridge.
    oppositeDirection = np.dot(r_vec, startDirection) &lt; 0
    in5InnermostRings = indices &lt;= start+5
    orbitting = np.logical_or(inInnerRings, 
        np.logical_and(oppositeDirection, in5InnermostRings))
    orbittingMask = particlesLeft[orbitting]
    r_vec = r_vec[np.logical_not(orbitting)]
    tailMask = particlesLeft[np.logical_not(orbitting)]

    if plot:
        axs[2].scatter(r_vec[:,0], r_vec[:,1], c=colour, alpha=0.1, s=0.1)
        axs[2].axis(&#39;equal&#39;)
        axs[2].scatter(0, 0, s=100, marker=&#34;*&#34;, c=&#39;black&#39;, alpha=.7)
        axs[2].axis(&#39;off&#39;)

    # Step 6: measure tail length and shape
    r = np.linalg.norm(r_vec, axis=-1)
    indices = np.digitize(r, edges)
    # Make list of barycenters
    points = [list(np.mean(r_vec[indices==i], axis=0))
        for i in range(1, nRings) if len(r_vec[indices==i])!=0]
    points = np.array(points)
    # Calculate total length
    lengths = np.sqrt(np.sum(np.diff(points, axis=0)**2, axis=1))
    length = np.sum(lengths)
    # Shape (for 2D only)
    angles = np.arctan2(points[:,1], points[:,0])
    distances = np.linalg.norm(points, axis=-1)
    shape = (distances, angles)

    if plot:
        axs[2].plot(points[:,0], points[:,1], c=&#39;black&#39;, linewidth=0.5, marker=&#39;+&#39;)

    if plot:
        plt.show()

    return (bridgeMask, stolenMask, orbittingMask, tailMask), shape, length


if __name__ == &#34;__main__&#34;:
    data = utils.loadData(&#39;200mass&#39;, 10400)
    segmentEncounter(data, plot=True)</code></pre>
</details>
</section>
<section>
</section>
<section>
</section>
<section>
<h2 class="section-title" id="header-functions">Functions</h2>
<dl>
<dt id="analysis.segmentation.segmentEncounter"><code class="name flex">
<span>def <span class="ident">segmentEncounter</span></span>(<span>data, plot=False, mode=&#39;all&#39;)</span>
</code></dt>
<dd>
<section class="desc"><p>Segment the encounter into tail, bridge, orbitting and
stolen particles as described in the report.</p>
<h2 id="parameters">Parameters</h2>
<dl>
<dt><strong><code>data</code></strong></dt>
<dd>A data instance as saved by the simulation to a pickle file</dd>
<dt><strong><code>plot</code></strong></dt>
<dd>If true the segmentation will be plotted and shown. Highly
useful for debugging.</dd>
<dt><strong><code>mode</code></strong> :&ensp;<code>string</code></dt>
<dd>If mode is 'all' all parts of the encounter will be
identified. If mode is 'bridge' only the bridge will be
identified. This is useful when there may be no tail.</dd>
</dl>
<h2 id="returns">Returns</h2>
<dl>
<dt><strong><code>masks</code></strong> :&ensp;<code>tupple</code></dt>
<dd>tupple of array corresponding to the masks of the
(bridge, stolen, orbitting, tail) particles. One can then use
e.g. data['r_vec'][bridgeMask].</dd>
<dt><strong><code>shape</code></strong> :&ensp;<code>tupple</code></dt>
<dd>tupple of (distances, angles) as measured from the
center of mass and with respect to the x axis. They define the
shape of the tail</dd>
<dt><strong><code>length</code></strong> :&ensp;<code>float</code></dt>
<dd>total length of the tail.</dd>
</dl></section>
<details class="source">
<summary>Source code</summary>
<pre><code class="python">def segmentEncounter(data, plot=False, mode=&#39;all&#39;):
    &#34;&#34;&#34;Segment the encounter into tail, bridge, orbitting and
    stolen particles as described in the report.

    Parameters:
        data: A data instance as saved by the simulation to a pickle file
        plot: If true the segmentation will be plotted and shown. Highly
            useful for debugging.
        mode (string): If mode is &#39;all&#39; all parts of the encounter will be
            identified. If mode is &#39;bridge&#39; only the bridge will be
            identified. This is useful when there may be no tail.

    Returns:
        masks (tupple): tupple of array corresponding to the masks of the 
            (bridge, stolen, orbitting, tail) particles. One can then use
            e.g. data[&#39;r_vec&#39;][bridgeMask].
        shape (tupple): tupple of (distances, angles) as measured from the
            center of mass and with respect to the x axis. They define the
            shape of the tail
        length (float): total length of the tail.
    &#34;&#34;&#34;
    nRings = 100 # number of rings to use when segmenting the data

    # Localize the central masses
    r_vec = data[&#39;r_vec&#39;]
    centers = r_vec[data[&#39;type&#39;][:,0]==&#39;center&#39;]
    rCenters_vec = centers[1] - centers[0]
    rCenters = np.linalg.norm(rCenters_vec)
    rCenters_unit = rCenters_vec/np.linalg.norm(rCenters_vec)
    # Take particles to be on the tail a priori and
    # remove them as they are found in other structures
    particlesLeft = np.arange(0, len(r_vec))


    if plot:
        colour = &#39;#c40f4c&#39;
        f, axs = plt.subplots(1, 3,  figsize=(9, 4), sharey=False)
        f.subplots_adjust(hspace=0, wspace=0)
        axs[0].scatter(r_vec[:,0], r_vec[:,1], c=colour, alpha=0.1, s=0.1)
        axs[0].axis(&#39;equal&#39;)
        utils.plotCenterMasses(axs[0], data)
        axs[0].axis(&#39;off&#39;)

    # Step 1: project points to see if they are part of the bridge
    parallelProjection = np.dot(r_vec - centers[0], rCenters_unit)
    perpendicularProjection = np.linalg.norm(r_vec - centers[0][np.newaxis] 
        - parallelProjection[:,np.newaxis] * rCenters_unit[np.newaxis], axis=-1)
    bridgeMask = np.logical_and(np.logical_and(0.3*rCenters &lt; parallelProjection,
        parallelProjection &lt; .7*rCenters), perpendicularProjection &lt; 2)

    # Remove the bridge
    notInBridge = np.logical_not(bridgeMask)
    r_vec = r_vec[notInBridge]
    particlesLeft = particlesLeft[notInBridge]

    if mode == &#39;bridge&#39;:
        return (bridgeMask, None, None, None), None, None

    # Step 2: select stolen particles by checking distance to centers
    stolenMask = (np.linalg.norm(r_vec - centers[0][np.newaxis], axis=-1) 
        &gt; np.linalg.norm(r_vec - centers[1][np.newaxis], axis=-1))
    # Remove the stolen part
    notStolen = np.logical_not(stolenMask)
    r_vec = r_vec[notStolen]
    particlesLeft, stolenMask = particlesLeft[notStolen], particlesLeft[stolenMask]

    # Step 3: segment data into concentric rings (spherical shells really)
    r_vec = r_vec - centers[0]
    r = np.linalg.norm(r_vec, axis=-1)
    edges = np.linspace(0, 30, nRings) # nRings concentric spheres 
    indices = np.digitize(r, edges) # Classify particles into shells

    if plot:
        axs[1].scatter(r_vec[:,0], r_vec[:,1], c=colour, alpha=.1, s=.1)
        axs[1].axis(&#39;equal&#39;)
        axs[1].scatter(0, 0, s=100, marker=&#34;*&#34;, c=&#39;black&#39;, alpha=.7)
        axs[1].axis(&#39;off&#39;)
    
    # Step 4: find start of tail
    start = None
    for i in range(1, nRings+1):
        rMean = np.mean(r[indices==i])
        rMean_vec = np.mean(r_vec[indices==i], axis=0)
        parameter = np.linalg.norm(rMean_vec)/rMean

        if plot:
            circ = patches.Circle((0,0), edges[i-1], linewidth=0.5,edgecolor=&#39;black&#39;,facecolor=&#39;none&#39;, alpha=.7)
            axs[1].add_patch(circ)
            txtxy = edges[i-1] * np.array([np.sin(i/13), np.cos(i/13)])
            axs[1].annotate(&#34;{:.2f}&#34;.format(parameter), xy=txtxy, backgroundcolor=&#39;#ffffff55&#39;)
        
        if start is None and parameter&gt;.8 : 
            start = i #Here starts the tail
            startDirection = rMean_vec/np.linalg.norm(rMean_vec)
            if not plot: break;

    if start is None: #abort if nothing found
        raise Exception(&#39;Could not identify tail&#39;)

    # Step 5: remove all circles before start
    inInnerRings = indices &lt; start
    # Remove particles on the opposite direction to startDirection. 
    # in the now innermost 5 rings. Likely traces of the bridge.
    oppositeDirection = np.dot(r_vec, startDirection) &lt; 0
    in5InnermostRings = indices &lt;= start+5
    orbitting = np.logical_or(inInnerRings, 
        np.logical_and(oppositeDirection, in5InnermostRings))
    orbittingMask = particlesLeft[orbitting]
    r_vec = r_vec[np.logical_not(orbitting)]
    tailMask = particlesLeft[np.logical_not(orbitting)]

    if plot:
        axs[2].scatter(r_vec[:,0], r_vec[:,1], c=colour, alpha=0.1, s=0.1)
        axs[2].axis(&#39;equal&#39;)
        axs[2].scatter(0, 0, s=100, marker=&#34;*&#34;, c=&#39;black&#39;, alpha=.7)
        axs[2].axis(&#39;off&#39;)

    # Step 6: measure tail length and shape
    r = np.linalg.norm(r_vec, axis=-1)
    indices = np.digitize(r, edges)
    # Make list of barycenters
    points = [list(np.mean(r_vec[indices==i], axis=0))
        for i in range(1, nRings) if len(r_vec[indices==i])!=0]
    points = np.array(points)
    # Calculate total length
    lengths = np.sqrt(np.sum(np.diff(points, axis=0)**2, axis=1))
    length = np.sum(lengths)
    # Shape (for 2D only)
    angles = np.arctan2(points[:,1], points[:,0])
    distances = np.linalg.norm(points, axis=-1)
    shape = (distances, angles)

    if plot:
        axs[2].plot(points[:,0], points[:,1], c=&#39;black&#39;, linewidth=0.5, marker=&#39;+&#39;)

    if plot:
        plt.show()

    return (bridgeMask, stolenMask, orbittingMask, tailMask), shape, length</code></pre>
</details>
</dd>
</dl>
</section>
<section>
</section>
</article>
<nav id="sidebar">
<h1>Index</h1>
<div class="toc">
<ul></ul>
</div>
<ul id="index">
<li><h3><a href="#header-functions">Functions</a></h3>
<ul class="">
<li><code><a title="analysis.segmentation.segmentEncounter" href="#analysis.segmentation.segmentEncounter">segmentEncounter</a></code></li>
</ul>
</li>
</ul>
</nav>
</main>
<footer id="footer">
<p>Generated by <a href="https://pdoc3.github.io/pdoc"><cite>pdoc</cite> 0.5.2</a>.</p>
</footer>
<script src="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.12.0/highlight.min.js"></script>
<script>hljs.initHighlightingOnLoad()</script>
</body>
</html>