<commit_before>import numpy as np
import scipy
import statsmodels.api as stats


def make_pixel_mask(image, saturation, input_mask=None):
    """Make a pixel mask that marks saturated pixels; optionally join with input_mask"""

    if input_mask is not None:
        new_mask = input_mask
    else:
        new_mask = np.zeros(image.shape)

    new_mask[image >= saturation] = 1
    return new_mask


def center_psf(psf):
    """Center psf at (0,0) based on max value"""

    peak = np.unravel_index(psf.argmax(), psf.shape)

    psf = np.roll(psf, -peak[0], 0)
    psf = np.roll(psf, -peak[1], 1)
    return psf


def fit_noise(data, n_stamps=1, mode='iqr'):
    """Find the standard deviation of the image background; returns standard deviation, median"""

    median_small = np.zeros([n_stamps, n_stamps])
    std_small = np.zeros([n_stamps, n_stamps])
    for y_stamp in range(n_stamps):
        for x_stamp in range(n_stamps):
            y_index = [y_stamp * data.shape[0] // n_stamps, (y_stamp + 1) * data.shape[0] // n_stamps]
            x_index = [x_stamp * data.shape[1] // n_stamps, (x_stamp + 1) * data.shape[1] // n_stamps]
            stamp_data = data[y_index[0]: y_index[1], x_index[0]: x_index[1]]
            if mode == 'gaussian':
                trimmed_stamp_data = stamp_data[stamp_data < np.percentile(stamp_data, 90)]
                trimmed_stamp_data = trimmed_stamp_data[trimmed_stamp_data != 0]
                histogram_data = np.histogram(trimmed_stamp_data, bins=100)
                x = histogram_data[1][:-1]
                y = histogram_data[0]
                guess = [np.max(y), np.median(trimmed_stamp_data), np.std(trimmed_stamp_data)]
                parameters, covariance = scipy.optimize.curve_fit(gauss, x, y, p0=guess, maxfev=1600)
                median_small[y_stamp, x_stamp] = parameters[1]
                std_small[y_stamp, x_stamp] = parameters[2]
            elif mode == 'iqr':
                quartile25, median, quartile75 = np.percentile(data, (25, 50, 75))
                median_small[y_stamp, x_stamp] = median
                # 0.741301109 is a tuning parameter that scales iqr to std
                std_small[y_stamp, x_stamp] = 0.741301109 * (quartile75 - quartile25)

    median = scipy.ndimage.zoom(median_small, [data.shape[0] / float(n_stamps), data.shape[1] / float(n_stamps)])
    std = scipy.ndimage.zoom(std_small, [data.shape[0] / float(n_stamps), data.shape[1] / float(n_stamps)])

    return std, median


def gauss(position, amplitude, median, std):
    """Return a gaussian function"""

    return amplitude * np.exp(-(position - median) ** 2 / (2 * std ** 2))


def interpolate_bad_pixels(image, mask, median_size=6):
    """Interpolate over bad pixels using a global median"""

    # from http://stackoverflow.com/questions/18951500/automatically-remove-hot-dead-pixels-from-an-image-in-python
    pix = np.transpose(np.where(mask == 1))
    blurred = scipy.ndimage.median_filter(image, size=median_size)
    interpolated_image = np.copy(image)
    for y, x in pix:
        interpolated_image[y, x] = blurred[y, x]

    return interpolated_image


def resize_psf(psf, shape):
    """Resize centered (0,0) psf to larger shape"""

    psf_extended = np.pad(psf, ((0, shape[0] - psf.shape[0]), (0, shape[1] - psf.shape[1])),
                          mode='constant', constant_values=0.)
    return psf_extended


def pad_to_power2(data):
    """Pad arrays to the nearest power of two"""

    n = 0
    defecit = [0, 0]
    while (data.shape[0] > (2 ** n)) or (data.shape[1] > (2 ** n)):
        n += 1
        defecit = [(2 ** n) - data.shape[0], (2 ** n) - data.shape[1]]
    padded_data = np.pad(data, ((0, defecit[0]), (0, defecit[1])), mode='constant', constant_values=0.)
    return padded_data


def solve_iteratively(science, reference,
                      mask_tolerance=10e-5, gain_tolerance=10e-6, max_iterations=5, sigma_cut=1):
    """Solve for linear fit iteratively"""

    gain = 1.
    gain0 = 10e5
    i = 0

    # pad image to power of two to speed fft
    old_size = science.image_data.shape
    science_image = pad_to_power2(science.image_data)
    reference_image = pad_to_power2(reference.image_data)

    science_psf = center_psf(resize_psf(science.raw_psf_data, science_image.shape))
    science_psf /= np.sum(science.raw_psf_data)
    reference_psf = center_psf(resize_psf(reference.raw_psf_data, reference_image.shape))
    reference_psf /= np.sum(reference.raw_psf_data)

    science_std = pad_to_power2(science.background_std)
    reference_std = pad_to_power2(reference.background_std)

    science_mask = pad_to_power2(science.pixel_mask)
    reference_mask = pad_to_power2(reference.pixel_mask)
<commit_msg>fixed divide by zero bug cause by padding to power of 2<commit_after>import numpy as np
import scipy
import statsmodels.api as stats


def make_pixel_mask(image, saturation, input_mask=None):
    """Make a pixel mask that marks saturated pixels; optionally join with input_mask"""

    if input_mask is not None:
        new_mask = input_mask
    else:
        new_mask = np.zeros(image.shape)

    new_mask[image >= saturation] = 1
    return new_mask


def center_psf(psf):
    """Center psf at (0,0) based on max value"""

    peak = np.unravel_index(psf.argmax(), psf.shape)

    psf = np.roll(psf, -peak[0], 0)
    psf = np.roll(psf, -peak[1], 1)
    return psf


def fit_noise(data, n_stamps=1, mode='iqr'):
    """Find the standard deviation of the image background; returns standard deviation, median"""

    median_small = np.zeros([n_stamps, n_stamps])
    std_small = np.zeros([n_stamps, n_stamps])
    for y_stamp in range(n_stamps):
        for x_stamp in range(n_stamps):
            y_index = [y_stamp * data.shape[0] // n_stamps, (y_stamp + 1) * data.shape[0] // n_stamps]
            x_index = [x_stamp * data.shape[1] // n_stamps, (x_stamp + 1) * data.shape[1] // n_stamps]
            stamp_data = data[y_index[0]: y_index[1], x_index[0]: x_index[1]]
            if mode == 'gaussian':
                trimmed_stamp_data = stamp_data[stamp_data < np.percentile(stamp_data, 90)]
                trimmed_stamp_data = trimmed_stamp_data[trimmed_stamp_data != 0]
                histogram_data = np.histogram(trimmed_stamp_data, bins=100)
                x = histogram_data[1][:-1]
                y = histogram_data[0]
                guess = [np.max(y), np.median(trimmed_stamp_data), np.std(trimmed_stamp_data)]
                parameters, covariance = scipy.optimize.curve_fit(gauss, x, y, p0=guess, maxfev=1600)
                median_small[y_stamp, x_stamp] = parameters[1]
                std_small[y_stamp, x_stamp] = parameters[2]
            elif mode == 'iqr':
                quartile25, median, quartile75 = np.percentile(data, (25, 50, 75))
                median_small[y_stamp, x_stamp] = median
                # 0.741301109 is a tuning parameter that scales iqr to std
                std_small[y_stamp, x_stamp] = 0.741301109 * (quartile75 - quartile25)

    median = scipy.ndimage.zoom(median_small, [data.shape[0] / float(n_stamps), data.shape[1] / float(n_stamps)])
    std = scipy.ndimage.zoom(std_small, [data.shape[0] / float(n_stamps), data.shape[1] / float(n_stamps)])

    return std, median


def gauss(position, amplitude, median, std):
    """Return a gaussian function"""

    return amplitude * np.exp(-(position - median) ** 2 / (2 * std ** 2))


def interpolate_bad_pixels(image, mask, median_size=6):
    """Interpolate over bad pixels using a global median"""

    # from http://stackoverflow.com/questions/18951500/automatically-remove-hot-dead-pixels-from-an-image-in-python
    pix = np.transpose(np.where(mask == 1))
    blurred = scipy.ndimage.median_filter(image, size=median_size)
    interpolated_image = np.copy(image)
    for y, x in pix:
        interpolated_image[y, x] = blurred[y, x]

    return interpolated_image


def resize_psf(psf, shape):
    """Resize centered (0,0) psf to larger shape"""

    psf_extended = np.pad(psf, ((0, shape[0] - psf.shape[0]), (0, shape[1] - psf.shape[1])),
                          mode='constant', constant_values=0.)
    return psf_extended


def pad_to_power2(data):
    """Pad arrays to the nearest power of two"""

    n = 0
    defecit = [0, 0]
    while (data.shape[0] > (2 ** n)) or (data.shape[1] > (2 ** n)):
        n += 1
        defecit = [(2 ** n) - data.shape[0], (2 ** n) - data.shape[1]]
    padded_data = np.pad(data, ((0, defecit[0]), (0, defecit[1])), mode='constant', constant_values=np.median(data))
    return padded_data


def solve_iteratively(science, reference,
                      mask_tolerance=10e-5, gain_tolerance=10e-6, max_iterations=5, sigma_cut=1):
    """Solve for linear fit iteratively"""

    gain = 1.
    gain0 = 10e5
    i = 0
    # pad image to power of two to speed fft
    old_size = science.image_data.shape
    science_image = pad_to_power2(science.image_data)
    reference_image = pad_to_power2(reference.image_data)

    science_psf = center_psf(resize_psf(science.raw_psf_data, science_image.shape))
    science_psf /= np.sum(science.raw_psf_data)
    reference_psf = center_psf(resize_psf(reference.raw_psf_data, reference_image.shape))
    reference_psf /= np.sum(reference.raw_psf_data)

    science_std = pad_to_power2(science.background_std)
    reference_std = pad_to_power2(reference.background_std)

    science_mask = pad_to_power2(science.pixel_mask)
    reference_mask = pad_to_power2(reference.pixel_mask)

{'fraction_non_alphanumeric': 0.1262381454162276, 'fraction_numerical': 0.027818756585879872, 'mean_word_length': 3.089185695820767, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1622815', 'n_tokens_mistral': 3591, 'n_tokens_neox': 3325, 'n_words': 725}