bd-j / forcepho

Generative modeling galaxy photometry for JWST
https://forcepho.readthedocs.io
19 stars 4 forks source link

Identifying Active sources for a Bulge+Disk galaxy #61

Open Gauri0002 opened 2 years ago

Gauri0002 commented 2 years ago

I am trying to perform Bulge-Disk Decomposition on generated galaxy images and was trying to use the color_fit.py and color_plot_together.py (from demo_color) on them.

I assumed that the galaxy only has 2 components, the bulge and the disk, and both the components have the same position while generating the image. I ran into a problem with forcepho when instead of detecting 2 active sources, only one active source was being detected (either the disk or the bulge).

I traced this issue back to the grow_scene and find_overlaps functions in the superscene.py file. I made a temporary fix in the the find_overlaps function which I have mentioned below,

    def find_overlaps(self, center, radius,seed_index=0, sort=False):
        """Find overlapping sources *not including at the provided center*
        Uses a metric based on |x_i - x_j|/|R_i + R_j|

        Parameters
        ----------
        center : ndarray of shape (2,)
            Floats representing the ra and dec of the center
            *in scene coordinates*

        radius : float
            Radius of influence of the `center`

        sort : bool, optional
            Whether to sort the output in increasing order of the metric

        Returns
        -------
        inds : ndarray of ints
            The indices of the sources that overlap, optionally sorted
            in order of distance from the center
        """
        kinds = self.kdt.query_ball_point(center, radius + self.boundary_radius)
        d = self.scene_coordinates[kinds] - center
        metric = np.hypot(*d.T) / (radius + self.roi[kinds])
        overlaps = (metric < 1) & (metric > 0)

#      LET US CONSIDER ONLY 2 (MAXIMUM) SOURCES
        if not np.any(overlaps)== True:     #IF BOTH VALUES IN OVERLAP ARE FALSE THEN THIS IS FOLLOWED
            if len(kinds)<=1:
                pass                #IF THERE IS ONLY ONE SOURCE THIS PART IS PASSED
            else:
                if seed_index == 0:
                    overlaps = [False,True]     #only defining for 2 sources
                else:
                    overlaps = [True,False]         #only defining for 2 sources

        if sort:
            order = np.argsort(metric[overlaps])
        else:
            order = slice(None)
        return np.array(kinds)[overlaps][order]

The "if- code" I wrote works only for a maximum of 2 sources. I believe it will fail if there are 2 sources at the centre and any other source at another position.

Is there another way of accepting active sources at the centre?

bd-j commented 2 years ago

Hi @Gauri0002, thanks for raising this issue! Indeed the assumption was distinct sources have different positions and this is not good for bulge+disk fits. I think optionally passing the seed index to find_overlaps is a good solution. You might be able to do something a little simpler (and that works for N>2 components), as

if seed_index:
    overlaps = (metric < 1) & (metric >= 0) & (kinds != seed_index)
else:
    overlaps = (metric < 1) & (metric > 0)

then seed_index should default to None or -1, and the calls to find_overlaps in grow_scene and grow_source shoufl pass in seed_index. If you agree I can add this to the code, or I'm happy to take a PR.

bd-j commented 2 years ago

I'd also note that for your use case you can probably get away with

region, active, fixed = sceneDB.checkout_region(seed_index=-1)
active = sceneDB.sourcecat.copy()

that is, just overwrite the checked out active source catalog (which erroneously has only 1 source) with the entire source catalog (which probably only has two sources, the bulge and the disk)

bd-j commented 2 years ago

@Gauri0002 I have pushed what I think is a fix for this in d01fccd but it would be nice to have a unit test.

Gauri0002 commented 2 years ago

@bd-j Thank you for the prompt response! And I apologise for a late one.

I ran the changed superscene.py code for the bulge+disk decomposition fits and it is working just fine!

I have a question regarding line 706 in the superscene.py code,

overlaps = (metric < 1) & (metric >= 0) & (kinds != seed_index)

seed_index by default is set to -1 and since it can only be an integer, the possibility of (kinds != seed_index) being False can occur only for cases in which there is one source right? So unless any user decides to pass a seed_index=0 for the find_overlaps function (seems highly unlikely), overlaps should return a True value and the active source should be registered.

bd-j commented 2 years ago

Great, glad it's working. I'm not sure I totally understand your question; In grow_source we pass the seed index to find_overlaps() so then kinds != seed_index will be an array of booleans that has a value of False for the initial source - the point is to find every overlapping source that is not the input source. However, find_overlaps is also used elsewhere in the code to find all overlapping sources even when the input coordinates do not necessarily correspond to a particular source, in which case we want kinds != seed_index to always be True. hope that helps

bd-j commented 2 years ago

And just to be clear, overlaps is an array and all elements where it is True will be chosen by grow_scene() and grow_source() as sources to include in addition to the initial seed source.

bd-j commented 2 years ago

oh wait, now I see what you mean, that's a bug.

bd-j commented 2 years ago

ok, just pushed a fix.