text:
import math

import numpy as np
from sklearn.cluster import MiniBatchKMeans

from probgf.methods_simple import SpatioTemporalBase

class NSPI(SpatioTemporalBase):
    """
    Implementation of the Neighborhood Similar Pixel Interpolator
    Most code is taken from the authors code ("NSPI_time_series_Python-updated_20200312.zip")
    More info:
    Jin Chena, Xiaolin Zhua, James E. Vogelmannc, Feng Gaod, and Suming Jinf
    "A simple and effective method forfilling gaps in Landsat ETM+ SLC-off images"
    (2011)
    https://xiaolinzhu.weebly.com/open-source-code.html
    https://doi.org/10.1016/j.rse.2010.12.010
    """

    @classmethod
    def default_config(cls):
        return '30,4,0.0,10000.0'

    def __init__(self, config, dates, console, emiters):
        super().__init__(config, dates, console, emiters)
        self.slice_shape = 'img'

    def name(self):
        return 'NSPI_MINSIM{}_NCLASS{}_DNMIN{}_DNMAX{}'.format(
            self.min_similar,
            self.num_class,
            self.DN_min,
            self.DN_max,
        ).replace('.', '_')

    def configure(self, config):
        "instead of parsing the parameters from the yaml config file, we here allow to set parameters via command line argument"
        parts = config.split(',')
        try:
            self.min_similar = int(parts[0])  # set the minimum sample size of similar pixels, 20~40 recommended
            self.num_class   = int(parts[1])  # set the number of classes (given empirically)
            # num_band = int(parts[2] not needed
            self.DN_min      = float(parts[2])  # set the range of DN value of the image,If byte, 0 and 255
            self.DN_max      = float(parts[3])
            # patch_long = param['patch_long'] not needed
            # DOY = param['DOY'] not needed
        except (ValueError, IndexError):
            raise RuntimeError('Invalid config "{}".\n{} needs a comma-seperated list of the following values:\n'.format(config, self.__class__.method_id()) + \
            '  min_similar : minimum sample size of similar pixels, 20~40 recommended (int)\n' + \
            '  num_class   : number of classes (given empirically) (int)\n' + \
            '  DN_min      : range min of DN value of the image (float)\n' + \
            '  DN_max      : range max of DN value of the image (float)')

    def run_training(self, data, obs, split, progr_train):
        super().run_training(data, obs, split, progr_train)

    def run_prediction(self, to_pred, obs, split, progr_pred):
        """
        NSPI logic, taken from the original "NSPI_TIMESERIES.py" script (removed boilerplate code with comments)
        Implementation was also changed to not subprocess smaller patches of the data, but instead predict the whole thing
        """
        # reshape to NSPI format+
        orig_shape = to_pred.shape
        mask = np.moveaxis(obs[:, :, :, 0], 2, 0)
        mask = np.invert(mask)
        to_pred = np.moveaxis(to_pred.reshape(obs.shape[0], obs.shape[1], obs.shape[2] * obs.shape[3]), 2, 0)

        # ******************************************************************************************************
        # *******************************Set parameters and read input data*************************************

        # starttime = datetime.datetime.now()  # initial time

        # root = tk.Tk()
        # root.withdraw()

        # Please set the appropriate parameters for your own data in parameters_timeseries.yaml before open this file!
        # f = open(filedialog.askopenfilename(title=u"Open the parameter settings file:"))
        # param = yaml.safe_load(f)
        min_similar = self.min_similar # param['min_similar']  # set the minimum sample size of similar pixels, 20~40 recommended
        num_class = self.num_class # param['num_class']  # set the number of classes (given empirically)
        num_band = to_pred.shape[0] // mask.shape[0] # param['num_band']  # set the number of bands in each image
        DN_min = self.DN_min # param['DN_min']  # set the range of DN value of the image,If byte, 0 and 255
        DN_max = self.DN_max # param['DN_max']
        DOY = self.dates # param['DOY']  # Day of year of each image in the time-series

        # # set the folder for storing temporary output results
        # temp_file = filedialog.askdirectory(title=u"Set the temporary folder")

        # # open the SLC-off ETM+ time-series image
        # path1 = filedialog.askopenfilename(title=u"Open the SLC_off ETM+ time-series image:")
        # suffix = os.path.splitext(path1)[-1]
        # nl, ns, data_off = read_raster(path1)
        # orig_ns = ns
        # orig_nl = nl
        # fp = gdal.Open(path1)
        # nb = fp.RasterCount

        # # divide the whole scene into blocks
        # n_nl = math.ceil(nl / patch_long)
        # n_ns = math.ceil(ns / patch_long)
        # n_image = math.ceil(nb / num_band)
        n_image = to_pred.shape[0] // num_band

        # ind_patch1 = np.zeros((n_nl * n_ns, 4), dtype=np.int)
        # ind_patch = np.zeros((n_nl * n_ns, 4), dtype=np.int)
        # location = np.zeros((n_nl * n_ns, 4), dtype=np.int)

        # for i_ns in range(0, n_ns):
        #     for i_nl in range(0, n_nl):
        #         ind_patch1[n_ns * i_nl + i_ns, 0] = i_ns * patch_long
        #         ind_patch[n_ns * i_nl + i_ns, 0] = np.max([0, ind_patch1[n_ns * i_nl + i_ns, 0] - 10])
        #         location[n_ns * i_nl + i_ns, 0] = ind_patch1[n_ns * i_nl + i_ns, 0] - ind_patch[n_ns * i_nl + i_ns, 0]

        #         ind_patch1[n_ns * i_nl + i_ns, 1] = np.min([ns - 1, (i_ns + 1) * patch_long - 1])
        #         ind_patch[n_ns * i_nl + i_ns, 1] = np.min([ns - 1, ind_patch1[n_ns * i_nl + i_ns, 1] + 10])
        #         location[n_ns * i_nl + i_ns, 1] = ind_patch1[n_ns * i_nl + i_ns, 1] - ind_patch1[n_ns * i_nl + i_ns, 0] + location[n_ns * i_nl + i_ns, 0]

        #         ind_patch1[n_ns * i_nl + i_ns, 2] = i_nl * patch_long
        #         ind_patch[n_ns * i_nl + i_ns, 2] = np.max([0, ind_patch1[n_ns * i_nl + i_ns, 2] - 10])
        #         location[n_ns * i_nl + i_ns, 2] = ind_patch1[n_ns * i_nl + i_ns, 2] - ind_patch[n_ns * i_nl + i_ns, 2]

        #         ind_patch1[n_ns * i_nl + i_ns, 3] = np.min([nl - 1, (i_nl + 1) * patch_long - 1])
        #         ind_patch[n_ns * i_nl + i_ns, 3] = np.min([nl - 1, ind_patch1[n_ns * i_nl + i_ns, 3] + 10])
        #         location[n_ns * i_nl + i_ns, 3] = ind_patch1[n_ns * i_nl + i_ns, 3] - ind_patch1[ n_ns * i_nl + i_ns, 2] + location[n_ns * i_nl + i_ns, 2]

        # tempoutname = temp_file + '\\temp_target'

        # for isub in range(0, n_nl * n_ns):
        #     col1 = ind_patch[isub, 0]
        #     col2 = ind_patch[isub, 1]
        #     row1 = ind_patch[isub, 2]
        #     row2 = ind_patch[isub, 3]
        #     data = data_off[:, row1:row2 + 1, col1:col2 + 1]
        #     out_name = tempoutname + str(isub + 1) + suffix
        #     fp = path1
        #     writeimage(data, out_name, fp)

        # # open the mask of SLC-off ETM+ image: 0 is good pixels, 1 gaps 2 shadows 3 clouds
        # path2 = filedialog.askopenfilename(title=u"Open the mask image):")
        # _, _, data_gap = read_raster(path2)

        # tempoutname = temp_file + '\\temp_mask'

        # for isub in range(0, n_nl * n_ns):
        #     col1 = ind_patch[isub, 0]
        #     col2 = ind_patch[isub, 1]
        #     row1 = ind_patch[isub, 2]
        #     row2 = ind_patch[isub, 3]
        #     data = data_gap[:, row1:row2 + 1, col1:col2 + 1]
        #     out_name = tempoutname + str(isub + 1) + suffix
        #     fp = path1
        #     writeimage(data, out_name, fp)

        # # **********************************************************************************************************************
        # # ***************************************begin process the gap for each block*******************************************

        # print('there are total', n_ns * n_nl, ' blocks')

        # for isub in range(0, n_nl * n_ns):

        # open each block image
        # FileName = temp_file + '\\temp_target' + str(isub + 1) + suffix
        # nl, ns, fine1_all = read_raster(FileName)
        nl = to_pred.shape[1]
        ns = to_pred.shape[2]
        fine1_all = to_pred
        # place the new image value
        # fine0 = fine1_all[:, location[isub, 2]:location[isub, 3] + 1, location[isub, 0]:location[isub, 1] + 1]
        fine0 = np.array(fine1_all)

        # FileName = temp_file + '\\temp_mask' + str(isub + 1) + suffix
        # _, _, mask = read_raster(FileName)

        # place the mark value
        mark = np.zeros(mask.shape, dtype=int)

        # classify all images by k-means
        classes = np.zeros([n_image, nl, ns]).astype(int)
        for i_img in range(0, n_image):
            imagei = fine1_all[i_img * num_band:(i_img * num_band + num_band), :, :]
            gapi = mask[i_img, :, :]
            ind_bad = np.where(gapi != 0)
            num_bad = int(int(np.size(ind_bad)) / len(ind_bad))
            if num_bad > 0:
                for iband in range(0, num_band):
                    temp = imagei[iband, :, :]
                    temp[ind_bad] = 0.0
                    imagei[iband, :, :] = temp

            # tempnamei = temp_file + '\\imagei' + suffix
            # fp = path1
            # writeimage(imagei, tempnamei, fp)

            # _, _, imagei_new = read_raster_new(tempnamei)
            imagei_new = np.moveaxis(np.array(imagei), 0, 2)
            imagei_new = np.maximum(imagei_new, 0)

            # parameter of kmeans
            CHANGE_THRESH = .02
            NUM_CLASSES = num_class + 1
            ITERATIONS = 20
            # out_bname = 'kmeans'
            # out_name = temp_file + '\\class_kmeans'
            new_shape = (imagei_new.shape[0] * imagei_new.shape[1], num_band)
            new_imagei = imagei_new[:, :, :num_band].reshape(new_shape)
            clf = MiniBatchKMeans(n_clusters=NUM_CLASSES, random_state=0, max_iter=ITERATIONS, reassignment_ratio=0.02)
            class_img = clf.fit(new_imagei).labels_
            labels = class_img.reshape(imagei_new[:, :, 0].shape)
            classes[i_img, :, :] = labels + 1

        # process each image
        for i_img in range(0, n_image):
            fine1 = fine1_all[i_img * num_band:(i_img * num_band + num_band), :, :]
            gap = mask[i_img, :, :]
            temp_dis = np.abs(np.array(DOY) - DOY[i_img])
            order = np.argsort(temp_dis)
            note_finish = 0
            i_input = 1

            while note_finish != 1 and i_input <= n_image - 1:

                fine2 = fine1_all[order[i_input] * num_band:(order[i_input] * num_band + num_band), :, :]
                mask2 = mask[order[i_input], :, :]
                class2 = classes[order[i_input], :, :]

                similar_th_band = np.zeros(num_band).astype(float)
                ind_nogrand = np.where(mask2 == 0)
                c_nogrand = int(np.size(ind_nogrand) / len(ind_nogrand))
                if c_nogrand > 0:
                    for iband in range(0, num_band):
                        # compute the threshold of similar pixel
                        fine2_ind = fine2[iband, :, :]
                        similar_th_band[iband] = np.std(fine2_ind[ind_nogrand]) * 2.0 / float(num_class)
                similar_th = np.mean(similar_th_band)

                # for i in range(location[isub, 0], location[isub, 1] + 1):
                #     for j in range(location[isub, 2], location[isub, 3] + 1):
                for i in range(to_pred.shape[2]):
                    for j in range(to_pred.shape[1]):

                        if 0 < gap[j, i] < 10:
                            if mask2[j, i] == 0:

                                classij = class2[j, i]
                                extent = np.ceil(0.5 * (np.sqrt(min_similar) - 1))  # compute the minimum window size
                                max_window = np.round(max([nl, ns]) * 0.25)  # maximum size of half seach window is 1/4 of image size

                                #  find the possible largest window
                                gap_profiles = np.hstack((gap[j, :], gap[:, i]))
                                input_profile = np.hstack((mask2[j, :], mask2[:, i]))
                                spatial_position = np.abs(np.hstack((np.arange(ns) - i, np.arange(nl) - j)))
                                ind_outcloud = np.logical_and(gap_profiles == 0, input_profile == 0)
                                num_outcloud = np.sum(ind_outcloud)
                                if num_outcloud > 0:
                                    max_window_p = np.min(spatial_position[ind_outcloud])  # find the possibel maximum half window
                                    max_window_p = np.min([max_window, max_window_p + 50])
                                else:
                                    max_window_p = max_window

                                start_w = np.ceil(0.5 * (np.sqrt(min_similar) - 1))
                                # check large cloud patches to avoid searching similar pixels for it
                                if max_window_p == max_window:
                                    a1 = int(np.max([0, i - max_window]))
                                    a2 = int(np.min([ns - 1, i + max_window]))
                                    b1 = int(np.max([0, j - max_window]))
                                    b2 = int(np.max([nl - 1, j + max_window]))

                                    sub_gap = gap[b1:b2 + 1, a1:a2 + 1]
                                    sub_class = class2[b1:b2 + 1, a1:a2 + 1]
                                    sub_mask2 = mask2[b1:b2 + 1, a1:a2 + 1]
                                    ind_class = np.logical_and(np.logical_and(sub_gap == 0, sub_mask2 == 0), sub_class == classij)
                                    c_common = np.sum(ind_class)
                                    if c_common < 2.0 * min_similar:
                                        start_w = max_window_p

                                # use binary search to find the best window size
                                end_w = max_window_p
                                ind_success = 0
                                while ind_success == 0 and start_w <= end_w:
                                    mid = np.floor(start_w + (end_w - start_w) / 2.0)
                                    a1 = int(np.max([0, i - mid]))
                                    a2 = int(np.min([ns - 1, i + mid]))
                                    b1 = int(np.max([0, j - mid]))
                                    b2 = int(np.min([nl - 1, j + mid]))
                                    sub_gap = gap[b1:b2 + 1, a1:a2 + 1]
                                    sub_mask2 = mask2[b1:b2 + 1, a1:a2 + 1]
                                    sub_class = class2[b1:b2 + 1, a1:a2 + 1]
                                    ind_common = np.logical_and(np.logical_and(sub_gap == 0, sub_mask2 == 0), sub_class == classij)
                                    c_common = np.sum(ind_common)
                                    if c_common > 3.0 * min_similar:  # search left range
                                        end_w = mid - 1
                                    else:
                                        if c_common < 2.0 * min_similar:  # search right range
                                            start_w = mid + 1
                                        else:
                                            ind_success = 1

                                extent = mid
                                sub_off = fine1[:, b1:b2 + 1, a1:a2 + 1]
                                sub_on = fine2[:, b1:b2 + 1, a1:a2 + 1]

                                col_temp = np.arange((a2 - a1 + 1)).astype(int)
                                col_num = int(b2 - b1 + 1)
                                col_wind = np.tile(col_temp, (col_num, 1))
                                row_temp = np.arange(b2 - b1 + 1).astype(int)
                                row_num = int(a2 - a1 + 1)
                                row_wind = np.repeat(row_temp, row_num).reshape(-1, row_num)
                                it = i - a1
                                jt = j - b1

                                if c_common >= 2.0 * min_similar:  # if same-class pixels more than required sample size

                                    disi = np.sqrt((it - col_wind[ind_common]) ** 2 + (jt - row_wind[ind_common]) ** 2)
                                    sub_on_common = np.zeros((num_band, c_common)).astype(float)
                                    sub_off_common = np.zeros((num_band, c_common)).astype(float)
                                    for iband in range(0, num_band):
                                        sub_on_common[iband, :] = (sub_on[iband, :, :])[ind_common] - fine2[iband, j, i]
                                        sub_off_common[iband, :] = (sub_on[iband, :, :])[ind_common] - sub_off[iband, :, :][ind_common]

                                    rmsei = np.sqrt(np.mean(sub_on_common ** 2, axis=0)) + 0.0001
                                    rmse12 = np.sqrt(np.mean(sub_off_common ** 2, axis=0)) + 0.0001

                                    order_rmse = np.argsort(rmsei)
                                    ind_similar = order_rmse[0:min_similar]   # find similar pixels

                                    similar_rmse = rmsei[ind_similar]
                                    similar_rmse12 = rmse12[ind_similar]
                                    similar_dis = disi[ind_similar]
                                    C_D = np.multiply(similar_rmse, similar_dis)
                                    weight = (1.0 / C_D) / np.sum(1.0 / C_D)
                                    T_1 = np.mean(similar_rmse)
                                    T_2 = np.mean(similar_rmse12)
                                    W_T1 = T_2 / (T_1 + T_2)
                                    W_T2 = T_1 / (T_1 + T_2)

                                    for iband in range(0, num_band):
                                        similar_on = (sub_on[iband, :, :])[ind_common][ind_similar]
                                        similar_off = (sub_off[iband, :, :])[ind_common][ind_similar]
                                        predict_1 = np.sum(np.multiply(similar_off, weight))
                                        predict_2 = fine2[iband, j, i] + np.sum(np.multiply((similar_off - similar_on), weight))
                                        if DN_min < predict_2 < DN_max:
                                            # fine0[i_img * num_band + iband, j - location[isub, 2], i - location[
                                            #     isub, 0]] = W_T1 * predict_1 + W_T2 * predict_2
                                            fine0[i_img * num_band + iband, j, i] = W_T1 * predict_1 + W_T2 * predict_2
                                        else:
                                            # fine0[i_img * num_band + iband, j - location[isub, 2], i - location[
                                            #     isub, 0]] = predict_1
                                            fine0[i_img * num_band + iband, j, i] = predict_1
                                    note = 1
                                    # mark[i_img, j - location[isub, 2], i - location[isub, 0]] = 1 + 10 * i_input
                                    mark[i_img, j, i] = 1 + 10 * i_input
                                    gap[j, i] = -2

                                else:  # if no enough same-class pixels exist
                                    # IF reach the maximum window, use all the same-class pixels
                                    if c_common >= 3 and extent >= max_window:
                                        disi = np.sqrt((it - col_wind[ind_common]) ** 2 + (jt - row_wind[ind_common]) ** 2)
                                        sub_on_common = np.zeros((num_band, c_common)).astype(float)
                                        sub_off_common = np.zeros((num_band, c_common)).astype(float)
                                        for iband in range(0, num_band):
                                            sub_on_common[iband, :] = (sub_on[iband, :, :])[ind_common] - fine2[iband, j, i]
                                            sub_off_common[iband, :] = (sub_on[iband, :, :])[ind_common] - sub_off[iband, :, :][ind_common]

                                        rmsei = np.sqrt(np.mean(sub_on_common ** 2, axis=0)) + 0.0001
                                        rmse12 = np.sqrt(np.mean(sub_off_common ** 2, axis=0)) + 0.0001

                                        similar_rmse = rmsei
                                        similar_rmse12 = rmse12
                                        similar_dis = disi
                                        C_D = np.multiply(similar_rmse, similar_dis)
                                        weight = (1.0 / C_D) / np.sum(1.0 / C_D)
                                        T_1 = np.mean(similar_rmse)
                                        T_2 = np.mean(similar_rmse12)
                                        W_T1 = T_2 / (T_1 + T_2)
                                        W_T2 = T_1 / (T_1 + T_2)
                                        for iband in range(0, num_band):
                                            similar_on = (sub_on[iband, :, :])[ind_common]
                                            similar_off = (sub_off[iband, :, :])[ind_common]
                                            predict_1 = np.sum(np.multiply(similar_off, weight))
                                            predict_2 = fine2[iband, j, i] + np.sum(np.multiply((similar_off - similar_on), weight))
                                            if DN_min < predict_2 < DN_max:
                                                # fine0[i_img * num_band + iband, j - location[isub, 2], i - location[
                                                #     isub, 0]] = W_T1 * predict_1 + W_T2 * predict_2
                                                fine0[i_img * num_band + iband, j, i] = W_T1 * predict_1 + W_T2 * predict_2
                                            else:
                                                # fine0[i_img * num_band + iband, j - location[isub, 2], i - location[
                                                #     isub, 0]] = predict_1
                                                fine0[i_img * num_band + iband, j, i] = predict_1
                                        # mark[i_img, j - location[isub, 2], i - location[isub, 0]] = 2 + 10 * i_input
                                        mark[i_img, j, i] = 2 + 10 * i_input
                                        note = 1
                                        gap[j, i] = -2
                                    else:
                                        # if no same-class pixels exist when reaching the maximum window, directly replace the pixel by input image
                                        # only use temporal interpolation if it have observations before and after
                                        first_half = np.max([0, i_img - 1])
                                        second_half = np.min([i_img + 1, n_image - 1])
                                        ind_good_before = np.where(mask[0:first_half+1, j, i] == 0)
                                        num_good_before = int(int(np.size(ind_good_before)) / len(ind_good_before))
                                        ind_good_after = np.where(mask[second_half:n_image, j, i] == 0)
                                        num_good_after = int(int(np.size(ind_good_after)) / len(ind_good_after))
                                        if num_good_before > 0 and num_good_after > 0:
                                            i_before = ind_good_before[0][num_good_before - 1]
                                            i_after = ind_good_after[0][0] + i_img + 1
                                            T1 = DOY[i_img] - DOY[int(i_before)]
                                            T2 = DOY[int(i_after)] - DOY[i_img]
                                            for iband in range(0, num_band):
                                                before_band = fine1_all[i_before * num_band + iband, j, i]
                                                after_band = fine1_all[i_after * num_band + iband, j, i]
                                                # fine0[i_img * num_band + iband, j - location[isub, 2], i - location[isub, 0]] = (T2 * before_band + T1 * after_band) / (T2 + T1)
                                                fine0[i_img * num_band + iband, j, i] = (T2 * before_band + T1 * after_band) / (T2 + T1)
                                        else:
                                            # fine0[i_img * num_band:((i_img + 1) * num_band), j - location[isub, 2], i - location[isub, 0]] = fine2[:, j, i]
                                            fine0[i_img * num_band:((i_img + 1) * num_band), j, i] = fine2[:, j, i]

                                        # mark[i_img, j - location[isub, 2], i - location[isub, 0]] = 3 + 10 * i_input
                                        mark[i_img, j, i] = 3 + 10 * i_input
                                        note = 1
                                        gap[j, i] = -2

                # ind_unfilled = np.where(gap[location[isub, 2]:location[isub, 3] + 1, location[isub, 0]:location[isub, 1] + 1] > 0)
                ind_unfilled = np.where(gap > 0)
                c_unfilled = int(int(np.size(ind_unfilled)) / len(ind_unfilled))
                if c_unfilled == 0:
                    note_finish = 1
                else:
                    i_input = i_input + 1

            # print('finish', i_img + 1, "image in the time-series")
            progr_pred[split] = float(i_img + 1) / n_image * 100
            self.cons.progress(progr_pred, split)
            # from probgf.visualization import visualize_rgb
            # img = np.moveaxis(fine0[i_img * num_band:((i_img + 1) * num_band)], 0, 2)
            # visualize_rgb(img, 'test_{}_{}.png'.format(split, i_img), np.array([0, 1, 2]), np.array([0, 0, 0]), np.array([3000, 3000, 3000]))

        # reshape to probgf format
        fine0 = np.moveaxis(fine0, 0, 2).reshape(orig_shape)
        super().run_prediction(to_pred, obs, split, progr_pred) # for final console output

        return fine0

        # print("finished", isub + 1, "block")
        # tempoutname1 = temp_file + '\\temp_filled' + str(isub + 1) + suffix
        # tempoutname2 = temp_file + '\\temp_mark' + str(isub + 1) + suffix
        # fp = path1
        # writeimage(fine0, tempoutname1, fp)
        # writeimage(mark, tempoutname2, fp)

        # # *******************************************************************************************
        # # *******************************mosaic all the filled patch*********************************

        # datalist = []
        # minx_list = []
        # maxX_list = []
        # minY_list = []
        # maxY_list = []

        # for isub in range(0, n_ns * n_nl):
        #     out_name = temp_file + '\\temp_filled' + str(isub + 1) + suffix
        #     datalist.append(out_name)
        #     col1 = ind_patch1[isub, 0]
        #     col2 = ind_patch1[isub, 1]
        #     row1 = ind_patch1[isub, 2]
        #     row2 = ind_patch1[isub, 3]

        #     minx_list.append(col1)
        #     maxX_list.append(col2)
        #     minY_list.append(row1)
        #     maxY_list.append(row2)

        # minX = min(minx_list)
        # maxX = max(maxX_list)
        # minY = min(minY_list)
        # maxY = max(maxY_list)

        # xOffset_list = []
        # yOffset_list = []
        # i = 0
        # for data in datalist:
        #     ds = gdal.Open(data)
        #     transform = ds.GetGeoTransform()
        #     pixelWidth = transform[1]
        #     pixelHeight = transform[5]
        #     xOffset = int((minx_list[i] - minX) / pixelWidth)
        #     yOffset = int((minY_list[i] - minY) / pixelHeight)
        #     xOffset_list.append(xOffset)
        #     yOffset_list.append(yOffset)
        #     i += 1

        # in_ds = gdal.Open(path1)
        # path = os.path.splitext(path1)[0] + "_fill_NSPI_all" + suffix
        # if suffix == '.tif':
        #     driver = gdal.GetDriverByName("GTiff")
        # elif suffix == "":
        #     driver = gdal.GetDriverByName("ENVI")
        # dataset = driver.Create(path, orig_ns, orig_nl, nb, gdal.GDT_Float32)
        # i = 0
        # for data in datalist:
        #     nl, ns, datavalue = read_raster(data)
        #     for j in range(0, nb):
        #         dataset.GetRasterBand(j + 1).WriteArray(datavalue[j, :, :], xOffset_list[i], yOffset_list[i])
        #     i += 1

        # geoTransform = in_ds.GetGeoTransform()
        # dataset.SetGeoTransform(geoTransform)
        # proj = in_ds.GetProjection()
        # dataset.SetProjection(proj)

        # # mosiac all the mask patch
        # masklist = []
        # for isub in range(0, n_ns * n_nl):
        #     out_name = temp_file + '\\temp_mark' + str(isub + 1) + suffix
        #     masklist.append(out_name)

        # n_mask = int(n_image)
        # in_ds = gdal.Open(path1)
        # path = os.path.splitext(path1)[0] + "_fill_quality_flag_all" + suffix
        # if suffix == '.tif':
        #     driver = gdal.GetDriverByName("GTiff")
        # elif suffix == "":
        #     driver = gdal.GetDriverByName("ENVI")
        # maskset = driver.Create(path, orig_ns, orig_nl, n_mask, gdal.GDT_Float32)
        # i = 0
        # for mask in masklist:
        #     nl, ns, maskvalue = read_raster(mask)
        #     for j in range(0, n_mask):
        #         maskset.GetRasterBand(j + 1).WriteArray(maskvalue[j, :, :], xOffset_list[i], yOffset_list[i])
        #     i += 1

        # geoTransform = in_ds.GetGeoTransform()
        # maskset.SetGeoTransform(geoTransform)
        # proj = in_ds.GetProjection()
        # maskset.SetProjection(proj)

        # endtime = datetime.datetime.now()
        # print("Finish mosaic image patches!")
        # print('Time used: ', (endtime - starttime).seconds, 'seconds')

max_stars_repo_name: raphischer/probgf

max_issues_repo_name: raphischer/probgf

max_forks_repo_name: raphischer/probgf


text:
import numpy as np

# takes in patterns that the neural network will learn
def hopfieldNetwork(patterns):
	# I make a weight array with length of a pattern * length of a pattern
	weights = np.zeros((len(patterns[0]), len(patterns[0])))

	# take the mean of all the patterns
	mean = np.mean(patterns)

	# I deduct mean from patterns array
	patterns -= mean

	for pattern in patterns:
		weights += np.outer(pattern, pattern) / len(patterns)

	return weights

def bipolarToBinary(weights):
	weights[weights < 0] = 0
	weights[weights > 0] = 1
	weights[weights == 0] = 0

	return weights

def binaryToBipolar(weights):
	weights[weights < 0] = -1
	weights[weights > 0] = 1
	weights[weights == 0] = -1

	return weights

def queryPattern(weights, testPatterns, iteration = 100):
	# we need to fill the diagonals with 0 first
	np.fill_diagonal(weights, 0)

	testPatterns = binaryToBipolar(testPatterns)

	results = []

	for pattern in testPatterns:

		newPattern = weights.dot(pattern) - iteration

		newPattern = bipolarToBinary(newPattern)

		results.append(newPattern)

	return np.array(results)

ext: py

max_stars_repo_name: tanvirtin/Hopfield-Network

max_issues_repo_name: tanvirtin/Hopfield-Network

max_forks_repo_name: tanvirtin/Hopfield-Network

reason: import numpy


text:
import argparse
import json
import os

import numpy as np
import torch
import SimpleITK as sitk
from torchvision.transforms import InterpolationMode
from scipy import ndimage
import torch.nn.functional as F
import cc3d
from scipy.ndimage.morphology import binary_opening, binary_dilation, binary_erosion

from config import ExperimentConfigs
from datasets.ct_dataset import get_transforms, LITS2017_VALSETS
from datasets.visualize_data import write_volume_slices
from models import get_model
from torchvision.transforms import Resize
from metrics import TverskyScore, compute_Recal

def get_model_from_dir(model_dir, ckpt_name):
    config = json.load(open(f'{model_dir}/exp_configs.json'))
    config = ExperimentConfigs(**config)
    ckpt = torch.load(f'{model_dir}/{ckpt_name}.pth')

    # get model
    model = get_model(config.get_model_config())
    model.load_state_dict(ckpt['model'])
    model.to(config.device)

    return model, config

def clean_liver_prediction(pred_volume):
    """
    Keep only largest blob in binary mask
    """
    # Convert to  binary-mask
    binary_mask = pred_volume.cpu().detach().bool().numpy()

    # erode to separate loosley connected blobs
    binary_mask = binary_erosion(binary_mask, iterations=5)

    if not binary_mask.any():
        return torch.from_numpy(binary_mask)
    # Find 3D blobs
    cc = cc3d.connected_components(binary_mask)
    clean_pred_volume = max([image for (label, image) in cc3d.each(cc, binary=True, in_place=False)],
                            key=lambda x: x.sum())

    # Restore blob to original
    clean_pred_volume = binary_dilation(clean_pred_volume, iterations=5)

    clean_pred_volume = torch.from_numpy(clean_pred_volume)

    return clean_pred_volume

def read_case(ct_dir, gt_dir, fname):
    ct_nii = sitk.ReadImage(os.path.join(ct_dir, fname))
    ct_volume = sitk.GetArrayFromImage(ct_nii)  # shape (S, 512, 512)
    gt_path = os.path.join(gt_dir, fname.replace('volume', 'segmentation'))

    gt_volume = sitk.GetArrayFromImage(sitk.ReadImage(gt_path, sitk.sitkUInt8)) if os.path.exists(gt_path) else None

    return ct_volume, gt_volume, ct_nii.GetSpacing()

class TwoStepsSegmentor:
    def __init__(self, liver_segmentation_model_dir, tumor_segmentation_model_dir):
        self.liver_segmentation_model, self.liver_segmentation_cfg = get_model_from_dir(liver_segmentation_model_dir, 'best')
        _, self.liver_segmentation_transforms = get_transforms(self.liver_segmentation_cfg.get_data_config())

        # Load Tumor model and cnfigs
        self.tumor_model, self.tumor_cfg = get_model_from_dir(tumor_segmentation_model_dir, 'best')
        _, self.tumor_transforms = get_transforms(self.tumor_cfg.get_data_config())

    def predict(self, cropped_ct):
        # Predict liver
        input_ct = self.liver_segmentation_transforms((cropped_ct, np.ones_like(cropped_ct)))[0].unsqueeze(0).float().cuda()
        predicted_liver_mask = self.liver_segmentation_model.predict_volume(input_ct)[0].argmax(0).cpu()
        predicted_liver_mask = clean_liver_prediction(predicted_liver_mask).long()

        # Predict tumors
        input_ct = self.tumor_transforms((cropped_ct, np.ones_like(cropped_ct)))[0].unsqueeze(0).float().cuda()  # shape (S, 128, 128)

        predicted_tumor_mask = self.tumor_model.predict_volume(input_ct)[0].argmax(0).cpu()  # shape (S, 128, 128)

        multiclass_mask = predicted_liver_mask
        multiclass_mask[torch.logical_and(multiclass_mask == 1, predicted_tumor_mask == 1)] = 2

        return multiclass_mask

class OneStepsSegmentor:
    def __init__(self, multiclass_model_dir):
        self.multiclass_model, self.multiclass_cfg = get_model_from_dir(multiclass_model_dir, 'best')
        _, self.multiclass_transforms = get_transforms(self.multiclass_cfg.get_data_config())

    def predict(self, cropped_ct):
        input_ct = self.multiclass_transforms((cropped_ct, np.ones_like(cropped_ct)))[0].unsqueeze(0).float().cuda()  # shape (S, 128, 128)
        multiclass_mask = self.multiclass_model.predict_volume(input_ct)[0].argmax(0).cpu()  # shape (S, 128, 128)

        return multiclass_mask

def inference(ct_path, gt_path, liver_localization_model_dir, multiclass_segmentation_model_dir,  normalized_mms=None, liver_crop_padding=(3, 20,20), dump_debug_images=False):
    with torch.no_grad():
        outputs_dir = os.path.join(os.path.dirname(ct_path), 'end2end_prediction')
        os.makedirs(outputs_dir, exist_ok=True)

        liver_localization_model, liver_localization_cfg = get_model_from_dir(liver_localization_model_dir, 'best')
        _, liver_localization_transforms = get_transforms(liver_localization_cfg.get_data_config())

        segmentor = OneStepsSegmentor(multiclass_segmentation_model_dir)

        liver_dice_scores = []
        tumor_dice_scores = []
        tumor_recalls = []
        for fname in os.listdir(ct_path):
            print(f"Case-{fname}")
            ct_volume, gt_volume, spacing = read_case(ct_path, gt_path, fname)

            # localize liver
            liver_input = liver_localization_transforms((ct_volume.copy(), np.ones_like(ct_volume).astype(np.uint8)))[0].unsqueeze(0).float().cuda()
            predicted_liver_mask = liver_localization_model.predict_volume(liver_input)[0].argmax(0).cpu()                                        # shape (S, 128, 128)

            # Clean liver and restore to origial size
            predicted_liver_mask = clean_liver_prediction(predicted_liver_mask)                                                      # shape (S, 128, 128)
            predicted_liver_mask = Resize(ct_volume.shape[-2:], interpolation=InterpolationMode.NEAREST)(predicted_liver_mask)       # shape (S, 512, 512)

            # Crop around liver for tumor segmentaiton                                                                           # shape (S, h, w)
            nwhere = np.where(predicted_liver_mask)
            liver_crop = tuple([slice(max(0, x.min() - liver_crop_padding[i]), x.max() + liver_crop_padding[i]) for i, x in enumerate(nwhere)])
            cropped_ct = ct_volume[liver_crop]                                                                                # shape (S, h, w)

            if normalized_mms is not None:
                cropped_ct = ndimage.zoom(cropped_ct, (spacing[-1] / normalized_mms, 1, 1), order=3)

            multiclass_mask = segmentor.predict(cropped_ct)

            # Restore input resolution
            multiclass_mask = F.interpolate(multiclass_mask.float().unsqueeze(0).unsqueeze(0), size=ct_volume[liver_crop].shape[-3:], mode='nearest')[0,0].long()

            # Create final prediction mask
            final_mask = torch.zeros(ct_volume.shape).long()
            final_mask[liver_crop] = multiclass_mask

            # write segmentation map
            new_seg = sitk.GetImageFromArray(final_mask, sitk.sitkInt8)
            new_seg.SetSpacing(spacing)
            sitk.WriteImage(new_seg, os.path.join(outputs_dir, fname.replace('volume', 'segmentation')))

            # # Compute dice
            if gt_volume is not None:
                gt_volume = torch.from_numpy(gt_volume)
                liver_score = TverskyScore(0.5,0.5)((final_mask != 0).long().unsqueeze(0), (gt_volume != 0).long().unsqueeze(0), torch.ones_like(gt_volume).unsqueeze(0))
                tumor_score = TverskyScore(0.5,0.5)((final_mask == 2).long().unsqueeze(0), (gt_volume == 2).long().unsqueeze(0), torch.ones_like(gt_volume).unsqueeze(0))
                tumor_recall = compute_Recal((final_mask == 2).long().unsqueeze(0), (gt_volume == 2).long().unsqueeze(0), torch.ones_like(gt_volume).unsqueeze(0))
                liver_dice_scores.append(liver_score)
                tumor_dice_scores.append(tumor_score)
                tumor_recalls.append(tumor_recall)

                # dump debug images
                if dump_debug_images:
                    ct_volume = torch.from_numpy(np.clip(ct_volume, -100, 400).astype(float))
                    write_volume_slices(ct_volume, [final_mask, gt_volume], os.path.join(outputs_dir, f"{os.path.splitext(fname)[0]}_{liver_score.item():.2f}_{tumor_score.item():.2f}"))

                print(f"AVG Dice per case: Liver: {np.mean(liver_dice_scores)}, Tumor: {np.mean(tumor_dice_scores)}, Tumor-Recall: {np.mean(tumor_recalls)}")

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Preprocess Lits2017 dataset')
    parser.add_argument('ct_dir')
    parser.add_argument('--gt_dir', default="", help='If GT is not specified no Dice score is computed')
    parser.add_argument('--localization_model_dir', default='trained_models/liver_localization/VGGUNet_Aug_Loss(0.0Dice+0.0WCE+1.0CE)_V-A')
    parser.add_argument('--segmentation_model_dir', default='trained_models/multiclass_segmentation/VGGUNet2_5D_Aug_FNE-0.5_Loss(0.0Dice+0.0WCE+1.0CE)_V-A')
    parser.add_argument('--debug', action='store_true')
    args = parser.parse_args()

    inference(args.ct_dir, args.gt_dir, args.localization_model_dir, args.segmentation_model_dir, dump_debug_images=args.debug)

reason: import numpy,from scipy


text:
from pathlib import Path

import numpy as np

from napari_lazy_openslide import napari_get_reader

fixture_dir = Path.cwd() / ".." / ".." / "fixtures"

# tmp_path is a pytest fixture
def test_reader(tmp_path):
    """An example of how you might test your plugin."""
    # TODO: create sample fixture
    # np.testing.assert_allclose(original_data, layer_data_tuple[0])

def test_get_reader_pass():
    assert napari_get_reader("fake.file") is None
    assert napari_get_reader(["d0.tif", "d1.tif"]) is None
    assert napari_get_reader(str(fixture_dir / "dataset.ome.tif")) is None

max_stars_repo_name: sdvillal/napari-lazy-openslide

max_issues_repo_name: sdvillal/napari-lazy-openslide

max_forks_repo_name: sdvillal/napari-lazy-openslide


text:
import glob
import os
import librosa
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation
from keras.layers import Dropout
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

def extract_feature(file_name):
  X, sample_rate = librosa.load(file_name)
  stft = np.abs(librosa.stft(X))
  mfccs = np.mean(librosa.feature.mfcc(y=X, sr=sample_rate, n_mfcc=40).T,axis=0)
  chroma = np.mean(librosa.feature.chroma_stft(S=stft, sr=sample_rate).T,axis=0)
  mel = np.mean(librosa.feature.melspectrogram(X, sr=sample_rate).T,axis=0)
  contrast = np.mean(librosa.feature.spectral_contrast(S=stft, sr=sample_rate).T,axis=0)
  tonnetz = np.mean(librosa.feature.tonnetz(y=librosa.effects.harmonic(X),
  sr=sample_rate).T,axis=0)
  return mfccs,chroma,mel,contrast,tonnetz;

def parse_audio_files(parent_dir,sub_dirs,file_ext="*.wav"):
    features, labels = np.empty((0,193)), np.empty(0)
    for label, sub_dir in enumerate(sub_dirs):
        for fn in glob.glob(os.path.join(parent_dir, sub_dir, file_ext)):
            try:
              mfccs, chroma, mel, contrast,tonnetz = extract_feature(fn)
            except Exception as e:
              print ("Error encountered while parsing file: ", fn)
              continue
            ext_features = np.hstack([mfccs,chroma,mel,contrast,tonnetz])
            features = np.vstack([features,ext_features])
            labels = np.append(labels, fn.split('/')[2].split('-')[2])
    return np.array(features), np.array(labels, dtype = np.int)

def one_hot_encode(labels):
  n_labels = len(labels)
  n_unique_labels = len(np.unique(labels))
  one_hot_encode = np.zeros((n_labels,n_unique_labels+1))
  one_hot_encode[np.arange(n_labels), labels] = 1
  one_hot_encode=np.delete(one_hot_encode, 0, axis=1)
  return one_hot_encode;

main_dir = 'Audio_Speech_Actors_01-24'
sub_dir=os.listdir(main_dir)
print ("\ncollecting features and labels...")
print("\nthis will take some time...")
features, labels = parse_audio_files(main_dir,sub_dir)
print("done")
#one hot encoding labels
labels = one_hot_encode(labels)
np.save('features',labels)
np.save('labels',features)

#prep train and test set
X=np.load('features.npy')
y=np.load('labels.npy')

train_x, test_x, train_y, test_y = train_test_split(X, y, test_size=0.33, random_state=42)

#build DNN
n_dim = train_x.shape[1]
n_classes = train_y.shape[1]
n_hidden_units_1 = n_dim
n_hidden_units_2 = 400 # approx n_dim * 2
n_hidden_units_3 = 200 # half of layer 2
n_hidden_units_4 = 100

def create_model(activation_function='relu', init_type='normal', optimiser='adam', dropout_rate=0.2):
  model = Sequential()
  # layer 1
  model.add(Dense(n_hidden_units_1, input_dim=n_hidden_units_1, kernel_initializer="normal", activation="relu"))
  # layer 2
  model.add(Dense(n_hidden_units_2, kernel_initializer="normal", activation="relu"))
  model.add(Dropout(dropout_rate))
  # layer 3
  model.add(Dense(n_hidden_units_3, kernel_initializer="normal", activation="relu"))
  model.add(Dropout(dropout_rate))
  #layer4
  model.add(Dense(n_hidden_units_4, kernel_initializer="normal", activation="relu"))
  model.add(Dropout(dropout_rate))
  # output layer
  model.add(Dense(n_classes, kernel_initializer="normal", activation="softmax"))
  #model compilation
  model.compile(loss='categorical_crossentropy', optimizer=optimiser, metrics=['accuracy'])
  return model;

#create the model
model = create_model()
model.summary()
#train the model
history = model.fit(train_x,train_y,epochs=200,batch_size=4)

#predicting from the model
predict=model.predict(test_x,batch_size=4)
np.save('prediction',predict)

emotions=['neutral', 'calm', 'happy', 'sad', 'angry', 'fearful', 'disgust', 'surprised']
#predicted emotions from the test set
y_pred = np.argmax(predict, 1)
np.save('y_pred',y_pred)
"""
predicted_emo=[]
for i in range(y_pred.shape[0]):
  emo=emotions[y_pred[i]]
  predicted_emo.append(emo)

actual_emo=[]
y_true=np.argmax(test_y, 1)
np.save('test_y',y_true)
l = np.load('test_y.npy')
print(l)

for i in range(0,test_y.shape[0]):
  emo=emotions[y_true[i]]
  actual_emo.append(emo)

np.save('emotion_array',actual_emo)
emotion_array = np.load('emotion_array.npy')
print(emotion_array)

#generate the confusion matrix
cm =confusion_matrix(actual_emo, actual_emo)
index = ['angry', 'calm', 'disgust', 'fearful', 'happy', 'neutral', 'sad', 'surprised']
columns = ['angry', 'calm', 'disgust', 'fearful', 'happy', 'neutral', 'sad', 'surprised']
cm_df = pd.DataFrame(cm,index,columns)
plt.figure(figsize=(10,6))
sns.heatmap(cm_df, annot=True)
"""

max_stars_repo_name: JumanaFM/emoviz

max_issues_repo_name: JumanaFM/emoviz

max_forks_repo_name: JumanaFM/emoviz


text:
import logging
from typing import List

import numpy as np
import torch
import torch.nn.functional as F
from torch.nn.functional import gumbel_softmax
from torch.distributions.categorical import Categorical
import wandb
from torch import nn, optim

import RL
from RL.agents.exp_buff_agent import ExperienceBuffer
from RL.utils.standard_models import FFModel
from RL.utils.util_fns import toNpFloat32

logger = logging.getLogger(__name__)
ldebug = logger.isEnabledFor(logging.DEBUG)
device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')

class Actor(nn.Module):
    def __init__(self, input_shape, convs, hidden_layers, n_actions, ln=False):
        super().__init__()
        self.logits = FFModel(input_shape, convs, list(
            hidden_layers) + [n_actions], ln=ln)
        self.logits.linear_layers[-1].bias.data.fill_(0)

    def forward(self, state):
        logits = Categorical(logits=self.logits(state)).logits
        return logits

    def sample_one_hot(self, state, return_logpi=False):
        act_one_hot = gumbel_softmax(self.logits(state), tau=1)
        if return_logpi:
            dist = Categorical(logits=self.logits(state))
            act = torch.argmax(act_one_hot, dim=-1)  # [n]
            logpi_act = dist.logits[torch.arange(len(state)), act]  # [batch]
            return act_one_hot, logpi_act
        return act_one_hot

    def sample(self, state):
        dist = Categorical(logits=self.logits(state))
        action = dist.sample()
        return action

    def entropy(self, state):
        dist = Categorical(logits=self.logits(state))
        return dist.entropy()

# class Critic(nn.Module):
#     def __init__(self, input_shape, convs, hidden_layers, n_actions, ln=False):
#         super().__init__()
#         self.convs = FFModel(input_shape, convs, [],
#                              flatten=True, apply_act_output=True, ln=ln)
#         self.q = FFModel(
#             (self.convs.output_shape[0] + n_actions,), [], linears=list(hidden_layers) + [1], ln=ln)
#         self.q.linear_layers[-1].bias.data.fill_(0)

#     def forward(self, state, action_one_hot):
#         convs = self.convs(state)
#         s_a_concat = torch.cat((convs, action_one_hot), dim=1)
#         q = self.q(s_a_concat)[:, 0]
#         return q

class CriticDQN(FFModel):
    '''Holds q value for every action'''

    def __init__(self, input_shape, convs, hidden_layers, n_actions, ln=False):
        super().__init__(input_shape, convs, list(
            hidden_layers) + [n_actions], ln=ln)
        self.linear_layers[-1].bias.data.fill_(0)

class SACDiscreteAgent(RL.Agent):
    def __init__(self, name: str, algo: RL.Algorithm, convs: List, hidden_layers: List, train_freq: int, sgd_steps: int, mb_size: int, dqn_mse_loss: bool, gamma: float, nsteps: int, td_clip: float, grad_clip: float, lr: float, a_lr: float, eval_mode: bool, no_train_for_steps: int, exp_buffer: ExperienceBuffer, alpha=0.2, fix_alpha=False):
        super().__init__(name, algo, supports_multiple_envs=False)
        self.convs = convs
        self.hidden_layers = hidden_layers
        self.train_freq = train_freq
        self.sgd_steps = sgd_steps
        self.mb_size = mb_size
        self.dqn_mse_loss = dqn_mse_loss
        self.gamma = gamma
        self.nsteps = nsteps
        self.td_clip = td_clip
        self.grad_clip = grad_clip
        self.lr = lr
        self.a_lr = a_lr
        self.no_train_for_steps = no_train_for_steps
        self.exp_buffer = exp_buffer
        assert alpha >= 0, "alpha must be >= 0"
        self.alpha = alpha
        self.fix_alpha = fix_alpha
        self.min_alpha = alpha
        self.derired_ent = np.log(self.env.action_space.n / 2)
        self.eval_mode = eval_mode
        self._loss, self._a_loss, self._mb_v, self._mb_ent = 0, 0, 0, 0

        logger.info(f'Creating main actor critic networks on device {device}')
        self.critic1 = CriticDQN(self.env.observation_space.shape,
                                 convs, hidden_layers, self.env.action_space.n)
        self.critic1.to(device)
        self.critic2 = CriticDQN(self.env.observation_space.shape,
                                 convs, hidden_layers, self.env.action_space.n)
        self.critic2.to(device)
        logger.info(str(self.critic1))
        self.actor = Actor(self.env.observation_space.shape,
                           convs, hidden_layers, self.env.action_space.n)
        self.actor.to(device)
        logger.info(str(self.actor))
        if ldebug:
            wandb.watch([self.critic1], log='all', log_freq=1000 // train_freq)
            wandb.watch([self.actor], log='all', log_freq=1000 // train_freq)
        if not eval_mode:
            logger.info(f'Creating target critic networks on device {device}')
            self.target_critic1 = CriticDQN(
                self.env.observation_space.shape, convs, hidden_layers, self.env.action_space.n)
            self.target_critic1.to(device)
            self.target_critic2 = CriticDQN(
                self.env.observation_space.shape, convs, hidden_layers, self.env.action_space.n)
            self.target_critic2.to(device)
            logger.info(f'Creating critic optimizer (lr={lr})')
            self.optim_critic = optim.Adam(
                list(self.critic1.parameters()) + list(self.critic2.parameters()), lr=self.lr)
            logger.info(f'Creating actor optimizer (lr={a_lr})')
            self.optim_actor = optim.Adam(
                self.actor.parameters(), lr=self.a_lr)

    def act(self):
        if self.manager.episode_type > 0:
            with torch.no_grad():
                obs = torch.from_numpy(toNpFloat32(
                    self.manager.obs, True)).to(device)
                action = self.actor.sample(obs).item()
            return action, {}
        else:
            ldebug and logger.debug('Behavior Policy')
            if self.manager.num_steps < self.no_train_for_steps:
                ldebug and logger.debug('Random action')
                return self.env.action_space.sample(), {}
            else:
                ldebug and logger.debug('Stochastic action')
                with torch.no_grad():
                    obs = torch.from_numpy(toNpFloat32(
                        self.manager.obs, True)).to(device)
                    action = self.actor.sample(obs).item()
                return action, {}

    def sgd_update(self):
        ldebug and logger.debug('Training')
        with torch.no_grad():
            states, actions, rewards, dones, info, next_states = self.exp_buffer.random_experiences_unzipped(
                self.mb_size)
            states = torch.from_numpy(toNpFloat32(states)).to(device)
            rewards = torch.from_numpy(toNpFloat32(rewards)).to(device)
            dones = torch.from_numpy(toNpFloat32(dones)).to(device)
            next_states = torch.from_numpy(
                toNpFloat32(next_states)).to(device)
            # next_actions, next_logpis = self.actor.sample_one_hot(
            #     next_states, return_logpi=True)
            # next_target_q = torch.min(self.target_critic1(
            #     next_states, next_actions), self.target_critic2(next_states, next_actions))
            # next_target_v = (next_target_q - self.alpha * next_logpis)

            # desired_q = rewards + (1 - dones) * \
            #     (self.gamma ** self.nsteps) * next_target_v
            next_logpis = self.actor(next_states)
            next_pis = torch.exp(next_logpis)
            next_target_q = torch.min(self.target_critic1(
                next_states), self.target_critic2(next_states))
            next_target_v = torch.sum(
                next_pis * (next_target_q - self.alpha * next_logpis), dim=-1)
            desired_q = (rewards + (1 - dones) * (self.gamma **
                                                  self.nsteps) * next_target_v)

        '''-------------------begin optimization----------------'''
        all_states_idx = torch.arange(self.mb_size).to(device)

        # optimize critic
        # self.optim_critic.zero_grad()
        # actions_one_hot = np.zeros((self.mb_size, self.env.action_space.n))
        # actions_one_hot[all_states_idx, actions] = 1
        # actions_one_hot = torch.from_numpy(
        #     toNpFloat32(actions_one_hot)).to(device)
        # critic_loss = 0.5 * (F.mse_loss(self.critic1(states, actions_one_hot),
        #                                 desired_q) + F.mse_loss(self.critic2(states, actions_one_hot), desired_q))
        # critic_loss.backward()
        # if self.grad_clip:
        #     nn.utils.clip_grad_norm_(list(self.critic1.parameters(
        #     )) + list(self.critic2.parameters()), self.grad_clip)
        # self.optim_critic.step()

        q1, q2 = self.critic1(states), self.critic2(states)
        q1_desired, q2_desired = q1.detach().clone(), q2.detach().clone()
        td1 = desired_q - q1_desired[all_states_idx, actions]
        td2 = desired_q - q2_desired[all_states_idx, actions]
        if self.td_clip is not None:
            ldebug and logger.debug('Doing TD Clipping')
            td1 = torch.clamp(td1, -self.td_clip, self.td_clip)
            td2 = torch.clamp(td2, -self.td_clip, self.td_clip)
        q1_desired[all_states_idx,
                   actions] = q1_desired[all_states_idx, actions] + td1
        q2_desired[all_states_idx,
                   actions] = q2_desired[all_states_idx, actions] + td2

        self.optim_critic.zero_grad()
        if self.dqn_mse_loss:
            critic_loss = 0.5 * (F.mse_loss(q1, q1_desired) +  # noqa
                                F.mse_loss(q2, q2_desired))
        else:
            critic_loss = 0.5 * (F.smooth_l1_loss(q1, q1_desired) +  # noqa
                                F.smooth_l1_loss(q2, q2_desired))
        critic_loss.backward()
        if self.grad_clip:
            nn.utils.clip_grad_norm_(list(self.critic1.parameters(
            )) + list(self.critic2.parameters()), self.grad_clip)
        self.optim_critic.step()

        # optimize actor
        # self.optim_actor.zero_grad()
        # actions_one_hot = self.actor.sample_one_hot(states)
        # q = torch.min(self.critic1(states, actions_one_hot),
        #               self.critic2(states, actions_one_hot))
        # entropy = torch.mean(self.actor.entropy(states))
        # mean_v = torch.mean(q) + self.alpha * entropy
        # actor_loss = -mean_v
        # actor_loss.backward()
        # if self.grad_clip:
        #     nn.utils.clip_grad_norm_(self.actor.parameters(), self.grad_clip)
        # self.optim_actor.step()

        self.optim_actor.zero_grad()
        logpis = self.actor(states)
        pis = torch.exp(logpis)
        q = torch.min(q1_desired, q2_desired)
        entropy = torch.mean(-torch.sum(pis * logpis, dim=-1))
        mean_v = torch.mean(torch.sum(pis * q, dim=-1)) + self.alpha * entropy
        actor_loss = -mean_v
        actor_loss.backward()
        if self.grad_clip:
            nn.utils.clip_grad_norm_(self.actor.parameters(), self.grad_clip)
        self.optim_actor.step()

        entropy = entropy.item()
        if not self.fix_alpha:
            logalpha = np.log(self.alpha)
            logalpha -= 0.001 * (entropy - self.derired_ent)
            self.alpha = max(np.exp(logalpha), self.min_alpha)

        '''====================optimization done================'''

        self._loss = critic_loss.item()
        self._a_loss = actor_loss.item()
        self._mb_v = mean_v.item()
        self._mb_ent = entropy

    def post_act(self):
        if self.manager.episode_type < 2 and self.manager.num_steps > self.no_train_for_steps and (self.manager.num_steps - 1) % self.train_freq == 0:
            for sgd_step in range(self.sgd_steps):
                self.sgd_update()

        if (self.manager.num_steps - 1) % 1000 == 0:
            wandb.log({
                'SAC/Loss': self._loss,
                'SAC/A_Loss': self._a_loss,
                'SAC/Value': self._mb_v,
                'SAC/Entropy': self._mb_ent,
                'SAC/Alpha': self.alpha
            }, step=self.manager.num_steps - 1)

    def post_episode(self):
        wandb.log({
            'SAC/Loss': self._loss,
            'SAC/A_Loss': self._a_loss,
            'SAC/Value': self._mb_v,
            'SAC/Entropy': self._mb_ent,
            'SAC/Alpha': self.alpha
        }, step=self.manager.num_steps - 1)


text:
import numpy as np
import gym
from copy import deepcopy as copy
import tensorflow as tf
from abc import ABC, abstractmethod
import os

def single_elem_support(func):
    """aop func"""
    type_list = (type([]), type(()), type(np.array(1)))

    def wrapper(*args, **kwargs):
        """wrapper func"""
        res = func(*args, **kwargs)
        if type(res) in type_list and len(res) == 1:
            return res[0]
        elif type(res[0]) in type_list and len(res[0]) == 1:
            return [x[0] for x in res]
        else:
            return res

    return wrapper

class RecState(ABC):
    def __init__(self, config, records):
        self.config = config
        self.records = records
        self._init_state = self.records_to_state(records)
        self._state = copy(self._init_state)

    @staticmethod
    def records_to_state(records):
        pass

    @property
    def state(self):
        return self._state

    @property
    @abstractmethod
    def user(self):
        pass

    @property
    @abstractmethod
    def info(self):
        pass

    @abstractmethod
    def act(self, actions):
        pass

    @abstractmethod
    def to_string(self):
        pass

class RecDataBase(object):
    '''
    file-based implementation of a RecommnedEnv's data source.

    Pulls data from file, preps for use by RecommnedEnv and then
    acts as data provider for each new episode.
    '''

    def __init__(self, config, state_cls):
        self.config = config
        self.sample_list = []
        self.state_cls = state_cls
        self.is_eval = config.get('is_eval', False)
        self.cache_size = config.get('cache_size', 2048)
        # sample file cache
        self.fp = open(config['sample_file'], 'r')
        # self.fp.readline()

    @staticmethod
    def seed(seed):
        np.random.seed(seed)

    def sample_cache(self, f, num):
        for i in range(num):
            tmp = f.readline().rstrip()
            if len(tmp) < 1:
                f.seek(0, 0)
                f.readline()
                self.sample_list.append(f.readline().rstrip())
            else:
                self.sample_list.append(tmp)

    def sample(self, batch_size):
        if self.is_eval:
            assert self.cache_size == batch_size
            assert len(self.sample_list) == batch_size
            records = self.sample_list[:batch_size]
        else:
            records = np.random.choice(self.sample_list, batch_size)
        samples = self.state_cls(self.config, records)
        return samples

    def reset(self, reset_file=False):
        # self.state_list = []
        self.sample_list = []
        # self.rawstate_cache(self.fs, 10000)
        if reset_file:
            self.fp.seek(0, 0)
        self.sample_cache(self.fp, self.cache_size)

class RecSimBase(ABC):
    """ Implemention of core recommendation simulator"""

    def __init__(self, config, state_cls):
        self.config = config
        self.max_steps = config['max_steps']
        self.batch_size = config['batch_size']
        model_file = config['model_file']
        self.graph = tf.Graph()
        with self.graph.as_default():
            self.model = self.get_model(config)
            if self.config.get('gpu', False):
                os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
                self.sess = tf.Session(graph=self.graph,
                                       config=tf.ConfigProto(device_count={"CPU": 4}))
            else:
                os.environ['TF_FORCE_GPU_ALLOW_GROWTH'] = 'true'
                self.sess = tf.Session(graph=self.graph)
            self.saver = tf.train.Saver()
        self.reload_model(model_file)
        self._recData = RecDataBase(config, state_cls)

    def reset(self, reset_file=False):
        self._recData.reset(reset_file)

    @abstractmethod
    def get_model(self, config):
        pass

    @abstractmethod
    def obs_fn(self, state):
        pass

    @abstractmethod
    def forward(self, model, samples):
        pass

    def reload_model(self, model_file):
        with self.sess.as_default():
            with self.graph.as_default():
                self.saver.restore(self.sess, model_file)

    def seed(self, sd=0):
        self._recData.seed(sd)
        np.random.seed(sd)

    def _step(self, samples, action, **kwargs):
        step = kwargs['step']
        samples.act(action)
        next_state = samples.state
        next_obs = self.obs_fn(next_state)
        reward = self.forward(self.model, samples)
        next_info = samples.info

        if step < self.max_steps - 1:
            done = [0] * self.batch_size
        else:
            done = [1] * self.batch_size

        return next_obs, reward, done, next_info

    def sample(self, batch_size):
        samples = self._recData.sample(batch_size)
        obs = self.obs_fn(samples.state)
        return samples, obs

class RecEnvBase(gym.Env):
    metadata = {'render.modes': ['human']}

    def __init__(self, recsim: RecSimBase):
        self.config = recsim.config
        self.batch_size = self.config['batch_size']
        self.cur_step = 0
        self.sim = recsim
        self.sim.reset()
        self.samples, self.obs = self.sim.sample(self.batch_size)
        if self.config.get("rawstate_as_obs", False):
            category_size = len(self.obs[0]['category_feature'])
            dense_size = len(self.obs[0]['dense_feature'])
            sequence_size = np.array(self.obs[0]['sequence_feature']).shape
            features = {
                "category_feature": gym.spaces.Box(-1000000.0, 1000000.0, shape=(category_size,)),
                "dense_feature": gym.spaces.Box(-1000000.0, 1000000.0, shape=(dense_size,)),
                "sequence_feature": gym.spaces.Box(-1000000.0, 1000000.0, shape=sequence_size),
            }
            if self.config.get("support_rllib_mask", False):
                action_feature_size = len(self.obs[0]['action_mask'])
                self.observation_space = gym.spaces.Dict({
                    "action_mask": gym.spaces.Box(0, 1, shape=(action_feature_size,)),
                    **features
                })
            else:
                self.observation_space = gym.spaces.Dict(features)
        else:
            if self.config.get("support_rllib_mask", False):
                action_feature_size = len(self.obs[0]['action_mask'])
                self.observation_space = gym.spaces.Dict({
                    "action_mask": gym.spaces.Box(0, 1, shape=(action_feature_size,)),
                    "obs": gym.spaces.Box(-1000.0, 1000.0, shape=(len(self.obs[0]["obs"]),))
                })
            else:
                self.observation_space = gym.spaces.Box(-1000.0, 1000.0, shape=(len(self.obs[0]),))
        if self.config.get("support_conti_env", False):
            self.action_space = gym.spaces.Box(-1, 1, shape=(self.config['action_emb_size'],))
        else:
            self.action_space = gym.spaces.Discrete(self.config['action_size'])
        # if self.config.get("support_rllib_mask", False):
        #     action_feature_size = len(self.obs[0]['action_mask'])
        #     # avail_actions_size = len(self.obs[0]['avail_actions'][0])
        #     # self.action_space = gym.spaces.Discrete(self.config['action_size'])
        #     self.observation_space = gym.spaces.Dict({
        #         "action_mask": gym.spaces.Box(0, 1, shape=(action_feature_size,)),
        #         "obs": self.observation_space,
        #     })
        # elif self.config.get("support_d3rl_mask", False):
        #     self.action_space = gym.spaces.Discrete(self.config['action_size'])
        # else:
        #     self.action_space = gym.spaces.Discrete(self.config['action_size'])
        self.reset()

    def seed(self, sd=0):
        self.sim.seed(sd)
        np.random.seed(sd)

    @property
    @single_elem_support
    def state(self):
        return self.obs

    @property
    @single_elem_support
    def user_id(self):
        return self.samples.user

    @property
    @single_elem_support
    def offline_action(self):
        return self.samples.offline_action

    @property
    @single_elem_support
    def offline_reward(self):
        return self.samples.offline_reward

    @single_elem_support
    def step(self, action):
        if not isinstance(action, (list, np.ndarray)):
            action = [action]
        obs, reward, done, info = \
            self.sim._step(self.samples, action, step=self.cur_step)
        self.cur_step += 1
        return obs, reward, done, info

    def reset(self, reset_file=False):
        self.cur_step = 0
        self.sim.reset(reset_file)
        self.samples, self.obs = self.sim.sample(self.batch_size)
        return self.state

    def render(self, mode='human', close=False):
        print('Current State:', '\n')
        print(self.samples.to_string())


text:
# coding=utf-8
# Copyright 2020 The Uncertainty Baselines Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Lint as: python3
"""Tests for speech_commands."""

from absl.testing import parameterized
import numpy as np
import tensorflow.compat.v2 as tf
import uncertainty_baselines as ub
from uncertainty_baselines.datasets import base

class SpeechCommandsDatasetTest(tf.test.TestCase, parameterized.TestCase):

  # Deleting the following because it was failing tests:
  #    ('RoomReverb12M', ('room_reverb', 12))
  #    ('LowPassFiltering4kHz', ('low_pass', 4))
  @parameterized.named_parameters(
      ('Train', base.Split.TRAIN), ('Validation', base.Split.VAL),
      ('Test', base.Split.TEST), ('WhiteNoiseNegative5DB', ('white_noise', -5)),
      ('PitchShift4Semitones', ('pitch_shift', 4)))
  def testSpeechCommandsDataset(self, split):
    batch_size = 35
    eval_batch_size = 25
    dataset_builder = ub.datasets.SpeechCommandsDataset(
        batch_size=batch_size,
        eval_batch_size=eval_batch_size,
        shuffle_buffer_size=20)
    dataset = dataset_builder.build(split).take(1)
    element = next(iter(dataset))
    features = element['features']
    labels = element['labels']
    if split == base.Split.TRAIN:
      self.assertEqual(
          features.shape,
          (batch_size, ub.datasets.SpeechCommandsDataset.AUDIO_LENGTH))
      self.assertEqual(labels.shape, (batch_size,))
    else:
      self.assertEqual(
          features.shape,
          (eval_batch_size, ub.datasets.SpeechCommandsDataset.AUDIO_LENGTH))
      self.assertEqual(labels.shape, (eval_batch_size,))
    # The in-distrubtions labels go up to 10 only.
    self.assertFalse(np.any(labels > 10))

    # Uncomment the code here to generate sample .wav files for inspection.
    # from scipy.io import wavfile
    # import tempfile
    # if isinstance(split, tuple):
    #   for i, label in enumerate(labels.numpy()):
    #     wav_file_path = tempfile.mktemp() + '_%s_label_%d.wav' % (
    #         split[0], label)
    #     wavfile.write(wav_file_path, 16000, features.numpy()[i] / 32768)
    #     print('Wrote to %s' % wav_file_path)

  def testSpeechCommandsDatasetSemanticShiftSplit(self):
    batch_size = 7
    eval_batch_size = 5
    dataset_builder = ub.datasets.SpeechCommandsDataset(
        batch_size=batch_size,
        eval_batch_size=eval_batch_size,
        shuffle_buffer_size=20)
    dataset = dataset_builder.build(('semantic_shift',)).take(1)
    element = next(iter(dataset))
    features = element['features']
    labels = element['labels']

    self.assertEqual(
        features.shape,
        (eval_batch_size, ub.datasets.SpeechCommandsDataset.AUDIO_LENGTH))
    self.assertEqual(labels.shape, (eval_batch_size,))
    # The semantic shift split has only the label 11, which differs from the
    # non-semantic-shifted labels: 0 - 10.
    self.assertFalse(np.any(labels != 11))

  def testInvalidSplitKeyNameLeadsToException(self):
    dataset_builder = ub.datasets.SpeechCommandsDataset(
        batch_size=7, eval_batch_size=5, shuffle_buffer_size=20)
    split = ('nonsensical_split', 42)
    with self.assertRaisesRegex(ValueError,
                                r'Unrecognized shift split: nonsensical_split'):
      dataset_builder.build(split)

if __name__ == '__main__':
  tf.test.main()

max_stars_repo_name: gpleiss/uncertainty-baselines

max_issues_repo_name: gpleiss/uncertainty-baselines

max_forks_repo_name: gpleiss/uncertainty-baselines


text:
###################3###################3###################3##################
# interpSphericalPotential.py: build spherical potential through interpolation
###################3###################3###################3##################
import numpy
from scipy import interpolate
from .SphericalPotential import SphericalPotential
from .Potential import _evaluateRforces, _evaluatePotentials
from ..util.conversion import physical_compatible, get_physical
class interpSphericalPotential(SphericalPotential):
    """__init__(self,rforce=None,rgrid=numpy.geomspace(0.01,20,101),Phi0=None,ro=None,vo=None)

Class that interpolates a spherical potential on a grid"""
    def __init__(self,rforce=None,rgrid=numpy.geomspace(0.01,20,101),Phi0=None,
                 ro=None,vo=None):
        """__init__(self,rforce=None,rgrid=numpy.geomspace(0.01,20,101),Phi0=None,ro=None,vo=None)

        NAME:

           __init__

        PURPOSE:

           initialize an interpolated, spherical potential

        INPUT:

           rforce= (None) Either a) function that gives the radial force as a function of r or b) a galpy Potential instance or list thereof

           rgrid= (numpy.geomspace(0.01,20,101)) radial grid on which to evaluate the potential for interpolation (note that beyond rgrid[-1], the potential is extrapolated as -GM(<rgrid[-1])/r)

           Phi0= (0.) value of the potential at rgrid[0] (only necessary when rforce is a function, for galpy potentials automatically determined)

           ro=, vo= distance and velocity scales for translation into internal units (default from configuration file)

        OUTPUT:

           (none)

        HISTORY:

           2020-07-13 - Written - Bovy (UofT)

        """
        SphericalPotential.__init__(self,amp=1.,ro=ro,vo=vo)
        self._rgrid= rgrid
        # Determine whether rforce is a galpy Potential or list thereof
        try:
            _evaluateRforces(rforce,1.,0.)
        except:
            _rforce= rforce
            Phi0= 0. if Phi0 is None else Phi0
        else:
            _rforce= lambda r: _evaluateRforces(rforce,r,0.)
            # Determine Phi0
            Phi0= _evaluatePotentials(rforce,rgrid[0],0.)
            # Also check that unit systems are compatible
            if not physical_compatible(self,rforce):
                raise RuntimeError('Unit conversion factors ro and vo incompatible between Potential to be interpolated and the factors given to interpSphericalPotential')
            # If set for the parent, set for the interpolated
            phys= get_physical(rforce,include_set=True)
            if phys['roSet']:
                self.turn_physical_on(ro=phys['ro'])
            if phys['voSet']:
                self.turn_physical_on(vo=phys['vo'])
        self._rforce_grid= numpy.array([_rforce(r) for r in rgrid])
        self._force_spline= interpolate.InterpolatedUnivariateSpline(
            self._rgrid,self._rforce_grid,k=3,ext=0)
        # Get potential and r2deriv as splines for the integral and derivative
        self._pot_spline= self._force_spline.antiderivative()
        self._Phi0= Phi0+self._pot_spline(self._rgrid[0])
        self._r2deriv_spline= self._force_spline.derivative()
        # Extrapolate as mass within rgrid[-1]
        self._rmin= rgrid[0]
        self._rmax= rgrid[-1]
        self._total_mass= -self._rmax**2.*self._force_spline(self._rmax)
        self._Phimax= -self._pot_spline(self._rmax)+self._Phi0\
            +self._total_mass/self._rmax
        self.hasC= True
        self.hasC_dxdv= True
        self.hasC_dens= True
        return None

    def _revaluate(self,r,t=0.):
        out= numpy.empty_like(r)
        out[r >= self._rmax]= -self._total_mass/r[r >= self._rmax]+self._Phimax
        out[r < self._rmax]= -self._pot_spline(r[r < self._rmax])+self._Phi0
        return out

    def _rforce(self,r,t=0.):
        out= numpy.empty_like(r)
        out[r >= self._rmax]= -self._total_mass/r[r >= self._rmax]**2.
        out[r < self._rmax]= self._force_spline(r[r < self._rmax])
        return out

    def _r2deriv(self,r,t=0.):
        out= numpy.empty_like(r)
        out[r >= self._rmax]= -2.*self._total_mass/r[r >= self._rmax]**3.
        out[r < self._rmax]= -self._r2deriv_spline(r[r < self._rmax])
        return out

    def _rdens(self,r,t=0.):
        out= numpy.empty_like(r)
        out[r >= self._rmax]= 0.
        # Fall back onto Poisson eqn., implemented in SphericalPotential
        out[r < self._rmax]= SphericalPotential._rdens(self,r[r < self._rmax])
        return out

max_stars_repo_name: gusbeane/galpy

max_issues_repo_name: gusbeane/galpy

max_forks_repo_name: gusbeane/galpy