text:
from numpy import array, matrix, inf, nan, ones, zeros, isinf, isfinite
from numpy import linspace, matmul, sign

import matplotlib

matplotlib.use('nbAgg')

import matplotlib.pyplot as plt

def get_image_pos(object_distance, focal_length):
    s = object_distance
    f = focal_length
    if s == f:
        # print('Image at infinity')
        return sign(s*f) * inf
    return 1 / (1 / f - 1 / s)

def get_image_size(object_size, object_distance, focal_length):
    A = object_size
    s = object_distance
    f = focal_length
    b = get_image_pos(s, f)

    return A * b / s

class OPE(object):
    def __init__(self, a=inf, d=None, f=None, name=''):
        """
        Optical path element used for ray tracing.

        Arguments
        ---------
        a : inf or float
            Aperture of the element.
        d : float or None
            Distance travelled through optical element.
        f : flaot or None
            Focal length of optical element.
        name : str
            Name of the element.
        """
        self.aperture = abs(a)
        self.name = name
        self.focal_length = f if f else inf
        self.distance = d if d else 0

        if f:
            self.c = -1 / f
        else:
            self.c = 0

        if d:
            self.b = d
        else:
            self.b = 0

    def copy(self, name=None):
        """
        Returns a copy of the OPE.
        """
        a = self.aperture
        name = name if name else f'{self.name}_copy'
        f = self.focal_length
        d = self.distance

        return OPE(a=a, d=d, f=f, name=name)

    def is_lens(self):
        if self.get_matrix()[1, 0] != 0:
            return True
        else:
            return False

    def get_matrix(self):
        return matrix([[1, self.b], [self.c, 1]])

    def passes_aperture(self, r, verbose=False):
        """
        Does an input ray r pass the optical element.
        """
        h = r.flatten()[0]
        if abs(h) > self.aperture:
            if verbose:
                print('{0:s}: stopped at aperture with height {1:1.3f}'.format(
                    self.name, h))
            return False
        else:
            if verbose:
                print('{0:s}: passed aperture with height {1:1.3f}'.format(
                    self.name, h))
            return True

    def pass_aperture(self, r, verbose=False):
        """ Return nan vector if ray cannot pass aperture."""
        if self.passes_aperture(r, verbose=verbose):
            return r
        else:
            return matrix(ones(r.shape) * nan)

    def transmit(self, r, verbose=False):
        """ Calculate output ray. """
        M = self.get_matrix()
        r_ = M * matrix(self.pass_aperture(r, verbose=verbose)).T
        return array(r_).flatten()

    def get_travel_length(self):
        """ Length the optical element takes. """
        return self.b

    def has_aperture(self):
        """ Return True is OPE has an aperture. """
        if self.aperture < inf:
            return True
        else:
            return False

def trace_ray(r, sequence):
    """Take r as input vector and trace the ray through the sequence of
    optical path elements OLEs.

    Arguments
    ---------
    r : 2-tuple
        Vector with the first element beeing the ray height and the
        second being the angle
    sequence : list of optical path elements (OPE)

    Returns
    -------
    (dist, ray)
    dist : array
        Array representing the position of the ray through the trace.
    ray : array
        Array representing the ray vectors at the corresponding
        distances.
    """
    ray0 = zeros((2, len(sequence) + 1))
    ray0[:, 0] = r

    # Do the calculations
    distances = [0.0]
    for idx, ope in enumerate(sequence):
        ray0[:, idx + 1] = ope.transmit(ray0[:, idx])
        distances.append(distances[idx] + ope.get_travel_length())

    dist = array(distances)

    return (dist, ray0)

def remove_apertures(sequence):
    """
    Return a copy of the sequence without apertures.
    """
    seq = []
    for ope in sequence:
        new_ope = ope.copy()
        new_ope.aperture = inf
        seq.append(new_ope)
    return seq

def get_first_aperture(sequence):
    """
    get aperture of first lens in path and calculate distance to it.

    Arguments
    ---------
    sequence : list of OPE
        List of optical path elements

    Returns
    -------
    tuple(idx, distance, aperture)
    """
    d = 0.0
    a = inf
    for idx, m in enumerate(sequence):
        if isinf(m.aperture):
            d += m.get_travel_length()
        else:
            a = m.aperture
            break
    return (idx, d, a)

def get_pos_at_idx(idx, sequence):
    """
    Return the position of the element at the given index.
    """
    d = sum([ope.get_travel_length() for ope in sequence[:idx + 1]])

    return d

def get_aperture_pos(sequence):
    """
    Calculate positions of apertures.

    Returns
    -------
    List of tuples with index and position of OPE in sequence.
    """
    d = 0.0
    d_ = []
    for idx, ope in enumerate(sequence):
        if ope.has_aperture():
            d_.append((idx, d))
        else:
            d += ope.get_travel_length()

    return d_

def get_lens_pos(sequence):
    """
    Calculate positions of lenses.

    Returns
    -------
    List of tuples with index and position of OPE in sequence.
    """
    d = 0.0
    d_ = []
    for idx, ope in enumerate(sequence):
        if ope.is_lens():
            d_.append((idx, d))
        else:
            d += ope.get_travel_length()

    return d_

def get_angle_lim(h, d, a):
    """
    Calculate the upper and lower angles of a source that is at
    distance d to an aperture a.
    """
    # return (arctan((a - h) / d), arctan(-(a + h) / d))
    return (-(a - h) / d, (a + h) / d)

def optical_system_parser(s):
    '''
    Convert a string into a seqeunce of optical path elements.

    Example
    -------
    s = 'd15 | l15/5.5 | d15'
    would result in a seqeunce of a distance element of 15 lu
    (lu=length units), followed by a Lens with focal length 15 lu
    with aperture 5.5 lu and another distance of 15 lu.
    '''
    if not isinstance(s, str):
        raise TypeError('Input string does not seem to be a string.')

    sequence = []

    for idx, si in enumerate(s.lower().replace(' ', '').split('|')):
        if si.startswith('d'):
            s_ = si.split('/')
            d = float(s_[0][1:])
            name = f'D_{int(round(d)):d}'
            if len(s_) == 2:
                a = float(s_[1])
                ope = OPE(d=d, a=a, name=name)
            else:
                ope = OPE(d=d, name=name)
        elif si.startswith('l'):
            s_ = si.split('/')
            f = float(s_[0][1:])
            name = f'L_{int(round(d)):d}'
            if len(s_) == 2:
                a = float(s_[1])
                ope = OPE(f=f, a=a, name=name)
            else:
                ope = OPE(f=f, name=name)
        elif si.startswith('a'):
            s_ = si.split('/')
            a = float(s_[0][1:])
            name = f'A_{int(round(d)):d}'
            ope = OPE(a=a, name=name)
        else:
            print('unknown element {}'.format(si))
        sequence.append(ope)

    return sequence

def get_max_aperture(sequence):
    """
    Return largest aperture in sequence.

    Returns None if no apertures are defined.
    """
    apertures = [abs(ope.aperture) for ope in sequence if ope.aperture != inf]

    if not apertures:
        out = None
    else:
        out = max(apertures)

    return out

class Ray(object):
    """
    Defines properties of a ray.
    """
    def __init__(self, height, angle, intensity=1.0):
        """
        Create a ray which is defined by its height and its angle,

        The angle is measured between the horizontal and teh ray,
        counter-clockwise.
        """
        self.height = height
        self.angle = angle
        self.intensity = intensity

    @property
    def vector(self):
        return (self.height, self.angle)

class OpticalSystem(object):
    """
    Manages an optical system.
    """
    def __init__(self, sequence):
        """
        Create optical system.

        Sequence can be a string or a list of OPEs.
        """
        try:
            seq = optical_system_parser(sequence)
        except TypeError:
            try:
                seq = sequence[:]
            except:
                raise TypeError('Sequence is not a proper string nor a sequence of OPEs.')

        self.sequence = seq
        self.max_y = 0.0
        self._plot_axis = None
        self._statics_drawn = None

    def _init_axis(self):
        self._statics_drawn = False
        axis = plt.figure().add_subplot(111)
        self._plot_axis = axis

    def get_matrix(self):
        """
        Return the tranfer matrix of the system.
        """
        M = matrix([[1, 0],
                    [0, 1]])
        mats = [ope.get_matrix() for ope in self.sequence]

        for mat in mats:
            M = matmul(mat, M)

        return M

    def add_ope(self, ope, index=None):
        """
        Add an optical path element to sequence of OPEs.

        Provide index to insert, leave to append.
        """
        if index:
            self.sequence.insert(index, ope)
        else:
            self.sequence.append(ope)

    @property
    def plot_axis(self):
        if not self._plot_axis:
            self._init_axis()
        return self._plot_axis

    @plot_axis.setter
    def plot_axis(self, axis):
        self._plot_axis = axis

    def set_max_y(self, y):
        if abs(y) > self.max_y:
            self.max_y = abs(y)

    def plot_statics(self, axis=None):
        """
        Plots statics, i.e. Lenses etc. and creates an axis if not
        given or present, yet.
        """
        if axis:
            self.plot_axis = axis
        else:
            axis = self.plot_axis

        if not self._statics_drawn:
            # draw optical axis
            axis.axhline(color='k', linewidth=0.5, linestyle='--')
            # draw lenses
            for idx, x in get_lens_pos(self.sequence):
                axis.axvline(x=x, ymin=0.02, ymax=0.98,
                             linewidth=0.5, linestyle='--')
            # draw apertures
            self.draw_apertures()

            self._statics_drawn = True

        return axis

    def draw_apertures(self):
        """
        Draw the apertures in the sequence
        """
        ax = self.plot_axis

        plt_kws = dict(linewidth=2, linestyle='-')

        a_max = get_max_aperture(self.sequence)

        if a_max:
            for idx, x in get_aperture_pos(self.sequence):
                a = self.sequence[idx].aperture
                if self.sequence[idx].is_lens():
                    plt_kws['color'] = 'darkslategrey'
                else:
                    plt_kws['color'] = 'teal'
                ax.plot([x, x], [-2*a_max, -a], **plt_kws)
                ax.plot([x, x], [a, 2*a_max], **plt_kws)
                self.set_max_y(2*a_max)

    def reset_plot(self):
        """ Reset internal plot axis."""
        self._statics_drawn = False

    def plot_rays(self, h, parallel=False, d=None, plot_fan=False,
                  axis=None, label=None, plot_statics=True, **pltkws):
        """
        Plot the ray trace through the sequence of OPEs.

        Arguments
        ---------
        h : float or iterable
            Height(s) of ray .
        sequence : list of OPE
            sequence of optical path elements.
        parallel : bool
            Whether the source ray is a parallel beam.
        d : float
            Diameter of an incoming parallel beam.
        plot_fan : bool or int
            Plot a fan of rays.
        axis : matplotlib.Axis
        label : str
            Label of the plotted line

        Keyword Arguments
        -----------------
        kws passed to plt.plot(**pltkws)

        Returns
        -------
        maplotlib.figure
        """
        try:
            heights = iter(h)
        except TypeError:
            heights = [h]

        if plot_statics:
            ax = self.plot_statics(axis=axis)
        else:
            ax = self.plot_axis

        if not any([a in pltkws for a in ['c', 'col', 'color']]):
            cycle_colors = True

        n = 1 * plot_fan

        for h in heights:
            if parallel:
                d = d or 1.0
                hs = list(linspace(h - d/2, h + d/2, n if n > 1 else 2))
                while True:
                    if cycle_colors:
                        pltkws['color'] = next(ax._get_lines.prop_cycler)['color']
                    try:
                        hin = hs.pop(0)
                        dist, r1 = trace_ray((hin, 0), self.sequence)
                        ax.plot(dist, r1[0, :], label=label or f'r = ({hin:1.2f}, 0)', **pltkws)
                        self.set_max_y(max(abs(r1[0, :])))
                    except:
                        # even
                        break
                    try:
                        dist, r2 = trace_ray((hs.pop(-1), 0), self.sequence)
                        ax.plot(dist, r2[0, :], **pltkws)
                        self.set_max_y(max(abs(r2[0, :])))
                    except:
                        # odd
                        break
            else:
                if self.has_aperture():
                    # get distance and aperture of first aperture
                    _, d, aperture = get_first_aperture(self.sequence)
                    a1, a2 = get_angle_lim(h, d, aperture)
                else:
                    try:
                        _, d = get_lens_pos(self.sequence)[0]
                    except:
                        d = inf
                    a1 = -h/d
                    a2 = 0

                if n > 0:
                    angles = list(linspace(a1, a2, n))
                else:
                    angles = [a1, a2]

                while True:
                    if cycle_colors:
                        pltkws['color'] = next(ax._get_lines.prop_cycler)['color']
                    try:
                        ang = angles.pop(0)
                        ray1 = (h, ang)
                        dist, r1 = trace_ray(ray1, self.sequence)
                        ax.plot(dist, r1[0, :], label=label or f'r =({h:1.2f}, {ang:1.2f})', **pltkws)
                        self.set_max_y(max(abs(r1[0, :])))
                    except:
                        # even
                        break
                    try:
                        ray2 = (h, angles.pop(-1))
                        dist, r2 = trace_ray(ray2, self.sequence)
                        ax.plot(dist, r2[0, :], **pltkws)
                        self.set_max_y(max(abs(r2[0, :])))
                    except:
                        # odd
                        break

        return ax

    def plot_ray(self, ray, axis=None, **pltkws):
        """
        Plot a ray.

        Arguments
        ---------
        ray : tuple(height, angle)

        Returns
        -------
        lines : plotted lines in axis
        """
        if not self._statics_drawn:
            ax = self.plot_statics(axis=axis)
        else:
            ax = self.plot_axis

        if not any([a in pltkws for a in ['c', 'col', 'color']]):
            pltkws['color'] = next(ax._get_lines.prop_cycler)['color']

        if not 'label' in pltkws:
            label = f'r =({ray[0]:1.2f}, {ray[1]:1.2f})'
        else:
            label = pltkws.pop('label')

        dist, r = trace_ray(ray, self.sequence)
        lines = ax.plot(dist, r[0, :], label=label, **pltkws)
        self.set_max_y(max(abs(r[0, :])))

        return lines

    @property
    def lens_positions(self):
        pos = list(zip(*get_lens_pos(self.sequence)))[1]
        return pos

    def print_ope_sequence(self):
        """
        Return a list of names of the OPEs in the sequence.
        """
        n = len(self.sequence)
        print('|'.join(f'{idx:^7d}' for idx in range(n)))
        print('|'.join(f'{name:^7s}' for name in self.names))

    @property
    def names(self):
        return [ope.name for ope in self.sequence]

    @property
    def aperture_positions(self):
        pos = list(zip(*get_aperture_pos(self.sequence)))[1]
        return pos

    @property
    def aperture_sizes(self):
        """ Return the aperture (half) sizes."""
        aps = [ope.aperture for ope in self.sequence if ope.has_aperture()]

        return aps

    def adjust_ylims(self):
        """
        Adjusts the y limits of the plot according to the apertures
        and rays.
        """
        max_y = self.max_y
        max_a = get_max_aperture(self.sequence)

        y = max([max_y, max_a])
        self.plot_axis.set_ylim([-0.9 * y, 0.9 * y])

    def get_idx_aperture_stop(self):
        """
        Find the optical element that defines the aperture stop of the
        system.

        Returns the index of the OPE in the sequence.
        """
        if get_max_aperture(self.sequence) is None:
            print('No apertures found.')
            out = None
        else:
            angle = 1e-10
            ctr = 0
            while True:
                ctr += 1
                _, rays = trace_ray((0, angle), self.sequence)
                if all(isfinite(rays[0])):
                    break
                else:
                    angle = angle / 2
                if ctr > 100:
                    raise Exception('Unable to trace ray through sequence.')

            ratio = 0.0
            for idx, ope in enumerate(self.sequence):
                ratio_ = abs(rays[0, idx+1]) / ope.aperture
                if ratio_ > ratio:
                    ratio = ratio_
                    out = idx

        return out

    def has_aperture(self):
        if any(ope.has_aperture() for ope in self.sequence):
            return True
        else:
            return False

    def get_aperture_stop_position(self, verbose=False):
        """
        Reduce sequence upto aperture and get distance from lens pos
        function.
        """
        _, d = get_aperture_pos(self.sequence[:self.get_idx_aperture_stop() + 1])[-1]

        if verbose:
            print(f'aperture stop position = {d:1.2f}')

        return d

    def get_aperture_stop_size(self, verbose=False):
        """
        Return the (half) aperture stop size.
        """
        a = self.sequence[self.get_idx_aperture_stop()].aperture
        if verbose:
            print(f'aperture half-diameter = {a:1.2f}')
        return a

    def indicate_aperture_stop(self, axis=None, color='orangered', verbose=False):
        """
        Draw the aperture stop in the ray tracing diagram.
        """
        if axis:
            ax = axis
        else:
            ax = self.plot_axis
        plt_kws = dict(linewidth=2.0, linestyle='-', color=color)

        x = self.get_aperture_stop_position(verbose=verbose)
        a = self.get_aperture_stop_size(verbose=verbose)
        y_max = self._get_y_max(verbose=verbose)

        ax.plot([x, x], [a, y_max], **plt_kws)
        ax.plot([x, x], [-y_max, -a], **plt_kws)

    def calc_entrance_pupil_position(self, verbose=False):
        """ sequence of OPEs preceeding the aperture stop """
        sequence_prec = self.sequence[:self.get_idx_aperture_stop()]
        d_ap = self.get_aperture_stop_position(verbose=verbose)
        x = d_ap
        mag = 1.0
        for idx, lens_pos in get_lens_pos(sequence_prec)[::-1]:
            # object distance
            d_obj = x - lens_pos
            # image distance
            d_img = get_image_pos(d_obj, sequence_prec[idx].focal_length)
            if verbose > 1:
                print(f'imaging lens position = {lens_pos:1.2f}')
                print(f'x_before = {x:1.2f}')
                print(f'd_obj = {d_obj:1.2f}')
                print(f'd_img = {d_img:1.2f}')

            x = lens_pos - d_img
            mag = mag * d_img / d_obj
            if verbose > 1:
                print(f'x_after = {x:1.2f}')
        return x

    def calc_entrance_pupil_size(self, verbose=False):
        """
        Return the size of the entrance pupil.
        """
        sequence_prec = self.sequence[:self.get_idx_aperture_stop()]
        d_ap = self.get_aperture_stop_position(verbose=verbose)
        x = d_ap
        mag = 1.0
        for idx, lens_pos in get_lens_pos(sequence_prec)[::-1]:
            # object distance
            d_obj = x - lens_pos
            # image distance
            d_img = get_image_pos(d_obj, sequence_prec[idx].focal_length)
            if verbose > 1:
                print(f'imaging lens position = {lens_pos:1.2f}')
                print(f'magnification = {mag:1.2f}')
                print(f'd_obj = {d_obj:1.2f}')
                print(f'd_img = {d_img:1.2f}')

            x = lens_pos - d_img
            mag = mag * d_img / d_obj
            if verbose > 1:
                print(f'magnification_after = {mag:1.2f}')

        en_pupil = self.get_aperture_stop_size() * abs(mag)

        return en_pupil

    def draw_entrance_pupil(self, axis=None, color='orangered', verbose=False):
        """
        Draw the apparent entrance pupil.
        """
        if axis:
            ax = axis
        else:
            ax = self.plot_axis
        x = self.calc_entrance_pupil_position(verbose=verbose)
        a = self.calc_entrance_pupil_size(verbose=verbose)
        y_max = self.max_y

        plt_kws = dict(linewidth=1.2, linestyle='-', color=color)
        ax.plot([x, x], [a, y_max], **plt_kws)
        ax.plot([x, x], [-y_max, -a], **plt_kws)

    def _get_y_max(self, verbose=False):
        """
        Return maximum y_value in plot.
        """
        if get_max_aperture(self.sequence):
            y_max = 2 * get_max_aperture(self.sequence)
        else:
            y_max = 2 * self.max_y

        if verbose > 1:
            print(f'y_max = {y_max:1.2f}')

        return y_max

    def get_NA(self):
        """
        Return the NA of the optical system.

        This is the entrance pupil size divided by the entrance pupil
        distance.
        """
        d = self.calc_entrance_pupil_position()
        a = self.calc_entrance_pupil_size()

        return a/d

    def trace_backward(self, idx, height, no_apertures=False, extra_distance=100.0, plot=True):
        """
        Trace a ray backward throught the optical system
        starting at the given index.

        Returns
        -------
        (positions, rays1, rays2)
        """
        # get only the sequence until the aperture and reverse it
        if no_apertures:
            seq_back = remove_apertures(self.sequence[:idx][::-1])
        else:
            seq_back = self.sequence[:idx][::-1]

        seq_back.append(OPE(d=extra_distance))

        # trace a ray from the given height or the edge of the
        # aperture with angle zero toward the next lense
        positions, rays1 = trace_ray((height, 0), seq_back)
        pos_back = get_pos_at_idx(idx, self.sequence) - positions

        if plot:
            pltkws = dict(linestyle='--', color='Grey')
            self.plot_axis.plot(pos_back, rays1[0], **pltkws)

        ## trace a ray from the given height or the edge of the
        ## aperture toward the center of the next lens
        # get index and position of the first element in the backward
        # sequence
        try:
            _, d = get_lens_pos(seq_back)[0]
        except IndexError:
            print('Nothing to track.')
            return None

        positions, rays2 = trace_ray((height, -height / d), seq_back)

        if plot:
            self.plot_axis.plot(pos_back, rays2[0], **pltkws)

        return (pos_back, rays1, rays2)

    def trace_forward(self, idx, height, no_apertures=False, extra_distance=100.0, plot=True):
        """
        Trace a ray throught the optical system starting at the given
        index.

        Returns
        -------
        (positions, rays1, rays2)
        """
        # get the sequence until the aperture and reverse it
        if no_apertures:
            seq = remove_apertures(self.sequence[idx:])
        else:
            seq = self.sequence[idx:]

        pos_offset = get_pos_at_idx(idx - 1, self.sequence)

        seq.append(OPE(d=extra_distance))

        # trace a ray from the given height or the edge of the
        # aperture with angle zero toward the next lens
        positions, rays1 = trace_ray((height, 0), seq)
        positions += pos_offset

        if plot:
            pltkws = dict(linestyle='--', color='Grey')
            self.plot_axis.plot(positions, rays1[0], **pltkws)

        ## trace a ray from the given height or the edge of the
        ## aperture toward the center of the next lens
        # get index and position of the first element in the backward
        # sequence
        try:
            _, d = get_lens_pos(seq)[0]
        except IndexError:
            print('Nothing to track.')
            return None

        positions, rays2 = trace_ray((height, -height / d), seq)
        positions += pos_offset

        if plot:
            self.plot_axis.plot(positions, rays2[0], **pltkws)

        return (positions, rays1, rays2)

    def _check_aperture_stop_rays_backward(self, extra_distance=100.0):
        """
        Re-trace the aperture stop.
        """
        idx = self.get_idx_aperture_stop()
        a = self.get_aperture_stop_size()
        self.trace_backward(idx, a, no_apertures=True,
                            extra_distance=extra_distance)

    def _check_aperture_stop_rays_forward(self, extra_distance=100.0):
        """
        Re-trace the aperture stop.
        """
        idx = self.get_idx_aperture_stop()
        a = self.get_aperture_stop_size()
        self.trace_forward(idx, a, no_apertures=True,
                           extra_distance=extra_distance)


text:
"""
Exercise: Linear Regression using Adagrad

@auth: Yu-Hsiang Fu
@date: 2018/09/18
"""
# --------------------------------------------------------------------------------
# 1.Import packages
# --------------------------------------------------------------------------------
import copy
import matplotlib.pyplot as plt
import numpy as np
import pickle
import sys

# --------------------------------------------------------------------------------
# 2.Const variables
# --------------------------------------------------------------------------------
# program variable
# EPSILON = 0.00000001
# ROUND_DIGITS = 8

# plot variable
PLOT_X_SIZE = 5
PLOT_Y_SIZE = 5
PLOT_DPI = 300
PLOT_FORMAT = "png"

# --------------------------------------------------------------------------------
# 3.Define function
# --------------------------------------------------------------------------------
def p_deepcopy(data, dumps=pickle.dumps, loads=pickle.loads):
    return loads(dumps(data, -1))

def f_theta(x, theta):
    return theta[0] + (x * theta[1])

def error(x, y):
    return 0.5 * (np.sum((y - x) ** 2))

def normalize_data(d):
    mean = np.mean(d)
    sigma = np.std(d)

    return (d - mean) / sigma

def gradient_descent(x, y, theta, rate_learning):
    counter_iter = 1
    error_diff = sys.maxsize
    error_prev = sys.maxsize
    error_stop = 0.0001
    max_iteration = 1000
    theta_history = [p_deepcopy(theta)]

    # DO gradient-descent
    theta0_grad = 0.0
    theta1_grad = 0.0

    while error_diff > error_stop:
        # calculate gradient
        # use Adagrad, the rate_learning should be at least larger, e.g., rate_learning = 1.0
        fx = f_theta(x, theta)
        grad0 = np.sum(fx - y)
        grad1 = np.sum((fx - y) * x)
        theta0_grad += np.sqrt(np.power(grad0, 2))
        theta1_grad += np.sqrt(np.power(grad1, 2))

        # update theta
        theta0 = theta[0] - (rate_learning / np.sqrt(theta0_grad)) * grad0
        theta1 = theta[1] - (rate_learning / np.sqrt(theta1_grad)) * grad1
        theta = [theta0, theta1]

        # update error_diff
        error_current = error(fx, y)
        error_diff = error_prev - error_current
        error_prev = copy.copy(error_current)

        # store updated-results
        theta_history.append(p_deepcopy(theta))

        # stop by max_iteration
        if counter_iter < max_iteration:
            pass
        else:
            break

        print(counter_iter, theta, error_current, error_diff, error_diff > error_stop)
        counter_iter += 1

    return theta_history

def draw_plot(x, y, theta_history):
    # contour-plot
    # plot-xy-range
    x_left, x_right = -100, 900
    y_left, y_right = -200, 400

    # plot-grid
    num_level = 50
    x_ls = np.linspace(x_left, x_right, num_level)
    y_ls = np.linspace(y_left, y_right, num_level)

    X_ls, Y_ls = np.meshgrid(x_ls, y_ls)
    Z_ls = np.zeros((len(x_ls), len(y_ls)))

    for i in range(len(x_ls)):
        for j in range(len(y_ls)):

            for k in range(len(x)):
                Z_ls[j][i] += pow(y[k] - (x_ls[i] + y_ls[j] * x[k]), 2)

            Z_ls[j][i] /= len(y_ls)

    # create a figure
    fig, ax = plt.subplots(figsize=(PLOT_X_SIZE, PLOT_Y_SIZE), facecolor="w")

    # draw plot
    theta0 = [theta0 for theta0, theta1, in theta_history]
    theta1 = [theta1 for theta0, theta1, in theta_history]
    ax.plot(theta0, theta1, "o-", ms=7, lw=1.2, mew=1.2, markevery=2, fillstyle="none", c="b")

    # plot setting
    ax.contourf(X_ls, Y_ls, Z_ls, num_level, alpha=0.5, cmap=plt.get_cmap("jet"))
    ax.grid(color="k", linestyle="dotted", linewidth=0.8, alpha=0.8)
    ax.set_xlabel(r"$\theta_0$", fontdict={"fontsize": 12})
    ax.set_ylabel(r"$\theta_1$", fontdict={"fontsize": 12})
    ax.set_xlim(x_left, x_right)
    ax.set_ylim(y_left, y_right)
    ax.set_title("Adagrad, iter={0}".format(len(theta_history)), fontdict={"fontsize": 12})
    ax.tick_params(axis="both", direction="in", which="major", labelsize=8)

    # save the figure
    plt.tight_layout()
    plt.savefig("gradient-descent_Adagrad.png", dpi=PLOT_DPI, format=PLOT_FORMAT, bbox_inches='tight', pad_inches=0.05)
    plt.close()

# --------------------------------------------------------------------------------
# 4.Main function
# --------------------------------------------------------------------------------
def main_function():
    # input data
    d = np.loadtxt("data.csv", delimiter=",")

    # input, output
    x = d[:, 0]
    y = d[:, 1]

    # normalize data
    x = normalize_data(x)

    # --------------------------------------------------
    # linear-regression setting (model, parameters)
    # model: f_theta = theta_0 + (theta_1 * x)
    # parameters: theta = (theta_0, theta_1)
    theta = np.random.rand(2)

    # gradient_descent, input: x, y, theta; output: theta_history
    rate_learning = 1
    theta_history = gradient_descent(x, y, theta, rate_learning)

    # --------------------------------------------------
    # draw scatter-line plot and contour-plot
    draw_plot(x, y, theta_history)

if __name__ == '__main__':
    main_function()

ext: py

max_stars_repo_name: yuhsiangfu/regression

max_issues_repo_name: yuhsiangfu/regression

max_forks_repo_name: yuhsiangfu/regression

reason: import numpy


text:
# Minhyuk Sung (mhsung@cs.stanford.edu)
# January 2017

import numpy as np
import random

class Walkerrandom:
    """ Walker's alias method for random objects with different probablities
    http://code.activestate.com/recipes/576564-walkers-alias-method-for-random-objects-with-diffe/
    """

    def __init__(self, weights, keys=None):
        """ builds the Walker tables prob and inx for calls to random().
        The weights (a list or tuple or iterable) can be in any order;
        they need not sum to 1.
        """

        n = self.n = len(weights)
        self.keys = keys
        sumw = sum(weights)
        prob = [w * n / sumw for w in weights]  # av 1
        inx = [-1] * n
        short = [j for j, p in enumerate(prob) if p < 1]
        long = [j for j, p in enumerate(prob) if p > 1]
        while short and long:
            j = short.pop()
            k = long[-1]
            # assert prob[j] <= 1 <= prob[k]
            inx[j] = k
            prob[k] -= (1 - prob[j])  # -= residual weight
            if prob[k] < 1:
                short.append(k)
                long.pop()
            #if Test:
            #    print "test Walkerrandom: j k pk: %d %d %.2g" % (j, k, prob[k])
        self.prob = prob
        self.inx = inx
        #if Test:
        #    print "test", self

    def __str__(self):
        """ e.g. "Walkerrandom prob: 0.4 0.8 1 0.8  inx: 3 3 -1 2" """

        probstr = " ".join(["%.2g" % x for x in self.prob])
        inxstr = " ".join(["%.2g" % x for x in self.inx])
        return "Walkerrandom prob: %s  inx: %s" % (probstr, inxstr)

    def random(self):
        """ each call -> a random int or key with the given probability
        fast: 1 randint(), 1 random.uniform(), table lookup
        """

        u = random.uniform(0, 1)
        j = random.randint(0, self.n - 1)  # or low bits of u
        randint = j if u <= self.prob[j] else self.inx[j]
        return self.keys[randint] if self.keys else randint

def resample_points(points, scales, n_out_points=None):
    assert (points.ndim == 3)
    assert (scales.ndim == 1)
    assert (points.shape[0] == scales.size)

    n_sets = points.shape[0]
    if n_sets == 1 and n_out_points is None:
        return points[0]

    n_in_points = points.shape[1]
    n_dim = points.shape[2]
    if n_out_points is None:
        n_out_points = points.shape[1]

    wrand = Walkerrandom(scales)
    set_idxs = np.empty(n_out_points, dtype='int')
    for i in range(n_out_points):
        set_idxs[i] = wrand.random()

    point_rand_idxs_in_sets = np.random.randint(n_in_points, size=n_out_points)

    sample_points = points[set_idxs, point_rand_idxs_in_sets, :]
    return sample_points

def axis_aligned_bbox_center(points):
    assert (points.ndim == 2)
    bb_min = np.nanmin(points, axis=0)
    bb_max = np.nanmax(points, axis=0)
    bb_center = 0.5 * (bb_min + bb_max)
    return bb_center

def centerize_points(points):
    assert (points.ndim == 2)
    bb_center = axis_aligned_bbox_center(points)
    centered_points = points - np.expand_dims(bb_center, 0)
    return centered_points, bb_center

max_stars_repo_name: mhsung/fuzzy-set-dual

max_issues_repo_name: mhsung/fuzzy-set-dual

max_forks_repo_name: mhsung/fuzzy-set-dual


text:
"""
Base class for a digital multimeter.
"""
import numpy as np
from msl.equipment import Backend
from msl.qt import (
    Signal,
    MICRO,
)

from . import BaseEquipment
from ..utils import ave_std

class DMM(BaseEquipment):

    FUNCTIONS = {
        'DCV': 'VOLTAGE',
        'VOLT': 'VOLTAGE',
        'VOLTAGE': 'VOLTAGE',
        '"VOLT:DC"': 'VOLTAGE',
        'CURR': 'CURRENT',
        'CURRENT': 'CURRENT',
        'DCI': 'CURRENT',
        '"CURR:DC"': 'CURRENT',
    }

    RANGES = {
        'AUTO': 'AUTO',
        'auto': 'AUTO',
        'DEF': 'DEFAULT',
        'def': 'DEFAULT',
        'DEFAULT': 'DEFAULT',
        'default': 'DEFAULT',
        'MAX': 'MAXIMUM',
        'max': 'MAXIMUM',
        'MAXIMUM': 'MAXIMUM',
        'maximum': 'MAXIMUM',
        'MIN': 'MINIMUM',
        'min': 'MINIMUM',
        'MINIMUM': 'MINIMUM',
        'minimum': 'MINIMUM',
        '100 mV': 0.1,
        '1 V': 1.0,
        '10 V': 10.0,
        '100 V': 100.0,
        '1000 V': 1000.0,
        f'1 {MICRO}A': 1e-6,
        f'10 {MICRO}A': 10e-6,
        f'100 {MICRO}A': 100e-6,
        '1 uA': 1e-6,
        '10 uA': 10e-6,
        '100 uA': 100e-6,
        '1 mA': 0.001,
        '10 mA': 0.01,
        '100 mA': 0.1,
        '1 A': 1.0,
        '3 A': 3.0,
        '10 A': 10.0,
    }

    NPLCS = {
        100.0: 100,
        10.0: 10,
        1.0: 1,
        0.2: 0.2,
        0.06: 0.06,
        0.02: 0.02,
    }

    EDGES = {
        'RISING': 'POSITIVE',
        'POSITIVE': 'POSITIVE',
        'POS': 'POSITIVE',
        'FALLING': 'NEGATIVE',
        'NEGATIVE': 'NEGATIVE',
        'NEG': 'NEGATIVE',
    }

    TRIGGERS = {
        'IMM': 'IMMEDIATE',
        'IMMEDIATE': 'IMMEDIATE',
        'EXT': 'EXTERNAL',
        'EXTERNAL': 'EXTERNAL',
        'BUS': 'BUS',
        'INT': 'INTERNAL',
        'INTERNAL': 'INTERNAL',
    }

    AUTO = {
        'OFF': 'OFF',
        'off': 'OFF',
        '0': 'OFF',
        False: 'OFF',
        0: 'OFF',
        'ON': 'ON',
        'on': 'ON',
        '1': 'ON',
        True: 'ON',
        1: 'ON',
        'ONCE': 'ONCE',
        'once': 'ONCE',
        '2': 'ONCE',
        2: 'ONCE',
    }

    fetched = Signal(float, float)  # average, stdev
    config_changed = Signal(dict)  # the result of info()

    def __init__(self, app, record, *, demo=None):
        """Base class for a digital multimeter.

        Parameters
        ----------
        app : :class:`photons.App`
            The main application entry point.
        record : :class:`~msl.equipment.record_types.EquipmentRecord`
            The equipment record.
        demo : :class:`bool`, optional
            Whether to simulate a connection to the equipment by opening
            a connection in demo mode.
        """
        if record.connection.backend == Backend.PyVISA:
            # `pop` this item from the properties to avoid the following if using PyVISA for the connection
            #   ValueError: 'clear' is not a valid attribute for type *Instrument
            #   ValueError: 'reset' is not a valid attribute for type *Instrument
            pop_or_get = record.connection.properties.pop
        else:
            # use `get` if the backend is MSL so that clear/reset remains in
            # the properties if it has been defined
            pop_or_get = record.connection.properties.get

        clear = pop_or_get('clear', False)
        reset = pop_or_get('reset', False)

        super(DMM, self).__init__(app, record, demo=demo)

        # suppress the warning that the following attributes cannot be made
        # available when starting the BaseEquipment as a Service
        self.ignore_attributes(['fetched', 'config_changed'])

        if clear:
            self.clear()
        if reset:
            self.reset()

    def reset(self) -> None:
        """Resets device to factory default state."""
        self.logger.info(f'reset {self.alias!r}')
        self._send_command_with_opc('*RST')

    def clear(self) -> None:
        """Clears the event registers in all register groups and the error queue."""
        self.logger.info(f'clear {self.alias!r}')
        self._send_command_with_opc('*CLS')

    def bus_trigger(self) -> None:
        """Send a software trigger."""
        self.logger.info(f'software trigger {self.alias!r}')
        self.write('INIT;*TRG')

    software_trigger = bus_trigger

    def fetch(self, initiate: bool = False) -> tuple:
        """Fetch the samples.

        Parameters
        ----------
        initiate : :class:`bool`
            Whether to send the ``'INIT'`` command before sending ``'FETCH?'``.

        Returns
        -------
        :class:`float`
           The average value.
        :class:`float`
           The standard deviation.
        """
        if initiate:
            self.logger.info(f'send INIT to {self.alias!r}')
            self.connection.write('INIT')
        samples = self.connection.query('FETCH?').rstrip()
        return self._average_and_emit(samples)

    def acquisition_time(self, *,
                         info: dict = None,
                         line_freq: float = 50.) -> float:
        """Get the approximate number of seconds it takes to acquire the data.

        Parameters
        ----------
        info : :class:`dict`, optional
            The value returned by :meth:`.info`. If not specified then it
            will be automatically retrieved.
        line_freq : :class:`float`, optional
            The line frequency, in Hz.

        Returns
        -------
        :class:`float`
            The number of seconds.
        """
        if not info:
            info = self.info()
        duration = info['nsamples'] * info['trigger_count'] * info['nplc'] / line_freq
        if info['auto_zero'] == 'ON':
            duration *= 2.
        return duration

    def check_errors(self):
        """Query the digital multimeter’s error queue.

        If there is an error then raise an exception.
        """
        raise NotImplementedError

    def info(self) -> dict:
        """Get the configuration information of the digital multimeter."""
        raise NotImplementedError

    def configure(self, *, function='voltage', range=10, nsamples=10, nplc=10, auto_zero=True,
                  trigger='bus', edge='falling', ntriggers=1, delay=None) -> dict:
        """Configure the digital multimeter.

        Parameters
        ----------
        function : :class:`str`, optional
            The function to measure. Can be any key in :attr:`.DMM.FUNCTIONS` (case insensitive).
        range : :class:`float` or :class:`str`, optional
            The range to use for the measurement. Can be any key in :attr:`.DMM.RANGES`.
        nsamples : :class:`int`, optional
            The number of samples to acquire after receiving a trigger.
        nplc : :class:`float`, optional
            The number of power line cycles.
        auto_zero : :class:`bool` or :class:`str`, optional
            The auto-zero mode. Can be any key in :attr:`.DMM.AUTO_ZEROS`.
        trigger : :class:`str`, optional
            The trigger mode. Can be any key in :attr:`.DMM.TRIGGERS` (case insensitive).
        edge : :class:`str` or :attr:`.DMM.TriggerEdge`, optional
            The edge to trigger on. Can be any key in :attr:`.DMM.EDGES` (case insensitive).
        ntriggers : :class:`int`, optional
            The number of triggers that are accepted by the digital multimeter
            before returning to the *wait-for-trigger* state.
        delay : :class:`float` or :data:`None`, optional
            The trigger delay in seconds. If :data:`None` then enables the auto-delay
            feature where the digital multimeter automatically determines the delay
            based on the function, range and NPLC.

        Returns
        -------
        :class:`dict`
            The result of :meth:`.info` after the settings have been written.
        """
        raise NotImplementedError

    def _send_command_with_opc(self, command: str):
        """Appends ``'*OPC?'`` to the end of a command.

        The *OPC? command guarantees that command's that were previously sent
        to the device have completed.
        """
        command += ';*OPC?'
        assert self.connection.query(command).startswith('1'), f'{command!r} did not return "1"'

    def _average_and_emit(self, samples) -> tuple:
        """Compute the average and emit the value.

        Parameters
        ----------
        samples : :class:`str` or :class:`list`
            A comma-separated string of readings or a list of readings.

        Returns
        -------
        :class:`float`
           The average value.
        :class:`float`
           The standard deviation.
        """
        if isinstance(samples, str):
            samples = samples.split(',')
        samples = [float(val) for val in samples]
        self.logger.info(f'fetch {self.alias!r} {samples}')
        ave, stdev = ave_std(np.asarray(samples))
        if abs(ave) > 1e30:  # overload
            ave, stdev = np.NaN, np.NaN
        self.fetched.emit(ave, stdev)
        self.emit_notification(ave, stdev)  # emit to all linked Clients
        return ave, stdev

max_stars_repo_name: MSLNZ/pr-single-photons

max_issues_repo_name: MSLNZ/pr-single-photons

max_forks_repo_name: MSLNZ/pr-single-photons


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

# In[12]:

import torch, os
import numpy as np
import matplotlib.pyplot as plt
from my_multi_main4 import Net,test
from torchvision import datasets, transforms

# In[13]:

model_file = 'mnist4_multi_97.4_0.0100_SGD_sparsemax_1_1.00_10.00_2.pt'
params = model_file.strip('.pt').split('_')
_,_,acc,lr,optim_type,max_type,norm_flag,lam,gamma,head_cnt =     params[0],params[1],float(params[2]),float(params[3]),params[4],    params[5],bool(int(params[6])),float(params[7]),float(params[8]),int(params[9])
model = Net(max_type,norm_flag,lam,gamma,head_cnt)
model.load_state_dict(torch.load(model_file,map_location='cpu'))

# In[14]:

model.eval()
use_cuda = torch.cuda.is_available()
use_cuda = False
device = torch.device("cuda" if use_cuda else "cpu")
model = model.to(device)
kwargs = {'num_workers': 1, 'pin_memory': True} if use_cuda else {}
test_loader = torch.utils.data.DataLoader(
        datasets.MNIST('../data', download=True, train=False, transform=transforms.Compose([
                           transforms.ToTensor(),
                           transforms.Normalize((0.1307,), (0.3081,))
                       ])),
        batch_size=1000, shuffle=True, **kwargs)
#test_acc = test(None, model, device, test_loader)
#print(test_acc)

# In[15]:

print(type(test_loader.dataset.test_data),test_loader.dataset.test_data.size())
test_data = test_loader.dataset.test_data

# In[16]:

import torch.nn as nn
import torch.nn.functional as F
def forward_att(att,x):
    N,C,H,W = x.size()
    x = torch.transpose(x,1,3) #[N,W,H,C]
    proj_x = att.proj_func(x) #[N,W,H,C']
    score_x = att.score_func(x) #[N,W,H,1]
    scores = torch.reshape(score_x,[N,-1,1]) #[N,W*H,1]
    projs = torch.reshape(proj_x,[N,-1,att.output_size]) #[N,W*H,C']
    weights = att.mapping_func(att.score_norm(scores),dim=-2) #[N,W*H,1]
    output = torch.sum(projs * weights,dim=-2) #[N,1]
    return proj_x,score_x,weights.reshape(score_x.size()).transpose(1,2),output
def forward(model,x):
    h1 = F.relu(model.conv1(x))
    h1_pooled = F.max_pool2d(h1, 2, 2)
    h2 = F.relu(model.conv2(h1_pooled))
    # h2_pooled = F.max_pool2d(h2, 2, 2)
    att_out = torch.cat([att(h2) for att in model.attention_pool],dim=1).view(-1,500)
    h3 = F.relu(att_out)
    h4 = model.fc2(h3)
    return h1,h1_pooled,h2,att_out,h3,h4
def weight_image(weight):
    image_size = 28
    recp_size = 8
    weight_size = 11
    stride = int((image_size-recp_size)/(weight_size-1))

    output = np.zeros([image_size,image_size])
    for i in range(weight_size):
        for j in range(weight_size):
            output[i*stride:(i*stride+recp_size),j*stride:(j*stride+recp_size)] += weight[i,j]
    return output
normalization = lambda x: (x-np.min(x))/(np.max(x)-np.min(x))
def weight_on_image(weight,image):
    output = np.zeros(weight.shape+(3,))
    output[:,:,2] = normalization(image)
    output[:,:,1] = normalization(weight)
    output = 1 - output
    return output
norm_image = lambda x: (x.type(torch.float)/255.-0.1307)/0.3081

# In[17]:

with torch.no_grad():
    index = np.random.choice(test_data.size()[0])
    plt.imshow(test_data[index])
    plt.show()
    h1,h1_pooled,h2,att_out,h3,h4 = forward(model,norm_image(test_data[index]).unsqueeze_(0).unsqueeze(0))
    plt.imshow(h1_pooled[0,10])
    plt.show()
    #plt.imshow(h2_pooled[0,10])
    #plt.show()
    plt.imshow(h4)
    plt.show()
    proj0,score0,weights0,att_out0 = forward_att(model.attention_pool[0],h2)
    print(proj0.shape,score0.shape,weights0.shape,att_out0.shape)
    plt.imshow(score0[0,:,:,0])
    plt.show()
    plt.imshow(weights0[0,:,:,0])
    plt.show()
    proj1,score1,weights1,att_out1 = forward_att(model.attention_pool[1],h2)
    print(proj1.shape,score1.shape,weights1.shape,att_out1.shape)
    plt.imshow(score1[0,:,:,0])
    plt.show()
    plt.imshow(weights1[0,:,:,0])
    plt.show()

# In[18]:

samples = 1
sample_per_row = 4
weight_sample_path = './weight_samples4_multi_%.1f'%(acc,)
if not os.path.exists(weight_sample_path):
    os.mkdir(weight_sample_path)
plt.axis('off')
plt.gcf().set_facecolor("white")
for _ in range(300):
    fig=plt.figure(figsize=(20, samples*10), dpi= 80)
    with torch.no_grad():
        for s in range(samples):
            index = np.random.choice(test_data.size()[0])
            h1,h1_pooled,h2,att_out,h3,h4 = forward(model,norm_image(test_data[index]).unsqueeze_(0).unsqueeze(0))
            proj0,score0,weights0,att_out0 = forward_att(model.attention_pool[0],h2)
            proj1,score1,weights1,att_out1 = forward_att(model.attention_pool[1],h2)
            image_weighted0 = weight_image(weights0[0,:,:,0].numpy())
            image_weighted1 = weight_image(weights1[0,:,:,0].numpy())
            image_with_weighted0 = weight_on_image(image_weighted0,test_data[index].numpy())
            image_with_weighted1 = weight_on_image(image_weighted1,test_data[index].numpy())

            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+1)
            plt.imshow(test_data[index])
            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+2)
            plt.imshow(image_with_weighted0)
            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+3)
            plt.imshow(score0[0,:,:,0])
            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+4)
            plt.imshow(weights0[0,:,:,0])
            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+5)
            plt.imshow(torch.norm(h2,dim=1)[0])
            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+6)
            plt.imshow(image_with_weighted1)
            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+7)
            plt.imshow(score1[0,:,:,0])
            plt.subplot(2*np.ceil(samples*4/sample_per_row),sample_per_row,8*s+8)
            plt.imshow(weights1[0,:,:,0])
    plt.savefig(os.path.join(weight_sample_path,'%d.png'%index))
    #plt.show()
    plt.close()

# In[21]:

h1,h1_pooled,h2,att_out,h3,h4 = forward(model,norm_image(test_data).unsqueeze_(1))
proj0,score0,weights0,att_out0 = forward_att(model.attention_pool[0],h2)
proj1,score1,weights1,att_out1 = forward_att(model.attention_pool[1],h2)

# In[22]:

all_weights0 = np.sum(np.abs(weights0[:,:,:,0].detach().numpy()),axis=0)
all_weights1 = np.sum(np.abs(weights1[:,:,:,0].detach().numpy()),axis=0)
print(all_weights0)
print(all_weights1)

# In[23]:

plt.imshow(all_weights0)
print(np.sum(all_weights0))
plt.show()
plt.imshow(all_weights1)
print(np.sum(all_weights1))
plt.show()

# In[17]:

print(model.attention_pool[0].score_func.weight)
print(model.attention_pool[0].score_func.bias)
print(model.attention_pool[1].score_func.weight)
print(model.attention_pool[1].score_func.bias)

# In[31]:

plt.imshow(h2[0,np.random.choice(50)])
plt.show()
plt.imshow(test_data[index])
plt.show()

# In[33]:

h2_trans = torch.transpose(h2,1,3)
scores_trans = model.attention_pool[0].score_func(h2_trans)
scores = torch.transpose(scores_trans,1,3)
plt.imshow(scores[0,0].detach().numpy())
plt.show()
scores_flat = scores.reshape(-1)
weights_flat = model.attention_pool[0].mapping_func(scores_flat,-1)
weights = weights_flat.reshape(scores.size())
plt.imshow(weights[0,0].detach().numpy())
plt.show()

# In[ ]:


text:
import math
import time
import torch
import numpy as np
from tqdm import tqdm
from pruning.function.helper import test

def load_sparse_model(net, path, fc_bits):
    '''load the model which is saved as sparse matrix

    Args:
        net:  the network object
        path: the path of the pruned model
        bits: the bits of each index in fc layer

    Returns:
        conv_layer_num:     the Number of convolutional layers
        nz_num:             the Number of non-zero value in each layers
        conv_diff:          the sparse index of each convolutional layers
        fc_diff:            the sparse index of each full-connect layers
        conv_value_array:   the sparse value of each convolutional layers
        fc_value_array:     the sparse value of each full-connect layers
    '''
    conv_layer_num = 0
    fc_layer_num = 0
    fin = open(path, 'rb')
    for name, x in net.named_parameters():
        if name.endswith('mask'):
            continue
        if name.startswith('conv'):
            conv_layer_num += 1
        elif name.startswith('fc'):
            fc_layer_num += 1
    nz_num = np.fromfile(fin, dtype=np.uint32, count=conv_layer_num + fc_layer_num)

    conv_diff_num = sum(nz_num[:conv_layer_num])
    conv_diff = np.fromfile(fin, dtype=np.uint8, count=conv_diff_num)

    fc_merge_num = math.floor((sum(nz_num[conv_layer_num:]) + 1) / 2)
    fc_merge_diff = np.fromfile(fin, dtype=np.uint8, count=fc_merge_num)

    conv_value_array = np.fromfile(fin, dtype=np.float32, count=sum(nz_num[:conv_layer_num]))
    fc_value_array = np.fromfile(fin, dtype=np.float32, count=sum(nz_num[conv_layer_num:]))

    # print(nz_num)
    # print(conv_diff.size, conv_diff[-10:])
    # print(len(fc_merge_diff), fc_merge_diff[-10:])
    # print(conv_value_array.size, conv_value_array[-10:])
    # print(fc_value_array.size, fc_value_array[-10:])

    # [  292    17  8213    15 77747    65   818     1]
    # 8537 [3 1 2 0 3 2 0 1 2 4]
    # 39316 [ 17 242  34  50 164  44  26   3   6 128]
    # 8537 [ 0.05500366 -0.0518913  -0.05787839  0.04747333 -0.07086759 -0.07142863
    #  -0.06043605 -0.06711546 -0.0698091  -0.06924898]
    # 78631 [ 0.13233908  0.16305041 -0.171971   -0.1353672   0.16033891 -0.19598335
    #  -0.11460102 -0.32042998 -0.12170218  0.14367148]

    # Split 8 bits index to 4 bits index
    fc_diff = []
    max_bits = (2 ** fc_bits) - 1
    for i in range(len(fc_merge_diff)):
        fc_diff.append(int(fc_merge_diff[i] >> fc_bits))  # first 4 bits
        fc_diff.append(fc_merge_diff[i] & max_bits)  # last 4 bits

    fc_num_sum = nz_num[conv_layer_num:].sum()
    if fc_num_sum % 2 != 0:
        fc_diff = fc_diff[:fc_num_sum]
    fc_diff = np.asarray(fc_diff, dtype=np.uint8)
    # print("if_error_more_15", (fc_diff > 15).sum())
    # print("if_error_less_0", (fc_diff < 0).sum())
    # print("fc_diff", (fc_diff).sum())
    # layer_index = fc_diff[0:0 + nz_num[4]]
    # print(sum(layer_index) + len(layer_index))

    # print(nz_num)
    # print(conv_diff.size, conv_diff[-10:])
    # print(len(fc_diff), fc_diff[-10:])
    # print(conv_value_array.size, conv_value_array[-10:])
    # print(fc_value_array.size, fc_value_array[-10:])

    # [  292    17  8213    15 77747    65   818     1]
    # 8537 [3 1 2 0 3 2 0 1 2 4]
    # 78631 [ 4  2 12  1 10  0  3  0  6  8]
    # 8537 [ 0.05500366 -0.0518913  -0.05787839  0.04747333 -0.07086759 -0.07142863
    #  -0.06043605 -0.06711546 -0.0698091  -0.06924898]
    # 78631 [ 0.13233908  0.16305041 -0.171971   -0.1353672   0.16033891 -0.19598335
    #  -0.11460102 -0.32042998 -0.12170218  0.14367148]

    return conv_layer_num, nz_num, conv_diff, fc_diff, conv_value_array, fc_value_array

def sparse_to_init(net, before_path, prune_fc_bits):
    conv_layer_length, nz_num, sparse_conv_diff, sparse_fc_diff, conv_value_array, fc_value_array \
        = load_sparse_model(net, before_path, prune_fc_bits)
    state_dict = net.state_dict()
    conv_layer_index = 0
    fc_layer_index = 0
    for i, (key, value) in enumerate(state_dict.items()):
        shape = value.shape
        value = value.view(-1)
        value.zero_()
        if i < conv_layer_length:
            layer_diff = sparse_conv_diff[conv_layer_index:conv_layer_index + nz_num[i]]
            layer_value = conv_value_array[conv_layer_index:conv_layer_index + nz_num[i]]
            conv_layer_index += nz_num[i]
        else:
            layer_diff = sparse_fc_diff[fc_layer_index:fc_layer_index + nz_num[i]]
            layer_value = fc_value_array[fc_layer_index:fc_layer_index + nz_num[i]]
            fc_layer_index += nz_num[i]
        dense_index = 0
        sparse_index = 0
        while sparse_index < len(layer_diff):
            dense_index += layer_diff[sparse_index]
            tmp = layer_value[sparse_index].item()
            value[dense_index] = tmp
            sparse_index += 1
            dense_index += 1
        value.reshape(shape)

def restructure_index(index_list, conv_layer_num, max_conv_bit, max_fc_bit):
    '''load the model which is saved as sparse matrix

    Args:
        index_list:         the index of the codebook
        conv_layer_num:     the Number of convolutional layers
        max_conv_bit:       the bits of each value in convolutional layer
        max_fc_bit:         the bits of each value in full-connect layer

    Returns:
        new_index_list:     Contains the index belonging to each value in codebook
        key_parameter:
    '''
    new_index_list = []

    for i in range(len(index_list)):
        num = max_conv_bit if i < conv_layer_num else max_fc_bit
        tmp_index = []
        for j in range(num):
            tmp_index.append(np.where(np.asarray(index_list[i]) == j)[0].tolist())
        new_index_list.append(tmp_index)

    key_parameter = []
    for j in range(int(len(index_list) / 2)):
        layer_index = np.concatenate((index_list[2 * j], index_list[2 * j + 1]))
        num = max_conv_bit if j < (conv_layer_num / 2) else max_fc_bit
        empty_parameter = [None] * num
        key_parameter.append(empty_parameter)
        for m in range(len(layer_index)):
            # print(m, layer_index[m])
            if layer_index[m] != -1 and key_parameter[j][layer_index[m]] is None:
                key_parameter[j][layer_index[m]] = m
    return new_index_list, key_parameter

def codebook_to_init(net, conv_layer_length, nz_num, sparse_conv_diff, sparse_fc_diff, codebook, max_conv_bit,
                     max_fc_bit):
    state_dict = net.state_dict()
    index_list = []
    conv_layer_index = 0
    fc_layer_index = 0
    for i, (key, value) in enumerate(state_dict.items()):
        shape = value.shape
        value = value.view(-1)

        index = np.empty_like(value, dtype=np.int16)
        index[:] = -1
        value.zero_()
        if i < conv_layer_length:
            layer_diff = sparse_conv_diff[conv_layer_index:conv_layer_index + nz_num[i]]
            conv_layer_index += nz_num[i]
        else:
            layer_diff = sparse_fc_diff[fc_layer_index:fc_layer_index + nz_num[i]]
            fc_layer_index += nz_num[i]
        dense_index = 0
        sparse_index = 0
        half_index = int(i / 2)
        codebook_index_array = codebook.codebook_index[half_index]
        while sparse_index < len(layer_diff):
            dense_index += layer_diff[sparse_index]
            tmp = codebook.codebook_value[half_index][codebook_index_array[sparse_index]].item()
            value[dense_index] = tmp
            index[dense_index] = int(codebook_index_array[sparse_index])
            sparse_index += 1
            dense_index += 1
        value.reshape(shape)
        index.reshape(shape)
        index_list.append(index)

    new_index_list, key_parameter = restructure_index(index_list, conv_layer_length, max_conv_bit,
                                                                  max_fc_bit)
    return new_index_list, key_parameter

def cluster_grad(net, index_list, max_conv_bit, max_fc_bit, conv_layer_length):
    params = list(net.parameters())
    # print('================')
    for i in range(0, len(params), 2):
        param = params[i]
        grad_shape = param.grad.shape
        grad = param.grad
        grad = grad.view(-1)
        index = index_list[i]

        bias = params[i + 1]
        bias_grad_shape = bias.grad.shape
        bias_grad = bias.grad
        bias_grad = bias_grad.view(-1)
        bias_index = index_list[i + 1]

        # start = time.clock()
        # Cluster grad using index, use mean of each class of grad to update weight
        cluster_bits = max_conv_bit if i < conv_layer_length else max_fc_bit
        for j in range(cluster_bits):
            sum_grad = grad[index[j]].sum()
            sum_grad += bias_grad[bias_index[j]].sum()
            grad[index[j]] = sum_grad
            bias_grad[bias_index[j]] = sum_grad
        # end = time.clock()
        # print(round(end - start, 5))

        grad = grad.view(grad_shape)
        params[i].grad = grad.clone()

        bias_grad = bias_grad.view(bias_grad_shape)
        params[i + 1].grad = bias_grad.clone()

def update_codebook(net, codebook, conv_layer_length, max_conv_bit, max_fc_bit, key_parameter):
    params = list(net.parameters())
    for i in range(0, len(params), 2):
        # start = time.clock()
        param = params[i]
        param = param.view(-1)
        bias_param = params[i + 1]
        bias_param = bias_param.view(-1)
        layer = torch.cat((param, bias_param))
        half_index = int(i / 2)
        cluster_bits = max_conv_bit if i < conv_layer_length else max_fc_bit
        codebook_centroids = codebook.codebook_value[half_index]
        for j in range(cluster_bits):
            if key_parameter[half_index][j] is not None:
                tmp = key_parameter[half_index][j]
                codebook_centroids[j] = layer[tmp]

def train_codebook(max_accuracy, nz_num, conv_diff, fc_diff, retrain_codebook_path, key_parameter, use_cuda, max_conv_bit, max_fc_bit, conv_layer_length,
                   codebook, index_list, testloader, net, trainloader, criterion, optimizer,
                   scheduler, epoch=1, top_5=False):
    for epoch in range(epoch):  # loop over the dataset multiple times
        train_loss = []
        net.train()
        # i = 0
        # for inputs, labels in tqdm(trainloader):
        for inputs, labels in trainloader:
            # get the inputs
            if use_cuda:
                inputs = inputs.cuda()
                labels = labels.cuda()

            # zero the parameter gradients
            optimizer.zero_grad()
            # forward + backward + optimize
            outputs = net(inputs)  # forward
            loss = criterion(outputs, labels)  # compute loss
            loss.backward()  # backward

            cluster_grad(net, index_list, max_conv_bit, max_fc_bit, conv_layer_length)

            optimizer.step()  # update weight
            train_loss.append(loss.item())

        # mean_train_loss = np.mean(train_loss)
        # print("Epoch:", epoch, "Training Loss: %5f" % mean_train_loss)
        accuracy = test(use_cuda, testloader, net, top_5)
        scheduler.step()
    update_codebook(net, codebook, conv_layer_length, max_conv_bit, max_fc_bit, key_parameter)
    if accuracy > max_accuracy:
       save_codebook(conv_layer_length, nz_num, conv_diff, fc_diff, codebook, retrain_codebook_path, net)
       max_accuracy = accuracy

def save_codebook(conv_layer_length, nz_num, conv_diff, fc_diff, codebook, path, net):
    fc_merge_diff = []

    # print(nz_num)
    # print(len(conv_diff), conv_diff[-10:])
    # print(len(fc_diff), fc_diff[-10:])
    # [   304     11   5353      1 400000    500   5000     10]
    # 5669 [ 0  2  0  1  1  1  0  9  8 44]
    # 405510 [0 0 0 0 0 0 0 0 0 0]

    length = len(fc_diff)
    fc_diff = list(fc_diff)
    if length % 2 != 0:
        fc_diff.append(0)
    for i in range(math.floor(len(fc_diff) / 2)):
        fc_merge_diff.append((fc_diff[2 * i] << 4) | fc_diff[2 * i + 1])
    nz_num = np.asarray(nz_num, dtype=np.uint32)
    conv_diff = np.asarray(conv_diff, dtype=np.uint8)
    fc_merge_diff = np.asarray(fc_merge_diff, dtype=np.uint8)

    conv_half_len = int(conv_layer_length / 2)
    conv_codebook_index = []
    for m in range(conv_half_len):
        conv_codebook_index.extend(codebook.codebook_index[m])

    fc_codebook_index = []
    for k in range(conv_half_len, len(codebook.codebook_index)):
        fc_codebook_index.extend(codebook.codebook_index[k])

    codebook_value = []
    for j in range(len(codebook.codebook_value)):
        codebook_value.extend(codebook.codebook_value[j])

    # print(len(conv_codebook_index), conv_codebook_index[-10:])
    # print(len(fc_codebook_index), fc_codebook_index[-10:])
    # print(len(codebook_value), codebook_value[-10:])
    # 5669 [2, 228, 211, 229, 76, 152, 23, 116, 111, 25]
    # 405510 [10, 11, 5, 6, 9, 7, 5, 7, 12, 5]
    # 544 [-0.11808116, -0.06328904, 0.1446653, 0.051914066, -0.03960273, -0.017428499, -0.017428499, 0.0050489083, 0.22879101, 0.051914066]

    length = len(fc_codebook_index)
    if length % 2 != 0:
        fc_codebook_index.append(0)

    fc_codebook_index = np.asarray(fc_codebook_index, dtype=np.uint8)
    fc_codebook_index_merge = []
    for i in range(math.floor((len(fc_codebook_index)) / 2)):
        fc_codebook_index_merge.append(
            (fc_codebook_index[2 * i] << 4) | fc_codebook_index[2 * i + 1])

    conv_codebook_index = np.asarray(conv_codebook_index, dtype=np.uint8)
    fc_codebook_index_merge = np.asarray(fc_codebook_index_merge, dtype=np.uint8)
    codebook_value = np.asarray(codebook_value, dtype=np.float32)

    # print(any(np.isnan(codebook_value)))

    # print(nz_num)
    # print(len(conv_diff), conv_diff[-10:])
    # print(len(fc_merge_diff), fc_merge_diff[-10:])
    # print(len(conv_codebook_index), conv_codebook_index[-10:])
    # print(len(fc_codebook_index_merge), fc_codebook_index_merge[-10:])
    # print(len(codebook_value), codebook_value[-10:])
    # [   304     11   5353      1 400000    500   5000     10]
    # 5669 [ 0  2  0  1  1  1  0  9  8 44]
    # 202755 [0 0 0 0 0 0 0 0 0 0]
    # 5669 [  2 228 211 229  76 152  23 116 111  25]
    # 202755 [200  66  71 152 140 171  86 151  87 197]
    # 544 [-0.11808116 -0.06328904  0.1446653   0.05191407 -0.03960273 -0.0174285
    #  -0.0174285   0.00504891  0.22879101  0.05191407]

    # Set to the same dtype uint8 to save
    nz_num.dtype = np.uint8
    codebook_value.dtype = np.uint8

    sparse_obj = np.concatenate((nz_num, conv_diff, fc_merge_diff, conv_codebook_index,
                                 fc_codebook_index_merge, codebook_value))
    sparse_obj.tofile(path)

max_stars_repo_name: edmontdants/DeepCompression-PyTorch

max_issues_repo_name: yangbingjie/DeepCompression-PyTorch

max_forks_repo_name: yangbingjie/DeepCompression-PyTorch


text:
from itertools import combinations

import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib.patches import Rectangle

import seaborn as sns
import numpy as np
import pandas as pd
import statsmodels.api as sm

from metroplot import metroplot

# Metroplot usage example

# Read Eysenck (1974) table from David Howell's website
data = pd.read_csv(
    "https://www.uvm.edu/~statdhtx/methods8/DataFiles/Tab13-2.dat", delimiter="\t"
)
data["Condition"] = [
    {1: "Counting", 2: "Rhyming", 3: "Adjective", 4: "Imagery", 5: "Intention"}[c]
    for c in data["Condition"]
]
data["Age"] = [{1: "Old", 2: "Young"}[c] for c in data["Age"]]

def pairwise_ttest(data, factor, dependent):
    """run t-tests between all levels of a factor and"""
    levels = data[factor].unique()

    # test all pairwise comparisons
    pairwise_comparisons = []
    for level1, level2 in combinations(levels, r=2):
        a = data.loc[data[factor] == level1, dependent]
        b = data.loc[data[factor] == level2, dependent]

        t, p_value, df = sm.stats.ttest_ind(a, b)

        pairwise_comparisons.append(
            {
                "level1": level1,
                "level2": level2,
                "effect_direction": np.sign(t),
                "p_value": p_value,
            }
        )
    pairwise_comparisons = pd.DataFrame(pairwise_comparisons)

    # multiple comparison correction
    (
        is_sig,
        corrected,
        _,
        _,
    ) = sm.stats.multipletests(pairwise_comparisons["p_value"], method="fdr_bh")

    pairwise_comparisons["is_sig"] = is_sig
    pairwise_comparisons["corrected_p_value"] = corrected  # this is not currently used

    return pairwise_comparisons

pairwise_comparisons = pairwise_ttest(data, "Condition", "Recall")

conditions = data["Condition"].unique()

print("pairwise comparisons results:")
print(pairwise_comparisons)

fig = plt.figure(figsize=(5, 4))
gs = GridSpec(
    2, 1, height_ratios=[0.75, 5], hspace=0.05, figure=fig
)  # we create two subplots, one for the main plot and one for the metroplot.
ax0 = fig.add_subplot(gs[1])
sns.barplot(data=data, x="Condition", y="Recall", order=conditions)
sns.despine(fig)
gs.tight_layout(fig)

# # grab colors from Seaborn's barplot.
bars = [r for r in ax0.get_children() if type(r) == Rectangle]
colors = [c.get_facecolor() for c in bars[:-1]]
level_palette = {cond: color for cond, color in zip(conditions, colors)}

level_to_location = {
    cond: i for i, cond in enumerate(conditions)
}  # map categories to x axis locations.
ax1 = fig.add_subplot(gs[0])
metroplot(
    pairwise_comparisons,
    level_to_location=level_to_location,
    metroplot_element_order=conditions,
    ax=ax1,
    level_axis="x",
    dominating_effect_direction=1,
    level_pallete=level_palette,
    level_axis_lim=ax0.get_xlim(),
)
# # note that level_axis_lim=ax0.get_xlim() aligns the metrplot with the main plot

# %%
# a slightly more complicated example

fig = plt.figure(figsize=(6, 6))
gs = GridSpec(1, 2, width_ratios=[4, 2.5], wspace=0.01, figure=fig)
ax0 = fig.add_subplot(gs[0])

data["Group"] = [cond + " " + age for age, cond in zip(data["Age"], data["Condition"])]
level_palette = {
    "Counting": "#a6cee3",
    "Rhyming": "#1f78b4",
    "Adjective": "#b2df8a",
    "Imagery": "#33a02c",
    "Intention": "#fb9a99",
}

pairwise_comparisons = pairwise_ttest(data, "Group", "Recall")

print("pairwise comparisons results:")
print(pairwise_comparisons)

conditions = data["Group"].unique()
conditions.sort()

sns.boxplot(
    data=data, y="Group", x="Recall", ax=ax0, color="w", order=conditions
)  # , hue='Condition', palette=level_palette)
gs.tight_layout(fig)

level_to_location = {
    cond: i for i, cond in enumerate(conditions)
}  # map categories to y axis locations.

ax1 = fig.add_subplot(gs[1])
metroplot(
    pairwise_comparisons,
    level_to_location=level_to_location,
    metroplot_element_order=conditions,
    ax=ax1,
    level_axis="y",
    dominating_effect_direction=1,
    level_pallete="k",
    level_axis_lim=ax0.get_ylim(),
)
# # note that level_axis_lim=ax0.get_ylim() aligns the metrplot with the main plot

plt.show()

max_stars_repo_name: Tal-Golan/metroplot

max_issues_repo_name: Tal-Golan/metroplot

max_forks_repo_name: Tal-Golan/metroplot

reason: import numpy,import statsmodels


text:
"""
Multivariate Wald-Wolfowitz test for two samples in separate CSV files.

See:
  Friedman, Jerome H., and Lawrence C. Rafsky.
  "Multivariate generalizations of the Wald-Wolfowitz and Smirnov two-sample tests."
  The Annals of Statistics (1979): 697-717.

Given multivariate sample X of length m and sample Y of length n, test the null hypothesis:

    H_0: The samples were generated by the same distribution

The algorithm uses a KD-tree to construct the minimum spanning tree, therefore is O(Nk log(Nk))
instead of O(N^3), where N = m + n is the total number of observations. Though approximate,
results are generally valid. See also:

Monaco, John V.
"Classification and authentication of one-dimensional behavioral biometrics."
Biometrics (IJCB), 2014 IEEE International Joint Conference on. IEEE, 2014.

The input files should be CSV files with no header and row observations, for example:

X.csv
-----
1,2
2,2
3,1

Y.csv
-----
1,1
2,4
3,2
4,2

Usage:

$ python X.csv Y.csv

> W = 0.485, 5 runs
> p = 0.6862
> Fail to reject H_0 at 0.05 significance level
> The samples appear to have similar distribution
"""
import sys
import numpy as np
import pandas as pd
import scipy.stats as stats
from sklearn.neighbors import kneighbors_graph
from scipy.sparse.csgraph import minimum_spanning_tree

SIGNIFICANCE = 0.05

def mst_edges(V, k):
    """
    Construct the approximate minimum spanning tree from vectors V

    :param: V: 2D array, sequence of vectors
    :param: k: int the number of neighbor to consider for each vector
    :return: V ndarray of edges forming the MST
    """

    # k = len(X)-1 gives the exact MST
    k = min(len(V) - 1, k)

    # generate a sparse graph using the k nearest neighbors of each point
    G = kneighbors_graph(V, n_neighbors=k, mode='distance')

    # Compute the minimum spanning tree of this graph
    full_tree = minimum_spanning_tree(G, overwrite=True)

    return np.array(full_tree.nonzero()).T

def ww_test(X, Y, k=10):
    """
    Multi-dimensional Wald-Wolfowitz test

    :param X: multivariate sample X as a numpy ndarray
    :param Y: multivariate sample Y as a numpy ndarray
    :param k: number of neighbors to consider for each vector
    :return: W the WW test statistic, R the number of runs
    """
    m, n = len(X), len(Y)
    N = m + n

    XY = np.concatenate([X, Y]).astype(np.float)

    # XY += np.random.normal(0, noise_scale, XY.shape)

    edges = mst_edges(XY, k)

    labels = np.array([0] * m + [1] * n)

    c = labels[edges]
    runs_edges = edges[c[:, 0] == c[:, 1]]

    # number of runs is the total number of observations minus edges within each run
    R = N - len(runs_edges)

    # expected value of R
    e_R = ((2.0 * m * n) / N) + 1

    # variance of R is _numer/_denom
    _numer = 2 * m * n * (2 * m * n - N)
    _denom = N ** 2 * (N - 1)

    # see Eq. 1 in Friedman 1979
    # W approaches a standard normal distribution
    W = (R - e_R) / np.sqrt(_numer/_denom)

    return W, R

if __name__ == '__main__':
    if len(sys.argv) != 3:
        print('Usage: $ python runs_test.py [X] [Y]')
        sys.exit(1)

    X = pd.read_csv(sys.argv[1], header=None).values
    Y = pd.read_csv(sys.argv[2], header=None).values

    W, R = ww_test(X, Y)
    pvalue = stats.norm.cdf(W)  # one sided test
    reject = pvalue <= SIGNIFICANCE

    print('W = %.3f, %d runs' % (W, R))
    print('p = %.4f' % pvalue)
    print('%s H_0 at 0.05 significance level' % ('Reject' if reject else 'Fail to reject'))
    print('The samples appear to have %s distribution' % ('different' if reject else 'similar'))

max_stars_repo_name: jjhenkel/dockerizeme

max_issues_repo_name: jjhenkel/dockerizeme

max_forks_repo_name: jjhenkel/dockerizeme

reason: import numpy,import scipy,from scipy


text:
"""
This module contains the loss classes.

Specific losses are used for regression, binary classification or multiclass
classification.
"""
from abc import ABC, abstractmethod

from scipy.special import expit, logsumexp
import numpy as np
from numba import njit, prange

from .utils import get_threads_chunks

@njit
def _logsumexp(a):
    """logsumexp(x) = log(sum(exp(x)))

    Custom logsumexp function with numerical stability, based on scipy's
    logsumexp which is unfortunately not supported (neither is
    np.logaddexp.reduce, which is equivalent). Only supports 1d arrays.
    """

    a_max = np.amax(a)
    if not np.isfinite(a_max):
        a_max = 0

    s = np.sum(np.exp(a - a_max))
    return np.log(s) + a_max

@njit(fastmath=True)
def _expit(x):
    # custom sigmoid because we cannot use that of scipy with numba
    return 1 / (1 + np.exp(-x))

class BaseLoss(ABC):
    """Base class for a loss."""

    def init_gradients_and_hessians(self, n_samples, prediction_dim):
        """Return initial gradients and hessians.

        Unless hessians are constant, arrays are initialized with undefined
        values.

        Parameters
        ----------
        n_samples : int
            The number of samples passed to `fit()`
        prediction_dim : int
            The dimension of a raw prediction, i.e. the number of trees
            built at each iteration. Equals 1 for regression and binary
            classification, or K where K is the number of classes for
            multiclass classification.

        Returns
        -------
        gradients : array-like, shape=(n_samples * prediction_dim)
        hessians : array-like, shape=(n_samples * prediction_dim).
            If hessians are constant (e.g. for ``LeastSquares`` loss, shape
            is (1,) and the array is initialized to ``1``.
        """
        shape = n_samples * prediction_dim
        gradients = np.empty(shape=shape, dtype=np.float32)
        if self.hessian_is_constant:
            hessians = np.ones(shape=1, dtype=np.float32)
        else:
            hessians = np.empty(shape=shape, dtype=np.float32)

        return gradients, hessians

    @abstractmethod
    def get_baseline_prediction(self, y_train, prediction_dim):
        """Return initial predictions (before the first iteration).

        Parameters
        ----------
        y_train : array-like, shape=(n_samples,)
            The target training values.
        prediction_dim : int
            The dimension of one prediction: 1 for binary classification and
            regression, n_classes for multiclass classification.

        Returns
        -------
        baseline_prediction: float or array of shape (1, prediction_dim)
            The baseline prediction.
        """
        pass

    @abstractmethod
    def update_gradients_and_hessians(self, gradients, hessians, y_true,
                                      raw_predictions):
        """Update gradients and hessians arrays, inplace.

        The gradients (resp. hessians) are the first (resp. second) order
        derivatives of the loss for each sample with respect to the
        predictions of model, evaluated at iteration ``i - 1``.

        Parameters
        ----------
        gradients : array-like, shape=(n_samples * prediction_dim)
            The gradients (treated as OUT array).
        hessians : array-like, shape=(n_samples * prediction_dim) or \
            (1,)
            The hessians (treated as OUT array).
        y_true : array-like, shape=(n_samples,)
            The true target values or each training sample.
        raw_predictions : array-like, shape=(n_samples, prediction_dim)
            The raw_predictions (i.e. values from the trees) of the tree
            ensemble at iteration ``i - 1``.
        """
        pass

class LeastSquares(BaseLoss):
    """Least squares loss, for regression.

    For a given sample x_i, least squares loss is defined as::

        loss(x_i) = (y_true_i - raw_pred_i)**2
    """

    hessian_is_constant = True

    def __call__(self, y_true, raw_predictions, average=True):
        # shape (n_samples, 1) --> (n_samples,). reshape(-1) is more likely to
        # return a view.
        raw_predictions = raw_predictions.reshape(-1)
        loss = np.power(y_true - raw_predictions, 2)
        return loss.mean() if average else loss

    def get_baseline_prediction(self, y_train, prediction_dim):
        return np.mean(y_train)

    def inverse_link_function(self, raw_predictions):
        return raw_predictions

    def update_gradients_and_hessians(self, gradients, hessians, y_true,
                                      raw_predictions):
        return _update_gradients_least_squares(gradients, y_true,
                                               raw_predictions)

@njit(parallel=True, fastmath=True)
def _update_gradients_least_squares(gradients, y_true, raw_predictions):
    # shape (n_samples, 1) --> (n_samples,). reshape(-1) is more likely to
    # return a view.
    raw_predictions = raw_predictions.reshape(-1)
    n_samples = raw_predictions.shape[0]
    starts, ends, n_threads = get_threads_chunks(total_size=n_samples)
    for thread_idx in prange(n_threads):
        for i in range(starts[thread_idx], ends[thread_idx]):
            # Note: a more correct exp is 2 * (raw_predictions - y_true) but
            # since we use 1 for the constant hessian value (and not 2) this
            # is strictly equivalent for the leaves values.
            gradients[i] = raw_predictions[i] - y_true[i]

class BinaryCrossEntropy(BaseLoss):
    """Binary cross-entropy loss, for binary classification.

    For a given sample x_i, the binary cross-entropy loss is defined as the
    negative log-likelihood of the model which can be expressed as::

        loss(x_i) = log(1 + exp(raw_pred_i)) - y_true_i * raw_pred_i

    See The Elements of Statistical Learning, by Hastie, Tibshirani, Friedman.
    """

    hessian_is_constant = False
    inverse_link_function = staticmethod(expit)

    def __call__(self, y_true, raw_predictions, average=True):
        # shape (n_samples, 1) --> (n_samples,). reshape(-1) is more likely to
        # return a view.
        raw_predictions = raw_predictions.reshape(-1)
        # logaddexp(0, x) = log(1 + exp(x))
        loss = np.logaddexp(0, raw_predictions) - y_true * raw_predictions
        return loss.mean() if average else loss

    def get_baseline_prediction(self, y_train, prediction_dim):
        proba_positive_class = np.mean(y_train)
        eps = np.finfo(y_train.dtype).eps
        proba_positive_class = np.clip(proba_positive_class, eps, 1 - eps)
        # log(x / 1 - x) is the anti function of sigmoid, or the link function
        # of the Binomial model.
        return np.log(proba_positive_class / (1 - proba_positive_class))

    def update_gradients_and_hessians(self, gradients, hessians, y_true,
                                      raw_predictions):
        return _update_gradients_hessians_binary_crossentropy(
            gradients, hessians, y_true, raw_predictions)

    def predict_proba(self, raw_predictions):
        # shape (n_samples, 1) --> (n_samples,). reshape(-1) is more likely to
        # return a view.
        raw_predictions = raw_predictions.reshape(-1)
        proba = np.empty((raw_predictions.shape[0], 2), dtype=np.float32)
        proba[:, 1] = expit(raw_predictions)
        proba[:, 0] = 1 - proba[:, 1]
        return proba

@njit(parallel=True, fastmath=True)
def _update_gradients_hessians_binary_crossentropy(gradients, hessians,
                                                   y_true, raw_predictions):
    # Note: using LightGBM version (first mapping {0, 1} into {-1, 1})
    # will cause overflow issues in the exponential as we're using float32
    # precision.

    # shape (n_samples, 1) --> (n_samples,). reshape(-1) is more likely to
    # return a view.
    raw_predictions = raw_predictions.reshape(-1)
    n_samples = raw_predictions.shape[0]
    starts, ends, n_threads = get_threads_chunks(total_size=n_samples)
    for thread_idx in prange(n_threads):
        for i in range(starts[thread_idx], ends[thread_idx]):
            gradients[i] = _expit(raw_predictions[i]) - y_true[i]
            gradient_abs = np.abs(gradients[i])
            hessians[i] = gradient_abs * (1. - gradient_abs)

class CategoricalCrossEntropy(BaseLoss):
    """Categorical cross-entropy loss, for multiclass classification.

    For a given sample x_i, the categorical cross-entropy loss is defined as
    the negative log-likelihood of the model and generalizes the binary
    cross-entropy to more than 2 classes.
    """

    hessian_is_constant = False

    def __call__(self, y_true, raw_predictions, average=True):
        one_hot_true = np.zeros_like(raw_predictions)
        prediction_dim = raw_predictions.shape[1]
        for k in range(prediction_dim):
            one_hot_true[:, k] = (y_true == k)

        loss = (logsumexp(raw_predictions, axis=1) -
                (one_hot_true * raw_predictions).sum(axis=1))
        return loss.mean() if average else loss

    def get_baseline_prediction(self, y_train, prediction_dim):
        init_value = np.zeros(
            shape=(1, prediction_dim),
            dtype=np.float32
        )
        eps = np.finfo(y_train.dtype).eps
        for k in range(prediction_dim):
            proba_kth_class = np.mean(y_train == k)
            proba_kth_class = np.clip(proba_kth_class, eps, 1 - eps)
            init_value[:, k] += np.log(proba_kth_class)

        return init_value

    def update_gradients_and_hessians(self, gradients, hessians, y_true,
                                      raw_predictions):
        return _update_gradients_hessians_categorical_crossentropy(
            gradients, hessians, y_true, raw_predictions)

    def predict_proba(self, raw_predictions):
        # TODO: This could be done in parallel
        # compute softmax (using exp(log(softmax)))
        return np.exp(raw_predictions -
                      logsumexp(raw_predictions, axis=1)[:, np.newaxis])

@njit(parallel=True)
def _update_gradients_hessians_categorical_crossentropy(
        gradients, hessians, y_true, raw_predictions):
    # Here gradients and hessians are of shape
    # (n_samples * prediction_dim,).
    # y_true is of shape (n_samples,).
    # raw_predictions is of shape (n_samples, raw_predictions)
    #
    # Instead of passing the whole gradients and hessians arrays and slicing
    # them here, we could instead do the update in the 'for k in ...' loop of
    # fit(), by passing gradients_at_k and hessians_at_k which are of size
    # (n_samples,).
    # That would however require to pass a copy of raw_predictions, so it does
    # not get partially overwritten at the end of the loop when
    # _update_y_pred() is called (see sklearn PR 12715)
    n_samples, prediction_dim = raw_predictions.shape
    starts, ends, n_threads = get_threads_chunks(total_size=n_samples)
    for k in range(prediction_dim):
        gradients_at_k = gradients[n_samples * k:n_samples * (k + 1)]
        hessians_at_k = hessians[n_samples * k:n_samples * (k + 1)]
        for thread_idx in prange(n_threads):
            for i in range(starts[thread_idx], ends[thread_idx]):
                # p_k is the probability that class(ith sample) == k.
                # This is a regular softmax.
                p_k = np.exp(raw_predictions[i, k] -
                             _logsumexp(raw_predictions[i, :]))
                gradients_at_k[i] = p_k - (y_true[i] == k)
                hessians_at_k[i] = p_k * (1. - p_k)
                # LightGBM uses 2 * p_k * (1 - p_k) which is not stricly
                # correct but equivalent to using half the learning rate.

_LOSSES = {'least_squares': LeastSquares,
           'binary_crossentropy': BinaryCrossEntropy,
           'categorical_crossentropy': CategoricalCrossEntropy}

max_stars_repo_name: P-Tillmann/pygbm

max_issues_repo_name: P-Tillmann/pygbm

max_forks_repo_name: P-Tillmann/pygbm

reason: import numpy,from scipy,from numba