text:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Thu Mar  1 11:20:18 2018

@author: dartoon

Propose:
    derive the SB/tot_flux profile when give a fit.
"""

import numpy as np
from regions import PixCoord, CirclePixelRegion
import astropy.io.fits as pyfits
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
from matplotlib.colors import LogNorm
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
import matplotlib as matt
matt.rcParams['font.family'] = 'STIXGeneral'

import copy, matplotlib
my_cmap = copy.copy(matplotlib.cm.get_cmap('gist_heat')) # copy the default cmap
my_cmap.set_bad('black')

def pix_region(center=([49,49]), radius=5):
    '''
    Creat a region file, with pixel units

    Parameter
    --------
        center: The center of the region, with ([reg_x, reg_y]);
        radius: The radius of the region.

    Return
    --------
        A region which is ds9-like.
    '''
    center= PixCoord(x=center[0],y=center[1])
    region = CirclePixelRegion(center, radius)
    return region

def flux_in_region(image,region,mode='exact'):
    '''
    Calculate the total flux inside a given region.

    Parameter
    --------
        image: 2-D array image;
        region: The region generated by pix_region;
        mode: mode type of the mask, 'exact', 'center', default is 'exact'.

    Returns
    --------
        Total flux
    '''
    mask = region.to_mask(mode=mode)
    data = mask.cutout(image)
    tot_flux= np.sum(mask.data * data)
    return tot_flux

def flux_profile(image, center, radius=35,start_p=1.5, grids=20, gridspace=None, ifplot=False,
                 fits_plot=False, mask_list=None, mask_cut=0., msk_image=None,
                 mask_plot = False):
    '''
    Derive the flux profile of one image start at the center.

    Parameters
    --------
        image: A 2-D array image;
        center: The center point of the profile;
        radius: The radius of the profile favourable with default equals to 35;
        grids: The number of points to sample the flux with default equals to 20;
        ifplot: if plot the profile
        fits_plot: if plot the fits file with the regions.
    Returns
    --------
        1. A 1-D array of the tot_flux value of each 'grids' in the profile sampled radius.
        2. The grids of each pixel radius.
        3. The region file for each radius.
    '''
    if gridspace == None:
        r_grids=(np.linspace(0,1,grids+1)*radius)[1:]
        diff = start_p - r_grids[0]
        r_grids += diff             #starts from pixel 0.5
    elif gridspace == 'log':
        r_grids=(np.logspace(-2,0,grids+1)*radius)[1:]
        diff =  start_p - r_grids[0]
        r_grids += diff             #starts from pixel 0.5
    r_flux = np.empty(grids)
    regions = []
    mask = np.ones(image.shape)
    if msk_image is not None:
        mask = msk_image
    elif mask_list != None:
        if msk_image is None:
            for i in range(len(mask_list)):
                mask *= cr_mask(image=image, filename=mask_list[i], mask_reg_cut = mask_cut)
            if mask_plot == True:
                plt.imshow(mask,origin='lower')
                plt.show()
    for i in range(len(r_grids)):
        region = pix_region(center, r_grids[i])
        r_flux[i] =flux_in_region(image*mask, region)
        regions.append(region)

    if fits_plot == True:
        ax=plt.subplot(1,1,1)
        cax=ax.imshow((image*mask),norm=LogNorm(),origin='lower')#,cmap=my_cmap)
        #ax.add_patch(mask.bbox.as_artist(facecolor='none', edgecolor='white'))
        for i in range(grids):
            ax.add_patch(regions[i].as_artist(facecolor='none', edgecolor='orange'))
        plt.colorbar(cax)
        plt.show()
    if ifplot == True:
        minorLocator = AutoMinorLocator()
        fig, ax = plt.subplots()
        plt.plot(r_grids, r_flux, 'x-')
        ax.xaxis.set_minor_locator(minorLocator)
        plt.tick_params(which='both', width=2)
        plt.tick_params(which='major', length=7)
        plt.tick_params(which='minor', length=4, color='r')
        plt.grid()
        ax.set_ylabel("Total Flux")
        ax.set_xlabel("Pixels")
        if gridspace == 'log':
            ax.set_xscale('log')
            plt.xlim(start_p*0.7, )
        plt.grid(which="minor")
        plt.show()
    return r_flux, r_grids, regions

def SB_profile(image, center, radius=35, start_p=1.5, grids=20, gridspace = None,
               ifplot=False, fits_plot=False,
               mask_list=None, mask_plot = False,
               mask_cut = 0., if_annuli= False, msk_image=None):
    '''
    Derive the SB profile of one image start at the center.

    Parameters
    --------
        image: A 2-D array image;
        center: The center point of the profile;
        radius: The radius of the profile favourable with default equals to 35;
        grids: The number of points to sample the flux with default equals to 20;
        ifplot: if plot the profile
        fits_plot: if plot the fits file with the regions.
        msk_image: if is not None, will use this image as mask.
    Returns
    --------
        A 1-D array of the SB value of each 'grids' in the profile with the sampled radius.
    '''
    if mask_list is None and msk_image is None:
        r_flux, r_grids, regions=flux_profile(image, center, radius=radius,start_p=start_p, grids=grids, gridspace=gridspace, ifplot=False, fits_plot=False)
        region_area = np.zeros([len(r_flux)])
        for i in range(len(r_flux)):
            circle=regions[i].to_mask(mode='exact')
            edge_mask = circle.cutout(np.ones(image.shape))
            region_area[i]=(circle.data * edge_mask).sum()
    else:
        mask = np.ones(image.shape)
        if mask_list is not None and msk_image is None:
            for i in range(len(mask_list)):
                mask *= cr_mask(image=image, filename=mask_list[i], mask_reg_cut = mask_cut)
        elif msk_image is not None:
            mask = msk_image
#        plt.imshow(mask, origin='low')
#        plt.show()
        r_flux, r_grids, regions=flux_profile(image*mask, center, radius=radius,start_p=start_p, grids=grids, gridspace=gridspace, ifplot=False, fits_plot=False)
        region_area = np.zeros([len(r_flux)])
        for i in range(len(r_flux)):
            circle=regions[i].to_mask(mode='exact')
            circle_mask =  circle.cutout(mask)
            if i ==len(r_flux)-1 and mask_plot == True:
#                print "plt circle_mask"
                plt.imshow((circle_mask),origin='lower')
                plt.show()
#                plt.imshow((circle.data),origin='lower') #The circle data is zero at outer circle area is zero
#                plt.show()
            region_area[i]=(circle.data * circle_mask).sum()
    if if_annuli ==False:
        r_SB= r_flux/region_area
    elif if_annuli == True:
        r_SB = np.zeros_like(r_flux)
        r_SB[0] = r_flux[0]/region_area[0]
        r_SB[1:] = (r_flux[1:]-r_flux[:-1]) / (region_area[1:]-region_area[:-1])
    if fits_plot == True:
        ax=plt.subplot(1,1,1)
        if mask_list is None and msk_image is None:
            cax=ax.imshow(image, norm=LogNorm(),origin='lower')
        else:
            cax=ax.imshow(image*mask,norm=LogNorm(),origin='lower')
        #ax.add_patch(mask.bbox.as_artist(facecolor='none', edgecolor='white'))
        for i in range(grids):
            ax.add_patch(regions[i].as_artist(facecolor='none', edgecolor='orange'))
        plt.colorbar(cax)
        plt.show()
    if ifplot == True:
        minorLocator = AutoMinorLocator()
        fig, ax = plt.subplots()
        plt.plot(r_grids, r_SB, 'x-')
        ax.xaxis.set_minor_locator(minorLocator)
        plt.tick_params(which='both', width=2)
        plt.tick_params(which='major', length=7)
        plt.tick_params(which='minor', length=4, color='r')
        plt.grid()
        ax.set_ylabel("Surface Brightness")
        ax.set_xlabel("Pixels")
        if gridspace == 'log':
            ax.set_xscale('log')
            plt.xlim(start_p*0.7, )
        plt.grid(which="minor")
        plt.show()
    return r_SB, r_grids

def concentration_profile(image, center, total_flux=None, radius=35,start_p=1.5,
                          grids=100, gridspace=None, ifplot=False,
                          fits_plot=False):
    '''
    Derive the flux profile of one image start at the center.

    Parameters
    --------
        image: A 2-D array image;
        center: The center point of the profile;
        radius: The radius of the profile favourable with default equals to 35;
        grids: The number of points to sample the flux with default equals to 20;
        ifplot: if plot the profile
        fits_plot: if plot the fits file with the regions.
    Returns
    --------
        1. A 1-D array of the tot_flux value of each 'grids' in the profile sampled radius.
        2. The grids of each pixel radius.
        3. The region file for each radius.
    '''
    if total_flux==None:
        total_flux = image.sum()
    if gridspace == None:
        r_grids=(np.linspace(0,1,grids+1)*radius)[1:]
        diff = start_p - r_grids[0]
        r_grids += diff             #starts from pixel 0.5
    elif gridspace == 'log':
        r_grids=(np.logspace(-2,0,grids+1)*radius)[1:]
        diff =  start_p - r_grids[0]
        r_grids += diff             #starts from pixel 0.5
    r80_flux = 0.8* total_flux
    r20_flux = 0.2* total_flux
    r100_flux = total_flux

    r_flux = np.empty(grids)
    regions = []
    for i in range(len(r_grids)):
        region = pix_region(center, r_grids[i])
        r_flux[i] =flux_in_region(image, region)
        regions.append(region)
    index = []
    index.append(np.where(abs(r20_flux - r_flux) == abs(r20_flux - r_flux).min())[0][0])
    index.append(np.where(abs(r80_flux - r_flux) == abs(r80_flux - r_flux).min())[0][0])
    index.append(np.where(abs(r100_flux - r_flux) == abs(r100_flux - r_flux).min())[0][0])
    print(index)
    if fits_plot == True:
        ax=plt.subplot(1,1,1)
        cax=ax.imshow((image),norm=LogNorm(),origin='lower')#,cmap=my_cmap)
        #ax.add_patch(mask.bbox.as_artist(facecolor='none', edgecolor='white'))
        for i in range(grids):
            if i in index:
                ax.add_patch(regions[i].as_artist(facecolor='none', edgecolor='red'))
            else:
                ax.add_patch(regions[i].as_artist(facecolor='none', edgecolor='orange'))
        plt.colorbar(cax)
        plt.show()
    if ifplot == True:
        minorLocator = AutoMinorLocator()
        fig, ax = plt.subplots()
        plt.plot(r_grids, r_flux, 'x-')
        plt.plot(r_grids[index[0]], r_flux[index[0]], 'ro')
        plt.text(r_grids[index[0]]+2, r_flux[index[0]], 'r_20',fontsize=14)
        plt.plot(r_grids[index[1]], r_flux[index[1]], 'ro')
        plt.text(r_grids[index[1]]+2, r_flux[index[1]], 'r_80',fontsize=14)
        plt.plot(r_grids[index[2]], r_flux[index[2]], 'ro')
        plt.text(r_grids[index[2]], r_flux[index[2]]*15/16, 'r model total',fontsize=14)
        ax.xaxis.set_minor_locator(minorLocator)
        plt.tick_params(which='both', width=2)
        plt.tick_params(which='major', length=7)
        plt.tick_params(which='minor', length=4, color='r')
        plt.grid()
        ax.set_ylabel("Total Flux")
        ax.set_xlabel("Pixels")
        if gridspace == 'log':
            ax.set_xscale('log')
            plt.xlim(start_p*0.7, )
        plt.grid(which="minor")
        plt.show()
    return r_grids[index[0]], r_grids[index[1]]

def text_in_string_list(text, string_list):
    '''
    Parameter
    --------
        text: a string that could existed in the string_list
        string_list: The list of name
    Return
    --------
        counts:
            how many items have the text,
            The remainds in the string_list which have text.
    '''
    counts = 0
    text_string=[]
    for i in range(len(string_list)):
        if text in string_list[i]:
            counts += 1
            text_string.append(string_list[i])
    return counts, text_string

def QSO_psfs_compare(QSO, psfs,QSO_msk=None, psf_mask_list=None, PSF_mask_img=None, plt_which_PSF=None,
                     include_QSO = True, gridspace = None , grids=30, norm_pix = 3.0,
                     if_annuli=False,plt_QSO=False, filt = None, not_plt=()):
    """
    Plot the QSO and PSFs SB compare.
    ------
    norm_pix : normalized position. If norm_pix = 'no', means no normalizaion.
    """
    psfs_not_none = [x for x in psfs if x is not None][0]
    if gridspace == None and filt != 'acs':
        radius = 6
    elif gridspace == None and filt == 'acs':
        radius = 6
    elif gridspace == 'log':
        radius = len(psfs_not_none)/6

    if filt != 'acs':
        start_p = 1.5
    elif filt == 'acs':
        start_p = 0.5
#    frm_qrt = int(len(QSO)/2.5)
    if include_QSO == True:
        if plt_QSO ==True:
            print("Plot for QSO:")
        center_QSO = np.array([len(QSO)/2,len(QSO)/2])
        r_SB_QSO, r_grids_QSO = SB_profile(QSO, center=center_QSO, start_p=start_p, radius=radius, grids=grids,
                                           fits_plot=plt_QSO, gridspace=gridspace, if_annuli=if_annuli,
                                           msk_image=QSO_msk)
        if isinstance(norm_pix,int) or isinstance(norm_pix,float):
            count = r_grids_QSO <= norm_pix
            idx = count.sum()-1
#            print "idx:",idx
            r_SB_QSO /= r_SB_QSO[idx]      #normalize the curves
    psfs_NO = len(psfs)
    center = np.array([len(psfs_not_none)/2,len(psfs_not_none)/2])
    if plt_which_PSF != None:
        for i in range(len(plt_which_PSF)):
            j = plt_which_PSF[i]
            if PSF_mask_img is None:
                msk_counts, mask_lists = text_in_string_list("PSF{0}_".format(j), psf_mask_list)
                print("Plot for fits: PSF{0}.fits".format(j))
                if msk_counts == 0:
                    SB_profile(psfs[j], center, start_p=start_p, radius=radius, grids=grids, fits_plot=True, gridspace=gridspace)
                elif msk_counts >0:
                    print(mask_lists)
                    SB_profile(psfs[j], center, start_p=start_p, radius=radius, grids=grids, fits_plot=True, gridspace=gridspace,
                                           mask_plot = False, mask_list=mask_lists)
            else:
                print("Plot for fits: PSF{0}.fits".format(j))
                SB_profile(psfs[j], center, start_p=start_p, radius=radius, grids=grids, fits_plot=True, gridspace=gridspace,
                                           mask_plot = False, msk_image=PSF_mask_img[j])
    minorLocator = AutoMinorLocator()
    fig, ax = plt.subplots(figsize=(10,7))
    for i in range(psfs_NO):
        if i ==0 and include_QSO == True:
            plt.plot(r_grids_QSO, r_SB_QSO, '-', color = 'red', label="QSO", linewidth=5, zorder = 10, alpha=0.5)
            plt.legend()
        if psfs[i] is not None:
            if PSF_mask_img is None:
                msk_counts, mask_lists = text_in_string_list("PSF{0}_".format(i), psf_mask_list)
                if msk_counts == 0:
                    r_SB, r_grids = SB_profile(psfs[i], center, start_p=start_p, radius=radius, grids=grids, gridspace=gridspace, if_annuli=if_annuli)
                elif msk_counts >0:
                    r_SB, r_grids = SB_profile(psfs[i], center, start_p=start_p, radius=radius, grids=grids, gridspace=gridspace, if_annuli=if_annuli,
                                               mask_list=mask_lists)
            else:
                    r_SB, r_grids = SB_profile(psfs[i], center, start_p=start_p, radius=radius, grids=grids, gridspace=gridspace, if_annuli=if_annuli,
                                               msk_image=PSF_mask_img[i])
            if isinstance(norm_pix,int) or isinstance(norm_pix,float):
                count = r_grids <= norm_pix
                idx = count.sum() -1
#                print "idx:",idx
                r_SB /= r_SB[idx]      #normalize the curves
            if gridspace == None:
                plt.text(start_p-0.25,r_SB[0],i)
            elif gridspace == 'log':
                plt.text(start_p-0.1,r_SB[0],i)
#            print r_grids[idx]
            if i not in not_plt:
                plt.plot(r_grids, r_SB, 'x-', label="PSF{0}".format(i))
                plt.legend()
    ax.xaxis.set_minor_locator(minorLocator)
    plt.tick_params(which='both', width=2)
    plt.tick_params(which='major', length=7)
    plt.tick_params(which='minor', length=4, color='r')
    plt.grid()
    if isinstance(norm_pix,int):
        ax.set_ylabel("Scaled Surface Brightness")
    else:
        ax.set_ylabel("Surface Brightness")
    ax.set_xlabel("Pixels")
    if gridspace == 'log':
        ax.set_xscale('log')
        plt.xlim(start_p*0.9, )
    plt.grid(which="minor")
    plt.show()
    return fig

def profiles_compare(prf_list, scal_list, prf_name_list = None, gridspace = None ,
                     grids = 20,  norm_pix = 3, if_annuli=False, filt = None,
                     y_log=False):
    '''
    Compare the profile between different profile (prf?). One can set the scal to uniformize the resolution.
    Note that the SB center from the center of the image; Not allow mask yet.
    Parameter
    --------
        prf_list: profile list of the arraies
        scal_list: a list for the scaled value for the resultion, default as zero, i.e. same resolution.
    Return
    --------
        A plot of SB comparison.
    '''
    if gridspace == None and filt != 'acs':
        radius = 6
    elif gridspace == None and filt == 'acs':
        radius = 6
    elif gridspace == 'log':
        radius = len(prf_list[1])/2

    from matplotlib.ticker import AutoMinorLocator
    minorLocator = AutoMinorLocator()
    fig, ax = plt.subplots(figsize=(10,7))
    prf_NO = len(prf_list)
    for i in range(prf_NO):
        b_c = int(len(prf_list[i])/2)
        b_r = int(len(prf_list[i])/6)
        center = np.reshape(np.asarray(np.where(prf_list[i]== prf_list[i][b_c-b_r:b_c+b_r,b_c-b_r:b_c+b_r].max())),(2))[::-1]
        scale = scal_list[i]
        r_SB, r_grids = SB_profile(prf_list[i], center, radius=radius*scale,
                                   grids=grids, gridspace=gridspace,if_annuli=if_annuli)

        if isinstance(norm_pix,int) or isinstance(norm_pix,float):
            count = r_grids <= norm_pix * scale
            idx = count.sum() -1
#            print "idx:",idx
            r_SB /= r_SB[idx]      #normalize the curves
        r_grids /= scale
        if prf_name_list == None:
            plt.plot(r_grids, r_SB, 'x-', label="prf_list{0}".format(i))
        elif len(prf_name_list)==len(prf_list) and y_log == False:
            plt.plot(r_grids, r_SB, 'x-', label=prf_name_list[i])
        elif len(prf_name_list)==len(prf_list) and y_log == True:
            plt.plot(r_grids, np.log10(r_SB), 'x-', label=prf_name_list[i])
            plt.ylim(0, 0.5)
        else:
            raise ValueError("The profile name is not in right length")
        plt.legend()

    ax.xaxis.set_minor_locator(minorLocator)
    plt.tick_params(which='both', width=2)
    plt.tick_params(which='major', length=7)
    plt.tick_params(which='minor', length=4, color='r')
    plt.grid()
    ax.set_ylabel("Scaled Surface Brightness")
    ax.set_xlabel("Pixels")
    if gridspace == 'log':
        ax.set_xscale('log')
        plt.xlim(1.3, )
    plt.grid(which="minor")
    plt.show()
    return fig

def string_find_between(s, first, last ):
    try:
        start = s.index( first ) + len( first )
        end = s.index( last, start )
        return s[start:end]
    except ValueError:
        return ""

def cr_mask_img(image, mask_list, mask_reg_cut = 0., cut_org=False):
    '''
    Creat a mask image given a mask_list
        if set cut_org=True, the input image don't need to be cutout (e.g. QSO[cut,cut...]).
    '''
    if cut_org ==False:
        mask = np.ones(image.shape)
        for i in range(len(mask_list)):
            mask *= cr_mask(image=image, filename=mask_list[i],mask_reg_cut=mask_reg_cut)
    elif cut_org ==True:
        cut = mask_reg_cut
        mask_org = np.ones(image.shape)
        for i in range(len(mask_list)):
            mask_org *= cr_mask(image=image, filename=mask_list[i])
        mask = np.ones(image[cut:-cut,cut:-cut].shape)
        mask = mask_org[cut:-cut,cut:-cut]
    return mask

def cr_mask(image, filename='test_circle.reg', mask_reg_cut = 0.):
    '''
    The creat a mask with a .reg file. The pixels in the region is 0, otherwise 1.

    Parameter
    --------
        filename: filename of the .reg
        image:

    Return
    --------
        A image.shape array. Pixels in the region is 0, otherwise 1.
    '''
    ##### Note the numpy starts from 0, especially in the center,
    ####!!!Need to check the center, shape of the boxtype, the size?
    with open(filename, 'r') as input_file:
        reg_string=input_file.read().replace('\n', '')
    if "physicalcircle" in reg_string:
        abc=string_find_between(reg_string, "(", ")")
        reg_info=np.fromstring(abc, dtype=float, sep=',')
        center, radius = reg_info[:2]-1 - mask_reg_cut, reg_info[2]
        region = pix_region(center, radius)
        box = 1-region.to_mask(mode='center').data
    elif "physicalbox" in reg_string:
        abc=string_find_between(reg_string, "(", ")")
        reg_info=np.fromstring(abc, dtype=float, sep=',')
        center = reg_info[:2] - 1 - mask_reg_cut
        x_r, y_r = reg_info[2:4]  # x_r is the length of the x, y_r is the length of the y
        box = np.zeros([np.int(x_r)+1, np.int(y_r)+1]).T
    else:
        print(reg_string)
        raise ValueError("The input reg is un-defined yet")
    frame_size = image.shape
    box_size = box.shape
    x_edge = np.int(center[1]-box_size[0]/2) #The position of the center is x-y switched.
    y_edge = np.int(center[0]-box_size[1]/2)
    mask = np.ones(frame_size)
    mask_box_part = mask[x_edge:x_edge+box_size[0],y_edge: y_edge + box_size[1]]
    mask_box_part *= box
    return mask

def total_compare(label_list, flux_list, zp=27.0, target_ID = 'target_ID',
                  add_background=0.0, data_mask_list = None, data_cut = 0.,plot_compare=False,
                  pix_sz = 'drz06', msk_image=None, if_annuli=False, arrows=False, host_comp_name=None, host_comp=None):
    if pix_sz == 'swarp':
        delatPixel = 0.127985
    elif pix_sz == 'drz06':
        delatPixel = 0.0642
    elif pix_sz == 'acs':
        delatPixel = 0.03
    elif isinstance(pix_sz,float):
        delatPixel = pix_sz
    else:
        delatPixel = 1.
        print("Warning: The pix_sz is inappropriate and thus delatPixel=1.0")

    norm = LogNorm() #ImageNormalize(stretch=SqrtStretch())
    f = plt.figure(0, figsize=(20.1,4))
    ax1 = plt.subplot2grid((6,5), (0,0), rowspan=6)
    ax2 = plt.subplot2grid((6,5), (0,1), rowspan=6)
    axE = plt.subplot2grid((6,5), (0,2), rowspan=6)
    ax3 = plt.subplot2grid((6,5), (0,3), rowspan=6)
    ax4 = plt.subplot2grid((6,5), (0,4), rowspan=5)
    ax5 = plt.subplot2grid((6,5), (5,4), rowspan=1)

    im1 = ax1.imshow(flux_list[0] + add_background,origin='lower',cmap=my_cmap, norm=norm, vmax = flux_list[0].max())
    clim=im1.properties()['clim']
    frame_size = len(flux_list[0])
    ax1.set_ylabel(target_ID, fontsize=15, weight='bold')
    # ax1.text(frame_size*0.05, frame_size*0.9, label_list[0],fontsize=20,weight='bold', color='white')
    ax1.get_xaxis().set_visible(False)
#    ax1.get_yaxis().set_visible(False)
    scale_bar(ax1, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(ax1, frame_size, arrow_size=0.03, color = 'white')
    cb1 = f.colorbar(im1, ax=ax1, shrink=0.48, pad=0.01,  orientation="horizontal", aspect=15, ticks= [1.e-4, 1.e-3, 1.e-2,1.e-1,0, 10])
    cb1.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1,10,100])
#    cb1.ax.()
#    print("lenfluxlist")
#    print(len(flux_list))
#    print("Flux List 1:")
#    print(len(flux_list[1]))
#    print("Flux List 1:")
#    print(flux_list[1][0][0])
#    print("BG")
#    print(add_background)
#    print("flux")
#    print(len(flux_list))
#    print("soom")
    soom=flux_list[1] + flux_list[2] + add_background

    djazz=soom[0][:][:]
#    print((soom[0][:][:]).shape)
#    print("S1")
#    print(soom[0][1][:])
#    print("S2")
#    print(soom[0][2][:])
    im2 = ax2.imshow(flux_list[1] + flux_list[2] + add_background,origin='lower',cmap="gist_heat", norm=norm, clim=clim)
#    pos2_o = ax2.get_position() # get the original position
#    pos2 = [pos2_o.x0 -0.03, pos2_o.y0, pos2_o.width, pos2_o.height]
#    ax2.set_position(pos2) # set a new position
    # ax2.text(frame_size*0.05, frame_size*0.9, label_list[-2],weight='bold',
         # fontsize=20, color='white')
    scale_bar(ax2, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(ax2, frame_size, arrow_size=0.03, color = 'white')
    ax2.get_xaxis().set_visible(False)
    ax2.get_yaxis().set_visible(False)
    cb2 = f.colorbar(im2, ax=ax2, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
    cb2.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1,10,100])

#    posE_o = axE.get_position() # get the original position
#    posE = [posE_o.x0 -0.1, pos3_o.y0 +0.025, pos3_o.width, pos3_o.height]
#    axE.set_position(posE) # set a new position
    diffjazz=flux_list[0]-flux_list[1]
    djazz1=diffjazz[0][:][:]
    #HERE GIRL
    imE = axE.imshow((flux_list[0]-flux_list[1])*msk_image, origin='lower',cmap=my_cmap, norm=norm, clim=clim)
    # axE.text(frame_size*0.05, frame_size*0.9, 'data - Point Source', weight='bold',
         # fontsize=17, color='white')
    scale_bar(axE, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(axE, frame_size, arrow_size=0.03, color = 'white')
    axE.get_xaxis().set_visible(False)
    axE.get_yaxis().set_visible(False)
    cbE =f.colorbar(imE, ax=axE, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
    cbE.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1,10,100])

    normjazz=(flux_list[0]-(flux_list[1]+flux_list[2]))/flux_list[3]
    djazz2=normjazz[0][:][:]
    norm_residual = (flux_list[0]-(flux_list[1]+flux_list[2]))/flux_list[3] * msk_image
#    pos3_o = ax3.get_position() # get the original position
#    pos3 = [pos3_o.x0 -0.1, pos3_o.y0 +0.025, pos3_o.width, pos3_o.height]
#    ax3.set_position(pos3) # set a new position
    im3 = ax3.imshow(norm_residual, origin='lower',cmap='bwr', vmin=-6, vmax=6)
    # ax3.text(frame_size*0.02, frame_size*0.9, label_list[-1],weight='bold',
    #      fontsize=17, color='black')
    scale_bar(ax3, frame_size, dist=1/delatPixel, text='1"')
    if arrows == True:
        coordinate_arrows(ax3, frame_size, arrow_size=0.03)
    ax3.get_xaxis().set_visible(False)
    ax3.get_yaxis().set_visible(False)
    f.colorbar(im3, ax=ax3, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
#    plt.tight_layout()
    plt.subplots_adjust(wspace=-0.5, hspace=0)
#    make_ticklabels_invisible(plt.gcf())
    for i in range(len(flux_list)-1):
        if i == 0:
            model_flux = flux_list[i+1] +0  # Don't share a same space
        else:
            model_flux += flux_list[2]
    model_flux = flux_list[1] + flux_list[2]

    label_SB_list = [label_list[0], label_list[-2], label_list[1], label_list[2]]
    flux_SB_list = [flux_list[0], model_flux, flux_list[1], flux_list[2]]
    radi = len(flux_list[0])/2
    if if_annuli == False:
        for i in range(len(label_SB_list)):
            center = len(flux_SB_list[i])/2, len(flux_SB_list[i])/2
            if label_SB_list[i] == 'data':
#                print "data_mask_lists:\t", data_mask_list
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log',
                                           radius= radi, grids = 50, mask_list=data_mask_list,
                                           mask_cut = data_cut, msk_image=msk_image, fits_plot=False)
            else:
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log', radius= radi,grids = 30, mask_list=None)
            r_mag = - 2.5 * np.log10(r_SB) + zp
            if label_SB_list[i] == 'data':
                ind = len(r_mag)-(r_mag == r_mag[-1]).sum()
                ax4.plot(r_grids[:ind], r_mag[:ind], 'o', color = 'whitesmoke',markeredgecolor="black", label=label_SB_list[i])
#                r_max =  r_grids[ind]
            else:
                ax4.plot(r_grids, r_mag, '-', label=label_SB_list[i])
        if host_comp_name is not None and host_comp is not None:
            for i in range(len(host_comp_name)):
                center = len(host_comp[i])/2, len(host_comp[i])/2
                r_SB, r_grids = SB_profile(host_comp[i], center, gridspace = 'log', radius= radi,grids = 30, mask_list=None)
                r_mag = - 2.5 * np.log10(r_SB) + zp
                ax4.plot(r_grids, r_mag, '-', label=host_comp_name[i])
        ax4.set_xlabel('pixel', fontsize=15)
        ax4.xaxis.set_label_position('top')
        ax4.xaxis.tick_top()
        ax4.set_xscale('log')
        ax4.set_xticks([2,4,6,10,15,20,30,50,100,150])
        from matplotlib.ticker import ScalarFormatter
        ax4.xaxis.set_major_formatter(ScalarFormatter())
        ax4.set_xlim([(r_grids).min()*0.85,r_grids.max()+6])

        ax4.invert_yaxis()
        ax4.set_ylabel('$\mu$(mag, pixel$^{-2}$)', fontsize=12)
        ax4.yaxis.set_label_position('right')
        ax4.yaxis.tick_right()
        ax4.yaxis.set_ticks_position('both')
        plt.gca().invert_yaxis()
        ax4.legend()
        x = np.linspace(1.e-4, 100, 2)
        y = x * 0
        r_mag_0 = 2.5 * np.log10(SB_profile(flux_SB_list[0], center, gridspace = 'log', radius= radi,
                                            mask_list=data_mask_list, mask_cut = data_cut,grids = 30,
                                            msk_image=msk_image)[0])
        r_mag_1 = 2.5 * np.log10(SB_profile(flux_SB_list[1], center, gridspace = 'log', grids = 30,radius= radi)[0])
        ind = len(r_mag_0)-(r_mag_0 == r_mag_0[-1]).sum()
        ax5.plot(r_grids[:ind]*delatPixel, (r_mag_0-r_mag_1)[:ind], 'ro')
        ax5.set_ylabel('$\Delta\mu$', fontsize=15)
        ax5.set_xlabel('arcsec', fontsize=15)
        ax5.set_xscale('log')
        ax5.set_xticks([0.1, 0.2, 0.5, 1, 2,5,10,20])
        ax5.set_yticks([-0.5,-0.25, 0., 0.25])
        import matplotlib
        ax5.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax5.plot(x, y, 'k--')
        ax5.yaxis.set_label_position('right')
        ax5.yaxis.tick_right()
        ax5.yaxis.set_ticks_position('both')
        ax5.set_xlim([(r_grids*delatPixel).min()*0.85, (r_grids.max()+6)*delatPixel])
        plt.ylim([-0.5,0.5])
    elif if_annuli == True:
        for i in range(len(label_SB_list)):
            center = len(flux_SB_list[i])/2, len(flux_SB_list[i])/2
            if label_SB_list[i] == 'data':
#                print "data_mask_lists:\t", data_mask_list
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log',
                                           radius= radi, grids = 50, mask_list=data_mask_list,
                                           mask_cut = data_cut, msk_image=msk_image, fits_plot=False, if_annuli = if_annuli)
                ax4.plot(r_grids, r_SB, 'o', color = 'whitesmoke',markeredgecolor="black", label=label_SB_list[i])
            else:
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log', radius=radi,grids = 30, mask_list=None, if_annuli = if_annuli)
                ax4.plot(r_grids, r_SB, '-', label=label_SB_list[i])
        if host_comp_name is not None and host_comp is not None:
            for i in range(len(host_comp_name)):
                center = len(host_comp[i])/2, len(host_comp[i])/2
                r_SB, r_grids = SB_profile(host_comp[i], center, gridspace = 'log', radius= radi,grids = 30, mask_list=None)
                ax4.plot(r_grids, r_SB, '-', label=host_comp_name[i])
        ax4.set_xlabel('pixel', fontsize=15)
        ax4.xaxis.set_label_position('top')
        ax4.xaxis.tick_top()
        ax4.set_xscale('log')
        ax4.set_xticks([2,4,6,10,15,20,30,50,100,150])
        from matplotlib.ticker import ScalarFormatter
        ax4.xaxis.set_major_formatter(ScalarFormatter())
        ax4.set_xlim([(r_grids).min()*0.85,(r_grids).max()+6])

        ax4.set_ylabel('$SB_{annuli}$(counts, pixel$^{-2}$)', fontsize=12)
        ax4.yaxis.set_label_position('right')
        ax4.yaxis.tick_right()
        ax4.yaxis.set_ticks_position('both')
        ax4.legend()

        x = np.linspace(1.e-4, 100, 2)
        y = x * 0
        r_SB_0 = (SB_profile(flux_SB_list[0], center, gridspace = 'log', radius= radi, if_annuli = if_annuli,
                                            mask_list=data_mask_list, mask_cut = data_cut,grids = 30,
                                            msk_image=msk_image)[0])
        r_SB_1 = (SB_profile(flux_SB_list[1], center, gridspace = 'log', grids = 30, if_annuli = if_annuli,radius= radi)[0])
    #    r_mag_diff = 2.5 * np.log10(SB_profile(flux_SB_list[1], center, gridspace = 'log', radius= radi)[0])
        ax5.plot(r_grids*delatPixel, (r_SB_0- r_SB_1), 'ro')
        ax5.set_ylabel('$\Delta SB$', fontsize=15)
        ax5.set_xlabel('arcsec', fontsize=15)
        ax5.set_xscale('log')
        ax5.set_xticks([0.1, 0.2, 0.5, 1, 2,5,10,20])
        import matplotlib
        ax5.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax5.plot(x, y, 'k--')
        ax5.yaxis.set_label_position('right')
        ax5.yaxis.tick_right()
        ax5.yaxis.set_ticks_position('both')
        ax5.set_xlim([(r_grids*delatPixel).min()*0.85, (r_grids.max()+6)*delatPixel])
        plt.ylim([-5,5])
        ax5.set_yticks([-5,-2.5, 0., 2.5])
    pos4_o = ax4.get_position() # get the original position
    pos4 = [pos4_o.x0+0.067, pos4_o.y0 + 0.10, pos4_o.width*0.72, pos4_o.height*0.8]
    ax4.set_position(pos4) # set a new position
    pos5_o = ax5.get_position() # get the original position
    pos5 = [pos5_o.x0+0.067, pos5_o.y0+0.08, pos5_o.width*0.72, pos5_o.height*1.1]
    ax5.set_position(pos5) # set a new position
    ax1.text(frame_size*0.05, frame_size*1.05, label_list[0],fontsize=20,weight='bold', color='black')
    ax2.text(frame_size*0.05, frame_size*1.05, label_list[-2],weight='bold',
      fontsize=20, color='black')
    ax3.text(frame_size*0.05, frame_size*1.05, label_list[-1],weight='bold',
          fontsize=17, color='black')
    axE.text(frame_size*0.05, frame_size*1.05, 'data - Point Source', weight='bold',
         fontsize=17, color='black')
    if plot_compare == True:
        plt.show()
    else:
        plt.close()
    return f

def total_compare_6(label_list, flux_list, zp=27.0, target_ID = 'target_ID',
                  add_background=0.0, data_mask_list = None, data_cut = 0.,plot_compare=False,
                  pix_sz = 'drz06', msk_image=None, if_annuli=False, arrows=False, host_comp_name=None, host_comp=None):
    if pix_sz == 'swarp':
        delatPixel = 0.127985
    elif pix_sz == 'drz06':
        delatPixel = 0.0642
    elif pix_sz == 'acs':
        delatPixel = 0.03
    elif isinstance(pix_sz,float):
        delatPixel = pix_sz
    else:
        delatPixel = 1.
        print("Warning: The pix_sz is inappropriate and thus delatPixel=1.0")

    norm = LogNorm() #ImageNormalize(stretch=SqrtStretch())
    f = plt.figure(0, figsize=(24.5,4))
    ax1 = plt.subplot2grid((6,6), (0,0), rowspan=6)
    ax2 = plt.subplot2grid((6,6), (0,1), rowspan=6)
    axE = plt.subplot2grid((6,6), (0,2), rowspan=6)
    axR = plt.subplot2grid((6,6), (0,3), rowspan=6)
    ax3 = plt.subplot2grid((6,6), (0,4), rowspan=6)
    ax4 = plt.subplot2grid((6,6), (0,5), rowspan=5)
    ax5 = plt.subplot2grid((6,6), (5,5), rowspan=1)

    im1 = ax1.imshow(flux_list[0] + add_background,origin='lower',cmap=my_cmap, norm=norm, vmax = flux_list[0].max())
    clim=im1.properties()['clim']
    frame_size = len(flux_list[0])
    ax1.set_ylabel(target_ID, fontsize=15, weight='bold')
    # ax1.text(frame_size*0.05, frame_size*0.9, label_list[0],fontsize=20,weight='bold', color='white')
    ax1.get_xaxis().set_visible(False)
#    ax1.get_yaxis().set_visible(False)
    scale_bar(ax1, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(ax1, frame_size, arrow_size=0.03, color = 'white')
    cb1 = f.colorbar(im1, ax=ax1, shrink=0.48, pad=0.01,  orientation="horizontal", aspect=15, ticks= [1.e-4, 1.e-3, 1.e-2,1.e-1,0, 10])
    cb1.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1,10,100])
#    cb1.ax.()

    im2 = ax2.imshow(flux_list[1] + flux_list[2] + add_background,origin='lower',cmap="gist_heat", norm=norm, clim=clim)
#    pos2_o = ax2.get_position() # get the original position
#    pos2 = [pos2_o.x0 -0.03, pos2_o.y0, pos2_o.width, pos2_o.height]
#    ax2.set_position(pos2) # set a new position
    # ax2.text(frame_size*0.05, frame_size*0.9, label_list[-2],weight='bold',
         # fontsize=20, color='white')
    scale_bar(ax2, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(ax2, frame_size, arrow_size=0.03, color = 'white')
    ax2.get_xaxis().set_visible(False)
    ax2.get_yaxis().set_visible(False)
    cb2 = f.colorbar(im2, ax=ax2, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
    cb2.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1,10,100])

#    posE_o = axE.get_position() # get the original position
#    posE = [posE_o.x0 -0.1, pos3_o.y0 +0.025, pos3_o.width, pos3_o.height]
#    axE.set_position(posE) # set a new position
    imE = axE.imshow((flux_list[0] - flux_list[1])*msk_image, origin='lower',cmap=my_cmap, norm=norm, clim=clim)
    # axE.text(frame_size*0.05, frame_size*0.9, 'data - Point Source', weight='bold',
         # fontsize=17, color='white')
    scale_bar(axE, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(axE, frame_size, arrow_size=0.03, color = 'white')
    axE.get_xaxis().set_visible(False)
    axE.get_yaxis().set_visible(False)
    cbE =f.colorbar(imE, ax=axE, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
    cbE.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1,10,100])

    imR = axR.imshow((flux_list[0] - flux_list[1]- flux_list[2])*msk_image, origin='lower',cmap=my_cmap, norm=norm, clim=clim)
    scale_bar(axR, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(axR, frame_size, arrow_size=0.03, color = 'white')
    axR.get_xaxis().set_visible(False)
    axR.get_yaxis().set_visible(False)
    cbR =f.colorbar(imR, ax=axR, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
    cbR.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1,10,100])

    norm_residual = (flux_list[0]-(flux_list[1]+flux_list[2]))/flux_list[3] * msk_image
#    pos3_o = ax3.get_position() # get the original position
#    pos3 = [pos3_o.x0 -0.1, pos3_o.y0 +0.025, pos3_o.width, pos3_o.height]
#    ax3.set_position(pos3) # set a new position
    im3 = ax3.imshow(norm_residual, origin='lower',cmap='bwr', vmin=-6, vmax=6)
    # ax3.text(frame_size*0.02, frame_size*0.9, label_list[-1],weight='bold',
    #      fontsize=17, color='black')
    scale_bar(ax3, frame_size, dist=1/delatPixel, text='1"')
    if arrows == True:
        coordinate_arrows(ax3, frame_size, arrow_size=0.03)
    ax3.get_xaxis().set_visible(False)
    ax3.get_yaxis().set_visible(False)
    f.colorbar(im3, ax=ax3, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
#    plt.tight_layout()
    plt.subplots_adjust(wspace=-0.5, hspace=0)
#    make_ticklabels_invisible(plt.gcf())
    for i in range(len(flux_list)-1):
        if i == 0:
            model_flux = flux_list[i+1] +0  # Don't share a same space
        else:
            model_flux += flux_list[2]
    model_flux = flux_list[1] + flux_list[2]

    label_SB_list = [label_list[0], label_list[-2], label_list[1], label_list[2]]
    flux_SB_list = [flux_list[0], model_flux, flux_list[1], flux_list[2]]
    radi = len(flux_list[0])/2
    if if_annuli == False:
        for i in range(len(label_SB_list)):
            center = len(flux_SB_list[i])/2, len(flux_SB_list[i])/2
            if label_SB_list[i] == 'data':
#                print "data_mask_lists:\t", data_mask_list
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log',
                                           radius= radi, grids = 50, mask_list=data_mask_list,
                                           mask_cut = data_cut, msk_image=msk_image, fits_plot=False)
            else:
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log', radius= radi,grids = 30, mask_list=None)
            r_mag = - 2.5 * np.log10(r_SB) + zp
            if label_SB_list[i] == 'data':
                ind = len(r_mag)-(r_mag == r_mag[-1]).sum()
                ax4.plot(r_grids[:ind], r_mag[:ind], 'o', color = 'whitesmoke',markeredgecolor="black", label=label_SB_list[i])
#                r_max =  r_grids[ind]
            else:
                ax4.plot(r_grids, r_mag, '-', label=label_SB_list[i])
        if host_comp_name is not None and host_comp is not None:
            for i in range(len(host_comp_name)):
                center = len(host_comp[i])/2, len(host_comp[i])/2
                r_SB, r_grids = SB_profile(host_comp[i], center, gridspace = 'log', radius= radi,grids = 30, mask_list=None)
                r_mag = - 2.5 * np.log10(r_SB) + zp
                ax4.plot(r_grids, r_mag, '-', label=host_comp_name[i])
        ax4.set_xlabel('pixel', fontsize=15)
        ax4.xaxis.set_label_position('top')
        ax4.xaxis.tick_top()
        ax4.set_xscale('log')
        ax4.set_xticks([2,4,6,10,15,20,30,50,100,150])
        from matplotlib.ticker import ScalarFormatter
        ax4.xaxis.set_major_formatter(ScalarFormatter())
        ax4.set_xlim([(r_grids).min()*0.85,r_grids.max()+6])

        ax4.invert_yaxis()
        ax4.set_ylabel('$\mu$(mag, pixel$^{-2}$)', fontsize=12)
        ax4.yaxis.set_label_position('right')
        ax4.yaxis.tick_right()
        ax4.yaxis.set_ticks_position('both')
        plt.gca().invert_yaxis()
        ax4.legend()
        x = np.linspace(1.e-4, 100, 2)
        y = x * 0
        r_mag_0 = 2.5 * np.log10(SB_profile(flux_SB_list[0], center, gridspace = 'log', radius= radi,
                                            mask_list=data_mask_list, mask_cut = data_cut,grids = 30,
                                            msk_image=msk_image)[0])
        r_mag_1 = 2.5 * np.log10(SB_profile(flux_SB_list[1], center, gridspace = 'log', grids = 30,radius= radi)[0])
        ind = len(r_mag_0)-(r_mag_0 == r_mag_0[-1]).sum()
        ax5.plot(r_grids[:ind]*delatPixel, (r_mag_0-r_mag_1)[:ind], 'ro')
        ax5.set_ylabel('$\Delta\mu$', fontsize=15)
        ax5.set_xlabel('arcsec', fontsize=15)
        ax5.set_xscale('log')
        ax5.set_xticks([0.1, 0.2, 0.5, 1, 2,5,10,20])
        ax5.set_yticks([-0.5,-0.25, 0., 0.25])
        import matplotlib
        ax5.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax5.plot(x, y, 'k--')
        ax5.yaxis.set_label_position('right')
        ax5.yaxis.tick_right()
        ax5.yaxis.set_ticks_position('both')
        ax5.set_xlim([(r_grids*delatPixel).min()*0.85, (r_grids.max()+6)*delatPixel])
        plt.ylim([-0.5,0.5])
    elif if_annuli == True:
        for i in range(len(label_SB_list)):
            center = len(flux_SB_list[i])/2, len(flux_SB_list[i])/2
            if label_SB_list[i] == 'data':
#                print "data_mask_lists:\t", data_mask_list
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log',
                                           radius= radi, grids = 50, mask_list=data_mask_list,
                                           mask_cut = data_cut, msk_image=msk_image, fits_plot=False, if_annuli = if_annuli)
                ax4.plot(r_grids, r_SB, 'o', color = 'whitesmoke',markeredgecolor="black", label=label_SB_list[i])
            else:
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log', radius=radi,grids = 30, mask_list=None, if_annuli = if_annuli)
                ax4.plot(r_grids, r_SB, '-', label=label_SB_list[i])
        if host_comp_name is not None and host_comp is not None:
            for i in range(len(host_comp_name)):
                center = len(host_comp[i])/2, len(host_comp[i])/2
                r_SB, r_grids = SB_profile(host_comp[i], center, gridspace = 'log', radius= radi,grids = 30, mask_list=None)
                ax4.plot(r_grids, r_SB, '-', label=host_comp_name[i])
        ax4.set_xlabel('pixel', fontsize=15)
        ax4.xaxis.set_label_position('top')
        ax4.xaxis.tick_top()
        ax4.set_xscale('log')
        ax4.set_xticks([2,4,6,10,15,20,30,50,100,150])
        from matplotlib.ticker import ScalarFormatter
        ax4.xaxis.set_major_formatter(ScalarFormatter())
        ax4.set_xlim([(r_grids).min()*0.85,(r_grids).max()+6])

        ax4.set_ylabel('$SB_{annuli}$(counts, pixel$^{-2}$)', fontsize=12)
        ax4.yaxis.set_label_position('right')
        ax4.yaxis.tick_right()
        ax4.yaxis.set_ticks_position('both')
        ax4.legend()

        x = np.linspace(1.e-4, 100, 2)
        y = x * 0
        r_SB_0 = (SB_profile(flux_SB_list[0], center, gridspace = 'log', radius= radi, if_annuli = if_annuli,
                                            mask_list=data_mask_list, mask_cut = data_cut,grids = 30,
                                            msk_image=msk_image)[0])
        r_SB_1 = (SB_profile(flux_SB_list[1], center, gridspace = 'log', grids = 30, if_annuli = if_annuli,radius= radi)[0])
    #    r_mag_diff = 2.5 * np.log10(SB_profile(flux_SB_list[1], center, gridspace = 'log', radius= radi)[0])
        ax5.plot(r_grids*delatPixel, (r_SB_0- r_SB_1), 'ro')
        ax5.set_ylabel('$\Delta SB$', fontsize=15)
        ax5.set_xlabel('arcsec', fontsize=15)
        ax5.set_xscale('log')
        ax5.set_xticks([0.1, 0.2, 0.5, 1, 2,5,10,20])
        import matplotlib
        ax5.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax5.plot(x, y, 'k--')
        ax5.yaxis.set_label_position('right')
        ax5.yaxis.tick_right()
        ax5.yaxis.set_ticks_position('both')
        ax5.set_xlim([(r_grids*delatPixel).min()*0.85, (r_grids.max()+6)*delatPixel])
        plt.ylim([-5,5])
        ax5.set_yticks([-5,-2.5, 0., 2.5])
    pos4_o = ax4.get_position() # get the original position
    pos4 = [pos4_o.x0+0.067, pos4_o.y0 + 0.10, pos4_o.width*0.72, pos4_o.height*0.8]
    ax4.set_position(pos4) # set a new position
    pos5_o = ax5.get_position() # get the original position
    pos5 = [pos5_o.x0+0.067, pos5_o.y0+0.08, pos5_o.width*0.72, pos5_o.height*1.1]
    ax5.set_position(pos5) # set a new position
    ax1.text(frame_size*0.05, frame_size*1.05, label_list[0],fontsize=20,weight='bold', color='black')
    ax2.text(frame_size*0.05, frame_size*1.05, label_list[-2],weight='bold',
      fontsize=20, color='black')
    ax3.text(frame_size*0.05, frame_size*1.05, label_list[-1],weight='bold',
          fontsize=17, color='black')
    axE.text(frame_size*0.05, frame_size*1.05, 'data - Point Sources', weight='bold',
         fontsize=17, color='black')
    axR.text(frame_size*0.05, frame_size*1.05, 'data - model', weight='bold',
         fontsize=17, color='black')
    if plot_compare == True:
        plt.show()
    else:
        plt.close()
    return f

def galaxy_total_compare(label_list, flux_list, zp=27.0, pix_sz=1., target_ID = 'target_ID',
                  add_background=0.0, data_mask_list = None, data_cut = 0.,plot_compare=False,
                  msk_image=None, if_annuli=False, arrows=False):
#    if facility == 'F140w':
#        zp = 26.4524
#    elif facility == 'F125w':
#        zp = 26.2303
#    elif facility == 'acs':
#        zp = 25.94333
    if pix_sz == 'swarp':
        delatPixel = 0.127985
    elif pix_sz == 'drz06':
        delatPixel = 0.0642
    elif pix_sz == 'acs':
        delatPixel = 0.03
    elif isinstance(pix_sz,float):
        delatPixel = pix_sz
    else:
        delatPixel = 1.
        print("Warning: The pix_sz is inappropriate and thus delatPixel=1.0")

    norm = LogNorm() #ImageNormalize(stretch=SqrtStretch())
    f = plt.figure(0, figsize=(17,4))
    ax1 = plt.subplot2grid((6,4), (0,0), rowspan=6)
    ax2 = plt.subplot2grid((6,4), (0,1), rowspan=6)
#    axE = plt.subplot2grid((6,4), (0,2), rowspan=6)
    ax3 = plt.subplot2grid((6,4), (0,2), rowspan=6)
    ax4 = plt.subplot2grid((6,4), (0,3), rowspan=5)
    ax5 = plt.subplot2grid((6,4), (5,3), rowspan=1)

    im1 = ax1.imshow(flux_list[0] + add_background,origin='lower',cmap="gist_heat", norm=norm, vmax = flux_list[0].max())
    clim=im1.properties()['clim']
    frame_size = len(flux_list[0])
    ax1.set_ylabel(target_ID, fontsize=15, weight='bold')
    ax1.text(frame_size*0.05, frame_size*0.9, label_list[0],fontsize=20,weight='bold', color='white')
    ax1.get_xaxis().set_visible(False)
#    ax1.get_yaxis().set_visible(False)
    scale_bar(ax1, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(ax1, frame_size, arrow_size=0.03, color = 'white')
    cb1 = f.colorbar(im1, ax=ax1, shrink=0.48, pad=0.01,  orientation="horizontal", aspect=15, ticks= [1.e-4, 1.e-3, 1.e-2,1.e-1,0, 10])
    cb1.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1])
#    cb1.ax.()

    im2 = ax2.imshow(flux_list[1] + flux_list[2] + add_background,origin='lower',cmap="gist_heat", norm=norm, clim=clim)
    ax2.text(frame_size*0.05, frame_size*0.9, label_list[-2],weight='bold',
         fontsize=20, color='white')
    scale_bar(ax2, frame_size, dist=1/delatPixel, text='1"', color = 'white')
    if arrows == True:
        coordinate_arrows(ax2, frame_size, arrow_size=0.03, color = 'white')
    ax2.get_xaxis().set_visible(False)
    ax2.get_yaxis().set_visible(False)
    cb2 = f.colorbar(im2, ax=ax2, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
    cb2.set_ticks([1.e-5, 1.e-4, 1.e-3, 1.e-2,1.e-1,0,1])

    norm_residual = (flux_list[0]-(flux_list[1]+flux_list[2]))/flux_list[3] * msk_image
#    pos3_o = ax3.get_position() # get the original position
#    pos3 = [pos3_o.x0 -0.1, pos3_o.y0 +0.025, pos3_o.width, pos3_o.height]
#    ax3.set_position(pos3) # set a new position
    im3 = ax3.imshow(norm_residual, origin='lower',cmap='bwr', vmin=-6, vmax=6)
    ax3.text(frame_size*0.02, frame_size*0.9, label_list[-1],weight='bold',
         fontsize=17, color='black')
    scale_bar(ax3, frame_size, dist=1/delatPixel, text='1"')
    if arrows == True:
        coordinate_arrows(ax3, frame_size, arrow_size=0.03)
    ax3.get_xaxis().set_visible(False)
    ax3.get_yaxis().set_visible(False)
    f.colorbar(im3, ax=ax3, shrink=0.48, pad=0.01,   orientation="horizontal", aspect=15)
#    plt.tight_layout()
    plt.subplots_adjust(wspace=-0.5, hspace=0)
#    make_ticklabels_invisible(plt.gcf())
    for i in range(len(flux_list)-1):
        if i == 0:
            model_flux = flux_list[i+1] +0  # Don't share a same space
        else:
            model_flux += flux_list[2]
    model_flux = flux_list[1] + flux_list[2]

    label_SB_list = [label_list[0],label_list[1]]
    flux_SB_list = [flux_list[0], model_flux]
    radi = len(flux_list[0])/2
    if if_annuli == False:
        for i in range(len(label_SB_list)):
            center = len(flux_SB_list[i])/2, len(flux_SB_list[i])/2
            if label_SB_list[i] == 'data':
#                print "data_mask_lists:\t", data_mask_list
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log',
                                           radius= radi, grids = 50, mask_list=data_mask_list,
                                           mask_cut = data_cut, msk_image=msk_image, fits_plot=False)
            else:
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log', radius= radi,grids = 30, mask_list=None)
            r_mag = - 2.5 * np.log10(r_SB) + zp
            if label_SB_list[i] == 'data':
                ind = len(r_mag)-(r_mag == r_mag[-1]).sum()
                ax4.plot(r_grids[:ind], r_mag[:ind], 'o', color = 'whitesmoke',markeredgecolor="black", label=label_SB_list[i])
#                r_max =  r_grids[ind]
            else:
                ax4.plot(r_grids, r_mag, '-', label=label_SB_list[i])
        ax4.set_xlabel('pixel', fontsize=15)
        ax4.xaxis.set_label_position('top')
        ax4.xaxis.tick_top()
        ax4.set_xscale('log')
        ax4.set_xticks([2,4,6,10,15,20,30])
        from matplotlib.ticker import ScalarFormatter
        ax4.xaxis.set_major_formatter(ScalarFormatter())
        ax4.set_xlim([(r_grids).min()*0.85,r_grids.max()+6])

        ax4.invert_yaxis()
        ax4.set_ylabel('$\mu$(mag, pixel$^{-2}$)', fontsize=12)
        ax4.yaxis.set_label_position('right')
        ax4.yaxis.tick_right()
        ax4.yaxis.set_ticks_position('both')
        plt.gca().invert_yaxis()

        ax4.legend()
        x = np.linspace(1.e-4, 100, 2)
        y = x * 0
        r_mag_0 = 2.5 * np.log10(SB_profile(flux_SB_list[0], center, gridspace = 'log', radius= radi,
                                            mask_list=data_mask_list, mask_cut = data_cut,grids = 30,
                                            msk_image=msk_image)[0])
        r_mag_1 = 2.5 * np.log10(SB_profile(flux_SB_list[1], center, gridspace = 'log', grids = 30,radius= radi)[0])
        ind = len(r_mag_0)-(r_mag_0 == r_mag_0[-1]).sum()
        ax5.plot(r_grids[:ind]*delatPixel, (r_mag_0-r_mag_1)[:ind], 'ro')
        ax5.set_ylabel('$\Delta\mu$', fontsize=15)
        ax5.set_xlabel('arcsec', fontsize=15)
        ax5.set_xscale('log')
        ax5.set_xticks([0.1, 0.2, 0.5, 1, 2])
        ax5.set_yticks([-0.5,-0.25, 0., 0.25])
        import matplotlib
        ax5.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax5.plot(x, y, 'k--')
        ax5.yaxis.set_label_position('right')
        ax5.yaxis.tick_right()
        ax5.yaxis.set_ticks_position('both')
        ax5.set_xlim([(r_grids*delatPixel).min()*0.85, (r_grids.max()+6)*delatPixel])
        plt.ylim([-0.5,0.5])
    elif if_annuli == True:
        for i in range(len(label_SB_list)):
            center = len(flux_SB_list[i])/2, len(flux_SB_list[i])/2
            if label_SB_list[i] == 'data':
                print("data_mask_lists:\t", data_mask_list)
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log',
                                           radius= radi, grids = 50, mask_list=data_mask_list,
                                           mask_cut = data_cut, msk_image=msk_image, fits_plot=False, if_annuli = if_annuli)
                ax4.plot(r_grids, r_SB, 'o', color = 'whitesmoke',markeredgecolor="black", label=label_SB_list[i])
            else:
                r_SB, r_grids = SB_profile(flux_SB_list[i], center, gridspace = 'log', radius=radi,grids = 30, mask_list=None, if_annuli = if_annuli)
                ax4.plot(r_grids, r_SB, '-', label=label_SB_list[i])
        ax4.set_xlabel('pixel', fontsize=15)
        ax4.xaxis.set_label_position('top')
        ax4.xaxis.tick_top()
        ax4.set_xscale('log')
        ax4.set_xticks([2,4,6,10,15,20,30])
        from matplotlib.ticker import ScalarFormatter
        ax4.xaxis.set_major_formatter(ScalarFormatter())
        ax4.set_xlim([(r_grids).min()*0.85,(r_grids).max()+6])

        ax4.set_ylabel('$SB_{annuli}$(counts, pixel$^{-2}$)', fontsize=12)
        ax4.yaxis.set_label_position('right')
        ax4.yaxis.tick_right()
        ax4.yaxis.set_ticks_position('both')
        ax4.legend()

        x = np.linspace(1.e-4, 100, 2)
        y = x * 0
        r_SB_0 = (SB_profile(flux_SB_list[0], center, gridspace = 'log', radius= radi, if_annuli = if_annuli,
                                            mask_list=data_mask_list, mask_cut = data_cut,grids = 30,
                                            msk_image=msk_image)[0])
        r_SB_1 = (SB_profile(flux_SB_list[1], center, gridspace = 'log', grids = 30, if_annuli = if_annuli,radius= radi)[0])
    #    r_mag_diff = 2.5 * np.log10(SB_profile(flux_SB_list[1], center, gridspace = 'log', radius= radi)[0])
        ax5.plot(r_grids*delatPixel, (r_SB_0- r_SB_1), 'ro')
        ax5.set_ylabel('$\Delta SB$', fontsize=15)
        ax5.set_xlabel('arcsec', fontsize=15)
        ax5.set_xscale('log')
        ax5.set_xticks([0.1, 0.2, 0.5, 1, 2])
        import matplotlib
        ax5.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax5.plot(x, y, 'k--')
        ax5.yaxis.set_label_position('right')
        ax5.yaxis.tick_right()
        ax5.yaxis.set_ticks_position('both')
        ax5.set_xlim([(r_grids*delatPixel).min()*0.85, (r_grids.max()+6)*delatPixel])
        plt.ylim([-5,5])
        ax5.set_yticks([-5,-2.5, 0., 2.5])
    pos4_o = ax4.get_position() # get the original position
    pos4 = [pos4_o.x0+0.085, pos4_o.y0 + 0.10, pos4_o.width*0.72, pos4_o.height*0.8]
    ax4.set_position(pos4) # set a new position
    pos5_o = ax5.get_position() # get the original position
    pos5 = [pos5_o.x0+0.085, pos5_o.y0+0.08, pos5_o.width*0.72, pos5_o.height*1.1]
    ax5.set_position(pos5) # set a new position
    if plot_compare == True:
        plt.show()
    else:
        plt.close()
    return f

def make_ticklabels_invisible(fig):
    for i, ax in enumerate(fig.axes):
#        ax.text(0.5, 0.5, "ax%d" % (i+1), va="center", ha="center")
        if i !=3 and i !=4:
            for tl in ax.get_xticklabels() + ax.get_yticklabels():
                tl.set_visible(False)

def coordinate_arrows(ax, d, color='black', arrow_size=0.02):
    d0 = d / 12.
    p0 = d / 12.
    pt = d / 7.
    deltaPix = 1
    ra0, dec0 = (d - d0) / deltaPix, d0 / deltaPix
    xx_, yy_ = ra0, dec0
    xx_ra, yy_ra = (ra0 - p0, dec0)
    xx_dec, yy_dec = (ra0, dec0 + p0)
    xx_ra_t, yy_ra_t = (ra0 - pt, dec0)
    xx_dec_t, yy_dec_t = (ra0, dec0 + pt)

    ax.arrow(xx_ * deltaPix, yy_ * deltaPix, (xx_ra - xx_) * deltaPix, (yy_ra - yy_) * deltaPix,
             head_width=arrow_size * d, head_length=arrow_size * d, fc=color, ec=color, linewidth=1.2)
    ax.text(xx_ra_t * deltaPix, yy_ra_t * deltaPix, "E", color=color, fontsize=12, ha='center')
    ax.arrow(xx_ * deltaPix, yy_ * deltaPix, (xx_dec - xx_) * deltaPix, (yy_dec - yy_) * deltaPix,
             head_width=arrow_size * d, head_length=arrow_size * d, fc
             =color, ec=color, linewidth=1.2)
    ax.text(xx_dec_t * deltaPix, yy_dec_t * deltaPix, "N", color=color, fontsize=12, ha='center')

def scale_bar(ax, d, dist=1/0.13, text='1"', color='black', flipped=False):
    if flipped:
        p0 = d - d / 15. - dist
        p1 = d / 15.
        ax.plot([p0, p0 + dist], [p1, p1], linewidth=2, color=color)
        ax.text(p0 + dist / 2., p1 + 0.02 * d, text, fontsize=15, color=color, ha='center')
    else:
        p0 = d / 15.
        ax.plot([p0, p0 + dist], [p0, p0], linewidth=2, color=color)
        ax.text(p0 + dist / 2., p0 + 0.02 * d, text, fontsize=15, color=color, ha='center')

def cal_bkg_var(sub, img, mask_list=None, test_n = 20):
    center = np.asarray(img.shape) /2
    t_r_flux, t_r_grids, regions=flux_profile(img-sub, center,
                                          radius=center.min(), grids=50,
                                          ifplot=False, fits_plot= False,
                                          mask_list=mask_list)
    value = (t_r_flux[-test_n:][1:]-t_r_flux[-test_n:][:-1]).sum()
    return value

#print cal_bkg_var(sub=0.0017, img=img, mask_list=mask_list)

def min_sub(ini, img, mask_list=None, test_n=20):
    sub_list1=np.linspace(ini-0.01,ini+0.01,20)
    values = np.zeros_like(sub_list1)
    for i in range(len(sub_list1)):
        values[i] = cal_bkg_var(sub=sub_list1[i], img=img, mask_list=mask_list)
    idx = np.where(abs(values).min() == abs(values))[0][0]
    ini2 = sub_list1[idx]
    sub_list2=np.linspace(ini2-0.001,ini2+0.001,40)
    values2 = np.zeros_like(sub_list2)
    for i in range(len(sub_list2)):
        values2[i] = cal_bkg_var(sub=sub_list2[i], img=img, mask_list=mask_list, test_n=test_n)
    idx2 = np.where(abs(values2).min() == abs(values2))[0][0]
    return sub_list2[idx2], values2[idx2]

def gaussian(x, mean, amplitude, standard_deviation):
    return amplitude * np.exp( - ((x - mean) / standard_deviation) ** 2)

def fit_data_as_gaussian(data, ifplot=1):
    bin_heights, bin_borders = np.histogram(data, bins='auto')
    bin_widths = np.diff(bin_borders)
    bin_centers = bin_borders[:-1] + bin_widths / 2
    plt.bar(bin_centers, bin_heights, width=bin_widths, label='histogram')
    from scipy.optimize import curve_fit
    popt, pcov = curve_fit(gaussian, bin_centers, bin_heights, p0=[0., bin_heights.max(), 0.01])
    x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 10000)
    gauss_grid = gaussian(x_interval_for_fit, *popt)
    plt.plot(x_interval_for_fit, gauss_grid, label='fit',c='red')
    fit_center_idx = np.where(gauss_grid==gauss_grid.max())[0][0]
    line = np.linspace(0,10000,10)
    peak_loc = x_interval_for_fit[fit_center_idx]
    plt.plot(peak_loc*np.ones_like(line), line, 'black')
    plt.plot(peak_loc*np.ones_like(line) + popt[2], line, 'black')
    plt.plot(peak_loc*np.ones_like(line) - popt[2], line, 'black')
    plt.ylim((0, bin_heights.max()*5./4.))
    plt.legend()
    if ifplot == 1:
        plt.show()
    else:
        plt.close()
    return peak_loc, popt[2]

def fit_bkg_as_gaussian(data, thre=None, ifplot=1, return_sig=0):
    if thre is None:
        thre = np.max(data)
    bin_heights, bin_borders = np.histogram(data[abs(data)<thre], bins='auto')
    bin_widths = np.diff(bin_borders)
    bin_centers = bin_borders[:-1] + bin_widths / 2
    plt.bar(bin_centers, bin_heights, width=bin_widths, label='histogram')
    from scipy.optimize import curve_fit
    popt, _ = curve_fit(gaussian, bin_centers, bin_heights, p0=[1., 0., 1.])
#    print popt
    x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 10000)
    gauss_grid = gaussian(x_interval_for_fit, *popt)
    plt.plot(x_interval_for_fit, gauss_grid, label='fit',c='red')
    fit_center_idx = np.where(gauss_grid==gauss_grid.max())[0][0]
    line = np.linspace(0,10000,10)
    QSO_peak_loc = x_interval_for_fit[fit_center_idx]
    plt.plot(QSO_peak_loc*np.ones_like(line), line, 'black')
    plt.ylim((0, bin_heights.max()*5./4.))
    plt.legend()
    if ifplot == 1:
        plt.show()
    else:
        plt.close()
    if return_sig == 0:
        return QSO_peak_loc
    elif return_sig == 1:
        return QSO_peak_loc, popt[2]

max_stars_repo_name: jasminwash/hostgal

max_issues_repo_name: jasminwash/hostgal

max_forks_repo_name: jasminwash/hostgal

reason: import numpy,from scipy,import astropy,from astropy


text:
__author__ = "Katri Leino"
__copyright__ = "Copyright (c) 2018, Aalto Speech Research"

import pybrain

import numpy as np
from scipy import *
import math
import copy

class UI():
    '''
    INPUT:
        NumStates : The number of the states
        Actionmatrix : A list of allowed actions (states which have not yet been visited)
        NumActions : The number of actions
        Goal : Has goal been achieved True/False
        SensorErrors : Sensor error percent for each modality
        ConfusionError : Sensor error percent for each modality
        Penalties : KLM penalties (in case 3 not needed because penalty is calculated from the distance)
        ItemGrid : x- and y-coordinates of the items [x_coordinates, y_coordinates]
        ItemDimensions : Dimension of items [Width, Height]
        InitPos : User's initial position [x coordinate, y coordinate]
        GoalPos : The last position user needs to press [x coordinate, y coordinate]
    '''
    def __init__(self, NumStates, Actionmatrix, NumActions, Goal, SensorErrors, ConfusionError, Penalties, ItemGrid, ItemDimensions, InitPos, GoalPos):
    	self.initState = 0
    	self.current_state = self.initState
        self.actionmatrix = Actionmatrix # allowed actions
        self.num_of_states = NumStates
        self.num_of_actions = NumActions
        self.visited_states = [0]*self.num_of_actions
        self.visited_states[self.num_of_actions-1] = 1
        self.used_actions = [0]*self.num_of_actions

        self.goal = Goal # True/False
        self.sensor_errors = SensorErrors
        self.confusion_error = ConfusionError
        self.penalties = Penalties
        self.sys_delay = 0 # System delay

        self.grid = ItemGrid
        self.dimensions = ItemDimensions
        self.initial_position = InitPos
        self.goal_position = GoalPos

        self.maxSteps = 30 # Maximum steps before the episode stops
        self.steps = 0

    def setGoal(self, Goal):
        self.goal = int(Goal)

    def setInitialState(self, initial_state):
        self.initState = initial_state
        self.current_state = initial_state

    def setModalities(self, Buttons):
        print "setModalities: Not available"

    def getSensors(self):
        """ the currently visible state of the world / Current state """
        return [float(self.current_state),]

    def getGoal(self):
    	return self.goal

    def performAction(self, action):
        # Performs given action and changes state
        self.steps += 1

        # Error model ::: NOT USED IN CASE3
        #modality = 0 # for remote controller there is only a single modality
        #action = self.errorModel(action, modality)
        #if action == -1: return

        self.current_state = action+1
        self.visited_states[action] = 1

        if np.sum(self.visited_states) == self.num_of_actions:
            self.visited_states[self.num_of_actions-1] = 0
            self.goal = True

    # INPUT:
        # action : performed action
        # modality : modality action belongs to
    def errorModel(self, action, modality):
        # User's error
        if random.uniform(0,1) < 0.04:
            return -1

        # Recognition error: Input not recognized
        if random.uniform(0,1) < self.sensor_errors[modality]:
            return -1

        # Confusion error: Confuses actions within the same modality.
        if random.uniform(0,1) < self.confusion_error[modality]:
            action_list = copy.deepcopy(self.mods[self.current_state])
            action_list.remove(float(action)) # Remove original action from the list
            if len(action_list) > 0:
                act = int(random.choice(action_list))
                return act

        return action

    def reset(self):
        self.current_state = self.initState
        self.visited_states = [0]*self.num_of_actions
        self.visited_states[self.num_of_actions-1] = 1
        self.used_actions = [0]*self.num_of_actions
        self.goal = False

    def getPenalty(self, action, prev_action):
        # Compute penalty
        penalty = 0

        if int(self.used_actions[action]) > 1: # Uses same action more than once. Should not be possible
            return 15

        # Calculating distance between items.
        if prev_action == -1: #Previous state was the initial state
            distance_x = np.abs(self.grid[0][int(action)]-self.initial_position[0])
            distance_y = np.abs(self.grid[1][int(action)]-self.initial_position[1])
            distance = math.sqrt(distance_x**2 + distance_y**2)
        elif np.sum(self.visited_states) == self.num_of_states: # Last command
            distance_x = np.abs(self.grid[0][int(action)]-self.goal_position[0])
            distance_y = np.abs(self.grid[1][int(action)]-self.goal_position[1])
            distance = math.sqrt(distance_x**2 + distance_y**2)
            penalty += 5
        else:
            distance_x = np.abs(self.grid[0][int(action)]-self.grid[0][int(prev_action)])
            distance_y = np.abs(self.grid[1][int(action)]-self.grid[1][int(prev_action)])
            distance = math.sqrt(distance_x**2 + distance_y**2)

        if distance_x== 0: # Approach from above or below
            width = self.dimensions[1]
        elif distance_y == 0: # Approach from either side
            width = self.dimensions[0]
        elif math.tan(distance_y/distance_x) < math.pi/4: # Approach from side
            width = self.dimensions[0]
        else:
            width = self.dimensions[1]

        #Fitts' Law ID value
        penalty += math.log((2.*distance/width), 2)
        penalty += self.sys_delay

        return penalty

# Defines the task
class UITask():
    discount = True

    # INPUT
        # environment : UI class object. Defines the environment

    def __init__(self, environment):
        """ All tasks are coupled to an environment. """
        self.env = environment
        self.cumreward = 0 # tracking cumulative reward
        self.samples = 0 # tracking the number of samples

        self.sensor_limits = None
        self.actor_limits = None
        self.clipping = True

        self.current_action = 0 # Saving current action
        self.prev_action = -1 # Saving previous action

    def reset(self):
        self.env.reset()
        self.cumreward = 0
        self.samples = 0
        self.prev_action = -1

    def getObservation(self):
        """ A filtered mapping to getSensors of the underlying environment. """
        sensors = self.env.getSensors()
        if self.sensor_limits:
            sensors = self.normalize(sensors)
        return sensors

    def isFinished(self):
        if self.samples >= self.env.maxSteps:
            return True

        # action list is empty (all states have been visited)
        if np.sum(self.env.visited_states) == self.env.num_of_actions:
            return True

        if self.env.goal:
            return True

        else: return False

    def performAction(self, action):
        act = int(action)
        self.current_action = act
        self.env.performAction(act)
        self.samples += 1

    def getReward(self):
        reward = -1*self.env.getPenalty(self.current_action, self.prev_action)
        return reward

    def addReward(self):
        if self.env.current_state == -1: return

        discount_value = 0.9
        if self.discount:
            self.cumreward += power(discount_value, self.samples) * self.getReward()
        else:
            self.cumreward += self.getReward()

    def getTotalReward(self):
        """ Return the accumulated reward since the start of the episode """
        return self.cumreward

ext: py

max_stars_repo_name: aalto-speech/rl-klm

max_issues_repo_name: aalto-speech/rl-klm

max_forks_repo_name: aalto-speech/rl-klm

reason: import numpy,from scipy


text:
#!/usr/bin/env python3

# This script produce appropriate SampleSheets for all samples requiring it from Sequencing_summary.xlsx
#if not os.path.isfile("../mw-sst/Sequencing_summary.csv"):

# Modified 03/09/2020
# if the data are single-cell RNA-seq bcl, we must not add 'I5_Index_ID' and 'index2' in the samplesheet or cellranger mkfastq throws an error.

# Modified 09/10/2020
# Merged scrna_bcl and bcl in the sequencing_summary after discussion with biologists.

import pandas
import os
import numpy as np
import hashlib
from pathlib import Path

if os.path.isfile("../mw-Madaci2021/Sequencing_summary.xlsx"):
    #def make_sample_sheets_from_sequencing_summary(seq_sum_path:
    samples = pandas.read_excel("../mw-Madaci2021/Sequencing_summary.xlsx", sheet_name="samples", engine="openpyxl")
    # Old line to get only line with bcl_Nextseq500 origin
    #samples_from_bcl = samples[samples['origin'] == 'bcl_NextSeq500']
    #samples_from_bcl = samples[samples['origin'] in ['bcl', 'bcl_no_mismatch']]
    #samples_from_bcl = samples[samples['origin'] == 'bcl']
    samples_from_bcl = samples[samples['origin'].isin(['bcl', 'bcl_no_mismatch', 'bcl_with_fastq_index'])]
    #samples_from_scrna_bcl = samples[samples['origin'] == 'scrna_bcl']

    #all_samples = [samples_from_bcl, samples_from_scrna_bcl]
    #all_samples_df = pandas.concat(all_samples)

    # Add adapter file reading and dictionnary creation to add adapter sequences in the samplesheet.
    if os.path.isfile("../mw-lib/src/snakemake/tables/adapter.tsv"):
        adapter_dict = {}
        with open("../mw-lib/src/snakemake/tables/adapter.tsv", "r") as adapter_file:
            # Skip first line containing header
            next(adapter_file)
            # Fill the dictionnary with all kits information:   Kit Adapter1    Adapter2
            for line in adapter_file:
                splitted_line = line.split("\t")
                adapter_dict[splitted_line[0]] = [splitted_line[1], splitted_line[2].replace("\n", "")]

    #experiments_from_bcl = all_samples_df.accession.unique()
    experiments_from_bcl = samples_from_bcl.accession.unique()

    mwconf['bcl2fastq_targets'] = []

    for experiment in experiments_from_bcl:
        # Variable to get the loop dataframe
        current_df = samples_from_bcl[samples['accession'] == experiment]
        kit_used = current_df.Kit.unique()[0]

        # Test to create symbolik link from bcl folder to out/ln/alias/sst/by_run/runX
        #run = str(int(current_df.run.unique()[0]))
        #print("out/ln/alias/sst/by_run/run" + run)
        #mwconf['bcl2fastq_targets'].append("out/ln/alias/sst/by_run/run" + run)

        #print(current_df.origin.unique())
        #print(current_df)

        # Condition to determine wether data are single or double indexed
        #if current_df[samples['index2']] != "-":
        #   current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','I5_Index_ID','index2','Sample_Project','Description']]
        #   data = current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','I5_Index_ID','index2','Sample_Project','Description']]
        #else:
        #   data = current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','Sample_Project','Description']]
        #current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','I5_Index_ID','index2','Sample_Project','Description']]
        # Special case for scrna cause cellranger seems to detect the second index even if it is empty.
        if(['scRNA', 'scRNA_HTO'] in current_df.type.unique()):
            #bcl2fastq_target = "out/cellranger/mkfastq/" + experiment + "/Run_" + str(int(current_df[['run']].iloc[0,0])) + "/outs/fastq_path/" + str(current_df[['Sample_Project']].iloc[0,0])  + "/Reports/html/tree.html"
            bcl2fastq_prefix = "out/cellranger/mkfastq/" + experiment
            bcl2fastq_target = bcl2fastq_prefix + "/Reports/html/tree.html"
            data_tmp = current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','Index_10X','I7_Index_ID','index', 'I5_Index_ID', 'index2', 'Sample_Project','Description']].fillna('')
            # 2021 04 07 Update to enable dual indexing demultiplexing. Cellranger can read index name and will replace it by the index sequence. It also detect dual index or single index.
            # 2021-04-23 Update the conditions to detect if samples have double index or not.
            for index, row in data_tmp.iterrows():
                if row['Index_10X'] == "yes":
                    mRNA_index=index
                    index_name=row["I7_Index_ID"]
                    data_tmp.at[mRNA_index, "index"] = index_name
                    if(str(row['I5_Index_ID']) != ''):
                        data_tmp.at[mRNA_index, "I5_Index_ID"] = index_name
                        data_tmp.at[mRNA_index, "index2"] = index_name
                        data_dropped = data_tmp.drop(columns=["Index_10X"])
                    else:
                        data_dropped = data_tmp.drop(columns=["Index_10X", "I5_Index_ID", "index2"])
                #data_dropped = data_tmp.drop(columns=["Index_10X"])
                elif(str(row['I5_Index_ID']) == ''):
                    data_dropped = data_tmp.drop(columns=["Index_10X", "I5_Index_ID", "index2"])
                else:
                    data_dropped = data_tmp.drop(columns=["Index_10X"])
            data = data_dropped
        elif(['bcl_no_mismatch'] in current_df.origin.unique()):
            bcl2fastq_prefix = "out/bcl2fastq/_--no-lane-splitting_--barcode-mismatches_0/" + experiment
            bcl2fastq_target = bcl2fastq_prefix + "/Reports/html/tree.html"
            data = current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','I5_Index_ID','index2','Sample_Project','Description']]
        elif(['bcl_index_generation'] in current_df.origin.unique()):
            bcl2fastq_prefix = "out/bcl2fastq/_--no-lane-splitting_--create-fastq-for-index-reads/" + experiment
            #bcl2fastq_prefix = "out/bcl2fastq/_--no-lane-splitting_--create-fastq-for-index-reads_--ignore-missing-bcls/" + experiment
            # uncomment this line to reproduce the error of run379
            #bcl2fastq_prefix = "out/bcl2fastq/_--no-lane-splitting_--create-fastq-for-index-reads/" + experiment
            bcl2fastq_target = bcl2fastq_prefix + "/Reports/html/tree.html"
            data = current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','I5_Index_ID','index2','Sample_Project','Description']]
        else:
            bcl2fastq_prefix = "out/bcl2fastq/_--no-lane-splitting/" + experiment
            bcl2fastq_target = bcl2fastq_prefix + "/Reports/html/tree.html"
            data = current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','I5_Index_ID','index2','Sample_Project','Description']]
        #data = current_df[['Sample_ID','Sample_Name','Sample_Plate','Sample_Well','I7_Index_ID','index','Sample_Project','Description']]
        #bcl2fastq_target = "out/bcl2fastq/_--no-lane-splitting/" + experiment + "/Reports/html/tree.html"

        #data = data.sort_values(by='Sample_ID')
        # Sort by index to prevent trouble with mix of ID and numbers
        data = data.sort_index()
        #filename = experiment + "/SampleSheet.csv"
        filename = bcl2fastq_prefix + "/SampleSheet.csv"
        tmp_filename = bcl2fastq_prefix + "/SampleSheet.csv.tmp"

        filename = filename.replace('//', '/')
        tmp_filename = tmp_filename.replace('//', '/')

        # SampleSheet should be created only if it does not exist
        # to ensure bcl2fastq is not redone only because the SampleSheet is touched
        # by the process.
        # Maybe replace this with a hash check to still write a new SampleSheet if the
        # infos from Sequencing_summary are different.
        #if not os.path.isfile(filename):
        os.makedirs(bcl2fastq_prefix, exist_ok = True)
        header=('[Header]\n'
                'IEMFileVersion,4\n'
                'Investigator Name,'+ current_df[['customer']].iloc[0,0]+'\n'
                'Experiment Name,'+ current_df[['Sample_Project']].iloc[0,0]+'\n'
                'Date,07/05/2020\n'
                'Workflow,GenerateFASTQ\n'
                'Application,NextSeq FASTQ Only\n'
                'Description,' + current_df[['Description']].iloc[0,0]+'\n'
                'Chemistry,Default\n\n'
                '[Reads]\n'
                #+current_df[['Reads']].iloc[0,0].replace('x','\n')+'\n\n'
                # Add str convertion to the read value in order to prevent error is SE (i.e x is absent)
                +str(current_df[['Reads']].iloc[0,0]).replace('x','\n')+'\n\n')
                # Add settings section containing adapter information
                #'[Settings]\n'
                #'Adapter,' + adapter_dict[kit_used][0] + '\n'
                #'AdapterRead2,' + adapter_dict[kit_used][1] + '\n\n'
                #'[Data]\n')

        if(['scRNA', 'scRNA_HTO'] not in current_df.type.unique()):
            adapter_content=('[Settings]\n'
         'Adapter,' + adapter_dict[kit_used][0] + '\n'
         'AdapterRead2,' + adapter_dict[kit_used][1] + '\n\n'
         '[Data]\n')
        else:
            adapter_content=("\n"
            '[Data]\n')
        if current_df[['process']].iloc[0,0] == "yes":
            myfile = open(filename, 'w')
            myfile.write(header)
            myfile.write(adapter_content)
            myfile.write(data.to_csv(index=False))
            myfile.close()

            #print(str(os.path.isfile(filename)) + " SampleSheet.csv " + filename)
            #print(str(os.path.isfile(tmp_filename)) + " SampleSheet.csv.tmp " + tmp_filename )
            #if(os.path.isfile(filename)):
                #print(hashlib.md5(open(filename,'rb').read()).hexdigest())
                #print(hashlib.md5(open(tmp_filename,'rb').read()).hexdigest())
                #print(hashlib.md5(open(filename,'rb').read()).hexdigest() == hashlib.md5(open(tmp_filename,'rb').read()).hexdigest())
                #print(hashlib.md5(open(filename,'rb').read()).hexdigest() != hashlib.md5(open(tmp_filename,'rb').read()).hexdigest())
                #if(hashlib.md5(open(filename,'rb').read()).hexdigest() != hashlib.md5(open(tmp_filename,'rb').read()).hexdigest()):
                    #os.replace(filename, tmp_filename)
                #else:
                    #if(os.path.isfile(tmp_filename)):
                        #print(tmp_filename)
                        #print(os.path.isfile(tmp_filename))
                        #os.remove(tmp_filename)

            #else:
                #if(os.path.isfile(tmp_filename)):
                    #os.rename(tmp_filename, filename)

    # Output line before the no lane splitting option
        #bcl2fastq_target = "out/bcl2fastq/_/" + experiment + "/Reports/html/tree.html"

        # Users may add '/' before and after accession in Sequencing_summary.xlsx
        # This leads to such Snakemake error:
        # File path XXX contains double '/'. This is likely unintended. It can also lead to inconsistent results
        # To avoid id, double slashes "//" are replaced by single one.
        bcl2fastq_target = bcl2fastq_target.replace('//', '/')

        # Add a condition to read the process column of the summary.
        # Sometimes, we may not want to run bcl2fastq.
        if current_df[['process']].iloc[0,0] == "yes":
            mwconf['bcl2fastq_targets'].append(bcl2fastq_target)

reason: import numpy


text:
import numpy as np

#create vector
vec = np.arange(9)

#reshape vector in 3x3 matrix
matrix = vec.reshape((3,3))
print(matrix)

#copy matrix
matrix2 = matrix.copy()
#change first column, first entry of matrix2
matrix2[:, 0][0] = 9
print(matrix)
print(matrix2)

#some sample for concatenate, vstack/hstack and split (vsplit, hsplit)
#those methods can also be found on numpy.org -> manual

#copy first column as new vector
vec2 = matrix[:, 0].copy()
print(vec2)

#concatenate on axis 0
print(np.concatenate([vec, vec2], axis=0))

#concatenate maxtrix and matrix2 on axis=1
print(np.concatenate([matrix, matrix2], axis=1))

# horizontally stack the arrays
y = np.array([[99],
              [99],
              [99]])
print(np.hstack([matrix, y]))

#finally we split the vector
x1, x2 = np.split(vec, [3])
print(x1)
print(x2)

#and the matrix
y1, y2 = np.vsplit(matrix, [2])
print(y1)
print(y2)


text:
# Light-Head R-CNN: In Defense of Two-Stage Object Detector  http://arxiv.org/abs/1711.07264
import chainer
import chainer.links as L
import chainer.functions as F
from chainer.links.model.vision.resnet import ResNet50Layers, BuildingBlock, _global_average_pooling_2d
import numpy as np
import copy
from chainer_maskrcnn.functions.roi_align_2d_yx import _roi_align_2d_yx

class LightRoIMaskHead(chainer.Chain):
    mask_size = 14

    def __init__(self,
                 n_class,
                 roi_size,
                 loc_initialW=None,
                 score_initialW=None,
                 mask_initialW=None):
        # n_class includes the background
        super().__init__()
        with self.init_scope():
            # Separable Convolution Layers 変数の名前はpaperに準拠してみた
            k = 15
            C_mid = 256
            C_out = 490
            # 初期値ワカラン・・ chainer.initializers.Normal(0.001)とかの方がいいか？
            # レイヤーの名前は論文の図を見たときに up left, bottom left, up right, bottom rightの４つw
            p = int(k / 2)
            self.conv_ul = L.Convolution2D(
                in_channels=None, out_channels=C_mid, ksize=(k, 1), pad=(p, 0))
            self.conv_bl = L.Convolution2D(
                in_channels=C_mid,
                out_channels=C_out,
                ksize=(1, k),
                pad=(0, p))
            self.conv_ur = L.Convolution2D(
                in_channels=None, out_channels=C_mid, ksize=(1, k), pad=(0, p))
            self.conv_br = L.Convolution2D(
                in_channels=C_mid,
                out_channels=C_out,
                ksize=(k, 1),
                pad=(p, 0))
            self.fc = L.Linear(None, 2048)
            self.cls_loc = L.Linear(2048, 4, initialW=loc_initialW)
            self.score = L.Linear(2048, n_class, initialW=score_initialW)
            # マスク推定ブランチへ
            self.conv2 = L.Convolution2D(
                in_channels=None,
                out_channels=256,
                ksize=3,
                stride=1,
                pad=1,
                initialW=mask_initialW)
            self.conv3_ = L.Convolution2D(
                in_channels=None,
                out_channels=256,
                ksize=3,
                stride=1,
                pad=1,
                initialW=mask_initialW)
            self.conv4 = L.Convolution2D(
                in_channels=None,
                out_channels=256,
                ksize=3,
                stride=1,
                pad=1,
                initialW=mask_initialW)
            self.deconv1_ = L.Deconvolution2D(
                in_channels=None,
                out_channels=n_class - 1,
                ksize=2,
                stride=2,
                pad=0,
                initialW=mask_initialW)

        self.n_class = n_class
        self.roi_size = roi_size

    def __call__(self, x, rois, roi_indices, spatial_scale):
        roi_indices = roi_indices.astype(np.float32)
        indices_and_rois = self.xp.concatenate(
            (roi_indices[:, None], rois), axis=1)

        # roi poolingをする前に、thin feature mapに変換します
        # activationしないっぽいことが書いてあるんだよなー
        left_path = self.conv_bl(self.conv_ul(x))
        right_path = self.conv_br(self.conv_ur(x))
        tfp = left_path + right_path

        pool = _roi_align_2d_yx(tfp, indices_and_rois, self.roi_size,
                                self.roi_size, spatial_scale)

        h = F.relu(self.fc(pool))
        roi_cls_locs = self.cls_loc(h)
        roi_scores = self.score(h)
        # at prediction time, we use two pass method.
        # at first path, we predict box location and class
        # at second path, we predict mask with accurate location from first path
        if chainer.config.train:
            mask = F.relu(self.conv2(pool))
            mask = F.relu(self.conv3_(mask))
            mask = F.relu(self.conv4(mask))
            mask = self.deconv1_(pool)
            #mask = self.conv2(self.deconv1(pool))
            return roi_cls_locs, roi_scores, mask
        else:
            # cache tfp for second path
            self.tfp = tfp
            return roi_cls_locs, roi_scores

    def predict_mask(self, rois, roi_indices, spatial_scale):
        roi_indices = roi_indices.astype(np.float32)
        indices_and_rois = self.xp.concatenate(
            (roi_indices[:, None], rois), axis=1)
        pool = _roi_align_2d_yx(self.tfp, indices_and_rois, self.roi_size,
                                self.roi_size, spatial_scale)

        mask = F.relu(self.conv2(pool))
        mask = F.relu(self.conv3_(mask))
        mask = F.relu(self.conv4(mask))
        mask = self.deconv1_(pool)
#        mask = self.deconv1(pool)
#        mask = self.conv2_(mask)
#        mask = self.conv3(mask)
#        mask = self.conv2(self.deconv1(pool))
        return mask

max_stars_repo_name: katotetsuro/chainer-maskrcnn

max_issues_repo_name: katotetsuro/chainer-maskrcnn

max_forks_repo_name: katotetsuro/chainer-maskrcnn