text:
import numpy as np
import torch
import argparse
import os
import math
import gym
import sys
import random
import time
import json
import dmc2gym
import copy

import utils
from logger import Logger
from video import VideoRecorder

from curl_sac import RadSacAgent
from multistep_replay import ReplayBufferPixelMultistep
import data_augs as rad

def parse_args():
    parser = argparse.ArgumentParser()
    # environment
    parser.add_argument('--domain_name', default='cartpole')
    parser.add_argument('--task_name', default='swingup')
    parser.add_argument('--pre_transform_image_size', default=100, type=int)
    parser.add_argument('--case', type=str, default='1') #cooresponding to overleaf cases

    parser.add_argument('--image_size', default=84, type=int)
    parser.add_argument('--action_repeat', default=1, type=int)
    parser.add_argument('--frame_stack', default=3, type=int)

    parser.add_argument('--resource_files', default=None, type=str)
    parser.add_argument('--eval_resource_files', type=str)
    parser.add_argument('--img_source', default=None, type=str, choices=['color', 'noise', 'black', 'images', 'video', 'none'])
    parser.add_argument('--total_frames', default=1000, type=int)

    # replay buffer
    parser.add_argument('--replay_type', default='multistep', type=str, choices=['multistep', 'singlestep'])
    parser.add_argument('--replay_buffer_capacity', default=100000, type=int)
    parser.add_argument('--normalize_obs', default=False, action='store_true')
    # train
    parser.add_argument('--init_steps', default=1000, type=int)
    parser.add_argument('--num_train_steps', default=1000000, type=int)
    parser.add_argument('--batch_size', default=32, type=int)
    parser.add_argument('--hidden_dim', default=1024, type=int)
    # eval
    parser.add_argument('--eval_freq', default=1000, type=int)
    parser.add_argument('--num_eval_episodes', default=10, type=int)
    # critic
    parser.add_argument('--critic_lr', default=1e-3, type=float)
    parser.add_argument('--critic_beta', default=0.9, type=float)
    parser.add_argument('--critic_tau', default=0.01, type=float) # try 0.05 or 0.1
    parser.add_argument('--critic_target_update_freq', default=2, type=int) # try to change it to 1 and retain 0.01 above
    # actor
    parser.add_argument('--actor_lr', default=1e-3, type=float)
    parser.add_argument('--actor_beta', default=0.9, type=float)
    parser.add_argument('--actor_log_std_min', default=-10, type=float)
    parser.add_argument('--actor_log_std_max', default=2, type=float)
    parser.add_argument('--actor_update_freq', default=2, type=int)
    # encoder
    parser.add_argument('--encoder_type', default='pixel', type=str)
    parser.add_argument('--encoder_feature_dim', default=50, type=int)
    parser.add_argument('--encoder_lr', default=1e-3, type=float)
    parser.add_argument('--encoder_tau', default=0.05, type=float)
    parser.add_argument('--num_layers', default=4, type=int)
    parser.add_argument('--num_filters', default=32, type=int)
    parser.add_argument('--latent_dim', default=128, type=int)
    parser.add_argument('--metric_loss', default=False, action='store_true')
    # sac
    parser.add_argument('--discount', default=0.99, type=float)
    parser.add_argument('--init_temperature', default=0.1, type=float)
    parser.add_argument('--alpha_lr', default=1e-4, type=float)
    parser.add_argument('--alpha_beta', default=0.5, type=float)
    # multistep
    parser.add_argument('--decoder_lr', default=1e-3, type=float)
    parser.add_argument('--decoder_weight_lambda', default=0.0, type=float)
    parser.add_argument('--horizon', default=3, type=int)
    # misc
    parser.add_argument('--seed', default=1, type=int)
    parser.add_argument('--work_dir', default='.', type=str)
    parser.add_argument('--save_tb', default=False, action='store_true')
    parser.add_argument('--save_buffer', default=False, action='store_true')
    parser.add_argument('--save_video', default=False, action='store_true')
    parser.add_argument('--save_model', default=False, action='store_true')
    parser.add_argument('--detach_encoder', default=False, action='store_true')
    parser.add_argument('--off_center', default=False, action='store_true')
    # data augs
    parser.add_argument('--data_augs', default='no_aug', type=str)

    parser.add_argument('--log_interval', default=100, type=int)
    args = parser.parse_args()
    return args

def evaluate(env, agent, video, num_episodes, L, step, args):
    all_ep_rewards = []

    def run_eval_loop(sample_stochastically=True):
        start_time = time.time()
        prefix = 'stochastic_' if sample_stochastically else ''
        for i in range(num_episodes):
            obs = env.reset()
            video.init(enabled=(i == 0))
            done = False
            episode_reward = 0
            while not done:
                # center crop image
                if args.encoder_type == 'pixel' and 'crop' in args.data_augs:
                    obs = utils.center_crop_image(obs,args.pre_transform_image_size)
                if args.encoder_type == 'pixel' and 'translate' in args.data_augs:
                    # first crop the center with pre_image_size
                    obs = utils.center_crop_image(obs, args.pre_transform_image_size)
                    # then translate cropped to center
                    obs = utils.center_translate(obs, args.image_size)

                with utils.eval_mode(agent):
                    if sample_stochastically:
                        action = agent.sample_action(obs / 255.)
                    else:
                        action = agent.select_action(obs / 255.)
                obs, reward, done, _ = env.step(action)
                video.record(env)
                episode_reward += reward

            video.save('%d.mp4' % step)
            L.log('eval/' + prefix + 'episode_reward', episode_reward, step)
            all_ep_rewards.append(episode_reward)

        L.log('eval/' + prefix + 'eval_time', time.time()-start_time , step)
        mean_ep_reward = np.mean(all_ep_rewards)
        best_ep_reward = np.max(all_ep_rewards)
        std_ep_reward = np.std(all_ep_rewards)
        L.log('eval/' + prefix + 'mean_episode_reward', mean_ep_reward, step)
        L.log('eval/' + prefix + 'best_episode_reward', best_ep_reward, step)

        filename = args.work_dir + '/' + args.domain_name + '--'+args.task_name + '-' + args.data_augs + '--s' + str(args.seed) + '--eval_scores.npy'
        key = args.domain_name + '-' + args.task_name + '-' + args.data_augs
        try:
            log_data = np.load(filename,allow_pickle=True)
            log_data = log_data.item()
        except:
            log_data = {}

        if key not in log_data:
            log_data[key] = {}

        log_data[key][step] = {}
        log_data[key][step]['step'] = step
        log_data[key][step]['mean_ep_reward'] = mean_ep_reward
        log_data[key][step]['max_ep_reward'] = best_ep_reward
        log_data[key][step]['std_ep_reward'] = std_ep_reward
        log_data[key][step]['env_step'] = step * args.action_repeat

        np.save(filename,log_data)

    run_eval_loop(sample_stochastically=False)
    L.dump(step)

def main():
    args = parse_args()
    if args.seed == -1:
        args.__dict__["seed"] = np.random.randint(1,1000000)
    utils.set_seed_everywhere(args.seed)

    if args.case == '0':
        from agents.agent_sac import PixelSacAgent
        assert args.horizon == 1, 'Horizon should be 1'
    elif args.case == '1':
        from agents.agent_sac_base import PixelSacAgent
        assert args.horizon == 1, 'Horizon should be 1'
    elif args.case == '2':
        from agents.agent_sac_noreward import PixelSacAgent
        assert args.horizon == 1, 'Horizon should be 1'
    elif args.case == '3':
        from agents.agent_sac_notransition import PixelSacAgent
        assert args.horizon == 1, 'Horizon should be 1'
    elif args.case == '4':
        from agents.agent_sac_norewtotransition import PixelSacAgent
        assert args.horizon == 1, 'Horizon should be 1'
    elif args.case == '5':
        from agents.agent_sac_value import PixelSacAgent
    elif args.case == '6':
        from agents.agent_sac_reconstruction import PixelSacAgent
    elif args.case == '7':
        from agents.agent_sac_contrastive import PixelSacAgent
    elif args.case == '8':
        from agents.agent_sac_SPR import PixelSacAgent
        args.horizon = 5

    pre_transform_image_size = args.pre_transform_image_size if 'crop' in args.data_augs else args.image_size
    pre_image_size = args.pre_transform_image_size # record the pre transform image size for translation

    env = dmc2gym.make(
        domain_name=args.domain_name,
        task_name=args.task_name,
        resource_files=args.resource_files,
        img_source=args.img_source,
        total_frames=args.total_frames,
        seed=args.seed,
        visualize_reward=False,
        from_pixels=(args.encoder_type == 'pixel'),
        height=pre_transform_image_size,
        width=pre_transform_image_size,
        frame_skip=args.action_repeat,
        off_center=args.off_center
    )
    env.seed(args.seed)

    eval_env = dmc2gym.make(
        domain_name=args.domain_name,
        task_name=args.task_name,
        resource_files=args.resource_files,
        img_source=args.img_source,
        total_frames=args.total_frames,
        seed=args.seed,
        visualize_reward=False,
        from_pixels=(args.encoder_type == 'pixel'),
        height=pre_transform_image_size,
        width=pre_transform_image_size,
        frame_skip=args.action_repeat,
        off_center=args.off_center
    )
    eval_env.seed(args.seed)

    # stack several consecutive frames together
    if args.encoder_type == 'pixel':
        env = utils.FrameStack(env, k=args.frame_stack)
        eval_env = utils.FrameStack(eval_env, k=args.frame_stack)

    # make directory
    ts = time.gmtime()
    ts = time.strftime("%m-%d", ts)
    env_name = args.domain_name + '-' + args.task_name
    exp_name = 'case-' + args.case + '-h' + str(args.horizon) + '-' + env_name + '-' + ts + '-im' + str(args.image_size) +'-b'  \
    + str(args.batch_size) + '-s' + str(args.seed)  + '-' + args.encoder_type
    args.work_dir = args.work_dir + '/'  + exp_name

    utils.make_dir(args.work_dir)
    video_dir = utils.make_dir(os.path.join(args.work_dir, 'video'))
    model_dir = utils.make_dir(os.path.join(args.work_dir, 'model'))
    buffer_dir = utils.make_dir(os.path.join(args.work_dir, 'buffer'))

    video = VideoRecorder(video_dir if args.save_video else None)

    utils.check_path(args.work_dir)

    with open(os.path.join(args.work_dir, 'args.json'), 'w') as f:
        json.dump(vars(args), f, sort_keys=True, indent=4)

    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

    action_shape = env.action_space.shape

    if args.encoder_type == 'pixel':
        obs_shape = (3*args.frame_stack, args.image_size, args.image_size)
        pre_aug_obs_shape = (3*args.frame_stack,pre_transform_image_size,pre_transform_image_size)

    replay_buffer = ReplayBufferPixelMultistep(
        obs_shape=pre_aug_obs_shape,
        action_shape=action_shape,
        capacity=args.replay_buffer_capacity,
        batch_size=args.batch_size,
        horizon=args.horizon+1,
        device=device,
        normalize_obs=args.normalize_obs
    )

    agent = PixelSacAgent(
        obs_shape=pre_aug_obs_shape,
        action_shape=action_shape,
        horizon=args.horizon,
        device=device,
        hidden_dim=args.hidden_dim,
        discount=args.discount,
        init_temperature=args.init_temperature,
        alpha_lr=args.alpha_lr,
        alpha_beta=args.alpha_beta,
        actor_lr=args.actor_lr,
        actor_beta=args.actor_beta,
        actor_log_std_min=args.actor_log_std_min,
        actor_log_std_max=args.actor_log_std_max,
        actor_update_freq=args.actor_update_freq,
        critic_lr=args.critic_lr,
        critic_beta=args.critic_beta,
        critic_tau=args.critic_tau,
        critic_target_update_freq=args.critic_target_update_freq,
        encoder_type=args.encoder_type,
        encoder_feature_dim=args.encoder_feature_dim,
        encoder_lr=args.encoder_lr,
        encoder_tau=args.encoder_tau,
        decoder_lr=args.decoder_lr,
        decoder_weight_lambda=args.decoder_weight_lambda,
        num_layers=args.num_layers,
        num_filters=args.num_filters,
        log_interval=args.log_interval,
        detach_encoder=args.detach_encoder,
        latent_dim=args.latent_dim,
        data_augs=args.data_augs,
        use_metric_loss=args.metric_loss
    )

    L = Logger(args.work_dir, use_tb=args.save_tb)

    episode, episode_reward, done = 0, 0, True
    start_time = time.time()

    for step in range(args.num_train_steps):
        # evaluate agent periodically

        if step % args.eval_freq == 0:
            L.log('eval/episode', episode, step)
            evaluate(eval_env, agent, video, args.num_eval_episodes, L, step,args)
            if args.save_model: # and args.agent == 'rad_sac':
                agent.save(model_dir, step)
            if args.save_buffer:
                replay_buffer.save(buffer_dir)

        if done:
            if step > 0:
                if step % args.log_interval == 0:
                    L.log('train/duration', time.time() - start_time, step)
                    L.dump(step)
                start_time = time.time()
            if step % args.log_interval == 0:
                L.log('train/episode_reward', episode_reward, step)

            obs = env.reset()
            done = False
            episode_reward = 0
            episode_step = 0
            episode += 1
            if step % args.log_interval == 0:
                L.log('train/episode', episode, step)

        # sample action for data collection
        if step < args.init_steps:
            action = env.action_space.sample()
        else:
            with utils.eval_mode(agent):
                action = agent.sample_action(obs / 255.)

        # run training update
        if step >= args.init_steps:
            num_updates = 1
            for _ in range(num_updates):
                agent.update(replay_buffer, L, step)

        next_obs, reward, done, _ = env.step(action)

        # allow infinit bootstrap
        done_bool = 0 if episode_step + 1 == env._max_episode_steps else float(
            done
        )
        episode_reward += reward
        replay_buffer.add(obs, action, reward, next_obs, done_bool)

        obs = next_obs
        episode_step += 1

if __name__ == '__main__':
    torch.multiprocessing.set_start_method('spawn')

    main()


text:
# Various 2D filters that can be applied to data.
from .filters1d import blockavg1d
import numpy as np
from astropy.modeling.models import Gaussian2D
from scipy.signal import gaussian, fftconvolve
############################################################################
# Master function that assigns operations based on chosen smoothing mode.
def smooth_master(data,smoothing_method,convolution_method,kernel_size_freq_chans,kernel_size_time_samples,freqs,times):
    if (smoothing_method=='Gaussian2D'):
        smooth_data = smooth_2DGaussian_zeroCCF(data,convolution_method,kernel_size_freq_chans,kernel_size_time_samples)
        return smooth_data, freqs, times
    if (smoothing_method=='Blockavg2D'):
        smooth_data, blkavg_freqs, blkavg_times = blockavg_ds(data,kernel_size_freq_chans,kernel_size_time_samples,freqs,times)
        return smooth_data, blkavg_freqs, blkavg_times
    if (smoothing_method in ['blackman','hanning','hamming']):
        smooth_data = pass_2dfilter(data,convolution_method,smoothing_method,kernel_size_freq_chans,kernel_size_time_samples)
        return smooth_data, freqs, times
############################################################################
# Function to recast a 2D array to a larger size by padding with zeros.
# Specify the start indices along the two axes. For i>start_index, non-zero values are permitted.
'''
Inputs:
array2D = 2D array of non-zeros values that needs to be recast into a larger shape by padding with zeroes
output_shape = Final shape of 2D array to be output
start_index_axis0 = Start index along axis 0 for permitting non-zero values
start_index_axis1 = Start index along axis 1 for permitting non-zero values
'''
def recast_2Darray(array2D, output_shape, start_index_axis0, start_index_axis1):
    output_array = np.zeros(output_shape)
    nonzero_length_axis0, nonzero_length_axis1 = array2D.shape
    stop_index_axis0 = start_index_axis0 + nonzero_length_axis0
    stop_index_axis1 = start_index_axis1 + nonzero_length_axis1
    output_array[start_index_axis0:stop_index_axis0,start_index_axis1:stop_index_axis1] = array2D
    return output_array
############################################################################
# Calculate 2D convolution using FFT.
'''
Inputs:
arr1, arr2 = Two 2D arrays of the same shape to be convolved.
'''
def convolve_fft2D(arr1,arr2):
    conv_product = np.fft.ifft2(np.fft.fft2(arr1)*np.fft.fft2(arr2)).real
    conv_product = np.fft.fftshift(conv_product)
    return conv_product
############################################################################
# Convole dynamic spectrum with a 2D Gaussian kernel assuming zero cross-correlation coefficient (CCF).
'''
Inputs:
data = dynamic spectrum, 2D array = [Frequency, Time]
convolution_method = 'fftconvolve' or 'fft2'
chan_fwhm = FWHM of 2D Gaussian along frequency axis measured in units of number of channels
t_sample_fwhm = Temporal FWHM of 2D Gaussian measured in units of number of time samples
'''
def smooth_2DGaussian_zeroCCF(data,convolution_method,chan_fwhm,t_sample_fwhm):
    chan_sigma = np.round(chan_fwhm/np.sqrt(8*np.log(2))).astype(int)
    t_sample_sigma = np.round(t_sample_fwhm/np.sqrt(8*np.log(2))).astype(int)
    print('Smoothing data with a 2D Gaussian kernel assuming zero CCF...')
    # 1D Gaussian along frequency axis
    freq_extent = 20*chan_sigma
    freq_Gaussian = gaussian(freq_extent,chan_sigma)
    # 1D Gaussian along time axis
    time_extent = 20*t_sample_sigma
    time_Gaussian = gaussian(time_extent,t_sample_sigma)
    # Construct 2D Gaussian from two 1D Gaussians assuming zero CCF.
    Gaussian_2D = np.outer(freq_Gaussian,time_Gaussian)
    Gaussian_2D = Gaussian_2D/np.sum(Gaussian_2D) # Normalize the 2D Gaussian to unit area.
    # Convolve 2D Gaussian with data.
    if (convolution_method=='fftconvolve'):
        output = fftconvolve(data,Gaussian_2D,mode='same')
    elif (convolution_method=='fft2'):
        # Recast 2D Gaussian to same shape as the data array.
        output_shape = data.shape
        start_index_axis0 = (data.shape[0]-freq_extent)//2
        start_index_axis1 = (data.shape[1]-time_extent)//2
        Gaussian_2D = recast_2Darray(Gaussian_2D, output_shape, start_index_axis0, start_index_axis1)
        output = convolve_fft2D(data,Gaussian_2D)
    return output
############################################################################
# Generate a general 2D Gaussian with user-supplied std. dev. along x-axis and y-axis.
# Angle between major axis of 2D gaussian and x-axis can also be specified.
# Also, a full covariance matrix can be provided.
'''
Inputs:
shape = Dimensions of 2D array to output = (No. of y-axis values, No. of x-axis samples)
x_fwhm = FWHM of 2D Gaussian along x-axis
y_fwhm = FWHM of 2D Gaussian along y-axis
theta = Rotation angle in radians. The rotation angle increases counterclockwise. Must be None if a covariance matrix (cov_matrix) is provided. If no cov_matrix is given, None means the default value (0).
cov_matrix = A 2x2 covariance matrix. If specified, overrides the x_stddev, y_stddev, and theta defaults.
'''
def gen_general_2DGaussian(shape,x_fwhm,y_fwhm,theta=None,cov_matrix=None):
    # Convert FWHM to sigma.
    x_stddev = int(np.round(x_fwhm/np.sqrt(8*np.log(2))))
    y_stddev = int(np.round(y_fwhm/np.sqrt(8*np.log(2))))
    # Initialize grid of x and y values.
    y_values = np.arange(shape[0])
    x_values = np.arange(shape[1])
    y_values, x_values = np.meshgrid(y_values,x_values)
    # Coordinates of the center.
    y_center = shape[0]//2
    x_center = shape[1]//2
    # Initialize Gaussian2D object.
    gauss2d_object = Gaussian2D(amplitude=1.0,x_mean=x_center, y_mean=y_center,x_stddev = x_stddev, y_stddev=y_stddev,theta=theta,cov_matrix=cov_matrix)
    gauss2d_array = gauss2d_object(freq_chans,time_samples)
    return gauss2d_array
############################################################################
# Perform a general 2D Gaussian smoothing of a dynamic spectrum.
'''
Inputs:
data = dynamic spectrum, 2D array = [Frequency, Time]
chan_fwhm = FWHM of 2D Gaussian along frequency axis measured in units of number of channels
t_sample_fwhm = Temporal FWHM of 2D Gaussian measured in units of number of time samples
theta = Rotation angle in radians. The rotation angle increases counterclockwise. Must be None if a covariance matrix (cov_matrix) is provided. If no cov_matrix is given, None means the default value (0).
cov_matrix = A 2x2 covariance matrix. If specified, overrides the x_stddev, y_stddev, and theta defaults.
'''
def smooth_2DGaussian(data,chan_fwhm,t_sample_fwhm,theta=None,cov_matrix=None):
    chan_sigma = np.round(chan_fwhm/np.sqrt(8*np.log(2))).astype(int)
    t_sample_sigma = np.round(t_sample_fwhm/np.sqrt(8*np.log(2))).astype(int)
    shape = data.shape
    # Generate 2D Gaussian kernel to smooth data.
    print('Generating 2D Gaussian kernel.')
    gauss2d_kernel = gen_general_2DGaussian(shape,t_sample_sigma,chan_sigma,theta=theta,cov_matrix=cov_matrix)
    # Normalize 2D Gaussian kernel.
    gauss2d_kernel = gauss2d_kernel/np.sum(gauss2d_kernel)
    # Smooth the data by convolving it with a 2D Gaussian.
    print('Smoothing data with a 2D Gaussian kernel...')
    smooth_data = convolve_fft2D(data,gauss2d_kernel)
    print('2D Gaussian smoothing complete.')
    return smooth_data
############################################################################
# Block average along the second axis of a 2D array.
'''
Inputs:
data = 2D array
R = Window length for performing block averaging.
'''
def blockavg_axis2(data,R):
    R = int(R) # Sanity check
    # Shape of the original data
    shape = data.shape
    ax1_length = shape[0]
    ax2_length = shape[1]
    # If ax2_length is not divisible by the length  R, pad the array with NaNs to a size divisible by R.
    pad_size = int(np.ceil(float(ax2_length)/R)*R - float(ax2_length))
    pad_array = np.zeros((ax1_length,pad_size))*np.NaN
    # Pad the data with pad_array.
    if np.ma.is_masked(data):
        padded_data = np.ma.concatenate((data,pad_array),axis=1)
    else:
        padded_data = np.concatenate((data,pad_array),axis=1)
    pad_ax2_length = padded_data.shape[1]
    blkavg_ax2_data = np.nanmean(padded_data.reshape((ax1_length,pad_ax2_length//R,R)),axis=2)
    return blkavg_ax2_data
############################################################################
# Block average independently along the frequency and time axes of a dynamic spectrum
'''
Inputs:
data = dynamic spectrum, 2D array = [Frequency, Time]
blkavg_factor_time = Window length (no. of time samples) for block averaging along time
blkavg_factor_freq = Window length (no. of channels) for block averaging along frequency
times = 1D array of observation times
freqs_GHz = 1D array of radio frequencies (usually GHz)
'''
def blockavg_ds(data, blkavg_factor_freq, blkavg_factor_time, freqs_GHz, times):
    # Ensure block averaging factors are atleast one.
    blkavg_factor_time = np.max([blkavg_factor_time,1])
    blkavg_factor_freq = np.max([blkavg_factor_freq,1])
    # Ensure block averaging factors are integers. If not, round them to nearest integer.
    blkavg_factor_time = int(np.round(blkavg_factor_time))
    blkavg_factor_freq = int(np.round(blkavg_factor_freq))
    # Block average data along time.
    print('Smoothing data by block averaging...')
    print('Block averaging along time by factor %d'% (blkavg_factor_time))
    if (blkavg_factor_time>=2):
        data = blockavg_axis2(data,blkavg_factor_time)
        blkavg_times = blockavg1d(times,blkavg_factor_time)
    else:
        blkavg_times = times
    print('Block averaging along time complete.')
    # Block average data along frequency.
    print('Block averaging along frequency by factor %d'% (blkavg_factor_freq))
    if (blkavg_factor_freq>=2):
        data = blockavg_axis2(data.T,blkavg_factor_freq).T
        blkavg_freqs = blockavg1d(freqs_GHz,blkavg_factor_freq)
    else:
        blkavg_freqs = freqs_GHz
    print('Block averaging along frequency complete.')
    return data,blkavg_freqs,blkavg_times
############################################################################
# Apply a 2D filter which is product of two 1D NumPy filters (np.hanning/np.hamming/np.blackman)
'''
Inputs:
data_2d = 2D data array
convolution_method = 'fftconvolve' or 'fft2'
smoothing_method = 'blackman','hanning','hamming'
kernel_size_freq_chans = No. of channels constituting a spectral kernel size (Gaussian std. dev. / window length)
kernel_size_time_samples = No. of channels constituting a spectral kernel size (Gaussian std. dev. / window length)
'''
def pass_2dfilter(data_2d,convolution_method,smoothing_method,kernel_size_freq_chans,kernel_size_time_samples):
    print('Smoothing data with a 2D %s filter'% (smoothing_method))
    if (smoothing_method=='blackman'):
        filter = np.blackman
    if (smoothing_method=='hanning'):
        filter = np.hanning
    if (smoothing_method=='hamming'):
        filter = np.hamming
    filter_time = filter(kernel_size_time_samples)
    filter_freq = filter(kernel_size_freq_chans)
    filter_2d = np.outer(filter_freq,filter_time)
    filter_2d = filter_2d/np.sum(filter_2d)
    # Convolve 2D filter with data.
    if (convolution_method=='fftconvolve'):
        output = fftconvolve(data_2d,filter_2d,mode='same')
    elif (convolution_method=='fft2'):
        # Recast filter_2d to the same shape as the data by padding with zeroes.
        output_shape = data_2d.shape
        start_index_axis0 = (data_2d.shape[0]-kernel_size_freq_chans)//2
        start_index_axis1 = (data_2d.shape[1]-kernel_size_time_samples)//2
        filter_2d = recast_2Darray(filter_2d, output_shape, start_index_axis0, start_index_axis1)
        output = convolve_fft2D(data_2d,filter_2d)
    return output
############################################################################

ext: py

reason: import numpy,from scipy,from astropy


text:
import math
import copy
import numpy as np
from mecode.main import G

class GMatrix(G):
    """This class passes points through a 2D transformation matrix before
    fowarding them to the G class.  A 2D transformation matrix was
    choosen over a 3D transformation matrix because GCode's ARC
    command cannot be arbitrary rotated in a 3 dimensions.

    This lets you write code like:

    def box(g, height, width):
        g.move(0, width)
        g.move(height, 0)
        g.move(0, -width)
        g.move(-height, 0)

    def boxes(g, height, width):
        g.push_matrix()
        box(g, height, width)
        g.rotate(math.pi/8)
        box(g, height, width)
        g.pop_matrix()

    To get two boxes at a 45 degree angle from each other.

    The 2D transformation matrices are arranged in a stack,
    similar to OpenGL.

    numpy is required.

    """
    def __init__(self, *args, **kwargs):
        super(GMatrix, self).__init__(*args, **kwargs)
        self._matrix_setup()
        self.position_savepoints = []

    # Position savepoints #####################################################
    def save_position(self):
        self.position_savepoints.append((self.current_position["x"],
                                         self.current_position["y"],
                                         self.current_position["z"]))

    def restore_position(self):
        return_position = self.position_savepoints.pop()
        self.abs_move(return_position[0], return_position[1], return_position[2])

    # Matrix manipulation #####################################################
    def _matrix_setup(self):
        " Create our matrix stack. "
        self.matrix_stack = [np.matrix([[1.0, 0], [0.0, 1.0]])]

    def push_matrix(self):
        " Push a copy of our current transformation matrix. "
        self.matrix_stack.append(copy.deepcopy(self.matrix_stack[-1]))

    def pop_matrix(self):
        " Pop the matrix stack. "
        self.matrix_stack.pop()

    def rotate(self, angle):
        """Rotate the current transformation matrix around the Z
        axis, in radians. """
        rotation_matrix = np.matrix([[math.cos(angle), -math.sin(angle)],
                                     [math.sin(angle), math.cos(angle)]])

        self.matrix_stack[-1] = rotation_matrix * self.matrix_stack[-1]

    def scale(self, scale):
        " Scale the current transformation matrix. "
        scale_matrix = np.identity(2) * scale
        self.matrix_stack[-1] = scale_matrix * self.matrix_stack[-1]

    def _matrix_transform(self, x, y, z):
        "Transform an x,y,z coordinate by our transformation matrix."
        matrix = self.matrix_stack[-1]

        if x is None: x = 0
        if y is None: y = 0

        transform = matrix * np.matrix([x, y]).T

        return (transform.item(0), transform.item(1), z)

    def _matrix_transform_length(self, length):
        (x,y,z) = self._matrix_transform(length, 0, 0)
        return math.sqrt(x**2 + y**2 + z**2)

    def abs_move(self, x=None, y=None, z=None, **kwargs):
        if x is None: x = self.current_position['x']
        if y is None: y = self.current_position['y']
        if z is None: z = self.current_position['z']
        # abs_move ends up invoking move, which means that
        # we don't need to do a matrix transform here.
        super(GMatrix, self).abs_move(x,y,z, **kwargs)

    def move(self, x=None, y=None, z=None, **kwargs):
        (x,y,z) = self._matrix_transform(x,y,z)
        super(GMatrix, self).move(x,y,z, **kwargs)

    def arc(self, x=None, y=None, z=None, direction='CW', radius='auto',
            helix_dim=None, helix_len=0, **kwargs):
        (x_prime,y_prime,z_prime) = self._matrix_transform(x,y,z)
        if x is None: x_prime = None
        if y is None: y_prime = None
        if z is None: z_prime = None
        if helix_len: helix_len = self._matrix_transform_length(helix_len)
        super(GMatrix, self).arc(x=x_prime,y=y_prime,z=z_prime,direction=direction,radius=radius,
                                 helix_dim=helix_dim, helix_len=helix_len,
                                 **kwargs)
    @property
    def current_position(self):
        x = self._current_position['x']
        y = self._current_position['y']
        z = self._current_position['z']
        if x is None: x = 0.0
        if y is None: y = 0.0

        matrix = self.matrix_stack[-1]
        transform = matrix.getI() * np.matrix([x, y]).T

        return { 'x':transform.item(0),
                 'y':transform.item(1),
                 'z':z }

reason: import numpy


text:
import numpy as np

import pennylane as qml
from pennylane.ops.qubit import QubitStateVector
import itertools as it

import tensorflow as tf

tf.enable_eager_execution()

class RockyModel(tf.keras.Model):
    """
    QML model template.
    """

    def __init__(self, nclasses: int, device="default.qubit"):
        """
        Initialize the keras model interface.

        Args:
            nclasses: The number of classes in the data, used the determine the required output qubits.
            device: name of Pennylane Device backend.
        """
        super(RockyModel, self).__init__()
        self.req_qub_out = int(np.ceil(np.log2(nclasses)))
        self.req_qub_in = None
        self.device = device
        self.data_dev = None
        self.model_dev = None
        self.nclasses = nclasses
        self.init = False
        self.circuit = None
        self.trainable_vars = []

    def __str__(self):
        return "RockyModel Template"

    def initialize(self, nfeatures: int):
        """
        Model initialization.

        Args:
            nfeatures: The number of features in X

        """
        self.init = True
        raise NotImplementedError

    def call(self, inputs: tf.Tensor, observable: tf.Tensor):
        """
        Given some obsersable, we calculate the output of the model.

        Args:
            inputs: N x d tf.Tensor of N samples and d features.
            observable: tf.Tensor with Hermitian matrix containing an observable

        Returns: N expectation values of the observable

        """

        raise NotImplementedError

class RaccoonWrapper:
    """
    QML model training
    """

    # TODO convert everything to complex128?
    measurements = {
        "sx": tf.constant(np.array([[0, 1], [1, 0]]), dtype=tf.float64),
        "sy": tf.constant(
            np.array([[0, complex(0, -1)], [complex(0, 1), 0]]), dtype=tf.float64
        ),
        "sz": tf.constant(np.array([[1, 0], [0, -1]]), dtype=tf.float64),
        "id": tf.constant(np.array([[1, 0], [0, 1]]), dtype=tf.float64),
    }

    def __init__(self, model: RockyModel):
        """
        Wrapper allows one to minimize the quantum log-likelihood for a given RockyModel

        Args:
            model: The QML model we want to use for learning
        """

        self.Q = None
        self.q_y_x = None
        self.q_x = None
        self.req_measurements = None
        self.model = model
        self.bias = True

        def circuit(x, obs=None):
            QubitStateVector(x, wires=list(range(self.model.req_qub_out)))
            return qml.expval.Hermitian(obs, wires=list(range(self.model.req_qub_out)))

        self.data_circuit = qml.QNode(circuit, device=self.model.data_dev, cache=True)

    def _get_discr_statistics(self, X: np.ndarray, y: np.ndarray) -> None:
        """
        Calculate the empirical distribution q(y|x) for data X and labels y.

        Args:
            X: N x d matrix of N samples and d features.
            y: Length N vector with labels

        """
        y = y.flatten()
        assert X.shape[0] == y.shape[0], "Data has {} samples with {} labels".format(
            X.shape[0], X.shape[1]
        )

        classes = np.unique(y)
        nsamples, nfeatures = X.shape

        def _q_y_x(c: int, sample: np.ndarray) -> np.ndarray:
            # find copies of our sample
            _idx = np.where((X == tuple(sample)).all(axis=1))[0]
            # return probablity q(y=1|x)
            return np.sum(y[_idx] == c) / len(_idx)

        def _q_x(sample: np.ndarray) -> np.ndarray:
            # find copies of our sample
            _idx = np.where((X == tuple(sample)).all(axis=1))[0]
            # calcualte emperical probability of sample
            return len(_idx) / X.shape[0]

        # Initialize arrays
        self.q_y_x = np.zeros((nsamples, len(classes)))
        self.q_x = np.zeros((nsamples, 1))
        # Get statistics per class
        for i in range(nsamples):
            for j, cl in enumerate(classes):
                self.q_y_x[i, j] = _q_y_x(cl, X[i, :])
            self.q_x[i] = _q_x(X[i, :])

    def _determine_req_measurements(self, complex=False) -> None:
        """
        Determine which measurements are required for constructing the density matrix.

        Args:
            complex: Boolean that determines whether we care about SigmaY observables.

        Returns: list of requirement measurements.

        """

        measurements = []
        if complex:
            for m in it.product(
                ["sx", "sy", "sz", "id"], repeat=self.model.req_qub_out
            ):
                measurements.append(m)
        else:
            for m in it.product(["sx", "sz", "id"], repeat=self.model.req_qub_out):
                measurements.append(m)
        # remove the idenity measurement
        measurements.remove(tuple("id" for _ in range(self.model.req_qub_out)))
        # create constant tensor containing the required measurements.
        self.req_measurements = []
        for m in measurements:
            obs = 1
            for ob in m:
                obs = np.kron(obs, RaccoonWrapper.measurements[ob])
            self.req_measurements.append(obs)
        self.req_measurements = tf.constant(self.req_measurements, dtype=tf.float64)

    def construct_density_matrix(self, phi: np.ndarray) -> np.ndarray:
        """
        Construct the pure density matrix belonging to a certain wavefunction using a quantum circuit.

        Args:
            phi: Amplitudes of wavefunction (should be normed).

        Returns: A density matrix as ndarray.

        """
        # Fill the full wavefunction with the conditional probabilities q(x|y)
        full_phi = np.zeros((2 ** self.model.req_qub_out))
        full_phi[: len(phi)] = phi

        assert np.isclose(
            sum(np.abs(full_phi) ** 2), 1
        ), "Wavefunction must be normed to 1"
        # obtain the required expectation values to construct the density matrix from the data circuit.
        expval = (
            self.req_measurements
            * tf.map_fn(
                lambda x: self.data_circuit(full_phi, obs=x.numpy()),
                self.req_measurements,
                dtype=tf.float64,
            )[:, tf.newaxis, tf.newaxis]
        )
        # combine the tomography measurements
        eta = tf.reduce_sum(expval, axis=0)
        eta += tf.eye(
            2 ** self.model.req_qub_out, 2 ** self.model.req_qub_out, dtype=tf.float64
        )
        # normalization,
        eta /= tf.trace(eta)
        return eta

    def loss(self, inputs: np.ndarray, eta: tf.Tensor) -> tf.float64:
        """
        Determine the quantum log-likelihood of a batch of inputs, eta is the
        corresponding data density matrix

        Args:
            inputs: Mxd matrix of a batch of M samples with d features
            eta: Data density matrices for each sample so a rank 3 tensor.

        Returns: Scalar with quantum log-likelihood

        """
        assert self.model.init, "Initialize the model before calculating the loss"

        # obtain the required expectation values to construct the density matrix from the data circuit.
        expval = tf.map_fn(
            lambda x: self.model(inputs, x), self.req_measurements, dtype=tf.float64
        )
        expval = tf.transpose(expval)
        rho = tf.einsum("no,oij->nij", expval, self.req_measurements)
        rho += tf.eye(
            2 ** self.model.req_qub_out, 2 ** self.model.req_qub_out, dtype=tf.float64
        )
        # normalization
        rho /= tf.trace(rho)[:, tf.newaxis, tf.newaxis]
        log_rho = self._matrix_log(rho)
        # calculate quantum log-likelihood
        lh = -tf.trace(eta @ log_rho)
        # aggregate loss over all samples
        return tf.reduce_mean(lh)

    def train(self, X: np.ndarray, y: np.ndarray, epsilon=0.1, maxiter=100, tol=0.1):
        """
        Train the QML model with gradient descent.

        Args:
            X: N x d matrix of N samples and d features.
            y: Length N vector with labels.
            epsilon: Learning rate of the optimizer.
            maxiter: Maximum number of iterations before convergence.
            tol: Likelihood tolerance threshhold.

        """
        # add bias to inputs
        if self.model.bias:
            X = self.add_bias(X)

        assert self.model.nclasses == len(
            np.unique(y)
        ), "Model expects at {} classes, when y contains {} classes".format(
            self.model.nclasses, len(np.unique(y))
        )
        # determine the required measurements and get empirical statistics
        self._determine_req_measurements(False)
        self._get_discr_statistics(X, y)
        data_states = tf.map_fn(
            self.construct_density_matrix, np.sqrt(self.q_y_x), dtype=tf.float64
        )
        # multiply with qx (prior)
        data_states *= tf.reshape(self.q_x, (-1, 1, 1))
        # intialize the model
        self.model.initialize(X.shape[1])

        optimizer = tf.compat.v1.train.GradientDescentOptimizer(learning_rate=epsilon)
        self.lh = []
        # print training info
        print("-------TRAINING-------")
        print(self.model)
        print("{} samples, {} features".format(*X.shape))
        print(
            "{} input qubits, {} output qubits".format(
                self.model.req_qub_in, self.model.req_qub_out
            )
        )
        # training loop
        for i in range(maxiter):
            # calculate and apply gradients
            with tf.GradientTape() as tape:
                loss_value = self.loss(X, data_states)
                grads = tape.gradient(loss_value, self.model.trainable_vars)
            optimizer.apply_gradients(
                zip(grads, self.model.trainable_vars),
                global_step=tf.compat.v1.train.get_or_create_global_step(),
            )
            self.lh.append(float(loss_value))
            if i % 20 == 0:
                print("Loss at step {:03d}: {:.6f}".format(i, loss_value))
                if abs(loss_value - self.loss(X, data_states)) < tol:
                    print("L<{} after {:03d} iterations".format(tol, i))

                    break
        print("Final loss: {:.6f}".format(self.loss(X, data_states)))

    def predict(self, inputs: np.ndarray) -> np.ndarray:
        """

        Args:
            inputs: Mxd matrix of a batch of M samples with d features

        Returns: array of size M x nclasses with probabilities

        """

        assert self.model.init, "Initialize the model before predicting data"
        # add bias to inputs
        if self.model.bias:
            inputs = self.add_bias(inputs)
        # obtain the required expectation values to construct the density matrix from the data circuit.
        expval = tf.map_fn(
            lambda x: self.model(inputs, x), self.req_measurements, dtype=tf.float64
        )
        expval = tf.transpose(expval)
        # combine the tomography measurements
        rho = tf.einsum("no,oij->nij", expval, self.req_measurements)
        rho += tf.eye(
            2 ** self.model.req_qub_out, 2 ** self.model.req_qub_out, dtype=tf.float64
        )
        # normalization
        rho /= tf.trace(rho)[:, tf.newaxis, tf.newaxis]
        # stack probabilties for their respective classes
        return tf.stack(
            [tf.to_float(rho[:, i, i]) for i in range(rho.shape[1])], axis=1
        ).numpy()

    @staticmethod
    def _matrix_log(matrix: tf.Tensor) -> tf.Tensor:
        """
        Calculate matrix log2 through diagonalization of Hermitian matrix  M = U^-1 D U

        Args:
            matrix: N x N matrix

        Returns: N x N matrix

        """
        rx, Ux = tf.linalg.eigh(matrix)
        Ux_inv = tf.linalg.adjoint(Ux)
        rx = tf.cast(
            tf.math.log(tf.clip_by_value(tf.math.real(rx), 1e-13, 1e13)), rx.dtype
        )
        tx = tf.linalg.LinearOperatorDiag(rx).to_dense()
        return tf.matmul(Ux, tf.matmul(tx, Ux_inv))

    @staticmethod
    def add_bias(X: np.ndarray) -> np.ndarray:
        """
        Add a bias to the input data by adding a column of ones.

        Args:
            X: Numpy array of size N x d

        Returns: Numpy array of size N x d+1

        """
        n_samples = X.shape[0]
        return np.hstack([X, np.ones((n_samples, 1))])

max_stars_repo_name: therooler/pennylane-qllh

max_issues_repo_name: therooler/pennylane-qllh

max_forks_repo_name: therooler/pennylane-qllh


text:
import pandas as pd
import netCDF4
import argparse
import dateutil
import numpy as np
from os import path
import datetime as dt

def can_be_float(string):
    """Returns True if a string can be converted to a float,
    False if not"""
    try:
        float(string)
        return True
    except ValueError:
        return False

def extract_notes_from_lat(lat_val):
    """Extracts notes from a latitude column"""
    if can_be_float(lat_val):
        return None
    else:
        return lat_val

def build_datetime_column(df):
    """Creates a date column from year, month, day and time columns"""
    hour_column = (df['Time (UTC)']
                   .fillna(0)
                   .map(lambda dec_hour : pd.Timedelta(hours=np.floor(dec_hour))))
    minute_column = (df['Time (UTC)']
                     .fillna(0)
                     .map(lambda dec_hour : pd.Timedelta(minutes=(dec_hour%1)*60)))
    date_column = pd.to_datetime(df[['Day','Month','Year']])
    return date_column + hour_column + minute_column

def calc_lat_lon_count(df):
    lat_min = df.Latitude.min()
    lat_max = df.Latitude.max()
    lon_min = df.Longitude.min()
    lon_max = df.Longitude.max()
    precision = 0.001   # Hard-coding for now

    lat_count = (lat_max - lat_min)/precision
    lon_count = (lon_max - lon_min)/precision

    return lat_count, lon_count

def calc_time_axis(df):
    """Returns the number of minutes between midnight of the earliest date
    and midnight of the day after the latest date in the dataset"""
    time_min = df.Timestamp.min()
    time_max = df.Timestamp.max()

    diff_in_days =(time_max.date() - time_min.date())
    diff_in_minutes = (diff_in_days.days + 1)*24*60
    return diff_in_minutes

def make_new_corads_netcdf(df, path):
    ds = netCDF4.Dataset(path, 'w', clobber="True")

    ds.createDimension('n', len(df.index))
    ds.createDimension('str_len', 10)

    for column_name, data in df.iteritems():
        short_name = column_name.partition('/')[0].strip()

        if data.dtype == object:
            new_variable = ds.createVariable(short_name, 'S1', ('n', 'str_len'))
            new_variable._Encoding = 'ascii'
            new_variable[:] = data.array.astype('S10')
        else:
            new_variable = ds.createVariable(short_name, 'f8', ('n'))
            new_variable[:] = data.array

        new_variable.long_name = column_name
        if column_name.endswith('kg'):
            new_variable.units = "kg"
        if column_name in ["Latitude", "Longitude"]:
            new_variable.units = "Degrees"

    ds.close()

def clean_dataframe(df):
    df.rename(columns=lambda x: x.strip(), inplace=True)
    df['Time (UTC)'] = df['Time (UTC)'].where(df['Time (UTC)'].map(can_be_float)).astype(float)
    df['Timestamp'] = build_datetime_column(df)
    df['Notes'] = df['Latitude'].map(extract_notes_from_lat)
    df['Latitude'] = df['Latitude'].where(df['Latitude'].map(can_be_float)).astype(float)
    return df

def main(excel_path, out_dir):
    debris_df = pd.read_excel(excel_path, sheet_name='Re-entries')
    debris_df = clean_dataframe(debris_df)
    make_new_corads_netcdf(debris_df, path.join(out_dir, "reentries.nc4"))

    launch_df = pd.read_excel(excel_path, sheet_name='Launch emissions')
    launch_total_df = launch_df.tail(1)
    launch_df.drop(launch_df.tail(1).index,inplace=True)
    launch_df = clean_dataframe(launch_df)
    make_new_corads_netcdf(launch_df, path.join(out_dir, "launches.nc4"))

if __name__ == "__main__":
    parser = argparse.ArgumentParser("Converts your rocketry spreadsheet to netCDF")
    parser.add_argument("--excel_path", default = "Launch-Re-entry-emis-only-03-08.xlsx")
    parser.add_argument("--out_dir", default = '.')
    args = parser.parse_args()
    main(args.excel_path, args.out_dir)

max_stars_repo_name: eamarais/eam-group

max_issues_repo_name: eamarais/eam-group

max_forks_repo_name: eamarais/eam-group


text:
# Copyright (c) 2009-2022 The Regents of the University of Michigan.
# Part of HOOMD-blue, released under the BSD 3-Clause License.

"""Test the C++ internal wall data structures."""

import numpy as np
import pytest

from hoomd.md import _md

class _TestCounter:

    def __init__(self):
        self.previous_cls = None
        self.count = 0

    def __call__(self, arg):
        # Return empty string for non class arguments
        if isinstance(arg, dict):
            return ""
        if arg == self.previous_cls:
            self.count += 1
        self.count = 1
        self.previous_cls = arg
        return f"{arg.__name__}-{self.count}"

_test_args = (
    (_md.SphereWall, ({
        "radius": 4.0,
        "origin": (1.0, 0, 0),
        "inside": True,
        "open": False
    }, {
        "radius": 1.0,
        "origin": (0.0, 0.0, 0.0),
        "inside": False,
        "open": False
    }, {
        "radius": 3.1415,
        "origin": (-1.0, -2.0, 2.0),
        "inside": False,
        "open": True
    })),
    (_md.CylinderWall, ({
        "radius": 4.0,
        "origin": (1.0, 0, 0),
        "axis": (1.0, 0.0, 0),
        "inside": True,
        "open": True
    }, {
        "radius": 1.0,
        "origin": (0.0, 0.0, 0.0),
        "axis": (0.0, 1.0, 0.0),
        "inside": False,
        "open": False
    }, {
        "radius": 3.1415,
        "origin": (-1.0, -2.0, 1.0),
        "axis": (0.0, 0.0, 1.0),
        "inside": False,
        "open": True
    })),
    (
        _md.PlaneWall,
        (
            # The normals have to be unit vectors for the equality checks to
            # hold.  The C++ class currently normalizes any input vector.
            {
                "origin": (1.0, 0, 0),
                "normal": (1.0, 0.0, 0),
                "open": False
            },
            {
                "origin": (0.0, 0.0, 0.0),
                "normal": (0.0, 1.0, 0.0),
                "open": True
            },
            {
                "origin": (-1.0, -2.0, 1.0),
                "normal": (0.0, 0.0, 1.0),
                "open": False
            })))

@pytest.mark.parametrize("cls, constructor_kwargs",
                         ((cls, constructor_kwargs)
                          for cls, arg_list in _test_args
                          for constructor_kwargs in arg_list),
                         ids=_TestCounter())
def test_valid_construction(cls, constructor_kwargs):
    obj = cls(**constructor_kwargs)
    for key, value in constructor_kwargs.items():
        assert np.allclose(getattr(obj, key), value)

@pytest.mark.parametrize("cls, constructor_kwargs",
                         ((cls, arg_list[0]) for cls, arg_list in _test_args),
                         ids=_TestCounter())
def test_immutability(cls, constructor_kwargs):
    obj = cls(**constructor_kwargs)
    for key, value in constructor_kwargs.items():
        with pytest.raises(AttributeError):
            setattr(obj, key, value)

max_stars_repo_name: USF-GT-Molecular-Modeling/hoomd-blue

max_issues_repo_name: USF-GT-Molecular-Modeling/hoomd-blue

max_forks_repo_name: USF-GT-Molecular-Modeling/hoomd-blue


text:
"""
Collatz conjecture
https://github.com/FMoller/Racing_the_Dragon
"""

__author__ = 'Frederico Moeller'

import time
import numpy as np

def r_collatz(n,c_list):
    c_list.append(n)
    if n>1:
        if n%2==0:
            r_collatz(n/2,c_list)
        else:
            r_collatz(3*n+1,c_list)

def r_collatz_2(n,c_list):
    c_list.append(n)
    if n>1:
        if n%2==1:
            n = 3*n+1
            c_list.append(n)
        r_collatz(n/2,c_list)

def i_collatz(n,c_list):
    while(n>1):
        c_list.append(n)
        if n%2==0:
            n=n/2
        else:
            n=3*n+1
    c_list.append(n)

def i_collatz2(n,c_list):
    while(n>1):
        c_list.append(n)
        if n%2==1:
            n=3*n+1
            c_list.append(n)
        n=n/2
    c_list.append(n)

test_val = np.array([744, 340, 480, 679, 227, 996, 984, 753, 931, 679, 325, 939, 153,
       529, 535, 914, 804, 746, 453,  58, 668, 406, 217, 254, 830, 530,
       867, 983, 149, 306, 852, 398, 229, 736,   9, 406, 184, 377, 322,
       367, 960, 583, 355, 534, 980, 906, 474, 921, 895, 504, 586, 281,
        32, 371, 866, 172, 478, 906, 844, 140])

rcltz = []
rcltz2 = []
icltz = []
icltz2 = []

for i in test_val:
    c_list = []
    st_time = time.perf_counter()
    rcltz.append(r_collatz(i,c_list))
    ed_time = time.perf_counter()
    rcltz.append(ed_time - st_time)


text:
# Building a neural network model in TensorFlow
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
import tensorflow_datasets as tfds

# Building a linear regression model
X_train = np.arange(10).reshape((10, 1))
y_train = np.array([1.0, 1.3, 3.1,
                    2.0, 5.0, 6.3,
                    6.6, 7.4, 8.0,
                    9.0])

plt.plot(X_train, y_train, 'o', markersize=10)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

X_train_norm = (X_train - np.mean(X_train))/np.std(X_train)

ds_train_orig = tf.data.Dataset.from_tensor_slices(
    (tf.cast(X_train_norm, tf.float32),
     tf.cast(y_train, tf.float32)))

class MyModel(tf.keras.Model):

    def __init__(self):
        super(MyModel, self).__init__()
        self.w = tf.Variable(0.0, name='weight')
        self.b = tf.Variable(0.0, name='bias')

    def call(self, x):
        return self.w*x + self.b

model = MyModel()

model.build(input_shape=(None, 1))
model.summary()

def loss_fn(y_true, y_pred):
    return tf.reduce_mean(tf.square(y_true - y_pred))

## testing the function:
yt = tf.convert_to_tensor([1.0])
yp = tf.convert_to_tensor([1.5])

loss_fn(yt, yp)

def train(model, inputs, outputs, learning_rate):
    with tf.GradientTape() as tape:
        current_loss = loss_fn(model(inputs), outputs)
    dW, db = tape.gradient(current_loss, [model.w, model.b])
    model.w.assign_sub(learning_rate * dW)
    model.b.assign_sub(learning_rate * db)

tf.random.set_seed(1)

num_epochs = 200
log_steps = 100
learning_rate = 0.001
batch_size = 1
steps_per_epoch = int(np.ceil(len(y_train) / batch_size))

ds_train = ds_train_orig.shuffle(buffer_size=len(y_train))
ds_train = ds_train.repeat(count=None)
ds_train = ds_train.batch(1)

Ws, bs = [], []

for i, batch in enumerate(ds_train):
    if i >= steps_per_epoch * num_epochs:
        break
    Ws.append(model.w.numpy())
    bs.append(model.b.numpy())

    bx, by = batch
    loss_val = loss_fn(model(bx), by)

    train(model, bx, by, learning_rate=learning_rate)
    if i%log_steps==0:
        print('Epoch {:4d} Step {:2d} Loss {:6.4f}'.format(
              int(i/steps_per_epoch), i, loss_val))

print('Final Parameters:', model.w.numpy(), model.b.numpy())

X_test = np.linspace(0, 9, num=100).reshape(-1, 1)
X_test_norm = (X_test - np.mean(X_train)) / np.std(X_train)
y_pred = model(tf.cast(X_test_norm, dtype=tf.float32))

fig = plt.figure(figsize=(13, 5))
ax = fig.add_subplot(1, 2, 1)
plt.plot(X_train_norm, y_train, 'o', markersize=10)
plt.plot(X_test_norm, y_pred, '--', lw=3)
plt.legend(['Training examples', 'Linear Reg.'], fontsize=15)
ax.set_xlabel('x', size=15)
ax.set_ylabel('y', size=15)
ax.tick_params(axis='both', which='major', labelsize=15)

ax = fig.add_subplot(1, 2, 2)
plt.plot(Ws, lw=3)
plt.plot(bs, lw=3)
plt.legend(['Weight w', 'Bias unit b'], fontsize=15)
ax.set_xlabel('Iteration', size=15)
ax.set_ylabel('Value', size=15)
ax.tick_params(axis='both', which='major', labelsize=15)
plt.show()

# Model training via the .compile() and .fit() methods

tf.random.set_seed(1)
model = MyModel()

model.compile(optimizer='sgd',
              loss=loss_fn,
              metrics=['mae', 'mse'])

model.fit(X_train_norm, y_train,
          epochs=num_epochs, batch_size=batch_size,
          verbose=1)

print(model.w.numpy(), model.b.numpy())

X_test = np.linspace(0, 9, num=100).reshape(-1, 1)
X_test_norm = (X_test - np.mean(X_train)) / np.std(X_train)

y_pred = model(tf.cast(X_test_norm, dtype=tf.float32))

fig = plt.figure(figsize=(13, 5))
ax = fig.add_subplot(1, 2, 1)
plt.plot(X_train_norm, y_train, 'o', markersize=10)
plt.plot(X_test_norm, y_pred, '--', lw=3)
plt.legend(['Training Samples', 'Linear Regression'], fontsize=15)

ax = fig.add_subplot(1, 2, 2)
plt.plot(Ws, lw=3)
plt.plot(bs, lw=3)
plt.legend(['W', 'bias'], fontsize=15)
plt.show()

# Building a multilayer perceptron for classifying flowers in the Iris dataset

iris, iris_info = tfds.load('iris', with_info=True)

print(iris_info)

tf.random.set_seed(1)

ds_orig = iris['train']
ds_orig = ds_orig.shuffle(150, reshuffle_each_iteration=False)

print(next(iter(ds_orig)))

ds_train_orig = ds_orig.take(100)
ds_test = ds_orig.skip(100)

## checking the number of examples:

n = 0
for example in ds_train_orig:
    n += 1
print(n)

n = 0
for example in ds_test:
    n += 1
print(n)

ds_train_orig = ds_train_orig.map(
    lambda x: (x['features'], x['label']))

ds_test = ds_test.map(
    lambda x: (x['features'], x['label']))

next(iter(ds_train_orig))

iris_model = tf.keras.Sequential([
    tf.keras.layers.Dense(16, activation='sigmoid',
                          name='fc1', input_shape=(4,)),
    tf.keras.layers.Dense(3, name='fc2', activation='softmax')])

iris_model.summary()

iris_model.compile(optimizer='adam',
                   loss='sparse_categorical_crossentropy',
                   metrics=['accuracy'])

num_epochs = 100
training_size = 100
batch_size = 2
steps_per_epoch = np.ceil(training_size / batch_size)

ds_train = ds_train_orig.shuffle(buffer_size=training_size)
ds_train = ds_train.repeat()
ds_train = ds_train.batch(batch_size=batch_size)
ds_train = ds_train.prefetch(buffer_size=1000)

history = iris_model.fit(ds_train, epochs=num_epochs,
                         steps_per_epoch=steps_per_epoch,
                         verbose=0)

hist = history.history

fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(1, 2, 1)
ax.plot(hist['loss'], lw=3)
ax.set_title('Training loss', size=15)
ax.set_xlabel('Epoch', size=15)
ax.tick_params(axis='both', which='major', labelsize=15)

ax = fig.add_subplot(1, 2, 2)
ax.plot(hist['accuracy'], lw=3)
ax.set_title('Training accuracy', size=15)
ax.set_xlabel('Epoch', size=15)
ax.tick_params(axis='both', which='major', labelsize=15)
plt.tight_layout()
plt.show()

# Evaluating the trained model on the test dataset
results = iris_model.evaluate(ds_test.batch(50), verbose=0)
print('Test loss: {:.4f}   Test Acc.: {:.4f}'.format(*results))

# Saving and reloading the trained model
iris_model.save('iris-classifier.h5',
                overwrite=True,
                include_optimizer=True,
                save_format='h5')

iris_model_new = tf.keras.models.load_model('iris-classifier.h5')

iris_model_new.summary()

results = iris_model_new.evaluate(ds_test.batch(50), verbose=0)
print('Test loss: {:.4f}   Test Acc.: {:.4f}'.format(*results))

labels_train = []
for i,item in enumerate(ds_train_orig):
    labels_train.append(item[1].numpy())

labels_test = []
for i,item in enumerate(ds_test):
    labels_test.append(item[1].numpy())
print('Training Set: ',len(labels_train), 'Test Set: ', len(labels_test))

iris_model_new.to_json()

max_stars_repo_path: Chapter 13/03 Building a NN model in TensorFlow/program.py


text:
# coding:utf-8
# 基金风险收益图
from __future__ import print_function
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from fund.fund_info import get_fund_acc_net_value_by_time, norm, \
    fund_corr_coef, max_drawdown, growth_ratio
from os.path import dirname

__ROOT__ = dirname(__file__)

df = get_fund_acc_net_value_by_time()[-200:]
fund_bound = ['易方达高等级信用债债券A', '易方达纯债债券A', '嘉实稳祥纯债债券A', '嘉实超短债债券C', '招商产业债券A',
        '易方达裕丰回报债券','鹏扬利泽债券C', '招商双债增强债券(LOF)C']
fund_consumer = ['易方达消费行业股票', '汇添富中证主要消费ETF', '易方达消费精选股票',
                 '富国中证消费50ETF', '招商中证白酒A', '天弘中证食品饮料指数A']
fund_tech = ['华安创业板50ETF', '易方达创业板ETF', '易方达中证科技50ETF', '易方达上证科创板50成份ETF']
fund = fund_tech
df = df[fund]
df_basic = pd.DataFrame({
    'max_drawdown': df.apply(max_drawdown, axis=0),
    'growth_ratio': df.apply(growth_ratio, axis=0)
})
print(df_basic)
x,y = 'max_drawdown', 'growth_ratio'
ax = df_basic.plot(kind='scatter', x=x,y=y, figsize=(10,6))
for name, row in df_basic.iterrows():
    ax.annotate(name.decode('utf-8'), (row[x], row[y]))
plt.grid()
plt.gcf().savefig(__ROOT__ + '/image/plot_risk_return_tech.svg')

max_stars_repo_name: tracholar/economics-data

max_issues_repo_name: tracholar/economics-data

max_forks_repo_name: tracholar/economics-data


text:
from PyQt5.QtWidgets import QWidget, QGridLayout, QLabel, QLineEdit, QPushButton, QHeaderView, QAbstractItemView, \
                            QComboBox, QTableWidget, QTableWidgetItem, QMessageBox
from PyQt5.QtGui import QFont, QColor
from PyQt5.QtCore import Qt
from data_processing.build_model_data import AMADataBuilder, SVMDataBuilder
from data_processing.load_data import load_svm_modeldata, load_ama_modeldata, load_rawdata
from data_processing.svm_preprocessing import preprocess
from uis.calculate_ama import calculate_ama
from sklearn.svm import SVC,classes
from sklearn.neural_network import MLPRegressor
import yahoo_fin.stock_info as si
import numpy as np
import stockstats
import pickle
import os

class TradingSignal(QWidget):
    def __init__(self):
        super(TradingSignal, self).__init__()

        self.setFont(QFont('arial', 12, QFont.Bold))

        # 采用栅格布局
        main_layout = QGridLayout(self)

        # 首先还是读取数据
        # 读取现在有的etfs
        try:
            self.etfs = ['^FCHI',]#np.load('./Data/etfs.npy').tolist()
        except FileNotFoundError:
            self.etfs = ['^FCHI',]#['SPY', 'QQQ', 'TLT', 'GLD', 'IWM', 'EFA', 'HYG', 'XLV']
        self.signals_names = ['Weekly Over EMA12', 'Monthly Over EMA20', 'Weekly Over AMA', 'AMA Direction',
                              'SVM Probability', 'Conclusion']

        # 读取参数
        try:
            self.model_parameters = np.load('./Data/model_parameters.npy', allow_pickle=True).item()
            print(self.model_parameters)
        except FileNotFoundError:
            self.model_parameters = {'SVM': {'hist_window': 4},
                                     'AMA': {'price_type': 'close', 'er_window': 4, 'slow_window': 22,
                                             'fast_window': 4}}

        # 创建数据集需要的参数
        svm_label = QLabel('SVM')
        main_layout.addWidget(svm_label, 0, 0, 2, 1)

        svm_hist_window_label = QLabel('Back Window')
        main_layout.addWidget(svm_hist_window_label, 0, 1, 2, 1)

        self.svm_hist_window_input = QLineEdit()
        self.svm_hist_window_input.setText(str(self.model_parameters['SVM']['hist_window']))
        main_layout.addWidget(self.svm_hist_window_input, 0, 2, 2, 1)

        self.svm_btn = QPushButton('Build Data')
        self.svm_btn.clicked.connect(self.build_svm_data)
        main_layout.addWidget(self.svm_btn, 0, 3, 2, 1)

        svm_label = QLabel('AMA')
        main_layout.addWidget(svm_label, 2, 0, 1, 1)

        price_type_label = QLabel('Price Type')
        main_layout.addWidget(price_type_label, 3, 0, 2, 1)

        self.price_type_input = QComboBox()
        self.price_type_input.addItems(['open', 'close'])
        self.price_type_input.setCurrentText('open')
        main_layout.addWidget(self.price_type_input, 3, 1, 2, 1)

        er_window_label = QLabel('ER Window')
        main_layout.addWidget(er_window_label, 3, 2, 2, 1)

        self.er_window_input = QLineEdit()
        self.er_window_input.setText(str(self.model_parameters['AMA']['er_window']))
        main_layout.addWidget(self.er_window_input, 3, 3, 2, 1)

        fast_window_label = QLabel('Fast Window')
        main_layout.addWidget(fast_window_label, 3, 4, 2, 1)

        self.fast_window_input = QLineEdit()
        self.fast_window_input.setText(str(self.model_parameters['AMA']['fast_window']))
        main_layout.addWidget(self.fast_window_input, 3, 5, 2, 1)

        slow_window_label = QLabel('Slow Window')
        main_layout.addWidget(slow_window_label, 3, 6, 2, 1)

        self.slow_window_input = QLineEdit()
        self.slow_window_input.setText(str(self.model_parameters['AMA']['slow_window']))
        main_layout.addWidget(self.slow_window_input, 3, 7, 2, 1)

        self.ama_btn = QPushButton('Build Data')
        self.ama_btn.clicked.connect(self.build_ama_data)
        main_layout.addWidget(self.ama_btn, 3, 8, 2, 1)

        # 创建表格1的名称
        signals_label = QLabel('Trading Signals Table')
        main_layout.addWidget(signals_label, 5, 0, 2, 2)

        self.update_models_btn = QPushButton('Update Models')
        self.update_models_btn.clicked.connect(self.update_models)
        main_layout.addWidget(self.update_models_btn, 5, 2, 2, 2)

        self.show_signals_btn = QPushButton('Show Signals')
        self.show_signals_btn.clicked.connect(self.calculate_signals)
        main_layout.addWidget(self.show_signals_btn, 5, 4, 2, 2)

        # 创建表格
        self.signals_table = QTableWidget()
        self.signals_table.setRowCount(len(self.signals_names))
        self.signals_table.setColumnCount(len(self.etfs))
        # 设置放缩方式
        self.signals_table.horizontalHeader().setMinimumHeight(60)
        self.signals_table.horizontalHeader().setDefaultAlignment(Qt.AlignCenter)
        self.signals_table.horizontalHeader().setSectionResizeMode(QHeaderView.Stretch)
        self.signals_table.verticalHeader().setSectionResizeMode(QHeaderView.Stretch)
        self.signals_table.verticalHeader().setDefaultAlignment(Qt.AlignCenter)
        # self.horizontalHeader().setFixedHeight(80)
        # 设置默认行高
        self.signals_table.verticalHeader().setDefaultSectionSize(30)
        self.signals_table.horizontalHeader().setFont(QFont('arial', 12, QFont.Bold))
        # 设置行和列的名称
        self.signals_table.setHorizontalHeaderLabels(self.etfs)
        self.signals_table.setVerticalHeaderLabels(self.signals_names)
        # 禁止编辑
        self.signals_table.setEditTriggers(QAbstractItemView.NoEditTriggers)

        main_layout.addWidget(self.signals_table, 7, 0, 12, 9)

        self.setLayout(main_layout)

    def build_ama_data(self):
        price_type = self.price_type_input.currentText()
        er_window = int(self.er_window_input.text())
        slow_window = int(self.slow_window_input.text())
        fast_window = int(self.fast_window_input.text())
        ama_builder = AMADataBuilder(price_type, er_window, slow_window, fast_window)
        for ticker in self.etfs:
            for tf in ['weekly', 'monthly']:
                ama_builder.build_model_data(ticker, tf)

        # 存储参数
        self.model_parameters['AMA']['price_type'] = price_type
        self.model_parameters['AMA']['er_window'] = er_window
        self.model_parameters['AMA']['slow_window'] = slow_window
        self.model_parameters['AMA']['fast_window'] = fast_window
        np.save('./Data/model_parameters.npy', np.array(self.model_parameters))
        QMessageBox.information(self, 'Success', 'AMA model data built successfully!')

    def build_svm_data(self):
        look_hist_window = int(self.svm_hist_window_input.text())
        svm_builder = SVMDataBuilder(look_hist_window=look_hist_window)
        for ticker in self.etfs:
            for tf in ['weekly', 'monthly']:
                svm_builder.build_model_data(ticker, tf)
        self.model_parameters['SVM']['hist_window'] = look_hist_window
        QMessageBox.information(self, 'Success', 'SVM model data built successfully!')

    def update_models(self):
        print('update_models')
        # 检查需要的目录是否已经存在
        if not os.path.exists('./models/'):
            os.makedirs('./models/')
        svm_path = './models/svm_models/'
        if not os.path.exists(svm_path):
            os.makedirs(svm_path)
        ama_path = './models/ama_models/'
        if not os.path.exists(ama_path):
            os.makedirs(ama_path)

        # svm模型的建立
        for ticker in self.etfs:
            print("Ticker ----- ",ticker)
            for tf in ['weekly', 'monthly']:
                # 加载数据
                data = load_svm_modeldata(ticker, tf)
                # 数据预处理，只采用最近的300个数据
                #Elona
                #if len(data) > 300:
                 #   X, Y = preprocess(data.iloc[-301:-2])

                #else:
                X, Y = preprocess(data)
                # 创建分类器
                classifier = SVC(C=100, kernel='rbf', degree=3, gamma='scale', probability=True,
                                 class_weight='balanced')
                # 拟合分类器
                classifier.fit(X, Y)
                # 保存模型
                file_path_name = svm_path + ticker + '_' + tf + '_svm_model'
                file = open(file_path_name, 'wb')
                pickle.dump(classifier, file)
                file.close()

        # # ama模型的建立
        # for ticker in self.etfs:
        #     for tf in ['weekly', 'monthly']:
        #         # 加载数据
        #         data = load_ama_modeldata(ticker, tf)
        #         # 去掉空值
        #         data.dropna(how='any', inplace=True)
        #         data = np.array(data)
        #         # 将最后一列作为目标
        #         X = data[:-1, :-1]
        #         Y = data[:-1, -1]
        #
        #         # 建立多层回归模型， 两个隐藏层，隐藏层节点个数分别为6和4
        #         regressor = MLPRegressor(hidden_layer_sizes=(4, 6, 2), max_iter=3000, learning_rate='adaptive',
        #                                  learning_rate_init=0.02, epsilon=0.0001)
        #         # 拟合模型
        #         regressor.fit(X, Y)
        #         print(regressor.loss_)
        #         print(regressor.loss_curve_)
        #
        #         # 保存模型
        #         file_path_name = ama_path + ticker + '_' + tf + '_ama_model'
        #         file = open(file_path_name, 'wb')
        #         pickle.dump(regressor, file)
        #         file.close()
        QMessageBox.information(self, 'Success', 'Models updated successfully!')

    def calculate_signals(self):
        self.over_ema12 = []
        self.over_ema20 = []
        self.over_ama = []
        self.ama_direction = []
        self.svm_probability = []
        for i in range(len(self.etfs)):
            # 用weekly数据计算ema12
            data = load_rawdata(self.etfs[i], 'weekly')
            data_stats = stockstats.StockDataFrame.retype(data.copy())
            self.over_ema12.append(data_stats['close_12_ema'].iloc[-1])

            # 用weekly数据计算ama
            # # 首先创建ama的模型数据
            # data = load_ama_modeldata(self.etfs[i], 'weekly')
            # input = np.array(data)[-2:, 0:-1]
            #
            # # 添加ama_direction
            # self.ama_direction.append(input[-1, 1])
            #
            # # 加载模型
            # file = open('./models/ama_models/' + self.etfs[i] + '_weekly_ama_model', 'rb')
            # ama_model = pickle.load(file)
            # # 计算ama
            # self.over_ama.append(ama_model.predict(input)[-1])
            # file.close()

            ama, ama_direction, _ = calculate_ama(data, data_stats, self.model_parameters['AMA']['price_type'],
                                      self.model_parameters['AMA']['er_window'],
                                      self.model_parameters['AMA']['slow_window'],
                                      self.model_parameters['AMA']['fast_window'])
            self.over_ama.append(ama[-1])
            self.ama_direction.append(ama_direction[-1])

            # 用weekly数据计算svm
            # 首先创建svm需要的数据
            data = load_svm_modeldata(self.etfs[i], 'weekly')
            if len(data) > 300:
                X, Y = preprocess(data.iloc[-301:-2])
            else:
                X, Y = preprocess(data)
            input = X[-2:]
            # 加载模型
            file = open('./models/svm_models/' + self.etfs[i] + '_weekly_svm_model', 'rb')
            svm_model = pickle.load(file)

            self.svm_probability.append(svm_model.predict_proba(input)[-1, 1])
            file.close()

            # 用月数据计算ema20
            data = load_rawdata(self.etfs[i], 'monthly')
            data_stats = stockstats.StockDataFrame.retype(data.copy())
            self.over_ema20.append(data_stats['close_20_ema'].iloc[-1])

        self.show_signals()

    def show_signals(self, live_price=None):
        if live_price is None:
            live_price = []
            for etf in self.etfs:
                price = si.get_live_price(etf)
                price = round(price, 2)
                live_price.append(price)

        self.conclusion = []

        # 计算conclusion
        # 比例分配， ema12=0.3， ema20=0.1， ama=0.3， ama_direction=0.15, svm_proba=0.5
        for i in range(len(self.etfs)):
            if live_price[i] > self.over_ema12[i]:
                self.over_ema12[i] = 1
            else:
                self.over_ema12[i] = -1

            if live_price[i] > self.over_ema20[i]:
                self.over_ema20[i] = 1
            else:
                self.over_ema20[i] = -1

            if live_price[i] > self.over_ama[i]:
                self.over_ama[i] = 1
            else:
                self.over_ama[i] = -1

            self.svm_probability[i] = self.svm_probability[i] * 2 - 1

            self.ama_direction[i] = round(self.ama_direction[i], 2)
            self.svm_probability[i] = round(self.svm_probability[i], 2)

            self.conclusion.append(round(0.25 * self.over_ema12[i]
                                         + 0.25 * self.over_ama[i]
                                         + 0.2 * (self.ama_direction[i] + self.svm_probability[i])
                                         + 0.1 * self.over_ema20[i], 2))
            print('Over Ema ' , self.over_ema12[i]," Over ama ", self.over_ama[i],"ama dir svm ", (self.ama_direction[i] + self.svm_probability[i]),"over ema 20",self.over_ema20[i], 2)
            # 保存conclusion
            np.save('./Data/conclusion.npy', np.array(self.conclusion))

        # 将这些数据显示到表格中
        for i in range(len(self.etfs)):

            # over ema12
            new_item = QTableWidgetItem()
            if self.over_ema12[i] == 1:
                new_item.setText('Yes')
                new_item.setForeground(QColor(0, 200, 0))
            else:
                new_item.setText('No')
                new_item.setForeground(QColor(200, 0, 0))
            new_item.setTextAlignment(Qt.AlignCenter)
            new_item.setFont(QFont('arial', 12, QFont.Bold))
            self.signals_table.setItem(0, i, new_item)

            # over ema20
            new_item = QTableWidgetItem()
            if self.over_ema20[i] == 1:
                new_item.setText('Yes')
                new_item.setForeground(QColor(0, 200, 0))
            else:
                new_item.setText('No')
                new_item.setForeground(QColor(200, 0, 0))
            new_item.setTextAlignment(Qt.AlignCenter)
            new_item.setFont(QFont('arial', 12, QFont.Bold))
            self.signals_table.setItem(1, i, new_item)

            # over ama
            new_item = QTableWidgetItem()
            if self.over_ama[i] == 1:
                new_item.setText('Yes')
                new_item.setForeground(QColor(0, 200, 0))
            else:
                new_item.setText('No')
                new_item.setForeground(QColor(200, 0, 0))
            new_item.setTextAlignment(Qt.AlignCenter)
            new_item.setFont(QFont('arial', 12, QFont.Bold))
            self.signals_table.setItem(2, i, new_item)

            # ama direction
            new_item = QTableWidgetItem()
            new_item.setText(str(self.ama_direction[i]))
            if self.ama_direction[i] > 0.2:
                new_item.setForeground(QColor(0, 200, 0))
            elif self.ama_direction[i] < -0.2:
                new_item.setForeground(QColor(200, 0, 0))
            else:
                new_item.setForeground(QColor(200, 200, 0))
            new_item.setTextAlignment(Qt.AlignCenter)
            new_item.setFont(QFont('arial', 12, QFont.Bold))
            self.signals_table.setItem(3, i, new_item)

            # SVM
            new_item = QTableWidgetItem()
            new_item.setText(str(self.svm_probability[i]))
            if self.svm_probability[i] > 0.2:
                new_item.setForeground(QColor(0, 200, 0))
            elif self.svm_probability[i] < -0.2:
                new_item.setForeground(QColor(200, 0, 0))
            else:
                new_item.setForeground(QColor(200, 200, 0))
            new_item.setTextAlignment(Qt.AlignCenter)
            new_item.setFont(QFont('arial', 12, QFont.Bold))
            self.signals_table.setItem(4, i, new_item)

            # conclusion
            new_item = QTableWidgetItem()
            new_item.setText(str(self.conclusion[i]))
            if self.conclusion[i] > 0.5:
                new_item.setBackground(QColor(0, 200, 0))
                new_item.setForeground(QColor(0, 0, 0))
            elif self.conclusion[i] < -0.5:
                new_item.setBackground(QColor(200, 0, 0))
            else:
                new_item.setBackground(QColor(200, 200, 0))
                new_item.setForeground(QColor(0, 0, 0))
            new_item.setTextAlignment(Qt.AlignCenter)
            new_item.setFont(QFont('arial', 12, QFont.Bold))
            self.signals_table.setItem(5, i, new_item)

max_stars_repo_name: HusseinLezzaik/Stock-Market-Prediction

max_issues_repo_name: HusseinLezzaik/Stock-Market-Prediction

max_forks_repo_name: HusseinLezzaik/Stock-Market-Prediction


text:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

X = pd.read_csv("./Training Data/Linear_X_Train.csv").values
Y = pd.read_csv("./Training Data/Linear_Y_Train.csv").values

theta = np.load("ThetaList.npy")
# 100, 2
T0  = theta[:,0]
T1 = theta[:,1]

plt.ion()
for i in range(0,50,3):
    y_ = T1[i]*X + T0
    #points
    plt.scatter(X,Y)
    # line
    plt.plot(X,y_,'red')
    plt.draw()
    plt.pause(1)
    plt.clf()

max_stars_repo_path: ml_repo/3. Linear Regression/Linear Regression/04_interactive_plots.py


text:
import numpy as np
import scipy
import matplotlib
import matplotlib.pyplot as plt
import networkx as nx
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import shortest_path

from .DDRtree import DDRTree_py

def remove_velocity_points(G, n):
    for nodeid in range(n,2*n):
        nb_ids = []
        for nb_id in range(len(G[0])):
            if G[nodeid][nb_id]!=0:
                nb_ids = nb_ids+[nb_id]
        num_nbs = len(nb_ids)

        if num_nbs == 1:
            G[nodeid][nb_ids[0]] =0
            G[nb_ids[0]][nodeid] = 0
        else:
            min_val = np.inf
            for i in range(len(G[0])):
                if G[nodeid][i]!=0:
                    if G[nodeid][i]<min_val:
                        min_val = G[nodeid][i]
                        min_ind = i
            for i in nb_ids:
                if i != min_ind:
                    new_weight = G[nodeid][i] + min_val
                    G[i][min_ind] = new_weight
                    G[min_ind][i] = new_weight
            # print('Add ege %s, %s\n',G.Nodes.Name {nb_ids(i)}, G.Nodes.Name {nb_ids(min_ind)});
            G[nodeid][nb_ids[0]] = 0
            G[nb_ids[0]][nodeid] = 0

    return G

def calculate_angle(o, y, x):
    yo = y - o
    norm_yo = yo/ scipy.linalg.norm(yo)
    xo = x - o
    norm_xo = xo/ scipy.linalg.norm(xo)
    angle = np.arccos(norm_yo.T * norm_xo)
    return angle

def construct_velocity_tree_py(X1,X2):

    n = X1.shape[1]

    # merge two data with a given time
    t = 0.5
    X_all = np.hstack((X1, X1 + t* X2))

    # parameter settings
    maxIter = 20
    eps = 1e-3
    sigma = 0.001
    gamma = 10

    # run DDRTree algorithm
    W, Z, stree, Y, R, history = DDRTree_py(X_all, maxIter = maxIter,eps = eps,sigma = sigma,gamma = gamma)

    # draw velocity figure

    # quiver(Z(1, 1: 100), Z(2, 1: 100), Z(1, 101: 200)-Z(1, 1: 100), Z(2, 101: 200)-Z(2, 1: 100));
    # plot(Z(1, 1: 100), Z(2, 1: 100), 'ob');
    # plot(Z(1, 101: 200), Z(2, 101: 200), 'sr');
    G = stree

    sG = remove_velocity_points(G, n)
    tree = sG
    row = []
    col = []
    val = []
    for i in range(sG.shape[0]):
        for j in range(sG.shape[1]):
            if sG[i][j]!=0:
                row = row+[i]
                col = col+[j]
                val = val + [sG[1][j]]
    tree_fname = 'tree.csv'
    # write sG data to tree.csv
    #######
    branch_fname = 'branch.txt'
    cmd = 'python extract_branches.py'+tree_fname + branch_fname

    branch_cell = []
    fid = open(branch_fname,'r')
    tline = next(fid)
    while isinstance(tline,str):
        path = regexp(tline, '\d*', 'Match')############
        branch_cell = branch_cell+[path]#################
        tline = next(fid)
    fid.close()

    dG = np.zeros((n, n))
    for p in range(len(branch_cell)):
        path = branch_cell[p]
        pos_direct = 0
        for bp in range(len(path)):
            u = path(bp)
            v = u + n

            # find the shorest path on graph G(works for trees)
            nodeid = u
            ve_nodeid = v
            shortest_mat = shortest_path(csgraph=G, directed=False, indices=nodeid, return_predecessors=True)
            velocity_path = []
            while ve_nodeid!=nodeid:
                velocity_path = [shortest_mat[nodeid][ve_nodeid]]+velocity_path
                ve_nodeid = shortest_mat[nodeid][ve_nodeid]
            velocity_path = [shortest_mat[nodeid][ve_nodeid]] + velocity_path
            ###v_path = G.Nodes.Name(velocity_path)

            # check direction consistency between path and v_path
            valid_idx = []
            for i in velocity_path:
                if i<=n:
                    valid_idx = valid_idx+[i]
            if len(valid_idx) == 1:
                # compute direction matching
                if bp < len(path):
                    tree_next_point = Z[:, path(bp)]
                    v_point = Z[:, v]
                    u_point = Z[:, u]
                    angle = calculate_angle(u_point, tree_next_point, v_point)
                    angle = angle / 3.14 * 180
                    if angle < 90:
                        pos_direct = pos_direct + 1

                else:
                    tree_pre_point = Z[:, path(bp - 1)]
                    v_point = Z[:, v]
                    u_point = Z[:, u]
                    angle = calculate_angle(u_point, tree_pre_point, v_point)
                    angle = angle / 3.14 * 180
                    if angle > 90:
                        pos_direct = pos_direct + 1

            else:

                if bp < len(path):
                    if path[bp + 1] == valid_idx[2]:
                        pos_direct = pos_direct + 1

                else:
                    if path[bp - 1] != valid_idx[2]:
                        pos_direct = pos_direct + 1

        neg_direct = len(path) - pos_direct
        print('branch='+ str(p) +', ('+path[0]+'->'+path[-1]+'), pos='+pos_direct+', neg='+neg_direct+'\n')
        print(path)
        print('\n')

        if pos_direct > neg_direct:
            for bp in range(len(path) - 1):
                dG[path[bp], path[bp + 1]] = 1

        else:
            for bp in range(len(path) - 1):
                dG[path(bp + 1), path(bp)] = 1

    # figure;
    # plot(digraph(dG));
    # title('directed graph') figure; hold on;
    row = []
    col = []
    for i in range(dG.shape[0]):
        for j in range(dG.shape[1]):
            if dG[i][j]!=0:
                row = row+[i]
                col = col+[j]
    for tn in range(len(row)):
        p1 = Y[:, row[tn]]
        p2 = Y[:, col[tn]]
        dp = p2 - p1
        h = quiver(p1(1), p1(2), dp(1), dp(2), 'LineWidth', 5)###############need to plot it
        set(h, 'MaxHeadSize', 1e3, 'AutoScaleFactor', 1)#############

    for i in range(n):
        text(Y(1, i), Y(2, i), str(i))##############
    savefig('./results/t01_figure3.fig')##################

max_stars_repo_name: softbear/dynamo-release

max_issues_repo_name: softbear/dynamo-release

max_forks_repo_name: softbear/dynamo-release

reason: import numpy,import scipy,from scipy,import networkx


text:
from convokit import Corpus, CorpusComponent
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from typing import Callable, List
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
import pandas as pd
from .classifier import Classifier
import numpy as np
from .util import extract_vector_feats_and_label

class VectorClassifier(Classifier):
    """
    Transformer that trains a classifier on the Corpus components' text vector representation
    (e.g. bag-of-words, TF-IDF, etc)

    Corpus must have a vector with the specified `vector_name`.

    Inherits from `Classifier` and has access to its methods.

    :param obj_type: "speaker", "utterance", or "conversation"
    :param vector_name: the metadata key where the Corpus object text vector is stored
    :param columns: list of column names of vector matrix to use; uses all columns by default.
    :param labeller: a (lambda) function that takes a Corpus object and returns True (y=1) or False (y=0) - i.e.
        labeller defines the y value of the object for fitting
    :param clf: a sklearn Classifier. By default, clf is a Pipeline with StandardScaler and LogisticRegression
    :param clf_attribute_name: the metadata attribute name to store the classifier prediction value under; default:
        "prediction"
    :param clf_prob_attribute_name: the metadata attribute name to store the classifier prediction score under;
        default: "pred_score"
    """
    def __init__(self, obj_type: str, vector_name: str, columns: List[str] = None,
                 labeller: Callable[[CorpusComponent], bool] = lambda x: True,
                 clf=None, clf_attribute_name: str = "prediction", clf_prob_attribute_name: str = "pred_score"):
        if clf is None:
            clf = Pipeline([("standardScaler", StandardScaler(with_mean=False)),
                            ("logreg", LogisticRegression(solver='liblinear'))])
            print("Initialized default classification model (standard scaled logistic regression).")

        super().__init__(obj_type=obj_type, pred_feats=[], labeller=labeller,
                         clf=clf, clf_attribute_name=clf_attribute_name, clf_prob_attribute_name=clf_prob_attribute_name)
        self.vector_name = vector_name
        self.columns = columns

    def fit(self, corpus: Corpus, selector: Callable[[CorpusComponent], bool] = lambda x: True, y=None):
        """
        Fit the Transformer's internal classifier model on the vector matrix that represents one of
        the Corpus components, with an optional selector that selects for objects to be fit on.

        :param corpus: the target Corpus
        :param selector: a (lambda) function that takes a Corpus object and returns True or False
            (i.e. include / exclude). By default, the selector includes all objects of the specified type in the Corpus.
        :return: the fitted VectorClassifier
        """
        # collect texts for vectorization
        obj_ids = []
        y = []
        for obj in corpus.iter_objs(self.obj_type, selector):
            obj_ids.append(obj.id)
            y.append(self.labeller(obj))
        X = corpus.get_vectors(self.vector_name, ids=obj_ids, columns=self.columns)
        y = np.array(y)
        # print(corpus.get_vector_matrix(self.vector_name).matrix.shape)
        # print(X.shape)
        # print(y.shape)
        self.clf.fit(X, y)
        return self

    def transform(self, corpus: Corpus, selector: Callable[[CorpusComponent], bool] = lambda x: True) -> Corpus:
        """
        Annotate the corpus components with the classifier prediction and prediction score, with an optional selector
        that selects for objects to be classified. Objects that are not selected will get a metadata value of 'None'
        instead of the classifier prediction.

        :param corpus: the target Corpus
        :param selector: a (lambda) function that takes a Corpus object and returns True or False
            (i.e. include / exclude). By default, the selector includes all objects of the specified type in the Corpus.

        :return: the target Corpus annotated
        """
        objs = []
        for obj in corpus.iter_objs(self.obj_type):
            if selector(obj):
                objs.append(obj)
            else:
                obj.add_meta(self.clf_attribute_name, None)
                obj.add_meta(self.clf_prob_attribute_name, None)

        obj_ids = [obj.id for obj in objs]
        X = corpus.get_vector_matrix(self.vector_name).get_vectors(obj_ids, self.columns)

        clfs, clfs_probs = self.clf.predict(X), self.clf.predict_proba(X)[:, 1]

        for idx, (clf, clf_prob) in enumerate(list(zip(clfs, clfs_probs))):
            obj = objs[idx]
            obj.add_meta(self.clf_attribute_name, clf)
            obj.add_meta(self.clf_prob_attribute_name, clf_prob)
        return corpus

    def transform_objs(self, objs: List[CorpusComponent]) -> List[CorpusComponent]:
        """
        Not implemented for VectorClassifier.
        """
        # """
        # Run classifier on list of Corpus component objects and annotate them with their predictions and
        # prediction scores.
        #
        # :param objs: list of Corpus objects
        # :return: list of annotated Corpus objects
        # """
        raise NotImplementedError("VectorClassifier can only be run on corpora, not arbitrary lists of corpus "
                                  "component objects.")

    def fit_transform(self, corpus: Corpus, selector: Callable[[CorpusComponent], bool] = lambda x: True) -> Corpus:
        """
        Runs the `fit()` and `transform()` steps in order, with the specified selector.

        :param corpus: the target Corpus
        :param selector: a (lambda) function that takes a Corpus object and returns True or False
            (i.e. include / exclude). By default, the selector includes all objects of the specified type in the Corpus.

        :return: the target Corpus annotated
        """
        self.fit(corpus, selector=selector)
        return self.transform(corpus, selector=selector)

    def summarize(self, corpus: Corpus, selector: Callable[[CorpusComponent], bool] = lambda x: True):
        """
        Generate a DataFrame indexed by object id with the classifier predictions and scores.

        :param corpus: the annotated Corpus
        :param selector: a (lambda) function that takes a Corpus object and returns True or False
            (i.e. include / exclude). By default, the selector includes all objects of the specified type in the Corpus.
        :return: a pandas DataFrame
        """
        objId_clf_prob = []

        for obj in corpus.iter_objs(self.obj_type, selector):
            objId_clf_prob.append((obj.id, obj.meta[self.clf_attribute_name], obj.meta[self.clf_prob_attribute_name]))

        return pd.DataFrame(list(objId_clf_prob),
                            columns=['id', self.clf_attribute_name, self.clf_prob_attribute_name])\
                            .set_index('id').sort_values(self.clf_prob_attribute_name, ascending=False)

    def summarize_objs(self, objs: List[CorpusComponent]):
        """
        Not implemented for VectorClassifier.
        """
        # Generate a pandas DataFrame (indexed by object id, with prediction and prediction score columns) of
        # classification results.
        #
        # Runs on a list of Corpus objects.
        #
        # :param objs: list of Corpus objects
        # :return: pandas DataFrame indexed by Corpus object id
        # objId_clf_prob = []
        # for obj in objs:
        #     objId_clf_prob.append((obj.id, obj.meta[self.clf_attribute_name], obj.meta[self.clf_prob_attribute_name]))
        #
        # return pd.DataFrame(list(objId_clf_prob),
        #                     columns=['id', self.clf_attribute_name, self.clf_prob_attribute_name]).set_index('id').sort_values(
        #                     self.clf_prob_attribute_name)
        raise NotImplementedError("VectorClassifier can only be run on corpora, not arbitrary lists of corpus "
                                  "component objects.")

    def evaluate_with_train_test_split(self, corpus: Corpus,
                                       selector: Callable[[CorpusComponent], bool] = lambda x: True,
                                       test_size: float = 0.2):
        """
        Evaluate the performance of predictive features (Classifier.pred_feats) in predicting for the label,
        using a train-test split.

        Run either on a Corpus (with Classifier labeller, selector, obj_type settings) or a list of Corpus objects

        :param corpus: target Corpus
        :param selector: a (lambda) function that takes a Corpus object and returns True or False (i.e. include /
            exclude). By default, the selector includes all objects of the specified type in the Corpus.
        :param test_size: size of test set
        :return: accuracy and confusion matrix
        """
        X, y = extract_vector_feats_and_label(corpus, self.obj_type, self.vector_name, self.columns,
                                              self.labeller, selector)

        print("Running a train-test-split evaluation...")
        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size)
        self.clf.fit(X_train, y_train)
        preds = self.clf.predict(X_test)
        accuracy = np.mean(preds == y_test)
        print("Done.")
        return accuracy, confusion_matrix(y_true=y_test, y_pred=preds)

    def evaluate_with_cv(self, corpus: Corpus, cv=KFold(n_splits=5, shuffle=True),
                         selector: Callable[[CorpusComponent], bool] = lambda x: True):
        """
        Evaluate the performance of predictive features (Classifier.pred_feats) in predicting for the label,
        using cross-validation for data splitting.

        :param corpus: target Corpus
        :param cv: cross-validation model to use: KFold(n_splits=5, shuffle=True) by default.
        :param selector: if running on a Corpus, this is a (lambda) function that takes a Corpus object and returns True or False (i.e. include / exclude). By default, the selector includes all objects of the specified type in the Corpus.

        :return: cross-validated accuracy score
        """

        X, y = extract_vector_feats_and_label(corpus, self.obj_type, self.vector_name,
                                              self.columns, self.labeller, selector)
        print("Running a cross-validated evaluation...", end="")
        score = cross_val_score(self.clf, X, y, cv=cv)
        print("Done.")
        return score


text:
# -*- encoding: utf-8 -*-
# test_coding.py - unit tests for the coding module
#
# Copyright 2015, 2016 NetworkX developers.
#
# This file is part of NetworkX.
#
# NetworkX is distributed under a BSD license; see LICENSE.txt for more
# information.
"""Unit tests for the :mod:`~networkx.algorithms.tree.coding` module."""
from itertools import product

from nose.tools import assert_equal
from nose.tools import assert_true
from nose.tools import raises

import networkx as nx

class TestPruferSequence(object):
    """Unit tests for the Prüfer sequence encoding and decoding
    functions.

    """

    @raises(nx.NotATree)
    def test_nontree(self):
        G = nx.cycle_graph(3)
        nx.to_prufer_sequence(G)

    @raises(nx.NetworkXPointlessConcept)
    def test_null_graph(self):
        nx.to_prufer_sequence(nx.null_graph())

    @raises(nx.NetworkXPointlessConcept)
    def test_trivial_graph(self):
        nx.to_prufer_sequence(nx.trivial_graph())

    @raises(KeyError)
    def test_bad_integer_labels(self):
        T = nx.Graph(nx.utils.pairwise('abc'))
        nx.to_prufer_sequence(T)

    def test_encoding(self):
        """Tests for encoding a tree as a Prüfer sequence using the
        iterative strategy.

        """
        # Example from Wikipedia.
        tree = nx.Graph([(0, 3), (1, 3), (2, 3), (3, 4), (4, 5)])
        sequence = nx.to_prufer_sequence(tree)
        assert_equal(sequence, [3, 3, 3, 4])

    def test_decoding(self):
        """Tests for decoding a tree from a Prüfer sequence."""
        # Example from Wikipedia.
        sequence = [3, 3, 3, 4]
        tree = nx.from_prufer_sequence(sequence)
        assert_equal(list(tree), list(range(6)))
        edges = [(0, 3), (1, 3), (2, 3), (3, 4), (4, 5)]
        assert_equal(list(tree.edges()), edges)

    def test_decoding2(self):
        # Example from "An Optimal Algorithm for Prufer Codes".
        sequence = [2, 4, 0, 1, 3, 3]
        tree = nx.from_prufer_sequence(sequence)
        assert_equal(list(tree), list(range(8)))
        edges = [(0, 1), (0, 4), (1, 3), (2, 4), (2, 5), (3, 6), (3, 7)]
        assert_equal(list(tree.edges()), edges)

    def test_inverse(self):
        """Tests that the encoding and decoding functions are inverses.

        """
        for T in nx.nonisomorphic_trees(4):
            T2 = nx.from_prufer_sequence(nx.to_prufer_sequence(T))
            assert_equal(list(T), list(T2))
            assert_equal(list(T.edges()), list(T2.edges()))

        for seq in product(range(4), repeat=2):
            seq2 = nx.to_prufer_sequence(nx.from_prufer_sequence(seq))
            assert_equal(list(seq), seq2)

class TestNestedTuple(object):
    """Unit tests for the nested tuple encoding and decoding functions.

    """

    @raises(nx.NotATree)
    def test_nontree(self):
        G = nx.cycle_graph(3)
        nx.to_nested_tuple(G, 0)

    @raises(nx.NodeNotFound)
    def test_unknown_root(self):
        G = nx.path_graph(2)
        nx.to_nested_tuple(G, 'bogus')

    def test_encoding(self):
        T = nx.full_rary_tree(2, 2 ** 3 - 1)
        expected = (((), ()), ((), ()))
        actual = nx.to_nested_tuple(T, 0)
        assert_equal(expected, actual)

    def test_canonical_form(self):
        T = nx.Graph()
        T.add_edges_from([(0, 1), (0, 2), (0, 3)])
        T.add_edges_from([(1, 4), (1, 5)])
        T.add_edges_from([(3, 6), (3, 7)])
        root = 0
        actual = nx.to_nested_tuple(T, root, canonical_form=True)
        expected = ((), ((), ()), ((), ()))
        assert_equal(actual, expected)

    def test_decoding(self):
        balanced = (((), ()), ((), ()))
        expected = nx.full_rary_tree(2, 2 ** 3 - 1)
        actual = nx.from_nested_tuple(balanced)
        assert_true(nx.is_isomorphic(expected, actual))

    def test_sensible_relabeling(self):
        balanced = (((), ()), ((), ()))
        T = nx.from_nested_tuple(balanced, sensible_relabeling=True)
        edges = [(0, 1), (0, 2), (1, 3), (1, 4), (2, 5), (2, 6)]
        assert_equal(list(T), list(range(2 ** 3 - 1)))
        assert_equal(list(T.edges()), edges)

max_stars_repo_name: KyleBenson/networkx

max_issues_repo_name: KyleBenson/networkx

max_forks_repo_name: KyleBenson/networkx

reason: import networkx