text:
import os
import time
# ignoring TF info messages
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'

import random
import numpy as np
import tensorflow as tf

import argparse

from data_loader import load_dataset
from train import train
from test import test

DATASETS = ['cardio', 'mammo', 'satellite', 'seismic', 'annthyroid', 'thyroid', 'vowels', 'yeast']

def get_execute_time(start_time, end_time):
    hours, rem = divmod(end_time - start_time, 3600)
    minutes, seconds = divmod(rem, 60)
    print("{:0>2}:{:0>2}:{:05.2f} ---".format(int(hours), int(minutes), seconds))

def seed_everything(seed):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    tf.random.set_seed(seed)

def run_experiment(args):
    """
    Running a specific experiment with the given args
    Parameters
    ----------
    args: argparse args. The args given to the program
    """

    # setting seed
    seed_everything(1234)

    # loading the preprocessed dataset
    print(f"--- Loading preprocessed {args.dataset_name.capitalize()} dataset ---")
    start_time = time.time()
    dataset_X, dataset_y, siamese_pairs, features_dim, folded_train_datasets_list, folded_test_datasets_list = load_dataset(args)
    end_time = time.time()
    print(f"--- {args.dataset_name.capitalize()} dataset ready after: ", end='')
    get_execute_time(start_time, end_time)

    # training
    print("--- Start training ---")
    start_time = time.time()
    trained_estimators_list, trained_siamese_network, euclidean_nn_model, siamese_nn_model = train(dataset_X, dataset_y, siamese_pairs, folded_train_datasets_list, features_dim, args)
    end_time = time.time()
    print("--- Training finished after: ", end='')
    get_execute_time(start_time, end_time)

    # testing
    print("--- Start testing ---")
    start_time = time.time()
    test(dataset_X, folded_test_datasets_list, trained_estimators_list, trained_siamese_network, euclidean_nn_model, siamese_nn_model, args)
    end_time = time.time()
    print("--- Testing finished after: ", end='')
    get_execute_time(start_time, end_time)
    print(f"--- Finished dataset {args.dataset_name.capitalize()} ---")

if __name__ == '__main__':

    # getting the dataset to preprocess
    parser = argparse.ArgumentParser()
    parser.add_argument('-d', "--dataset", required=True, dest='dataset_name', type=str, help='the dataset to preprocess and save to disk for later use')
    parser.add_argument("-n", "--neighbors", dest="num_neighbors", default=10, type=int, help="The number of neighbors to retrieve from the Neareset Neighbors model")
    parser.add_argument("-a", "--augmentations", dest="num_augmentations", default=7, type=int, help="The number test-time augmentations to apply on every test sample")
    parser.add_argument("-c", "--cuml", dest="with_cuml", default=True, type=bool, help='Whether of not to use cuML')
    parser.add_argument("-s", "--siamesebatchsize", dest="siamese_batch_size", default=64, type=int, help="Batch size used to train the Siamese network")
    parser.add_argument("-b", "--siameseecpohs", dest="siamese_n_epochs", default=10, type=int, help="Number of epochs to train with the Siamese netowrk")
    parser.add_argument("-e", "--epochs", dest="num_epochs", default=500, type=int, help="The numberof epochs to train with the anomaly detector model")
    parser.add_argument("-f", "--folds", dest="n_folds", default=10, type=int, help="The number of folds in the K-Fold cross-validation")
    args = parser.parse_args()

    # adding support for runnning all of the available datasets
    if args.dataset_name == "all":
        dataset_index = int(os.environ['SLURM_ARRAY_TASK_ID'])-1
        args.dataset_name = DATASETS[dataset_index]

    # running experiments with the given program arguments
    run_experiment(args)

max_stars_repo_name: nivgold/TTAD

max_issues_repo_name: nivgold/TTAD

max_forks_repo_name: nivgold/TTAD


text:
#!/usr/bin/env python
''' Simulation class for the HEMS System

Inputs:
Outputs:
'''

__author__ = "D. Esteban M. Bondy, based on work by Tue V. Jensen"
__copyright__ = "Copyright 2018, DTU"
__license__ = "MIT"
__version__ = "0.1"
__maintainer__ = "Esteban Bondy"
__email__ = "bondy@elektro.dtu.dk"
__status__ = "Development"

import mosaik
import mosaik.util
import pandas as pd
import numpy as np
import os
# from dtu_mosaik import *

data_path = 'temp_files/'

def run_simulation(database_name, run_data, verbose=True):
    SIM_CONFIG = {
        'DemandModel': {
            'python': 'dtu_mosaik.mosaik_demand:DemandModel',
        },
        'SimpleGridModel': {
            'python': 'dtu_mosaik.mosaik_grid:SimpleGridModel',
        },
        'CollectorSim': {
            'python': 'dtu_mosaik.collector:Collector',
        },
        'PVModel': {
            'python': 'dtu_mosaik.mosaik_pv:PVModel'
        },
        'BatteryModel': {
            'python': 'dtu_mosaik.mosaik_battery:BatteryModel'
        },
        'Control': {
            'python': 'dtu_mosaik.mosaik_control:Control'
        },
        #'WTModel': {
        #    'python': 'dtu_mosaik.mosaik_wt:WTModel'
        #}
        'NoiseGenerator': {
            'python': 'dtu_mosaik.mosaik_noise:NoiseGenerator'
        }
    }
    run_id = run_data['ID']
    pv1_cap = run_data['pv1_scaling']
    battery_cap = run_data['batt_storage_capacity']
    battery_rate = run_data['batt_charge_capacity']
    change_rate = run_data['controller_change_rate']
    day_type = run_data['climate_conditions']
    random_weather = run_data['random_weather']
    stochastic = run_data['stochastic']
    noise_scale = run_data['noise_scale']
    season = run_data['season']

    seasons = {'summer': 1, 'autumn': 3, 'winter': 5, 'spring': 2}
    demand = seasons[run_data['season']]

    weather_base = {'cloudy': ['/PV715_20180125', '/PV715_20180126', '/PV715_20180127', '/PV715_20180130'],
                    'intermittent': ['/PV715_20180423', '/PV715_20180430', '/PV715_20180820', '/PV715_20180722'],
                    'sunny': ['/PV715_20180730', '/PV715_20180728', '/PV715_20180729', '/PV715_20180721']}
    if random_weather:
        day = weather_base[day_type][np.random.randint(0, 4)]
    else:
        day = weather_base[day_type][0]

    def init_entities(world, pv1_rated_capacity=pv1_cap, batt1_rated_cap=battery_cap,
                      batt1_rate=battery_rate, con_change_rate=change_rate, weather=day, base_demand=demand,
                      noise_scale=noise_scale, filename=data_path+database_name):
        sim_dict = {}
        entity_dict = {}

        # quick check if filepath exists
        directory = os.path.dirname(filename)
        if not os.path.exists(directory):
            os.makedirs(directory)

        ## Demand
        demand_sim = world.start(
            'DemandModel',
            eid_prefix='demand_',
            step_size=1)
        demand_entity_1 = demand_sim.DemandModel(rated_capacity=base_demand, seriesname='/flexhouse_20180219')
        sim_dict['demand'] = demand_sim
        entity_dict['demand1'] = demand_entity_1

        ## Grid model
        grid_sim = world.start(
            'SimpleGridModel',
            eid_prefix='grid_',
            step_size=1)
        grid_entity_1 = grid_sim.SimpleGridModel(V0=240, droop=0.1)
        sim_dict['grid'] = grid_sim
        entity_dict['grid1'] = grid_entity_1

        ## Collector
        collector_sim = world.start(
            'CollectorSim',
            step_size=60,
            save_h5=True,
            h5_storename='{}_data.h5'.format(filename),
            h5_framename='timeseries/sim_{}'.format(run_id),
            print_results=False)
        collector_entity = collector_sim.Collector()
        sim_dict['collector'] = collector_sim
        entity_dict['collector'] = collector_entity

        ##  PV
        pv_sim = world.start(
            'PVModel',
            eid_prefix='pv_',
            step_size=1)
        sim_dict['pv1'] = pv_sim

        pv_entity_1 = pv_sim.PVModel(
            rated_capacity=pv1_rated_capacity,
            series_name=weather)
        entity_dict['pv1'] = pv_entity_1

        ## WT
        #wt_sim = world.start(
        #    'WTModel',
        #    eid_prefix='wt_',
        #    step_size=1)
        #sim_dict['wt1'] = wt_sim
        #wt_entity_1 = wt_sim.WTModel(
        #    rated_capacity=1,
        #    series_name='/Gaia_20180510')
        #entity_dict['wt1'] = wt_entity_1

        ## Battery
        batt_sim = world.start(
            'BatteryModel',
            eid_prefix='batt_',
            step_size=1)
        batt_entity_1 = batt_sim.BatteryModel(
            rated_capacity=batt1_rated_cap,
            rated_discharge_capacity=batt1_rate,
            rated_charge_capacity=batt1_rate,
            initial_charge_rel=0.5,
            charge_change_rate=0.90,
        )
        sim_dict['batt'] = batt_sim
        entity_dict['batt1'] = batt_entity_1

        ## Controller
        control_sim = world.start(
            'Control',
            eid_prefix='demand_',
            step_size=1)
        control_entity_1 = control_sim.Control(setpoint_change_rate=con_change_rate)

        sim_dict['control'] = control_sim
        entity_dict['control1'] = control_entity_1

        ## "Noisifier"
        noise_injector = world.start(
            'NoiseGenerator',
            eid_prefix='ng_',
            step_size=1)
        noise_entity_1 = noise_injector.NoiseGenerator(scale=noise_scale)

        sim_dict['noise'] = noise_injector
        entity_dict['noise1'] = noise_entity_1

        return sim_dict, entity_dict

    world = mosaik.World(SIM_CONFIG)
    sim_dict, entity_dict = init_entities(world)

    # Connect units to grid busbar
    world.connect(entity_dict['demand1'], entity_dict['grid1'], ('P', 'P'))
    world.connect(entity_dict['pv1'], entity_dict['grid1'], ('P', 'P'))
    world.connect(entity_dict['batt1'], entity_dict['grid1'], ('P', 'P'))

    # Connect controller, note that the grid measurement is passed by a noise injection
    if stochastic:
        world.connect(entity_dict['grid1'], entity_dict['noise1'], ('Pgrid','input'))
        world.connect(entity_dict['noise1'], entity_dict['control1'], ('output', 'Pgrid'))
    else:
        world.connect(entity_dict['grid1'], entity_dict['control1'], ('Pgrid', 'Pgrid'))
    world.connect(entity_dict['control1'], entity_dict['batt1'], ('Pset', 'Pset'), time_shifted=True,
                  initial_data={'Pset': 0.0})
    world.connect(entity_dict['batt1'], entity_dict['control1'], ('relSoC', 'relSoC'))

    # Connect to Collector
    world.connect(entity_dict['batt1'], entity_dict['collector'], ('P', 'BattP'))
    world.connect(entity_dict['batt1'], entity_dict['collector'], ('SoC', 'BattSoC'))
    world.connect(entity_dict['demand1'], entity_dict['collector'], ('P', 'DemP'))
    if stochastic:
        world.connect(entity_dict['noise1'], entity_dict['collector'], ('output', 'GridP'))
    else:
        world.connect(entity_dict['grid1'], entity_dict['collector'], ('Pgrid', 'GridP'))
    world.connect(entity_dict['pv1'], entity_dict['collector'], ('P', 'SolarP'))
    #world.connect(entity_dict['wt1'], entity_dict['collector'], ('P', 'WindP'))

    END = 24*60*60-1 # 24 hours, 1 MosaikTime = 1 second
    world.run(END)
    ## End of simulation

    # # Data processing
    sim_store = pd.HDFStore('temp_files/{}_data.h5'.format(database_name))
    sim_data = sim_store['/timeseries/sim_{}'.format(run_data['ID'])]
    sim_store.close()
    # The measurement at PCC
    # print(sim_data.keys())
    if stochastic:
        grid_balance = sim_data['NoiseGenerator-0.ng__0']
    else:
        grid_balance = sim_data['SimpleGridModel-0.grid__0']
    cost = 0.39  # DKK/kWh
    energy_imported = grid_balance.apply(lambda x: x[x > 0].sum()) / 60
    energy_exported = grid_balance.apply(lambda x: x[x < 0].sum()) / 60
    energy_bill = energy_imported*cost
    maximum_infeed = np.abs(grid_balance).max()
    pv_generated = sim_data['PVModel-0.pv__0']['SolarP'].sum()/60

    self_consumption_index = (pv_generated+energy_exported)/pv_generated

    sim_data = {'ID': [run_id],
                'Energy bill [DKK/kWh]': [energy_bill['GridP']],
                'Max. in-feed [kW]': [maximum_infeed['GridP']],
                'Energy Imported [kWh]': [energy_imported['GridP']],
                'Energy exported [kWh]': [energy_exported['GridP']],
                'Self consumption index': [self_consumption_index['GridP']],
                'pv1_capacity [kW]': [pv1_cap],
                'battery storage capacity [kWh]': [battery_cap],
                'battery charge capacity[kW]': [battery_rate],
                'climate_conditions': [day_type],
                'controller_change_rate': [change_rate],
                'season': [season],
                'File ID/dataframe': ['{}'.format(database_name) + '/' + 'timeseries/sim_{}'.format(run_data['ID'])]}

    run_store = pd.HDFStore('temp_files/runs_summary_{}.h5'.format(database_name))
    run_df = pd.DataFrame(data=sim_data)
    print(run_df)
    run_store['run_{}'.format(run_data['ID'])] = run_df
    run_store.close()

ext: py

reason: import numpy


text:
# convert ensembl gff3 to gtf

import numpy as np
import pandas as pd
import csv
import argparse
import os

def main():
    parser = argparse.ArgumentParser(description='Convert a GFF3 file to basic GTF format')
    parser.add_argument('GFF', metavar='gff3_file', type=str, help='gff3 file to convert')
    parser.add_argument('-o', '--output', type=str, help='file name to write gtf')

    args = parser.parse_args()
    input_file = args.GFF

    # check if output file name specified / sub gff for gtf if no output filename
    if args.output is None:
        output_path = os.path.splitext(input_file)[0] + ".gtf"
    else:
        output_path = args.output

    # check if file already exists
    if os.path.isfile(output_path):
        print(output_path + " already exists")

        overwrite = input("Do you want to overwrite these files? ([Y]/n): ").lower().strip()[:1]
        if not (overwrite == "y" or overwrite == ""):
            sys.exit(1)
        else:
            # remove old files so you don't append new data to old files
            if os.path.isfile(output_path):
                os.remove(output_path)

    # read in gff
    gff3 = pd.read_csv(input_file, sep='\t', comment="#", header=None)
    gff3.columns = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes']

    # remove rows that are chromosomes, or CDS/UTR annotations. We don't need these ATM
    gff3 = gff3[~gff3.type.isin(['CDS', 'cds', 'chromosome', 'five_prime_utr', 'three_prime_utr', 'five_prime_UTR', 'three_prime_UTR'])]

    # what does the attributes col start with?
    starts_with = gff3.attributes.str.split(":").str[0]

    # replace any non-standard (i.e. noncoding transcripts) with gene/transcript
    gff3.loc[(starts_with == 'ID=gene'),'type'] = 'gene'
    gff3.loc[gff3.type.str.contains('gene'),'type'] = 'gene'

    gff3.loc[(starts_with == 'ID=transcript'),'type'] = 'transcript'
    gff3.loc[gff3.type.str.contains('RNA'),'type'] = 'transcript'
    gff3.loc[gff3.type.str.contains('exon'),'type'] = 'exon'

    # find gene_id
    gff3['gene_id'] = ''
    gff3.loc[gff3.type=='gene','gene_id'] = gff3.loc[gff3.type=='gene','attributes'].str.split("ID=").str[1].str.split(";").str[0]
    gff3.loc[gff3.type=='transcript','gene_id'] = gff3.loc[gff3.type=='transcript','attributes'].str.split("Parent=").str[1].str.split(";").str[0]

    # find transcript_id
    gff3['transcript_id'] = ''
    gff3.loc[gff3.type=='transcript','transcript_id'] = gff3.loc[gff3.type=='transcript','attributes'].str.split("ID=").str[1].str.split(";").str[0]
    gff3.loc[gff3.type=='exon','transcript_id'] = gff3.loc[gff3.type=='exon','attributes'].str.split("Parent=").str[1].str.split(";").str[0]

    # find exon_id
    gff3['exon_id'] = ''
    # check if ID= or Name= or exon_id=
    replace_with_exon_id = sum(gff3.loc[gff3.type=='exon','attributes'].str.contains('exon_id='))
    replace_with_Name = sum(gff3.loc[gff3.type=='exon','attributes'].str.contains('Name='))
    replace_with_ID = sum(gff3.loc[gff3.type=='exon','attributes'].str.contains('ID='))

    if replace_with_exon_id == len(gff3.loc[gff3.type=='exon','attributes']):
        gff3.loc[gff3.type=='exon','exon_id'] = gff3.loc[gff3.type=='exon','attributes'].str.split("exon_id=").str[1].str.split(";").str[0]
    elif replace_with_Name == len(gff3.loc[gff3.type=='exon','attributes']):
        gff3.loc[gff3.type=='exon','exon_id'] = gff3.loc[gff3.type=='exon','attributes'].str.split("Name=").str[1].str.split(";").str[0]
    elif replace_with_ID == len(gff3.loc[gff3.type=='exon','attributes']):
        gff3.loc[gff3.type=='exon','exon_id'] = gff3.loc[gff3.type=='exon','attributes'].str.split("ID=").str[1].str.split(";").str[0]
    else:
        gff3.loc[gff3.type=='exon','exon_id'] = gff3.loc[gff3.type=='exon','attributes']

    # to find gene_id for exons, copy gene_id from parent transcript
    transcript_ids = gff3.copy()
    transcript_ids = transcript_ids.loc[gff3.type=='transcript',['gene_id','transcript_id']]
    gff3 = pd.merge(gff3, transcript_ids, left_on='transcript_id', right_on='transcript_id', how='left')
    gff3.loc[gff3.type=='exon','gene_id_x'] = gff3.loc[gff3.type=='exon','gene_id_y']

    # add a gtf formatted description/attributes column
    gff3['gtf_desc'] = ''
    gff3.loc[gff3.type=='gene','gtf_desc'] = 'gene_id "' + gff3.loc[gff3.type=='gene','gene_id_x'].astype(str) + '";'
    gff3.loc[gff3.type=='transcript','gtf_desc'] = 'gene_id "' + gff3.loc[gff3.type=='transcript','gene_id_x'].astype(str) + '";' + ' transcript_id "' + gff3.loc[gff3.type=='transcript','transcript_id'].astype(str) + '";'
    gff3.loc[gff3.type=='exon','gtf_desc'] = 'gene_id "' + gff3.loc[gff3.type=='exon','gene_id_x'].astype(str) + '";' + ' transcript_id "' + gff3.loc[gff3.type=='exon','transcript_id'].astype(str) + '";' + ' exon_id "' + gff3.loc[gff3.type=='exon','exon_id'].astype(str) + '";'

    # write file
    gff3.to_csv(output_path, sep='\t', columns = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase','gtf_desc'], header=False, index=False, quoting=csv.QUOTE_NONE)


text:
from jax.config import config
config.update("jax_enable_x64", True)
from jaxns.likelihood_samplers.ellipsoid_utils import minimum_volume_enclosing_ellipsoid, sample_ellipsoid
from jax import numpy as jnp, random
from jax.lax import while_loop, dynamic_update_slice, scan

def debug_triangular_solve():
    from jax.lax_linalg import triangular_solve
    from scipy.linalg import solve_triangular
    from jax.scipy.linalg import solve_triangular as solve_triangular_jax
    import jax.numpy as jnp
    ndims = 2
    A = jnp.diag(jnp.ones(ndims))
    A = jnp.where(A == 0., 0.95, A)
    b = jnp.ones(ndims)
    L = jnp.linalg.cholesky(A)
    assert jnp.all(L @ L.T == A)

    x = jnp.linalg.solve(L, b)
    print("Solving L.x  = b with scipy")
    print("x should be {}".format(x))
    scipy_x = solve_triangular(L, b, lower=True)
    assert jnp.all(scipy_x == x)
    print("Works as expected!")

    print("Now note JAX's solution to L^T.x = b corresponds to scipy's L.x = b")
    jax_x = triangular_solve(L, b, lower=True, transpose_a=True)
    assert jnp.all(jax_x==scipy_x)
    print("Likewise, JAX's solution to L.x = b corresponds to scipy's L^T.x = b")
    assert jnp.all(triangular_solve(L, b, lower=True) == solve_triangular(L, b, lower=True, trans=1))

    print("Note, I have not tested for the L^H.x=b case.")

    jax_x = solve_triangular_jax(L, b, lower=True)
    assert jnp.all(scipy_x == jax_x)

def debug_mvee():
    import pylab as plt

    n = random.normal(random.PRNGKey(0), (10000,2))
    n = n /jnp.linalg.norm(n, axis=1, keepdims=True)
    angle = jnp.arctan2(n[:,1], n[:,0])
    plt.hist(angle, bins=100)
    plt.show()
    N = 120
    D = 2
    points = random.uniform(random.PRNGKey(0), (N, D))

    from jax import disable_jit
    with disable_jit():
        center, radii, rotation = minimum_volume_enclosing_ellipsoid(points, 0.01)

    plt.hist(jnp.linalg.norm((rotation.T @ (points.T - center[:, None])) / radii[:, None], axis=0))
    plt.show()
    print(center, radii, rotation)
    plt.scatter(points[:, 0], points[:, 1])
    theta = jnp.linspace(0., jnp.pi*2, 100)
    ellipsis = center[:, None] + rotation @ jnp.stack([radii[0]*jnp.cos(theta), radii[1]*jnp.sin(theta)], axis=0)
    plt.plot(ellipsis[0,:], ellipsis[1,:])

    for i in range(1000):
        y = sample_ellipsoid(random.PRNGKey(i), center, radii, rotation)
        plt.scatter(y[0], y[1])

    C = jnp.linalg.pinv(jnp.cov(points, rowvar=False, bias=True))
    p = (N - D - 1)/N
    def q(p):
        return p + p**2/(4.*(D-1))
    C = C / q(p)
    c = jnp.mean(points, axis=0)
    W, Q, Vh = jnp.linalg.svd(C)
    radii = jnp.reciprocal(jnp.sqrt(Q))
    rotation = Vh.conj().T
    ellipsis = c[:, None] + rotation @ jnp.stack([radii[0] * jnp.cos(theta), radii[1] * jnp.sin(theta)], axis=0)
    plt.plot(ellipsis[0, :], ellipsis[1, :])

    plt.show()

def debug_riddle():
    import numpy as np
    import pylab as plt

    t = np.linspace(0., 12., 10000)
    m = (t - np.floor(t))*2*np.pi
    h = t/12.*2.*np.pi

    plt.plot(t, h)
    plt.show()
    plt.plot(t, m)
    plt.show()

    def time(m, h):
        return h/(2.*np.pi)*12.

    plt.plot(t, time(m,h))
    plt.plot(t, time(h,m))
    plt.show()

def debug_vmap_bfgs():
    import jax.numpy as jnp
    from jax import jit, config
    from jax.scipy.optimize import minimize
    import os
    config.enable_omnistaging()
    ncpu=2
    os.environ['XLA_FLAGS'] = f"--xla_force_host_platform_device_count={ncpu}"

    def cost_fn(x):
        return -jnp.sum(x**2)

    x = random.uniform(random.PRNGKey(0), (3,), minval=-1, maxval=1)
    result = jit(lambda x: minimize(cost_fn, x, method='BFGS'))(x)
    print(result)

def debug_random_angles():
    from jax import random
    import jax.numpy as jnp
    import pylab as plt
    S = 10000
    p1 = random.uniform(random.PRNGKey(0),shape=(S,2))
    p1 /= jnp.linalg.norm(p1, axis=-1, keepdims=True)
    plt.hist(jnp.arctan2(p1[:,1],p1[:,0]), bins='auto', alpha=0.5)
    p2 = random.normal(random.PRNGKey(0), shape=(S, 2))
    p2 /= jnp.linalg.norm(p2, axis=-1, keepdims=True)
    plt.hist(jnp.arctan2(p2[:,1],p2[:,0]), bins='auto', alpha=0.5)
    plt.show()

def debug_simplex_sampler():
    from jax import random
    import jax.numpy as jnp
    import pylab as plt
    from jaxns.gaussian_process.utils import squared_norm
    S = 10
    D = 2
    u = random.uniform(random.PRNGKey(0),shape=(S,2))
    plt.scatter(u[:,0],u[:,1])
    inter_point_distance = squared_norm(u, u)
    inter_point_distance = jnp.where(inter_point_distance < 1e-20, jnp.inf, inter_point_distance)
    knn_indices = jnp.argsort(inter_point_distance, axis=-1)[:, :D + 1]
    i = 2
    plt.scatter(u[i,0],u[i,1],c='red')
    for j in range(S):
        for k in range(D+1):
            plt.plot([u[i,0], u[knn_indices[i,k],0]],[u[i,1], u[knn_indices[i,k],1]])
    plt.show()

def debug_fft():
    from jax.config import config

    config.update("jax_enable_x64", True)

    import time
    import numpy as np

    import jax
    from jax import numpy as jnp

    np.random.seed(0)

    signal = np.random.randn(2 ** 20)
    signal_jax = jnp.array(signal)

    jfft = jax.jit(jnp.fft.fft)

    import tensorflow as tf
    signal_tf = tf.constant(signal, dtype=tf.complex128)

    def tffft(x):
        return tf.signal.fft(x).numpy()

    X_np = np.fft.fft(signal)
    X_jax = jfft(signal_jax)
    X_tf = tffft(signal_tf)

    print(np.mean(np.abs(X_np)))
    print("With JAX:")
    print('max:\t', jnp.max(jnp.abs(X_np - X_jax)))
    print('mean:\t', jnp.mean(jnp.abs(X_np - X_jax)))
    print('min:\t', jnp.min(jnp.abs(X_np - X_jax)))

    print("With Tensorflow:")
    print('max:\t', jnp.max(jnp.abs(X_np - X_tf)))
    print('mean:\t', jnp.mean(jnp.abs(X_np - X_tf)))
    print('min:\t', jnp.min(jnp.abs(X_np - X_tf)))

    ### CPU
    # 907.3490574884647
    # max:	 2.8773885332210747
    # mean:	 0.3903197564919141
    # min:	 2.4697454729898156e-05

    ### GPU
    # 907.3490574884647
    # max:	 0.001166179716824765
    # mean:	 0.00020841654559267488
    # min:	 2.741492442122853e-07

    R = 100
    ts = time.time()
    for i in range(R):
        _ = np.fft.fft(signal)
    print('numpy fft execution time [ms]:\t', (time.time() - ts) / R * 1000)

    # Compile
    _ = jfft(signal_jax).block_until_ready()

    ts = time.time()
    for i in range(R):
        _ = jfft(signal_jax).block_until_ready()
    print('jax fft execution time [ms]:\t', (time.time() - ts) / R * 1000)

    ts = time.time()
    for i in range(R):
        _ = tffft(signal_tf)
    print('tensorflow fft execution time [ms]:\t', (time.time() - ts) / R * 1000)

def delay_accuracy():
    import astropy.units as au
    from jax import vmap
    V = 30.*au.km/au.s
    c = 299792458.*au.m/au.s
    def delay_perfect(b):
        return -b/c * (1. - 0.5*V**2/c**2) - b*V/c**2 * (1. + 0.5*V/c)/(1. + V/c)
    def delay_nonperfect(b):
        return -b/c
    def delay(b,s, V):
        return (-(b@s)/c * (1. - 0.5*V@V/c**2) - (b @ V)/c**2 * (1. + 0.5*(s@V)/c))/(1.+(s@V)/c)

    for _b in 10**jnp.linspace(0., 4, 10):
        delay_diff = []
        for i in range(1000):
            keys = random.split(random.PRNGKey(i), 3)
            b = random.normal(keys[0], shape=(3,))
            b /= jnp.linalg.norm(b)
            b = b * _b * au.km
            s = random.normal(keys[1], shape=(3,))
            s /= jnp.linalg.norm(s)
            V = random.normal(keys[2], shape=(3,))
            V /= jnp.linalg.norm(V)
            V = V * 30.*au.km/au.s
            delay_diff.append(jnp.abs(delay(b, s, V).to(au.ns).value - delay(b, s, 0*V).to(au.ns).value))
        delay_diff = jnp.array(delay_diff)
        print("Baseline={:.2f}, max diff={:.2f}ns, mean diff={:.2f}ns".format(
            _b*au.km, jnp.max(delay_diff), jnp.mean(delay_diff)
        ))
        # print('For baseline of ',b,'magnitude of relativistic corrections', -delay_perfect(b*au.km).to(au.ns)+delay_nonperfect(b*au.km).to(au.ns))

def quick_select(k, a):
    """
    Performs iterative quick select.
    Args:
        k: the kth smallest element to select
        a: the 1D array to select from

    Returns: the kth smallest element in `a`
    """

    def partition(arr, low, high):
        """
        Lomuto partition function.
        """
        if len(arr.shape) > 1:
            raise ValueError("Partition works on 1D arrays. Use vmap.")
        pivot = arr[high]
        def body(state):
            (j, i,arr) = state
            do_swap = arr[j] <= pivot
            i = jnp.where(do_swap, i+1, i)
            ai = arr[i, None]
            aj = arr[j, None]
            arr_swapped = dynamic_update_slice(arr, aj, [i])
            arr_swapped = dynamic_update_slice(arr_swapped, ai, [j])
            arr_swapped = jnp.where(do_swap, arr_swapped, arr)
            return (j + 1, i, arr_swapped)

        (j, i, arr) = while_loop(lambda state: state[0] < high,
                              body,(low, low - 1, arr))

        ai = arr[i+1, None]
        aj = arr[high, None]
        arr_swapped = dynamic_update_slice(arr, aj, [i+1])
        arr_swapped = dynamic_update_slice(arr_swapped, ai, [high])
        return (i + 1, arr_swapped)

    def body(state):
        (left, right, _, _, arr) = state
        pivot_index, arr = partition(arr, left, right)
        done = pivot_index == k-1
        kth_smallest = arr[pivot_index]
        right = jnp.where(pivot_index > k-1, pivot_index - 1, right)
        left = jnp.where(pivot_index < k-1, pivot_index + 1, left)
        return (left, right, done, kth_smallest, arr)

    (left, right, _, kth_smallest, arr) = while_loop(lambda state: (state[0]<=state[1]) & (~state[2]),
                               body,
                               (jnp.asarray(0), jnp.asarray(len(a) - 1), jnp.array(False), a[0], a))
    return kth_smallest

def quick_sort_median(a):
    """
    Performs a quickselect median which is equivalent to,
        jnp.percentile(a, 50, interpolation='higher')
    Thus, for odd sized arrays it's the same as jnp.median.
    Args:
        a: 1D array. To do ND arrays use vmap.

    Returns: median of array choosing the 'higher' of the two middle values
    when `a` has an even number of points.
    """
    if len(a.shape) > 1:
        raise ValueError("Only works on 1D arrays. Use vmap.")
    k = a.size//2 + 1
    median = quick_select(k, a)
    return median

def test_quick_select():

    arr = jnp.asarray([10, 4, 5, 8, 11, 6, 26, 7]) # even number of points
    assert quick_sort_median(arr) == jnp.percentile(arr, 50, interpolation='higher')
    arr = jnp.asarray([10, 4, 5, 8, 11, 6, 26]) # odd number of points
    assert quick_sort_median(arr) == jnp.percentile(arr, 50, interpolation='higher')
    arr = jnp.asarray([10, 4, 5, 8, 11, 6, 26])  # odd number of points
    assert quick_sort_median(arr) == jnp.median(arr)

    arr = jnp.asarray([10, 4, 5, 8, 11, 6, 26, 7])  # even number of points
    try:
        assert quick_sort_median(arr) == jnp.median(arr)
    except AssertionError:
        print("Not equivalent to median when array has an even size.")

def test_median_speed():
    from jax import jit, random
    import numpy as np
    from timeit import default_timer
    import pylab as plt

    jax_median = jit(jnp.median)
    quick_median = jit(quick_sort_median)

    n = []
    jax_time = []
    quick_time = []
    numpy_time = []
    for b in range(1,15):
        random_array = random.uniform(random.PRNGKey(b), shape=(2**b + 1,))
        jax_median(random_array)
        quick_median(random_array)
        for i in range(10):
            n.append(2**b + 1)
            # random_array = random.uniform(random.PRNGKey(b+1), shape=(2 ** b + 1,))
            t0 = default_timer()
            print(jax_median(random_array))
            jax_time.append(default_timer() - t0)
            t0 = default_timer()
            print(quick_median(random_array))
            quick_time.append(default_timer() - t0)
            random_array = np.array(random_array)
            t0 = default_timer()
            print(np.median(random_array))
            numpy_time.append(default_timer() - t0)
    plt.scatter(n, jax_time, label='jax median')
    plt.scatter(n, quick_time, label='quick select median')
    plt.scatter(n, numpy_time, label='numpy median')
    plt.legend()
    plt.yscale('log')
    plt.ylabel('time (s)')
    plt.show()

def test_speed_jax():

    import numpy as np
    from jax import jit,value_and_grad, random
    from jaxns.optimize import minimize as minimize_jax
    from scipy.optimize import minimize as minimize_np
    import pylab as plt
    from timeit import default_timer

    S = 3
    t_scipy_halfjax,t_scipy_jax,t_jax,t_numpy = [],[],[],[]
    N_array = [2,10,50, 100, 200, 400]
    for N in N_array:
        print("Working on N={}".format(N))

        A = random.normal(random.PRNGKey(0), shape=(N, N))
        u = jnp.ones(N)
        x0 = -2. * jnp.ones(N)

        def f_prescale(x, u):
            y = A @ x
            dx = u - y
            return jnp.sum(dx**2) + 0.1*jnp.sum(jnp.abs(x))

        # Due to https://github.com/google/jax/issues/4594 we scale the loss
        # so that scipy and jax linesearch perform similarly.
        jac_norm = jnp.linalg.norm(value_and_grad(f_prescale)(x0, u)[1])
        jac_norm_np = np.array(jac_norm)

        def f(x, u):
            y = A @ x
            dx = u - y
            return (jnp.sum(dx**2) + 0.1*jnp.sum(jnp.abs(x)))/jac_norm

        def f_np(x, u):
            y = A @ x[:,None]
            dx = u - y
            return (np.sum(dx**2) + 0.1*np.sum(np.abs(x)))/jac_norm_np

        print("Testing scipy+numpy")
        t0 = default_timer()
        args= (np.array(x0), (np.array(u),))
        results_np = minimize_np(f_np, *args, method='BFGS')
        for _ in range(S):
            results_np = minimize_np(f_np, *args, method='BFGS')
        t_numpy.append((default_timer() - t0) / S)
        print("nfev",results_np.nfev, "njev", results_np.njev)
        print("Time for scipy + numpy", t_numpy[-1])

        print("Testing scipy + jitted function and numeric grad")
        @jit
        def _f(x0, u):
            return f(x0, u)

        _f(x0, u).block_until_ready()
        t0 = default_timer()
        for _ in range(S):
            results_np = minimize_np(_f, x0, (u,), method='BFGS')
        t_scipy_halfjax.append((default_timer() - t0) / S)
        print("nfev",results_np.nfev, "njev", results_np.njev)
        print("Time for scipy + jitted function and numeric grad", t_scipy_halfjax[-1])

        print("Testing scipy + jitted function and grad")
        @jit
        def _f(x0, u):
            v, g = value_and_grad(f)(x0, u)
            return v, g
        _f(x0, u)[1].block_until_ready()
        t0 = default_timer()
        for _ in range(S):
            results_np = minimize_np(_f, x0, (u,), method='BFGS', jac=True)
        t_scipy_jax.append((default_timer() - t0) / S)
        print("nfev",results_np.nfev, "njev", results_np.njev)
        print("Time for scipy + jitted function and grad", t_scipy_jax[-1])

        print("Testing pure JAX implementation")
        @jit
        def do_minimize_jax(x0, u):
            results = minimize_jax(f, x0, args=(u,),method='BFGS')
            return results.x

        results_jax = minimize_jax(f, x0, args=(u,),method='BFGS')
        do_minimize_jax(x0, u).block_until_ready()

        t0 = default_timer()
        for _ in range(S):
            do_minimize_jax(x0, u).block_until_ready()
        t_jax.append((default_timer() - t0)/S)
        print("nfev", results_jax.nfev, "njev", results_jax.njev)
        print("Time for pure JAX implementation", t_jax[-1])

    plt.figure(figsize=(8,5))
    plt.plot(N_array,t_scipy_jax,label='scipy+jitted(func and grad)')
    plt.plot(N_array,t_scipy_halfjax,label='scipy+jitted(func)')
    plt.plot(N_array,t_jax,label='pure JAX')
    plt.plot(N_array,t_numpy,label='scipy+numpy')
    plt.yscale('log')
    plt.legend()
    plt.title("Run time of BFGS on N-D Least squares + L1 regularisation.")
    plt.ylabel('Time [s]')
    plt.xlabel("N")
    plt.show()

def test_speed_MM_jax():
    from jax import jit,vmap
    import numpy as np
    import pylab as plt
    from timeit import default_timer
    S = 10
    plt.figure(figsize=(8, 5))
    t1,t2 = [],[]
    for N in [100,200,400,800,1600,3200]:

        M1 = jnp.ones((N,N))
        M2 = jnp.ones((N,N))
        M3 = jnp.ones((N,N))

        def mm(M1,M2,M3):
            return (M1**2 @ M2**2 + M3**2).sum()**2

        _M1 = np.array(M1)
        _M2 = np.array(M2)
        _M3 = np.array(M3)
        t0 = default_timer()
        for _ in range(S):
            results_jax = mm(_M1, _M2, _M3)
        t1.append((default_timer() - t0) / S)
        print("Time for numpy", t1)

        _mm = jit(mm)
        _mm(M1,M2,M3).block_until_ready()
        t0 = default_timer()
        for _ in range(S):
            results_jax = _mm(M1,M2,M3).block_until_ready()
        t2.append((default_timer() - t0) / S)
        print("Time for jax", t2)

    plt.plot([100,200,400,800,1600,3200], t1, label='numpy')
    plt.plot([100,200,400,800,1600,3200], t2, label='jax')
    plt.legend()
    plt.yscale('log')
    plt.xlabel("N")
    plt.title("GEMM: M1 @ M2 + M3")
    plt.ylabel('Time [s]')
    plt.show()

def test_example_omni():
    import jax.numpy as jnp
    from jax import jit

    @jit
    def select_tril(x):
        mask = jnp.arange(x.shape[0])[:, None] > jnp.arange(x.shape[1])
        return jnp.where(mask, x, jnp.zeros_like(x))

    import numpy as np
    x = np.arange(12).reshape((3, 4)) # mix in numpy as desired
    select_tril(x)

if __name__ == '__main__':
    # test_example_omni()
    test_speed_jax()
    # delay_accuracy()
    # debug_vmap_bfgs()
    # debug_random_angles()
    # debug_simplex_sampler()
    # debug_triangular_solve()
    # debug_fft()
    # debug_nestest_sampler()
    # debug_tec_clock_inference()
    # debug_riddle()
    # debug_mvee()

max_stars_repo_name: fehiepsi/jaxns

max_issues_repo_name: fehiepsi/jaxns

max_forks_repo_name: fehiepsi/jaxns

reason: import numpy,from scipy,import astropy,import jax,from jax


text:
# encoding: utf-8
"""
@author: jemmy li
@contact: zengarden2009@gmail.com
"""

from IPython import embed
from .light_config import cfg
import tensorflow as tf
import tensorflow.contrib.slim as slim
# from tensorflow.contrib.slim.python.slim.nets import resnet_utils, resnet_v1
import numpy as np

from utils.tf_utils.basemodel import resnet_utils, resnet_v1
from tensorflow.contrib.slim import arg_scope
from tensorflow.python.framework import ops
from tensorflow.python.ops import nn_ops
from tensorflow.contrib.layers.python.layers import regularizers, \
    initializers, layers

# todo : delete snippets
from detection_opr.rpn.snippets import generate_anchors_opr
from detection_opr.rpn_batched.proposal_target_layer import proposal_target_layer
from detection_opr.rpn_batched.anchor_target_layer_without_boxweight import \
    anchor_target_layer
from detection_opr.rpn_batched.proposal_opr import proposal_opr

from detection_opr.utils import loss_opr
from lib_kernel.lib_psroi_pooling import psroi_pooling_op, psroi_pooling_op_grad
from lib_kernel.lib_psalign_pooling import psalign_pooling_op, psalign_pooling_op_grad
from detection_opr.rfcn_plus_plus import rfcn_plus_plus_opr
from collections import OrderedDict as dict

def resnet_arg_scope(
        is_training=True, weight_decay=cfg.weight_decay, batch_norm_decay=0.997,
        batch_norm_epsilon=1e-5, batch_norm_scale=True):
    batch_norm_params = {
        'is_training': False, 'decay': batch_norm_decay,
        'epsilon': batch_norm_epsilon, 'scale': batch_norm_scale,
        'trainable': cfg.bn_training,
        'updates_collections': ops.GraphKeys.UPDATE_OPS
    }

    with arg_scope(
            [slim.conv2d],
            weights_regularizer=regularizers.l2_regularizer(weight_decay),
            weights_initializer=initializers.variance_scaling_initializer(),
            trainable=is_training,
            activation_fn=nn_ops.relu,
            normalizer_fn=layers.batch_norm,
            normalizer_params=batch_norm_params):
        with arg_scope([layers.batch_norm], **batch_norm_params) as arg_sc:
            return arg_sc

class Network(object):
    def __init__(self):
        pass

    def get_inputs(self, mode=0):
        """gets network inputs
        Returns:
             mode:
                0: return inputs tensor_holder
                1: return inputs name
        """
        if mode == 0:
            inputs = []
            inputs.append(tf.placeholder(tf.float32, shape=[None, None, None, 3]))
            inputs.append(tf.placeholder(tf.float32, shape=[None, 6]))
            inputs.append(tf.placeholder(tf.float32, shape=[None, None, 5]))
            return inputs
        elif mode == 1:
            inputs_names = ['data', 'im_info', 'boxes']
            return inputs_names

    def inference(self, mode, inputs):
        is_training = mode == 'TRAIN'

        ###decode your inputs
        [image, im_info, gt_boxes] = inputs

        image.set_shape([None, None, None, 3])
        im_info.set_shape([None, cfg.nr_info_dim])
        if mode == 'TRAIN':
            gt_boxes.set_shape([None, None, 5])
        ##end of decode

        num_anchors = len(cfg.anchor_scales) * len(cfg.anchor_ratios)
        bottleneck = resnet_v1.bottleneck

        blocks = [
            resnet_utils.Block('block1', bottleneck,
                               [(256, 64, 1, 1)] * 2 + [(256, 64, 1, 1)]),
            resnet_utils.Block('block2', bottleneck,
                               [(512, 128, 2, 1)] + [(512, 128, 1, 1)] * 3),
            resnet_utils.Block('block3', bottleneck,
                               [(1024, 256, 2, 1)] + [(1024, 256, 1, 1)] * 22),
            resnet_utils.Block('block4', bottleneck,
                               [(2048, 512, 1, 2)] + [(2048, 512, 1, 2)] * 2)
        ]

        with slim.arg_scope(resnet_arg_scope(is_training=False)):
            with tf.variable_scope('light_resnet_v1_101', 'light_resnet_v1_101'):
                net = resnet_utils.conv2d_same(
                    image, 64, 7, stride=2, scope='conv1')
                net = slim.max_pool2d(
                    net, [3, 3], stride=2, padding='SAME', scope='pool1')
            net, _ = resnet_v1.resnet_v1(
                net, blocks[0:1], global_pool=False, include_root_block=False,
                scope='light_resnet_v1_101')

        with slim.arg_scope(resnet_arg_scope(is_training=is_training)):
            net_conv3, _ = resnet_v1.resnet_v1(
                net, blocks[1:2], global_pool=False, include_root_block=False,
                scope='light_resnet_v1_101')

        with slim.arg_scope(resnet_arg_scope(is_training=is_training)):
            net_conv4, _ = resnet_v1.resnet_v1(
                net_conv3, blocks[2:3], global_pool=False,
                include_root_block=False, scope='light_resnet_v1_101')

        with slim.arg_scope(resnet_arg_scope(is_training=is_training)):
            net_conv5, _ = resnet_v1.resnet_v1(
                net_conv4, blocks[-1:], global_pool=False,
                include_root_block=False, scope='light_resnet_v1_101')

        initializer = tf.random_normal_initializer(mean=0.0, stddev=0.01)
        initializer_bbox = tf.random_normal_initializer(mean=0.0, stddev=0.001)

        with tf.variable_scope(
                'light_resnet_v1_101', 'light_resnet_v1_101',
                regularizer=tf.contrib.layers.l2_regularizer(
                    cfg.weight_decay)):
            # rpn
            rpn = slim.conv2d(
                net_conv4, 512, [3, 3], trainable=is_training,
                weights_initializer=initializer, activation_fn=nn_ops.relu,
                scope="rpn_conv/3x3")
            rpn_cls_score = slim.conv2d(
                rpn, num_anchors * 2, [1, 1], trainable=is_training,
                weights_initializer=initializer, padding='VALID',
                activation_fn=None, scope='rpn_cls_score')
            rpn_bbox_pred = slim.conv2d(
                rpn, num_anchors * 4, [1, 1], trainable=is_training,
                weights_initializer=initializer, padding='VALID',
                activation_fn=None, scope='rpn_bbox_pred')

            # generate anchor
            height = tf.cast(tf.shape(rpn)[1], tf.float32)
            width = tf.cast(tf.shape(rpn)[2], tf.float32)
            anchors = generate_anchors_opr(
                height, width, cfg.stride[0], cfg.anchor_scales,
                cfg.anchor_ratios)
            # change it so that the score has 2 as its channel size
            rpn_cls_prob = tf.reshape(rpn_cls_score, [-1, 2])
            rpn_cls_prob = tf.nn.softmax(rpn_cls_prob, name='rpn_cls_prob')
            rpn_cls_prob = tf.reshape(rpn_cls_prob, tf.shape(rpn_cls_score))

            rois, roi_scores = proposal_opr(
                rpn_cls_prob, rpn_bbox_pred, im_info, mode, cfg.stride,
                anchors, num_anchors, is_tfchannel=True, is_tfnms=True)

            if is_training:
                with tf.variable_scope('anchor') as scope:
                    rpn_labels, rpn_bbox_targets = \
                        tf.py_func(
                            anchor_target_layer,
                            [gt_boxes, im_info, cfg.stride, anchors,
                             num_anchors],
                            [tf.float32, tf.float32])
                    rpn_labels = tf.to_int32(rpn_labels, name="to_int32")

                with tf.control_dependencies([rpn_labels]):
                    with tf.variable_scope('rpn_rois') as scope:
                        rois, labels, bbox_targets = \
                            tf.py_func(
                                proposal_target_layer,
                                [rois, gt_boxes, im_info],
                                [tf.float32, tf.float32, tf.float32])
                        labels = tf.to_int32(labels, name="to_int32")

        with tf.variable_scope(
                'light_resnet_v1_101', 'light_resnet_v1_101',
                regularizer=tf.contrib.layers.l2_regularizer(
                    cfg.weight_decay)):

            ps_chl = 7 * 7 * 10
            ps_fm = rfcn_plus_plus_opr.global_context_module(
                net_conv5, prefix='conv_new_1',
                ks=15, chl_mid=256, chl_out=ps_chl)
            ps_fm = nn_ops.relu(ps_fm)

            [psroipooled_rois, _, _] =  psalign_pooling_op.psalign_pool(
                ps_fm, rois, group_size=7,
                sample_height=2, sample_width=2, spatial_scale=1.0/16.0)

            #[psroipooled_rois, _] = psroi_pooling_op.psroi_pool(
            #    ps_fm, rois, group_size=7, spatial_scale=1.0 / 16.0)
            psroipooled_rois = slim.flatten(psroipooled_rois)
            ps_fc_1 = slim.fully_connected(
                psroipooled_rois, 2048, weights_initializer=initializer,
                activation_fn=nn_ops.relu, trainable=is_training, scope='ps_fc_1')
            cls_score = slim.fully_connected(
                ps_fc_1, cfg.num_classes, weights_initializer=initializer,
                activation_fn=None, trainable=is_training, scope='cls_fc')
            bbox_pred = slim.fully_connected(
                ps_fc_1, 4 * cfg.num_classes, weights_initializer=initializer_bbox,
                activation_fn=None, trainable=is_training, scope='bbox_fc')

            cls_prob = loss_opr.softmax_layer(cls_score, "cls_prob")

            #conv_new_1 = slim.conv2d(
            #    net_conv5, 1024, [1, 1], trainable=is_training,
            #    weights_initializer=initializer, activation_fn=nn_ops.relu,
            #    scope="conv_new_1")
            #rfcn_cls = slim.conv2d(
            #    conv_new_1, 7 * 7 * cfg.num_classes, [1, 1],
            #    trainable=is_training, weights_initializer=initializer,
            #    activation_fn=None, scope="rfcn_cls")
            #rfcn_bbox = slim.conv2d(
            #    conv_new_1, 7 * 7 * 4, [1, 1], trainable=is_training,
            #    weights_initializer=initializer,
            #    activation_fn=None, scope="rfcn_bbox")

            #[psroipooled_cls_rois, _] = psroi_pooling_op.psroi_pool(
            #    rfcn_cls, rois, group_size=7, spatial_scale=1.0 / 16.0)
            #[psroipooled_loc_rois, _] = psroi_pooling_op.psroi_pool(
            #    rfcn_bbox, rois, group_size=7, spatial_scale=1.0 / 16.0)

            #cls_score = tf.reduce_mean(psroipooled_cls_rois, axis=[1, 2])
            #bbox_pred = tf.reduce_mean(psroipooled_loc_rois, axis=[1, 2])
            #cls_prob = loss_opr.softmax_layer(cls_score, "cls_prob")
            # cls_prob = tf.nn.softmax(cls_score, name="cls_prob")
            #bbox_pred = tf.tile(bbox_pred, [1, cfg.num_classes])

        if not is_training:
            stds = np.tile(
                np.array(cfg.TRAIN.BBOX_NORMALIZE_STDS), (cfg.num_classes))
            means = np.tile(
                np.array(cfg.TRAIN.BBOX_NORMALIZE_MEANS), (cfg.num_classes))
            bbox_pred *= stds
            bbox_pred += means

            ##############add prediction#####################
            tf.add_to_collection("rpn_cls_score", rpn_cls_score)
            tf.add_to_collection("rpn_cls_prob", rpn_cls_prob)
            tf.add_to_collection("rpn_bbox_pred", rpn_bbox_pred)
            tf.add_to_collection("cls_score", cls_score)
            tf.add_to_collection("cls_prob", cls_prob)
            tf.add_to_collection("bbox_pred", bbox_pred)
            tf.add_to_collection("rois", rois)

        else:
            #--------------------  rpn loss ---------------------------------#
            from detection_opr.utils import loss_opr_without_box_weight
            rpn_loss_box = loss_opr_without_box_weight.smooth_l1_loss_rpn(
                tf.reshape(rpn_bbox_pred, [-1, 4]),
                tf.reshape(rpn_bbox_targets, [-1, 4]),
                tf.reshape(rpn_labels, [-1]), sigma=cfg.simga_rpn)

            rpn_cls_score = tf.reshape(rpn_cls_score, [-1, 2])
            rpn_label = tf.reshape(rpn_labels, [-1])
            rpn_select = tf.where(tf.not_equal(rpn_label, -1))
            rpn_cls_score = tf.reshape(
                tf.gather(rpn_cls_score, rpn_select), [-1, 2])
            rpn_label = tf.reshape(tf.gather(rpn_label, rpn_select), [-1])
            rpn_cross_entropy = tf.reduce_mean(
                tf.nn.sparse_softmax_cross_entropy_with_logits(
                    logits=rpn_cls_score, labels=rpn_label))

            #-------------------- rcnn loss  --------------------------------#
            label = tf.reshape(labels, [-1])
            cross_entropy, loss_box = loss_opr_without_box_weight.sum_ohem_loss(
                tf.reshape(cls_score, [-1, cfg.num_classes]), label,
                bbox_pred, bbox_targets, cfg.TRAIN.nr_ohem_sampling,
                cfg.num_classes)
            loss_box *= 2

            #--------------------add to colloection ------------------------#
            tf.add_to_collection('loss_cls', cross_entropy)
            tf.add_to_collection('loss_box', loss_box)
            tf.add_to_collection('rpn_loss_cls', rpn_cross_entropy)
            tf.add_to_collection('rpn_loss_box', rpn_loss_box)
            loss = cross_entropy + loss_box + rpn_cross_entropy + rpn_loss_box
            tf.add_to_collection('losses', loss)
            return loss

    def get_train_collection(self):
        ret = dict()
        ret['rpn_loss_cls'] = tf.add_n(tf.get_collection('rpn_loss_cls'))
        ret['rpn_loss_box'] = tf.add_n(tf.get_collection('rpn_loss_box'))
        ret['loss_cls'] = tf.add_n(tf.get_collection('loss_cls'))
        ret['loss_box'] = tf.add_n(tf.get_collection('loss_box'))
        ret['tot_losses'] = tf.add_n(tf.get_collection('losses'))
        return ret

    def get_test_collection(self):
        ret = dict()
        ret['cls_score'] = tf.get_collection('cls_score')[0]
        ret['cls_prob'] = tf.get_collection('cls_prob')[0]
        ret['bbox_pred'] = tf.get_collection('bbox_pred')[0]
        ret['rois'] = tf.get_collection('rois')[0]
        return ret


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
This module deals with hexagonal twinboundary structure.
"""

from copy import deepcopy
import math
import numpy as np
from twinpy.structure.base import _BaseTwinStructure
from twinpy.structure.shear import ShearStructure

class TwinBoundaryStructure(_BaseTwinStructure):
    """
    Twinboundary structure class.
    """

    def __init__(
           self,
           lattice:np.array,
           symbol:str,
           twinmode:str,
           twintype:int,
           wyckoff:str='c',
           ):
        """
        Args:
            twintype: Twin type choose 1 or 2.

        Note:
            To see detail, visit _BaseTwinStructure class.
        """
        super().__init__(lattice=lattice,
                         symbol=symbol,
                         twinmode=twinmode,
                         wyckoff=wyckoff)
        self._shear_strain_ratio = None
        self._twintype = twintype
        self._rotation_matrix = None
        self._set_rotation_matrix()
        self._dichromatic_operation = None
        self._set_dichromatic_operation()
        self._layers = None
        self._delta = None

    def _set_rotation_matrix(self):
        """
        Set rotation matrix.
        """
        indices = self._indices.indices
        rotation_matrix = np.array([
                indices['m'].get_cartesian(normalize=True),
                indices['eta1'].get_cartesian(normalize=True),
                indices['k1'].get_cartesian(normalize=True),
                ])
        self._rotation_matrix = rotation_matrix

    @property
    def rotation_matrix(self):
        """
        Rotation matrix.
        """
        return self._rotation_matrix

    @property
    def shear_strain_ratio(self):
        """
        Shear twinboundary ratio.
        """
        return self._shear_strain_ratio

    @property
    def twintype(self):
        """
        Twin type.
        """
        return self._twintype

    def _set_dichromatic_operation(self):
        """
        Set dichromatic operation.

        Raises:
            ValueError: Twintype is neither equal 1 nor 2.
        """
        twintype = self._twintype
        if twintype == 1:
            W = np.array([[ 1, 0, 0],
                          [ 0, 1, 0],
                          [ 0, 0,-1]])
        elif twintype == 2:
            W = np.array([[-1, 0, 0],
                          [ 0, 1, 0],
                          [ 0, 0,-1]])
        else:
            msg = "Twintype must be 1 or 2."
            raise ValueError(msg)
        self._dichromatic_operation = W

    @property
    def dichromatic_operation(self):
        """
        Dichromatic operation.
        """
        return self._dichromatic_operation

    @property
    def layers(self):
        """
        The number of layers.
        """
        return self._layers

    @property
    def delta(self):
        """
        Delta.
        """
        return self._delta

    def _get_parent_structure(self,
                              dim:np.array) -> dict:
        """
        Get parent.
        """
        parent = ShearStructure(lattice=self._hexagonal_lattice,
                                symbol=self._symbol,
                                shear_strain_ratio=0.,
                                twinmode=self._twinmode,
                                wyckoff=self._wyckoff)
        parent.run(dim=dim)
        output_structure = parent.output_structure
        return output_structure

    def _get_shear_twinboundary_lattice(self,
                                        tb_lattice:np.array,
                                        shear_strain_ratio:float) -> np.array:
        """
        Get shear twinboudnary lattice.
        """
        lat = deepcopy(tb_lattice)
        e_b = lat[1] / np.linalg.norm(lat[1])
        shear_func = self._indices.get_shear_strain_function()
        lat[2] += np.linalg.norm(lat[2]) \
                  * shear_func(self._r) \
                  * shear_strain_ratio \
                  * e_b
        return lat

    def get_twinboundary_lattice(self,
                                 layers,
                                 delta,
                                 xshift:float=0.,
                                 yshift:float=0.):
        """
        Get twinboundary lattice.

        Args:
            layers: The number of layers in bulk.
            delta: Additional interval both sites of twin boundary.
            xshift: x shift.
            yshift: y shift.
        """
        def _add_delta(lat_cell:tuple,
                       delta:float):
            """
            Add delta interval to the twin boundary lattce.
            """
            ld = np.array([0., 0., delta])
            frame = lat_cell[0].copy()
            frame[2,:] = frame[2,:] + 4 * ld

            multi = {
                    'white_tb': 0,
                    'white': 1,
                    'black_tb': 2,
                    'black': 3,
                    }
            orig_cart_coords = np.dot(lat_cell[0].T, lat_cell[1].T).T.tolist()
            frac_coords = []
            for i in range(len(lat_cell[1])):
                cart_coord = orig_cart_coords[i] + multi[lat_cell[2][i]] * ld
                frac_coord = np.dot(cart_coord, np.linalg.inv(frame))
                frac_coords.append(frac_coord)

            return (frame, np.array(frac_coords), lat_cell[2])

        tot_layers = layers + 1
        prim_layers = self._indices.layers
        zdim = math.ceil(tot_layers / prim_layers)
        multi_dim = np.array([1,1,2*zdim])
        zratio = tot_layers / (zdim * prim_layers)
        p_orig_structure = self._get_parent_structure(dim=multi_dim)

        p_frame = np.dot(self._rotation_matrix,
                         p_orig_structure['lattice'].T).T
        t_frame = np.dot(self._dichromatic_operation,
                         p_frame.T).T
        tb_frame = np.array([
            p_frame[0],
            p_frame[1],
            (p_frame[2] - t_frame[2]) / 2
            ])
        crop_tb_frame = tb_frame * np.array([1., 1., zratio])

        p_cart_points = np.dot(p_frame.T,
                               p_orig_structure['lattice_points']['white'].T).T
        p_frac_points = np.round(np.dot(np.linalg.inv(crop_tb_frame.T),
                                        p_cart_points.T).T, decimals=8)
        t_cart_points = np.dot(self._dichromatic_operation,
                               p_cart_points.T).T
        t_frac_points = np.round(np.dot(np.linalg.inv(crop_tb_frame.T),
                                        t_cart_points.T).T, decimals=8)

        crop_p_tb_frac_points = np.array([ frac for frac in p_frac_points
                                           if np.allclose(frac[2], 0.) ]) % 1
        crop_p_frac_points = np.array([ frac for frac in p_frac_points
                                        if 0. < frac[2] < 0.5 ]) % 1
        crop_t_tb_frac_points = np.array([ frac for frac in t_frac_points
                                           if np.allclose(frac[2], -0.5) ]) % 1
        crop_t_frac_points = np.array([ frac for frac in t_frac_points
                                        if -0.5 < frac[2] < 0. ]) % 1

        p_shift = np.array([-xshift/2, -yshift/2, 0.])
        t_shift = np.array([ xshift/2,  yshift/2, 0.])
        scaled_positions = np.vstack([crop_p_tb_frac_points+p_shift,
                                      crop_p_frac_points+p_shift,
                                      crop_t_tb_frac_points+t_shift,
                                      crop_t_frac_points+t_shift])
        scaled_positions = scaled_positions[np.argsort(scaled_positions[:,2])]
        symbols = ['white_tb'] * len(crop_p_tb_frac_points) \
                + ['white'] * len(crop_p_frac_points) \
                + ['black_tb'] * len(crop_t_tb_frac_points) \
                + ['black'] * len(crop_t_frac_points)

        lat_cell = (crop_tb_frame, scaled_positions, symbols)
        delta_lat_cell = _add_delta(lat_cell, delta)

        return delta_lat_cell

    def run(self,
            layers,
            delta=0.,
            xshift:float=0.,
            yshift:float=0.,
            shear_strain_ratio:float=0.,
            make_tb_flat:bool=True,
            ):
        """
        Build structure.

        Args:
            layers: The number of layers.
            delta: Additional interval both sites of twin boundary.
            xshift: x shift.
            yshift: y shift.
            shear_strain_ratio: Shear strain ratio.
            make_tb_flat: If True, atoms on the twin boundary plane are
                          projected to twin boundary.

        Note:
            The structure built is set self.output_structure.
        """
        tb_frame, lat_points, dichs = self.get_twinboundary_lattice(
                layers=layers, delta=delta, xshift=xshift, yshift=yshift)

        W = self._dichromatic_operation
        R = self._rotation_matrix
        orig_cart_atoms = np.dot(self._hexagonal_lattice.T,
                                 self._atoms_from_lattice_points.T).T
        parent_cart_atoms = np.dot(R, orig_cart_atoms.T).T
        twin_cart_atoms = np.dot(W, parent_cart_atoms.T).T
        parent_frac_atoms = np.dot(np.linalg.inv(tb_frame.T),
                                   parent_cart_atoms.T).T
        twin_frac_atoms = np.dot(np.linalg.inv(tb_frame.T),
                                 twin_cart_atoms.T).T

        white_ix = \
                [ i for i in range(len(dichs)) if dichs[i] == 'white' ]
        white_tb_ix = \
                [ i for i in range(len(dichs)) if dichs[i] == 'white_tb' ]
        black_ix = \
                [ i for i in range(len(dichs)) if dichs[i] == 'black' ]
        black_tb_ix = \
                [ i for i in range(len(dichs)) if dichs[i] == 'black_tb' ]
        symbols = [ self._symbol ] * len(lat_points) * len(parent_frac_atoms)
        tb_shear_frame = self._get_shear_twinboundary_lattice(
                tb_lattice=tb_frame,
                shear_strain_ratio=shear_strain_ratio)

        lattice_points = {
                'white_tb': lat_points[white_tb_ix],
                'white': lat_points[white_ix],
                'black_tb': lat_points[black_tb_ix],
                'black': lat_points[black_ix]
                }

        atoms_from_lp = {
                'white': parent_frac_atoms,
                'white_tb': parent_frac_atoms.copy(),
                'black': twin_frac_atoms,
                'black_tb': twin_frac_atoms.copy(),
                }

        if make_tb_flat:
            atoms_from_lp['white_tb'] *= np.array([1., 1., 0.])
            atoms_from_lp['black_tb'] *= np.array([1., 1., 0.])

        output_structure = {
                'lattice': tb_shear_frame,
                'lattice_points': lattice_points,
                'atoms_from_lattice_points': atoms_from_lp,
                'symbols': symbols,
                }

        self._output_structure = output_structure
        self._natoms = len(self._output_structure['symbols'])
        self._xshift = xshift
        self._yshift = yshift
        self._shear_strain_ratio = shear_strain_ratio
        self._layers = layers
        self._delta = delta

    def get_translation_vectors(self) -> tuple:
        """
        Get translation vectors for white and black lattices.

        Return:
            tuple: (white_tl_vec, black_tl_vec).
        """
        tl_vecs = []
        lps = self._output_structure['lattice_points']
        w_tl_vec = lps['white'][0] - lps['white'][1]
        b_tl_vec = lps['black'][1] - lps['black'][0]
        tl_vecs = (w_tl_vec, b_tl_vec)

        return tl_vecs

def get_twinboundary(lattice:np.array,
                     symbol:str,
                     twinmode:str,
                     layers:int,
                     wyckoff:str='c',
                     delta:float=0.,
                     twintype:int=1,
                     xshift:float=0.,
                     yshift:float=0.,
                     shear_strain_ratio:float=0.,
                     expansion_ratios:np.array=np.ones(3),
                     make_tb_flat:bool=True,
                     ) -> TwinBoundaryStructure:
    """
    Get twinboudnary structure object.

    Args:
        lattice: Lattice matrix.
        symbol: Element symbol.
        twinmode: Twinmode.
        layers: The number of layers.
        wyckoff: No.194 Wycoff position ('c' or 'd').
        delta: Additional interval both sites of twin boundary.
        twintype: Twintype, choose from 1 and 2.
        xshift: x shift.
        yshift: y shift.
        shear_strain_ratio (float): Shear twinboundary ratio.
        expansion_ratios: Expansion ratios.
        make_tb_flat: If True, atoms on the twin boundary plane are
                      projected to twin boundary.
    """
    tb = TwinBoundaryStructure(lattice=lattice,
                               symbol=symbol,
                               twinmode=twinmode,
                               twintype=twintype,
                               wyckoff=wyckoff)
    tb.set_expansion_ratios(expansion_ratios)
    tb.run(layers=layers,
           delta=delta,
           xshift=xshift,
           yshift=yshift,
           shear_strain_ratio=shear_strain_ratio,
           make_tb_flat=make_tb_flat)
    return tb

def plot_nearest_atomic_distance_of_twinboundary(
        lattice:np.array,
        symbol:str,
        twinmode:str,
        layers:int,
        wyckoff:str='c',
        delta:float=0.,
        twintype:int=1,
        xshift:float=0.,
        yshift:float=0.,
        shear_strain_ratio:float=0.,
        expansion_ratios:np.array=np.ones(3),
        make_tb_flat:bool=True,
        ):
    """
    Show nearest atomic distance of twinboundary
    by changing xshift and yshift.

    Args:
        lattice: Lattice matrix.
        symbol: Element symbol.
        twinmode: Twinmode.
        layers: The number of layers.
        wyckoff: No.194 Wycoff position ('c' or 'd').
        delta: Additional interval both sites of twin boundary.
        twintype: Twintype, choose from 1 and 2.
        shear_strain_ratio (float): Shear twinboundary ratio.
        expansion_ratios: Expansion ratios.
        make_tb_flat: If True, atoms on the twin boundary plane are
                      projected to twin boundary.
    """
    import matplotlib.pyplot as plt
    from twinpy.structure.bonding import get_nearest_atomic_distance
    from twinpy.structure.lattice import CrystalLattice

    tb = TwinBoundaryStructure(lattice=lattice,
                               symbol=symbol,
                               twinmode=twinmode,
                               twintype=twintype,
                               wyckoff=wyckoff)
    tb.set_expansion_ratios(expansion_ratios)

    x = np.arange(0, 1.025, 0.025)
    y = np.arange(0, 1.025, 0.025)

    X, Y = np.meshgrid(x, y)
    shape = X.shape
    Z = np.zeros(shape)

    for i in range(shape[0]):
        for j in range(shape[1]):
            tb.run(layers=layers,
                   delta=delta,
                   xshift=X[i,j],
                   yshift=Y[i,j],
                   shear_strain_ratio=shear_strain_ratio,
                   make_tb_flat=make_tb_flat)
            cell = tb.get_cell_for_export()
            Z[i,j] = get_nearest_atomic_distance(cell)

    a, b, _ = CrystalLattice(lattice=cell[0]).abc
    fig = plt.figure(figsize=(4,4*b/a))
    cont=plt.contour(X,Y,Z,  5, vmin=0,vmax=4, colors=['black'])
    cont.clabel(fmt='%1.1f')

    plt.xlabel('xshift')
    plt.ylabel('yshift')

    plt.pcolormesh(X,Y,Z, cmap='cool')
    pp=plt.colorbar (orientation="vertical")
    pp.set_label("Nearest Atomic Distance")

    plt.show()


text:
#!/usr/bin/env python2.7
# encoding: utf-8
"""
uniform_order_statistics.py

Created by Jakub Konka on 2011-04-18.
Copyright (c) 2011 University of Strathclyde. All rights reserved.
"""

import sys
import os
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def plot_pdf(order, N, iterations):
	order_stats = []
	for it in range(iterations):
		numbers = [np.random.uniform(0,1) for i in range(N)]
		numbers.sort()
		order_stats.append(numbers[order-1])
	plt.figure()
	n, bins, patches = plt.hist(order_stats, iterations/20, normed=1, facecolor='green')
	y = lambda x: int(sp.factorial(N))/(int(sp.factorial(N-order))*int(sp.factorial(order-1))) * x**(order-1) * (1-x)**(N-order)
	plt.plot(bins, y(bins), 'r--', linewidth=3)

if __name__ == '__main__':
	iterate = 10000
	N = 5
	for i in range(N):
		plot_pdf(i+1,N,iterate)
	plt.show()

max_stars_repo_name: kubkon/Phd-python

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: kubkon/Phd-python

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: kubkon/Phd-python

max_forks_repo_licenses: MIT Unlicense

reason: import numpy,import scipy


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

import os, sys
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt

try:
    from vafnet import VAFnet, Visualizer
except:
    print("vafnet is not installed yet! Trying to call directly from source...")
    from sys import path
    path.append("../../")
    from vafnet import VAFnet, Visualizer
    print("done.")

BEST_MODEL_PATH = 'best_model.pkl'

############################ Prior Setting ############################
archits = [[30, 20, 10], [35, 15], [35]]
reps = 5
plot_metric = 'mse'
select_params_metric = 'mse'
select_model_metric = 'score'
visualizer = None
# fig = plt.figure(figsize=(8, 6), facecolor='white')
# visualizer = Visualizer(fig, plot_metric)
algo = {
    'algo': 'adadelta',
    'algo_params': {
        'learning_rate':1e-2,
        'rho':0.75,
    }
}
opt_params = {
    'obj': select_params_metric,
    'algo': algo,
    'cv_nfolds': 3,
    'cvrg_tol': 1e-5,
    'max_cvrg': 50,
    'max_iter': 1000
}
evals = {
    'score': [
        'Model Selection Score',
        {"("+",".join([str(D) for D in archit])+")": [] for archit in archits}
    ],
    'obj': [
        'Params Optimization Objective',
        {"("+",".join([str(D) for D in archit])+")": [] for archit in archits}
    ],
    'mse': [
        'Mean Square Error',
        {"("+",".join([str(D) for D in archit])+")": [] for archit in archits}
    ],
    'rmse': [
        'Root Mean Square Error',
        {"("+",".join([str(D) for D in archit])+")": [] for archit in archits}
    ],
    'nmse': [
        'Normalized Mean Square Error',
        {"("+",".join([str(D) for D in archit])+")": [] for archit in archits}
    ],
    'mnlp': [
        'Mean Negative Log Probability',
        {"("+",".join([str(D) for D in archit])+")": [] for archit in archits}
    ],
    'time': [
        'Training Time(s)',
        {"("+",".join([str(D) for D in archit])+")": [] for archit in archits}
    ],
}

############################ General Methods ############################
def debug(local):
    locals().update(local)
    print('Debug Commands:')
    while True:
        cmd = input('>>> ')
        if(cmd == ''):
            break
        try:
            exec(cmd)
        except Exception as e:
            import traceback
            traceback.print_tb(e.__traceback__)

def plot_dist(*args):
    import seaborn as sns
    for x in args:
        plt.figure()
        sns.distplot(x)
    plt.show()

def load_data(prop=0.9):
    from sklearn import datasets
    boston = datasets.load_boston()
    X, y = boston.data, boston.target
    y = y[:, None]
    ntrain = y.shape[0]
    train_inds = npr.choice(range(ntrain), int(prop*ntrain), replace=False)
    test_inds = np.setdiff1d(range(ntrain), train_inds)
    X_train, y_train = X[train_inds].copy(), y[train_inds].copy()
    X_test, y_test = X[test_inds].copy(), y[test_inds].copy()
    return X_train.T, y_train.T, X_test.T, y_test.T

############################ Training Phase ############################
for archit in archits:
    funcs = None
    results = {en:[] for en in evals.keys()}
    for round in range(reps):
        X_train, y_train, X_test, y_test = load_data()
        model = VAFnet(archit)
        model.optimize(X_train, y_train, funcs, visualizer, **opt_params)
        if(funcs is None):
            funcs = model.get_compiled_funcs()
        if(not os.path.exists(BEST_MODEL_PATH)):
            model.save(BEST_MODEL_PATH)
        else:
            best_model = VAFnet(archit).load(BEST_MODEL_PATH)
            best_model.fit(X_train, y_train)
            best_model.score(X_test, y_test)
            model.score(X_test, y_test)
            best_model._print_evals_comparison(model.evals)
            if(model.evals[select_model_metric][1][-1] <
                best_model.evals[select_model_metric][1][-1]):
                model.save(BEST_MODEL_PATH)
                print("!"*80)
                print("!"*30, "NEW BEST PREDICTOR", "!"*30)
                print("!"*80)
        if(round >= reps//2):
            for res in evals.keys():
                results[res].append(model.evals[res][1][-1])
    for en in evals.keys():
        eval = (np.mean(results[en]), np.std(results[en]))
        evals[en][1]["("+",".join([str(D) for D in archit])+")"].append(eval)

############################ Plot Performances ############################
    # import os
    # if not os.path.exists('plots'):
    #     os.mkdir('plots')
    # fig = plt.figure(facecolor='white', dpi=120)
    # ax = fig.add_subplot(111)
    # for en, (metric_name, metric_result) in evals.items():
    #     maxv, minv = 0, 1e5
    #     for archit in metric_result.keys():
    #         for j in range(i+1):
    #             maxv = max(maxv, metric_result[archit][j][0]+\
    #                 metric_result[archit][j][1])
    #             minv = min(minv, metric_result[archit][j][0]-\
    #                 metric_result[archit][j][1])
    #         line = ax.errorbar(nfeats_choice[:i+1],
    #             [metric_result[archit][j][0] for j in range(i+1)],
    #             yerr=[metric_result[archit][j][1] for j in range(i+1)],
    #             fmt='o', capsize=6, label=archit, alpha=0.6)
    #     ax.set_xlim([min(nfeats_choice)-10, max(nfeats_choice)+10])
    #     ax.set_ylim([minv-abs(maxv-minv)*0.2,maxv+abs(maxv-minv)*0.2])
    #     ax.grid(True)
    #     plt.title(metric_name, fontsize=18)
    #     handles, labels = ax.get_legend_handles_labels()
    #     ax.legend(handles, labels, loc='upper right', ncol=3, fancybox=True)
    #     plt.xlabel('Optimized Features', fontsize=13)
    #     plt.ylabel(en, fontsize=13)
    #     plt.savefig('plots/'+en.lower()+'.png')
    #     ax.cla()

max_stars_repo_name: MaxInGaussian/VAFnet

max_issues_repo_name: MaxInGaussian/VAFnet

max_forks_repo_name: MaxInGaussian/VAFnet


text:
import numpy as np
import coord_transforms as ct
import constants
import scipy.interpolate as scint

magn = 2.16
a = 1/10 # paraboloid scaling factor
#a = 1/10*magn
xpixels = 1024;
ypixels = 1024;
pixelsize = 13*10**(-3)
pixelarea = pixelsize**2
xmax = 10.75/magn
focusdistance = 0.5/magn
xcut = 10.75/magn
ni = 1;

def parabola_position(direction):
    '''
    direction: emission direction vector in Cartesian coordinates, a numpy
    array of shape (N by 3)
    '''
    direction = direction/np.expand_dims(ct.field_magnitude(direction), -1)
    r, theta, phi = ct.cartesian_to_spherical_coords(direction)
    c = np.empty(np.shape(direction))
    positive_x_condition = np.logical_not(
        np.logical_and(theta==np.pi/2, phi==0)
        )
    c[np.logical_not(positive_x_condition)] = np.nan
    c[positive_x_condition] = np.expand_dims(
            np.array(
                1/(
                    2*a*(
                        1-np.sin(theta[positive_x_condition])
                        *np.cos(phi[positive_x_condition])
                        )
                    )
                ), -1
        )
    parabola = direction * c
    return parabola

def parabola_normals(parabola_positions):
    '''
    parabola_positions: positions on the parabola, N by 3 numpy array
    R = (a*r**2, r*cos(theta), r*sin(theta))
    '''
#    r, theta = ct.cartesian_to_polar(parabola_positions[:, 1], parabola_positions[:, 2])

#    ## One way to calculate the normals
#    normalizing_factor = 1/np.sqrt(1 + 4 * a**2 * r**2)
#    normal_x = np.ones(np.shape(r))
#    normal_y = -2 * a * r * np.cos(theta)
#    normal_z = -2 * a * r * np.sin(theta)
#    normals = np.transpose(normalizing_factor * np.vstack((normal_x, normal_y, normal_z)))
#
    # parameterize x in terms of y and z: x = (a*(y**2 + z**2)-2.5, y, z)
    normalizing_factor = 1/np.sqrt(
        1 + 4 * a**2 * (
            parabola_positions[:, 1]**2 + parabola_positions[:, 2]**2
            )
        )
    normal_x = np.ones(np.shape(parabola_positions)[0])
    normal_y = -2 * a * parabola_positions[:, 1]
    normal_z = -2 * a * parabola_positions[:, 2]
    normals = normalizing_factor * np.transpose(
        np.vstack((normal_x, normal_y, normal_z))
        )

    ## Another way to calculate the normals
#    # derivative of parabola along r
#    r_gradient_x = 2*a*r
#    r_gradient_y = np.cos(theta)
#    r_gradient_z = np.sin(theta)
#    r_gradient = np.vstack((r_gradient_x, r_gradient_y, r_gradient_z)).transpose()
#
#    # derivative of parabola along theta
#    theta_gradient_x = np.zeros(np.shape(r_gradient_x))
#    theta_gradient_y = -r * np.sin(theta)
#    theta_gradient_z = r * np.cos(theta)
#    theta_gradient = np.vstack((theta_gradient_x, theta_gradient_y, theta_gradient_z)).transpose()
#
#    normals = np.cross(r_gradient, theta_gradient) #cross product of r and theta to get the surface normal
#    normals = normals/ct.field_magnitude(normals)[:, None] # normalize the normal vectors
    return normals

def surface_polarization_directions(theta, phi):
    '''
    calculate direction vectors for p- and s- polarized light on the mirror surface
    p-polarized: polarized along theta (is this true for the paraboloid??)
    s-polarized: polarized along phi (is this true for the paraboloid??)
    '''
    p_direction = np.vstack(ct.spherical_to_cartesian_vector_field(
        theta,
        phi,
        np.zeros(np.shape(theta)),
        np.ones(np.shape(theta)),
        np.zeros(np.shape(theta)),
        ))

    s_direction = np.vstack(ct.spherical_to_cartesian_vector_field(
        theta,
        phi,
        np.zeros(np.shape(theta)),
        np.zeros(np.shape(theta)),
        np.ones(np.shape(theta)),
        ))
    return p_direction, s_direction

def fresnel_reflection_coefficients(normals, e_incident_direction, n_mirror, n_environment=1):
    '''
    calculate s and p fresnel reflection coefficients for the parabolic mirror surface
    normals: surface normals of the mirror
    e_incident_direction: direction vector of the incident electromagnetic wave
    n_environment: refractive index of the environment at the desired wavelength
    n_mirror: refractive index of the mirror at the desired wavelength
    '''
    incidence_angle = np.arccos(
        np.dot(normals, e_incident_direction)
        ) / (
            ct.field_magnitude(normals)*ct.field_magnitude(e_incident_direction)
        )
    n_factor = np.sqrt(1-np.square(n_environment/n_mirror * np.sin(incidence_angle)))
    R_s = np.square(
        (
            n_environment * np.cos(incidence_angle) -
            n_mirror * n_factor
        ) / (
            n_environment * np.cos(incidence_angle) +
            n_mirror * n_factor
            )
        )
    R_p = np.square(
        (
            n_environment * n_factor -
            n_mirror * np.cos(incidence_angle)
        ) / (
            n_environment * n_factor +
            n_mirror * np.cos(incidence_angle)
            )
        )
    return R_s, r_p

def aluminium_refractive_index(wavelength, filename):
    '''
    import data file containing aluminium's dielectric function
    wavelength: wavelength at which it is desired to extract the refractive index (in m)
    filename: file containing the dielectric function (should point to al_pa.mat)
    returns complex refractive index at the given wavelength
    '''
    load_data = np.loadtxt(filename, skiprows=15, max_rows=141-15) # must fix that to allow other files than al_pa.mat
    index_factor = np.sqrt(load_data[:, 1]**2 + load_data[:, 2]**2)
    n = np.sqrt(0.5 * (index_factor + load_data[:, 1])) + 1j * np.sqrt(0.5 * (index_factor - load_data[:, 1]))
    wavelength_list = constants.PLANCK * constants.LIGHTSPEED / (load_data[:, 0] * constants.COULOMB)
    interp_function_n = scint.interp1d(wavelength_list, n, kind='cubic')
    n_wavelength = interp_function_n(wavelength)
    return n_wavelength

max_stars_repo_name: icbicket/CLFields

max_issues_repo_name: icbicket/CLFields

max_forks_repo_name: icbicket/CLFields


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Jan 04 13:23:27 2013

@author: Joey
"""
import numpy as np

import spectroscopy.Spacetime as Spacetime
import spectroscopy.ElectronicOperator as ElectronicOperator
import spectroscopy.NuclearWavefunction as NuclearWavefunction
import spectroscopy.NuclearOperator as NuclearOperator
import spectroscopy.ElectronicWavefunction as ElectronicWavefunction

omega_1_g_wavenumbers = 640.0
omega_1_e_wavenumbers = omega_1_g_wavenumbers * (.90)
HR_1 = 1.0
mass_1_amu = 63.466

c_natural = .20
c_min_log = -3. #np.log10(2.06913808e-01) #
c_max_log = 1.
c_number = 20
c_values = np.logspace(c_min_log, c_max_log, c_number)
c_values
pulse_carrier_frequency_wavenumbers = 0.0

##all units from here on out are defined by mySpace
mySpace = Spacetime.Spacetime(xMax_tuple = (35.0,),
             numberOfNuclearDimenions = 1,
             numberOfElectronicDimensions = 2,
             numberOfSimulationSpacePointsPerNuclearDimension_tuple = (28,),
             dt_SECONDS = .300E-15,
             UnityMassInElectronMasses = 100.0)

omega_1_g = mySpace.unitHandler.energyUnitsFromWavenumbers(omega_1_g_wavenumbers)
omega_1_e = mySpace.unitHandler.energyUnitsFromWavenumbers(omega_1_e_wavenumbers)

m_1 = mySpace.unitHandler.massUnitsFromAmu(mass_1_amu)

pulse_carrier_frequency = mySpace.unitHandler.energyUnitsFromWavenumbers(pulse_carrier_frequency_wavenumbers)

dx_1 = np.sqrt(HR_1 * 2.0 * mySpace.hbar / (m_1 * omega_1_e))

x_0 = 0.0
groundCenter = x_0
excitedCenter = x_0 + dx_1

nuclear_ground_1 = NuclearOperator.harmonicOscillator(mySpace,
                                                 omega = omega_1_g,
                                                 mass = m_1,
                                                 center = groundCenter)
nuclear_excited_1 = NuclearOperator.harmonicOscillator(mySpace,
                                                 omega = omega_1_e,
                                                 mass = m_1,
                                                 center = excitedCenter)

x_max_needed = 2 * dx_1 + 4.0 *  nuclear_ground_1.sigma
print "xmax we need: ", x_max_needed, "  Versus what we got: ", mySpace.xMax_values[0]
dx_needed_2_interactions = (mySpace.hbar * np.pi) / (nuclear_ground_1.omega * nuclear_ground_1.mass * (2.0 * dx_1 + 4.0 *  nuclear_ground_1.sigma))
print "dx needed for 2 interactions: ", dx_needed_2_interactions, "but we have: ", mySpace.Dx_values[0]

electronic_ground = NuclearOperator.nuclearHamiltonian(mySpace, listOfOneDimensionalHamiltonians = [nuclear_ground_1] )

electronic_excited = NuclearOperator.nuclearHamiltonian(mySpace, listOfOneDimensionalHamiltonians = [nuclear_excited_1 ] )

groundStateNuclearWF = NuclearWavefunction.nuclearWavefunction(mySpace, groundStateNuclearHamiltonian = electronic_ground )

ElectronicHamiltonian = ElectronicOperator.ElectronicHamiltonian(mySpace, [(0,0, electronic_ground),
                                                            (1,1, electronic_excited)])

initialEWF = ElectronicWavefunction.electronicWavefunction(mySpace,
                       listOfNuclearWavefunctions = [groundStateNuclearWF, 0],
                       Normalize=True)

mu_0 = 1.0

x_max = mySpace.xMax_values[0]

constantMu = NuclearOperator.constantPositionNuclearOperator(mySpace, mu_0)

xTransitionDipole_FC = ElectronicOperator.ElectronicPositionOperator(mySpace, [(0, 1, constantMu), (1, 0, constantMu)])
yTransitionDipole_FC = ElectronicOperator.ElectronicPositionOperator(mySpace, [])
zTransitionDipole_FC = ElectronicOperator.ElectronicPositionOperator(mySpace, [])

transitionDipoleTuple_FC = (xTransitionDipole_FC, yTransitionDipole_FC, zTransitionDipole_FC)

#linearMu = NuclearOperator.functionalPositionNuclearOperator(mySpace, lambda x: mu_0 + c_natural * (x- x_max)/HR_1 )
linearMu = NuclearOperator.functionalPositionNuclearOperator(mySpace, lambda x: mu_0 + c_natural * (x)/dx_1 )

xTransitionDipole_linear = ElectronicOperator.ElectronicPositionOperator(mySpace, [(0, 1, linearMu ), (1, 0, linearMu)])
yTransitionDipole_linear = ElectronicOperator.ElectronicPositionOperator(mySpace, [])
zTransitionDipole_linear = ElectronicOperator.ElectronicPositionOperator(mySpace, [])

transitionDipoleTuple_linear = (xTransitionDipole_linear, yTransitionDipole_linear, zTransitionDipole_linear)

def transitionDipoleTuple_gaussian_noisy(center, stdev, seed = 42):
    np.random.seed(seed)
    noisyMu = NuclearOperator.constantPositionNuclearOperator(mySpace, center)
    noisyMu.surface = noisyMu.surface + stdev * np.random.randn(*noisyMu.surface.shape)
    x = ElectronicOperator.ElectronicPositionOperator(mySpace, [(0, 1, noisyMu ), (1, 0, noisyMu)])
    y = ElectronicOperator.ElectronicPositionOperator(mySpace, [])
    z = ElectronicOperator.ElectronicPositionOperator(mySpace, [])
    return (x, y, z)

def transitionDipoleTuple_sine_groundCenter(mu_0, amplitude, frequency):
    sineMu = NuclearOperator.functionalPositionNuclearOperator(mySpace, lambda x: mu_0 + np.sin( frequency * (x - groundCenter) ) )
    x = ElectronicOperator.ElectronicPositionOperator(mySpace, [(0, 1, sineMu ), (1, 0, sineMu)])
    y = ElectronicOperator.ElectronicPositionOperator(mySpace, [])
    z = ElectronicOperator.ElectronicPositionOperator(mySpace, [])
    return (x, y, z)

frequency_for_half_cycle_between_wells = np.pi / dx_1

c_dipole_tuples = []
for c in c_values:
    lin_mu = NuclearOperator.functionalPositionNuclearOperator(mySpace, lambda x: mu_0 + c * (x)/dx_1 )

    xD = ElectronicOperator.ElectronicPositionOperator(mySpace, [(0, 1, lin_mu ), (1, 0, lin_mu)])
    yD = ElectronicOperator.ElectronicPositionOperator(mySpace, [])
    zD = ElectronicOperator.ElectronicPositionOperator(mySpace, [])

    tupleTD = (xD, yD, zD)
    c_dipole_tuples.append(tupleTD)

max_stars_repo_name: jgoodknight/spectroscopy

max_issues_repo_name: jgoodknight/spectroscopy

max_forks_repo_name: jgoodknight/spectroscopy


text:
# Author: Jintao Huang
# Email: hjt_study@qq.com
# Date:

import numpy as np
from sklearn.linear_model import Ridge as _Ridge
import time
from dev.linear_model import Ridge

# In[0]: test run time
# In[1]: test accuracy
# In[2]: test shape
rng = np.random.default_rng(42)
X = rng.random((2000, 500), dtype=np.float32)
y = rng.random((2000, 50), dtype=np.float32)

# sklearn warm up
t = time.time()
reg = _Ridge(solver='svd')
reg.fit(X, y)
#
y_sk = reg.predict(X)
s_sk = reg.score(X, y)
c_sk = reg.coef_
i_sk = reg.intercept_
print("time: %.6f" % (time.time() - t))

# sklearn
t = time.time()
reg = _Ridge(solver='svd')
reg.fit(X, y)
reg.predict(X)
print("time: %.6f" % (time.time() - t))

# torch warm up gpu
t = time.time()
reg = Ridge(device='cuda')
reg.fit(X, y)
#
y_t = reg.predict(X).cpu().numpy()
s_t = reg.score(X, y).cpu().numpy()
c_t = reg.coef_.cpu().numpy()
i_t = reg.intercept_.cpu().numpy()
print("time: %.6f" % (time.time() - t))

# torch gpu
t = time.time()
reg = Ridge(device='cuda')
reg.fit(X, y)
reg.predict(X)
print("time: %.6f" % (time.time() - t))

# torch cpu warm up
t = time.time()
reg = Ridge()
reg.fit(X, y)
reg.predict(X)
print("time: %.6f" % (time.time() - t))

# torch cpu
t = time.time()
reg = Ridge()
reg.fit(X, y)
reg.predict(X)
print("time: %.6f" % (time.time() - t))
"""Out[0]
time: 0.388993
time: 0.371973
time: 2.902385
time: 0.062803
time: 0.056849
time: 0.049865
"""
# judge correctness
print(np.allclose(y_sk, y_t, rtol=1e-4, atol=1e-5))
print(np.allclose(s_sk, s_t, rtol=1e-4, atol=1e-6))
print(np.allclose(c_sk, c_t, rtol=1e-4, atol=1e-6))
print(np.allclose(i_sk, i_t, rtol=1e-4, atol=1e-5))
"""Out[1]
True
True
True
True
"""
# test shape
print(y_sk.shape, y_t.shape)
print(c_sk.shape, c_t.shape)
print(i_sk.shape, i_t.shape)
"""Out[2]
(20000, 50) (20000, 50)
(50, 500) (50, 500)
(50,) (50,)
"""

# In[3]: Differences from SkLearn: shape
rng = np.random.default_rng(42)
X = rng.random((2, 1), dtype=np.float32)
y = rng.random((2,), dtype=np.float32)

# sklearn
reg = _Ridge()
reg.fit(X, y)
#
y_sk = reg.predict(X)
s_sk = reg.score(X, y)
c_sk = reg.coef_
i_sk = reg.intercept_

# torch gpu
reg = Ridge(device='cuda')
reg.fit(X, y)
#
y_t = reg.predict(X).cpu().numpy()
s_t = reg.score(X, y).cpu().numpy()
c_t = reg.coef_.cpu().numpy()
i_t = reg.intercept_.cpu().numpy()
print(y_sk, y_t)
print(s_sk, s_t)
print(c_sk, c_t)
print(i_sk, i_t)
"""Out[3]
[0.5672046 0.5262452] [[0.5672046]
 [0.5262452]]
0.34373288682296077 0.3437329
[-0.0598205] [[-0.05982049]]
0.5725436 [0.5725436]
"""

max_stars_repo_name: ustcml/TorchML

max_issues_repo_name: ustcml/TorchML

max_forks_repo_name: ustcml/TorchML