text:
from __future__ import print_function
from __future__ import absolute_import

import unittest
import qgate
from qgate.script import *
import numpy as np

class TestSampling(unittest.TestCase) :

    def test_observation(self) :
        refs = new_references(3)
        obs = qgate.simulator.observation.Observation(refs, 0b100, 0b001)
        self.assertEqual(repr(obs), '10*')
        self.assertEqual(obs.int, 4)

        self.assertIs(obs(refs[0]), None)
        self.assertEqual(obs(refs[1]), 0)
        self.assertEqual(obs(refs[2]), 1)

    def test_observation_eq(self) :
        refs = new_references(3)
        obs0 = qgate.simulator.observation.Observation(refs, 0b100, 0b001)
        self.assertTrue(obs0 == obs0)

        obs1 = qgate.simulator.observation.Observation(refs, 0b100, 0b000)
        self.assertFalse(obs0 == obs1)

        obs1 = qgate.simulator.observation.Observation(refs[0:1], 0b100, 0b001)
        self.assertFalse(obs0 == obs1)

        obs1 = qgate.simulator.observation.Observation(refs, 0b110, 0b001)
        self.assertFalse(obs0 == obs1)

    def test_observation_lt(self) :
        refs = new_references(3)
        obs = qgate.simulator.observation.Observation(refs, 0b100, 0b001)
        with self.assertRaises(TypeError) :
            v = obs < obs
        with self.assertRaises(TypeError) :
            v = obs < self

    def test_observation_list(self) :
        refs = new_references(3)
        obs = qgate.simulator.observation.Observation(refs, 0b100, 0b001)
        with self.assertRaises(TypeError) :
            v = obs < obs
        with self.assertRaises(TypeError) :
            v = obs < self

    def test_sampling_one_shot(self) :
        qregs = new_qregs(4)
        cregs = new_references(4)
        circuit = [I(qreg) for qreg in qregs]

        sim = qgate.simulator.py()
        sim.run(circuit)
        obs = sim.obs(cregs)
        self.assertEqual(int(obs), 0)

    def test_sampling_128_shots(self) :
        qregs = new_qregs(4)
        cregs = new_references(4)
        circuit = [[I(qreg) for qreg in qregs],
                   X(qregs[0]),
                   [measure(creg, qreg) for creg, qreg in zip(cregs[0:-1], qregs[0:-1])]]

        sim = qgate.simulator.cpu()
        obslist = sim.sample(circuit, cregs, 128)
        # sample list length
        self.assertEqual(len(obslist), 128)
        # sampled value
        obs = obslist[0]
        self.assertTrue(isinstance(obs, qgate.simulator.observation.Observation))
        self.assertEqual(repr(obs), "*001")

        # intarray
        values = obslist.intarray
        self.assertEqual(len(values), 128)
        self.assertTrue(np.all(values == 1))

        # __call__ to extract
        # creg[0] -> 1
        extracted = obslist(cregs[0])
        self.assertEqual(len(extracted), 128)
        b = [v == 1 for v in extracted]
        self.assertTrue(np.all(b))
        # creg[1] -> 0
        extracted = obslist(cregs[1])
        b = [v == 0 for v in extracted]
        self.assertTrue(np.all(b))
        # creg[3] -> None
        extracted = obslist(cregs[3])
        b = [v is None for v in extracted]
        self.assertTrue(np.all(b))

        # __getitem__
        half = obslist[0:64]
        self.assertEqual(len(half), 64)

    def test_observation_list_getitem(self) :
        cregs = new_references(4)
        values = np.random.random((64))
        obslist = qgate.simulator.observation.ObservationList(cregs, values, 0)
        obslist32 = obslist[:32]
        self.assertEqual(obslist32._reflist, obslist._reflist)
        self.assertTrue(np.allclose(obslist32._values, obslist._values[:32]))

    def test_observation_histgram(self) :
        qregs = new_qregs(4)
        cregs = new_references(4)
        circuit = [[I(qreg) for qreg in qregs],
                   X(qregs[0]),
                   [measure(creg, qreg) for creg, qreg in zip(cregs[0:-1], qregs[0:-1])]]

        sim = qgate.simulator.cpu()
        obslist = sim.sample(circuit, cregs, 128)

        hist = obslist.histgram()
        self.assertEqual(hist.n_samples, 128)
        self.assertEqual(len(hist), 1)
        self.assertEqual(hist[0], 0)
        self.assertEqual(hist[1], 128)

        keys = [key for key in hist.keys()]
        self.assertEqual(len(keys), 1)
        self.assertEqual(keys[0], 1)

        values = [value for value in hist.values()]
        self.assertEqual(len(values), 1)
        self.assertEqual(values[0], 128)

        copied = dict(hist.items())
        self.assertEqual(copied, hist._hist)

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

max_stars_repo_name: rickyHong/Gate-Quantum-gpu-repl

max_issues_repo_name: rickyHong/Gate-Quantum-gpu-repl


text:
import rosbag
import argparse
import numpy as np
parser = argparse.ArgumentParser()
parser.add_argument('--bag_file', type=str)

opt = parser.parse_args()

ODOM_TOPIC = '/airsim_node/PhysXCar/odom_local_ned'
CMD_TOPIC = '/airsim_node/PhysXCar/cmd_vel'

bag = rosbag.Bag(opt.bag_file)

active_cmd = None
cmd_start_time = 0
from matplotlib import pyplot as plt

command_responses = []

for topic, msg, t in bag.read_messages(topics=[ODOM_TOPIC, CMD_TOPIC]):
  if topic == CMD_TOPIC:
    active_cmd = msg
    cmd_start_time = t
    command_responses.append([])
  elif topic == ODOM_TOPIC and active_cmd:
    x_vel = np.linalg.norm([msg.twist.twist.linear.x, msg.twist.twist.linear.y, msg.twist.twist.linear.z])
    target_vel = active_cmd.linear.x
    command_responses[-1].append((t.to_sec(), x_vel, target_vel))

plt.figure()
for response_info in command_responses:
  data = np.array(response_info)
  if len(data):
    plt.scatter(data[:, 0], data[:, 1])
    plt.plot(data[:, 0], data[:, 2])
plt.show()

bag.close()

ext: py

max_stars_repo_name: ut-amrl/AirSim

max_issues_repo_name: ut-amrl/AirSim

max_forks_repo_name: ut-amrl/AirSim

reason: import numpy


text:
# Copyright (C) 2013 Joshua L. Adelman, Carsen A. Stringer and Daniel M. Zuckerman
#
# This file is part of WESTPA.
#
# WESTPA is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# WESTPA is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with WESTPA.  If not, see <http://www.gnu.org/licenses/>.

import numpy as np
import numpy.ma as ma
import itertools

TOL=1.0E-6

class UncertContainer(object):
    """ Container to hold uncertainty measurements. Data is convert to np masked arrays
        to avoid possible numerical problems
    """
    def __init__(self,vals,vals_dmin,vals_dmax,mask=ma.nomask):
        super(UncertContainer, self).__init__()

        # If input data already masked arrays extract unmasked data
        if ma.isMaskedArray(vals):
            vals = vals.data
        if ma.isMaskedArray(vals_dmin):
            vals_dmin = vals_dmin.data
        if ma.isMaskedArray(vals_dmax):
            vals_dmax = vals_dmax.data

        # Adjust negative values
        ineg = np.where(vals_dmin <= 0.0)
        vals_dmin[ineg] = TOL*vals[ineg]

        # Calculate weight based on fractional uncertainty
        diff = vals_dmax - vals_dmin
        diff_m = ma.masked_where(vals_dmax == vals_dmin,diff)

        self.vals = ma.masked_where(vals == 0.0,vals)

        self.wt = (self.vals/diff_m)**2
        self.uncert = diff_m/self.vals

        self.wt.fill_value = np.inf
        self.uncert.fill_vaule = np.inf

        assert np.all(self.wt.mask == self.uncert.mask)

        # Mask data if uncertainty is not finite or if any of the inputs were
        # already masked

        mm = ma.mask_or(self.wt.mask,mask)

        self.vals.mask = mm
        self.wt.mask = mm
        self.uncert.mask = mm
        self.dmin = ma.array(vals_dmin,mask=mm,fill_value=np.inf)
        self.dmax = ma.array(vals_dmax,mask=mm,fill_value=np.inf)

        self.mask = ma.getmaskarray(self.vals)

    def __getitem__(self,indx):
        vals = self.vals[indx]
        dmin = self.dmin[indx]
        dmax = self.dmax[indx]

        if type(vals) is ma.core.MaskedConstant:
            dum = np.zeros((1,))
            return UncertContainer(dum.copy(),dum.copy(),dum.copy())
        elif isinstance(vals,(float,int,np.float,np.int)):
            return UncertContainer(np.array([vals,]),np.array([dmin,]),np.array([dmax,]))
        elif isinstance(vals,np.ndarray):
            return UncertContainer(vals,dmin,dmax,mask=vals.mask)
        else:
            raise TypeError

    def __setitem__(self,indx,value):
        if isinstance(value,UncertContainer):
            self.vals[indx] = value.vals
            self.dmin[indx] = value.dmin
            self.dmax[indx] = value.dmax
            self.wt[indx] = value.wt
            self.uncert[indx] = value.uncert
        else:
            raise TypeError('Can only set values with an UncertContainer object')

    def __repr__(self):
        return 'shape={} vals={} dmin={} dmax={} vals.mask={}'.format(self.vals.shape,self.vals,self.dmin,self.dmax,self.vals.mask)

    def __add__(self,value):
        if isinstance(value,UncertContainer):
            vals = self.vals + value.vals
            dmin = self.dmin + value.dmin
            dmax = self.dmax + value.dmax

            return UncertContainer(vals,dmin,dmax,mask=vals.mask)
        elif isinstance(value,(float,int,np.float,np.int)):
            vals = self.vals + value
            dmin = self.dmin + value
            dmax = self.dmax + value

            return UncertContainer(vals,dmin,dmax,mask=vals.mask)
        else:
            raise TypeError('Attempt to add value of unsupported type')

    def __sub__(self,value):
        if isinstance(value,UncertContainer):
            vals = self.vals - value.vals
            dmin = self.dmin - value.dmin
            dmax = self.dmax - value.dmax

            return UncertContainer(vals,dmin,dmax,mask=vals.mask)
        else:
            raise TypeError ('Attempted to subtract by value of unsupported type')

    def __mul__(self,value):
        if isinstance(value,UncertContainer):
            vals = self.vals * value.vals
            dmin = self.dmin * value.dmin
            dmax = self.dmax * value.dmax

            return UncertContainer(vals,dmin,dmax,mask=vals.mask)

        elif isinstance(value,(float,int,np.float,np.int)):
            vals = self.vals * value
            dmin = self.dmin * value
            dmax = self.dmax * value

            return UncertContainer(vals,dmin,dmax,mask=vals.mask)
        else:
            raise TypeError('Attempted to multiply by value of unsupported type')

    def __div__(self,value):
        if isinstance(value,UncertContainer):
            vals = self.vals / value.vals
            dmin = self.dmin / value.dmax
            dmax = self.dmax / value.dmin

            return UncertContainer(vals,dmin,dmax,mask=vals.mask)
        else:
            raise TypeError('Attempted to divide by unsupported type')

    def transpose(self):
        vals = self.vals.T
        dmin = self.dmin.T
        dmax = self.dmax.T

        return UncertContainer(vals,dmin,dmax,mask=vals.mask)

    def recip(self):
        vals = 1.0 / self.vals
        dmin = 1.0 / self.dmax
        dmax = 1.0 / self.dmin

        return UncertContainer(vals,dmin,dmax,mask=vals.mask)

    def update_mask(self):
        self.vals.mask = self.mask
        self.dmin.mask = self.mask
        self.dmax.mask = self.mask
        self.wt.mask = self.mask
        self.uncert.mask = self.mask

    def concatenate(self,value,axis=0):
        """ Concatentate UncertContainer value to self.
            Assumes that if dimensions of self and value do not match, to
            add a np.newaxis along axis of value
        """

        if isinstance(value,UncertContainer):
            if value.vals.ndim == self.vals.ndim:
                vals = value.vals
                dmin = value.dmin
                dmax = value.dmax
                wt = value.wt
                uncert = value.uncert
                mask = value.mask
            elif (value.vals.ndim + 1) == self.vals.ndim:
                vals =  ma.expand_dims(value.vals,axis)
                dmin =  ma.expand_dims(value.dmin,axis)
                dmax =  ma.expand_dims(value.dmax,axis)
                wt =  ma.expand_dims(value.wt,axis)
                uncert =  ma.expand_dims(value.uncert,axis)
                mask =  np.expand_dims(value.mask,axis)
            else:
                raise ValueError('Could not propery match dimensionality')

            self.vals = ma.concatenate((self.vals,vals),axis=axis)
            self.dmin = ma.concatenate((self.dmin,dmin),axis=axis)
            self.dmax = ma.concatenate((self.dmax,dmax),axis=axis)
            self.wt = ma.concatenate((self.wt,wt),axis=axis)
            self.uncert = ma.concatenate((self.uncert,uncert),axis=axis)

            self.mask = np.concatenate((self.mask,mask),axis=axis)
        else:
            raise ValueError('Can only concatenate with an UncertContainer object')

    def weighted_average(self,axis=0,expaxis=None):
        """ Calculate weighted average of data along axis
            after optionally inserting a new dimension into the
            shape array at position expaxis
        """

        if expaxis is not None:
            vals = ma.expand_dims(self.vals,expaxis)
            dmin = ma.expand_dims(self.dmin,expaxis)
            dmax = ma.expand_dims(self.dmax,expaxis)
            wt = ma.expand_dims(self.wt,expaxis)
        else:
            vals = self.vals
            wt = self.wt
            dmin = self.dmin
            dmax = self.dmax

        # Get average value
        avg,norm = ma.average(vals,axis=axis,weights=wt,returned=True)
        avg_ex = ma.expand_dims(avg,0)

        # Calculate weighted uncertainty
        wtmax = ma.max(wt,axis=axis)
        neff = norm/wtmax       # Effective number of samples based on uncertainties

        # Seeking max deviation from the average; if above avg use max, if below use min
        term = np.empty_like(vals)

        indices = np.where(vals > avg_ex)
        i0 = indices[0]
        irest = indices[1:]
        ii = tuple(x for x in itertools.chain([i0],irest))
        jj = tuple(x for x in itertools.chain([np.zeros_like(i0)],irest))
        term[ii] = (dmax[ii] - avg_ex[jj])**2

        indices = np.where(vals <= avg_ex)
        i0 = indices[0]
        irest = indices[1:]
        ii = tuple(x for x in itertools.chain([i0],irest))
        jj = tuple(x for x in itertools.chain([np.zeros_like(i0)],irest))
        term[ii] = (avg_ex[jj] - dmin[ii])**2

        dsum = ma.sum(term*wt,axis=0)     # Sum for weighted average of deviations

        dev = 0.5*np.sqrt(dsum/(norm*neff))

        if isinstance(avg,(float,np.float)):
            avg = avg_ex

        tmp_min = avg - dev
        ii = np.where(tmp_min < 0)
        tmp_min[ii] = TOL*avg[ii]

        return UncertContainer(avg,tmp_min,avg+dev)

max_stars_repo_name: ajoshpratt/westpa

max_issues_repo_name: astatide/westpa

max_forks_repo_name: ajoshpratt/westpa


text:
import numpy as np
import json
import math

def fov2length(angle):
    return math.tan(math.radians(angle) / 2) * 2

def _convert_camera_params(input_camera_params,
                           view_res):
    """
    Check and convert camera parameters in config file to pixel-space

    :param cam_params: { ["fx", "fy" | "fov"], "cx", "cy", ["normalized"] },
        the parameters of camera
    :return: camera parameters
    """
    input_is_normalized = bool(input_camera_params.get('normalized'))
    camera_params = {}
    if 'fov' in input_camera_params:
        if input_is_normalized:
            camera_params['fy'] = 1 / fov2length(input_camera_params['fov'])
            camera_params['fx'] = camera_params['fy'] / view_res["y"] * view_res["x"]
        else:
            camera_params['fx'] = camera_params['fy'] = view_res["x"] / \
                                                        fov2length(input_camera_params['fov'])
        camera_params['fy'] *= -1
    else:
        camera_params['fx'] = input_camera_params['fx']
        camera_params['fy'] = input_camera_params['fy']
    camera_params['cx'] = input_camera_params['cx']
    camera_params['cy'] = input_camera_params['cy']
    if input_is_normalized:
        camera_params['fx'] *= view_res["y"]
        camera_params['fy'] *= view_res["x"]
        camera_params['cx'] *= view_res["y"]
        camera_params['cy'] *= view_res["x"]
    return camera_params

def load_pose_json(filename,near,far):
    info = None
    with open(filename, encoding='utf-8') as file:
        info = json.loads(file.read())
    poses = []
    for i in range(len(info["view_centers"])):
        pose = np.zeros((3, 5), dtype=np.float32)
        Rt = np.eye(4, dtype=np.float32)
        if len(info["view_rots"][i]) == 9:
            for j in range(9):
                Rt[j // 3, j % 3] = info["view_rots"][i][j]
        #     right up backward -> down right backward
        #     temp = Rt[:,0]
        #     Rt[:,0] = Rt[:,1]
        #     Rt[:,1] = temp
        #     Rt[:, 0] *= -1
        elif len(info["view_rots"][i]) == 2:
            Rt[1, 1] = math.cos(math.radians(info["view_rots"][i][0]))
            Rt[2, 1] = math.sin(math.radians(info["view_rots"][i][0]))
            Rt[1, 2] = -Rt[2,1]
            Rt[2, 2] = Rt[1,1]
            #     flip rots row3 and col3
            Rt[2] *= -1
            Rt[:, 2] *= -1
        for j in range(3):
            Rt[j, 3] = info["view_centers"][i][j]
        # Rt = np.linalg.inv(Rt)
        Rt = Rt[:3,:4]
        Rt = np.concatenate([-Rt[:, 1:2], Rt[:, 0:1], -Rt[:, 2:3], Rt[:,3:4]],1)

        # OGL to DX coordinate
        if len(info["view_rots"][i]) == 2:
            #     flip fy
            # cam_par["fy"] *= -1
            #     flip view center z
            Rt[2, 3] *= -1

        pose[0,:4] = Rt[0]
        pose[1,:4] = Rt[2]
        pose[2,:4] = Rt[1]

        pose[0, 4] = info["view_res"]["x"]
        pose[1, 4] = info["view_res"]["y"]
        # print(info["view_res"])
        if "fx" in info["cam_params"]:
            pose[2, 4] = info["cam_params"]["fx"]
        else:
            cam_par = _convert_camera_params(info["cam_params"],info["view_res"])
            # print("cam_par", cam_par)
            pose[2,4] = cam_par["fx"]
            pose[2,4] = 560.166104085872
            print("fx", cam_par["fx"])
        # OGL to DX coordinate
        if len(info["view_rots"][i]) == 2:
            #     flip fy
            cam_par["fy"] *= -1
        pose = pose.flatten()
        pose = np.concatenate((pose, [near,far]))
        ## better to know the near and far plane from Unity-DengNianchen.
        ## Hard to understand the exact meaning (the depth bound of scene geometry as I can see)
        # print(pose.shape)
        poses.append(pose)
    poses = np.stack(poses)

    print(poses.shape)
    return poses

import argparse
parser = argparse.ArgumentParser()
parser.add_argument('train_json', type=str,
                    help='input training set info')
parser.add_argument('near', type=float,
                    help='scene depth range near')
parser.add_argument('far', type=float,
                    help='scene depth range far')
args = parser.parse_args()

if __name__ == "__main__":
    jsonfile = args.train_json
    near = args.near
    far = args.far
    poses = load_pose_json(jsonfile,near,far)
    np.save(args.train_json+"poses_bounds.npy", poses)

max_stars_repo_name: snowymo/nerf

max_issues_repo_name: snowymo/nerf

max_forks_repo_name: snowymo/nerf


text:
import numpy as np
import matplotlib.pyplot as plt
import imports.arrival_process_simulation as arrival
import imports.stochastic_plots as plots
import imports.BM_simulators as bm_sim

from scipy.special import factorial, i0
from scipy import stats
from collections import Counter
from sklearn.neighbors import KernelDensity

from matplotlib.animation import FuncAnimation
from celluloid import Camera

COLOR1 = '#1f77b4'
COLOR2 = '#ff7f0e'

# ---------------------------------------- EXERCISE 1 ----------------------------------------

def exercise_1(t=2, max_n=40, lamb=10, n_samples=10**4):
    """
        Implements the whole exercise 1, in which we compare the empirical
        distribution obtained by simulating a Poisson process and the
        theoretical distribution.

        Parameters
        ----------
        t : float
            Current time at the simulation
        max_n : int
            Maximum number of
        lamb : float
            lambda parameter of the Poisson process
        n_samples : int
            Number of samples to generate for the empirical distribution
        Returns
        -------
        No returns

        Example
        -------
        >>> exercise_1()
    """
    # Theoretical
    ns = np.arange(max_n+1) * 1.0

    y_theoretical = (lamb*t)**ns * np.exp(-lamb*t) / factorial(ns)

    # Simulation
    counter = Counter([ len(arrival_times) for arrival_times in arrival.simulate_poisson(0, t, lamb, n_samples)])
    total_count = sum(counter.values())
    y_simulated = [counter[i]/total_count for i in ns]

    # Plotting
    width = 0.35
    plt.figure(figsize=(10, 6))
    plt.bar(ns - width/2, y_theoretical, width, label='Theoretical')
    plt.bar(ns + width/2, y_simulated, width, label='Simulated')

    plt.xlabel('n')
    plt.ylabel('P[ N(2) = n ]')
    plt.title('Empirical comparisson of a Poisson process with λ=10, t=2')
    plt.legend()

# ---------------------------------------- EXERCISE 2 ----------------------------------------

def plot_kde_in_axis(axis, ts, pdf, pdf_kde, samples, lamb, n):
    """
        Function to plot kernel estimation of the n-th arrival of a Poisson
        process andtheoretical density, along with generated samples
        for the estimation.

        Parameters
        ----------
        axis : matplotlib.axes
            Matplotlib axes of a figure to plot the charts
        ts : numpy.ndarray
            Array of times uses for the plot
        pdf: function
            Theoretical pdf values of the n-th arrival of a Poisson process
        pdf_kde: function
            Estimated pdf values of the n-th arrival of a Poisson process
        samples : numpy.ndarray
            Array with the generated arrival times
        lamb : float
            Lambda parameter of the Poisson process
        n : int
            Number of arrivals estimated

        Returns
        -------
        No returns, it fills the axis

        Example
        -------
        See function exercise_2() for an example.
    """
    # Plot pdf and estimation
    axis.fill_between(ts, pdf, alpha=0.3, color='C0', label='Theoretical')
    axis.plot(ts, pdf_kde, color='C1', label='Kernel density estimation')

    # Plot little x's near the X axis
    axis.plot(samples[:, 0], -0.06 - 0.05 * np.random.random(samples.shape[0]),
              "+k", label='Samples')

    axis.set_xlabel('t')
    axis.set_ylabel('P[ N(t) = {} ]'.format(n))
    axis.set_title('Empirical comparisson of a Poisson process with λ={}, n={}'.format(lamb, n))
    axis.legend()

def exercise_2(ns=[1,2,5,10], lamb=5, max_t=7, n_samples_kde=10**4, kde_bandwidth=0.1):
    """
        Implements the whole exercise 2, in which we try to find
        the differences between the kernel density estimation and the
        theoretical. Function to plot kernel estimation of the
        n-th arrival of a Poisson process and theoretical density,
        along with generated samples for the estimation.

        Parameters
        ----------
        axis : matplotlib.axes
            Matplotlib axes of a figure to plot the charts
        ts : numpy.ndarray
            Array of the times uses for the plot
        pdf: function
            Theoretical pdf of the n-th arrival of a Poisson process
        pdf_kde: function
            Estimated pdf of the n-th arrival of a Poisson process
        samples : numpy.ndarray
            Array with the generated arrival times
        lamb : float
            Lambda parameter of the Poisson process
        n : int
            Number of arrivals estimated

        Returns
        -------
        No returns, it fills the axis

        Example
        -------
        >>> exercise_2()
    """
    ts = np.arange(max_t, step=0.01)
    _, axis = plt.subplots(2, 2, figsize=(18, 12))

    for ax, n in zip(axis.flatten(), ns):
        # Theoretical
        sn_pdf = stats.erlang.pdf(ts, a=n, scale=1/lamb)

        # Kernel Density Estimation
        sn_empirical = [ arrival_times[n-1] for arrival_times in arrival.simulate_poisson(0, max_t, lamb, n_samples_kde)]
        sn_empirical = np.array(sn_empirical).reshape(-1, 1)
        kernel_density_estimatior = KernelDensity(kernel='gaussian', bandwidth=kde_bandwidth).fit(sn_empirical)
        sn_pdf_estimation = np.exp(kernel_density_estimatior.score_samples(np.array(ts).reshape(-1, 1)))

        # Plotting
        plot_kde_in_axis(ax, ts, sn_pdf, sn_pdf_estimation, sn_empirical, lamb, n)

# ---------------------------------------- EXERCISE 4 ----------------------------------------

def simulate_team_scores(t=90, lamb1=0.02, lamb2=0.03, n_samples=10**5):
    """

    Simulates n_samples games between the two teams.

    Parameters
    ----------
    t : int, optional
        Array of the times uses for the plot
    lamb1 : float, optional
        Lambda parameter of the first Poisson process
    lamb2 : float, optional
        Lambda parameter of the second Poisson process
    n_samples : int, optional
        Number of games estimated estimated

    Returns
    -------
    numpy.ndarray
        2d-array where each row is a game with two values, the scores for
        teams A and B respectively

    Example
    -------
    >>> simuations = simulate_team_scores()
    """
    return [
        [ len(arrival_times1), len(arrival_times2) ]
        for arrival_times1, arrival_times2
        in zip(arrival.simulate_poisson(0, t, lamb1, n_samples),
               arrival.simulate_poisson(0, t, lamb2, n_samples))
    ]

def estimate_prob(condition, t=90, lamb1=0.02, lamb2=0.03, n_samples=10**5, seed=123):
    """
    Estimates the probability of the given condition being fullfilled
    in a random match between the two teams. It simulates n_samples games
    and return the proportion on games that passed the condition.

    Parameters
    ----------
    condition : lambda function
        Condition to apply to the games, recieves 2 values (the scores)
        and returns a boolean whether the condition was fullfilled.
    t : int, optional
        Array of the times uses for the plot
    lamb1 : float, optional
        Lambda parameter of the first Poisson process
    lamb2 : float, optional
        Lambda parameter of the second Poisson process
    n_samples : int, optional
        Number of games estimated estimated
    seed : int, optional
        seed to fix on numpy.random

    Returns
    -------
    float
        An estimation of the probability of the given condition being
        fullfilled.

    Example
    -------
    >>> tie_on_0 = lambda x: x[0] == 0 and x[1] == 0
    >>> estimate_prob(tie_on_0)
    """
    np.random.seed(seed)
    scores = simulate_team_scores(t=t, lamb1=lamb1, lamb2=lamb2, n_samples=n_samples)
    fullfil_condition = np.sum([ condition(score) for score in scores ])
    return 1.0*fullfil_condition / len(scores)

def team_B_scores_first_prob(t=90, lamb1=0.02, lamb2=0.03, n_samples=10**5):
    """
    Estimates the probability of team B scoring first by simulating
    n_samples games and returning the proportion of games that satisfy
    the condition.

    Parameters
    ----------
    t : int, optional
        Array of the times uses for the plot
    lamb1 : float, optional
        Lambda parameter of the first Poisson process
    lamb2 : float, optional
        Lambda parameter of the second Poisson process
    n_samples : int, optional
        Number of games estimated estimated

    Returns
    -------
    float
        An estimation of the probability of the given condition being
        fullfilled.

    Example
    -------
    >>> team_B_scores_first_prob()
    """
    count = 0
    for arrival_times1, arrival_times2 \
        in zip(arrival.simulate_poisson(0, t, lamb1, n_samples),
               arrival.simulate_poisson(0, t, lamb2, n_samples)):
        if len(arrival_times2) == 0:
            # Team B didn't score a single goal.
            continue
        elif len(arrival_times1) == 0:
            # Team A didn't score a single goal but B did.
            count += 1
        else:
            # Both teams scored, compare when.
            count += arrival_times1[0] > arrival_times2[0]
    return 1.0 * count / n_samples

# ---------------------------------------- EXERCISE 6 ----------------------------------------

def simulate_wiener_process(t0=0, t1=1, delta_t=0.001, n_processes=100):
    """
    Simulates n_processes wiener processes.

    Parameters
    ----------
    t0 : float, optional
        Initial time
    t1 : float, optional
        Final time
    delta_t : float, optional
        Time increment
    n_processes : int, optional
        Number of processes to be simulated

    Returns
    -------
    numpy.ndarray
        The array of times when the simulations took place.
    numpy.ndarray
        A 2d-array containing the simulations. Each row corresponds
        to a simulated processed, and has the same length and the
        array of times also returned.

    Example
    -------
    >>> ts, trajectories = simulate_wiener_process(n_processes=1000)
    """
    n_steps = int((t1-t0) / delta_t)
    return bm_sim.simulate_arithmetic_BM(t0, B0=0, T=t1-t0, mu=0, sigma=1,
                                         M=n_processes, N=n_steps)

def subplot_mean_and_std(axis, x, mean, std, color=COLOR1,
                         xlims=None, ylims=None, xlabel=None,
                         ylabel=None, title=None, alpha_std=.3):
    """
    Plots the passed mean and std.

    Parameters
    ----------
    axis : matplotlib.pyplot.axis
        axis to plot on
    x : numpy.ndarray
        x-component to plot
    mean : numpy.ndarray
        mean of the y-component to plot
    std : numpy.ndarray
        std of the y-component to plot
    color : string, optional
        Color to plot the mean and std on
    xlims : numpy.ndarray, optional
        xlims for the plot
    ylims : numpy.ndarray, optional
        xlims for the plot
    xlabel : string, optional
        xlabel for the plot
    ylabel : string, optional
        ylabel for the plot
    title : string, optional
        Title for the plot
    alpha_std : float, optional
        Alpha of the std-filling color

    Returns
    -------
    No returns, it fills the axis

    Example
    -------
    >>> simulate_wiener_process(n_processes=1000)
    >>> mean, std = np.mean(trajectories, axis=0), np.std(trajectories, axis=0)
    >>> fig, axis = plt.figure(figsize=(12, 8))
    >>> subplot_mean_and_std(axis, ts, mean, 2*std)
    """
    axis.plot(x, mean, color='b')
    axis.fill_between(x, mean-std, mean+std, color=color, alpha=alpha_std)
    if xlims is not None: axis.set_xlim(xlims)
    if ylims is not None: axis.set_ylim(ylims)
    if xlabel is not None: axis.set_xlabel(xlabel)
    if ylabel is not None: axis.set_ylabel(ylabel)
    if title is not None: axis.set_title(title)

def plot_trajectories(ts, trajectories, axis=None, max_trajectories=50):
    """
    Plots at most max_trajectories of the given trajectories,
    along with the mean and std of all the trajectories.

    Parameters
    ----------
    ts : numpy.ndarray
        Vector of times where the trajectories happened
    trajectories : numpy.ndarray
        2d-array containing all the trajectories
    axis : matplotlib.pyplot.axis, optional
        Axis to plot on
    max_trajectories : int, optional
        Max number of trajectories to plot

    Returns
    -------
    No returns, it fills the axis

    Example
    -------
    >>> ts, trajectories = simulate_wiener_process(n_processes=1000)
    >>> plot_trajectories(axis, ts)
    """
    mean, std = np.mean(trajectories, axis=0), np.std(trajectories, axis=0)

    if axis is None:
        plt.figure(figsize=(12, 8))
        axis = plt.gca()

    for t in trajectories[0:max_trajectories]:
        axis.plot(ts, t, linewidth=0.7, label='_nolegend_')
    subplot_mean_and_std(axis, ts, mean, 2*std)

    axis.legend(['Mean trajectory', '$\pm$ 2 * standard deviation'])
    axis.set_xlabel("t")
    axis.set_ylabel("W(t)")

def estimate_covariance(fixed_t=0.25, t0=0, t1=1, delta_t=0.001, n_processes=100):
    """
    Estimates the covariance of a wiener process between the fixed time
    and times from t0 to t1.

    Parameters
    ----------
    fixed_t : float, optional
        Fixed time to compute the covariance against
    t0 : float, optional
        Initial time
    t1 : float, optional
        Final time
    delta_t : float, optional
        Time increment
    n_processes : int, optional
        Number of trajectories to simulate

    Returns
    -------
    numpy.ndarray
        The array of times when the simulations took place.
    numpy.ndarray
        An array containing in position i Cov( W(t_i), W(fixed_t) )

    Example
    -------
    >>> estimate_covariance()
    """
    ts, trajectories = simulate_wiener_process(
        t0=t0, t1=t1, delta_t=delta_t, n_processes=n_processes)

    fixed_index = np.where( np.abs(ts-fixed_t) < 1e-10 )[0][0]
    fixed_time_samples = trajectories.T[fixed_index]
    fixed_time_samples -= np.mean(fixed_time_samples)
    cov = [ np.mean( fixed_time_samples * (time_sample - np.mean(time_sample)))
            for time_sample in trajectories.T]

    return ts, cov

def plot_estimated_covariance(fixed_t=0.25, t0=0, t1=1,
                              delta_t=0.001, n_processes=100, n_estimations=100):
    """
    Plots a estimation (mean and std) of the covariance of a wiener process.
    For time t in [t0, t1], plots the estimated covariance between W(t) and
    W(fixed_t).

    Parameters
    ----------
    fixed_t : float, optional
        Fixed time to compute the covariance against
    t0 : float, optional
        Initial time
    t1 : float, optional
        Final time
    delta_t : float, optional
        Time increment
    n_processes : int, optional
        Number of trajectories to simulate
    n_estimations : int, optional
        Number of times the n_processes are simulated for
        the estimation

    Returns
    -------
    No returns, it fills the axis

    Example
    -------
    >>> plot_estimated_covariance()
    """
    t, _ = estimate_covariance(fixed_t=fixed_t, t0=t0, t1=t1, delta_t=delta_t, n_processes=1)
    cov_estimations = [
        estimate_covariance(fixed_t=fixed_t,
                            t0=t0, t1=t1, delta_t=delta_t, n_processes=n_processes)[1]
        for _ in range(n_estimations)
    ]
    cov_mean, cov_std = np.mean(cov_estimations, axis=0), np.std(cov_estimations, axis=0)

    # Plotting
    plt.figure(figsize=(12, 8))
    plt.plot([0, fixed_t, 1], [0, fixed_t, fixed_t], color='red')
    subplot_mean_and_std(plt.gca(), t, cov_mean, cov_std, xlims=[0,1],
                         xlabel='t', ylabel='Cov[ W(t), W({}) ]'.format(fixed_t))
    plt.legend(['Theoretical covariance', 'Mean empirical covariance',
                '$\pm$ standard deviation'], loc='lower right')

# ---------------------------------------- EXERCISE 7 ----------------------------------------

def plot_hist_and_pdf(X, pdf, axis=None, max_bins=50,
                      xlims=None, ylims=None):
    """
    Plots the histogam and the theorical pdf of the given data
    in the given axis.

    Parameters
    ----------
    X : numpy.ndarray
        An array containing the data to histogram
    pdf : function
        pdf function to evaluate
    axis : matplotlib.pyplot.axis, optional
        Axis to plot on
    max_bins : int, optional
        Max number of bins to plot
    xlims : numpy.ndarray, optional
        xlims for the plot
    ylims : numpy.ndarray, optional
        xlims for the plot

    Returns
    -------
    No returns, it fills the axis

    Example
    -------
    See function plot_trajectories_and_hist() for an example
    """
    if axis is None:
        plt.figure(figsize=(12, 8))
        axis = plt.gca()

    # Plot histogram
    n_bins = np.min((np.int(np.round(np.sqrt(len(X)))), max_bins))
    axis.hist(X, bins=n_bins, density=True, color=COLOR1, alpha=0.3)
    axis.set_xlabel('x')
    axis.set_ylabel('pdf(x)')
    if xlims is not None:
        axis.set_xlim(xlims)
    if ylims is not None:
        axis.set_ylim(ylims)

    # Compare with exact distribution
    n_plot = 1000
    X_from, X_to = np.min(X), np.max(X)
    if xlims is not None:
        X_from, X_to = xlims

    x_plot = np.linspace(X_from, X_to, n_plot)
    y_plot = pdf(x_plot)

    axis.plot(x_plot, y_plot, linewidth=2, color=COLOR2)
    axis.legend(['Theoretical distribution', 'Empirical histogram'])

def plot_trajectories_and_hist(ts, trajectories, fixed_t_index, fig=None, axis=None, pdf=None,
                                B0=0, BM_mu=0, BM_sigma=1, hist_xlims=None, hist_ylims=None):
    """
    Graphs a 2-plots subplot with the trajectories and the histogram
    at a given fixed_t_index time

    Parameters
    ----------
    ts : numpy.ndarray
        Vector of times where the trajectories happened
    trajectories : numpy.ndarray
        2d-array containing all the trajectories
    fixed_t_index : int
        Index refering to the fixed time whose histogram
        will be plotted.
    fig : matplotlib.figure.Figure, optional
        Figure to plot on
    axis : matplotlib.pyplot.axis, optional
        Axis to plot on
    pdf : function
        pdf function to evaluate for the histogram
    B0 : float, optional
        B0 value for the brownian motion, in case we want to
        represent that
    BM_mu : float, optional
        mu value for the brownian motion, in case we want to
        represent that
    BM_sigma : float, optional
        sigma value for the brownian motion, in case we want to
        represent that
    hist_xlims : numpy.ndarray, optional
        xlims for the plot
    hist_ylims : numpy.ndarray, optional
        xlims for the plot

    Returns
    -------
    No returns, it fills the axis

    Example
    ------
    >>> ts, trajectories = ex.simulate_wiener_process(n_processes=1000)
    >>> plot_trajectories_and_hist(ts, trajectories, 25)
    """
    if fig is None or axis is None:
        fig, axis = plt.subplots(1, 2, figsize=(20, 6))
    fixed_time = ts[fixed_t_index]

    # Plot the trajectories
    plot_trajectories(ts, trajectories, axis=axis[0])
    axis[0].axvline(x=fixed_time, color=COLOR2)
    axis[0].legend(['Mean trajectory', 't = {:.2f}'.format(fixed_time),
                    '$\pm$ 2 * standard deviation'])

    # Define the pdf for the Wiener and Brownian motions
    if pdf is None:
        pdf_mu, pdf_std = B0 + BM_mu*(fixed_time-ts[0]), BM_sigma*np.sqrt(fixed_time-ts[0])
        if pdf_std == 0:
            # Dirac's delta centered on mu
            pdf = lambda x: x == pdf_mu
        else:
            pdf = lambda x: stats.norm(pdf_mu, scale=pdf_std).pdf(x)

    # Plot the histogram
    plot_hist_and_pdf(trajectories[:, fixed_t_index], pdf,
                        axis=axis[1], xlims=hist_xlims, ylims=hist_ylims)
    fig.suptitle('Trajectories and histogram for t={:.2f}'.format(fixed_time))

# ---------------------------------------- EXERCISE 8 ----------------------------------------

def plot_trajectories_animation(ts, trajectories, n_frames=10,
                                hist_xlims=[-2,2], hist_ylims=[0,2], gif_name='',
                                B0=0, BM_mu=0, BM_sigma=1, save_gif=False):
    """
    Creates and saves an animation of n_frames with the given trayectories.
    It plots using plot_trajectories_and_hist at different times.

    Drawback of using celluloid library:
    - Lims must be fixed for the whole animation, the will be the lims of the last
        figure plotted.
    - Legends are acummulative. This issue, summed up with the matplotlib.fill_between
        legend beeing always the last legend to be added means the legend of
        matplotlib.fill_between is broken in this animation.

    Parameters
    ----------
    ts : numpy.ndarray
        Vector of times where the trajectories happened
    trajectories : numpy.ndarray
        2d-array containing all the trajectories
    n_frames : int, optional
        Number of frames in the animation
    hist_pdf : function
        pdf function to evaluate for the histogram
    hist_xlims : numpy.ndarray, optional
        xlims for the plot
    hist_ylims : numpy.ndarray, optional
        xlims for the plot
    gif_name : string, optional
        Name for the gif to be save
    B0 : float, optional
        B0 value for the brownian motion, in case we want to
        represent that
    BM_mu : float, optional
        mu value for the brownian motion, in case we want to
        represent that
    BM_sigma : float, optional
        sigma value for the brownian motion, in case we want to
    save_gif : bool, optional
        A boolean checking whether to save the animation or not

    Returns
    -------
    No returns, it fills the axis

    Example
    -------
    >>> ts, trajectories = ex.simulate_wiener_process(n_processes=1000)
    >>> plot_trajectories_animation(ts, trajectories, n_frames=20, gif_name='wiener')
    """
    plt.rc('figure', figsize=(20, 6))
    fig, ax = plt.subplots(1, 2)
    camera = Camera(fig)
    time_indexes = np.arange(0, len(ts), int(len(ts)/n_frames))

    for t_index in time_indexes:
        plot_trajectories_and_hist(ts, trajectories, fixed_t_index=t_index,
                                   fig=fig, axis=ax,
                                   B0=B0, BM_mu=BM_mu, BM_sigma=BM_sigma,
                                   hist_xlims=hist_xlims, hist_ylims=hist_ylims)
        camera.snap()
    fig.suptitle('Trajectories and histogram for t $\in$[0, {:.1f}]'.format(ts[-1]))

    animation = camera.animate()
    if save_gif:
        animation.save('{}.gif'.format(gif_name), writer='pillow')
    plt.close()

max_stars_repo_name: Ocete/mcd-pe-continuo

max_issues_repo_name: Ocete/mcd-pe-continuo

max_forks_repo_name: Ocete/mcd-pe-continuo

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
""" The IMED policy of [Honda & Takemura, JMLR 2015].

- Reference: [["Non-asymptotic analysis of a new bandit algorithm for semi-bounded rewards", J. Honda and A. Takemura, JMLR, 2015](http://jmlr.csail.mit.edu/papers/volume16/honda15a/honda15a.pdf)].
"""
from __future__ import division, print_function  # Python 2 compatibility

__author__ = "Lilian Besson"
__version__ = "0.9"

import numpy as np
np.seterr(divide='ignore')  # XXX dangerous in general, controlled here!

try:
    from .DMED import DMED
    from .usenumba import jit
    from .kullback import klBern
except ImportError:
    from DMED import DMED
    from usenumba import jit
    from kullback import klBern

# --- Utilitary functions Dinf

from scipy.optimize import minimize_scalar
from warnings import filterwarnings
filterwarnings("ignore", message="Method 'bounded' does not support relative tolerance in x; defaulting to absolute tolerance", category=RuntimeWarning)

def Dinf(x=None, mu=None, kl=klBern,
        lowerbound=0, upperbound=1,
        precision=1e-6, max_iterations=50
    ):
    r""" The generic Dinf index computation.

    - ``x``: value of the cum reward,
    - ``mu``: upperbound on the mean ``y``,
    - ``kl``: the KL divergence to be used (:func:`klBern`, :func:`klGauss`, etc),
    - ``lowerbound``, ``upperbound=1``: the known bound of the values ``y`` and ``x``,
    - ``precision=1e-6``: the threshold from where to stop the research,
    - ``max_iterations``: max number of iterations of the loop (safer to bound it to reduce time complexity).

    .. math::

        D_{\inf}(x, d) \simeq \inf_{\max(\mu, \mathrm{lowerbound}) \leq y \leq \mathrm{upperbound}} \mathrm{kl}(x, y).

    .. note:: It uses a call the :func:`scipy.optimize.minimize_scalar`. If this fails, it uses a **bisection search**, and one call to ``kl`` for each step of the bisection search.
    """
    # lower and upper bounds
    l = max(lowerbound, mu)
    u = upperbound
    def f(i): return kl(x, i)  # objective function!
    res = minimize_scalar(f, bounds=[l, u],
        method="bounded",
        tol=precision,
        options={"maxiter": max_iterations, "disp": False}
    )
    if hasattr(res, "x"):
        return res.x
    else:
        print("Warning: the call to scipy.optimize.minimize_scalar failed, using hand-written bisection search instead...")  # DEBUG
    # start in the middle
    value = (l + u) / 2.
    current_kl = kl(x, value)
    _count_iteration = 0
    while _count_iteration < max_iterations and u - value > precision:
        _count_iteration += 1
        # try to see right
        value = (l + u) / 2.
        new_kl = kl(x, value)
        if new_kl > current_kl:  # need to go left
            u = value
        else:  # need to go right
            current_kl = new_kl
            l = value
    return current_kl

# --- IMED

class IMED(DMED):
    """ The IMED policy of [Honda & Takemura, JMLR 2015].

    - Reference: [["Non-asymptotic analysis of a new bandit algorithm for semi-bounded rewards", J. Honda and A. Takemura, JMLR, 2015](http://jmlr.csail.mit.edu/papers/volume16/honda15a/honda15a.pdf)].
    """

    def __init__(self, nbArms, tolerance=1e-4, kl=klBern, lower=0., amplitude=1.):
        super(IMED, self).__init__(nbArms, tolerance=tolerance, kl=kl, lower=lower, amplitude=amplitude)

    def __str__(self):
        return r"IMED({})".format(self.kl.__name__[2:])

    def one_Dinf(self, x, mu):
        r""" Compute the :math:`D_{\inf}` solution, for one value of ``x``, and one value for ``mu``."""
        return Dinf(x=x, mu=mu, kl=self.kl, lowerbound=self.lower, upperbound=self.lower + self.amplitude, precision=self.tolerance)

    # XXX Use this hack to vectorize one_Dinf ?
    # Dinf = np.vectorize(one_Dinf) ?
    def Dinf(self, xs, mu):
        r""" Compute the :math:`D_{\inf}` solution, for a vector of value of ``xs``, and one value for ``mu``."""
        return np.array([ self.one_Dinf(x, mu) for x in xs ])

    def choice(self):
        r""" Choose an arm with **minimal** index (uniformly at random):

        .. math:: A(t) \sim U(\arg\min_{1 \leq k \leq K} I_k(t)).

        Where the indexes are:

        .. math:: I_k(t) = N_k(t) D_{\inf}(\hat{\mu_{k}}(t), \max_{k'} \hat{\mu_{k'}}(t)) + \log(N_k(t)).
        """
        empiricalMeans = self.rewards / self.pulls
        bestEmpiricalMean = np.max(empiricalMeans)
        values_Dinf = self.Dinf(empiricalMeans, bestEmpiricalMean)

        # now compute the indexes
        indexes_to_minimize = self.pulls * values_Dinf + np.log(self.pulls)
        indexes_to_minimize[self.pulls < 1] = float('-inf')

        # then do as IndexPolicy but with a min instead
        # try:
        return np.random.choice(np.nonzero(indexes_to_minimize == np.min(indexes_to_minimize))[0])
        # except ValueError:
        #     if not np.all(np.isnan(indexes_to_minimize)):
        #         raise ValueError("Error: unknown error in IMED.choice(): the indexes were {} but couldn't be used to select an arm.".format(indexes_to_minimize))
        #     return np.random.randint(self.nbArms)


text:
import numpy as np
from cobra import Model, Reaction, Metabolite
#from optlang import Variable, Constraint, Objective, Model
from cobra.sampling import sample
from scipy.sparse import diags

EPSILON = 1e-7

def sample_initial_concentrations(kmodel,
                                  reference_concentrations,
                                  lower_bound=0.8,
                                  upper_bound=1.2,
                                  n_samples=10,
                                  absolute_bounds = False):

    concentrations = np.array([reference_concentrations[k] for k in kmodel.variables])

    if kmodel.conservation_relation is None:
        N = len(kmodel.variables)
        rand = np.random.uniform(low=lower_bound,
                                 high=upper_bound,
                                 size=(N,n_samples))

    #
    else:
        # Todo Also allow
        lower_bound = {k: v * lower_bound for k,v in  reference_concentrations.items()}
        upper_bound = {k: v * upper_bound for k, v in reference_concentrations.items()}

        # Account for volume differences in compartments
        if kmodel.volume_ratio_func is None:
            effective_conservation_relations = kmodel.conservation_relation.todense()
        else:
            param_values = {p.symbol: p.value for p in kmodel.parameters.values()}
            volume_ratios = kmodel.volume_ratio_func(param_values)
            inv_volume_ratio_matrix =  diags(1./np.array(volume_ratios)).todense()
            effective_conservation_relations = kmodel.conservation_relation.todense()\
                                               .dot(inv_volume_ratio_matrix)

        # Create a linear problem and sample it
        rhs = effective_conservation_relations.dot(concentrations)
        linmodel = create_linear_model(effective_conservation_relations,
                                       rhs,
                                       kmodel.reactants,
                                       lower_bound=lower_bound,
                                       upper_bound=upper_bound)

        return sample(linmodel, n_samples, )

def create_linear_model(A, rhs, variables, lower_bound=None, upper_bound=None):
    # This is a bit retarded but this way we can use the sampling functions from cobra

    model = Model('lin_model')
    # TODO have cases for lower and upper bound  None
    rxns = np.array([ Reaction(v, lower_bound=lower_bound[v],
                                  upper_bound=upper_bound[v]) for v in variables])

    model.add_reactions(rxns)

    lin_variables = np.array([r.forward_variable for r in rxns])
    NR, NC = A.shape
    constraints = [model.problem.Constraint( (A[i,:].dot(lin_variables))[0,0] ,
                                             lb=np.float(rhs[:,i]),
                                             ub=np.float(rhs[:,i])) for i in range(NR)]
    model.add_cons_vars(constraints)

    return model


text:
import numpy as np
import pandas as pd
import matplotlib.pylab as plt
import gradio as gr
import PIL.Image as Image
import tensorflow as tf
import tensorflow_hub as hub

TF_MODEL_URL = 'https://tfhub.dev/google/on_device_vision/classifier/landmarks_classifier_asia_V1/1'
LABEL_MAP_URL = 'https://www.gstatic.com/aihub/tfhub/labelmaps/landmarks_classifier_asia_V1_label_map.csv'
IMAGE_SHAPE = (321, 321)

classifier = tf.keras.Sequential([hub.KerasLayer(TF_MODEL_URL,
                                                 input_shape=IMAGE_SHAPE+(3,),
                                                 output_key="predictions:logits")])

df = pd.read_csv(LABEL_MAP_URL)
label_map = dict(zip(df.id, df.name))
img_loc = "Image.jpeg"
img = Image.open(img_loc).resize(IMAGE_SHAPE)
img
img = np.array(img)/255.0
img.shape

img = img[np.newaxis, ...]
img.shape
result = classifier.predict(img)
result
label_map[np.argmax(result)]

class_names=list(label_map.values())

def classify_image(image):
    img = np.array(image)/255.0
    img = img[np.newaxis, ...]
    prediction = classifier.predict(img)
    return label_map[np.argmax(prediction)]

image = gr.inputs.Image(shape=(321, 321))
label = gr.outputs.Label(num_top_classes=1)

gr.Interface(
    classify_image,
    image,
    label,
    capture_session=True).launch(debug=True);


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Minimal test script to check for modules needed in python workshop.

Execute this code at the command line by typing:

python workshop_checklist.py

How to get a command line:

- On OSX run this with the Terminal application.

- On Windows, go to the Start menu, select 'Run' and type 'cmd' (without the
quotes) to run the cmd.exe Windows Command Prompt.

Run the script and follow the instructions it prints at the end.
"""

##############################################################################
# Config here

# Modules whose imports we're going to validate
MODULES = ['IPython',
           'numpy',
           'scipy',
           'scipy.io',
           'matplotlib','pylab',
           # 'sympy', 'Cython', 'networkx', 'mayavi.mlab',
           # 'setuptools',
           ]

##############################################################################
# Code begins

# Standard library imports
import glob
import os
import platform
import sys

from StringIO import StringIO

# Third-party imports
import nose
import nose.tools as nt

#-----------------------------------------------------------------------------
# Generic utility functions
def sys_info():
    """Summarize some info about the system"""

    print '=================='
    print 'System information'
    print '=================='
    print 'os.name      :',os.name
    try:
        print 'os.uname     :',os.uname()
    except:
        pass
    print 'platform     :',sys.platform
    print 'platform+    :',platform.platform()
    print 'prefix       :',sys.prefix
    print 'exec_prefix  :',sys.exec_prefix
    print 'executable   :',sys.executable
    print 'version_info :',sys.version_info
    print 'version      :',sys.version
    print '=================='

#-----------------------------------------------------------------------------
# Tests

def check_import(mname):
    "Check that the given name imports correctly"
    exec "import %s as m" % mname

    if mname == 'matplotlib':
        m.use('Agg')
        m.rcParams['figure.subplot.top']= 0.85

    try:
        vinfo = m.__version__
    except AttributeError:
        vinfo = '*no info*'

    print 'MOD: %s, version: %s' % (mname,vinfo)

# Test generators are best written without docstrings, because nose can then
# show the parameters being used.
def test_imports():
    for mname in MODULES:
        yield check_import, mname

# Test generator, don't put a docstring in it
def test_loadtxt():
    import numpy as np
    import numpy.testing as npt

    # Examples taken from the loadtxt docstring
    array = np.array

    c = StringIO("0 1\n2 3")
    a1 = np.loadtxt(c)
    a2 = np.array([[ 0.,  1.],
                   [ 2.,  3.]])
    yield npt.assert_array_equal,a1,a2

    d = StringIO("M 21 72\nF 35 58")
    a1 = np.loadtxt(d, dtype={'names': ('gender', 'age', 'weight'),
                         'formats': ('S1', 'i4', 'f4')})

    a2 = np.array([('M', 21, 72.0), ('F', 35, 58.0)],
                  dtype=[('gender', '|S1'), ('age', '<i4'), ('weight', '<f4')])
    yield npt.assert_array_equal,a1,a2

    c = StringIO("1,0,2\n3,0,4")
    x,y = np.loadtxt(c, delimiter=',', usecols=(0,2), unpack=True)
    yield npt.assert_array_equal,x,np.array([ 1.,  3.])
    yield npt.assert_array_equal,y,np.array([ 2.,  4.])

def test_plot():
    "Simple plot generation."
    from matplotlib import pyplot as plt
    plt.figure()
    plt.plot([1,2,3])
    plt.xlabel('some numbers')
    plt.savefig('tmp_test_plot.png')

def test_plot_math():
    "Plots with math"
    from matplotlib import pyplot as plt
    plt.figure()
    plt.plot([1,2,3],label='data')
    t=(r'And X is $\sum_{i=0}^\infty \gamma_i + \frac{\alpha^{i^2}}{\gamma}'
       r'+ \cos(2 \theta^2)$')
    plt.title(t)
    plt.legend()
    plt.grid()
    plt.savefig('tmp_test_plot_math.png')

def cleanup_pngs():
    """Remove temporary pngs made by our plotting tests"""

    for f in glob.glob('tmp_test_plot*.png'):
        try:
            os.remove(f)
        except OSError:
            print '*** Error: could not remove file',f

#-----------------------------------------------------------------------------
# Main routine, executed when this file is run as a script
#
if __name__ == '__main__':
    print "Running tests:"
    # This call form is ipython-friendly
    nose.runmodule(argv=[__file__,'-vvs'],
                   exit=False)
    print """
***************************************************************************
                           TESTS FINISHED
***************************************************************************

If the printout above did not finish in 'OK' but instead says 'FAILED', copy
and send the *entire* output (all info above and summary below) to the
instructor for help.
"""
    sys_info()


text:
from IMLearn.utils import split_train_test
from IMLearn.learners.regressors import LinearRegression

from typing import NoReturn
import numpy as np
import pandas as pd
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio

pio.templates.default = "simple_white"

def load_data(filename: str) -> (pd.DataFrame, pd.Series):
    """
    Load house prices dataset and preprocess data.
    Parameters
    ----------
    filename: str
        Path to house prices dataset

    Returns
    -------
    Design matrix and response vector (prices) - either as a single
    DataFrame or a Tuple[DataFrame, Series]
    """
    x = pd.DataFrame(pd.read_csv(filename, usecols=np.arange(2,17))
                     .drop_duplicates())
    x = x.where(x>=0, other=pd.NaT).dropna()
    x = x[x['zipcode'] > 0]
    x = pd.get_dummies(x, columns=['zipcode'])
    x.loc[x['yr_renovated'] > 0, 'yr_since_renovated'] = x['yr_renovated']
    x.loc[x['yr_renovated'] <= 0, 'yr_since_renovated'] = x['yr_built']
    del x['yr_built']
    del x['yr_renovated']
    x = x[x['yr_since_renovated'] > 0]
    x = x[x['sqft_living'] > 0]
    x = x[x['sqft_lot'] > 0]
    x = x.astype(float)

    x.to_csv(r"C:\Users\romyb\Documents\iml\IML.HUJI\exercises\temp_labels.csv")
    y = x['price']
    del x['price']
    return x,y

def feature_evaluation(X: pd.DataFrame, y: pd.Series, output_path: str = ".") \
        -> NoReturn:
    """
    Create scatter plot between each feature and the response.
        - Plot title specifies feature name
        - Plot title specifies Pearson Correlation between feature and response
        - Plot saved under given folder with file name including feature name
    Parameters
    ----------
    X : DataFrame of shape (n_samples, n_features)
        Design matrix of regression problem

    y : array-like of shape (n_samples, )
        Response vector to evaluate against

    output_path: str (default ".")
        Path to folder in which plots are saved
    """
    X.insert(0,'price', y)
    cov_XY = X.cov()['price']
    std_y = y.std()
    X.sort_values(by=['price'])

    for col in X.columns:
        corr = float(cov_XY[col]/(std_y*(X[col].std())))
        print(col)
        print(corr)
        # figi = go.Figure([go.Scatter(x=X['price'], y= X[col], mode='markers')],
        #           layout=go.Layout(title=f"feature is: {col} and the pearson "
        #                                  f"correlation between is: {corr}"
        #                            , xaxis_title = "price"
        #                            , yaxis_title = f"number of {col}", width= 1000))
        # pio.write_image(figi, output_path+'/'+ col + '.png')
    del X['price']

def fit_test_model(x_train: pd.DataFrame, y_train: pd.Series, x_test, y_test):
    # For every percentage p in 10%, ..., 100%, repeat the following 10 times:
    #   1) Sample p% of the overall training data
    #   2) Fit linear model (including intercept) over sampled set
    #   3) Test fitted model over test set
    #   4) Store average and variance of loss over test set
    # Then plot average loss as function of training size with error ribbon of

    sample_num = 10
    lini_model = LinearRegression(True)
    mean_loss, std_loss = np.zeros(91), np.zeros(91)
    x_test, y_test = np.array(x_test), np.array(y_test)
    a = np.arange(10,100,1)
    for p in a:
        cur_loss = np.zeros(sample_num)
        for j in range(sample_num):
            cur_x = x_train.sample(frac=p/100)
            cur_y = y_train.loc[cur_x.index]
            cur_x, cur_y = np.array(cur_x), np.array(cur_y)
            lini_model.fit(cur_x, cur_y)
            cur_loss[j] = lini_model.loss(x_test, y_test)

        mean_loss[p-10] = np.mean(cur_loss)
        std_loss[p-10] = np.std(cur_loss)

    up_func = mean_loss + (2 * std_loss)
    down_func = mean_loss - (2 * std_loss)
    a = a/100
    fig = go.Figure(
        (go.Scatter(x=a, y=mean_loss, mode="markers+lines",
                    name="Mean Prediction", line=dict(dash="dash"),
                    marker=dict(color="green", opacity=.7),
                    ),
         go.Scatter(x=a, y=up_func, fill=None,
                    mode="lines", line=dict(color="lightgrey"),
                    showlegend=False),
         go.Scatter(x=a, y=down_func, fill='tonexty',
                    mode="lines", line=dict(color="lightgrey"),
                    showlegend=False),),
          layout=go.Layout(title="mean loss as a function of p% and condifence"
                                 " interval"
                           , xaxis_title = "p% of train set"
                           , yaxis_title = "mean loss"))
    fig.show()

if __name__ == '__main__':

    np.random.seed(0)

    # Question 1 - Load and preprocessing of housing prices dataset
    X, y_true = load_data(r"C:\Users\romyb\Documents\iml\IML.HUJI\datasets\house_prices.csv")

    # Question 2 - Feature evaluation with respect to response
    feature_evaluation(X, y_true, r"C:\Users\romyb\Documents\iml\IML.HUJI\exercises\house_price_predic_plot2")

    # Question 3 - Split samples into training- and testing sets.
    train_x, train_y, test_x, test_y = split_train_test(X,y_true)

    # Question 4
    fit_test_model(train_x, train_y, test_x, test_y)


text:
# Code from Chapter 10 of Machine Learning: An Algorithmic Perspective (2nd Edition)
# by Stephen Marsland (http://stephenmonika.net)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008, 2014

# An exhaustive search to solve the Knapsack problem
import numpy as np

def exhaustive():
    maxSize = 500
    sizes = np.array(
        [109.60, 125.48, 52.16, 195.55, 58.67, 61.87, 92.95, 93.14, 155.05, 110.89, 13.34, 132.49, 194.03, 121.29,
         179.33, 139.02, 198.78, 192.57, 81.66, 128.90])

    best = 0

    twos = np.arange(-len(sizes), 0, 1)
    twos = 2.0 ** twos

    for i in range(2 ** len(sizes) - 1):
        string = np.remainder(np.floor(i * twos), 2)
        fitness = np.sum(string * sizes)
        if fitness > best and fitness < 500:
            best = fitness
            bestString = string
    print(best)
    print(bestString)

exhaustive()

max_stars_repo_name: quietcoolwu/MLCode

max_issues_repo_name: quietcoolwu/MLCode

max_forks_repo_name: quietcoolwu/MLCode


text:
# Python3 file
# Created by Marissa Bennett

import math, sys, csv, ast, time, re, warnings
import numpy as np
import pandas as pd
import rpyca as rp
from helperFiles.logger import *
from helperFiles.matrixOp import *
from helperFiles.fileHandler import *
from helperFiles.plotter import *
from helperFiles.models import *
from helperFiles.configParser import *

# main function
if __name__ == '__main__':
    start_time = time.time()
    # Ask for configuration to use
    configType, con = setConfig()
    # Set log for debugging or other purposes (can be overridden)
    setLog(con['LogFile'])
    logMsg(1,"CONFIGURATION USED: %s" % str(configType))
    # Set all other configuration variables
    fileName = con['CSVFile']
    logMsg(1,"File Name: %s" % str(fileName))
    labelsName = re.sub(r'[^\w]', '', con['Labels'])
    #####
    ## TODO need to make these better and more reliable
#    onehot = toList(con['OneHot'], integer=False)
#    skip = toList(con['Skip'], integer=False)
    ######
    seed = (0 if (con['RandomSeed'] == 0) else con['RandomSeed'])
    sample = (0 if (con['SampleSize'] == 0) else con['SampleSize'])
    ratioTrain, ratioValid = con['RatioTrainData'], con['RatioValidData']
    # Set ML model to run
    toRun = [con['Models']]
    if "all" == con['Models']:
        # NOTE these are not all the models in the model.py file
        toRun = ['rf','knn','logreg','svm','dtree','nb','kmeans','gb','pynn']
    # Set Looping actions
    howToRun = []
    mode = con['Mode']
    if mode == 1:
        howToRun = [con['LambdaStartValue']]
    elif mode == 2: # this is used for plotting
        howToRun = [con['LambdaStartValue']] * 10
    else:           # default for finding a good lambda
        howToRun = frange(con['LambdaStartValue'], con['LambdaEndValue'], con['LambdaIncrValue'])

    # ensures preprocessing happens at least once
    # TODO look into if I could just randomize the random data again instead???
    pre = True

    Xlis,LSlis,XLSlis = [], [], []

    # main loop
    # TODO normalize each matrix with X1 things (see paper)
    for l in howToRun:
        if not mode == 0 or pre:
            if "ISCX" in fileName:  # TODO these should be changed into a function or something in the future
                skip = ['FlowID', 'SourceIP', 'Timestamp', 'Label']
                [X1, X2, X3], ymat = preproc(fileName, labelsName, sample, seed, ratioTrain, ratioValid, skip=skip)# onehot, skip)
            elif "LLS_DDOS" in fileName:
                skip = ['No.', 'Label']
                [X1, X2, X3], ymat = preprocLLSDOS(fileName, labelsName, sample, seed, ratioTrain, ratioValid, skip=skip)# onehot, skip)
#            [X1, X2, X3], ymat = preprocKaggle(fileName, labelsName, sample, seed, ratioTrain, ratioValid, onehot, skip)
            pre = False     # done preprocessing for mode 0 only!

        # XXX
        #plotU(X1, ymat[0])

        logMsg(1, "Lambda: %s" % (str(l)))
        print("\n\nLAMBDA: ", l)

        # runs RPCA
        [LS1, LS2], [XLS1, XLS2] = rp.rpca(X1, X2, l)
        # XXX Future Work: see if lambda can be tuned outside of using ML models??

        # ML/AI loop
        for m in toRun:
            print("running ML")
            Xmat, LSmat, XLSmat, ymatX12 = [X1, X2], [LS1, LS2], [XLS1, XLS2], [ymat[0], ymat[1]]
            res, dall = runModels(Xmat, LSmat, XLSmat, ymatX12, code=m)

            # Validates ONLY if a good f1 score occurred
            if res:
                print("Validating...")
                logMsg(1, "Validating GOOD Lambda: %s" % (str(l)))

                # validate
                [LS1, LS3], [XLS1, XLS3] = rp.project(X1, X3)

                # ML/AI
                Xmat, LSmat, XLSmat, ymatX13 = [X1, X3], [LS1, LS3], [XLS1, XLS3], [ymat[0], ymat[2]]
                res, dgood = runModels(Xmat, LSmat, XLSmat, ymatX13, code=m)

                Xlis.append(dgood[0])
                LSlis.append(dgood[1])
                XLSlis.append(dgood[2])

    generateResults(toRun[0],l,Xlis,LSlis,XLSlis)
    logMsg(1, "Time to complete: %s" % str(time.time() - start_time))


text:
# Copyright (c) 2020, NVIDIA CORPORATION.
#
# 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 cudf
import cupy
import pandas as pd
import numpy as np
import torch
import s3fs
import transformers
from clx.analytics.cybert import Cybert

S3_BASE_PATH = "models.huggingface.co/bert/raykallen/cybert_apache_parser"
CONFIG_FILENAME = "config.json"
MODEL_FILENAME = "pytorch_model.bin"

fs = s3fs.S3FileSystem(anon=True)
fs.get(S3_BASE_PATH + "/" + MODEL_FILENAME, MODEL_FILENAME)
fs.get(S3_BASE_PATH + "/" + CONFIG_FILENAME, CONFIG_FILENAME)

cyparse = Cybert()

input_logs = cudf.Series(['109.169.248.247 - -',
                          'POST /administrator/index.php HTTP/1.1 200 4494'])

def get_expected_preprocess():
    tokens = torch.tensor(
        [[11523, 119, 20065, 119, 27672, 119, 26049, 118, 118, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
         [153, 9025, 1942, 120, 11065, 120, 7448, 119, 185, 16194, 145, 20174,
          2101, 120, 122, 119, 122, 2363, 3140, 1580, 1527, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0]], device='cuda:0'
    )

    masks = torch.tensor(
        [[1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0],
         [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 0, 0, 0, 0]], device='cuda:0'
    )

    metadata = cupy.array([[0, 0, 8], [1, 0, 20]], dtype='uint32')
    return tokens, masks, metadata

def get_expected_inference():
    expected_parsed_df = pd.DataFrame({
        'remote_host': ['109.169.248.247', np.NaN],
        'other': ['-', np.NaN],
        'request_method': [np.NaN, 'POST'],
        'request_url': [np.NaN, "/administrator/index.php"],
        'request_http_ver': [np.NaN, 'HTTP/1.1'],
        'status': [np.NaN, '200'],
        'response_bytes_clf': [np.NaN, '449']
    })

    expected_confidence_df = pd.DataFrame({
        'remote_host': [0.999628, np.NaN], 'other': [0.999579, np.NaN],
        'request_method': [np.NaN, 0.99822], 'request_url': [np.NaN, 0.999629],
        'request_http_ver': [np.NaN, 0.999936], 'status': [np.NaN, 0.999866],
        'response_bytes_clf': [np.NaN, 0.999751]
    })
    return expected_parsed_df, expected_confidence_df

def test_load_model():
    cyparse.load_model(MODEL_FILENAME, CONFIG_FILENAME)
    assert isinstance(cyparse._label_map, dict)
    assert isinstance(cyparse._model,
                      transformers.modeling_bert.BertForTokenClassification)

def test_preprocess():
    expected_tokens, expected_masks, expected_metadata = get_expected_preprocess()
    actual_tokens, actual_masks, actual_metadata = cyparse.preprocess(input_logs)
    assert actual_tokens.equal(expected_tokens)
    assert actual_masks.equal(expected_masks)
    assert cupy.equal(actual_metadata, expected_metadata).all()

def test_inference():
    if torch.cuda.is_available():
        expected_parsed_df, expected_confidence_df = get_expected_inference()
        actual_parsed_df, actual_confidence_df = cyparse.inference(input_logs)
        pd._testing.assert_frame_equal(actual_parsed_df, expected_parsed_df)
        pd._testing.assert_frame_equal(actual_confidence_df, expected_confidence_df)

max_stars_repo_name: mdemoret-nv/clx

max_issues_repo_name: mdemoret-nv/clx

max_forks_repo_name: mdemoret-nv/clx

reason: import numpy,import cupy


text:
"""
=======================================================
Lasso path computation on Leukemia dataset
=======================================================

The example runs the Celer algorithm on the Leukemia
dataset which is a dense dataset.

Running time is compared with the scikit-learn implementation.
"""

import time

import numpy as np
import matplotlib.pyplot as plt

from sklearn.linear_model import lasso_path
from sklearn.datasets import fetch_openml

from celer import celer_path
from celer.plot_utils import plot_path_hist

print(__doc__)

print("Loading data...")
dataset = fetch_openml("leukemia")
X = np.asfortranarray(dataset.data.astype(float))
y = 2 * ((dataset.target != "AML") - 0.5)
n_samples = len(y)

y -= np.mean(y)
y /= np.std(y)

print("Starting path computation...")
alpha_max = np.max(np.abs(X.T.dot(y))) / n_samples

fine = True  # fine or coarse grid
n_alphas = 100 if fine else 10
alphas = alpha_max * np.logspace(0, -2, n_alphas)

gap_freq = 10
prune = 1
verbose = 0
verbose_inner = 0
tols = [1e-2, 1e-4, 1e-6, 1e-8]
results = np.zeros([2, len(tols)])
for tol_ix, tol in enumerate(tols):
    t0 = time.time()
    res = celer_path(X, y, pb='lasso', alphas=alphas, max_iter=100, gap_freq=gap_freq,
                     p0=100, verbose=verbose, verbose_inner=verbose_inner,
                     tol=tol, prune=prune, return_thetas=True)
    results[0, tol_ix] = time.time() - t0
    print('Celer time: %.2f s' % results[0, tol_ix])
    _, coefs, gaps, thetas = res

    t0 = time.time()
    _, coefs, dual_gaps = lasso_path(X, y, tol=tol, alphas=alphas)
    results[1, tol_ix] = time.time() - t0
    coefs = coefs.T

labels = [r"\sc{Celer}", "scikit-learn"]
figsize = (7.1, 4.3)
fig = plot_path_hist(results, labels, tols, figsize, ylim=None)
plt.show()

max_stars_repo_name: yjkimnada/celer

max_issues_repo_name: yjkimnada/celer

max_forks_repo_name: yjkimnada/celer


text:
from __future__ import division
from __future__ import print_function

from util import LoadData, Load, Save, DisplayPlot
import sys
import numpy as np
import matplotlib.pyplot as plt

def InitNN(num_inputs, num_hiddens, num_outputs):
    """Initializes NN parameters.

    Args:
        num_inputs:    Number of input units.
        num_hiddens:   List of two elements, hidden size for each layer.
        num_outputs:   Number of output units.

    Returns:
        model:         Randomly initialized network weights.
    """
    W1 = 0.1 * np.random.randn(num_inputs, num_hiddens[0])
    W2 = 0.1 * np.random.randn(num_hiddens[0], num_hiddens[1])
    W3 = 0.01 * np.random.randn(num_hiddens[1], num_outputs)
    b1 = np.zeros((num_hiddens[0]))
    b2 = np.zeros((num_hiddens[1]))
    b3 = np.zeros((num_outputs))
    v1 = np.zeros_like(W1)
    v2 = np.zeros_like(W2)
    v3 = np.zeros_like(W3)
    v4 = np.zeros_like(b1)
    v5 = np.zeros_like(b2)
    v6 = np.zeros_like(b3)

    model = {
        'W1': W1,
        'W2': W2,
        'W3': W3,
        'b1': b1,
        'b2': b2,
        'b3': b3,
        'v1': v1,
        'v2': v2,
        'v3': v3,
        'v4': v4,
        'v5': v5,
        'v6': v6
    }
    return model

def Affine(x, w, b):
    """Computes the affine transformation.

    Args:
        x: Inputs
        w: Weights
        b: Bias

    Returns:
        y: Outputs
    """
    # y = np.dot(w.T, x) + b
    y = x.dot(w) + b
    return y

def AffineBackward(grad_y, x, w):
    """Computes gradients of affine transformation.

    Args:
        grad_y: gradient from last layer
        x: inputs
        w: weights

    Returns:
        grad_x: Gradients wrt. the inputs.
        grad_w: Gradients wrt. the weights.
        grad_b: Gradients wrt. the biases.
    """
    grad_x = grad_y.dot(w.T)
    grad_w = x.T.dot(grad_y)
    grad_b = np.sum(grad_y, axis=0)
    return grad_x, grad_w, grad_b
    ###########################
    # Insert your code here.
    # grad_x = ...
    # grad_w = ...
    # grad_b = ...
    # return grad_x, grad_w, grad_b
    ###########################
    raise Exception('Not implemented')

def ReLU(x):
    """Computes the ReLU activation function.

    Args:
        x: Inputs

    Returns:
        y: Activation
    """
    return np.maximum(x, 0.0)

def ReLUBackward(grad_y, x, y):
    """Computes gradients of the ReLU activation function.

    Returns:
        grad_x: Gradients wrt. the inputs.
    """
    grad_x= None
    grad_x= grad_y
    grad_x[x<=0]=0
    return grad_x
    ###########################
    # Insert your code here.
    # grad_x = ...
    # return grad_x
    ###########################
    raise Exception('Not implemented')

def Softmax(x):
    """Computes the softmax activation function.

    Args:
        x: Inputs

    Returns:
        y: Activation
    """
    return np.exp(x) / np.exp(x).sum(axis=1, keepdims=True)

def NNForward(model, x):
    """Runs the forward pass.

    Args:
        model: Dictionary of all the weights.
        x:     Input to the network.

    Returns:
        var:   Dictionary of all intermediate variables.
    """
    h1 = Affine(x, model['W1'], model['b1'])
    h1r = ReLU(h1)
    h2 = Affine(h1r, model['W2'], model['b2'])
    h2r = ReLU(h2)
    y = Affine(h2r, model['W3'], model['b3'])
    var = {
        'x': x,
        'h1': h1,
        'h1r': h1r,
        'h2': h2,
        'h2r': h2r,
        'y': y
    }
    return var

def NNBackward(model, err, var):
    """Runs the backward pass.

    Args:
        model:    Dictionary of all the weights.
        err:      Gradients to the output of the network.
        var:      Intermediate variables from the forward pass.
    """
    dE_dh2r, dE_dW3, dE_db3 = AffineBackward(err, var['h2r'], model['W3'])
    dE_dh2 = ReLUBackward(dE_dh2r, var['h2'], var['h2r'])
    dE_dh1r, dE_dW2, dE_db2 = AffineBackward(dE_dh2, var['h1r'], model['W2'])
    dE_dh1 = ReLUBackward(dE_dh1r, var['h1'], var['h1r'])
    _, dE_dW1, dE_db1 = AffineBackward(dE_dh1, var['x'], model['W1'])
    model['dE_dW1'] = dE_dW1
    model['dE_dW2'] = dE_dW2
    model['dE_dW3'] = dE_dW3
    model['dE_db1'] = dE_db1
    model['dE_db2'] = dE_db2
    model['dE_db3'] = dE_db3
    pass

def NNUpdate(model, eps, momentum):
    """Update NN weights.

    Args:
        model:    Dictionary of all the weights.
        eps:      Learning rate.
        momentum: Momentum.
    """
    model['v1'] = momentum*model['v1']+eps*model['dE_dW1']
    model['v2'] = momentum*model['v2']+eps*model['dE_dW2']
    model['v3'] = momentum*model['v3']+eps*model['dE_dW3']
    model['v4'] = momentum*model['v4']+eps*model['dE_db1']
    model['v5'] = momentum*model['v5']+eps*model['dE_db2']
    model['v6'] = momentum*model['v6']+eps*model['dE_db3']
    model['W1'] = model['W1']-model['v1']
    model['W2'] = model['W2']-model['v2']
    model['W3'] = model['W3']-model['v3']
    model['b1'] = model['b1']-model['v4']
    model['b2'] = model['b2']-model['v5']
    model['b3'] = model['b3']-model['v6']
 _

def Train(model, forward, backward, update, eps, momentum, num_epochs,
          batch_size):
    """Trains a simple MLP.

    Args:
        model:           Dictionary of model weights.
        forward:         Forward prop function.
        backward:        Backward prop function.
        update:          Update weights function.
        eps:             Learning rate.
        momentum:        Momentum.
        num_epochs:      Number of epochs to run training for.
        batch_size:      Mini-batch size, -1 for full batch.

    Returns:
        stats:           Dictionary of training statistics.
            - train_ce:       Training cross entropy.
            - valid_ce:       Validation cross entropy.
            - train_acc:      Training accuracy.
            - valid_acc:      Validation accuracy.
    """
    inputs_train, inputs_valid, inputs_test, target_train, target_valid, \
        target_test = LoadData('toronto_face.npz')
    rnd_idx = np.arange(inputs_train.shape[0])
    train_ce_list = []
    valid_ce_list = []
    train_acc_list = []
    valid_acc_list = []
    num_train_cases = inputs_train.shape[0]
    if batch_size == -1:
        batch_size = num_train_cases
    num_steps = int(np.ceil(num_train_cases / batch_size))
    for epoch in range(num_epochs):
        np.random.shuffle(rnd_idx)
        inputs_train = inputs_train[rnd_idx]
        target_train = target_train[rnd_idx]
        for step in range(num_steps):
            # Forward prop.
            start = step * batch_size
            end = min(num_train_cases, (step + 1) * batch_size)
            x = inputs_train[start: end]
            t = target_train[start: end]

            var = forward(model, x)
            prediction = Softmax(var['y'])

            train_ce = -np.sum(t * np.log(prediction)) / x.shape[0]
            train_acc = (np.argmax(prediction, axis=1) ==
                         np.argmax(t, axis=1)).astype('float').mean()
            print(('Epoch {:3d} Step {:2d} Train CE {:.5f} '
                   'Train Acc {:.5f}').format(
                epoch, step, train_ce, train_acc))

            # Compute error.
            error = (prediction - t) / x.shape[0]

            # Backward prop.
            backward(model, error, var)

            # Update weights.
            update(model, eps, momentum)

        valid_ce, valid_acc = Evaluate(
            inputs_valid, target_valid, model, forward, batch_size=batch_size)
        print(('Epoch {:3d} '
               'Validation CE {:.5f} '
               'Validation Acc {:.5f}\n').format(
            epoch, valid_ce, valid_acc))
        train_ce_list.append((epoch, train_ce))
        train_acc_list.append((epoch, train_acc))
        valid_ce_list.append((epoch, valid_ce))
        valid_acc_list.append((epoch, valid_acc))
        DisplayPlot(train_ce_list, valid_ce_list, 'Cross Entropy', number=0)
        DisplayPlot(train_acc_list, valid_acc_list, 'Accuracy', number=1)

    print()
    train_ce, train_acc = Evaluate(
        inputs_train, target_train, model, forward, batch_size=batch_size)
    valid_ce, valid_acc = Evaluate(
        inputs_valid, target_valid, model, forward, batch_size=batch_size)
    test_ce, test_acc = Evaluate(
        inputs_test, target_test, model, forward, batch_size=batch_size)
    print('CE: Train %.5f Validation %.5f Test %.5f' %
          (train_ce, valid_ce, test_ce))
    print('Acc: Train {:.5f} Validation {:.5f} Test {:.5f}'.format(
        train_acc, valid_acc, test_acc))

    stats = {
        'train_ce': train_ce_list,
        'valid_ce': valid_ce_list,
        'train_acc': train_acc_list,
        'valid_acc': valid_acc_list
    }

    return model, stats

def Evaluate(inputs, target, model, forward, batch_size=-1):
    """Evaluates the model on inputs and target.

    Args:
        inputs: Inputs to the network.
        target: Target of the inputs.
        model:  Dictionary of network weights.
    """
    num_cases = inputs.shape[0]
    if batch_size == -1:
        batch_size = num_cases
    num_steps = int(np.ceil(num_cases / batch_size))
    ce = 0.0
    acc = 0.0
    for step in range(num_steps):
        start = step * batch_size
        end = min(num_cases, (step + 1) * batch_size)
        x = inputs[start: end]
        t = target[start: end]
        prediction = Softmax(forward(model, x)['y'])
        ce += -np.sum(t * np.log(prediction))
        acc += (np.argmax(prediction, axis=1) == np.argmax(
            t, axis=1)).astype('float').sum()
    ce /= num_cases
    acc /= num_cases
    return ce, acc

def CheckGrad(model, forward, backward, name, x):
    """Check the gradients

    Args:
        model: Dictionary of network weights.
        name: Weights name to check.
        x: Fake input.
    """
    np.random.seed(0)
    var = forward(model, x)
    loss = lambda y: 0.5 * (y ** 2).sum()
    grad_y = var['y']
    backward(model, grad_y, var)
    grad_w = model['dE_d' + name].ravel()
    w_ = model[name].ravel()
    eps = 1e-7
    grad_w_2 = np.zeros(w_.shape)
    check_elem = np.arange(w_.size)
    np.random.shuffle(check_elem)
    # Randomly check 20 elements.
    check_elem = check_elem[:20]
    for ii in check_elem:
        w_[ii] += eps
        err_plus = loss(forward(model, x)['y'])
        w_[ii] -= 2 * eps
        err_minus = loss(forward(model, x)['y'])
        w_[ii] += eps
        grad_w_2[ii] = (err_plus - err_minus) / 2 / eps
    np.testing.assert_almost_equal(grad_w[check_elem], grad_w_2[check_elem],
                                   decimal=3)

def main():
    """Trains a NN."""
    model_fname = 'nn_model.npz'
    stats_fname = 'nn_stats.npz'

    # Hyper-parameters. Modify them if needed.
    num_hiddens = [16, 32]
    eps = 0.01
    momentum = 0.3
    num_epochs = 1000
    batch_size = 200

    # Input-output dimensions.
    num_inputs = 2304
    num_outputs = 7

    # Initialize model.
    model = InitNN(num_inputs, num_hiddens, num_outputs)

    # Uncomment to reload trained model here.
    # model = Load(model_fname)

    # Check gradient implementation.
    print('Checking gradients...')
    x = np.random.rand(10, 48 * 48) * 0.1
    CheckGrad(model, NNForward, NNBackward, 'W3', x)
    CheckGrad(model, NNForward, NNBackward, 'b3', x)
    CheckGrad(model, NNForward, NNBackward, 'W2', x)
    CheckGrad(model, NNForward, NNBackward, 'b2', x)
    CheckGrad(model, NNForward, NNBackward, 'W1', x)
    CheckGrad(model, NNForward, NNBackward, 'b1', x)

    # Train model.
    stats = Train(model, NNForward, NNBackward, NNUpdate, eps,
                  momentum, num_epochs, batch_size)

    # Uncomment if you wish to save the model.
    # Save(model_fname, model)

    # Uncomment if you wish to save the training statistics.
    # Save(stats_fname, stats)

    #plot the first layer weight
    #plt.figure(2)
    #plt.clf()
    #plt.subplot(1, model['W1'].shape[1], 2)
    #plt.imshow(model['W1'][:, 1].reshape(48, 48), cmap=plt.cm.gray)
    #plt.draw()
    #plt.pause(1)

if __name__ == '__main__':
    main()

max_stars_repo_path: Nerual Network/nn.py

max_issues_repo_path: Nerual Network/nn.py

max_forks_repo_path: Nerual Network/nn.py


text:
import sys
sys.path.append("..") #this is to add the avobe folder to the package directory
import geopandas as gpd
import pandas as pd
import numpy as np
import nexustool.weap_tools as wp
import re
from functools import reduce
import os
from shutil import copyfile

def extract_points_data(results, spacial_data, points, melt_on, regex):
    df = pd.DataFrame()
    for key, row in spacial_data.groupby(points):
        df_temp = results.filter(regex=regex.format(key)).copy()
        df_temp.columns = [i.split('\\')[-1] for i in df_temp.columns]
        df_temp[points] = key
        for i in melt_on:
            df_temp[i] = results[i]
        melt_vars = melt_on.copy()
        melt_vars.append(points)
        df = df.append(df_temp.melt(id_vars=melt_vars), sort=False)
    return df.reset_index(drop=True)

def extract_within(s, left, rigth):
    mo = re.search(f'{left}(.*){rigth}', s)
    if mo:
        return mo.group(1)
    return ''

def special_conditions(df, sheet_name):
    if sheet_name == 'Pipelines':
        df["PL_ZaraMain2SZ09SZ06 0 \ Reach"] = df["PL_ZaraMain2SZ09SZ06 0 \ Headflow"]

def get_data(sheet_names, data, spatial_data, variable, regex,
             rename={'Unnamed: 0': 'Date'}, melt_on=['Date', 'Year', 'Month'],
             look_in_sector=False):
    merged_data = pd.DataFrame()
    for sheet_name, value in sheet_names.items():
        data_temp = data.parse(sheet_name, skiprows=3)
        data_temp = data_temp.loc[data_temp['Unnamed: 0']!='Sum']
        data_temp.rename(columns=rename, inplace=True)
        if 'Date' in data_temp.columns:
            data_temp['Date'] = pd.to_datetime(data_temp['Date'])
            data_temp['Year'] = data_temp['Date'].dt.year
            data_temp['Month'] = data_temp['Date'].dt.month
        data_temp.columns = data_temp.columns.str.replace('"', '')
        data_temp.columns = data_temp.columns.str.replace('  ', ' ')
        data_temp.columns = data_temp.columns.str.strip()
        special_conditions(data_temp, sheet_name)
        if look_in_sector:
            merged_data_temp = extract_points_data(data_temp,
                                                   spatial_data.loc[spatial_data['type']==value],
                                                   variable, melt_on, regex)
        else:
            merged_data_temp = extract_points_data(data_temp,
                                                   spatial_data,
                                                   variable, melt_on, regex)
        merged_data_temp['variable'] = [i.split('[')[0].strip() for i in merged_data_temp.variable]
        merged_data_temp['type'] = value
        merged_data_temp.loc[merged_data_temp['type']=='variable', 'type'] = merged_data_temp.loc[merged_data_temp['type']=='variable', 'variable']

        merged_data = merged_data.append(merged_data_temp, sort=False)
    return merged_data

def enumerate_segments(df):
    x = len(df.Year.unique())*12
    df['n'] = 0
    for key, group in df.groupby('pipeline'):
        n = [[i]*x for i in range(int(group.shape[0]/x))]
        n = reduce(lambda a, b: a + b, n)
        df.loc[df.pipeline==key,'n'] = n

def get_elevation_delta(df):
    df['elevation_delta'] = 0
    for pipeline in df.pipeline.unique():
        for n in range(df.loc[df.pipeline==pipeline, 'n'].max()):
            elevation_1 = df.loc[(df.pipeline==pipeline) & (df.n==n), 'elevation']
            elevation_2 = df.loc[(df.pipeline==pipeline) & (df.n==(n+1)), 'elevation']
            df.loc[(df.pipeline==pipeline) & (df.n==n), 'elevation_delta'] = np.array(elevation_2) - np.array(elevation_1)

def save_dataframe(df, init_year, end_year, output_folder, file_name):
    df.loc[(df.Year>=init_year) & (df.Year<=end_year)].to_csv(os.path.join(output_folder, file_name), index=False)

max_stars_repo_name: KTH-dESA/FAO

max_issues_repo_name: KTH-dESA/FAO

max_forks_repo_name: KTH-dESA/FAO


text:
#!/usr/bin/python
#
# Copyright 2019 Google LLC
#
# 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 tensorflow as tf
import numpy as np
import logging
import pickle
import os
from annoy import AnnoyIndex

VECTOR_LENGTH = 512
METRIC = 'angular'

def build_index(embedding_files_pattern, index_filename,
                num_trees=100):

  annoy_index = AnnoyIndex(VECTOR_LENGTH, metric=METRIC)
  mapping = {}

  embed_files = tf.gfile.Glob(embedding_files_pattern)[:250]
  logging.info('{} embedding files are found.'.format(len(embed_files)))

  item_counter = 0
  for f, embed_file in enumerate(embed_files):
    logging.info('Loading embeddings in file {} of {}...'.format(
      f, len(embed_files)))
    record_iterator = tf.python_io.tf_record_iterator(
      path=embed_file)

    for string_record in record_iterator:
      example = tf.train.Example()
      example.ParseFromString(string_record)
      string_identifier = example.features.feature['id'].bytes_list.value[0]
      mapping[item_counter] = string_identifier
      embedding = np.array(
        example.features.feature['embedding'].float_list.value)
      annoy_index.add_item(item_counter, embedding)
      item_counter += 1

    logging.info('Loaded {} items to the index'.format(item_counter))

  logging.info('Start building the index with {} trees...'.format(num_trees))
  annoy_index.build(n_trees=num_trees)
  logging.info('Index is successfully built.')
  logging.info('Saving index to disk...')
  annoy_index.save(index_filename)
  logging.info('Index is saved to disk.')
  logging.info("Index file size: {} GB".format(
    round(os.path.getsize(index_filename) / float(1024 ** 3), 2)))
  annoy_index.unload()
  logging.info('Saving mapping to disk...')
  with open(index_filename + '.mapping', 'wb') as handle:
    pickle.dump(mapping, handle, protocol=pickle.HIGHEST_PROTOCOL)
  logging.info('Mapping is saved to disk.')
  logging.info("Mapping file size: {} MB".format(
    round(os.path.getsize(index_filename + '.mapping') / float(1024 ** 2), 2)))