text:
import numpy as np
from skimage.transform import resize
from keras.optimizers import Adam
from keras.models import Model, Sequential
from keras.layers import Input, InputLayer, Dense, BatchNormalization, Dropout, Flatten, Reshape

class Autoencoder:
    def __init__(self, img_shape=(90, 90, 3), code_size=128, hidden_size=512,
                 encoder_weights_path = None, decoder_weights_path = None,
                 happiness_code_path = None):
        # Basic input parameters
        self.img_shape = img_shape
        self.code_size = code_size
        self.hidden_size = hidden_size
        self.encoder_weights_path = encoder_weights_path
        self.decoder_weights_path = decoder_weights_path
        self.happiness_code_path = happiness_code_path

        # Compile model
        self.encoder, self.decoder = self.build_autoencoder(img_shape, code_size, hidden_size)
        inp = Input(self.img_shape)
        code = self.encoder(inp)
        reconstruction = self.decoder(code)
        self.autoencoder = Model(inputs=inp, outputs=reconstruction)
        self.autoencoder.compile(optimizer=Adam(lr=0.0001), loss='binary_crossentropy')

        # Load weights and happiness code if needed
        if self.encoder_weights_path is not None and self.decoder_weights_path is not None:
            # load weights
            self.load_weights(self.encoder_weights_path, self.decoder_weights_path)
            if self.happiness_code_path is not None:
              self.happiness_code = np.load(self.happiness_code_path)

    # \brief Fit method wrapper
    # Note that y_train == x_train, y_val == x_val, because this is unsupervised learning model
    def fit(self, x_train, x_val, epochs, batch_size=128, verbose=1):
        self.autoencoder.fit(x=x_train, y=x_train, epochs=epochs, batch_size=batch_size,
                             validation_data=[x_val, x_val], verbose=verbose)

    # \brief Predict method wrapper
    def predict(self, data):
        assert np.min(data) >= 0 and np.max(data) <= 1
        return self.autoencoder.predict(data)

    # \brief Single image predict wrapper
    def predict_img(self, img):
        assert np.min(img) >= 0 and np.max(img) <= 1
        # img[None] is the same as img[np.newaxis, :]
        return self.autoencoder.predict(img[None])[0]

    # \brief Single code predict wrapper
    def predict_img_code(self, img):
        assert np.min(img) >= 0 and np.max(img) <= 1
        # img[None] is the same as img[np.newaxis, :]
        return self.encoder.predict(img[None])[0]

    # \brief Single code reconstruction predict wrapper
    def _predict_code_reconstruction(self, code):
        # img[None] is the same as img[np.newaxis, :]
        return self.decoder.predict(code[None])[0]

    # \brief Load pretrained model weights
    def load_weights(self, encoder_weights_path, decoder_weights_path):
        self.encoder.load_weights(encoder_weights_path)
        self.decoder.load_weights(decoder_weights_path)

    # \brief Generate happy code from img code
    def _get_happy_img_code(self, img_code, inverse):
        if inverse:
          happy_img_code = img_code - self.happiness_code
        else:
          happy_img_code = img_code + self.happiness_code
        return happy_img_code

    # \brief Add happiness to image
    def _add_happiness(self, img, inverse):
        img_code = self.encoder.predict(img[None])[0]
        happy_img_code = self._get_happy_img_code(img_code, inverse)
        happy_img = self.decoder.predict(happy_img_code[None])[0]
        return happy_img

    # \brief Add happiness wrapper
    def add_happiness(self, raw_img, inverse=False):
        formatted_img = self.prepare_photo_before_feeding(raw_img)
        reconstr_img = self._add_happiness(formatted_img, inverse)
        return reconstr_img

    # \brief Photo must suit the input format of network
    def prepare_photo_before_feeding(self, raw_jpg):
        dx, dy = self.img_shape[0], self.img_shape[1]
        img = resize(raw_jpg, [dx, dy], mode='reflect', anti_aliasing=True)
        # this method automatically normalizes by 255.
        # dimx, dimy = 3 * dx, 3 * dy
        #img = resize(raw_jpg, [dimx, dimy], mode='reflect', anti_aliasing=True)
        # crop the center part
        #img = img[dy:-dy, dx:-dx]
        # convert to float32, maybe redundant
        img = img.astype('float32')
        return img

    # \brief Feed network with raw_jpg
    def feed_photo(self, raw_jpg):
        formatted_img = self.prepare_photo_before_feeding(raw_jpg)
        reconstr_img = self.predict_img(formatted_img)
        return formatted_img, reconstr_img

    # \brief Deep autoencoder network construction
    # \return encoder -- encoding part, which compress image to the code
    # \return decoder -- decoding part, which reconstruct image for input code
    @staticmethod
    def build_autoencoder(img_shape=(90, 90, 3), code_size=128, hidden_size=512):
        # encoder
        encoder = Sequential()
        encoder.add(InputLayer(img_shape))
        encoder.add(Flatten())
        encoder.add(Dense(hidden_size * 2, activation='relu'))
        encoder.add(BatchNormalization())
        encoder.add(Dropout(0.3))
        encoder.add(Dense(hidden_size, activation='relu'))
        encoder.add(BatchNormalization())
        encoder.add(Dropout(0.3))
        encoder.add(Dense(code_size, activation='relu'))
        # decoder
        decoder = Sequential()
        decoder.add(InputLayer((code_size,)))
        decoder.add(Dense(hidden_size, activation='relu'))
        decoder.add(BatchNormalization())
        decoder.add(Dropout(0.3))
        decoder.add(Dense(hidden_size * 2, activation='relu'))
        decoder.add(BatchNormalization())
        decoder.add(Dropout(0.3))
        decoder.add(Dense(img_shape[0] * img_shape[1] * img_shape[2], activation='sigmoid'))
        decoder.add(Reshape(img_shape))
        return encoder, decoder

max_stars_repo_name: GimmeDanger/autoencoders-bot

max_issues_repo_name: GimmeDanger/deep-autoencoder-bot

max_forks_repo_name: GimmeDanger/deep-autoencoder-bot


text:
import unittest
import numpy

from qupulse.pulses.plotting import PlottingNotPossibleException, render, iter_waveforms, iter_instruction_block, plot
from qupulse._program.instructions import InstructionBlock
from qupulse.pulses.table_pulse_template import TablePulseTemplate
from qupulse.pulses.sequence_pulse_template import SequencePulseTemplate
from qupulse.pulses.sequencing import Sequencer
from qupulse._program._loop import MultiChannelProgram, Loop

from tests.pulses.sequencing_dummies import DummyWaveform, DummyInstruction, DummyPulseTemplate

class PlotterTests(unittest.TestCase):

    def test_render_unsupported_instructions(self) -> None:
        block = InstructionBlock()
        block.add_instruction(DummyInstruction())

        with self.assertRaises(NotImplementedError):
            render(block)

    def test_render_no_waveforms(self) -> None:
        time, channel_data = render(InstructionBlock())
        self.assertEqual(channel_data, dict())
        numpy.testing.assert_equal(time, numpy.empty(0))

    def test_iter_waveforms(self) -> None:
        wf1 = DummyWaveform(duration=7)
        wf2 = DummyWaveform(duration=5)
        wf3 = DummyWaveform(duration=3)

        repeated_block = InstructionBlock()
        repeated_block.add_instruction_meas([('m', 1, 2)])
        repeated_block.add_instruction_exec(wf2)
        repeated_block.add_instruction_exec(wf1)

        main_block = InstructionBlock()
        main_block.add_instruction_exec(wf1)
        main_block.add_instruction_repj(2, repeated_block)
        main_block.add_instruction_exec(wf3)

        for idx, (expected, received) in enumerate(zip([wf1, wf2, wf1, wf2, wf1, wf3], iter_waveforms(main_block))):
            self.assertIs(expected, received, msg="Waveform {} is wrong".format(idx))

    def test_iter_waveform_exceptions(self) -> None:
        wf1 = DummyWaveform(duration=7)
        wf2 = DummyWaveform(duration=5)
        wf3 = DummyWaveform(duration=3)

        repeated_block = InstructionBlock()
        repeated_block.add_instruction_meas([('m', 1, 2)])
        repeated_block.add_instruction_exec(wf2)
        repeated_block.add_instruction_exec(wf1)

        main_block = InstructionBlock()
        main_block.add_instruction_exec(wf1)
        main_block.add_instruction_repj(2, repeated_block)
        main_block.add_instruction_exec(wf3)
        main_block.add_instruction_goto(repeated_block)

        with self.assertRaises(NotImplementedError):
            list(iter_waveforms(main_block))

        repeated_block.add_instruction(DummyInstruction())
        with self.assertRaises(NotImplementedError):
            list(iter_waveforms(main_block))

        main_block = InstructionBlock()
        main_block.add_instruction_stop()

        with self.assertRaises(StopIteration):
            next(iter_waveforms(main_block))

    def test_iter_instruction_block(self) -> None:
        wf1 = DummyWaveform(duration=7)
        wf2 = DummyWaveform(duration=5)
        wf3 = DummyWaveform(duration=3)

        repeated_block = InstructionBlock()
        repeated_block.add_instruction_meas([('m', 1, 2)])
        repeated_block.add_instruction_exec(wf2)
        repeated_block.add_instruction_exec(wf1)

        main_block = InstructionBlock()
        main_block.add_instruction_exec(wf1)
        main_block.add_instruction_repj(2, repeated_block)
        main_block.add_instruction_exec(wf3)

        waveforms, measurements, total_time = iter_instruction_block(main_block, True)

        for idx, (expected, received) in enumerate(zip([wf1, wf2, wf1, wf2, wf1, wf3], waveforms)):
            self.assertIs(expected, received, msg="Waveform {} is wrong".format(idx))
        self.assertEqual([('m', 8, 2), ('m', 20, 2)], measurements)
        self.assertEqual(total_time, 34)

    def test_iter_instruction_block_exceptions(self) -> None:
        wf1 = DummyWaveform(duration=7)
        wf2 = DummyWaveform(duration=5)
        wf3 = DummyWaveform(duration=3)

        repeated_block = InstructionBlock()
        repeated_block.add_instruction_meas([('m', 1, 2)])
        repeated_block.add_instruction_exec(wf2)

        main_block = InstructionBlock()
        main_block.add_instruction_exec(wf1)
        main_block.add_instruction_repj(2, repeated_block)
        main_block.add_instruction_exec(wf3)

        repeated_block.add_instruction_goto(main_block)

        with self.assertRaises(NotImplementedError):
            iter_instruction_block(main_block, False)

        repeated_block = InstructionBlock()
        repeated_block.add_instruction_meas([('m', 1, 2)])
        repeated_block.add_instruction_exec(wf2)
        repeated_block.add_instruction(DummyInstruction())

        main_block = InstructionBlock()
        main_block.add_instruction_exec(wf1)
        main_block.add_instruction_repj(2, repeated_block)
        main_block.add_instruction_exec(wf3)

        with self.assertRaises(NotImplementedError):
            iter_instruction_block(main_block, False)

    def test_render(self) -> None:
        with self.assertWarnsRegex(DeprecationWarning, ".*InstructionBlock.*"):
            wf1 = DummyWaveform(duration=19)
            wf2 = DummyWaveform(duration=21)

            block = InstructionBlock()
            block.add_instruction_exec(wf1)
            block.add_instruction_meas([('asd', 0, 1)])
            block.add_instruction_exec(wf2)

            wf1_expected = ('A', [0, 2, 4, 6, 8, 10, 12, 14, 16, 18])
            wf2_expected = ('A', [x-19 for x in [20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40]])
            wf1_output_array_len_expected = len(wf1_expected[1])
            wf2_output_array_len_expected = len(wf2_expected[1])

            wf1.sample_output = numpy.linspace(start=4, stop=5, num=len(wf1_expected[1]))
            wf2.sample_output = numpy.linspace(6, 7, num=len(wf2_expected[1]))

            expected_times = numpy.arange(start=0, stop=42, step=2)
            expected_result = numpy.concatenate((wf1.sample_output, wf2.sample_output))

            times, voltages, _ = render(block, sample_rate=0.5)

            self.assertEqual(len(wf1.sample_calls), 1)
            self.assertEqual(len(wf2.sample_calls), 1)

            self.assertEqual(wf1_expected[0], wf1.sample_calls[0][0])
            self.assertEqual(wf2_expected[0], wf2.sample_calls[0][0])

            numpy.testing.assert_almost_equal(wf1_expected[1], wf1.sample_calls[0][1])
            numpy.testing.assert_almost_equal(wf2_expected[1], wf2.sample_calls[0][1])

            self.assertEqual(wf1_output_array_len_expected, len(wf1.sample_calls[0][2]))
            self.assertEqual(wf2_output_array_len_expected, len(wf2.sample_calls[0][2]))

            self.assertEqual(voltages.keys(), dict(A=0).keys())

            numpy.testing.assert_almost_equal(expected_times, times)
            numpy.testing.assert_almost_equal(expected_result, voltages['A'])
            self.assertEqual(expected_result.shape, voltages['A'].shape)

            times, voltages, measurements = render(block, sample_rate=0.5, render_measurements=True)
            self.assertEqual(voltages.keys(), dict(A=0).keys())

            numpy.testing.assert_almost_equal(expected_times, times)
            numpy.testing.assert_almost_equal(expected_result, voltages['A'])
            self.assertEqual(expected_result.shape, voltages['A'].shape)

            self.assertEqual(measurements, [('asd', 19, 1)])

    def test_render_block_time_slice(self) -> None:
        with self.assertWarnsRegex(DeprecationWarning, ".*InstructionBlock.*"):
            with self.assertRaises(ValueError):
                wf1 = DummyWaveform(duration=19)
                wf2 = DummyWaveform(duration=21)

                block = InstructionBlock()
                block.add_instruction_exec(wf1)
                block.add_instruction_exec(wf2)

                times, voltages, _ = render(block, sample_rate=0.5, time_slice=(1, 16))

    def test_render_loop_compare(self) -> None:
        wf1 = DummyWaveform(duration=19)
        wf2 = DummyWaveform(duration=21)

        block = InstructionBlock()
        block.add_instruction_exec(wf1)
        block.add_instruction_meas([('asd', 0, 1), ('asd', 1, 1)])
        block.add_instruction_exec(wf2)

        mcp = MultiChannelProgram(block)
        loop = next(iter(mcp.programs.values()))

        block_times, block_voltages, _ = render(block, sample_rate=0.5)
        loop_times, loop_voltages, _ = render(loop, sample_rate=0.5)

        numpy.testing.assert_equal(block_times, loop_times)
        numpy.testing.assert_equal(block_voltages, loop_voltages)

        block_times, block_voltages, block_measurements = render(block, sample_rate=0.5, render_measurements=True)
        loop_times, loop_voltages, loop_measurements = render(loop, sample_rate=0.5, render_measurements=True)
        numpy.testing.assert_equal(block_times, loop_times)
        numpy.testing.assert_equal(block_voltages, loop_voltages)
        self.assertEqual(block_measurements, loop_measurements)

    def test_render_loop_sliced(self) -> None:
        wf = DummyWaveform(duration=19)

        full_times = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
        wf.sample_output = numpy.linspace(start=4, stop=5, num=len(full_times))[3:7]

        full_measurements = [('foo', 1, 3), ('bar', 5, 2), ('foo', 7, 3), ('bar', 11, 2), ('foo', 15, 3)]
        loop = Loop(waveform=wf, measurements=full_measurements)

        expected_times = full_times[3:7]
        expected_measurements = full_measurements[1:4]

        times, voltages, measurements = render(loop, sample_rate=0.5, render_measurements=True, time_slice=(6, 12))

        numpy.testing.assert_almost_equal(expected_times, times)
        numpy.testing.assert_almost_equal(wf.sample_output, voltages['A'])
        self.assertEqual(expected_measurements, measurements)

    def test_render_loop_invalid_slice(self) -> None:
        with self.assertRaises(ValueError):
            render(Loop(waveform=DummyWaveform()), time_slice=(5, 1))
            render(Loop(waveform=DummyWaveform()), time_slice=(-1, 4))
            render(Loop(waveform=DummyWaveform()), time_slice=(-7, -3))

    def test_render_warning(self) -> None:
        wf1 = DummyWaveform(duration=19)
        wf2 = DummyWaveform(duration=21)

        block = InstructionBlock()
        block.add_instruction_exec(wf1)
        block.add_instruction_meas([('asd', 0, 1)])
        block.add_instruction_exec(wf2)

        with self.assertWarns(UserWarning):
            render(block, sample_rate=0.51314323423)

    def integrated_test_with_sequencer_and_pulse_templates(self) -> None:
        # Setup test data
        square = TablePulseTemplate()
        square.add_entry('up', 'v', 'hold')
        square.add_entry('down', 0, 'hold')
        square.add_entry('length', 0)

        mapping1 = {
            'up': 'uptime',
            'down': 'uptime + length',
            'v': 'voltage',
            'length': '0.5 * pulse_length'
        }

        outer_parameters = ['uptime', 'length', 'pulse_length', 'voltage']

        parameters = {}
        parameters['uptime'] = 5
        parameters['length'] = 10
        parameters['pulse_length'] = 100
        parameters['voltage'] = 10

        sequence = SequencePulseTemplate([(square, mapping1), (square, mapping1)], outer_parameters)

        # run the sequencer and render the plot
        sample_rate = 20
        sequencer = Sequencer()
        sequencer.push(sequence, parameters)
        block = sequencer.build()
        times, voltages = render(block, sample_rate=sample_rate)

        # compute expected values
        expected_times = numpy.linspace(0, 100, sample_rate)
        expected_voltages = numpy.zeros_like(expected_times)
        expected_voltages[100:300] = numpy.ones(200) * parameters['voltage']

        # compare
        self.assertEqual(expected_times, times)
        self.assertEqual(expected_voltages, voltages)

    def test_plot_empty_pulse(self) -> None:
        import matplotlib
        matplotlib.use('svg') # use non-interactive backend so that test does not fail on travis

        pt = DummyPulseTemplate()
        with self.assertWarnsRegex(UserWarning, "empty", msg="plot() did not issue a warning for an empty pulse"):
            plot(pt, dict(), show=False)

class PlottingNotPossibleExceptionTests(unittest.TestCase):

    def test(self) -> None:
        t = DummyPulseTemplate()
        exception = PlottingNotPossibleException(t)
        self.assertIs(t, exception.pulse)
        self.assertIsInstance(str(exception), str)

ext: py

max_stars_repo_name: lankes-fzj/qupulse

max_issues_repo_name: lankes-fzj/qupulse

max_forks_repo_name: lankes-fzj/qupulse

reason: import numpy


text:
import os,sys
# import tensorflow as tf
import numpy as np
import csv

def generate(shape, rank):
    U = np.random.randn(shape[0], rank)
    V = np.random.randn(shape[1], rank)
    return U, V

def centralization(mat):
    tmp = np.matmul(np.ones((mat.shape[0], mat.shape[0])), mat)/mat.shape[0]
    centermat = mat-tmp
    return centermat

def subspace(shape, rank, mode=None):
    U, V = generate(shape, rank)
    tmp = np.matmul(U, np.transpose(V))
    if mode == 'B':
        Psb = centralization(tmp)
        return Psb
    if mode == 'C':
        Psc = centralization(tmp)
        return Psc
    elif mode == 'A':
        row = tmp.shape[0]
        col = tmp.shape[1]
        vec1 = np.ones((row, 1))
        vec1_t = np.transpose(vec1)
        vec2 = np.ones((col, 1))
        vec2_t = np.transpose(vec2)
        Psa = centralization(tmp) + np.matmul(np.matmul(vec1_t, tmp), vec2) * (vec1 * vec2_t) / (row * col)
        '''
        ss = centralization(tmp)
        a = np.matmul(vec1_t, tmp)
        b = np.matmul(a, vec2)
        c = b * (vec1 * vec2_t)
        d = c / (row * col)
        e = ss+d
        return e
        '''
        return Psa
    return False

def sample_rule4mat(tensor_shape, ra, rb, rc, sample_num):
    A = subspace((tensor_shape[0], tensor_shape[1]), ra, 'A')
    B = subspace((tensor_shape[1], tensor_shape[2]), rb, 'B')
    C = subspace((tensor_shape[2], tensor_shape[0]), rc, 'C')
    return A, B, C

def sample3D_rule(tensor_shape, sample_num):
    len_1d = tensor_shape[0]
    len_2d = tensor_shape[1]
    len_3d = tensor_shape[2]
    a = np.random.randint(0, len_1d, sample_num)
    b = np.random.randint(0, len_2d, sample_num)
    c = np.random.randint(0, len_3d, sample_num)
    return a, b, c

def Pomega_mat(spl, mat, tensor_shape, m, dim=None):
    # delta_vec = np.eye(m, dtype=int)
    #sum_vec = 0
    sum_vec = []
    # a, b, c = sample3D_rule(tensor_shape, m)
    if dim == 0:
        X = mat
        for i in range(m):
            # sum_vec += X[a[i]][b[i]] * delta_vec[:, i]
            sum_vec.append(X[spl[0][i]][spl[1][i]])
        sum_vec = np.asarray(sum_vec)
        Pomega_AX = sum_vec/np.sqrt(tensor_shape[2])
        return Pomega_AX
    elif dim == 1:
        Y = mat
        for i in range(m):
            # sum_vec += Y[b[i]][c[i]] * delta_vec[:, i]
            sum_vec.append(Y[spl[1][i]][spl[2][i]])
        sum_vec = np.asarray(sum_vec)
        Pomega_BY = sum_vec/np.sqrt(tensor_shape[0])
        return Pomega_BY
    elif dim == 2:
        Z = mat
        for i in range(m):
            #sum_vec += Z[c[i]][a[i]] * delta_vec[:, i]
            sum_vec.append(Z[spl[2][i]][spl[0][i]])
        sum_vec = np.asarray(sum_vec)
        Pomega_CZ = sum_vec/np.sqrt(tensor_shape[1])
        return Pomega_CZ

def adjoint_operator(spl, sample_vec, tensor_shape,m ,dim=None):# sample_vec[i]????  OR sample_vec[i][0]???
    # a, b, c = sample3D_rule(tensor_shape, m)
    if dim == 0:
        # sample_vec = Pomega_mat(a, b, c, mat, tensor_shape, m, 0)
        X = np.zeros((tensor_shape[0], tensor_shape[1]))
        for i in range(m):
            X[spl[0][i]][spl[1][i]] = sample_vec[i]
        return X
    elif dim == 1:
        # sample_vec = Pomega_mat(a, b, c, mat, tensor_shape, m, 1)
        Y = np.zeros((tensor_shape[1], tensor_shape[2]))
        for i in range(m):
            Y[spl[1][i]][spl[2][i]] = sample_vec[i]
        return Y
    elif dim == 2:
        # sample_vec = Pomega_mat(a, b, c, mat, tensor_shape, m, 2)
        Z = np.zeros((tensor_shape[2], tensor_shape[0]))
        for i in range(m):
            Z[spl[2][i]][spl[0][i]] = sample_vec[i]
        return Z
    else:
        print('error')
        return -1

def Pomega_tensor(spl, tensor, tensor_shape, sample_number):
    # a, b, c = sample3D_rule(tensor_shape, sample_number)
    # sample_t = np.zeros((sample_number, 1))
    sample_t = np.zeros(sample_number)
    for i in range(sample_number):
        sample_t[i] = tensor[spl[0][i]][spl[1][i]][spl[2][i]]
    # print(sample_t)
    return sample_t

def Pomega_Pair(spl, X, Y, Z, tensor_shape, m):
    Pomega_A = Pomega_mat(spl, X, tensor_shape, m, 0)
    Pomega_B = Pomega_mat(spl, Y, tensor_shape, m, 1)
    Pomega_C = Pomega_mat(spl, Z, tensor_shape, m, 2)
    Pomega_pair = Pomega_A + Pomega_B + Pomega_C
    # print(Pomega_pair)
    return Pomega_pair

def cone_projection_operator(x, t):
    norm_x = np.linalg.norm(x, ord=2)
    if norm_x <= t:
        return x, t
    if t <= np.negative(norm_x):
        return 0, 0
    if t >= np.negative(norm_x) and t <= norm_x:
        tmp = (norm_x + t)/(2 * norm_x)
        return tmp*x, tmp*norm_x

def SVT(mat, tao):
    u, s, v = np.linalg.svd(centralization(mat), full_matrices=True, compute_uv=True)
    for i in range(len(s)):
        s[i]=max(0, s[i]-tao)
    # print(u, s, v)
    return u, s, v

'''
def shrink(mat, mode='normal'):
        u,s,v = SVT(mat)
        if mode == 'normal':
            return np.matmul(np.matmul(u,s),v)
        if mode == 'complicated':
            delta = np.matmul(np.matmul(np.transpose(vecr1), mat),vecr2)/np.sqrt(row*col)
            tmp1 = np.matmul(np.matmul(u,s),v)
            tmp2 = (max(0,delta-tao)+min(0,delta+tao))*np.ones((mat.shape))/np.sqrt(row*col)
            return tmp1+tmp2
'''

def shrink(mat, tao, mode='normal'):
        u, s, v = SVT(mat, tao)
        sm = np.zeros((u.shape[0], v.shape[0]))
        for i in range(len(s)):
            sm[i][i] = s[i]
        if mode == 'normal':
            return np.matmul(np.matmul(u, sm), v)
        if mode == 'complicated':
            vecr1 = np.ones((mat.shape[0], 1))
            vecr2 = np.ones((mat.shape[1], 1))
            delta = np.matmul(np.matmul(np.transpose(vecr1), mat), vecr2)/np.sqrt(mat.shape[0]*mat.shape[1])
            tmp1 = np.matmul(np.matmul(u,sm), v)
            tmp2 = (max(0, delta-tao)+min(0, delta+tao))*np.ones(mat.shape)/np.sqrt(mat.shape[0]*mat.shape[1])
            return tmp1+tmp2

def shrinkageBorC(X_hat, tao, r):
    sum = 0
    s = r + 1
    U, S, V = np.linalg.svd(centralization(X_hat))
    '''
    while True:
        if (s + 5 < len(S)):
            s = s + 5
            if (S[s-5] <= tao):
                break
        else:
            if (S[s] <= tao):
                break
            s = s + 1
            if(s>=len(S)):

    for j in range(s-5, s):
        if(S[j] > tao):
            r = j
                 break
    '''
    for i in range(s, len(S)):
        if(S[i]<tao):
            r = i-1
            break

    for j in range(r):
        shape1 = U.shape
        shape2 = V.shape
        m = np.matmul(np.reshape(U[j,:],(shape1[0], 1)), np.reshape(np.transpose(V[j, :]), (1,shape2[0])))
        sum = sum +(S[j]-tao)*m
    X = sum
    return X, r

def shrinkageA(X_hat, tao, r):
    X, r = shrinkageBorC(X_hat, tao, r)
    # delta = np.inner(np.ones(X_hat.shape), X_hat) # elementwise sum of X_hat
    delta = np.sum(X_hat)
    gamma = (max(0, delta-tao)+min(0, delta+tao))/(X_hat.shape[0]*X_hat.shape[1])
    ones_vec1 = np.reshape(np.ones(X_hat.shape[0]), (X_hat.shape[0], 1))
    ones_vec2 = np.reshape(np.ones(X_hat.shape[1]), (1, X_hat.shape[1]))
    result = X + gamma*np.matmul(ones_vec1, ones_vec2)

    return result, r


text:
from models.duplo import DuPLO
import torch
from train import prepare_dataset
from argparse import Namespace
from tqdm import tqdm
import numpy as np
import sklearn.metrics
import pandas as pd
import os
import argparse

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        'mode', type=str, default="12classes", help='either 12classes or 23classes')
    parser.add_argument(
        'dataset', type=str, default="BavarianCrops", help='either BavarianCrops or GAFv2')
    args, _ = parser.parse_known_args()

    return args

def main(args):

    traindataloader, testdataloader, model, args, device = setup(args.dataset, args.mode)

    optimizer = torch.optim.Adam(
                filter(lambda x: x.requires_grad, model.parameters()),
                betas=(0.9, 0.999), eps=1e-08, weight_decay=args.weight_decay, lr=args.learning_rate)

    stats = list()
    for epoch in range(args.epochs):
        trainlosses = train_epoch(traindataloader, optimizer, model, device)
        testlosses, predictions, labels = test_epoch(testdataloader, model, device)

        trainloss = np.array(trainlosses).mean()
        testloss = np.array(testlosses).mean()
        stat = metrics(labels, predictions)
        stat["epoch"] = epoch
        stat["trainloss"] = trainloss
        stat["testloss"] = testloss
        stats.append(stat)
        print(f"epoch {epoch}: trainloss: {trainloss:.2f}, testloss: {testloss:.2f}, kappa: {stat['kappa']:.2f}, accuracy: {stat['accuracy']:.2f}")

        stat_df = pd.DataFrame(stats).set_index("epoch")

        os.makedirs(os.path.join(args.store, args.experiment), exist_ok=True)
        stat_df.to_csv(os.path.join(args.store,args.experiment, "log.csv"))

        better = not (stat_df["testloss"] < stat["testloss"]).any()
        if better:
            model.save(os.path.join(args.store,args.experiment, "model.pth"))

def setup(dataset, mode, dataroot="../data", store='/tmp/'):
    if mode == "12classes":
        classmapping = os.path.join(dataroot, "BavarianCrops", 'classmapping12.csv')
    elif mode == "23classes":
        classmapping = os.path.join(dataroot, "BavarianCrops", 'classmapping23.csv')

    args = Namespace(batchsize=256,epochs=1500,
                     classmapping=classmapping,
                     dataroot=dataroot, dataset=dataset,
                     model='duplo', mode=None,weight_decay=0, learning_rate=1e-3,
                     seed=0, store=store, workers=0)

    if dataset == "BavarianCrops":
        args = merge([args, TUM_dataset])
        exp = "isprs_tum_duplo"
    elif dataset == "GAFv2":
        args = merge([args, GAF_dataset])
        exp = "isprs_gaf_duplo"
    args.experiment = exp
    args.store = f"/tmp/{mode}"

    args.train_on = "train"
    args.test_on = "valid"
    traindataloader, testdataloader = prepare_dataset(args)

    input_dim = traindataloader.dataset.datasets[0].ndims
    nclasses = len(traindataloader.dataset.datasets[0].classes)

    device = torch.device("cuda")
    model = DuPLO(input_dim=input_dim, nclasses=nclasses, sequencelength=args.samplet, dropout=0.4)

    model.to(device)

    return traindataloader, testdataloader, model, args, device

def train_epoch(dataloader, optimizer, model, device):
    losses = list()
    model.train()
    for iteration, data in tqdm(enumerate(dataloader), total=len(dataloader), leave=False):
        optimizer.zero_grad()

        inputs, targets, _ = data

        inputs = inputs.to(device)
        targets = targets.to(device)

        logprobabilities, logprobabilities_cnn, logprobabilities_rnn = model.forward(inputs.transpose(1,2))

        loss = torch.nn.functional.nll_loss(logprobabilities, targets[:, 0])
        loss_cnn = torch.nn.functional.nll_loss(logprobabilities_cnn, targets[:, 0])
        loss_rnn = torch.nn.functional.nll_loss(logprobabilities_rnn, targets[:, 0])

        loss += 0.3 * loss_cnn
        loss += 0.3 * loss_rnn

        loss.backward()
        optimizer.step()

        losses.append(loss.detach().cpu().numpy())
        return losses

def test_epoch(dataloader, model, device):
    predictions = list()
    losses = list()
    labels = list()
    model.eval()

    with torch.no_grad():
        for iteration, data in tqdm(enumerate(dataloader), total=len(dataloader), leave=False):

            inputs, targets, _ = data

            inputs = inputs.to(device)
            targets = targets.to(device)

            logprobabilities, logprobabilities_cnn, logprobabilities_rnn = model.forward(inputs.transpose(1, 2))

            loss = torch.nn.functional.nll_loss(logprobabilities, targets[:, 0])
            loss_cnn = torch.nn.functional.nll_loss(logprobabilities_cnn, targets[:, 0])
            loss_rnn = torch.nn.functional.nll_loss(logprobabilities_rnn, targets[:, 0])

            loss += 0.3 * loss_cnn
            loss += 0.3 * loss_rnn

            losses.append(loss.detach().cpu().numpy())
            y_pred = logprobabilities.argmax(-1)

            predictions.append(y_pred.detach().cpu().numpy())
            labels.append(targets[:, 0].detach().cpu().numpy())

    return losses, np.hstack(predictions), np.hstack(labels)

def metrics(y_true, y_pred):
    accuracy = sklearn.metrics.accuracy_score(y_true, y_pred)
    kappa = sklearn.metrics.cohen_kappa_score(y_true, y_pred)
    f1_micro = sklearn.metrics.f1_score(y_true, y_pred, average="micro")
    f1_macro = sklearn.metrics.f1_score(y_true, y_pred, average="macro")
    f1_weighted = sklearn.metrics.f1_score(y_true, y_pred, average="weighted")
    recall_micro = sklearn.metrics.recall_score(y_true, y_pred, average="micro")
    recall_macro = sklearn.metrics.recall_score(y_true, y_pred, average="macro")
    recall_weighted = sklearn.metrics.recall_score(y_true, y_pred, average="weighted")
    precision_micro = sklearn.metrics.precision_score(y_true, y_pred, average="micro")
    precision_macro = sklearn.metrics.precision_score(y_true, y_pred, average="macro")
    precision_weighted = sklearn.metrics.precision_score(y_true, y_pred, average="weighted")

    return dict(
        accuracy=accuracy,
        kappa=kappa,
        f1_micro=f1_micro,
        f1_macro=f1_macro,
        f1_weighted=f1_weighted,
        recall_micro=recall_micro,
        recall_macro=recall_macro,
        recall_weighted=recall_weighted,
        precision_micro=precision_micro,
        precision_macro=precision_macro,
        precision_weighted=precision_weighted,
    )

def merge(namespaces):
    merged = dict()

    for n in namespaces:
        d = n.__dict__
        for k,v in d.items():
            merged[k]=v

    return Namespace(**merged)

TUM_dataset = Namespace(
    dataset = "BavarianCrops",
    trainregions = ["holl","nowa","krum"],
    testregions = ["holl","nowa","krum"],
    scheme="blocks",
    test_on = "test",
    train_on = "trainvalid",
    samplet = 70
)

GAF_dataset = Namespace(
    dataset = "GAFv2",
    trainregions = ["holl","nowa","krum"],
    testregions = ["holl","nowa","krum"],
    features = "optical",
    scheme="blocks",
    test_on="test",
    train_on="train",
    samplet = 23
)

if __name__ == '__main__':
    args = parse_args()
    main(args)


text:
from __future__ import absolute_import, division, print_function

import os
import warnings
import numpy as np
from multiprocessing import cpu_count

from astropy.wcs import WCS
from astropy import units as u
from astropy.nddata import StdDevUncertainty
from astropy.utils.console import ProgressBar

from .nikamap import NikaMap

__all__ = ["HalfDifference", "Jackknife", "Bootstrap"]

def compare_header(header_ref, header_target):
    """Crude comparison of two header

    Parameters
    ----------
    header_ref : astropy.io.fits.Header
        the reference header
    header_target : astropy.io.fits.Header
        the target header to check

    Notes
    -----
    This will raise assertion error if the two header are not equivalent
    """
    wcs_ref = WCS(header_ref)
    wcs_target = WCS(header_target)

    assert wcs_ref.wcs == wcs_target.wcs, "Different header found"
    for key in ["UNIT", "NAXIS1", "NAXIS2"]:
        if key in header_ref:
            assert header_ref[key] == header_target[key], "Different key found"

def check_filenames(filenames):
    """check filenames existence

    Parameters
    ----------
    filenames : list of str
        filenames list to be checked

    Returns
    -------
    list of str
        curated list of files
    """
    _filenames = []
    for filename in filenames:
        if os.path.isfile(filename):
            _filenames.append(filename)
        else:
            warnings.warn("{} does not exist, removing from list".format(filename), UserWarning)
    return _filenames

class MultiScans(object):
    """A class to hold multi single scans from a list of fits files.

    This acts as a python lazy iterator and/or a callable

    Parameters
    ----------
    filenames : list or `~MultiScans` object
        the list of fits files to produce the Jackknifes or an already filled object
    ipython_widget : bool, optional
        If True, the progress bar will display as an IPython notebook widget.
    ignore_header : bool, optional
        if True, the check on header is ignored
    n : int
        the number of iteration for the iterator

    Notes
    -----
    A crude check is made on the wcs of each map when instanciated
    """

    def __init__(self, filenames, n=None, ipython_widget=False, ignore_header=False, **kwd):

        self.i = 0
        self.n = n
        self.kwargs = kwd
        self.ipython_widget = ipython_widget

        if isinstance(filenames, MultiScans):
            data = filenames

            self.filenames = data.filenames
            self.primary_header = data.primary_header
            self.header = data.header
            self.unit = data.unit
            self.shape = data.shape
            self.datas = data.datas
            self.weights = data.weights
            self.time = data.time
            self.mask = data.mask
        else:
            self.filenames = check_filenames(filenames)

            nm = NikaMap.read(self.filenames[0], **kwd)

            self.primary_header = nm.meta.get("primary_header", None)
            self.header = nm.meta["header"]
            self.unit = nm.unit
            self.shape = nm.shape

            # This is a low_mem=False case ...
            # TODO: How to refactor that for low_mem=True ?
            datas = np.zeros((len(self.filenames),) + self.shape)
            weights = np.zeros((len(self.filenames),) + self.shape)
            time = np.zeros(self.shape) * u.h

            for i, filename in enumerate(ProgressBar(self.filenames, ipython_widget=self.ipython_widget)):

                nm = NikaMap.read(filename, **kwd)
                try:
                    compare_header(self.header, nm.meta["header"])
                except AssertionError as e:
                    if ignore_header:
                        warnings.warn("{} for {}".format(e, filename), UserWarning)
                    else:
                        raise ValueError("{} for {}".format(e, filename))

                datas[i, :, :] = nm.data
                with np.errstate(invalid="ignore", divide="ignore"):
                    weights[i, :, :] = nm.uncertainty.array ** -2
                time += nm.time

                # make sure that we do not have nans in the data
                datas[i, nm.time == 0] = 0
                weights[i, nm.time == 0] = 0

            unobserved = time == 0

            self.datas = datas
            self.weights = weights
            self.time = time
            self.mask = unobserved

    def __len__(self):
        # to retrieve the legnth of the iterator, enable ProgressBar on it
        return self.n

    def __iter__(self):
        # Iterators are iterables too.
        # Adding this functions to make them so.
        return self

    def __call__(self):
        """The main method which should be overrided

        should return a  :class:`nikamap.NikaMap`
        """
        pass

    def __next__(self):
        """Iterator on the objects"""
        if self.n is None or self.i < self.n:
            # Produce data until last iter
            self.i += 1
            data = self.__call__()
        else:
            raise StopIteration()

        return data

class HalfDifference(MultiScans):
    """A class to create weighted half differences uncertainty maps from a list of scans.

    This acts as a python lazy iterator and/or a callable

    Parameters
    ----------
    filenames : list
        the list of fits files to produce the Jackknifes
    ipython_widget : bool, optional
        If True, the progress bar will display as an IPython notebook widget.
    n : int
        the number of Jackknifes maps to be produced in the iterator

            if set to `None`, produce only one weighted average of the maps

    parity_threshold : float
        mask threshold between 0 and 1 to keep partially jackknifed area
        * 1 pure jackknifed
        * 0 partially jackknifed, keep all

    Notes
    -----
    A crude check is made on the wcs of each map when instanciated
    """

    def __init__(self, filenames, parity_threshold=1, **kwd):

        super(HalfDifference, self).__init__(filenames, **kwd)
        self.parity_threshold = parity_threshold

        # Create weights for Half differences
        jk_weights = np.ones(len(self.filenames))

        if self.n is not None:
            jk_weights[::2] *= -1

        if self.n is not None and len(self.filenames) % 2:
            warnings.warn("Even number of files, dropping a random file", UserWarning)
            jk_weights[-1] = 0

        assert np.sum(jk_weights != 0), "Less than 2 existing files in filenames"

        self.jk_weights = jk_weights

    @property
    def parity_threshold(self):
        return self._parity

    @parity_threshold.setter
    def parity_threshold(self, value):
        if value is not None and isinstance(value, (int, float)) and 0 <= value <= 1:
            self._parity = value
        else:
            raise TypeError("parity must be between 0 and 1")

    def __call__(self):
        """Compute a jackknifed dataset

        Returns
        -------
        :class:`nikamap.NikaMap`
            a jackknifed data set
        """
        np.random.shuffle(self.jk_weights)

        with np.errstate(invalid="ignore", divide="ignore"):
            e_data = 1 / np.sqrt(np.sum(self.weights, axis=0))
            data = np.sum(self.datas * self.weights * self.jk_weights[:, np.newaxis, np.newaxis], axis=0) * e_data ** 2
            parity = np.mean((self.weights != 0) * self.jk_weights[:, np.newaxis, np.newaxis], axis=0)
            # TBC: In principle we should use a weighted parity to avoid different scans/weights problems
            # weighted_parity = np.sum(self.weights * self.jk_weights[:, np.newaxis, np.newaxis], axis=0) * e_data ** 2

        if self.n is not None:
            mask = (1 - np.abs(parity)) < self.parity_threshold
        else:
            mask = parity < self.parity_threshold

        mask = mask | self.mask

        data[mask] = np.nan
        e_data[mask] = np.nan

        # TBC: time will have a different mask here....
        data = NikaMap(
            data,
            mask=mask,
            uncertainty=StdDevUncertainty(e_data),
            unit=self.unit,
            wcs=WCS(self.header),
            meta={"header": self.header, "primary_header": self.primary_header},
            time=self.time,
        )

        return data  # , weighted_parity

class Jackknife(MultiScans):
    """A class to create weighted Jackknife maps from a list of scans.

    This acts as a python lazy iterator and/or a callable

    Parameters
    ----------
    filenames : list
        the list of fits files to produce the Jackknifes
    n_samples : int
        The number of (sub) samples to use (from 2 to len(filenames))
    parity_threshold : float
        mask threshold between 0 and 1 to keep partially jackknifed area
        * 1 pure jackknifed
        * 0 partially jackknifed, keep all
    ipython_widget : bool, optional
        If True, the progress bar will display as an IPython notebook widget.
    n : int
        the number of Jackknifes maps to be produced by the iterator

    Notes
    -----
    A crude check is made on the wcs of each map when instanciated
    """

    def __init__(self, filenames, n_samples=None, parity_threshold=1, **kwd):

        super(Jackknife, self).__init__(filenames, **kwd)

        assert len(self.filenames) > 1, "Less than 2 existing files in filenames"

        self.n_samples = n_samples  # Will create the indexes for the sub-samples
        self.parity_threshold = parity_threshold

    @property
    def parity_threshold(self):
        return self._parity

    @parity_threshold.setter
    def parity_threshold(self, value):
        if value is not None and isinstance(value, (int, float)) and 0 <= value <= 1:
            self._parity = value
        else:
            raise TypeError("parity must be between 0 and 1")

    @property
    def n_samples(self):
        return self._n_samples

    @n_samples.setter
    def n_samples(self, value):
        if value is None:
            value = len(self.filenames)

        assert (2 <= value) and (value <= len(self.filenames)), "n_samples must be between 2 and the number of scans"

        self._n_samples = value

        # Check compatibility between n_samples and filenames length
        n_filenames = len(self.filenames)
        remainder = n_filenames % value

        if remainder:
            warnings.warn("Remainder in number of files for {} samples, dropping the last {}".format(value, remainder), UserWarning)
            n_filenames -= remainder

        assert n_filenames, "Less than 2 existing files in filenames"

        # Create the indexes for the sub-samples
        indexes = np.repeat(np.arange(value), n_filenames // value)

        if remainder:
            indexes = np.concatenate([indexes, np.full(remainder, np.nan)])

        self.indexes = indexes

    def __call__(self):
        """Compute a jackknifed dataset

        Returns
        -------
        :class:`nikamap.NikaMap`
            a jackknifed data set
        """
        np.random.shuffle(self.indexes)

        with np.errstate(invalid="ignore", divide="ignore"):

            # Compute sub-samples
            sub_datas = []
            sub_weights = []
            for idx in range(self.n_samples):
                mask = self.indexes == idx
                data, weight = np.ma.average(self.datas[mask], weights=self.weights[mask], axis=0, returned=True)
                sub_datas.append(data)
                sub_weights.append(weight)

            sub_datas = np.ma.array(sub_datas)
            sub_weights = np.ma.array(sub_weights)

            data = np.ma.average(sub_datas, weights=sub_weights, axis=0)
            # unweighted sample variance
            V1 = self.n_samples
            e_data = np.sqrt(np.sum((sub_datas - data)**2, axis=0) / (V1 * (V1 - 1)))
            # TODO : weighted sample variance (NOT WORKING !!!)
            # V1 = np.sum(sub_weights, axis=0)
            # V2 = np.sum(sub_weights**2, axis=0)
            # e_data = np.sqrt(np.sum(sub_weights * (sub_datas - data)**2, axis=0)  / (V1 - V2 / V1) )
            # e_data = e_data.filled(np.nan)

            parity = np.mean(sub_weights != 0, axis=0)

            # TBC: In principle we should use a weighted parity to avoid different scans/weights problems

            mask = parity < self.parity_threshold

        mask = mask | self.mask

        data[mask] = np.nan
        e_data[mask] = np.nan

        # TBC: time will have a different mask here....
        data = NikaMap(
            data,
            mask=mask,
            uncertainty=StdDevUncertainty(e_data),
            unit=self.unit,
            wcs=WCS(self.header),
            meta={"header": self.header, "primary_header": self.primary_header},
            time=self.time,
        )

        return data  # , weighted_parity

class Bootstrap(MultiScans):
    """A class to create bootstraped maps from a list of IDL fits files.

    This acts as a python lazy iterator and/or a callable

    Parameters
    ----------
    filenames : list
        the list of fits files to produce the Jackknifes
    n_bootstrap : int
        the number of realization to produce a bootsrapped map, by default 20 times the length of the input filename list
    ipython_widget : bool, optional
        If True, the progress bar will display as an IPython notebook widget.
    n : int
        the number of bootstrap maps to be produced by the iterator

    Notes
    -----
    A crude check is made on the wcs of each map when instanciated
    """

    def __init__(self, filenames, n_bootstrap=None, **kwd):
        super(Bootstrap, self).__init__(filenames, **kwd)

        if n_bootstrap is None:
            n_bootstrap = 20 * len(self.filenames)

        self.n_bootstrap = n_bootstrap

    def shuffled_average(self, *args):
        """Actually do one shuffled average."""
        if len(args) > 0:
            n_shuffle = len(args[0])
        else:
            n_shuffle = 1
        n_scans = self.datas.shape[0]
        outputs = []
        for _ in range(n_shuffle):
            shuffled_index = np.floor(np.random.uniform(0, n_scans, n_scans)).astype(np.int)
            # np.ma.average is needed as some of the pixels have zero weights (should be masked)
            outputs.append(np.ma.average(self.datas[shuffled_index], weights=self.weights[shuffled_index], axis=0, returned=False))
        return outputs

    def __call__(self):
        """Compute a bootstraped map

        Returns
        -------
        :class:`nikamap.NikaMap`
            a bootstraped data set
        """

        bs_array = np.concatenate(
            ProgressBar.map(self.shuffled_average, np.array_split(np.arange(self.n_bootstrap), cpu_count()), ipython_widget=self.ipython_widget, multiprocess=True,)
        )

        data = np.mean(bs_array, axis=0)
        e_data = np.std(bs_array, axis=0, ddof=1)

        # Mask unobserved regions
        unobserved = self.time == 0
        data[unobserved] = np.nan
        e_data[unobserved] = np.nan

        data = NikaMap(
            data,
            mask=unobserved,
            uncertainty=StdDevUncertainty(e_data),
            unit=self.unit,
            wcs=WCS(self.header),
            meta={"header": self.header, "primary_header": self.primary_header},
            time=self.time,
        )

        return data

max_stars_repo_name: abeelen/nikamap

max_issues_repo_name: abeelen/nikamap

max_forks_repo_name: abeelen/nikamap

reason: import numpy,from astropy


text:
# Copyright 2021 The Deluca 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.

"""run controller."""
import datetime
import functools
import os
import pickle

from deluca.lung.controllers._expiratory import Expiratory
from deluca.lung.core import BreathWaveform
from deluca.lung.envs._balloon_lung import BalloonLung
import jax
import jax.numpy as jnp
import tqdm

# pylint: disable=invalid-name
# pylint: disable=unused-argument

def run_controller(
    controller,
    T=1000,
    dt=0.03,
    abort=60,
    env=None,
    waveform=None,
    use_tqdm=False,
    directory=None,
):
  """run controller over horizon of length T."""
  env = env or BalloonLung()
  waveform = waveform or BreathWaveform.create()
  expiratory = Expiratory.create(waveform=waveform)

  result = {}

  controller_state = controller.init()
  expiratory_state = expiratory.init()

  tt = range(T)
  if use_tqdm:
    tt = tqdm.tqdm(tt, leave=False)

  timestamps = jnp.zeros(T)
  pressures = jnp.zeros(T)
  flows = jnp.zeros(T)
  u_ins = jnp.zeros(T)
  u_outs = jnp.zeros(T)

  state, obs = env.reset()

  try:
    for i, _ in enumerate(tt):
      pressure = obs.predicted_pressure
      if env.should_abort():
        break

      controller_state, u_in = controller.__call__(controller_state, obs)
      expiratory_state, u_out = expiratory.__call__(expiratory_state, obs)

      u_in = u_in.squeeze()
      state, obs = env(state, (u_in, u_out))

      timestamps = timestamps.at[i].set(env.time(state) - dt)
      u_ins = u_ins.at[i].set(u_in)
      u_outs = u_outs.at[i].set(u_out)
      pressures = pressures.at[i].set(pressure)
      flows = flows.at[i].set(env.flow)

      env.wait(max(dt - env.dt, 0))

  finally:
    env.cleanup()

  timeseries = {
      "timestamp": jnp.array(timestamps),
      "pressure": jnp.array(pressures),
      "flow": jnp.array(flows),
      "target": waveform.at(timestamps),
      "u_in": jnp.array(u_ins),
      "u_out": jnp.array(u_outs),
  }

  for key, val in timeseries.items():
    timeseries[key] = val[:T + 1]

  result["timeseries"] = timeseries
  result["waveform"] = waveform

  if directory is not None:
    if not os.path.exists(directory):
      os.makedirs(directory)
    timestamp = datetime.datetime.now().strftime("%Y-%m-%dT%H-%M-%S")
    pickle.dump(result, open(f"{directory}/{timestamp}.pkl", "wb"))

  return result

def loop_over_tt(envState_obs_ctrlState_ExpState_i, dummy_data, controller,
                 expiratory, env, dt):
  """loop function over time points."""
  state, obs, controller_state, expiratory_state, i = envState_obs_ctrlState_ExpState_i
  pressure = obs.predicted_pressure

  # if env.should_abort(): # TODO(dsuo): how to handle break in scan
  #     break

  controller_state, u_in = controller.__call__(controller_state, obs)
  expiratory_state, u_out = expiratory.__call__(expiratory_state, obs)

  state, obs = env(state, (u_in, u_out))

  timestamps_i = env.time(state) - dt
  u_ins_i = u_in
  u_outs_i = u_out
  pressures_i = pressure
  flows_i = env.flow

  env.wait(max(dt - env.dt, 0))
  return (state, obs, controller_state, expiratory_state,
          i + 1), (timestamps_i, u_ins_i, u_outs_i, pressures_i, flows_i)

def run_controller_scan(
    controller,
    T=1000,
    dt=0.03,
    abort=60,
    env=None,
    waveform=None,
    use_tqdm=False,
    directory=None,
    init_controller=False,
):
  """run controller scan version."""
  env = env or BalloonLung()
  waveform = waveform or BreathWaveform.create()
  expiratory = Expiratory.create(waveform=waveform)

  result = {}

  if init_controller:
    controller_state = controller.init(waveform)
  else:
    controller_state = controller.init()
  expiratory_state = expiratory.init()

  tt = range(T)
  if use_tqdm:
    tt = tqdm.tqdm(tt, leave=False)

  timestamps = jnp.zeros(T)
  pressures = jnp.zeros(T)
  flows = jnp.zeros(T)
  u_ins = jnp.zeros(T)
  u_outs = jnp.zeros(T)

  state, obs = env.reset()
  # xp = jnp.array(waveform.xp)
  # fp = jnp.array(waveform.fp)
  # period = waveform.period
  # dtype = waveform.dtype

  jit_loop_over_tt = jax.jit(functools.partial(loop_over_tt,
                                               controller=controller,
                                               expiratory=expiratory,
                                               env=env,
                                               dt=dt))
  try:
    _, (timestamps, u_ins, u_outs, pressures, flows) = jax.lax.scan(
        jit_loop_over_tt, (state, obs, controller_state, expiratory_state, 0),
        jnp.arange(T))

  finally:
    env.cleanup()

  timeseries = {
      "timestamp": jnp.array(timestamps),
      "pressure": jnp.array(pressures),
      "flow": jnp.array(flows),
      "target": waveform.at(timestamps),
      "u_in": jnp.array(u_ins),
      "u_out": jnp.array(u_outs),
  }

  for key, val in timeseries.items():
    timeseries[key] = val[:T + 1]

  result["timeseries"] = timeseries
  result["waveform"] = waveform

  if directory is not None:
    if not os.path.exists(directory):
      os.makedirs(directory)
    timestamp = datetime.datetime.now().strftime("%Y-%m-%dT%H-%M-%S")
    pickle.dump(result, open(f"{directory}/{timestamp}.pkl", "wb"))

  return result

max_stars_repo_name: google/deluca

max_issues_repo_name: google/deluca

max_forks_repo_name: google/deluca

reason: import jax


text:
# Copyright 2021-present Kensho Technologies, LLC.
"""Create a CSV that contains pre-wikitext parsing article metadata."""
import logging
import os

import networkx as nx
import pandas as pd

from kwnlp_preprocessor import argconfig

logger = logging.getLogger(__name__)

# TODO: centralize this in argconfig
# have to change in task 39p1 create article too
ROOT_NQIDS = [
    17442446,  # Wikimedia internal item
    14795564,  # point in time with respect to recurrent timeframe
    18340514,  # events in a specific year or time period
    5,  # human
    2221906,  # geographic location
    43229,  # organization
    4830453,  # business
]

def main(
    wp_yyyymmdd: str,
    wd_yyyymmdd: str,
    data_path: str = argconfig.DEFAULT_KWNLP_DATA_PATH,
    wiki: str = argconfig.DEFAULT_KWNLP_WIKI,
) -> None:

    wp_dump_path = os.path.join(data_path, f"wikipedia-derived-{wp_yyyymmdd}")
    wd_dump_path = os.path.join(data_path, f"wikidata-derived-{wd_yyyymmdd}")

    # read title mapper CSV as base
    # ====================================================================
    file_path = os.path.join(
        wp_dump_path,
        "kwnlp-sql",
        f"kwnlp-{wiki}-{wp_yyyymmdd}-title-mapper.csv",
    )
    logger.info(f"reading {file_path}")
    df = pd.read_csv(
        file_path,
        keep_default_na=False,
        usecols=["source_id", "source_title", "is_redirect"],
    )

    # get base information from title mapper
    # ====================================================================
    df = df[~df["is_redirect"]]
    df = df[["source_id", "source_title"]]
    df = df.rename(columns={"source_id": "page_id", "source_title": "page_title"})

    # read item id from page props CSV
    # ====================================================================
    file_path = os.path.join(
        wp_dump_path,
        "kwnlp-sql",
        f"kwnlp-{wiki}-{wp_yyyymmdd}-page-props.csv",
    )
    logger.info(f"reading {file_path}")
    df_pp = pd.read_csv(file_path, usecols=["page_id", "wikibase_item"])

    # add wikidata item id
    # ====================================================================
    df = pd.merge(df, df_pp, on="page_id", how="left")

    df["wikibase_item"] = df["wikibase_item"].fillna("Q-1").apply(lambda x: int(x[1:]))

    df = df.rename(columns={"wikibase_item": "item_id"})

    # read subclass of info
    # ====================================================================
    file_path = os.path.join(
        wd_dump_path,
        "p279-claim",
        f"kwnlp-wikidata-{wd_yyyymmdd}-p279-claim.csv",
    )
    logger.info(f"reading {file_path}")
    df_p279 = pd.read_csv(file_path, usecols=["source_id", "target_id"])
    logger.info("building p279 graph")
    g_p279 = nx.DiGraph()
    g_p279.add_edges_from(df_p279.values)

    # read instance of info
    # ====================================================================
    file_path = os.path.join(
        wd_dump_path,
        "p31-claim",
        f"kwnlp-wikidata-{wd_yyyymmdd}-p31-claim.csv",
    )
    logger.info(f"reading {file_path}")
    df_p31 = pd.read_csv(file_path, usecols=["source_id", "target_id"])

    # add root qid tags
    # ====================================================================
    for root_nqid in ROOT_NQIDS:

        logger.info(f"tagging items in subclass tree of {root_nqid}")
        col_name = f"isa_Q{root_nqid}"
        subclass_qids = set(nx.ancestors(g_p279, root_nqid)).union(set([root_nqid]))
        mask = df_p31["target_id"].isin(subclass_qids)
        tmp_df = df_p31[mask].copy()
        tmp_df = tmp_df.sort_values(by=["source_id", "target_id"])
        tmp_df = tmp_df.drop_duplicates(subset=["source_id"], keep="first")
        tmp_df = tmp_df.rename(columns={"source_id": "item_id", "target_id": col_name})

        df = pd.merge(df, tmp_df, on="item_id", how="left")
        df[col_name] = df[col_name].fillna(False).astype(int)

    # read views CSV
    # ====================================================================
    file_path = os.path.join(
        wp_dump_path,
        "kwnlp-sql",
        f"kwnlp-{wiki}-{wp_yyyymmdd}-prior-month-pageviews-complete.csv",
    )
    logger.info(f"reading {file_path}")
    df_views = pd.read_csv(file_path, keep_default_na=False)

    # add views
    # ====================================================================
    df = pd.merge(df, df_views, on="page_title", how="left")
    df["views"] = df["views"].fillna(0).astype("int")

    # sort and write output
    # ====================================================================
    df = df[["page_id", "item_id", "page_title", "views"] + [f"isa_Q{nqid}" for nqid in ROOT_NQIDS]]
    df = df.sort_values("page_id")
    file_path = os.path.join(
        wp_dump_path,
        "kwnlp-sql",
        f"kwnlp-{wiki}-{wp_yyyymmdd}-article-pre.csv",
    )
    os.makedirs(os.path.dirname(file_path), exist_ok=True)
    logger.info(f"writing {file_path}")
    df.to_csv(file_path, index=False)
    return df

if __name__ == "__main__":

    description = "create kwnlp article pre"
    arg_names = ["wp_yyyymmdd", "wd_yyyymmdd", "data_path", "wiki", "loglevel"]
    parser = argconfig.get_argparser(description, arg_names)

    args = parser.parse_args()
    logging.basicConfig(level=args.loglevel)
    logger.info(f"args={args}")

    main(
        args.wp_yyyymmdd,
        args.wd_yyyymmdd,
        data_path=args.data_path,
        wiki=args.wiki,
    )

max_stars_repo_name: kensho-technologies/kwnlp-preprocessor

max_issues_repo_name: kensho-technologies/kwnlp-preprocessor

max_forks_repo_name: kensho-technologies/kwnlp-preprocessor

reason: import networkx


text:
"""pytest initialization."""
import numpy as np
import scipy
from distutils.version import LooseVersion
import warnings

""" Set numpy print options to "legacy" for new versions of numpy
 If imported into a file, pytest will run this before any doctests.
 References
-----------
https://github.com/numpy/numpy/commit/710e0327687b9f7653e5ac02d222ba62c657a718
https://github.com/numpy/numpy/commit/734b907fc2f7af6e40ec989ca49ee6d87e21c495
https://github.com/nipy/nibabel/pull/556
"""
if LooseVersion(np.__version__) >= LooseVersion('1.14'):
    np.set_printoptions(legacy='1.13')

warnings.simplefilter(action="default", category=FutureWarning)
warnings.simplefilter("always", category=UserWarning)
# List of files that pytest should ignore
collect_ignore = ["testing/decorators.py", ]

max_stars_repo_name: Raniac/NEURO-LEARN

max_forks_repo_name: Raniac/NEURO-LEARN

reason: import numpy,import scipy


text:
"""self_attention
PyTorch Specialized Self Attention Modules
"""

import torch
import torch.nn as nn
import torch.nn.functional as F
import numpy as np

######################################################################
# Base Class for Self Attention Mechanisms
# =========================================
#

class BaseSelfAttention(nn.Module):

    def __init__(self):
        super(BaseSelfAttention, self).__init__()

    def init_linear(self, input_linear):
        """Initialize linear transformation"""
        bias = np.sqrt(6.0 / (input_linear.weight.size(0) + \
            input_linear.weight.size(1)))
        nn.init.uniform_(input_linear.weight, -bias, bias)
        if input_linear.bias is not None:
            input_linear.bias.data.zero_()

    def initialize_layers(self):
        raise NotImplementedError

    def forward(self, X):
        raise NotImplementedError

    def score(self, a, b):
        raise NotImplementedError

######################################################################
# Structured Self-Attentive Sentence Embedding
# =============================================
#
# Paper: https://arxiv.org/pdf/1703.03130.pdf
#
# Given a sentence of words, turn the sequence of word embeddings into
# a single sentence embedding that attends to important part of the
# sentence. The number of output sentence embeddings is determined by
# the number of annotations specified at module initialization.
#

class SentenceEmbedding(BaseSelfAttention):

    def __init__(self, embedding_dim, hidden_dim, num_annotations):
        super(SentenceEmbedding, self).__init__()

        # Save Parameters for reference
        self.embedding_dim = embedding_dim
        self.hidden_dim = hidden_dim
        self.num_annotations = num_annotations

        # Define layers
        self.initialize_layers()

    def initialize_layers(self):
        self.Ws1 = nn.Linear(self.embedding_dim, self.hidden_dim)
        self.Ws2 = nn.Linear(self.hidden_dim, self.num_annotations)
        # Initialize Layers
        self.init_linear(self.Ws1)
        self.init_linear(self.Ws2)

    def forward(self, word_embeddings):
        """
        Args:
            word_embeddings:
                (batch_size, doc_maxlen, embedding_dim)
        Output:
            sentence_embedding:
                (batch_size, num_annotations, embedding_dim)
        """
        # dim: (batch_size, doc_maxlen, embedding_dim)
        hidden = F.tanh(self.Ws1(word_embeddings))
        # dim: (batch_size, doc_maxlen, hidden_dim)
        atten_weights = F.softmax(self.Ws2(hidden), dim=2)
        # dim: (batch_size, doc_maxlen, num_annotations)
        atten_weights = atten_weights.transpose(1, 2)
        # dim: (batch_size, num_annotations, doc_maxlen)
        sentence_embedding = atten_weights.bmm(word_embeddings)
        # dim: (batch_size, num_annotations, embedding_dim)
        return sentence_embedding

######################################################################
# Self Attention using Luong Attention
# =====================================
#
# Paper: https://arxiv.org/pdf/1508.04025.pdf
#
# Luong Attention Scoring Function
#
#   score(h_t, h_s_bar)
#      = h_t^T * h_s_bar (dot)
#      = h_t^T * W_a * h_s_bar (general)
#      = v_a^T * W_a * [h_t; h_s_bar] (concat)
#
# This does not implement the local attention proposed by the original paper.
#
# ADVICE
# =======
# It is not recommended to use `concat` scoring for self attention
# when document length is too long (>20-100), because it will VERY VERY SLOW and
# memory intensive, because of the amount of individual scoring we need to do
#

class SelfAttention(BaseSelfAttention):

    def __init__(self, hidden_dim, scoring="general"):
        super(SelfAttention, self).__init__()

        # Save parameters for reference
        self.scoring = scoring
        self.hidden_dim = hidden_dim

        ### Define Layers ###
        self.initialize_layers()

    def initialize_layers(self):
        if self.scoring == 'general':
            self.W = nn.Linear(self.hidden_dim, self.hidden_dim)
            # initialization
            self.init_linear(self.W)
        elif self.scoring == 'concat':
            self.W = nn.Linear(2*self.hidden_dim, self.hidden_dim)
            self.v = nn.Linear(self.hidden_dim, 1)
            # initialization
            self.init_linear(self.W)
            self.init_linear(self.v)
        elif self.scoring == "dot":
            pass # don't need to do anything
        else:
            raise RuntimeError("Unrecognized attention scoring method: %s" % self.scoring)

    def forward(self, hidden_outputs):
        # hidden_outputs: (batch_size, doc_maxlen, hidden_dim)
        scores = self.score(hidden_outputs)
        # scores: (batch_size, doc_maxlen, doc_maxlen)
        context = scores.bmm(hidden_outputs)
        # context: (batch_size, doc_maxlen, hidden_dim)
        return context

    def score(self, hidden_outputs):
        # hidden_outputs: (batch_size, doc_maxlen, hidden_dim)
        if self.scoring == "dot":
            H = hidden_outputs.transpose(1,2)
            # H: (batch_size, hidden_dim, doc_maxlen)
            attention_energies = hidden_outputs.bmm(H)
            # (batch_size, doc_maxlen, doc_maxlen)
            scores = F.softmax(attention_energies, dim=2)
            return scores
        elif self.scoring == "general":
            H = self.W(hidden_outputs)
            # H: (batch_size, doc_maxlen, hidden_dim) with new hidden values
            H = H.transpose(1, 2)
            # H: (batch_size, hidden_dim, doc_maxlen)
            attention_energies = hidden_outputs.bmm(H)
            # (batch_size, doc_maxlen, doc_maxlen)
            scores = F.softmax(attention_energies, dim=2)
            return scores
        elif self.scoring == "concat":
            # hidden_outputs: (batch_size, doc_maxlen, hidden_dim)
            H = hidden_outputs.transpose(1,2)
            # H: (batch_size, hidden_dim, doc_maxlen)
            scores = []
            batch_size, doc_maxlen, hidden_dim = hidden_outputs.shape
            for doc_idx in range(H.shape[-1]):
                h_t = hidden_outputs[:,doc_idx, :]
                # h_t: (batch_size, hidden_dim)
                h_t = h_t.unsqueeze(1)
                # h_t: (batch_size, 1, hidden_dim)
                h_t = h_t.repeat(1, doc_maxlen, 1)
                # h_t: (batch_size, doc_maxlen, hidden_dim)
                # hidden_outputs: (batch_size, doc_maxlen, hidden_dim)
                H_t = torch.cat((h_t, hidden_outputs), dim=2)
                # H_t: (batch_size, doc_maxlen, 2 * hidden_dim)
                H_t = self.W(H_t)
                # H_t: (batch_size, doc_maxlen, hidden_dim) with new hidden values
                H_t = torch.nn.functional.tanh(H_t)
                # H_t: (batch_size, doc_maxlen, hidden_dim)
                H_t = self.v(H_t)
                # H_t: (batch_size, doc_maxlen, 1)
                H_t = H_t.view(batch_size, doc_maxlen)
                # H_t: (batch_size, doc_maxlen)
                scores.append(H_t)
            scores = torch.stack(scores)
            # scores: (doc_maxlen, batch_size, doc_maxlen)
            scores = scores.transpose(0, 1)
            # scaling trick: scaling 1/sqrt(d)
            scores = scores / torch.sqrt(torch.Tensor([hidden_dim]))
            # scores: (batch_size, doc_maxlen, doc_maxlen)
            scores = F.softmax(scores, dim=2)
            return scores
        else:
            raise RuntimeError("Unrecognized scoring method: %s" % self.scoring)

######################################################################
# Self Attention using Key, Value, and Query
# ==========================================
#
# Paper: https://arxiv.org/pdf/1706.03762.pdf
#
# This method uses the KVQ attention technique from
# the Transformer architecture's multi-head self-attention
#

class KVQSelfAttention(SelfAttention):

    def __init__(self, hidden_dim):
        super(KVQSelfAttention, self).__init__(hidden_dim, "KVQ")

    def initialize_layers(self):
        self.Wk = nn.Linear(self.hidden_dim, self.hidden_dim)
        self.Wv = nn.Linear(self.hidden_dim, self.hidden_dim)
        self.Wq = nn.Linear(self.hidden_dim, self.hidden_dim)
        # initialization
        self.init_linear(self.Wk)
        self.init_linear(self.Wv)
        self.init_linear(self.Wq)

    def forward(self, hidden_outputs):
        # hidden_outputs: (batch_size, doc_maxlen, hidden_dim)
        scores = self.score(hidden_outputs)
        # scores: (batch_size, doc_maxlen, doc_maxlen)
        hidden_outputs = self.Wv(hidden_outputs)
        # hidden_outputs: (batch_size, doc_maxlen, hidden_dim)
        context = scores.bmm(hidden_outputs)
        # context: (batch_size, doc_maxlen, hidden_dim)
        return context

    def score(self, hidden_outputs):
        K = self.Wk(hidden_outputs)
        Q = self.Wv(hidden_outputs)
        # K: (batch_size, doc_maxlen, hidden_dim) with new hidden values #1
        # Q: (batch_size, doc_maxlen, hidden_dim) with new hidden values #2
        Q = Q.transpose(1, 2)
        # K: (batch_size, doc_maxlen, hidden_dim) with new hidden values #1
        # Q: (batch_size, hidden_dim, doc_maxlen) with new hidden values #2
        attention_energies = K.bmm(Q)
        # (batch_size, doc_maxlen, doc_maxlen)
        scores = F.softmax(attention_energies, dim=2)
        # (batch_size, doc_maxlen, doc_maxlen)
        return scores

######################################################################
#  Multi-head Self Attention Mechanisms
# =============================================
#
# Paper: https://arxiv.org/pdf/1706.03762.pdf
#
# This method uses the multi-head self attention technique from
# the Transformer architecture.
#

class MultiHeadSelfAttention(SelfAttention):

    def __init__(self, hidden_dim, num_heads):
        self.num_heads = num_heads
        super(MultiHeadSelfAttention, self).__init__(hidden_dim, "multi-head")

    def initialize_layers(self):
        self.Wk = nn.ModuleList()
        self.Wv = nn.ModuleList()
        self.Wq = nn.ModuleList()
        for i in range(self.num_heads):
            self.Wk.append(nn.Linear(self.hidden_dim, self.hidden_dim))
            self.Wv.append(nn.Linear(self.hidden_dim, self.hidden_dim))
            self.Wq.append(nn.Linear(self.hidden_dim, self.hidden_dim))
        self.linear = nn.Linear(self.num_heads*self.hidden_dim, self.hidden_dim)
        # initialization
        for i in range(self.num_heads):
            self.init_linear(self.Wk[i])
            self.init_linear(self.Wv[i])
            self.init_linear(self.Wq[i])
        self.init_linear(self.linear)

    def forward(self, hidden_outputs):
        # hidden_outputs: (batch_size, doc_maxlen, hidden_dim)
        scores = self.score(hidden_outputs)
        # scores: (num_heads, batch_size, doc_maxlen, doc_maxlen)
        attention_outputs = []
        for i in range(self.num_heads):
            outputs_i = self.Wv[i](hidden_outputs)
            # outputs_i: (batch_size, doc_maxlen, hidden_dim)
            context_i = scores[i].bmm(outputs_i)
            # context_i: (batch_size, doc_maxlen, hidden_dim)
            attention_outputs.append(context_i)
        # (num_heads, batch_size, doc_maxlen, hidden_dim)
        attention_outputs = torch.cat(attention_outputs, dim=2)
        # (batch_size, doc_maxlen, hidden_dim * num_heads)
        context = self.linear(attention_outputs)
        # (batch_size, doc_maxlen, hidden_dim)
        return context

    def score(self, hidden_outputs):
        scores = []
        # technique from paper: scaled KVQ attention
        scale = torch.sqrt(torch.Tensor([hidden_outputs.shape[2]]))
        # move new tensor to the correct device, in case GPU is used
        scale = scale.to(hidden_outputs.device)
        # scoring for each heads
        for i in range(self.num_heads):
            K = self.Wk[i](hidden_outputs)
            # K: (batch_size, doc_maxlen, hidden_dim) with new hidden values #1
            Q = self.Wv[i](hidden_outputs)
            # Q: (batch_size, doc_maxlen, hidden_dim) with new hidden values #2
            Q = Q.transpose(1, 2)
            # K: (batch_size, doc_maxlen, hidden_dim) with new hidden values #1
            # Q: (batch_size, hidden_dim, doc_maxlen) with new hidden values #2
            attention_energies = K.bmm(Q)
            # (batch_size, doc_maxlen, doc_maxlen)
            scores.append(F.softmax(attention_energies / scale, dim=2))
            # (batch_size, doc_maxlen, doc_maxlen)
        return scores # (num_heads, batch_size, doc_maxlen, doc_maxlen)

##########################################################
# Basic Test
# ===========
#
# This only test if a model runs without any exceptions.
# It does NOT test for accuracy of the implementation.
#

def test_modules(batch_size=64, doc_maxlen=1000, hidden_dim=256,
                    num_heads=5, num_annotations=10):
    # dummy data
    hidden_outputs = torch.randn(size=(batch_size, doc_maxlen, hidden_dim))
    # test 0
    print("Testing Sentence Self Embedding")
    model = SentenceEmbedding(hidden_dim, hidden_dim, num_annotations)
    _ = model(hidden_outputs)
    # test 1
    print("Testing Self Attention with `general` scoring")
    model = SelfAttention(hidden_dim, scoring="general")
    _ = model(hidden_outputs)
    # test 2
    print("Testing Self Attention with `dot` scoring")
    model = SelfAttention(hidden_dim, scoring="dot")
    _ = model(hidden_outputs)
    # test 3
    print("Testing Self Attention with `concat` scoring")
    x = doc_maxlen
    if x >= 20:
        x = np.random.randint(2, 10)
        print("MODIFICATION: doc_maxlen is too large!")
        print("-- to make the test slower, a new doc_maxlen %d is used" % x)
    tmp = torch.randn(size=(batch_size, x, hidden_dim))
    model = SelfAttention(hidden_dim, scoring="concat")
    _ = model(tmp)
    # test 4
    print("Testing Self Attention with `other` scoring (expects exception)")
    exception_raised = False
    try:
        model = SelfAttention(hidden_dim, scoring="other")
    except RuntimeError as e:
        exception_raised = True
    assert exception_raised, "no exception is raised for unknown scoring"
    # test 5
    print("Testing KVQ Self Attention")
    model = KVQSelfAttention(hidden_dim)
    _ = model(hidden_outputs)
    # test 5
    print("Testing Multi-Head Self Attention")
    model = MultiHeadSelfAttention(hidden_dim, num_heads)
    _ = model(hidden_outputs)
    # end
    print("Test passed.")

max_stars_repo_name: Gan-Tu/ganutils

max_issues_repo_name: Gan-Tu/ganutils

max_forks_repo_name: Gan-Tu/ganutils


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 26 20:21:33 2019

@author: hc3
"""

import numpy as np

entradas = np.array([ [0,0], [0,1], [1,0], [1,1] ])
saidas = np.array([0,0,0,1])
pesos = np.array([0.0, 0.0])
taxaAprendizagem = 0.1

def stepFunction(soma):
    if(soma >= 1):
        return 1
    return 0

def calculaSaida(registro):
    s = registro.dot(pesos)
    return stepFunction(s)

def treinar():
    erroTotal = 1
    while(erroTotal != 0):
        erroTotal = 0
        for i in range(len(saidas)):
            saidaCalculada = calculaSaida(np.asarray(entradas[i]))
            erro = abs(saidas[i] - saidaCalculada)
            erroTotal += erro
            for j in range(len(pesos)):
                pesos[j] = pesos[j] + (taxaAprendizagem * entradas[i][j])
                print('Peso atualizado: '+ str(pesos[j]))

        print('Total de erros: '+ str(erroTotal))

treinar()


text:
# Copyright 2020 MONAI Consortium
# 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.

import unittest
import numpy as np
import torch
from parameterized import parameterized

from monai.inferers import sliding_window_inference

TEST_CASE_1 = [(1, 3, 16, 15, 7), (4, 10, 7), 3, 0.25, "constant"]  # 3D small roi

TEST_CASE_2 = [(1, 3, 16, 15, 7), (20, 22, 23), 10, 0.25, "constant"]  # 3D large roi

TEST_CASE_3 = [(1, 3, 15, 7), (2, 6), 1000, 0.25, "constant"]  # 2D small roi, large batch

TEST_CASE_4 = [(1, 3, 16, 7), (80, 50), 7, 0.25, "constant"]  # 2D large roi

TEST_CASE_5 = [(1, 3, 16, 15, 7), (20, 22, 23), 10, 0.5, "constant"]  # 3D large overlap

TEST_CASE_6 = [(1, 3, 16, 7), (80, 50), 7, 0.5, "gaussian"]  # 2D large overlap, gaussian

TEST_CASE_7 = [(1, 3, 16, 15, 7), (4, 10, 7), 3, 0.25, "gaussian"]  # 3D small roi, gaussian

class TestSlidingWindowInference(unittest.TestCase):
    @parameterized.expand([TEST_CASE_1, TEST_CASE_2, TEST_CASE_3, TEST_CASE_4, TEST_CASE_5, TEST_CASE_6, TEST_CASE_7])
    def test_sliding_window_default(self, image_shape, roi_shape, sw_batch_size, overlap, mode):
        inputs = torch.ones(*image_shape)
        device = torch.device("cpu:0")

        def compute(data):
            return data + 1

        result = sliding_window_inference(
            inputs.to(device), roi_shape, sw_batch_size, compute, overlap, blend_mode=mode
        )
        expected_val = np.ones(image_shape, dtype=np.float32) + 1
        self.assertTrue(np.allclose(result.numpy(), expected_val))

if __name__ == "__main__":
    unittest.main()

max_stars_repo_name: xidchen/MONAI

max_issues_repo_name: xidchen/MONAI

max_forks_repo_name: xidchen/MONAI


text:
from abc import ABC

import numpy as np

class Activation(object):
    def __call__(self, x):
        raise NotImplementedError

    def gradient(self, x):
        raise NotImplementedError

class Sigmoid(Activation):
    def __call__(self, x):
        return 1 / (1 + np.exp(-x))

    def gradient(self, x):
        return self.__call__(x) * (1 - self.__call__(x))

class ReLU(Activation):
    def __call__(self, x):
        return np.where(x >= 0, x, 0)

    def gradient(self, x):
        return np.where(x >= 0, 1, 0)

class TanH(Activation):
    def __call__(self, x):
        return (2 / (1 + np.exp(-2 * x))) - 1

    def gradient(self, x):
        return 1 - self.__call__(x) ** 2

class LeakyReLU(Activation):
    def __init__(self):
        self.alpha = .01

    def __call__(self, x):
        return np.where(x > 0, x, self.alpha * x)

    def gradient(self, x):
        return np.where(x > 0, 1, self.alpha)

class ELU(Activation):
    """ Exponential linear unit """

    def __init__(self, alpha=1.):
        self.alpha = alpha

    def __call__(self, x):
        return np.where(x > 0, x, self.alpha * (np.exp(x) - 1))

    def gradient(self, x):
        return np.where(x > 0, 1, self.__call__(x) + self.alpha)

class SELU(Activation):
    def __init__(self):
        self.alpha = 1.6732632423543772848170429916717
        self.l = 1.0507009873554804934193349852946

    def __call__(self, x):
        return self.l * np.where(x > 0, x, self.alpha * np.exp(x) - self.alpha)

    def gradient(self, x):
        return self.l * np.where(x > 0, 1, self.alpha * np.exp(x))

class Softmax(Activation):
    def __call__(self, x):
        e_x = np.exp(x)
        return e_x / np.sum(e_x, axis=1, keepdims=True)

    def gradient(self, x):
        p = self.__call__(x)
        return p * (1 - p)

class LogSoftmax(Activation, ABC):
    def __call__(self, x):
        e_x = np.exp(x)
        return np.log(e_x / np.sum(e_x, axis=1, keepdims=True))

class Linear(Activation):
    def __call__(self, x):
        return x

    def gradient(self, x):
        return np.ones((x.shape[0], x.shape[1]))


text:
"""
:Author: Pierre Barbier de Reuille <pierre.barbierdereuille@gmail.com>

This module specifically implement the curve fitting, wrapping the default
scipy.optimize.leastsq function. It allows for parameter value fixing,
different kind of residual and added constraints function.
"""

from __future__ import division, print_function, absolute_import
from scipy import optimize
from .compat import lrange
import numpy as np

class CurveFitting(object):
    r"""
    Fit a curve using the :py:func:`scipy.optimize.leastsq` function

    :type  xdata: ndarray
    :param xdata: Explaining values

    :type  ydata: ndarray
    :param ydata: Target values

    Once fitted, the following variables contain the result of
    the fitting:

    :ivar ndarray popt: The solution (or the result of the last iteration for
        an unsuccessful call)
    :ivar ndarray pcov: The estimated covariance of popt.  The diagonals
        provide the variance of the parameter estimate.
    :ivar ndarray res: Final residuals

    :ivar dict infodict: a dictionary of outputs with the keys:

            ``nfev``
                the number of function calls
            ``fvec``
                the function evaluated at the output
            ``fjac``
                A permutation of the R matrix of a QR factorization of
                the final approximate Jacobian matrix, stored column wise.
                Together with ipvt, the covariance of the estimate can be
                approximated.
            ``ipvt``
                an integer array of length N which defines a permutation
                matrix, ``p``, such that ``fjac*p = q*r``, where ``r`` is upper
                triangular with diagonal elements of nonincreasing
                magnitude. Column ``j`` of ``p`` is column ``ipvt(j)`` of the
                identity matrix.
            ``qtf``
                the vector ``(transpose(q) * fvec)``
            ``CI``
                list of tuple of parameters, each being the lower and
                upper bounds for the confidence interval in the CI
                argument at the same position.
            ``est_jacobian``
                True if the jacobian is estimated, false if the
                user-provided functions have been used

    .. note::

        In this implementation, residuals are supposed to be a generalisation
        of the notion of difference. In the end, the mathematical expression
        of this minimisation is:

        .. math::

            \hat{\theta} = \argmin_{\theta\in \mathbb{R}^p}
                \sum_i r(y_i, f(\theta, x_i))^2

        Where :math:`\theta` is the vector of :math:`p` parameters to optimise,
        :math:`r` is the residual function and :math:`f` is the function being
        fitted.
    """

    def __init__(self, xdata, ydata, **kwords):
        self._fct = None
        self._Dfun = None
        self._residuals = None
        self._Dres = None
        self._col_deriv = True
        self._constraints = None
        self._lsq_args = ()
        self._lsq_kwords = {}
        self._xdata = None
        self._ydata = None
        self._p0 = None
        self._fix_params = None

        self.xdata = xdata
        self.ydata = ydata

        self._fitted = False

        for n in kwords:
            setattr(self, n, kwords[n])

        if self._residuals is None:
            self._residuals = lambda x, y: (x - y)
            self._Dres = lambda y1, y0: -1

    def need_fit(self):
        """
        Function to be called if the object need to be fitted again
        """
        self._fitted = False

    @property
    def fitted(self):
        """
        Check if the object has been fitted or not
        """
        return self._fitted

    @property
    def function(self):
        """
        Function to be fitted. The call of the function will be::

            function(params, xs)
        """
        return self._fct

    @function.setter
    def function(self, f):
        self.need_fit()
        self._fct = f

    @property
    def Dfun(self):
        """
        Jacobian of the function with respect to its parameters.

        :Note: col_deriv defines if the derivative with respect to a given parameter is in column or row

        If not provided, a numerical approximation will be used instead.
        """
        return self._Dfun

    @Dfun.setter
    def Dfun(self, df):
        self.need_fit()
        self._Dfun = df

    @Dfun.deleter
    def Dfun(self):
        self.need_fit()
        self._Dfun = None

    @property
    def col_deriv(self):
        """
        Define if Dfun returns the derivatives by row or column.

        If ``col_deriv`` is ``True``, each line correspond to a parameter and each column to a point.
        """
        return self._col_deriv

    @col_deriv.setter
    def col_deriv(self, value):
        self._col_deriv = bool(value)
        self.need_fit()

    @property
    def residuals(self):
        """
        Residual function to use. The call will be::

            residuals(y_measured, y_est)

        where ``y_measured`` are the estimated values and ``y_est`` the measured ones.

        :Default: the defauls is ``y_measured - y_est``
        """
        return self._residuals

    @residuals.setter
    def residuals(self, f):
        self.need_fit()
        self._residuals = f

    @property
    def Dres(self):
        """
        Derivative of the residual function with respec to the estimated values. The call will be:

            Dres(y_measured, y_est)

        :Default: as the default residual is ``y_measured - y_est``, the default derivative is ``-1``
        """
        return self._Dres

    @Dres.setter
    def Dres(self, df):
        self.need_fit()
        self._Dres = df

    @Dres.deleter
    def Dres(self):
        self.need_fit()
        self._Dres = None

    @property
    def lsq_args(self):
        """
        Extra arguments to give to the least-square algorithm.

        See :py:func:`scipy.optimize.leastsq` for details
        """
        return self._lsq_args

    @lsq_args.setter
    def lsq_args(self, val):
        self.need_fit()
        self._lsq_args = tuple(val)

    @lsq_args.deleter
    def lsq_args(self):
        self._lsq_args = ()

    @property
    def lsq_kwords(self):
        """
        Extra named arguments to give to the least-square algorithm.

        See :py:func:`scipy.optimize.leastsq` for details
        """
        return self._lsq_kwords

    @lsq_kwords.setter
    def lsq_kwords(self, val):
        self.need_fit()
        self._lsq_kwords = dict(val)

    @lsq_kwords.deleter
    def lsq_kwords(self):
        self._lsq_kwords = {}

    @property
    def xdata(self):
        """
        Explaining values.
        """
        return self._xdata

    @xdata.setter
    def xdata(self, value):
        value = np.atleast_1d(value).squeeze()
        assert len(value.shape) < 3, "Error, xdata must be at most a 2D array"
        self._xdata = value
        self.need_fit()

    @property
    def ydata(self):
        """
        Target values.
        """
        return self._ydata

    @ydata.setter
    def ydata(self, value):
        value = np.atleast_1d(value).squeeze()
        assert len(value.shape) == 1, "Error, ydata must be at most a 1D array"
        self._ydata = value
        self.need_fit()

    @property
    def p0(self):
        """
        Initial fitting parameters
        """
        return self._p0

    @p0.setter
    def p0(self, value):
        value = np.atleast_1d(value)
        assert len(value.shape) == 1, "Error, p0 must be at most a 1D array"
        self._p0 = value

    @property
    def constraints(self):
        """
        Function returning additional constraints to the problem
        """
        return self._constraints

    @constraints.setter
    def constraints(self, value):
        assert callable(value), "Error, constraints must be a callable returning a 1d array"
        self._constraints = value

    @constraints.deleter
    def constraints(self):
        self._constraints = None

    @property
    def fix_params(self):
        """
        Index of parameters that shouldn't be touched by the algorithm
        """
        return self._fix_params

    @fix_params.setter
    def fix_params(self, value):
        self._fix_params = tuple(value)

    @fix_params.deleter
    def fix_params(self):
        self._fix_params = None

    def fit(self):
        """
        Fit the curve
        """
        Dres = self.Dres
        Dfun = self.Dfun
        fct = self.function
        residuals = self.residuals
        col_deriv = self.col_deriv
        p0 = self.p0
        xdata = self.xdata
        ydata = self.ydata
        fix_params = self.fix_params

        use_derivs = (Dres is not None) and (Dfun is not None)
        df = None
        f = None

        if fix_params:
            p_save = np.array(p0, dtype=float)
            change_params = lrange(len(p0))
            try:
                for i in fix_params:
                    change_params.remove(i)
            except ValueError:
                raise ValueError("List of parameters to fix is incorrect: "
                                 "contains either duplicates or values "
                                 "out of range.")
            p0 = p_save[change_params]

            def f_fixed(p):
                p1 = np.array(p_save)
                p1[change_params] = p
                y0 = fct(p1, xdata)
                return residuals(ydata, y0)
            f = f_fixed
            if use_derivs:
                def df_fixed(p):
                    p1 = np.array(p_save)
                    p1[change_params] = p
                    y0 = fct(p1, xdata)
                    dfct = Dfun(p1, xdata)
                    dr = Dres(ydata, y0)
                    if col_deriv:
                        return dfct[change_params]*dr
                    return dfct[:,change_params]*dr[:, np.newaxis]
                df = df_fixed
        else:
            def f_free(p):
                y0 = fct(p, xdata)
                return residuals(ydata, y0)
            f = f_free
            if use_derivs:
                def df_free(p):
                    dfct = Dfun(p, xdata)
                    y0 = fct(p, xdata)
                    dr = np.atleast_1d(Dres(ydata, y0))
                    if col_deriv:
                        return dfct*dr
                    return dfct*dr[:, np.newaxis]
                df = df_free

        if use_derivs:
            self.df = df

        cd = 1 if col_deriv else 0
        optim = optimize.leastsq(f, p0, full_output=1, Dfun=df,
                                 col_deriv=cd, *self.lsq_args, **self.lsq_kwords)
        popt, pcov, infodict, mesg, ier = optim
        #infodict['est_jacobian'] = not use_derivs

        if fix_params:
            p_save[change_params] = popt
            popt = p_save

        if not ier in [1, 2, 3, 4]:
            raise RuntimeError("Unable to determine number of fit parameters. "
                               "Error returned by scipy.optimize.leastsq:\n%s"
                               % (mesg,))

        res = residuals(ydata, fct(popt, xdata))
        if (len(res) > len(p0)) and pcov is not None:
            s_sq = (res ** 2).sum() / (len(ydata) - len(p0))
            pcov = pcov * s_sq
        else:
            pcov = np.inf

        self.popt = popt
        self.pcov = pcov
        self.res = res
        self.infodict = infodict
        self._fitted = True

    def __call__(self, xdata):
        """
        Return the value of the fitted function for each of the points in
        ``xdata``
        """
        if not self.fitted:
            self.fit()
        return self.function(self.popt, xdata)

max_stars_repo_name: GPrathap/OpenBCIPython

max_issues_repo_name: GPrathap/OpenBCIPython

max_forks_repo_name: GPrathap/OpenBCIPython

reason: import numpy,from scipy


text:
import rasterio
from rasterio import warp
from rasterio.mask import mask
from rasterio.plot import show, reshape_as_raster, reshape_as_image, adjust_band
from shapely.geometry import box
from fiona.crs import from_epsg
import geopandas as gpd
import numpy as np
import cv2 as cv
import pycrs
from optparse import OptionParser
import matplotlib.pyplot as plt
import uuid
import os

'''
    Overview:
        raster(R, G, B, NIR, null) ---->  | naip2map | ----> raster(R, G, B, NIR, Binary Occupancy)
'''

def window_from_extent(xmin, xmax, ymin, ymax, aff):
    # Source: Loic Dutrieux
    # https://gis.stackexchange.com/a/244727
    col_start, row_start = ~aff * (xmin, ymax)
    col_stop, row_stop = ~aff * (xmax, ymin)
    return ((int(row_start), int(row_stop)), (int(col_start), int(col_stop)))

# Options
parser = OptionParser()
parser.add_option("-i", "--in_file",
                  help = "File with raster of bands (red, green, blue, NIR, null)",
                  default = "in/naip18-nc-cir-60cm_2797141_20181210.jp2")
parser.add_option("-o", "--out_file",
                  help = "File to store output raster (fifth band changed to occupancy grid)",
                  default = "out/naip18-nc-cir-60cm_2797141_20181210_occ.jp2")
parser.add_option("-p", "--out_file_plot",
                  help = "File to store output occupancy grid plot",
                  default = "out/naip18-nc-cir-60cm_2797141_20181210_occ.png")
parser.add_option("-b", "--bounds",
                  help = "Comma-separated coordinates to crop raster. Ex: 'ymin,xmin,ymax,xmax'",
                  default = "27.849235,-97.361426,27.875268,-97.325000")
parser.add_option("-n", "--band",
                  help = "Index of raster band with NIR data", type = "int",
                  default = 4)
parser.add_option("-s", "--start_coords",
                  help = "Comma-separated pair of lat,lon coords to determine accessable map area. Ex: '27.870539,-97.334531'",
                  default = None)

(options, args) = parser.parse_args()
inFile = options.in_file
outFile = options.out_file
outFilePlotOcc = options.out_file_plot
bounds = startCoords = None
nirBand = options.band
if options.bounds is not None:
    bounds = [float(b) for b in options.bounds.split(",")]
if options.start_coords is not None:
    startCoords = [float(s) for s in options.start_coords.split(",")]
    start_latlon = { "lat" : startCoords[0], "lon": startCoords[1] }

# Open NAIP raster
aRaster = rasterio.open(inFile, driver = "JP2OpenJPEG")

# Crop raster
bbox = box(bounds[1] + 0.5, bounds[0] + 0.5, bounds[3] - 0.5, bounds[2] - 0.5)
geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(4326))
geo = geo.to_crs(crs=aRaster.crs.data)

def getFeatures(gdf):
    import json
    return [json.loads(gdf.to_json())['features'][0]['geometry']]

coords = getFeatures(geo)
out_img, out_transform = mask(aRaster, shapes = coords, crop = True, nodata = 255)
out_meta = aRaster.meta.copy()
epsg_code = int(aRaster.crs.data['init'][5:])

out_meta.update({"driver": "JP2OpenJPEG",
                 "height": out_img.shape[1],
                 "width": out_img.shape[2],
                 "transform": out_transform,
                 "crs": aRaster.crs})

tempFile = str(uuid.uuid4()) + ".jp2"
with rasterio.open(tempFile, "w", **out_meta) as dest:
    dest.write(out_img)
aRaster.close()

aRaster = rasterio.open(tempFile)
aRows, aCols = aRaster.shape

# Raster information
print("Number of bands: {n}".format(n = aRaster.count))
print("Image size: {r} rows x {c} columns".format(r = aRows, c = aCols))
print("Raster description: {desc}".format(desc = aRaster.descriptions))
print("Image projection: {proj}".format(proj = aRaster.crs))
print("Image geo-transform: {gt}".format(gt = aRaster.transform))
print("Image extent: {bounds}".format(bounds = aRaster.bounds))

# Open image
img = aRaster.read().astype(np.uint8)

# Select NIR band
nir = img[nirBand - 1, :, :]
nirRows, nirCols = nir.shape

# Filter
nir[nir > 70] = 255
nir[nir < 50] = 0

# Dilation
kernel = np.ones((5, 5), np.uint8)
dilation = cv.dilate(nir, kernel, iterations = 1)

# Final map
aMap = dilation.copy()

if startCoords is None:
    # Make new raster
    with rasterio.open(outFile, 'w', **aRaster.meta) as dst:
        # Copy existing bands
        for ID, b in enumerate([x + 1 for x in list(range(aRaster.count))], 1):
            dst.write(aRaster.read(b), ID)
            dst.set_band_description(ID, aRaster.descriptions[ID - 1])

        # Add occupancy grid band
        dst.write(aMap, aRaster.count)
        dst.set_band_description(5, 'Occupancy')

    # Save plot of occupancy grid
    if outFilePlotOcc is not None:
        fig, ax = plt.subplots(1, 1)
        ax.imshow(aMap, cmap = 'Blues_r', alpha = 0.5)
        plt.title("Occupancy Grid")
        plt.tight_layout()
        fig.savefig(outFilePlotOcc)
        plt.close()

    os.remove(tempFile)
    exit(0)

#################################################
# Limit to region accessible from a start point #
#################################################

# Convert coordinates
start_local = warp.transform({'init' : 'epsg:4326'}, aRaster.crs, [start_latlon["lon"]], [start_latlon["lat"]])
py, px = aRaster.index(start_local[0], start_local[1])

# Determine accessible region via flood fill
filled = aMap.copy()
regionMask = np.zeros((nirRows + 2, nirCols + 2), np.uint8)
cv.floodFill(filled, regionMask, (px[0], py[0]), 200, loDiff = (50, 50, 50, 50), upDiff = (50, 50, 50, 50))
regionMask = regionMask[1:1+nirRows, 1:1+nirCols]

# Convert all inaccessible areas to land
regionMaskRev = regionMask.copy()
regionMaskRev[regionMaskRev > 0] = 100
regionMaskRev[regionMaskRev == 0] = 255
regionMaskRev[regionMaskRev != 255] = aMap[regionMaskRev != 255]

# Make new raster
with rasterio.open(outFile, 'w', **aRaster.meta) as dst:
    # Copy existing bands
    for ID, b in enumerate([x + 1 for x in list(range(aRaster.count))], 1):
        dst.write(aRaster.read(b), ID)
        dst.set_band_description(ID, aRaster.descriptions[ID - 1])

    # Add occupancy grid band
    dst.write(regionMaskRev, aRaster.count)
    dst.set_band_description(5, 'Occupancy')

# Save plot
if outFilePlotOcc is not None:
    fig, axs = plt.subplots(1, 2)
    axs[0].imshow(aMap, cmap = 'Blues_r', alpha = 0.5)
    axs[0].set_title("Occupancy grid")
    axs[1].imshow(regionMaskRev, cmap = 'Blues_r', alpha = 0.5)
    axs[1].scatter(px, py, s = 5, c = 'red', marker = 'x')
    axs[1].set_title("Accessible grid")
    plt.tight_layout()
    fig.savefig(outFilePlotOcc)

os.remove(tempFile)