text:
import os
import tensorflow as tf
import numpy as np

from modules.modules_ambient_pmd import sampler, \
        generator, pitch_discriminator, mfc_discriminator
from utils.model_utils import l1_loss
from utils.tf_forward_tan import lddmm

class VariationalCycleGAN(object):

    def __init__(self, dim_pitch=1, dim_mfc=23, n_frames=128,
            pitch_discriminator=pitch_discriminator,
            mfc_discriminator=mfc_discriminator,
            generator=generator, sampler=sampler,
            lddmm=lddmm, mode='train',
            log_file_name='no_name', pre_train=None):

        self.n_frames = n_frames
        self.pitch_shape = [None, dim_pitch, None] #[batch_size, num_features, num_frames]
        self.mfc_shape = [None, dim_mfc, None]

        self.first_order_diff_mat = np.eye(self.n_frames, dtype=np.float32)
        for i in range(1, self.n_frames):
            self.first_order_diff_mat[i-1,i] = -1

        # Create the kernel for lddmm
        self.kernel = tf.expand_dims(tf.constant([6,50],
            dtype=tf.float32), axis=0)

        self.sampler = sampler
        self.generator = generator
        self.pitch_discriminator = pitch_discriminator
        self.mfc_discriminator = mfc_discriminator
        self.lddmm = lddmm
        self.mode = mode

        self.build_model()
        self.optimizer_initializer()

        self.saver = tf.train.Saver()
        self.sess = tf.Session()
        if pre_train is not None:
            self.saver.restore(self.sess, pre_train)
        else:
            self.sess.run(tf.global_variables_initializer())

        if self.mode == 'train':
            self.train_step = 0
            self.writer = tf.summary.FileWriter('./tensorboard_log/'+log_file_name,
                    tf.get_default_graph())
            self.generator_summaries, self.discriminator_summaries = self.summary()

    def build_model(self):

        # Placeholders for real training samples
        self.pitch_A_real = tf.placeholder(tf.float32, shape=self.pitch_shape,
                name='pitch_A_real')
        self.pitch_B_real = tf.placeholder(tf.float32, shape=self.pitch_shape,
                name='pitch_B_real')

        self.mfc_A_real = tf.placeholder(tf.float32, shape=self.mfc_shape,
                name='mfc_A_real')
        self.mfc_B_real = tf.placeholder(tf.float32, shape=self.mfc_shape,
                name='mfc_B_real')

        # Placeholders for fake generated samples
        self.pitch_A_fake = tf.placeholder(tf.float32, shape=self.pitch_shape,
                name='pitch_A_fake')
        self.pitch_B_fake = tf.placeholder(tf.float32, shape=self.pitch_shape,
                name='pitch_B_fake')

        self.mfc_A_fake = tf.placeholder(tf.float32, shape=self.mfc_shape,
                name='mfc_A_fake')
        self.mfc_B_fake = tf.placeholder(tf.float32, shape=self.mfc_shape,
                name='mfc_B_fake')

        # Placeholder for test samples
        self.pitch_A_test = tf.placeholder(tf.float32, shape=self.pitch_shape,
                name='pitch_A_test')
        self.mfc_A_test = tf.placeholder(tf.float32, shape=self.mfc_shape,
                name='mfc_A_test')

        self.pitch_B_test = tf.placeholder(tf.float32, shape=self.pitch_shape,
                name='pitch_B_test')
        self.mfc_B_test = tf.placeholder(tf.float32, shape=self.mfc_shape,
                name='mfc_B_test')

        # Place holder for lambda_cycle and lambda_identity
        self.lambda_cycle_pitch = tf.placeholder(tf.float32, None,
                name='lambda_cycle_pitch')
        self.lambda_cycle_mfc = tf.placeholder(tf.float32, None,
                name='lambda_cycle_mfc')
        self.lambda_momenta = tf.placeholder(tf.float32, None,
                name='lambda_momenta')
        self.lambda_identity_mfc = tf.placeholder(tf.float32, None,
                name='lambda_identity_mfc')

        '''
        Generator A
        '''
        # Generate pitch from A to B
        self.momenta_generation_A2B = self.sampler(input_pitch=self.pitch_A_real,
                input_mfc=self.mfc_A_real, reuse=False, scope_name='sampler_A2B')
        self.pitch_generation_A2B = self.lddmm(x=self.pitch_A_real,
                p=self.momenta_generation_A2B, kernel=self.kernel, reuse=False, scope_name='lddmm')
        self.mfc_generation_A2B = self.generator(input_pitch=self.pitch_generation_A2B,
                input_mfc=self.mfc_A_real, reuse=False, scope_name='generator_A2B')
        # Cyclic generation
        self.momenta_cycle_A2A = self.sampler(input_pitch=self.pitch_generation_A2B,
                input_mfc=self.mfc_generation_A2B, reuse=False, scope_name='sampler_B2A')
        self.pitch_cycle_A2A = self.lddmm(x=self.pitch_generation_A2B,
                p=self.momenta_cycle_A2A, kernel=self.kernel, reuse=True, scope_name='lddmm')
        self.mfc_cycle_A2A = self.generator(input_pitch=self.pitch_cycle_A2A,
                input_mfc=self.mfc_generation_A2B, reuse=False, scope_name='generator_B2A')
        self.mfc_identity_A2B = self.generator(input_pitch=self.pitch_B_real,
                input_mfc=self.mfc_B_real, reuse=True, scope_name='generator_A2B')

        '''
        Generator B
        '''
        # Generate pitch from B to A
        self.momenta_generation_B2A = self.sampler(input_pitch=self.pitch_B_real,
                input_mfc=self.mfc_B_real, reuse=True, scope_name='sampler_B2A')
        self.pitch_generation_B2A = self.lddmm(x=self.pitch_B_real,
                p=self.momenta_generation_B2A, kernel=self.kernel, reuse=True, scope_name='lddmm')
        self.mfc_generation_B2A = self.generator(input_pitch=self.pitch_generation_B2A,
                input_mfc=self.mfc_B_real, reuse=True, scope_name='generator_B2A')
        # Cyclic generation
        self.momenta_cycle_B2B = self.sampler(input_pitch=self.pitch_generation_B2A,
                input_mfc=self.mfc_generation_B2A, reuse=True, scope_name='sampler_A2B')
        self.pitch_cycle_B2B = self.lddmm(x=self.pitch_generation_B2A,
                p=self.momenta_cycle_B2B, kernel=self.kernel, reuse=True, scope_name='lddmm')
        self.mfc_cycle_B2B = self.generator(input_pitch=self.pitch_cycle_B2B,
                input_mfc=self.mfc_generation_B2A, reuse=True, scope_name='generator_A2B')
        self.mfc_identity_B2A = self.generator(input_pitch=self.pitch_A_real,
                input_mfc=self.mfc_A_real, reuse=True, scope_name='generator_B2A')

        '''
        Initialize the pitch discriminators
        '''
        # Discriminator initialized to keep parameters in memory
        self.pitch_discrimination_B_fake = self.pitch_discriminator(input_pitch=tf.concat([self.pitch_A_real,
                self.pitch_generation_A2B], axis=1), reuse=False, scope_name='pitch_discriminator_A')

        self.pitch_discrimination_A_fake = self.pitch_discriminator(input_pitch=tf.concat([self.pitch_B_real,
                self.pitch_generation_B2A], axis=1), reuse=False, scope_name='pitch_discriminator_B')

        '''
        Initialize the mfc discriminators
        '''
        # Discriminator initialized to keep parameters in memory
        self.mfc_discrimination_B_fake = self.mfc_discriminator(input_mfc=tf.concat([self.mfc_A_real,
                self.mfc_generation_A2B], axis=1), reuse=False, scope_name='mfc_discriminator_A')

        self.mfc_discrimination_A_fake = self.mfc_discriminator(input_mfc=tf.concat([self.mfc_B_real,
                self.mfc_generation_B2A], axis=1), reuse=False, scope_name='mfc_discriminator_B')

        '''
        Computing loss for generators
        '''
        # Cycle loss
        self.cycle_loss_pitch = (l1_loss(y=self.pitch_A_real,
            y_hat=self.pitch_cycle_A2A) + l1_loss(y=self.pitch_B_real,
                y_hat=self.pitch_cycle_B2B)) / 2.0

        self.cycle_loss_mfc = (l1_loss(y=self.mfc_A_real,
            y_hat=self.mfc_cycle_A2A) + l1_loss(y=self.mfc_B_real,
                y_hat=self.mfc_cycle_B2B)) / 2.0

        # Identity Loss
        self.identity_loss_mfc = (l1_loss(y=self.mfc_identity_A2B,
            y_hat=self.mfc_B_real) + l1_loss(y=self.mfc_identity_B2A,
                    y_hat=self.mfc_A_real)) / 2.0

        # Sampler-Generator loss
        # Sampler-Generator wants to fool discriminator
        self.generator_loss_A2B = (l1_loss(y=tf.ones_like(self.pitch_discrimination_B_fake),
                y_hat=self.pitch_discrimination_B_fake) \
                + l1_loss(y=tf.ones_like(self.mfc_discrimination_B_fake),
                          y_hat=self.mfc_discrimination_B_fake)) / 2.0
        self.generator_loss_B2A = (l1_loss(y=tf.ones_like(self.pitch_discrimination_A_fake),
                y_hat=self.pitch_discrimination_A_fake) \
                + l1_loss(y=tf.ones_like(self.mfc_discrimination_A_fake),
                          y_hat=self.mfc_discrimination_A_fake)) / 2.0

        self.gen_disc_loss = (self.generator_loss_A2B + self.generator_loss_B2A) / 2.0

        self.momenta_loss_A2B = tf.reduce_sum(tf.square(tf.matmul(self.first_order_diff_mat,
            tf.reshape(self.momenta_generation_A2B, [-1,1])))) \
                    + tf.reduce_sum(tf.square(tf.matmul(self.first_order_diff_mat,
                        tf.reshape(self.momenta_cycle_A2A, [-1,1]))))

        self.momenta_loss_B2A = tf.reduce_sum(tf.square(tf.matmul(self.first_order_diff_mat,
            tf.reshape(self.momenta_generation_B2A, [-1,1])))) \
                    + tf.reduce_sum(tf.square(tf.matmul(self.first_order_diff_mat,
                        tf.reshape(self.momenta_cycle_B2B, [-1,1]))))

        self.momenta_loss = (self.momenta_loss_A2B + self.momenta_loss_B2A) / 2.0

        # Merge the two sampler-generator, the cycle loss and momenta prior
        self.generator_loss \
            = self.gen_disc_loss + self.lambda_cycle_pitch * self.cycle_loss_pitch \
                + self.lambda_cycle_mfc * self.cycle_loss_mfc \
                + self.lambda_momenta * self.momenta_loss \
                + self.lambda_identity_mfc * self.identity_loss_mfc

        # Compute the pitch discriminator probability for pair of inputs
        self.pitch_discrimination_input_A_real_B_fake \
            = self.pitch_discriminator(input_pitch=tf.concat([self.pitch_A_real, self.pitch_B_fake], axis=1),
                    reuse=True, scope_name='pitch_discriminator_A')
        self.pitch_discrimination_input_A_fake_B_real \
            = self.pitch_discriminator(input_pitch=tf.concat([self.pitch_A_fake, self.pitch_B_real], axis=1),
                    reuse=True, scope_name='pitch_discriminator_A')

        self.pitch_discrimination_input_B_real_A_fake \
            = self.pitch_discriminator(input_pitch=tf.concat([self.pitch_B_real, self.pitch_A_fake], axis=1),
                    reuse=True, scope_name='pitch_discriminator_B')
        self.pitch_discrimination_input_B_fake_A_real \
            = self.pitch_discriminator(input_pitch=tf.concat([self.pitch_B_fake, self.pitch_A_real], axis=1),
                    reuse=True, scope_name='pitch_discriminator_B')

        # Compute pitch discriminator loss for backprop
        self.pitch_discriminator_loss_input_A_real \
            = l1_loss(y=tf.zeros_like(self.pitch_discrimination_input_A_real_B_fake),
                    y_hat=self.pitch_discrimination_input_A_real_B_fake)
        self.pitch_discriminator_loss_input_A_fake \
            = l1_loss(y=tf.ones_like(self.pitch_discrimination_input_A_fake_B_real),
                    y_hat=self.pitch_discrimination_input_A_fake_B_real)
        self.pitch_discriminator_loss_A = (self.pitch_discriminator_loss_input_A_real \
                                     + self.pitch_discriminator_loss_input_A_fake) / 2.0

        self.pitch_discriminator_loss_input_B_real \
            = l1_loss(y=tf.zeros_like(self.pitch_discrimination_input_B_real_A_fake),
                    y_hat=self.pitch_discrimination_input_B_real_A_fake)
        self.pitch_discriminator_loss_input_B_fake \
            = l1_loss(y=tf.ones_like(self.pitch_discrimination_input_B_fake_A_real),
                    y_hat=self.pitch_discrimination_input_B_fake_A_real)
        self.pitch_discriminator_loss_B = (self.pitch_discriminator_loss_input_B_real \
                                     + self.pitch_discriminator_loss_input_B_fake) / 2.0

        # Merge the two pitch discriminators into one
        self.pitch_discriminator_loss = (self.pitch_discriminator_loss_A \
                + self.pitch_discriminator_loss_B) / 2.0

        # Compute the mfcc discriminator probability for mfcc
        self.mfc_discrimination_input_A_real_B_fake \
            = self.mfc_discriminator(input_mfc=tf.concat([self.mfc_A_real, self.mfc_B_fake], axis=1),
                    reuse=True, scope_name='mfc_discriminator_A')
        self.mfc_discrimination_input_A_fake_B_real \
            = self.mfc_discriminator(input_mfc=tf.concat([self.mfc_A_fake, self.mfc_B_real], axis=1),
                    reuse=True, scope_name='mfc_discriminator_A')

        self.mfc_discrimination_input_B_real_A_fake \
            = self.mfc_discriminator(input_mfc=tf.concat([self.mfc_B_real, self.mfc_A_fake], axis=1),
                    reuse=True, scope_name='mfc_discriminator_B')
        self.mfc_discrimination_input_B_fake_A_real \
            = self.mfc_discriminator(input_mfc=tf.concat([self.mfc_B_fake, self.mfc_A_real], axis=1),
                    reuse=True, scope_name='mfc_discriminator_B')

        # Compute mfcc discriminator loss for backprop
        self.mfc_discriminator_loss_input_A_real \
            = l1_loss(y=tf.zeros_like(self.mfc_discrimination_input_A_real_B_fake),
                    y_hat=self.mfc_discrimination_input_A_real_B_fake)
        self.mfc_discriminator_loss_input_A_fake \
            = l1_loss(y=tf.ones_like(self.mfc_discrimination_input_A_fake_B_real),
                    y_hat=self.mfc_discrimination_input_A_fake_B_real)
        self.mfc_discriminator_loss_A = (self.mfc_discriminator_loss_input_A_real \
                                     + self.mfc_discriminator_loss_input_A_fake) / 2.0

        self.mfc_discriminator_loss_input_B_real \
            = l1_loss(y=tf.zeros_like(self.mfc_discrimination_input_B_real_A_fake),
                    y_hat=self.mfc_discrimination_input_B_real_A_fake)
        self.mfc_discriminator_loss_input_B_fake \
            = l1_loss(y=tf.ones_like(self.mfc_discrimination_input_B_fake_A_real),
                    y_hat=self.mfc_discrimination_input_B_fake_A_real)
        self.mfc_discriminator_loss_B = (self.mfc_discriminator_loss_input_B_real \
                                     + self.mfc_discriminator_loss_input_B_fake) / 2.0

        # Merge the two mfc discriminators into one
        self.mfc_discriminator_loss = (self.mfc_discriminator_loss_A \
                + self.mfc_discriminator_loss_B) / 2.0

        # Final merging of pitch and mfc discriminators
        self.discriminator_loss = (self.pitch_discriminator_loss + self.mfc_discriminator_loss) / 2.0

        # Categorize variables to optimize the two sets separately
        trainable_variables = tf.trainable_variables()
        self.discriminator_vars = [var for var in trainable_variables if 'discriminator' in var.name]
        self.generator_vars = [var for var in trainable_variables if 'generator' in var.name \
                                                                    or 'sampler' in var.name]

        # Reserved for test
        self.momenta_A2B_test = self.sampler(input_pitch=self.pitch_A_test,
                input_mfc=self.mfc_A_test, reuse=True, scope_name='sampler_A2B')
        self.pitch_A2B_test = self.lddmm(x=self.pitch_A_test,
                p=self.momenta_A2B_test, kernel=self.kernel, reuse=True, scope_name='lddmm')
        self.mfc_A2B_test = self.generator(input_pitch=self.pitch_A2B_test,
                input_mfc=self.mfc_A_test, reuse=True, scope_name='generator_A2B')

        self.momenta_B2A_test = self.sampler(input_pitch=self.pitch_B_test,
                input_mfc=self.mfc_B_test, reuse=True, scope_name='sampler_B2A')
        self.pitch_B2A_test = self.lddmm(x=self.pitch_B_test,
                p=self.momenta_B2A_test, kernel=self.kernel, reuse=True, scope_name='lddmm')
        self.mfc_B2A_test = self.generator(input_pitch=self.pitch_B2A_test,
                input_mfc=self.mfc_B_test, reuse=True, scope_name='generator_B2A')

    def optimizer_initializer(self):

        self.generator_learning_rate = tf.placeholder(tf.float32, None,
                name='generator_learning_rate')
        self.discriminator_learning_rate = tf.placeholder(tf.float32, None,
                name='discriminator_learning_rate')

        self.discriminator_optimizer \
                = tf.train.AdamOptimizer(learning_rate=self.discriminator_learning_rate,
                        beta1=0.5).minimize(self.discriminator_loss, var_list=self.discriminator_vars)

        self.generator_optimizer \
                = tf.train.AdamOptimizer(learning_rate=self.generator_learning_rate,
                        beta1=0.5).minimize(self.generator_loss, var_list=self.generator_vars)

    def train(self, pitch_A, mfc_A, pitch_B, mfc_B, lambda_cycle_pitch, lambda_cycle_mfc,
              lambda_momenta, lambda_identity_mfc, generator_learning_rate,
              discriminator_learning_rate):

        generated_momenta_B, generated_pitch_B, generated_mfc_B, \
                generated_momenta_A, generated_pitch_A, generated_mfc_A, \
                generator_loss, _, generator_summaries \
                = self.sess.run([self.momenta_generation_A2B, self.pitch_generation_A2B,
                    self.mfc_generation_A2B, self.momenta_generation_B2A, self.pitch_generation_B2A,
                    self.mfc_generation_B2A, self.gen_disc_loss, self.generator_optimizer,
                    self.generator_summaries],
                    feed_dict = {self.lambda_cycle_pitch:lambda_cycle_pitch,
                        self.lambda_cycle_mfc:lambda_cycle_mfc,
                        self.lambda_momenta:lambda_momenta,
                        self.lambda_identity_mfc:lambda_identity_mfc,
                        self.pitch_A_real:pitch_A, self.mfc_A_real:mfc_A,
                        self.pitch_B_real:pitch_B, self.mfc_B_real:mfc_B,
                        self.generator_learning_rate:generator_learning_rate})

        discriminator_loss, _, discriminator_summaries \
            = self.sess.run([self.discriminator_loss, self.discriminator_optimizer,
                self.discriminator_summaries],
                    feed_dict = {self.pitch_A_real:pitch_A, self.pitch_B_real:pitch_B,
                        self.mfc_A_real:mfc_A, self.mfc_B_real:mfc_B,
                        self.discriminator_learning_rate:discriminator_learning_rate,
                        self.pitch_A_fake:generated_pitch_A, self.pitch_B_fake:generated_pitch_B,
                        self.mfc_A_fake:generated_mfc_A, self.mfc_B_fake:generated_mfc_B})

        self.writer.add_summary(generator_summaries, self.train_step)
        self.writer.add_summary(discriminator_summaries, self.train_step)

        self.train_step += 1

        return generator_loss, discriminator_loss, generated_pitch_A, generated_mfc_A, \
                generated_pitch_B, generated_mfc_B, generated_momenta_A, generated_momenta_B

    def test_gen(self, mfc_A, pitch_A, mfc_B, pitch_B):
        gen_mom_B, gen_pitch_B, gen_mfc_B, = self.sess.run([self.momenta_A2B_test, \
                                    self.pitch_A2B_test, self.mfc_A2B_test], \
                                    feed_dict={self.pitch_A_test:pitch_A, \
                                        self.mfc_A_test:mfc_A})

        gen_mom_A, gen_pitch_A, gen_mfc_A = self.sess.run([self.momenta_B2A_test, \
                                    self.pitch_B2A_test, self.mfc_B2A_test], \
                                    feed_dict={self.pitch_B_test:pitch_B, \
                                        self.mfc_B_test:mfc_B})

        return gen_mom_A, gen_pitch_A, gen_mfc_A, gen_mom_B, gen_pitch_B, gen_mfc_B

    def test(self, input_pitch, input_mfc, direction):

        if direction == 'A2B':
            generated_pitch, generated_mfc = self.sess.run([self.pitch_A2B_test,
                self.mfc_A2B_test], feed_dict = {self.pitch_A_test:input_pitch,
                    self.mfc_A_test:input_mfc})

        elif direction == 'B2A':
            generated_pitch, generated_mfc = self.sess.run([self.pitch_B2A_test,
                self.mfc_B2A_test], feed_dict = {self.pitch_B_test:input_pitch,
                    self.mfc_B_test:input_mfc})
        else:
            raise Exception('Conversion direction must be specified.')

        return generated_pitch, generated_mfc

    def save(self, directory, filename):

        if not os.path.exists(directory):
            os.makedirs(directory)
        self.saver.save(self.sess, \
                        os.path.join(directory, filename))

    def load(self, filepath):

        self.saver.restore(self.sess, filepath)

    def summary(self):

        with tf.name_scope('generator_summaries'):
            cycle_loss_pitch_summary = tf.summary.scalar('cycle_loss_pitch', \
                                    self.cycle_loss_pitch)
            cycle_loss_mfc_summary = tf.summary.scalar('cycle_loss_mfc', \
                                    self.cycle_loss_mfc)
            identity_loss_summary = tf.summary.scalar('identity_loss_mfc', \
                                    self.identity_loss_mfc)
            generator_loss_A2B_summary = tf.summary.scalar('generator_loss_A2B', \
                                    self.generator_loss_A2B)
            generator_loss_B2A_summary = tf.summary.scalar('generator_loss_B2A', \
                                    self.generator_loss_B2A)
            generator_loss_summary = tf.summary.scalar('generator_loss', \
                                    self.generator_loss)
            generator_summaries = tf.summary.merge([cycle_loss_pitch_summary,
                                    cycle_loss_mfc_summary, identity_loss_summary,
                                    generator_loss_A2B_summary, generator_loss_B2A_summary,
                                    generator_loss_summary])

        with tf.name_scope('discriminator_summaries'):
            discriminator_loss_pitch_summary \
                = tf.summary.scalar('discriminator_loss_pitch', \
                        self.pitch_discriminator_loss)
            discriminator_loss_mfc_summary \
                = tf.summary.scalar('discriminator_loss_mfc', \
                        self.mfc_discriminator_loss)
            discriminator_loss_summary \
                = tf.summary.scalar('discriminator_loss', \
                        self.discriminator_loss)
            discriminator_summaries \
                = tf.summary.merge([discriminator_loss_pitch_summary, \
                        discriminator_loss_mfc_summary, \
                        discriminator_loss_summary])

        return generator_summaries, discriminator_summaries

if __name__ == '__main__':

    model = VariationalCycleGAN(num_features = 23)
    print('Graph Compile Successful.')


text:
# -*- coding: utf-8 -*-

import os
import numpy as np

from PIL import Image
from PIL import ImageFile
from torchvision import transforms
import matplotlib.pyplot as plt

import torch
import torch.utils.data as data

ImageFile.LOAD_TRUNCATED_IMAGES = True

def get_loader(transform,
               mode='train',
               batch_size=1,
               img_folder="/",
               label_file=[],
               num_workers=0,
               fourier=0):
    """
    Args:
        transform : Image tranform
        mode: One of 'train' or 'test', this arguments are not used at here
        batch_size: Batch Size
        img_folder: the folder path that images are stored
        label_file: The file that stores annotation about the image and label
        num_workers: Number of Subprocesses to use data loading

    [Detailed explantion about the label file]
    =========================================================================
    Annotation label_file data structure

    list(dict)
    [{'image': image_id(str), 'label', label(int)}]

    Ex: [{'image': 'im_1' 'label': 1},
         {'image': 'im_2' 'label': 3},
         {'image': 'im_3' 'label': 7}]

    Image id is another expression of image path.
    The loader will load the image by it's path.

    label number starts from 0.
    For this label it has 18 different categories.
    So it will start at 0 and end at 17.
    =========================================================================

    There is no need to distinguish the train and the test.
    If your going to split to train and test image,

    Split the label_file to for example 8:2
    """

    assert len(label_file) != 0, "Please fill the label file QUARK!"
    if fourier == 0:
        dataset = aesthetics_dataset(transform=transform,
                                  mode=mode,
                                  batch_size=batch_size,
                                  label_file=label_file,
                                  img_folder=img_folder)
    elif fourier == 1:
        dataset = aesthetics_dataset_f(transform=transform,
                                  mode=mode,
                                  batch_size=batch_size,
                                  label_file=label_file,
                                  img_folder=img_folder)
    elif fourier == 2:
        dataset = aesthetics_dataset_f2(transform=transform,
                                  mode=mode,
                                  batch_size=batch_size,
                                  label_file=label_file,
                                  img_folder=img_folder)

    data_loader = data.DataLoader(dataset=dataset,
                                  batch_size=dataset.batch_size,
                                  shuffle=True,
                                  num_workers=num_workers)

    return data_loader

class aesthetics_dataset(data.Dataset):
    '''
    Dataset to train the clothes category classification model

    The __getitem__ will return 2 variables label and image.

    image = numpy array with size (3, 224, 224)
    label = Number label file with size (1)
    '''
    def __init__(self, transform, mode, batch_size, label_file, img_folder):
        self.transform = transform
        self.mode = mode
        self.batch_size = batch_size
        self.label_file = label_file
        self.img_folder = img_folder

    def __getitem__(self, index):
        label_num = self.label_file[index]['label']
        label = torch.tensor(label_num)

        #open image
        img_path = os.path.join(self.img_folder,
                                self.label_file[index]['image'])
        PIL_image = Image.open(img_path).convert('RGB')
        orig_image = PIL_image
        image = self.transform(orig_image)

        return image, label

    def __len__(self):
        return len(self.label_file)

class aesthetics_dataset_f(data.Dataset):
    '''
    Dataset to train the clothes category classification model

    The __getitem__ will return 2 variables label and image.

    image = numpy array with size (3, 224, 224)
    label = Number label file with size (1)
    '''
    def __init__(self, transform, mode, batch_size, label_file, img_folder):
        self.transform = transform
        self.mode = mode
        self.batch_size = batch_size
        self.label_file = label_file
        self.img_folder = img_folder
        self.transform_1 = transforms.Compose([
                transforms.Resize(256),
                transforms.RandomCrop(224)])
        self.transform_2 = transforms.Normalize((0.485, 0.456, 0.406),
                                                (0.229, 0.224, 0.225))

    def __getitem__(self, index):
        label_num = self.label_file[index]['label']
        label = torch.tensor(label_num)

        #open image
        img_path = os.path.join(self.img_folder,
                                self.label_file[index]['image'])
        PIL_image = Image.open(img_path).convert('RGB')
        PIL_image = self.transform_1(PIL_image)
        grey_im = PIL_image.convert('L')
        im = np.array(grey_im)

        # Fast Fourier Transform
        im = np.fft.fft2(im)
        im = np.fft.fftshift(im)
        im = np.log(np.abs(im+1))
        im = im / 255.0

        im = np.array(im, np.float32, copy=False)
        im = torch.from_numpy(im)
        im = im.unsqueeze(0)
        image = self.transform_2(im)

        return image, label

    def __len__(self):
        return len(self.label_file)

class aesthetics_dataset_f2(data.Dataset):
    '''
    Dataset to train the clothes category classification model

    The __getitem__ will return 2 variables label and image.

    image = numpy array with size (3, 224, 224)
    label = Number label file with size (1)
    '''
    def __init__(self, transform, mode, batch_size, label_file, img_folder):
        self.transform = transform
        self.mode = mode
        self.batch_size = batch_size
        self.label_file = label_file
        self.img_folder = img_folder
        self.transform_1 = transforms.Compose([
                transforms.Resize(256),
                transforms.RandomCrop(224)])
        self.transform_2 = transforms.Compose([
                transforms.ToTensor(),
                transforms.Normalize((0.485, 0.456, 0.406),
                                     (0.229, 0.224, 0.225))])
        self.transform_3 = transforms.Normalize((0.485, 0.456, 0.406),
                                                (0.229, 0.224, 0.225))

    def __getitem__(self, index):
        label_num = self.label_file[index]['label']
        label = torch.tensor(label_num)

        #open image
        img_path = os.path.join(self.img_folder,
                                self.label_file[index]['image'])
        PIL_image = Image.open(img_path).convert('RGB')
        PIL_image = self.transform_1(PIL_image)
        orig_im = np.array(PIL_image)
        image = self.transform_2(PIL_image)

        grey_im = PIL_image.convert('L')
        im = np.array(grey_im)

        # Fast Fourier Transform
        im = np.fft.fft2(im)
        im = np.fft.fftshift(im)
        im = np.log(np.abs(im+1))
        im = im / 255.0
        f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
        print("This image is : ", "GOOD" if label_num == 1 else "BAD")
        ax1.imshow(orig_im)
        ax2.imshow(im, cmap='gray')

        im = np.array(im, np.float32, copy=False)
        im = torch.from_numpy(im)
        im = im.unsqueeze(0)
        image_f = self.transform_3(im)

        return image, image_f, label

    def __len__(self):
        return len(self.label_file)

ext: py

max_stars_repo_name: IllgamhoDuck/ResTR

max_issues_repo_name: IllgamhoDuck/ResTR

max_forks_repo_name: IllgamhoDuck/ResTR

reason: import numpy


text:
import torch
import gpytorch
import time

import numpy as np

from meta_learn.models import AffineTransformedDistribution, EqualWeightedMixtureDist
from meta_learn.random_gp import RandomGPMeta
from meta_learn.util import _handle_input_dimensionality, DummyLRScheduler
from meta_learn.svgd import SVGD, RBF_Kernel, IMQSteinKernel
from meta_learn.abstract import RegressionModelMetaLearned
from config import device

class GPRegressionMetaLearnedSVGD(RegressionModelMetaLearned):

    def __init__(self, meta_train_data, num_iter_fit=10000, feature_dim=1,
                 prior_factor=0.01, weight_prior_std=0.5, bias_prior_std=3.0,
                 covar_module='NN', mean_module='NN', mean_nn_layers=(32, 32), kernel_nn_layers=(32, 32),
                 optimizer='Adam', lr=1e-3, lr_decay=1.0, kernel='RBF', bandwidth=None, num_particles=10,
                 task_batch_size=-1, normalize_data=True, random_seed=None):
        """
        PACOH-SVGD: Stein Variational Gradient Descent on PAC-optimal hyper-posterior.
        Meta-learns a set of GP priors (i.e. mean and kernel function)

        Args:
            meta_train_data: list of tuples of ndarrays[(train_x_1, train_t_1), ..., (train_x_n, train_t_n)]
            num_iter_fit: (int) number of gradient steps for fitting the parameters
            feature_dim: (int) output dimensionality of NN feature map for kernel function
            prior_factor: (float) weighting of the hyper-prior (--> meta-regularization parameter)
            weight_prior_std (float): std of Gaussian hyper-prior on weights
            bias_prior_std (float): std of Gaussian hyper-prior on biases
            covar_module: (gpytorch.mean.Kernel) optional kernel module, default: RBF kernel
            mean_module: (gpytorch.mean.Mean) optional mean module, default: ZeroMean
            mean_nn_layers: (tuple) hidden layer sizes of mean NN
            kernel_nn_layers: (tuple) hidden layer sizes of kernel NN
            optimizer: (str) type of optimizer to use - must be either 'Adam' or 'SGD'
            lr: (float) learning rate for prior parameters
            lr_decay: (float) lr rate decay multiplier applied after every 1000 steps
            kernel (std): SVGD kernel, either 'RBF' or 'IMQ'
            bandwidth (float): bandwidth of kernel, if None the bandwidth is chosen via heuristic
            num_particles: (int) number particles to approximate the hyper-posterior
            task_batch_size: (int) mini-batch size of tasks for estimating gradients
            normalize_data: (bool) whether the data should be normalized
            random_seed: (int) seed for pytorch
        """
        super().__init__(normalize_data, random_seed)

        assert mean_module in ['NN', 'constant', 'zero'] or isinstance(mean_module, gpytorch.means.Mean)
        assert covar_module in ['NN', 'SE'] or isinstance(covar_module, gpytorch.kernels.Kernel)
        assert optimizer in ['Adam', 'SGD']

        self.num_iter_fit, self.prior_factor, self.feature_dim = num_iter_fit, prior_factor, feature_dim
        self.weight_prior_std, self.bias_prior_std = weight_prior_std, bias_prior_std
        self.num_particles = num_particles
        if task_batch_size < 1:
            self.task_batch_size = len(meta_train_data)
        else:
            self.task_batch_size = min(task_batch_size, len(meta_train_data))

        # Check that data all has the same size
        self._check_meta_data_shapes(meta_train_data)
        self._compute_normalization_stats(meta_train_data)

        """ --- Setup model & inference --- """
        self._setup_model_inference(mean_module, covar_module, mean_nn_layers, kernel_nn_layers,
                                    kernel, bandwidth, optimizer, lr, lr_decay)

        # Setup components that are different across tasks
        self.task_dicts = []

        for train_x, train_y in meta_train_data:
            task_dict = {}

            # a) prepare data
            x_tensor, y_tensor = self._prepare_data_per_task(train_x, train_y)
            task_dict['train_x'], task_dict['train_y'] = x_tensor, y_tensor
            self.task_dicts.append(task_dict)

        self.fitted = False

    def meta_fit(self, valid_tuples=None, verbose=True, log_period=500, n_iter=None):
        """
        fits the hyper-posterior particles with SVGD

        Args:
            valid_tuples: list of valid tuples, i.e. [(test_context_x_1, test_context_t_1, test_x_1, test_t_1), ...]
            verbose: (boolean) whether to print training progress
            log_period (int) number of steps after which to print stats
            n_iter: (int) number of gradient descent iterations
        """

        assert (valid_tuples is None) or (all([len(valid_tuple) == 4 for valid_tuple in valid_tuples]))

        t = time.time()

        if n_iter is None:
            n_iter = self.num_iter_fit

        for itr in range(1, n_iter + 1):

            task_dict_batch = self.rds_numpy.choice(self.task_dicts, size=self.task_batch_size)
            self.svgd_step(task_dict_batch)
            self.lr_scheduler.step()

            # print training stats stats
            if itr == 1 or itr % log_period == 0:
                duration = time.time() - t
                t = time.time()

                message = 'Iter %d/%d - Time %.2f sec' % (itr, self.num_iter_fit, duration)

                # if validation data is provided  -> compute the valid log-likelihood
                if valid_tuples is not None:
                    valid_ll, valid_rmse, calibr_err = self.eval_datasets(valid_tuples)
                    message += ' - Valid-LL: %.3f - Valid-RMSE: %.3f - Calib-Err %.3f' % (valid_ll, valid_rmse, calibr_err)

                if verbose:
                    self.logger.info(message)

        self.fitted = True

    def predict(self, context_x, context_y, test_x, return_density=False):
        """
        Performs posterior inference (target training) with (context_x, context_y) as training data and then
        computes the predictive distribution of the targets p(y|test_x, test_context_x, context_y) in the test points

        Args:
            context_x: (ndarray) context input data for which to compute the posterior
            context_y: (ndarray) context targets for which to compute the posterior
            test_x: (ndarray) query input data of shape (n_samples, ndim_x)
            return_density: (bool) whether to return result as mean and std ndarray or as MultivariateNormal pytorch object

        Returns:
            (pred_mean, pred_std) predicted mean and standard deviation corresponding to p(t|test_x, test_context_x, context_y)
        """

        context_x, context_y = _handle_input_dimensionality(context_x, context_y)
        test_x = _handle_input_dimensionality(test_x)
        assert test_x.shape[1] == context_x.shape[1]

        # normalize data and convert to tensor
        context_x, context_y = self._prepare_data_per_task(context_x, context_y)

        test_x = self._normalize_data(X=test_x, Y=None)
        test_x = torch.from_numpy(test_x).float().to(device)

        with torch.no_grad():
            pred_dist = self.get_pred_dist(context_x, context_y, test_x)
            pred_dist = AffineTransformedDistribution(pred_dist, normalization_mean=self.y_mean,
                                                      normalization_std=self.y_std)
            pred_dist = EqualWeightedMixtureDist(pred_dist, batched=True)

            if return_density:
                return pred_dist
            else:
                pred_mean = pred_dist.mean.cpu().numpy()
                pred_std = pred_dist.stddev.cpu().numpy()
                return pred_mean, pred_std

    def _setup_model_inference(self, mean_module_str, covar_module_str, mean_nn_layers, kernel_nn_layers,
                               kernel, bandwidth, optimizer, lr, lr_decay):
        assert mean_module_str in ['NN', 'constant']
        assert covar_module_str in ['NN', 'SE']

        """ random gp model """
        self.random_gp = RandomGPMeta(size_in=self.input_dim, prior_factor=self.prior_factor,
                                  weight_prior_std=self.weight_prior_std, bias_prior_std=self.bias_prior_std,
                                  covar_module_str=covar_module_str, mean_module_str=mean_module_str,
                                  mean_nn_layers=mean_nn_layers, kernel_nn_layers=kernel_nn_layers)

        """ Setup SVGD inference"""

        if kernel == 'RBF':
            kernel = RBF_Kernel(bandwidth=bandwidth)
        elif kernel == 'IMQ':
            kernel = IMQSteinKernel(bandwidth=bandwidth)
        else:
            raise NotImplemented

        # sample initial particle locations from prior
        self.particles = self.random_gp.sample_params_from_prior(shape=(self.num_particles,))

        self._setup_optimizer(optimizer, lr, lr_decay)

        self.svgd = SVGD(self.random_gp, kernel, optimizer=self.optimizer)

        """ define svgd step """

        def svgd_step(tasks_dicts):
            # tile data to svi_batch_shape
            train_data_tuples_tiled = []
            for task_dict in tasks_dicts:
                x_data, y_data = task_dict['train_x'], task_dict['train_y']
                x_data = x_data.view(torch.Size((1,)) + x_data.shape).repeat(self.num_particles, 1, 1)
                y_data = y_data.view(torch.Size((1,)) + y_data.shape).repeat(self.num_particles, 1)
                train_data_tuples_tiled.append((x_data, y_data))

            self.svgd.step(self.particles, train_data_tuples_tiled)

        """ define predictive dist """

        def get_pred_dist(x_context, y_context, x_valid):
            with torch.no_grad():
                x_context = x_context.view(torch.Size((1,)) + x_context.shape).repeat(self.num_particles, 1, 1)
                y_context = y_context.view(torch.Size((1,)) + y_context.shape).repeat(self.num_particles, 1)
                x_valid = x_valid.view(torch.Size((1,)) + x_valid.shape).repeat(self.num_particles, 1, 1)

                gp_fn = self.random_gp.get_forward_fn(self.particles)
                gp, likelihood = gp_fn(x_context, y_context, train=False)
                pred_dist = likelihood(gp(x_valid))
            return pred_dist

        self.svgd_step = svgd_step
        self.get_pred_dist = get_pred_dist

    def _setup_optimizer(self, optimizer, lr, lr_decay):
        assert hasattr(self, 'particles'), "SVGD must be initialized before setting up optimizer"

        if optimizer == 'Adam':
            self.optimizer = torch.optim.Adam([self.particles], lr=lr)
        elif optimizer == 'SGD':
            self.optimizer = torch.optim.SGD([self.particles], lr=lr)
        else:
            raise NotImplementedError('Optimizer must be Adam or SGD')

        if lr_decay < 1.0:
            self.lr_scheduler = torch.optim.lr_scheduler.StepLR(self.optimizer, 1000, gamma=lr_decay)
        else:
            self.lr_scheduler = DummyLRScheduler()

    def _vectorize_pred_dist(self, pred_dist):
        multiv_normal_batched = pred_dist.dists
        normal_batched = torch.distributions.Normal(multiv_normal_batched.mean, multiv_normal_batched.stddev)
        return EqualWeightedMixtureDist(normal_batched, batched=True, num_dists=multiv_normal_batched.batch_shape[0])

if __name__ == "__main__":
    """ 1) Generate some training data from GP prior """
    from experiments.data_sim import GPFunctionsDataset

    data_sim = GPFunctionsDataset(random_state=np.random.RandomState(26))

    # meta_train_data = data_sim.generate_meta_train_data(n_tasks=5, n_samples=40)
    meta_test_data = data_sim.generate_meta_test_data(n_tasks=5, n_samples_context=40, n_samples_test=160)
    meta_train_data = [(context_x, context_y) for context_x, context_y, _, _ in meta_test_data]

    NN_LAYERS = (16, 16)

    plot = False
    from matplotlib import pyplot as plt

    if plot:
        for x_train, y_train in meta_train_data:
            plt.scatter(x_train, y_train)
        plt.title('sample from the GP prior')
        plt.show()

    for prior_factor in [1e-3]:
        gp_model = GPRegressionMetaLearnedSVGD(meta_train_data, num_iter_fit=2000, prior_factor=prior_factor, num_particles=10,
                                             covar_module='SE', mean_module='NN', mean_nn_layers=NN_LAYERS, kernel_nn_layers=NN_LAYERS,
                                             bandwidth=0.5, task_batch_size=2)

        for i in range(10):
            itrs = 0
            gp_model.meta_fit(valid_tuples=meta_test_data, log_period=500, n_iter=1000)
            itrs += 10

            x_test = np.linspace(-5, 5, num=150)
            x_context, t_context, _, _ = meta_test_data[0]
            pred_mean, pred_std = gp_model.predict(x_context, t_context, x_test)
            ucb, lcb = gp_model.confidence_intervals(x_context, t_context, x_test, confidence=0.9)

            plt.scatter(x_context, t_context)
            plt.plot(x_test, pred_mean)
            plt.fill_between(x_test, lcb, ucb, alpha=0.2)
            plt.title('GPR meta SVGD (prior-factor =  %.4f) itrs = %i'%(prior_factor, itrs))
            plt.show()


text:
"""
Runs frozen vortex wake solver // lifting line theory
"""

import numpy as np
from scipy.interpolate import interp1d
from create_geometry import BladeGeometry

class LiftingLineSolver:

    def __init__(self, geo, r_rotor, weight=0.3, tol=1e-6, n_iter=1000):
        """
        :param geo: BladeGeometry class
        :param r_rotor: rotor radius [m]
        :param weight: weighting factor for next step
        :param tol: stopping criteria
        :param n_iter: number of iterations
        """

        # rotor discretization
        self.geo = geo
        # rotor properties
        self.u_rot = np.array([self.geo.tsr / r_rotor * geo.v_inf, 0, 0])
        self.u_inf = geo.v_inf
        self.r_rotor = r_rotor
        self.double_rotor = self.geo.double_rotor

        # solver settings
        self.weight = weight
        self.tol = tol
        self.n_iter = n_iter

        # airfoil variables
        self._polarAirfoil()  # computes minimum/maximum alpha, polynomials for CL and CD alpha approximations

    def _compute_circ(self, gamma, weight):
        self.geo.filaments[-1] = self.geo.filaments[-1] * (1 - weight) + (weight * gamma)

    def _velocity_3D_from_vortex_filament(self, cp_i, core):
        """
        Computes the induced velocity at a control point due to ALL rings (vectorized part)
        :param cp_i: single control point
        :param core: minimum value for circulation for stability
        :return: matrix of size (3, (n_blades x (n_span-1))) containing the induced velocities on the control point
        in u, v, w.
        """

        r_gamma = self.geo.filaments[-1]  # vortex strength of each filament
        x1 = self.geo.filaments[0]
        y1 = self.geo.filaments[1]
        z1 = self.geo.filaments[2]
        x2 = self.geo.filaments[3]
        y2 = self.geo.filaments[4]
        z2 = self.geo.filaments[5]
        xc, yc, zc = cp_i

        R1 = np.sqrt((xc - x1) ** 2 + (yc - y1) ** 2 + (zc - z1) ** 2)
        R2 = np.sqrt((xc - x2) ** 2 + (yc - y2) ** 2 + (zc - z2) ** 2)

        R12_xx = ((yc - y1) * (zc - z2)) - ((zc - z1) * (yc - y2))
        R12_xy = -((xc - x1) * (zc - z2)) + ((zc - z1) * (xc - x2))
        R12_xz = ((xc - x1) * (yc - y2)) - ((yc - y1) * (xc - x2))

        R12_sq = (R12_xx ** 2) + (R12_xy ** 2) + (R12_xz ** 2)

        R01 = ((x2 - x1) * (xc - x1)) + ((y2 - y1) * (yc - y1)) + ((z2 - z1) * (zc - z1))
        R02 = ((x2 - x1) * (xc - x2)) + ((y2 - y1) * (yc - y2)) + ((z2 - z1) * (zc - z2))

        # check if target point is in the vortex filament core,
        # and modify to solid body rotation
        R12_sq = np.where(R12_sq < core ** 2, core ** 2, R12_sq)
        R1 = np.where(R1 < core, core, R1)
        R2 = np.where(R1 < core, core, R2)

        K = (r_gamma / (4 * np.pi * R12_sq)) * ((R01 / R1) - (R02 / R2))

        U = np.sum(K * R12_xx, axis=1)
        V = np.sum(K * R12_xy, axis=1)
        W = np.sum(K * R12_xz, axis=1)

        return np.array([U, V, W])

    def _compute_induced_velocity(self):
        # create u, v, w matrix to store induced velocities. The matrix consists of elements in which the rows
        # represent the control points and columns the net effect of a single "ring".
        uvw_mat = np.zeros((3, self.geo.cp.shape[0], self.geo.filaments.shape[1]))  # three square matrices
        core = 0.00001  # no idea what this is
        for i, cp_i in enumerate(self.geo.cp[:, :3]):  # only loop over the coordinates of the control points
            temp_v = self._velocity_3D_from_vortex_filament(cp_i, core)
            uvw_mat[:, i, :] = temp_v

        return uvw_mat

    @staticmethod
    def _geo_blade(r_R):
        pitch = 2  # in deg
        chord = 3 * (1 - r_R) + 1
        twist = -14 * (1 - r_R)  # in deg
        result = [chord.flatten(), np.radians(twist + pitch).flatten()]

        return result

    def _polarAirfoil(self):
        data = np.loadtxt("polar_DU95W180.csv", delimiter=';')
        data = data[3:]
        data = np.array(data, dtype=float)

        alphaRad = np.radians(data[:, 0])
        self.amax = max(alphaRad)
        self.amin = min(alphaRad)
        self.fcl = interp1d(alphaRad, data[:, 1], fill_value=(data[0, 1], data[-1, 1]), bounds_error=False,
                            kind='cubic')
        self.fcd = interp1d(alphaRad, data[:, 2], fill_value=(data[0, 2], data[-1, 2]), bounds_error=False,
                            kind='cubic')

    def _compute_loads_blade(self, v_norm, v_tan, r_R):
        V_mag2 = (v_norm ** 2 + v_tan ** 2)  # Velocity magnitude squared
        phi = np.arctan(v_norm / v_tan)  # Inflow angle

        # Get chord and twist
        [chord, twist] = self._geo_blade(r_R)

        alpha = twist + phi

        cl = self.fcl(alpha)
        cd = self.fcd(alpha)

        L = 0.5 * V_mag2 * cl * chord
        D = 0.5 * V_mag2 * cd * chord

        F_norm = L * np.cos(phi) + D * np.sin(phi)
        F_tan = L * np.sin(phi) - D * np.cos(phi)

        Gamma = 0.5 * np.sqrt(V_mag2) * cl * chord

        return np.array([F_norm, F_tan, Gamma, alpha, phi])

    def _initialize_solver(self):

        # update Gamma given Gamma matrix, weight, and new Gamma
        self._compute_circ(gamma=1, weight=1)  # updates self.geo itself

        # compute [ui, vi, wi] based on vortex strength and distance
        # between control point and vortex
        v_induced = self._compute_induced_velocity()

        return v_induced

    def run_solver(self):
        # determine radial position of control point
        pos_radial = np.sqrt(np.sum(self.geo.cp[:(self.geo.n_blades * (self.geo.n_span - 1)), :3] ** 2,
                                    axis=1)).reshape(-1, 1)
        r_R = pos_radial / self.r_rotor

        if self.double_rotor:  # copy r/R for second rotor
            pos_radial = np.tile(pos_radial, 2).reshape((-1, 1), order='F')
            cp = self.geo.cp[:, :3].copy()
            cp[int(len(cp) / 2):, :3] = self.geo._compute_cp(self.geo.phase_diff)[:, :3]

        # initialize gamma vectors new and old
        gamma_new = np.ones((len(self.geo.cp), 1))
        # initialize output variables
        a = np.ones(len(self.geo.cp)) * 0.33
        aline = np.ones(len(self.geo.cp))
        # r_R = np.ones(len(self.geo.cp))
        f_norm = np.ones(len(self.geo.cp))
        f_tan = np.ones(len(self.geo.cp))
        gamma = np.ones(len(self.geo.cp))
        alpha = np.ones(len(self.geo.cp))
        phi = np.ones(len(self.geo.cp))
        # initial error
        err = 1.0
        error_log = []

        for i in range(self.n_iter):

            # re-discretise wake sheet based on new induction factor
            self.geo.a = np.mean(a[(self.geo.n_span - 1):])  # take only the values of the first rotor

            if self.double_rotor:  # shift filament coords of the 1st rotor
                self.geo.doubleRotorUpdate()
            else:
                self.geo.filaments = self.geo.compute_ring()

            # compute system of linear eqn. (influence of each filament)
            uvw_mat = self._initialize_solver()
            # update circulation
            gamma_curr = gamma_new

            # calculate velocity, circulation, control points
            # directly compute total velocity at each control point by mat. vec. product
            u = uvw_mat[0] @ gamma_curr
            v = uvw_mat[1] @ gamma_curr
            w = uvw_mat[2] @ gamma_curr

            if self.double_rotor:
                # compute perceived velocity by blade element
                vel_rot = np.cross(-self.u_rot, cp)
                vel_per = np.hstack([self.u_inf + u + vel_rot[:, 0].reshape(-1, 1),
                                     v + vel_rot[:, 1].reshape(-1, 1),
                                     w + vel_rot[:, 2].reshape(-1, 1)])

                # calculate azimuthal and axial velocity
                inv_pos_radial = np.hstack([-1 / pos_radial, np.zeros(pos_radial.shape), np.zeros(pos_radial.shape)])
                azim_dir = np.cross(inv_pos_radial, cp)
                u_azim = np.array([azim @ vel for azim, vel in zip(azim_dir, vel_per)])
                u_axial = vel_per @ np.array([1, 0, 0])  # should be the same as [1, 0, 0] @ vel_per (dot product)
            else:
                # compute perceived velocity by blade element
                vel_rot = np.cross(-self.u_rot, self.geo.cp[:, :3])
                vel_per = np.hstack([self.u_inf + u + vel_rot[:, 0].reshape(-1, 1),
                                     v + vel_rot[:, 1].reshape(-1, 1),
                                     w + vel_rot[:, 2].reshape(-1, 1)])

                # calculate azimuthal and axial velocity
                inv_pos_radial = np.hstack([-1 / pos_radial, np.zeros(pos_radial.shape), np.zeros(pos_radial.shape)])
                azim_dir = np.cross(inv_pos_radial, self.geo.cp[:, :3])
                u_azim = np.array([azim @ vel for azim, vel in zip(azim_dir, vel_per)])
                u_axial = vel_per @ np.array([1, 0, 0])  # should be the same as [1, 0, 0] @ vel_per (dot product)

            # calculate loads using BEM
            blade_loads = self._compute_loads_blade(u_axial, u_azim, pos_radial / self.r_rotor)

            # update loads and circulation
            gamma_new = blade_loads[2].reshape(-1, 1)
            a = -(u + vel_rot[:, 0].reshape(-1, 1)) / self.u_inf
            aline = u_azim / (pos_radial.flatten() * self.u_rot[0]) - 1
            f_norm = blade_loads[0]
            f_tan = blade_loads[1]
            gamma = blade_loads[2]
            alpha = blade_loads[3]
            phi = blade_loads[4]

            # check convergence
            err_ref = max(0.001, np.max(np.abs(gamma_new)))  # choose highest value for reference error
            err = np.max(np.abs(gamma_new - gamma_curr)) / err_ref
            error_log.append(err)
            print("iteration: {} | current error: {}".format(i, err))
            if err < self.tol:
                print("solution converged")
                break

            # set new estimate of bound circulation
            gamma_new = (1 - self.weight) * gamma_curr + self.weight * gamma_new
        print("solution unconverged error: {}".format(err)) if err > self.tol else None
        return [a, aline, r_R, f_norm, f_tan, gamma, alpha, phi, error_log]

    def CP_and_CT(self, a, aline, r_R, f_norm, f_tan, v_inf, omega, radius, nblades):

        CT_LLM = []
        CP_LLM = []
        CP_flow = []

        for i in range(len(r_R) - 1):
            r_R_temp = (r_R[i] + r_R[i + 1]) / 2
            drtemp = (-r_R[i] + r_R[i + 1])
            # Prof. Ferreira
            CT_LLM.append((drtemp * f_norm[i] * nblades) / (0.5 * (v_inf ** 2) * np.pi * radius))
            CP_LLM.append((drtemp * f_tan[i] * r_R_temp * omega * nblades) / (0.5 * (v_inf ** 3) * np.pi))
            # another method
            U_tan = r_R_temp * radius * omega * aline[i]
            U_norm = v_inf * (1 - a[i])
            CP_flow.append((drtemp * (f_norm[i] * U_norm - f_tan[i] * U_tan) * nblades) / (0.5 * (v_inf ** 3) *
                                                                                np.pi * radius))

        # Based on induction factor
        CT_LLM2 = 4 * np.mean(a) * (1 - np.mean(a))
        CP_LLM2 = 4 * np.mean(a) * (1 - np.mean(a))**2

        return [CP_LLM, CT_LLM, CP_LLM2, CT_LLM2]

reason: import numpy,from scipy


text:
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import cPickle

from keras.layers import Input, Dense, Lambda, Flatten, Reshape, Layer
from keras.layers import Conv2D, Conv2DTranspose
from keras.models import Model
from keras import backend as K
from keras import metrics
from keras import optimizers

from full_params_conv_101 import *

if K.image_data_format() == 'channels_first':
    original_img_size = (img_chns, img_rows, img_cols)
else:
    original_img_size = (img_rows, img_cols, img_chns)

x = Input(shape=original_img_size)
conv_1 = Conv2D(32,
                kernel_size=(4, 4),
                strides=(2, 2),
                padding='same', activation='relu')(x)
conv_2 = Conv2D(64,
                kernel_size=(4, 4),
                padding='same', activation='relu',
                strides=(2, 2))(conv_1)
conv_3 = Conv2D(128,
                kernel_size=(4, 4),
                padding='same', activation='relu',
                strides=(2, 2))(conv_2)
conv_4 = Conv2D(256,
                kernel_size=(4, 4),
                padding='same', activation='relu',
                strides=(2, 2))(conv_3)
flat = Flatten()(conv_4)
hidden = Dense(intermediate_dim, activation='relu')(flat)

z_mean = Dense(latent_dim)(hidden)
z_log_var = Dense(latent_dim)(hidden)

def sampling(args):
    z_mean, z_log_var = args
    epsilon = K.random_normal(shape=(K.shape(z_mean)[0], latent_dim),
                              mean=0., stddev=epsilon_std)
    return z_mean + K.exp(z_log_var) * epsilon

# note that "output_shape" isn't necessary with the TensorFlow backend
# so you could write `Lambda(sampling)([z_mean, z_log_var])`
z = Lambda(sampling, output_shape=(latent_dim,))([z_mean, z_log_var])

# we instantiate these layers separately so as to reuse them later
decoder_hid = Dense(intermediate_dim, activation='relu')
decoder_upsample = Dense(16384, activation='relu')

if K.image_data_format() == 'channels_first':
    output_shape = (batch_size, 256, 8, 8)
else:
    output_shape = (batch_size, 8, 8, 256)

decoder_reshape = Reshape(output_shape[1:])
decoder_deconv_1 = Conv2DTranspose(128,
                                   kernel_size=(4, 4),
                                   padding='same',
                                   strides=(2, 2),
                                   activation='relu')
decoder_deconv_2 = Conv2DTranspose(64,
                                   kernel_size=(4, 4),
                                   padding='same',
                                   strides=(2, 2),
                                   activation='relu')
decoder_deconv_3_upsamp = Conv2DTranspose(32,
                                          kernel_size=(4, 4),
                                          strides=(2, 2),
                                          padding='same',
                                          activation='relu')
decoder_mean_squash = Conv2DTranspose(3,
                             kernel_size=(4, 4),
                             strides=(2, 2),
                             padding='same',
                             activation='relu')

hid_decoded = decoder_hid(z)
up_decoded = decoder_upsample(hid_decoded)
reshape_decoded = decoder_reshape(up_decoded)
deconv_1_decoded = decoder_deconv_1(reshape_decoded)
deconv_2_decoded = decoder_deconv_2(deconv_1_decoded)
x_decoded_relu = decoder_deconv_3_upsamp(deconv_2_decoded)
x_decoded_mean_squash = decoder_mean_squash(x_decoded_relu)

# Custom loss layer
class CustomVariationalLayer(Layer):
    def __init__(self, **kwargs):
        self.is_placeholder = True
        super(CustomVariationalLayer, self).__init__(**kwargs)

    def vae_loss(self, x, x_decoded_mean_squash):
        x = K.flatten(x)
        x_decoded_mean_squash = K.flatten(x_decoded_mean_squash)
        xent_loss = img_rows * img_cols * metrics.binary_crossentropy(x, x_decoded_mean_squash)
        kl_loss = - 0.5 * K.mean(1 + z_log_var - K.square(z_mean) - K.exp(z_log_var), axis=-1)
        return K.mean(xent_loss + kl_loss)

    def call(self, inputs):
        x = inputs[0]
        x_decoded_mean_squash = inputs[1]
        loss = self.vae_loss(x, x_decoded_mean_squash)
        self.add_loss(loss, inputs=inputs)
        # We don't use this output.
        return x

y = CustomVariationalLayer()([x, x_decoded_mean_squash])
vae = Model(x, y)

sgd = optimizers.SGD(lr=0.01)
vae.compile(optimizer=sgd, loss=None)
vae.summary()

"""
with open('../datasets/101_ObjectCategories.pkl') as f:
    dic = cPickle.load(f)
    x_train = dic['all_images']
"""
x_train = np.load('../datasets/full_x.npy')
print "dataset loaded"

history = vae.fit(x_train,
        shuffle=True,
        epochs=epochs,
        batch_size=batch_size,
        )

# build a model to project inputs on the latent space
encoder = Model(x, z_mean)

"""
# display a 2D plot of the digit classes in the latent space
x_test_encoded = encoder.predict(x_test, batch_size=batch_size)
plt.figure(figsize=(6, 6))
plt.scatter(x_test_encoded[:, 0], x_test_encoded[:, 1], c=y_test)
plt.colorbar()
plt.show()
"""

# build a digit generator that can sample from the learned distribution
decoder_input = Input(shape=(latent_dim,))
_hid_decoded = decoder_hid(decoder_input)
_up_decoded = decoder_upsample(_hid_decoded)
_reshape_decoded = decoder_reshape(_up_decoded)
_deconv_1_decoded = decoder_deconv_1(_reshape_decoded)
_deconv_2_decoded = decoder_deconv_2(_deconv_1_decoded)
_x_decoded_relu = decoder_deconv_3_upsamp(_deconv_2_decoded)
_x_decoded_mean_squash = decoder_mean_squash(_x_decoded_relu)
generator = Model(decoder_input, _x_decoded_mean_squash)

vae.save('../models/object101_ld_%d_conv_%d_id_%d_e_%d_vae.h5' % (latent_dim, num_conv, intermediate_dim, epochs))
encoder.save('../models/object101_ld_%d_conv_%d_id_%d_e_%d_encoder.h5' % (latent_dim, num_conv, intermediate_dim, epochs))
generator.save('../models/object101_ld_%d_conv_%d_id_%d_e_%d_generator.h5' % (latent_dim, num_conv, intermediate_dim, epochs))

fname = '../models/object101_ld_%d_conv_%d_id_%d_e_%d_history.pkl' % (latent_dim, num_conv, intermediate_dim, epochs)
with open(fname, 'wb') as file_pi:
    cPickle.dump(history.history, file_pi)

"""
# display a 2D manifold of the digits
n = 15  # figure with 15x15 digits
digit_size = 28
figure = np.zeros((digit_size * n, digit_size * n))
# linearly spaced coordinates on the unit square were transformed through the inverse CDF (ppf) of the Gaussian
# to produce values of the latent variables z, since the prior of the latent space is Gaussian
grid_x = norm.ppf(np.linspace(0.05, 0.95, n))
grid_y = norm.ppf(np.linspace(0.05, 0.95, n))

for i, yi in enumerate(grid_x):
    for j, xi in enumerate(grid_y):
        z_sample = np.array([[xi, yi]])
        z_sample = np.tile(z_sample, batch_size).reshape(batch_size, 2)
        x_decoded = generator.predict(z_sample, batch_size=batch_size)
        digit = x_decoded[0].reshape(digit_size, digit_size)
        figure[i * digit_size: (i + 1) * digit_size,
               j * digit_size: (j + 1) * digit_size] = digit

plt.figure(figsize=(10, 10))
plt.imshow(figure, cmap='Greys_r')
plt.show()
"""


text:
"""
Perona Malik Diffusion
"""

# Author: Khoa Nguyen <khoa.v18nguyen@gmail.com>

import copy
import numpy as np

def _div(v, K, esf):
    """Divergence operator"""

    grad_x = np.gradient(v)[0]
    grad_y = np.gradient(v)[1]

    norm2 = grad_x**2 + grad_y**2

    if esf == 'c1':
        factor = 1 - np.exp(-norm2/K**2)
    elif esf == 'c2':
        factor = 1/np.sqrt(1 + norm2/K**2)
    elif esf == 'c3':
        factor = 1/(1 + norm2/K**2)
    else:
        raise ValueError(' Edge Stopping Function Errors !')

    return np.gradient(factor*grad_x)[0] + np.gradient(factor*grad_y)[1]

def _pm(u, N=50, lambda_=0.07, K=0.2, esf='c1'):
    """

    :param u: ndarray
        observed image

    :param N: int
            number of iterations

    :param lambda_: float
            learning rate

    :param K: int
            edge stopping function's constant

    :param esf: string
        edge stopping function:
            + c1
            + c2
            + c3

    :return:
        uf: ndarray
            de-noised image result

    """
    uf = copy.deepcopy(u)  # initialization

    for k in range(N):
        uf = uf + lambda_*_div(uf, K, esf)

    return uf


text:
import scipy as sp
import numpy as np

def compute_perm_amg(parents):
    """
    Return a list of indices to reorder the adjacency and data matrices so
    that the union of two neighbors from layer to layer forms a binary tree.
    """

    # Order of last layer is random (chosen by the clustering algorithm).
    indices = []
    if len(parents) > 0:
        M_last = max(parents[-1]) + 1
        indices.append(list(range(M_last)))

    for parent in parents[::-1]:
        indices_layer = []
        for i in indices[-1]:
            indices_node = list(np.where(parent == i)[0])
            indices_layer.extend(indices_node)
        indices.append(indices_layer)

    return indices[::-1]

def perm_adjacency(A, indices):
    """
    Permute adjacency matrix, i.e. exchange node ids,
    so that binary unions form the clustering tree.
    """
    if indices is None:
        return A

    M, M = A.shape
    Mnew = len(indices)
    assert Mnew >= M
    A = A.tocoo()

    # Add Mnew - M isolated vertices.
    if Mnew > M:
        rows = sp.sparse.coo_matrix((Mnew-M,    M), dtype=np.float32)
        cols = sp.sparse.coo_matrix((Mnew, Mnew-M), dtype=np.float32)
        A = sp.sparse.vstack([A, rows])
        A = sp.sparse.hstack([A, cols])

    # Permute the rows and the columns.
    perm = np.argsort(indices)
    A.row = np.array(perm)[A.row]
    A.col = np.array(perm)[A.col]

    # assert np.abs(A - A.T).mean() < 1e-9
    assert type(A) is sp.sparse.coo.coo_matrix
    return A

max_stars_repo_name: zifeo/Titanic

max_issues_repo_name: zifeo/Titanic

max_forks_repo_name: zifeo/Titanic

reason: import numpy,import scipy


text:
import json
import math as m
import psycopg2
import sys
import numpy as np
from urllib.parse import urlparse
import socket
import asyncio
import ast
import config
from __init__ import app

paths = config.Path('matt')
sys.path.insert(0, paths.markov_chain())
sys.path.insert(0, paths.personal_info())
from hmm import *
from obs import *

CHUNK_TYPE_MOUSE = 0
CHUNK_TYPE_KEYBOARD = 1
WEBPAGE_TIME = 2
MAX_EVENTS_SIZE = 3
MAX_OBS_SIZE = 2
LEARN = 1
PREDICT = 2
LEARNED_ENOUGH = 50

class HMM(DB, hmm):
	def __init__(self, db, user, password, host, port, user_id):
		DB.__init__(self, db, user, password, host, port)
		self.get_connection()
		self.get_cursor()
		self.user_id = user_id

	def hmm_get_(self):
		SELECT_HMM = '''SELECT * FROM "hmm" WHERE "user_id" = {} '''.format(self.user_id)
		self.cursor.execute(SELECT_HMM)
		[(A, B, pi, self.status, tmp, ),] = self.cursor.fetchall()
		hmm.__init__(self, np.array(pi), np.array(A), np.array(B))

	def hmm_get_status(self):
		if self.status > LEARNED_ENOUGH:
			return PREDICT
		return LEARN

	def hmm_init_(self, unique_obs_len):
		a = (m.sqrt(unique_obs_len) * np.array([[0.5, 0.5] , [0.5, 0.5]])).tolist()
		r_1 = np.random.random(unique_obs_len)
		r_1 = (m.sqrt(unique_obs_len) * r_1 / np.sum(r_1)).tolist()
		r_2 = np.random.random(unique_obs_len)
		r_2 = (m.sqrt(unique_obs_len) * r_2 / sum(r_2)).tolist()
		b = np.array([r_1, r_2]).tolist()
		pi = (m.sqrt(unique_obs_len) * np.array([0.1, 0.9])).tolist()
		self.status = 0
		INSERT_HMM = '''INSERT INTO "hmm" (transition, emission, distribution, status, user_id) VALUES (array {}, array {}, array {}, 0, {});'''.format(a, b, pi, self.user_id)
		try:
			self.cursor.execute(INSERT_HMM)
			self.conn.commit()
		except:
			self.conn.rollback()
		hmm.__init__(self, np.array(pi), np.array(a), np.array(b))

	def hmm_update(self, obs_len):
		print("A", self.A, "B", self.B, "pi", self.pi)

		UPDATE_HMM = '''UPDATE "hmm" SET "transition" = array {}, "emission" = array {}, "distribution" = array {}, "status" = {} WHERE "user_id" = {}'''.format(self.A.tolist(), self.B.tolist(), self.pi.tolist(),  self.status + obs_len, self.user_id)
		try:
			self.cursor.execute(UPDATE_HMM)
			self.conn.commit()
		except:
			self.conn.rollback()
		self.hmm_get_()

	# def hmm_set(self, A = [[0.5, 0.5], [0.5, 0.5]], B = [[0.5, 0.5] for _ in range(100)], pi = [0.5, 0.5], obs_num = 0):
	# 	SET_HMM = '''UPDATE "hmm" SET "transition" = array {}, "emission" = array {}, "distribution" = array {}, "status" = {} WHERE "user_id" = {}'''.format(A, B, pi, self.status + obs_num, self.user_id)
	# 	try:
	# 		self.cursor.execute(SET_HMM)
	# 		self.conn.commit()
	# 	except:
	# 		self.conn.rollback()
	# 	self.hmm_get_()
	# 	print("A", self.A, "B", self.B, "pi", self.pi)

	def hmm_learn(self, obs):
		print(obs)
		self.hmm_init_(len(set(obs)))
		self.learn(obs)
		self.hmm_update(len(obs))
		print("updated")

	def hmm_predict(self, obs):
		print(obs)
		self.hmm_get_()
		path, prob = self.viterbi(obs)
		print(prob)

dictionary = {}

class ClientServerProtocol(asyncio.Protocol):
	def connection_made(self, transport):
		self.transport = transport

	def process_data(self, data):
		status, user_id, data = data[:-1].split(" ", 2)
		print(status, user_id)
		try:
			user_id = int(user_id)
			data = ast.literal_eval(data)
			if status not in ['learn', 'predict']:
				return "error\nwrong data\n\n"
		except:
			return "error\nwrong data\n\n"
		res = ""
		markov_model = HMM(app.config['DB'], app.config['USER'], app.config['PASSWORD'], app.config['HOST'], app.config['PORT'], user_id)
		print(data)
		if status == 'predict':
			markov_model.hmm_predict(data)
			print("PREDICT")
		if status == 'learn':
			markov_model.hmm_learn(data)
			print("LEARN")

		# markov_model.hmm_predict(data)
		res += "ok\n\n"
		return res

	def data_received(self, data):
		resp = self.process_data(data.decode())
		self.transport.write(resp.encode())

def run_server(host, port):
    loop = asyncio.get_event_loop()
    coro = loop.create_server(ClientServerProtocol, host, port)
    server = loop.run_until_complete(coro)

    try:
        loop.run_forever()
    except KeyboardInterrupt:
        pass

    server.close()
    loop.run_until_complete(server.wait_closed())
    loop.close()

run_server("127.0.0.1", 8181)


text:
import pandas as pd
import seaborn as sns
import numpy as np
import altair as alt
import matplotlib.pyplot as plt
import os
import sys
from pathlib import Path
from ocp_table_tpot.globals import Globals as gd
from tpot import TPOTRegressor
from src.models.model import HistoricalMedian,XGBoost,LinearModel,RF,KNN,SVM,mase,TimeSeriesSplitImproved,my_custom_scorer
from sklearn.metrics import make_scorer

root = Path(__file__).resolve().parents[2]
print(root)

def mase_error(y_true,y_pred):
    return mase(y_pred,y_true)

df_tsfresh = pd.read_pickle(f'{root}/data/processed/train_test_tsfresh.pkl').reset_index(level = 0)
data_dict = pd.read_pickle(f'{root}/data/processed/data_dict_all.pkl')

year = 2016
tgt = 'rougher.output.recovery'
X = data_dict[year]['X_train']
y = data_dict[year]['y_train'][tgt].dropna()

Nmonths_total = 8
Nspl = int(Nmonths_total * 30 / 15)
Nmonths_test = 4
Nmonths_min_train = 2.5
cv = TimeSeriesSplitImproved(n_splits=Nspl)

X = X.sample(frac=0.3).sort_index()
inds = X.index.intersection(y.index)

X = X.loc[inds]
y = y.loc[inds]

gens = 50
pops = 50
seed  = 123

sc= my_custom_scorer()

tpot = TPOTRegressor(generations=50,
                                     verbosity=2,
                                     cv=cv,
                                     scoring=gd.SCORING_LOOKUP_DICT['mae'],
                                     config_dict=gd.CONFIG_LOOKUP_DICT['rf_boost'],
                                     population_size=50,
                                     random_state=seed,
                                     periodic_checkpoint_folder='../results/',
                                     n_jobs=10,
                                     memory=None)
tpot.fit(X,y)
tpot.export(f'../results/tpot-pipeline_{year}_{tgt}.py')

max_stars_repo_name: AntonBiryukovUofC/unearthed-zinc

max_issues_repo_name: AntonBiryukovUofC/unearthed-zinc

max_forks_repo_name: AntonBiryukovUofC/unearthed-zinc


text:
import os
import cv2 as cv
import numpy as np

# people = ['Ben Affleck','Elton John','Jerry Seinfield','Madonna','Mindy Kaling']
people = []
for i in os.listdir('Faces/train'):         #append all the names in the folders into a list
    people.append(i)

DIR = r'/home/maitidevi/Documents/opencv-course/Resources/Faces/train'

haar_classifier = cv.CascadeClassifier('Face Detection/haar_cascade.xml')

features = []  #the image arrays of the faces
labels = [] #for every face in the features list, points to the corresponding people

def create_train():
    for person in people:
        path = os.path.join(DIR,person)
        label = people.index(person)

        for img in os.listdir(path):
            img_path = os.path.join(path,img)

            img_array = cv.imread(img_path)
            gray = cv.cvtColor(img_array,cv.COLOR_BGR2GRAY)

            faces_rect = haar_classifier.detectMultiScale(gray, scaleFactor=1.1, minNeighbors=4)

            for (x,y,w,h) in faces_rect:
                faces_roi = gray[y:y+h,x:x+h]   #regions of interest --> faces detection
                features.append(faces_roi)
                labels.append(label)

create_train()
print("~~~~~~~~~~~~~~Training done~~~~~~~~~~~~~")

features = np.array(features,dtype=np.object)
labels = np.array(labels)

face_recogniser = cv.face.LBPHFaceRecognizer_create()
face_recogniser.train(features,labels)

face_recogniser.save('Face Detection/face_trained.yml')   #train this same face recogniser model into another file/directory
# np.save('Face Detection/features.npy',features)     #save the features of the trained image
# np.save('Face Detection/labels.npy',labels)         #save the labels of the trained image

max_stars_repo_path: Face Detection/face_train.py

max_stars_repo_name: sanzeshstha/opencv

max_issues_repo_path: Face Detection/face_train.py

max_issues_repo_name: sanzeshstha/opencv

max_forks_repo_path: Face Detection/face_train.py

max_forks_repo_name: sanzeshstha/opencv


text:
"""Do CPT using Gibbs sampling.

Uses Gensim fucntionality.

Papers:
- Finding scientific topics
- A Theoretical and Practical Implementation Tutorial on Topic Modeling and
Gibbs Sampling
- Mining contrastive opinions on political texts using cross-perspective topic
model
"""

from __future__ import division
import numpy as np
import CPTCorpus
import glob
import logging
import time
import pandas as pd
import os

from gibbs_inner import gibbs_inner

logger = logging.getLogger(__name__)

class GibbsSampler():
    PARAMETER_DIR = '{}/parameter_samples'
    PHI_TOPIC = 'phi_topic'
    PHI_OPINION = 'phi_opinion_{}'
    THETA = 'theta'
    NKS = 'nks'

    def __init__(self, corpus, nTopics=10, alpha=0.02, beta=0.02, beta_o=0.02,
                 nIter=2, out_dir=None, sample_interval=10, initialize=True):
        self.corpus = corpus
        self.nTopics = nTopics
        self.alpha = alpha
        self.beta = beta
        self.nPerspectives = len(self.corpus.perspectives)
        self.beta_o = beta_o
        self.nIter = nIter
        self.sampleInterval = sample_interval

        self.VT = len(self.corpus.topicDictionary)
        self.VO = len(self.corpus.opinionDictionary)

        self.out_dir = out_dir
        if self.out_dir:
            if not os.path.exists(self.out_dir):
                os.makedirs(out_dir)
        logger.debug('output directory: {}'.format(self.out_dir))
        parameter_dir = self.get_parameter_dir_name()
        if not os.path.exists(parameter_dir):
            os.makedirs(parameter_dir)

        if initialize:
            self._initialize()

    def _initialize(self, phi_topic=None):
        """Initializes the Gibbs sampler."""
        logger.info('started initialization ({})'.format(str(self)))

        if not isinstance(phi_topic, np.ndarray):
            logger.debug('working with train set')
            self.documents = self.corpus
        else:
            logger.debug('working with test set')
            self.documents = self.corpus.testSet()

        self.DT = sum([len(p.corpus(phi_topic).topicCorpus)
                       for p in self.corpus.perspectives])
        self.DO = max([len(p.corpus(phi_topic).opinionCorpus)
                       for p in self.corpus.perspectives])
        self.maxDocLengthT = max([p.corpus(phi_topic).topicCorpus.maxDocLength
                                 for p in self.corpus.perspectives])
        self.maxDocLengthO = max([p.corpus(phi_topic).opinionCorpus.
                                  maxDocLength
                                  for p in self.corpus.perspectives])

        logger.debug('{} documents in the corpus (DT)'.format(self.DT))
        msg = '{} documents in the biggest perspective (DO)'.format(self.DO)
        logger.debug(msg)
        msg = 'maximum document lengths found: {} (topic) {} (opinion)'
        logger.debug(msg.format(self.maxDocLengthT, self.maxDocLengthO))

        # topics
        self.z = np.zeros((self.DT, self.maxDocLengthT), dtype=np.int)
        self.ndk = np.zeros((self.DT, self.nTopics), dtype=np.int)
        self.nkw = np.zeros((self.nTopics, self.VT), dtype=np.int)
        self.nk = np.zeros(self.nTopics, dtype=np.int)
        self.ntd = np.zeros(self.DT, dtype=np.int)

        # opinions
        self.x = np.zeros((self.nPerspectives, self.DO, self.maxDocLengthO),
                          dtype=np.int)
        self.nrs = np.zeros((self.nPerspectives, self.nTopics, self.VO),
                            dtype=np.int)
        self.ns = np.zeros((self.nPerspectives, self.nTopics), dtype=np.int)

        if phi_topic is not None:
            phi_topic_sum = np.sum(phi_topic, axis=0)

        logger.debug('Starting loop over documents/words')

        # loop over the words in the corpus
        for d, persp, d_p, doc in self.documents:
            for w_id, i in self.corpus.words_in_document(doc, 'topic'):
                #print d, persp, d_p, w_id, i
                if phi_topic is None:
                    topic = np.random.randint(0, self.nTopics)
                else:
                    p = phi_topic[:, w_id] / phi_topic_sum[w_id]
                    topic = self.sample_from(p)
                    #print len(phi_topic[:, w_id])
                    #print np.sum(phi_topic[:, w_id])
                #print topic
                self.z[d, i] = topic
                self.ndk[d, topic] += 1
                self.nkw[topic, w_id] += 1
                self.nk[topic] += 1
                self.ntd[d] += 1

            for w_id, i in self.corpus.words_in_document(doc, 'opinion'):
                opinion = np.random.randint(0, self.nTopics)
                self.x[persp][d_p, i] = opinion
                self.nrs[persp, opinion, w_id] += 1
                self.ns[persp, opinion] += 1
        logger.info('finished initialization ({})'.format(str(self)))

    def sample_from(self, p):
        """Sample (new) topic from multinomial distribution p.
        Returns a word's the topic index based on p_z.

        The searchsorted method is used instead of
        np.random.multinomial(1,p).argmax(), because despite normalizing the
        probabilities, sometimes the sum of the probabilities > 1.0, which
        causes the multinomial method to crash. This probably has to do with
        machine precision.
        """
        return np.searchsorted(np.cumsum(p), np.random.rand())

    def calc_theta_topic(self):
        """Calculate theta based on the current word/topic assignments.
        """
        f1 = self.ndk+self.alpha
        f2 = np.sum(self.ndk, axis=1, keepdims=True)+self.nTopics*self.alpha
        return f1/f2

    def calc_phi_topic(self):
        """Calculate phi based on the current word/topic assignments.
        """
        f1 = self.nkw+self.beta
        f2 = np.sum(self.nkw, axis=1, keepdims=True)+self.VT*self.beta
        return f1/f2

    def calc_phi_opinion(self, persp):
        """Calculate phi based on the current word/topic assignments.
        """
        f1 = self.nrs[persp]+float(self.beta_o)
        f2 = np.sum(self.nrs[persp], axis=1, keepdims=True)+self.VO*self.beta_o
        return f1/f2

    def run(self):
        logger.info('started sampling ({})'.format(str(self)))
        if not self.out_dir:
            # store all parameter samples in memory
            logger.debug('storing parameter samples in memory')
            self.theta_topic = np.zeros((self.nIter, self.DT, self.nTopics))
            self.phi_topic = np.zeros((self.nIter, self.nTopics, self.VT))

            self.phi_opinion = [np.zeros((self.nIter, self.nTopics, self.VO))
                                for p in self.corpus.perspectives]

        # variable to store nk for each iteration
        # nk is used for Contrastive Opinion Modeling
        self.nks = np.zeros((self.nIter, self.nTopics), dtype=np.int)

        for t in range(self.nIter):
            t1 = time.clock()
            logger.debug('iteration {} of {}'.format(t+1, self.nIter))

            gibbs_inner(self)

            # calculate theta and phi
            logger.debug('calculate theta and phi')
            if not self.out_dir:
                self.theta_topic[t] = self.calc_theta_topic()
                self.phi_topic[t] = self.calc_phi_topic()

                for p in range(self.nPerspectives):
                    self.phi_opinion[p][t] = self.calc_phi_opinion(p)
            else:
                if t == 0 or (t+1) % self.sampleInterval == 0:
                    fname = self.get_theta_file_name(t)
                    np.save(fname, self.calc_theta_topic())
                    logger.debug('writing to disk {}'.format(fname))

                    fname = self.get_phi_topic_file_name(t)
                    np.save(fname, self.calc_phi_topic())
                    logger.debug('writing to disk {}'.format(fname))

                    for p in range(self.nPerspectives):
                        fname = self.get_phi_opinion_file_name(p, t)
                        np.save(fname, self.calc_phi_opinion(p))
                        logger.debug('writing to disk {}'.format(fname))

                # save nk (for Contrastive Opinion Mining)
                self.nks[t] = np.copy(self.nk)
                np.save(self.get_nks_file_name(), self.nks)

            t2 = time.clock()
            logger.debug('time elapsed: {}'.format(t2-t1))

        logger.info('finished sampling ({})'.format(str(self)))

        self.estimate_parameters()

    def estimate_parameters(self, index=None, start=None, end=None):
        """Default: return single point estimate of the last iteration"""
        logger.debug('estimating parameters (index: {}, start: {}, end: {})'
                     .format(index, start, end))
        self.theta = self.get_theta(index, start, end)
        self.topics = self.get_phi_topic(index, start, end)
        self.opinions = self.get_phi_opinion(index, start, end)

    def get_theta(self, index=None, start=None, end=None):
        index = self._check_index(index)
        start, end = self._check_start_and_end(start=start, end=end)

        if hasattr(self, 'theta_topic'):
            logger.debug('retrieving theta from memory')
            if start and end:
                return np.mean(self.theta_topic[start:end], axis=0)
            return self.theta_topic[index]
        elif hasattr(self, 'out_dir'):
            logger.debug('retrieving theta from file')
            return self.load_parameters(self.THETA, index=index, start=start,
                                        end=end)
        else:
            # error: no parameter samples.
            # TODO: properly handle this case
            logger.error('no parameter samples found. Please run the Gibbs '
                         'sampler before trying to retrieve paramters.')

    def get_phi_topic(self, index=None, start=None, end=None):
        index = self._check_index(index)
        start, end = self._check_start_and_end(start=start, end=end)

        if hasattr(self, 'phi_topic'):
            logger.debug('retrieving phi topic from memory')
            if start and end:
                return np.mean(self.phi_topic[start:end], axis=0)
            return self.phi_topic[index]
        elif hasattr(self, 'out_dir'):
            logger.debug('retrieving phi topic from file')
            return self.load_parameters(self.PHI_TOPIC, index=index,
                                        start=start, end=end)
        else:
            # error: no parameter samples.
            # TODO: properly handle this case
            logger.error('no parameter samples found. Please run the Gibbs '
                         'sampler before trying to retrieve paramters.')

    def get_phi_opinion(self, index=None, start=None, end=None):
        index = self._check_index(index)
        start, end = self._check_start_and_end(start=start, end=end)

        phi_opinion = [np.zeros((self.nPerspectives, self.nTopics, self.VO))
                       for p in self.corpus.perspectives]

        if hasattr(self, 'phi_opinion'):
            logger.debug('retrieving phi opinion from memory')
            if start and end:
                for p in range(self.nPerspectives):
                    phi_opinion[p] = np.mean(self.phi_opinion[p][start:end],
                                             axis=0)
                return phi_opinion
            for p in range(self.nPerspectives):
                phi_opinion[p] = self.phi_opinion[p][index]
            return phi_opinion
        elif hasattr(self, 'out_dir'):
            logger.debug('retrieving phi opinion from from file')
            for p in range(self.nPerspectives):
                phi_opinion[p] = self.load_parameters(self.PHI_OPINION.
                                                      format(p),
                                                      index=index, start=start,
                                                      end=end)
            return phi_opinion
        else:
            # error: no parameter samples.
            # TODO: properly handle this case
            logger.error('no parameter samples found. Please run the Gibbs '
                         'sampler before trying to retrieve paramters.')

    def _check_index(self, index=None):
        if index is None:
            index = self.nIter-1
        if index == self.nIter:
            logger.warn('requested index {}; setting it to {}'.format(index,
                        self.nIter-1))
            index = self.nIter-1
        return index

    def _check_start_and_end(self, start=None, end=None):
        if start is None and end is None:
            return start, end
        if start < 0 or start > self.nIter or start > end:
            start = 0

        if end < 0 or end > self.nIter:
            end = self.nIter

        return start, end

    def perplexity(self, index=None, phi_topic=None, phi_opinion=None):
        """Calculate topic word and opinion word perplexity of the test set.

        Parameters:
            index : int (optional)
                index of the parameter sample of phi topic to use for a single
                point parameter estimate of phi topic (default nIter-1)
            phi_topic : ndarray (optional)
                phi_topic is used to initialize the Gibbs sampler to determine
                theta for the test set
            phi_opinion : ndarray (optional)
                phi_opinion provides estimates of p(o_i|z_i=k)

        Returns:
            topic word perplexity : float
            opinion word perplexity : float

        Topic word perplexity is calculated using importance sampling, as in
        [Griffiths and Steyvers, 2004]. However, according to
        [Wallach et al. 2009], this does not always result in an accurate
        estimate of the perplexity.

        Opinion word perplexity is calculated as described in [Fang et al.,
        2012] section 5.1.1.
        """
        # TODO: implement more accurate estimate of perplexity
        logger.info('calculating perplexity ({})'.format(str(self)))

        # load parameters
        if phi_topic is None:
            phi_topic = self.get_phi_topic(index)

        if phi_opinion is None:
            phi_opinion = self.get_phi_opinion(index)

        # run Gibbs sampler to find estimates for theta of the test set
        s = GibbsSampler(self.corpus, nTopics=self.nTopics, nIter=self.nIter,
                         alpha=self.alpha, beta=self.beta, beta_o=self.beta_o,
                         initialize=False)
        s._initialize(phi_topic=phi_topic)
        s.run()

        # topic word perplexity
        total_topic_words_in_test_documents = 0
        log_p_w = 0.0

        # opinion word perplexity
        total_opinion_words_in_test_documents = 0
        log_p_od = 0.0

        for d, persp, d_p, doc in self.corpus.testSet():
            for w_id, freq in doc['topic']:
                total_topic_words_in_test_documents += freq
                log_p_w += freq * np.log(np.sum(s.theta[d]*phi_topic[:, w_id]))

            for w_id, freq in doc['opinion']:
                total_opinion_words_in_test_documents += freq
                f1 = np.log(np.sum(s.theta[d]*phi_opinion[persp][:, w_id]))
                log_p_od += freq * f1

        tw_perp = np.exp(-log_p_w/total_topic_words_in_test_documents)
        ow_perp = np.exp(-log_p_od/total_opinion_words_in_test_documents)

        return tw_perp, ow_perp

    def opinion_word_perplexity_per_document(self, index=None, phi_topic=None,
                                             phi_opinion=None):
        """Calculate opinion word perplexity per document

        Opinion word perplexity is calculated per document in the
        test set.
        """
        logger.info('calculating opinion wor perplexity per document '
                    '({})'.format(str(self)))

        # load parameters
        if phi_topic is None:
            phi_topic = self.get_phi_topic(index)

        if phi_opinion is None:
            phi_opinion = self.get_phi_opinion(index)

        # run Gibbs sampler to find estimates for theta and phi_opinions
        # of the test set
        out_dir = os.path.join(self.out_dir, 'per_document_ow_perplexity')
        s = GibbsSampler(self.corpus, nTopics=self.nTopics, nIter=self.nIter,
                         alpha=self.alpha, beta=self.beta, beta_o=self.beta_o,
                         out_dir=out_dir, initialize=False)
        s._initialize(phi_topic=phi_topic)
        s.run()

        file_dict = s.corpus.get_files_in_train_and_test_sets()

        results = {}
        for d, persp, d_p, doc in self.corpus.testSet():
            opinion_words_in_document = 0
            log_p_od = 0.0

            for w_id, freq in doc['opinion']:
                opinion_words_in_document += freq
                f1 = np.log(np.sum(s.theta[d] * phi_opinion[persp][:, w_id]))
                log_p_od += freq * f1

            ow_perp = np.exp(-log_p_od / opinion_words_in_document)
            p_name = self.corpus.perspectives[persp].name

            if not results.get(p_name):
                results[p_name] = []
            results[p_name].append(ow_perp)

        df_index = [os.path.basename(f) for f in file_dict['0']['test']]

        # get theta
        documents = []
        for persp in self.corpus.perspectives:
            print persp
            documents += df_index
        print documents
        theta = s.theta_to_df(s.theta, documents)

        return pd.DataFrame(results, index=df_index), theta

    def load_parameters(self, name, index=None, start=None, end=None):
        index = self._check_index(index)
        start, end = self._check_start_and_end(start, end)

        if start and end:
            logger.debug('loading parameter files {} - {}'.
                         format(self.get_parameter_file_name(name, start),
                                self.get_parameter_file_name(name, end)))
            data = None
            for i in range(start, end):
                fName = self.get_parameter_file_name(name, i)
                if os.path.isfile(fName):
                    ar = np.load(fName)
                    if data is None:
                        data = np.array([ar])
                    else:
                        data = np.append(data, np.array([ar]), axis=0)
            return np.mean(data, axis=0)

        if not index is None:
            logger.debug('loading parameter file {}'.
                         format(self.get_parameter_file_name(name, index)))
            return np.load(self.get_parameter_file_name(name, index))

    def get_phi_topic_file_name(self, number):
        return self.get_parameter_file_name(self.PHI_TOPIC, number)

    def get_phi_opinion_file_name(self, persp, number):
        return self.get_parameter_file_name(self.PHI_OPINION.format(persp),
                                            number)

    def get_theta_file_name(self, number):
        return self.get_parameter_file_name(self.THETA, number)

    def get_nks_file_name(self):
        return '{}/{}.npy'.format(self.get_parameter_dir_name(), self.NKS)

    def get_parameter_file_name(self, name, number):
        return '{}/{}_{:04d}.npy'.format(self.get_parameter_dir_name(),
                                         name, number)

    def get_parameter_dir_name(self):
        return self.PARAMETER_DIR.format(self.out_dir)

    def print_topics_and_opinions(self, topics, opinions, top=10,
                                  threshold=0.0):
        """Print topics and associated opinions.

        Parameters:
            topics : pandas DataFrame
            opinions : list of pandas DataFrames

        The <top> top words and weights are printed.
        """
        for i in range(self.nTopics):
            print u'Topic #{}: {}'. \
                  format(i, self.print_topic(topics.loc[:, i].copy(), top=top,
                                             threshold=threshold))
            print
            for p in range(self.nPerspectives):
                print u'Opinion {}: {}'. \
                      format(self.corpus.perspectives[p].name,
                             self.print_topic(opinions[p].loc[:, i].copy(),
                                              top=top, threshold=threshold))
            print '-----'
            print

    def print_topic(self, series, top=10, threshold=0.0):
        """Prints the top 10 words in the topic/opinion on a single line."""
        s = series.copy()
        s.sort(ascending=False)
        t = [u'{} ({:.4f})'.format(word, p)
             for word, p in s[0:top].iteritems() if p > threshold]
        return u' - '.join(t)

    def topics_to_df(self, phi, words):
        """Returns a pandas DataFrame containing words x probabilities.

        The DataFrame is used as input for print_topics_and_opinions() and
        print_topic().

        Parameters:
            phi : numpy array
                Array containing the word probabilities for the topics.
            words : list of strings
                One of CPTCorpus.topic_words() and CPTCorpus.opinion_words().

        Returns:
            pandas DataFrame
                containing words x probabilities
        """
        logger.debug('creating dataframe with topics')
        df = pd.DataFrame(phi, columns=words)
        df = df.transpose()
        return df

    def theta_to_df(self, theta, documents):
        """Returns a pandas DataFrame containing doc ids x topic probabilities.

        Parameters:
        theta : numpy array
            Array containing the word probabilities for the topics.
        documents : list of strings
            list of document ids.

        Returns:
        pandas DataFrame
            containing document ids x topic probabilities.
        """
        logger.debug('creating dataframe with theta')
        df = pd.DataFrame(theta, index=documents)
        return df

    def __str__(self):
        return 'CPT GibbsSampler: {} perspectives, {} topics,  {} ' \
               'iterations, alpha: {}, beta: {}, beta_o: {}'. \
               format(self.nPerspectives, self.nTopics, self.nIter,
                      self.alpha, self.beta, self.beta_o)

if __name__ == '__main__':
    #logger.setLevel(logging.DEBUG)
    #logger.setLevel(logging.INFO)

    #files = glob.glob('/home/jvdzwaan/data/tmp/dilipad/gov_opp/*')
    files = glob.glob('/home/jvdzwaan/data/tmp/test/*')
    #files = glob.glob('/home/jvdzwaan/data/dilipad/perspectives/*')
    #out_dir = '/home/jvdzwaan/data/tmp/generated/test/'
    out_dir = '/home/jvdzwaan/data/tmp/dilipad/test_perplexity/'

    corpus = CPTCorpus.CPTCorpus(files, testSplit=20)
    #corpus.filter_dictionaries(minFreq=5, removeTopTF=100, removeTopDF=100)
    #corpus.save_dictionaries(directory=out_dir)
    #corpus.save('{}corpus.json'.format(out_dir))
    #corpus = CPTCorpus.CPTCorpus.load('{}corpus.json'.format(out_dir),
    #                                  topicDict='{}/topicDict.dict'.format(out_dir),
    #                                  opinionDict='{}/opinionDict.dict'.format(out_dir))
    #sampler = GibbsSampler(corpus, nTopics=30, nIter=100, out_dir=out_dir)
    sampler = GibbsSampler(corpus, nTopics=3, nIter=10, out_dir=out_dir)
    #sampler = GibbsSampler(corpus, nTopics=100, nIter=2)
    sampler._initialize()
    sampler.run()

    print sampler

    #topics_df = sampler.topics_to_df(sampler.get_phi_topic(),
    #                                 sampler.corpus.topic_words())
    #opinion_df = [sampler.topics_to_df(sampler.get_phi_opinion()[p],
    #                                   sampler.corpus.opinion_words())
    #              for p in range(sampler.nPerspectives)]
    #sampler.print_topics_and_opinions(topics_df, opinion_df, threshold=0.09)
    ps = []
    for i in range(0, 11):
        ps.append(sampler.perplexity(index=i))
    print ps

max_stars_repo_name: egpbos/cptm

max_issues_repo_name: egpbos/cptm

max_forks_repo_name: egpbos/cptm


text:
# use time() instead on unix
import sys
if sys.platform=='win32':
    from time import clock
else:
    from time import time as clock

from sympycore import profile_expr

def time1(n=500):
    import sympycore as sympy
    w = sympy.mpq((3,4))
    x = sympy.polynomials.poly([0, 1, 1])
    b = (x-1)*(x-2)*(x-w)
    a = (x-1)*(x-2)*(x-3)*(x-4)*(x-5)# + (x-1)*(x-2)*x**10

    t1 = clock()
    while n:
        divmod(a, b)
        n -= 1
    t2 = clock()
    return 100 / (t2-t1)

def time2(n=500):
    import sympycore as sympy
    w = sympy.mpq((3,4))
    x = sympy.polynomials.PolynomialRing[1]([0, 1, 1])
    b = (x-1)*(x-2)*(x-w)
    a = (x-1)*(x-2)*(x-3)*(x-4)*(x-5)# + (x-1)*(x-2)*x**10
    t1 = clock()
    while n:
        divmod(a, b)
        n -= 1
    t2 = clock()
    return 100 / (t2-t1)

def timing():
    t1 = time1()
    t2 = time2()
    return t1, t2, t1/t2

if 1:
    print timing()
    print timing()
    print timing()

profile_expr('time2(100)')

max_stars_repo_name: pearu/sympycore

max_issues_repo_name: pearu/sympycore

max_forks_repo_name: pearu/sympycore

reason: import sympy,from sympy