text:
# -*- coding: utf-8 -*-
from __future__ import division
import time
from collections import OrderedDict
import itertools
import ast
import numpy as np
import numpy.ma as ma
#import Splines
#from Splines import spline1d
from scipy.interpolate import interp1d
from scipy.optimize import minimize, curve_fit, leastsq
from scipy.stats import gaussian_kde  # do this for some lists
import pandas as pd
#import matplotlib.pyplot as plt

# TODO: Use Bayes to refine offset estimates given slip rate constraints

def flatten(nested_iterator):
    return list(itertools.chain(*nested_iterator))

def tspline_interpolate():
    pass

def fit_history_spline(age_array, offset_array):

    return interp1d(age_array, offset_array)

def sample_slip_history(age_array, offset_array, time_array,
                        extend_time=False):

    history_spline = fit_history_spline(age_array, offset_array)

    if extend_time == False:
        time_array = time_array[time_array <= np.max(age_array)]
    elif extend_time == True:
        raise Exception('extrapolating not supported yet')

    return history_spline(time_array)

def inverse_transform_sample(vals, probs, n_samps, n_interp=1000, seed=False,
                            seed_val=69):
    pdf_range, pdf_probs = make_pdf(vals, probs, n_interp)
    cdf_range, cdf_probs = make_cdf(pdf_range, pdf_probs)

    cdf_interp = interp1d(cdf_probs, cdf_range, bounds_error=False,
                          fill_value=0.)

    if seed == True:
        np.random.seed(seed_val)

    samps = np.random.rand(n_samps)

    return cdf_interp(samps)

def make_pdf(vals, probs, n_interp=1000):
    val_min = np.min(vals)
    val_max = np.max(vals)

    pdf_range = np.linspace(val_min, val_max, n_interp)

    pmf = interp1d(vals, probs)
    pmf_samples = pmf(pdf_range)
    pdf_probs = pmf_samples / np.sum(pmf_samples)

    return pdf_range, pdf_probs

def make_cdf(pdf_range, pdf_probs):
    return (pdf_range, np.cumsum(pdf_probs))

class OffsetMarker:
    """Represents an offset geologic marker.

    Attributes:
        offsets: list of possible offset distances for the given marker.
            If offset_type = normal, offsets = [mean, sd]
        offset_probs: list of probabilities of corresponding offset distances
        offset_dist_type: offset prob distribution (normal, uniform, arbitrary)
        ages: list of possible ages for the given marker
        age_probs: list of probabilities of corresponding ages
        age_dist_type: age prob. distribution (normal, uniform, arbitrary)
        source: Source for information (e.g., what article, field campaign)

    """
    # TODO: Need to make a random.choice setting for large arrays of vals

    def __init__(self, offsets=np.array([]), offset_probs=None,
                 offset_vals=None, offset_mean=None, offset_median=None,
                 offset_sd=None, offset_mad=None,
                 offset_min=None, offset_max=None,
                 offset_seed=None,
                 offset_dist_type='unspecified', offset_units='unspecified',
                 ages=np.array([]),
                 age_probs=None, age_vals=None,
                 age_mean=None, age_median=None, age_sd=None, age_mad=None,
                 age_min=None, age_max=None,
                 age_seed=None,
                 age_dist_type='unspecified', age_units='unspecified',
                 source='None'):

            self.offsets = offsets
            self.offset_probs = offset_probs
            self.offset_vals = offset_vals
            self.offset_mean = offset_mean
            self.offset_median = offset_median
            self.offset_sd = offset_sd
            self.offset_mad = offset_mad
            self.offset_min = offset_min
            self.offset_max = offset_max
            self.offset_units = offset_units

            if offset_dist_type != 'unspecified':
                self.offset_dist_type = offset_dist_type
            elif offset_dist_type == 'unspecified':
                if offset_mean is not None and offset_sd is not None:
                    self.offset_dist_type = 'normal'
                elif (offset_min is not None and offset_max is not None
                      and offset_sd == None):
                    self.offset_dist_type = 'uniform'
                elif offset_probs is not None and offset_vals is not None:
                    self.offset_dist_type = 'arbitrary'

            self.ages = ages
            self.age_probs = age_probs
            self.age_vals = age_vals
            self.age_mean = age_mean
            self.age_median = age_median
            self.age_sd = age_sd
            self.age_mad = age_mad
            self.age_min = age_min
            self.age_max = age_max
            self.age_units = age_units

            if age_dist_type != 'unspecified':
                self.age_dist_type = age_dist_type
            elif age_dist_type == 'unspecified':
                if age_mean is not None and age_sd is not None:
                    self.age_dist_type = 'normal'
                elif (age_min is not None and age_max is not None
                      and age_sd == None):
                    self.age_dist_type = 'uniform'
                elif age_probs is not None and age_vals is not None:
                    self.age_dist_type = 'arbitrary'

            self.source = source

    def sample_offset_from_normal(self, n):
        """Generates n-length sample from normal distribution of offsets"""

        return sample_from_bounded_normal(self.offset_mean, self.offset_sd, n,
                                          self.offset_min, self.offset_max)

    def sample_offset_from_uniform(self, n):
        """Generates n-length sample from uniform distribution of ages"""

        return np.random.uniform(self.offset_min, self.offset_max, n)

    def sample_offset_from_arbitrary(self, n):
        """not supported yet"""
        offset_sample = inverse_transform_sample(self.offset_vals,
                                                 self.offset_probs, n)
        return offset_sample

    def sample_offset(self, n):
        """Generates n-length array of samples from distribution"""
        if self.offset_dist_type == 'normal':
            offset_sample = self.sample_offset_from_normal(n)

        elif self.offset_dist_type == 'uniform':
            offset_sample = self.sample_offset_from_uniform(n)

        elif self.offset_dist_type == 'arbitrary':
            offset_sample = self.sample_offset_from_arbitrary(n)

        else:
            print('What is the offset distribution type?')

        return offset_sample

    def sample_age_from_normal(self, n):
        """Generates n-length sample from normal distribution of ages"""
        if self.age_min:
            age_min = self.age_min
        else:
            age_min = 0.

        age_sample = sample_from_bounded_normal(self.age_mean, self.age_sd, n,
                                           age_min, self.age_max)

        return age_sample

    def sample_age_from_uniform(self, n):
        """Generates n-length sample from uniform distribution of ages"""
        return np.random.uniform(self.age_min, self.age_max, n)

    def sample_age_from_arbitrary(self, n):
        """not supported yet"""
        return inverse_transform_sample(self.age_vals, self.age_probs, n)

    def sample_age(self, n):
        """Generates n-length array of samples from distribution"""
        if self.age_dist_type == 'normal':
            age_sample = self.sample_age_from_normal(n)

        elif self.age_dist_type == 'uniform':
            age_sample = self.sample_age_from_uniform(n)

        elif self.age_dist_type == 'arbitrary':
            age_sample = self.sample_age_from_arbitrary(n)

        else:
            print('What is the age distribution type?')

        return age_sample

    def sample(self, n):
        age_sample = self.sample_age(n)
        offset_sample = self.sample_offset(n)

        asl = len(age_sample)
        osl = len(offset_sample)

        if asl > osl:
            age_sample = age_sample[0:osl]
        elif osl > asl:
            offset_sample = offset_sample[0:asl]

        return age_sample, offset_sample

def offset_list_from_gui(tabledata, table_header):
    offsets_d = offset_markers_from_gui(tabledata, table_header)

    return list(offsets_d.values())

def offset_markers_from_gui(tabledata, table_header):
    offsets_d = OrderedDict()

    for row in tabledata:
        off_mark_d = offset_marker_dict_from_row(row, table_header)
        offsets_d[off_mark_d['Name']] = offset_marker_from_dict(off_mark_d)

    return offsets_d

def offset_marker_dict_from_row(row, table_header):
    # header_table should be passed from gui
    off_mark_d = OrderedDict()

    for i, key in enumerate(table_header):
        try:
            off_mark_d[key] = ast.literal_eval(row[i])
        except ValueError:
            off_mark_d[key] = row[i]

    for key, val in off_mark_d.items():
        if key in ['Age', 'Age_Err', 'Offset', 'Offset_Err']:
            if not isinstance(val, (list, tuple, np.ndarray)):
                if not isinstance( val, (int, float, complex)):
                    raise Exception(
                        ('Error in {}: value for {} is not numeric. Maybe '
                         +'a string?').format(off_mark_d['Name'], key) )
            else:
                for item in val:
                    if not isinstance( item, (int, float, complex)):
                        raise Exception(
                            ('Error in {}: value for {} is not numeric. Maybe '
                            +'a string?').format(off_mark_d['Name'], key) )
    return off_mark_d

def offset_marker_from_dict(off_row_d):
    or_d = off_row_d

    args = {'offset_units': or_d['Offset_Units'],
            'age_units': or_d['Age_Units']}

    # get offset arguments
    if or_d['Offset_Type'] == 'mean':
        if not np.isscalar(or_d['Offset']):
            raise Exception('Mean Offset has to be a scalar!')
        args['offset_mean'] = or_d['Offset']

    elif or_d['Offset_Type'] == 'median':
        if not np.isscalar(or_d['Offset']):
            raise Exception('Median Offset has to be a scalar!')
        args['offset_median'] = or_d['Offset']

    elif or_d['Offset_Type'] == 'list':
        if len(or_d['Offset']) < 2:
            raise Exception('List Offsets have to be longer than 1!')
        args['offset_vals'] = or_d['Offset']

    else:
        raise Exception('Offset_Type must be mean, median or list!')

    # get offset err arguments
    # TODO: More consistency checking between arg types
    if or_d['Offset_Err_Type'] == 'sd':
        if not np.isscalar(or_d['Offset_Err']):
            raise Exception('sd Offset_Err must be a scalar!')
        args['offset_sd'] = or_d['Offset_Err']

    elif or_d['Offset_Err_Type'] == 'mad':
        if not np.isscalar(or_d['Offset_Err']):
            raise Exception('mad Offset_Err must be a scalar!')
        args['offset_mad'] = or_d['Offset_Err']

    elif or_d['Offset_Err_Type'] == 'minmax':
        if not np.isscalar(or_d['Offset_Err']):
            raise Exception('minmax Offset_Err must be a scalar!')
        if not np.isscalar(or_d['Offset']):
            raise Exception('Mean Offset has to be a scalar!')
        args['offset_min'] = or_d['Offset'] - or_d['Offset_Err']
        args['offset_max'] = or_d['Offset'] + or_d['Offset_Err']
        args['offset_sd'] = None # just to make sure the class inits right

    elif or_d['Offset_Err_Type'] == 'probs':
        if len(or_d['Offset_Err']) < 2:
            raise Exception('probs Offset_Err have to be longer than 1!')
        args['offset_probs'] = or_d['Offset_Err']
        # check to make sure offset vals are set too?

    elif or_d['Offset_Err_Type'] == 'kde':
        if len(or_d['Offset_Err']) < 2:
            raise Exception('kde Offset_Err have to be longer than 1!')
        args['offset_probs'] = kde(or_d['Offset'])

    else:
        raise Exception('Offset_Err_Type must be sd, mad, minmax, probs, '
                        +'or kde!')

    # get age arguments
    if or_d['Age_Type'] == 'mean':
        if not np.isscalar(or_d['Age']):
            raise Exception('Mean Age has to be a scalar!')
        args['age_mean'] = or_d['Age']

    elif or_d['Age_Type'] == 'median':
        if not np.isscalar(or_d['Age']):
            raise Exception('Median Age has to be a scalar!')
        args['age_median'] = or_d['Age']

    elif or_d['Age_Type'] == 'list':
        if len(or_d['Age']) < 2:
            raise Exception('List Ages have to be longer than 1!')
        args['age_vals'] = or_d['Age']

    else:
        raise Exception('Age_Type must be mean, median or list!')

    # get age err arguments
    # TODO: More consistency checking between arg types
    if or_d['Age_Err_Type'] == 'sd':
        if not np.isscalar(or_d['Age_Err']):
            raise Exception('sd Age_Err must be a scalar!')
        args['age_sd'] = or_d['Age_Err']

    elif or_d['Age_Err_Type'] == 'mad':
        if not np.isscalar(or_d['Age_Err']):
            raise Exception('mad Age_Err must be a scalar!')
        args['age_mad'] = or_d['Age_Err']

    elif or_d['Age_Err_Type'] == 'minmax':
        if not np.isscalar(or_d['Age_Err']):
            raise Exception('minmax Age_Err must be a scalar!')
        if not np.isscalar(or_d['Age']):
            raise Exception('Mean Age has to be a scalar!')
        args['age_min'] = or_d['Age'] - or_d['Age_Err']
        args['age_max'] = or_d['Age'] + or_d['Age_Err']
        args['age_sd'] = None # just to make sure the class inits right

    elif or_d['Age_Err_Type'] == 'probs':
        if len(or_d['Age_Err']) < 2:
            raise Exception('probs Age_Err have to be longer than 1!')
        args['age_probs'] = or_d['Age_Err']
        # check to make sure age vals are set too?

    elif or_d['Age_Err_Type'] == 'kde':
        if len(or_d['Age_Err']) < 2:
            raise Exception('kde Age_Err have to be longer than 1!')
        args['age_probs'] = kde(or_d['Age'])

    else:
        raise Exception('Age_Err_Type must be sd, mad, minmax, probs, '
                        +'or kde!')

    return OffsetMarker(**args)

def kde(vals):
    # not sure how to do this yet
    # need to match input length?  or just resample? pass resampling to class?
    # will need to re-do vals too!

    raise Exception('Not Implemented Yet')

def sample_from_bounded_normal(mean, sd, n, sample_min=None, sample_max=None):

    sample = np.random.normal(mean, sd, n)
    sample = trim_distribution(sample, sample_min=sample_min,
                                       sample_max=sample_max)

    while len(sample) < n:
        next_sample = np.random.normal(mean, sd, n)
        next_sample = trim_distribution(next_sample, sample_min, sample_max)
        sample = np.hstack([sample, next_sample])

    return sample[:n]

def trim_distribution(sample, sample_min=None, sample_max=None):

    if sample_min is not None and sample_max is not None:
        if sample_min >= sample_max:
            raise Exception('min must be less than max!')

    if sample_min is not None:
        sample = sample[sample >= sample_min]

    if sample_max is not None:
        sample = sample[sample <= sample_max]

    return sample

def check_monot_increasing(in_array):
    """Checks to see if array is monotonically increasing, returns bool value
    """
    dx = np.diff(in_array)

    return np.all(dx >= 0)

def check_unit_consistency(offset_list):
    off_unit_list = [om.offset_units for om in offset_list]
    age_unit_list = [om.age_units for om in offset_list]

    for off_u in off_unit_list:
        if off_u != off_unit_list[0]:
            raise Exception('OffsetMarker units not consistent.')

    for age_u in age_unit_list:
        if age_u != age_unit_list[0]:
            raise Exception('OffsetMarker units not consistent.')

    return

def get_log_pts(p_min, p_max, n_pts=50, base=np.e):
    """Generates n_pts length array of logarithmically spaced points"""
    if p_min == 0:
        pts_array = np.hstack([0, np.logspace(np.log(1e-5), np.log(p_max),
                                                num=n_pts-1, base=base)])
    else:
        pts_array = np.logspace(p_min, p_max, num=n_pts, base=base)

    return pts_array

def make_age_offset_arrays(offset_list, n, force_increasing=False,
                           zero_offset_age=0., seed=False, seed_value=None,
                           sample_chunks=1):

    # TODO: implement sample chunking (using n samples per marker per fit)

    if seed == True:
        np.random.seed(seed_value)

    age_array = np.zeros((n, len(offset_list)+1 * sample_chunks))
    off_array = np.zeros((n, len(offset_list)+1 * sample_chunks))

    age_array[:,0] = zero_offset_age

    for i, off_mark in enumerate(offset_list):
        col = i+1
        age_array[:,col], off_array[:,col] = off_mark.sample(n)

    if force_increasing == True:

        def make_inc_bool(age_array, off_array, n):

            inc_bool = np.ones((age_array.shape[0]), dtype=int)
            for row in range(n):
                age_inc = check_monot_increasing(age_array[row,:])
                off_inc = check_monot_increasing(off_array[row,:])

                if not (age_inc and off_inc):
                    inc_bool[row] = 0

            inc_bool = np.array(inc_bool, dtype=bool)

            return inc_bool

        inc_bool = make_inc_bool(age_array, off_array, n)

        age_array = age_array[inc_bool, :]
        off_array = off_array[inc_bool, :]

        while age_array.shape[0] < n:

            next_age_array, next_off_array = make_age_offset_arrays(
                                                offset_list, n,
                                                force_increasing=False,
                                                zero_offset_age=zero_offset_age)

            next_inc_bool = make_inc_bool(next_age_array, next_off_array, n)

            next_age_array = next_age_array[next_inc_bool, :]
            next_off_array = next_off_array[next_inc_bool, :]

            off_array = np.vstack([off_array, next_off_array])
            age_array = np.vstack([age_array, next_age_array])

    return age_array[:n,:], off_array[:n,:]

####
# Piecewise linear fitting.  Multiple methods here, pick your poison.
###

def fit_piecewise_linear_w_breakpts(x_data, y_data, breakpts):
    '''
    Modified from an email by Josef Perktold on the
    StatsModels mailing list.
    '''

    # make breakpts into list, so we can prepend 0
    if hasattr(breakpts, 'shape'):
        breakpts = breakpts.tolist()
    else:
        breakpts = list(breakpts)
    breakpts.insert(0,0) # slope over entire array

    # make exog array
    A = np.column_stack([np.maximum(0, x_data - knot) for knot in breakpts])

    # returned slopes are in difference from last slope where slope1 is from 0
    # don't know how to make exog array otherwise
    slopes, sum_sq_err = np.linalg.lstsq(A, y_data)[0:2]

    return np.cumsum(slopes), sum_sq_err # cumsum makes each slope the real one

def piecewise_linear_breakpt_search(x_data, y_data, n_pieces=2, n_iters=20,
                                    penalize_rate_changes=False,
                                    weight_pen=0.2,
                                    allow_slip_reversals=False):
    '''
    docs
    '''

    x_d = x_data - x_data[0] # adjust for zero_offset_age

    n_breaks = n_pieces - 1
    breakpt_samples = np.random.uniform(0., x_d.max(), (n_iters, n_breaks))
    breakpt_samples = np.sort(breakpt_samples, axis=1)
    slopes = {}
    if allow_slip_reversals == False: # 1 means no reversal, 0 means reversal
        monotonic_index = np.zeros(len(breakpt_samples), dtype=int)

    # make this huge so failures won't be selected as min
    sum_sq_errs = np.ones(n_iters) * np.inf
    if penalize_rate_changes == True:
        pen_sum_sq = sum_sq_errs.copy()

    for i, breakpt in enumerate(breakpt_samples):
        try:
            slopes[i], sum_sq_errs[i] = fit_piecewise_linear_w_breakpts(x_d,
                                                                        y_data,
                                                                       breakpt)
            if penalize_rate_changes == True:
                pen_sum_sq[i] = sum_sq_errs[i] * rate_change_penalization(
                                                                    slopes[i],
                                                                    weight_pen)
        except ValueError: # returned when least_sqs problem ill-conditioned
            pass

        if allow_slip_reversals == False:
            monotonic_index[i] = check_slip_monotonicity(slopes[i])
            rev_index = np.bool_(1 - monotonic_index)

    if allow_slip_reversals == False: # give slip reversals inf err, keep inds
        sum_sq_errs[rev_index] = np.inf
        if penalize_rate_changes == True:
            pen_sum_sq[rev_index] = np.inf

    if penalize_rate_changes == True:
        min_i = np.argmin(pen_sum_sq)
    else:
        min_i = np.argmin(sum_sq_errs)

    return flatten([slopes[min_i], breakpt_samples[min_i] + x_data[0],
                    [sum_sq_errs[min_i]]])

def check_slip_monotonicity(rates):
    '''
    Arguments: 'rates', a sequence of slip rates.

    Checks for monotonic slip, i.e. no changes in sign of slip rates.
    If no slip reversals are found, returns True.
    '''
    if all(rate >= 0. for rate in rates) or all(rate <= 0. for rate in rates):
        return 1
    else:
        return 0

def rate_change_penalization(slopes, weight_pen):
    return 1 + np.abs(np.diff(slopes)) * weight_pen

###
# Older, not currently used piecewise fitting stuff
###

def piece_lin_objective(params, x_data, y_data):
    '''docs

    Modified from a function by Andreas Hillboll on the StatsModels
    mailing list.
    '''
    y1 = 0.
    y2, y3, x2 = params
    x1, x3 = x_data[0], x_data[-1]
    Xbefore = y1 + (x_data - x1) * (y2 - y1) / (x2 - x1)
    Xafter = y2 + (x_data - x2) * (y3 - y2) / (x3 - x2)
    Xbreak = np.where(x_data <= x2, Xbefore, Xafter)
    return (ma.masked_invalid(Xbreak - y_data)**2).sum()

def piece_lin_opt(x_data, y_data):

    init_guesses = (np.mean(y_data), np.mean(y_data), np.mean(x_data))
    bounds = ((0, np.max(y_data)), (0., np.max(y_data)), (0., np.max(y_data)))

    res = minimize(piece_lin_objective, init_guesses, (x_data, y_data),
                   method="TNC", bounds=bounds)

    sum_sq_err = piece_lin_objective(res.x, x_data, y_data)

    y2, y3, x2 = res.x

    slope1 = y2 / x2
    slope2 = ((y3 - y2) / (np.max(x_data) - x2))
    breakpoint = x2

    return slope1, slope2, breakpoint, sum_sq_err

def piecewise_linear(x, breakpt, m1, m2):
    return np.piecewise(x, [x < breakpt], [lambda x: m1 * x,
                             lambda x: m2 * x + (m1 * breakpt) - m2 * breakpt])

def piecewise_linear_objective(params, x_data, y_data):
    return ( (y_data - piecewise_linear(x_data, *params))**2).sum()

def penalized_piecewise_linear_objective(params, x_data, y_data, weight=0.1):

    breakpt, m1, m2 = params

    resids = np.array( (y_data - piecewise_linear(x_data, *params)) )

    rate_change_penalization = np.sum(np.abs(resids)) * np.abs(m1 - m2) * weight

    new_resids = np.append(resids, rate_change_penalization)

    return new_resids

def piecewise_linear_opt(x_data, y_data):

    breakpt_guess = np.median(x_data)
    m1_guess = x_data.max() / y_data.max()
    m2_guess = x_data.max() / y_data.max()
    init_vals = [breakpt_guess, m1_guess, m2_guess]

    try:
        params, cov_matrix = curve_fit(piecewise_linear, x_data, y_data,
                                       init_vals)
    except RuntimeError:
        results = minimize(piecewise_linear_objective, init_vals,
                           (x_data, y_data), method='SLSQP')
        #print('slsqp')
        params = results.x

    # params =

    breakpt, m1, m2 = params

    errs = y_data - piecewise_linear(x_data, breakpt, m1, m2)

    sum_sq_err = np.sum(errs**2)

    return m1, m2, breakpt, sum_sq_err

def penalized_piecewise_linear_opt(x_data, y_data, weight=0.3):
    breakpt_guess = np.median(x_data)
    m1_guess = x_data.max() / y_data.max()
    m2_guess = x_data.max() / y_data.max()
    init_vals = (breakpt_guess, m1_guess, m2_guess)

    params, success = leastsq(penalized_piecewise_linear_objective, init_vals,
                              args=(x_data, y_data, weight))

    breakpt, m1, m2 = params

    errs = y_data - piecewise_linear(x_data, breakpt, m1, m2)

    sum_sq_err = np.sum(errs**2)

    return m1, m2, breakpt, sum_sq_err

####
# Other fitting stuff
####

def lin_fit(x_data, y_data):
    x = x_data[:,np.newaxis] - x_data[0]

    # to solve for y = mx + b:
    #x = np.vstack([x_data, np.ones(len(x_data))]).T

    m, _, _, _ = np.linalg.lstsq(x, y_data)
    m = m[0]

    sum_sq_err = ((y_data - (m * x))**2).sum()

    return m, sum_sq_err

def make_linear_results_columns(fit_type=None, n_linear_pieces=None):

    # TODO: fix for arbitrary breakpts

    results_columns = ['m', 'sumsq1']

    if fit_type == 'piecewise':
        m_cols = ['m{}'.format(num + 1) for num in range(n_linear_pieces)]
        results_columns = flatten([m_cols, ['breakpt', 'sumsq2'],
                                   results_columns])
    return results_columns

def do_linear_fits(age_arr, off_arr, fit_type=None, trim_results=True,
                   n_linear_pieces=None, allow_slip_reversals=False):

    n_iters = age_arr.shape[0]
    results_columns = make_linear_results_columns(fit_type, n_linear_pieces)
    results_arr = np.zeros( (n_iters, len(results_columns) ) )

    if fit_type == 'linear':
        for i in range(n_iters):
            xd = age_arr[i,:]
            yd = off_arr[i,:]

            results_arr[i,:] = lin_fit(xd, yd)

    elif fit_type == 'piecewise':
        for i in range(n_iters):
            xd = age_arr[i,:]
            yd = off_arr[i,:]

            results_arr[i, 4:6] = lin_fit(xd, yd)

            #results_arr[i, 0:4] = piece_lin_opt(xd, yd)
            #results_arr[i, 0:4] = piecewise_linear_opt(xd, yd)
            #results_arr[i, 0:4] = penalized_piecewise_linear_opt(xd, yd)
            results_arr[i, 0:4] = piecewise_linear_breakpt_search(xd, yd,
                                    n_pieces=n_linear_pieces,
                                    penalize_rate_changes=True, weight_pen=0.2,
                                    allow_slip_reversals=allow_slip_reversals)

        if allow_slip_reversals==False: #extra filter to catch wiley minnows
            mon_inds = np.bool_([check_slip_monotonicity((results_arr[i,0:2]))
                                for i in range(n_iters)])

    results_df = pd.DataFrame(results_arr, columns=results_columns)

    if fit_type == 'piecewise':
       if trim_results==True:
           # option will be set in the GUI

            results_df = trim_results_df(results_df, age_arr,
                                    allow_slip_reversals=allow_slip_reversals)

    return results_df

def trim_results_df(results_df, age_arr, trim_mag=5,
                    allow_slip_reversals=False):

    results_df = results_df[(results_df.breakpt > age_arr[:,0])
                            &(results_df.breakpt < age_arr[:,-1])]

    m1_75 = results_df.m1.describe()['75%']
    m2_75 = results_df.m2.describe()['75%']
    m1_25 = results_df.m1.describe()['25%']
    m2_25 = results_df.m2.describe()['25%']

    m1_inter_quart_range = m1_75 - m1_25
    m2_inter_quart_range = m2_75 - m2_25

    m1_range = trim_mag * m1_inter_quart_range
    m2_range = trim_mag * m2_inter_quart_range

    results_df = results_df[(np.abs(results_df.m1 - results_df.m1.median())
                             < m1_range)]

    results_df = results_df[(np.abs(results_df.m2 - results_df.m2.median())
                             < m2_range)]

    if allow_slip_reversals == False:
        pos_slip = ((results_df.m1 >= 0.) & (results_df.m2 >= 0.))
        neg_slip = ((results_df.m1 <= 0.) & (results_df.m2 <= 0.))

        results_df = results_df[(pos_slip) ^ (neg_slip)]

    return results_df

def log_likelihood(sum_sq, n):

    return -n / 2 * np.log(sum_sq)

def BIC(log_likelihood, n, p):

    return log_likelihood - ( 0.5 * p * np.log(n / 2 * np.pi))

def AIC(log_likelihood, n, p):
    '''Akaiki's Information Criterion.  Uses same function call as BIC(),
       though *n* is not used.'''

    return 2 * p - 2 * log_likelihood
    return 2 * p - 2 * log_likelihood

def AICc(log_likelihood, n, p):

    aic = AIC(log_likelihood, n, p)

    correction_numerator = 2 * p * (p + 1)
    correction_denominator = (n - p - 1)

    if correction_denominator == 0:
        correction = np.inf
    else:
        correction = correction_numerator / correction_denominator

    return aic + correction

def find_nearest_index(array, value):
    idx = (np.abs(array-value)).argmin()
    return idx

def rate_change_test(results_df, n_offsets, print_res=False):
    results_df['log_like_2'] = log_likelihood(results_df.sumsq2, n_offsets)
    n_iters_out = results_df.shape[0]

    # pn = num params, incl. sum_sq_err and fixed intercept
    p1 = 3 #1 # number of parameters for single linear fit
    p2 = 5 #3 # number of parameters for 2 part piecewise fit

    #if n_offsets > 46:
    #    results_df['bic_1'] = BIC(results_df.log_like_1, n_offsets, p1)
    #    results_df['bic_2'] = BIC(results_df.log_like_2, n_offsets, p2)

    #else:
    #    results_df['bic_1'] = AICc(results_df.log_like_1, n_offsets, p1)
    #    results_df['bic_2'] = AICc(results_df.log_like_2, n_offsets, p2)

    #num_1_count = results_df[results_df.bic_1 > results_df.bic_2].shape[0]
    #num_2_count = n_iters_out - num_1_count
    #num_1_odds = num_1_count / n_iters_out
    #num_2_odds = num_2_count / n_iters_out

    #if num_1_odds > num_2_odds:
    #    n_pieces_best = 1
    #else:
    #    n_pieces_best = 2

    results_df['bic_1'] = AIC(results_df.log_like_1, n_offsets, p1)
    results_df['bic_2'] = AIC(results_df.log_like_2, n_offsets, p2)

    num_1_count = results_df[results_df.bic_1 < results_df.bic_2].shape[0]
    num_2_count = n_iters_out - num_1_count

    num_1_odds = num_1_count / n_iters_out
    num_2_odds = num_2_count / n_iters_out

    if num_1_count > num_2_count:
        n_pieces_best = 1
    else:
        n_pieces_best = 2

    if print_res==True:
        if n_pieces_best == 1:
            print('1 line fits best.  {}/{} ({}% chance)'.format(num_1_count,
                                                               n_iters_out,
                                                               num_1_odds*100))
            print('\nbest fit slip rate results:')
            print(results_df.m.describe())

        else:
            print('2 lines fit best.  {}/{} ({}% chance)'.format(num_2_count,
                                                               n_iters_out,
                                                               num_2_odds*100))
            print('\nbest fit slip rate results:')
            print('rate 1 (younger):')
            print(results_df.m1.describe())
            print('rate change timing:')
            print(results_df.breakpt.describe())
            print('rate 2 (older):')
            print(results_df.m2.describe())
            print('rate_change:')
            print((results_df.m2 - results_df.m1).describe())

    return n_pieces_best

def linear_rate_interp(rate, run_time_max, sim_time_max, zero_offset_age=0.,
                       num_pts=1000):
    ''' Makes a history array of slip rates.  In this case, the slip
    rate is a constant from zero_offset_age to run_time_max, and is
    zero outside of those boundaries.  Returns a Pandas Series.

    Arguments:
    rate (float): slip rate.
    run_time_max (float): Maximum age of slip rate for this MC iteration,
                          i.e. age of oldest offset feature. Times older
                          than this will have zero slip rate.
    sim_time_max (float): Maximum age of oldest feature in the whole MC
                          simulation. This determines the length of the
                          array.
    zero_offset_age (float): Youngest age of faulting.  Times younger than
                             this time will have zero rate.
    num_pts (int): Number of points in the array.
    '''

    times = np.linspace(zero_offset_age, sim_time_max, num_pts)
    slip_rate_history = pd.Series(index=times, data=np.zeros(num_pts))

    slip_rate_history.ix[zero_offset_age : run_time_max] = rate

    return slip_rate_history

def piecewise_rate_interp(rate1, rate2, breakpt, run_time_max, sim_time_max,
                          zero_offset_age=0., num_pts=1000):

    times = np.linspace(zero_offset_age, sim_time_max, num_pts)
    slip_rate_history = np.zeros(num_pts)

    zero_offset_idx = find_nearest_index(times, zero_offset_age)
    run_time_max_idx = find_nearest_index(times, run_time_max)
    breakpt_idx = find_nearest_index(times, breakpt)

    slip_rate_history[zero_offset_idx : breakpt_idx] = rate1
    slip_rate_history[breakpt_idx : run_time_max_idx] = rate2

    return slip_rate_history

def make_rate_hist_array(results_df, age_arr, n_segments=1, num_pts=1000,
                         zero_offset_age=0., return_array=False,
                         sim_time_max='mc_age_max'):

    if sim_time_max == 'mc_age_max':
        sim_time_max = np.max(age_arr)

    times = np.linspace(zero_offset_age, sim_time_max, num_pts)

    rate_hist_df = pd.DataFrame(columns=times, index=results_df.index)
    rate_hist_ar = np.zeros((len(results_df.index), num_pts))

    if n_segments == 1:
        for i in rate_hist_df.index:
            rate = results_df.ix[i, 'm']
            run_time_max = age_arr[i, -1]
            rate_hist_df.ix[i, :] = linear_rate_interp(rate, run_time_max,
                                                       sim_time_max,
                                                       zero_offset_age,
                                                       num_pts)
    elif n_segments == 2:
        for i, row in enumerate(results_df.index):
            rate1 = results_df.ix[row, 'm1']
            rate2 = results_df.ix[row, 'm2']
            breakpt = results_df.ix[row, 'breakpt']
            run_time_max = age_arr[row, -1]
            rate_hist_ar[i, :] = piecewise_rate_interp(rate1, rate2,
                                                       breakpt, run_time_max,
                                                       sim_time_max,
                                                       zero_offset_age,
                                                       num_pts)
    else:
        raise Exception('Only 1 or 2 rates supported now.')

    #return rate_hist_df if return_array == True else rate_hist_df.values
    return rate_hist_ar

def make_cum_hist_array(rate_hist_array):

    return np.cumsum(rate_hist_array, axis=0)

def run_interp_from_gui(offset_list, run_config_dict):
    t0 = time.time()

    rc = run_config_dict

    check_unit_consistency(offset_list)

    n_offsets = len(offset_list) + 1

    print('sampling offset markers')
    age_arr, off_arr = make_age_offset_arrays(offset_list, rc['n_iters'],
                                       force_increasing=rc['force_increasing'],
                                       zero_offset_age=rc['zero_offset_age'],
                                       seed=rc['random_seed'],
                                       seed_value=rc['random_seed_value'])

    print('doing fits')
    if rc['fit_type'] in ['linear', 'piecewise']:
        results_df = do_linear_fits(age_arr, off_arr, fit_type=rc['fit_type'],
                                    n_linear_pieces=rc['n_linear_pieces'],
                                    allow_slip_reversals=rc['slip_reversals'])
    else:
        raise Exception('fit type not implemented yet')

    results_df['log_like_1'] = log_likelihood(results_df.sumsq1, n_offsets)

    if rc['fit_type']  == 'linear':
        print(results_df.m.describe())
        n_pieces_best = 1

    elif rc['fit_type'] == 'piecewise':
        n_pieces_best = rate_change_test(results_df, n_offsets, print_res=True)

    print("\ndone in {:.2f} seconds".format(time.time() - t0))

    return results_df, age_arr, off_arr, n_pieces_best

def trim_age_offset_arrays(res_df, age_arr, off_arr=None):
    """
    Trims age and offset arrays based on retained values from the results_df.
    """
    good_inds = res_df.index.values

    age_arr_trim = age_arr[good_inds, :]

    if off_arr is not None:
        off_arr_trim = off_arr[good_inds, :]

        return age_arr_trim, off_arr_trim

    else:
        return age_arr_trim

def cumulative_offsets(prev_age, prev_rate, new_age, new_rate):
    return prev_age * prev_rate + (new_age - prev_age) * new_rate

reason: import numpy,from scipy


text:
#!/usr/bin/env python2
"""
    Simple interface for drawing the ground truth versions of shadow
    images.

    Shows the original image with the ground truth overlaid.

    - Holding down left click paints onto the ground truth layer.
    - Clicking middle mouse button switches the class you are currently
      painting with (e.g. switches from 'penumbra' to 'shadow' mode).
    - Holding down right-click erases things.
    - Scrolling up/down makes the painting radius larger/smaller.
    - Holding down shift whilst scrolling up/down zooms the image in/out.
    - Holding down shift + left click whilst moving mouse pans the image.
    - Holding down space hides the image (shows only ground truth).
    - Holding down return key hides the ground truth (shows only image).
    - Holding down x/y locks the cursor on that axis (for straight lines).
    - Left/Right keys: Switch to the previous/next image in the set.
    - 's' key saves the ground truth to the output directory.
    - 'f' key flood-fills an area with your currently selected class.
    - Holding down 'e' key shows detected edges in input image, locks
      cursor to nearest edge.

    If previous ground truth images exist, the program will attempt to
    load them when displaying the input images.

    Original Author: Owain Jones [odj@aber.ac.uk]
"""

from __future__ import print_function
import argparse
import pygame
import numpy
import glob
import cv2
import sys
import os

# Pixel classes of the ground truth images.
CLASSES = {
    'shadow': 0,
    'penumbra': 25,
    'object': 153,
    'unknown': 250,
    'background': 255
}
DEFAULT_CLASS = 'background'
zoom = False
image = None
ground_truth = None
pan = (0, 0)
output_fmt = '.png'
output_dir = os.getcwd()
images = []
ground_truths = {}
image_idx = 0
clear = False
ground = True
line = False
alpha = 128
flood = False
axis = None
edges = None
show_edges = False

def class2rgb(pixel_class=DEFAULT_CLASS, alpha=255):
    """Takes the name of a pixel class and returns a pygame.Color RGB
       object."""
    return pygame.Color(*tuple([CLASSES[pixel_class]] * 3 + [alpha]))

class Cursor:
    """Handles storing of data for, and drawing of, a mouse cursor.
       Used for both the program's graphical cursor, and also for
       painting directly onto ground truth images.
    """
    def __init__(self):
        self.x, self.y = (0, 0)
        self.switch_class(DEFAULT_CLASS)
        self.old_class = DEFAULT_CLASS
        self.radius = 16
        self.active = False

    def draw(self, surface, pan=(0, 0)):
        """Paint cursor onto a specified surface. If the cursor is
           active, paints a filled circle. Otherwise, paints a circle
           with an outline.

           Arguments:
           surface: pygame.Surface onto which to draw this cursor.
           pan: If the above surface has been panned, use this tuple to
                specify how much it has been panned by, to offset the
                cursor when drawn to it.

           Class properties affect how it's drawn:
           cursor.active: if True, paint filled circle,
                          if False, paint ring
           cursor.radius: controls the radius of the circle that's drawn
           cursor.current_class_str: Name of the ground truth class that
                                     we're getting the cursor's colour
                                     from (e.g. 'shadow' or 'object').

           Returns: Nothing. Operates directly on the surface passed
                    to it.
        """
        if self.active is False:
            self.lastx, self.lasty = (None, None)
            pygame.draw.circle(
                surface,
                pygame.Color(255 - self.current_class,
                             (255 - self.current_class) * 2 % 255,
                             (255 - self.current_class) * 3 % 255),
                (self.x, self.y), self.radius + 1, 3
            )
            pygame.draw.circle(
                surface,
                class2rgb(self.current_class_str),
                (self.x, self.y), self.radius, 1
            )

        else:
            color = class2rgb(self.current_class_str)
            if (self.lastx, self.lasty) != (None, None):
                pygame.draw.line(
                    surface, color,
                    (self.lastx - pan[0], self.lasty - pan[1]),
                    (self.x - pan[0], self.y - pan[1]),
                    self.radius * 2
                )
            pygame.draw.circle(
                surface,
                class2rgb(self.current_class_str),
                (self.x - pan[0], self.y - pan[1]), self.radius, 0
            )
            self.lastx, self.lasty = (self.x, self.y)

    def switch_class(self, new_class=None):
        """Selects the next class in global list CLASSES, makes this
           cursor use that class."""
        if new_class is None:
            new_class_idx = (self.current_class_idx + 1) % len(CLASSES)
            new_class = list(CLASSES.keys())[new_class_idx]

        self.current_class_str = new_class
        self.current_class_idx = list(CLASSES.keys()).index(new_class)
        self.current_class = CLASSES[new_class]
        pygame.display.set_caption('mode: ' + self.current_class_str)

    def change_radius(self, radius):
        """Change the cursor's radius, as long as it's larger than 2.
           'radius' argument is a delta change to the cursor's actual
           current radius.
        """
        self.radius += radius
        if self.radius < 2:
            self.radius = 2

    def follow_closest_edge(self, edges, pos, pan=(0, 0)):
        """Snap cursor X,Y coordinates to closest white pixel + cursor
           radius."""
        edges_size = edges.get_size()
        black = pygame.Color(0x000000)
        surface_pos = (pos[0] + pan[0], pos[1] + pan[1])
        close = []
        half_radius = int(self.radius / 2)
        for x in range(-half_radius, half_radius + 1):
            for y in range(-half_radius, half_radius + 1):
                corner = (surface_pos[0] + x, surface_pos[1] + y)
                if corner[0] < 0 or corner[1] < 0:
                    continue
                if corner[0] >= edges_size[0] or corner[1] >= edges_size[1]:
                    continue
                if edges.get_at(corner) != black:
                    close.append((half_radius + x, half_radius + y))
        if len(close) > 0:
            self.x = pos[0] + (min(close)[0] - half_radius)
            self.y = pos[1] + (min(close)[1] - half_radius)
        else:
            self.x = pos[0]
            self.y = pos[1]

cursor = Cursor()

def flood_fill(surface, seed_point, color, pan=(0, 0)):
    """Flood fills a pygame surface, starting off at specific seed point.
       Returns the original surface with that area filled in.

       Thanks to wonderfully concise example of (non-recursive)
       flood-fill algorithm in Python:
       http://stackoverflow.com/a/11747550
    """
    seed_point = (seed_point[0] - pan[0], seed_point[1] - pan[1])

    if seed_point[0] > surface.get_width()\
    or seed_point[1] > surface.get_height()\
    or seed_point[0] < 0 or seed_point[1] < 0:
        return surface

    to_fill = set()
    to_fill.add(seed_point)
    background = surface.get_at(seed_point)
    while len(to_fill) > 0:
        x, y = to_fill.pop()
        surface.set_at((x, y), color)
        for i in range(x - 1, x + 2):
            if i < 0 or i > surface.get_width() - 1:
                continue
            for j in range(y - 1, y + 2):
                if j < 0 or j > surface.get_height() - 1:
                    continue
                if color != background\
                and surface.get_at((i, j)) == background:
                    to_fill.add((i, j))

    return surface

def get_ground_truth(input_file, output_dir, output_fmt):
    """Given a path to an input file, the output directory and output
       format (filename suffix), returns what the ground truth file's
       path should be."""
    name = os.path.splitext(os.path.basename(input_file))[0]
    output = os.path.join(output_dir, '%s.%s' % (name, output_fmt))
    return output

def handle_event(event):
    # TODO: This function is far too complex. Bits of it are repeated and could
    # probably be moved into their own functions with some thought.
    """Handle pygame events."""
    if event.type == pygame.QUIT:
        sys.exit(0)

    elif event.type == pygame.KEYDOWN:
        if event.key in [pygame.K_LSHIFT, pygame.K_RSHIFT]:
            global zoom
            zoom = True

        elif event.key == pygame.K_SPACE:
            global clear
            clear = True

        elif event.key in[pygame.K_LCTRL, pygame.K_RCTRL]:
            global line
            line = True

        elif event.key == pygame.K_RETURN:
            global ground
            ground = False

        elif event.key == pygame.K_x:
            global axis
            axis = 'x'
        elif event.key == pygame.K_y:
            global axis
            axis = 'y'

        elif event.key == pygame.K_e:
            global show_edges
            show_edges = True

    elif event.type == pygame.KEYUP:
        if event.key in [pygame.K_LSHIFT, pygame.K_RSHIFT]:
            global zoom
            zoom = False

        elif event.key == pygame.K_s:
            global ground_truth, images, image_idx, output_dir, output_fmt,\
                   ground_truths
            output = get_ground_truth(
                images[image_idx], output_dir, output_fmt
            )
            ground_truths[images[image_idx]] = output
            pygame.image.save(ground_truth, output)
            print('Saved', output)

        elif event.key == pygame.K_LEFT:
            global ground_truth, image, images, image_idx, pan, ground_truths,\
                   edges
            if image_idx > 0:
                image_idx -= 1
                pan = (0, 0)
                if images[image_idx] in ground_truths.keys():
                    ground_truth = pygame.image.load(
                        ground_truths[images[image_idx]]
                    )
                    ground_truth.set_colorkey(class2rgb())
                else:
                    ground_truth.fill(class2rgb())
                image = pygame.image.load(images[image_idx])
                edges = cv_edge_detect(images[image_idx])

        elif event.key == pygame.K_RIGHT:
            global ground_truth, image, images, image_idx, pan, ground_truths
            if image_idx < len(images) - 1:
                image_idx += 1
                pan = (0, 0)
                if images[image_idx] in ground_truths.keys():
                    ground_truth = pygame.image.load(
                        ground_truths[images[image_idx]]
                    )
                    ground_truth.set_colorkey(class2rgb())
                else:
                    ground_truth.fill(class2rgb())
                image = pygame.image.load(images[image_idx])
                edges = cv_edge_detect(images[image_idx])

        elif event.key == pygame.K_SPACE:
            global clear
            clear = False

        elif event.key in [pygame.K_LCTRL, pygame.K_RCTRL]:
            global line
            line = False

        elif event.key == pygame.K_RETURN:
            global ground
            ground = True

        elif event.key == pygame.K_f:
            global ground_truth, cursor, pan
            ground_truth = flood_fill(
                ground_truth, (cursor.x, cursor.y),
                class2rgb(cursor.current_class_str),
                pan
            )

        elif event.key in [pygame.K_x, pygame.K_y]:
            global axis
            axis = None

        elif event.key == pygame.K_e:
            global show_edges
            show_edges = False

    elif event.type == pygame.MOUSEMOTION:
        global cursor, zoom, axis, show_edges, edges
        if zoom is True and cursor.active is True:
            global pan
            pan = (pan[0] + event.rel[0], pan[1] + event.rel[1])
        if show_edges is True:
            global pan
            cursor.follow_closest_edge(edges, event.pos, pan)
        else:
            if axis is not 'y':
                cursor.x = event.pos[0]
            if axis is not 'x':
                cursor.y = event.pos[1]

    elif event.type == pygame.MOUSEBUTTONDOWN:
        global cursor
        if event.button == 1:
            cursor.active = True

        elif event.button == 3:
            cursor.active = True
            cursor.old_class = cursor.current_class_str
            cursor.switch_class(DEFAULT_CLASS)

    elif event.type == pygame.MOUSEBUTTONUP:
        global cursor
        if event.button == 1:
            cursor.active = False

        elif event.button == 2:
            cursor.switch_class()

        elif event.button == 3:
            cursor.active = False
            cursor.switch_class(cursor.old_class)

        elif event.button == 4:
            global zoom
            if zoom is True:
                global image, ground_truth, pan, edges
                image_w, image_h = image.get_size()
                gt_w, gt_h = ground_truth.get_size()
                edges_w, edges_h = edges.get_size()
                image = pygame.transform.scale(
                    image, (image_w * 2, image_h * 2)
                )
                ground_truth = pygame.transform.scale(
                    ground_truth, (gt_w * 2, gt_h * 2)
                )
                edges = pygame.transform.scale(
                    edges, (edges_w * 2, edges_h * 2)
                )
                pan = (pan[0] * 2, pan[1] * 2)

            else:
                cursor.change_radius(1)

        elif event.button == 5:
            global zoom
            if zoom is True:
                global image, ground_truth, pan, edges
                image_w, image_h = image.get_size()
                gt_w, gt_h = ground_truth.get_size()
                edges_w, edges_h = edges.get_size()
                image = pygame.transform.scale(
                    image, (int(image_w / 2), int(image_h / 2))
                )
                ground_truth = pygame.transform.scale(
                    ground_truth, (int(gt_w / 2), int(gt_h / 2))
                )
                edges = pygame.transform.scale(
                    edges, (int(edges_w / 2), int(edges_h / 2))
                )
                pan = (int(pan[0] / 2), int(pan[1] / 2))

            else:
                cursor.change_radius(-1)

    elif event.type == pygame.VIDEORESIZE:
        global window
        window = pygame.display.set_mode(event.size, pygame.RESIZABLE)

def find_input_images(working_dir=os.getcwd(), image_format='jpg'):
    """Find the input images. Defaults to looking for all .jpg files in
       current working directory."""
    images = glob.glob(os.path.join(working_dir, '*.%s' % image_format))
    if images is None:
        return []
    return images

def cv_edge_detect(image):
    """Use OpenCV2 to perform Canny Edge detection on an input image.

       Returns as a new pygame surface."""
    cv_surface = cv2.imread(image, cv2.CV_LOAD_IMAGE_COLOR)
    cv_surface = cv2.cvtColor(cv_surface, cv2.COLOR_BGR2GRAY)
    cv_surface = cv2.fastNlMeansDenoising(cv_surface)
    edges = cv2.Canny(cv_surface, 100, 100, apertureSize=5)
    edges = cv2.cvtColor(edges, cv2.COLOR_GRAY2RGB)
    new_surface = pygame.image.frombuffer(
        edges.tostring(), edges.shape[1::-1], 'RGB'
    )
    new_surface.set_colorkey(0x000000)
    return new_surface

def parse_program_arguments():
    """Parses command line arguments and returns a Namespace."""
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.RawDescriptionHelpFormatter
    )
    parser.add_argument('-o', '--output-dir', default=os.getcwd())
    parser.add_argument('-i', '--input-images', nargs='+', type=str,
                        default=find_input_images())
    parser.add_argument('-f', '--output-format', default='png')
    return parser.parse_args()

def main():
    args = parse_program_arguments()
    if args.input_images is None or len(args.input_images) == 0:
        sys.exit('No input images were found.')

    global window, cursor, image, ground_truth, images, output_fmt, output_dir
    global clear, line, ground, ground_truths, edges, show_edges

    output_dir = args.output_dir
    output_fmt = args.output_format
    images = args.input_images

    # Look for existing ground truth images for input images.
    for image_file in images:
        ground_truth_file = get_ground_truth(
            image_file, output_dir, output_fmt
        )
        if os.path.isfile(ground_truth_file):
            ground_truths[image_file] = ground_truth_file

    # Load the first input image and its ground truth.
    image = pygame.image.load(images[0])
    edges = cv_edge_detect(images[0])
    edges.set_colorkey(0x000000)
    if images[0] in ground_truths.keys():
        ground_truth = pygame.image.load(ground_truths[images[0]])
    else:
        ground_truth = pygame.Surface(image.get_size(), pygame.SRCCOLORKEY)
        ground_truth.fill(class2rgb())
    ground_truth.set_colorkey(class2rgb())

    # Open pygame window, set to the same size as the input image.
    window = pygame.display.set_mode(
        image.get_size(), pygame.RESIZABLE | pygame.SRCALPHA
    )
    cursor = Cursor()
    pygame.mouse.set_visible(False)

    while True:
        window.fill(class2rgb())  # clear previous frame from window
        ground_truth.set_alpha(alpha)  # make ground truth transparent
        handle_event(pygame.event.wait())  # wait for an event, handle it
        for event in pygame.event.get():  # handle any other queued events
            handle_event(event)

        # If left mouse button is down, paint onto the ground truth surface
        if cursor.active is True and zoom is False:
            cursor.draw(ground_truth, pan)

        # Draw the surfaces onto the window.
        # Image is first to be blitted, then ground truth, then the cursor.
        if clear is False:
            window.blit(image, pan)
        if ground is True:
            window.blit(ground_truth, pan)
        if show_edges is True:
            window.blit(edges, pan)
        if cursor.active is False:
            cursor.draw(window)

        pygame.display.flip()

if __name__ == '__main__':
    main()

# vim: set tabstop=4 shiftwidth=4 textwidth=79 cc=72,79:

ext: py

max_stars_repo_name: erinaceous/shadows

max_issues_repo_name: erinaceous/shadows

max_forks_repo_name: erinaceous/shadows

reason: import numpy


text:
import numpy as np
import csv
from math import log

fn="S-cases-age-stp-week-nhser-labeled-1446-PHE.csv"

np.set_printoptions(precision=6,linewidth=120)

d={}
e={}
weeks=set()
nhs=set()
locs=set()
with open(fn,"r") as fp:
  r=csv.reader(fp)
  headings=next(r)
  for row in r:
    week=int(row[0]);weeks.add(week)
    loc=row[1];locs.add(loc)
    nhser=row[2];nhs.add(nhser)
    ageband=row[3]
    agenum=int(row[3].split('-')[0])//10
    neg=float(row[4])
    pos=float(row[5])
    tot=int(row[9])
    d[week,loc,agenum]=(neg,pos)
    if (week,agenum) not in e: e[week,agenum]=[0,0]
    e[week,agenum][0]+=neg
    e[week,agenum][1]+=pos

weeks=sorted(list(weeks))
locs=sorted(list(locs))
nhs=sorted(list(nhs))

for t in range(1,8):
  w0,w1=weeks[0],weeks[-1]
  group0=range(0,t)
  group1=range(t,8)
  acc=k=0
  for loc in locs:
    if 0:
      ok=[(w,loc,a) in d for w in [w0,w1] for a in range(8)]
      if not all(ok): continue
      r0=sum(d[w0,loc,a][0] for a in group0)/sum(d[w0,loc,a][1] for a in group0)
      r1=sum(d[w0,loc,a][0] for a in group1)/sum(d[w0,loc,a][1] for a in group1)
      r2=sum(d[w1,loc,a][0] for a in group0)/sum(d[w1,loc,a][1] for a in group0)
      r3=sum(d[w1,loc,a][0] for a in group1)/sum(d[w1,loc,a][1] for a in group1)
      if r0*r3>0 or r1*r2>0:
        acc+=1
        # r=r0/r1/(r2/r3)
        # print(log(r))
        if r0*r3<=r1*r2: k+=1
    if 1:
      ok=[(w,loc,a) in d for w in [w1] for a in range(8)]
      if not all(ok): continue
      r2=sum(d[w1,loc,a][0] for a in group0)/sum(d[w1,loc,a][1] for a in group0)
      r3=sum(d[w1,loc,a][0] for a in group1)/sum(d[w1,loc,a][1] for a in group1)
      acc+=1
      if r2<=r3: k+=1
    if 0:
      ok=[(w,loc,a) in d for w in weeks for a in range(8)]
      if not all(ok): continue
      r2=sum(d[w,loc,a][0] for a in group0 for w in weeks)/sum(d[w,loc,a][1] for a in group0 for w in weeks)
      r3=sum(d[w,loc,a][0] for a in group1 for w in weeks)/sum(d[w,loc,a][1] for a in group1 for w in weeks)
      acc+=1
      if r2<=r3: k+=1

  rej=len(locs)-acc
  print()
  print("Separating into age groups %d-%d and %d-%d"%(0,t*10,t*10,80))
  print("Locations accepted:",acc)
  print("Locations rejected:",rej)
  print("Number of locations where lower age group has more SGTF:",acc-k,"out of",acc)
  totp=0;p=1;n=acc
  for r in range(k+1):
    totp+=p
    p*=(n-r)/(r+1)
  print("p = %.3g"%(totp/2**n))


text:
import numpy as np
from .translation import translation

def create_points():
    # test data: square of side 1 with lower left corner at (1, 1)
    return np.array([[1.0, 1.0], [1.0, 2.0], [2.0, 2.0], [2.0, 1.0]])

def test_translates_and_returns_translation_factor_with_translation_factor():
    points = create_points()
    translation_factor = np.array([-1.0, 2.0])

    expected_image = np.array([[0.0, 3.0], [0.0, 4.0], [1.0, 4.0], [1.0, 3.0]])

    image, returned_translation_factor = translation(points, translation_factor)

    np.testing.assert_array_equal(image, expected_image)
    np.testing.assert_array_equal(returned_translation_factor, translation_factor)

def test_translates_with_center_of_mass_and_returns_translation_factor_with_no_translation_factor():
    points = create_points()
    expected_image = np.array([[-0.5, -0.5], [-0.5, 0.5], [0.5, 0.5], [0.5, -0.5]])
    expected_translation_factor = -np.mean(points, axis=0)

    image, translation_factor = translation(points)

    np.testing.assert_array_equal(image, expected_image)
    np.testing.assert_array_equal(translation_factor, expected_translation_factor)

max_stars_repo_name: algrs/polliwog

max_issues_repo_name: algrs/polliwog

max_forks_repo_name: algrs/polliwog


text:
# -*- coding: utf-8 -*-
# @Time    : 2019/8/1
# @Author  : Godder
# @Github  : https://github.com/WangGodder
import os
import torch
import sys
from torch.optim.lr_scheduler import StepLR
from torch.utils.data import DataLoader
from torchcmh.utils import get_plotter, AverageMeter
from torchcmh.utils import calc_map_k
from tqdm import tqdm

class TrainBase(object):
    def __init__(self, name, data_name, bit, batch_size, visdom=True, cuda=True):
        self.parameters = {}
        self.loss_store = []
        self.name = name
        self.bit = bit
        self.batch_size = batch_size
        self.cuda = cuda
        self.plotter = get_plotter(self.name) if visdom else None
        self.data_name = data_name
        self.img_model = self.txt_model = None
        self.train_data = self.valid_data = None
        self.save_code_freq = 10
        self.checkpoint_dir = os.path.join('..', 'checkpoints', self.name, data_name)
        if not os.path.exists(self.checkpoint_dir):
            os.makedirs(self.checkpoint_dir)
        self.max_epoch = 500
        self.lr = {}
        self.lr_decay = 0.9
        self.lr_decay_freq = 1
        self.optimizers = []
        self.schedulers = None
        self.max_mapi2t = self.max_mapt2i = self.max_mapi2i = self.max_mapt2t = 0.
        self.best_epoch = 0
        self.qB_img = self.qB_txt = self.rB_img = self.rB_txt = None
        self.best_train_img = self.best_train_txt = None
        self.valid_frequency = 1

    def _init(self):
        self.schedulers = self.__lr_scheduler()
        self.loss_store = self._loss_store_init(self.loss_store)
        print("training %s for %d bit. hyper-paramter list:" % (self.name, self.bit))
        for key, value in self.parameters.items():
            print("%8s: %4.3f" % (key, value))
        print("learning rate decay: %3.2f, decay frequency %2d, learning rate:" % (self.lr_decay, self.lr_decay_freq))
        for key, value in self.lr.items():
            print("%8s: %5.5f" % (key, value))
        print("img net: %s, txt net: %s" % (self.img_model.module_name, self.txt_model.module_name))

    def __lr_scheduler(self):
        """
        create optimizers scheduler for all optimizers.
        default use StepLR scheduler base on lr_decay_freq and lr_decay
        :return:
        """
        schedulers = []
        for optimizer in self.optimizers:
            scheduler = StepLR(optimizer, self.lr_decay_freq, self.lr_decay)
            schedulers.append(scheduler)
        return schedulers

    def lr_schedule(self):
        """
        adjust the learning rate of all optimizers by schedulers.
        :return:
        """
        for scheduler in self.schedulers:
            scheduler.step()

    def plot_loss(self, title, loss_store=None):
        """
        plot loss in loss_store at a figure
        :param title: the title of figure name
        :param loss_store: the loss store to plot, if none, the default loss store will plot
        """
        if loss_store is None:
            loss_store = self.loss_store
        if self.plotter:
            for name, loss in loss_store.items():
                self.plotter.plot(title, name, loss.avg)

    def print_loss(self, epoch, loss_store=None):
        loss_str = "epoch: [%3d/%3d], " % (epoch+1, self.max_epoch)
        if loss_store is None:
            loss_store = self.loss_store
        for name, value in loss_store.items():
            loss_str += name + " {:4.3f}".format(value.avg) + "\t"
        print(loss_str)
        sys.stdout.flush()

    def reset_loss(self, loss_store=None):
        if loss_store is None:
            loss_store = self.loss_store
        for store in loss_store.values():
            store.reset()

    def train(self, num_works=4):
        pass

    def object_function(self, *args):
        pass

    def remark_loss(self, *args):
        """
        store loss into loss store by order
        :param args: loss to store
        :return:
        """
        for i, loss_name in enumerate(self.loss_store.keys()):
            if isinstance(args[i], torch.Tensor):
                self.loss_store[loss_name].update(args[i].item())
            else:
                self.loss_store[loss_name].update(args[i])

    def load_model(self):
        """
        load models from default url.
        :return:
        """
        self.img_model.load_dict(os.path.join(self.checkpoint_dir, str(self.bit) + '-' + self.img_model.module_name + '.pth'))
        self.txt_model.load_dict(os.path.join(self.checkpoint_dir, str(self.bit) + '-' + self.txt_model.module_name + '.pth'))
        if self.cuda:
            self.img_model = self.img_model.cuda()
            self.txt_model = self.txt_model.cuda()

    @staticmethod
    def to_cuda(*args):
        """
        chagne all tensor from cpu tensor to cuda tensor
        :param args: tensors or models
        :return: the tensors or models in cuda by input order
        """
        cuda_args = []
        for arg in args:
            cuda_args.append(arg.cuda())
        return cuda_args

    @staticmethod
    def _loss_store_init(loss_store):
        """
        initialize loss store, transform list to dict by (loss name -> loss register)
        :param loss_store: the list with name of loss
        :return: the dict of loss store
        """
        dict_store = {}
        for loss_name in loss_store:
            dict_store[loss_name] = AverageMeter()
        loss_store = dict_store
        return loss_store

    def valid(self, epoch):
        """
        valid current training model, and save the best model and hash code.
        :param epoch: current epoch
        :return:
        """
        mapi2t, mapt2i, mapi2i, mapt2t, qB_img, qB_txt, rB_img, rB_txt = \
            self.valid_calc(self.img_model, self.txt_model, self.valid_data, self.bit, self.batch_size, return_hash=True, cuda=self.cuda)
        self.max_mapi2i = max(mapi2i, self.max_mapi2i)
        self.max_mapt2t = max(mapt2t, self.max_mapt2t)
        if mapt2i + mapi2t >= self.max_mapi2t + self.max_mapt2i:
            self.max_mapi2t = mapi2t
            self.max_mapt2i = mapt2i
            self.best_epoch = epoch
            self.img_model.save_dict(os.path.join(self.checkpoint_dir, str(self.bit) + '-' + self.img_model.module_name + '.pth'))
            self.txt_model.save_dict(os.path.join(self.checkpoint_dir, str(self.bit) + '-' + self.txt_model.module_name + '.pth'))
            self.qB_img = qB_img.cpu()
            self.qB_txt = qB_txt.cpu()
            self.rB_img = rB_img.cpu()
            self.rB_txt = rB_txt.cpu()
        print('epoch: [%3d/%3d], valid MAP: MAP(i->t): %3.4f, MAP(t->i): %3.4f, MAP(i->i): %3.4f, MAP(t->t): %3.4f \n '
              'max MAP: MAP(i->t): %3.4f, MAP(t->i): %3.4f, MAP(i->i): %3.4f, MAP(t->t): %3.4f. the best epoch is %d' %
              (epoch + 1, self.max_epoch, mapi2t, mapt2i, mapi2i, mapt2t,
               self.max_mapi2t, self.max_mapt2i, self.max_mapi2i, self.max_mapt2t, self.best_epoch + 1))
        if self.plotter:
            self.plotter.plot("mAP", 'i->t', mapi2t.item())
            self.plotter.plot("mAP", "t->i", mapt2i.item())
            self.plotter.plot("mAP", "i->i", mapi2i.item())
            self.plotter.plot("mAP", "t->t", mapt2t.item())
        self.save_code(epoch)

    @staticmethod
    def valid_calc(img_model, txt_model, dataset, bit, batch_size, drop_integer=False, return_hash=False, cuda=True):
        """
        get valid data hash code and calculate mAP
        :param img_model: the image model
        :param txt_model: the txt model
        :param dataset: the valid dataset
        :param bit: the length of hash code
        :param batch_size: the batch size of valid
        :param drop_integer: if true, the excrescent data will be drop
        :param return_hash: if true, the hash codes will be returned
        :param cuda: if use cuda
        :return: mAP and hash codes(if return_hash = True)
        """
        # get query img and txt binary code
        dataset.query()
        qB_img, qB_txt = TrainBase.get_codes(img_model, txt_model, dataset, bit, batch_size, cuda=cuda)
        query_label = dataset.get_all_label()
        # get retrieval img and txt binary code
        dataset.retrieval()
        rB_img, rB_txt = TrainBase.get_codes(img_model, txt_model, dataset, bit, batch_size, cuda=cuda)
        retrieval_label = dataset.get_all_label()
        mAPi2t = calc_map_k(qB_img, rB_txt, query_label, retrieval_label)
        mAPt2i = calc_map_k(qB_txt, rB_img, query_label, retrieval_label)
        mAPi2i = calc_map_k(qB_img, rB_img, query_label, retrieval_label)
        mApt2t = calc_map_k(qB_txt, rB_txt, query_label, retrieval_label)
        if return_hash:
            return mAPi2t, mAPt2i, mAPi2i, mApt2t, qB_img.cpu(), qB_txt.cpu(), rB_img.cpu(), rB_txt.cpu()
        return mAPi2t, mAPt2i, mAPi2i, mApt2t

    def get_train_hash(self):
        img_hash = self.get_img_code(self.img_model, self.train_data, self.bit, self.batch_size)
        txt_hash = self.get_txt_code(self.txt_model, self.train_data, self.bit, self.batch_size)
        return img_hash, txt_hash

    @staticmethod
    def calc_map(qB_img, qB_txt, rB_img, rB_txt, query_label, retrieval_label):
        mAPi2t = calc_map_k(qB_img, rB_txt, query_label, retrieval_label)
        mAPt2i = calc_map_k(qB_txt, rB_img, query_label, retrieval_label)
        return mAPi2t, mAPt2i

    @staticmethod
    def get_codes(img_model, txt_model, dataset, bit, batch_size, cuda=True):
        dataset.both_load()
        dataloader = DataLoader(dataset=dataset, batch_size=batch_size, num_workers=4, pin_memory=True)
        img_buffer = txt_buffer = torch.empty(len(dataset), bit, dtype=torch.float)
        if cuda:
            img_buffer, txt_buffer = TrainBase.to_cuda(img_buffer, txt_buffer)
        img_model.eval()
        txt_model.eval()
        for data in tqdm(dataloader):
            index = data['index'].numpy()
            img = data['img']
            txt = data['txt']
            if cuda:
                img, txt = TrainBase.to_cuda(img, txt)
            img_hash = torch.tanh(img_model(img))
            txt_hash = torch.tanh(txt_model(txt))
            img_buffer[index, :] = img_hash.data
            txt_buffer[index, :] = txt_hash.data
        return img_buffer, txt_buffer

    @staticmethod
    def get_img_code(img_model, dataset, bit, batch_size, drop_integer=False, cuda=True):
        dataset.img_load()
        dataloader = DataLoader(dataset=dataset, batch_size=batch_size, num_workers=4, pin_memory=True, drop_last=drop_integer)
        B_img = torch.zeros(len(dataset), bit, dtype=torch.float)
        if cuda:
            B_img = B_img.cuda()
        img_model.eval()
        for data in tqdm(dataloader):
            index = data['index'].numpy()  # type: np.ndarray
            img = data['img']  # type: torch.Tensor
            if cuda:
                img = img.cuda()

            f = img_model(img)
            B_img[index, :] = f.data
        return B_img

    @staticmethod
    def get_txt_code(txt_model, dataset, bit, batch_size, drop_integer=False, cuda=True):
        dataset.txt_load()
        dataloader = DataLoader(dataset=dataset, batch_size=batch_size, num_workers=4, pin_memory=True, drop_last=drop_integer)
        B_txt = torch.zeros(len(dataset), bit, dtype=torch.float)
        if cuda:
            B_txt = B_txt.cuda()
        txt_model.eval()
        for data in tqdm(dataloader):
            index = data['index'].numpy()  # type: np.ndarray
            txt = data['txt']  # type: torch.Tensor
            if cuda:
                txt = txt.cuda()

            g = txt_model(txt)
            B_txt[index, :] = g.data
        return B_txt

    def save_code(self, epoch):
        """
        save valid hash code in default url.
        :param epoch: current epoch
        :return:
        """
        if (epoch+1) % self.save_code_freq != 0 or epoch - self.best_epoch < 10 is None:
            return
        qB_img = torch.sign(self.qB_img).numpy().astype(int)
        qB_txt = torch.sign(self.qB_txt).numpy().astype(int)
        rB_img = torch.sign(self.rB_img).numpy().astype(int)
        rB_txt = torch.sign(self.rB_txt).numpy().astype(int)
        import scipy.io as sio
        path = os.path.join(self.checkpoint_dir, str(self.bit) + '-hash.mat')
        sio.savemat(path, {"q_img": qB_img, "q_txt": qB_txt, "r_img": rB_img, "r_txt": rB_txt})
        print("save hash code in %s" % path)

    def resume_train(self, state_path):
        pass

max_stars_repo_name: chrisbyd/deep-cross-modal-hashing

max_issues_repo_name: silencelzx/deep-cross-modal-hashing

max_forks_repo_name: silencelzx/deep-cross-modal-hashing

reason: import scipy


text:
# Copyright (c) Facebook, Inc. and its affiliates.
#
# 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.
"""Gap"""
"""Template task in which you prevent something from falling so ball can roll into container."""
import numpy as np
import phyre.creator as creator_lib
import phyre.virtual_tools as vt

@creator_lib.define_task_template(
    seed=range(1000),
    version="2",
    search_params=dict(required_flags=['BALL:GOOD_STABLE'],
                       excluded_flags=['BALL:TRIVIAL'],
                       diversify_tier='ball'),
)
def build_task(C, seed):
    rng = np.random.RandomState(seed=seed)
    slopeLeft = [300, 450]
    slopeRight = [200, 350]
    slopeWidth = [200, 300]
    ballXPos = [50, 150]
    ballYPos = [15, 50]
    gapWidth = [60, 200]
    supportX = [5, 25]
    supportY = [-5, 80]
    lidThick = [15, 40]

    flip_lr = rng.uniform(0, 1) < 0.5
    ## Make the slope
    sL = rng.uniform(slopeLeft[0], slopeLeft[1])
    sR = rng.uniform(slopeRight[0], slopeRight[1])
    sW = rng.uniform(slopeWidth[0], slopeWidth[1])
    slopeVerts = [[0.0, 0.0], [0.0, sL], [sW, sR], [sW, 0]]

    ## Make the gap
    gW = rng.uniform(gapWidth[0], gapWidth[1])
    strutSupportBox = [sW, 0, sW + gW, sR - 150]

    ## Make the support
    jitterX = rng.uniform(supportX[0], supportX[1])
    jitterY = rng.uniform(supportY[0], supportY[1])
    lidT = rng.uniform(lidThick[0], lidThick[1])
    supportBox = [
        sW + jitterX, sR + jitterY, sW + gW - jitterX, sR + jitterY + lidT
    ]

    ## Make the goal
    goalSegs = [[sW + gW, sR - 50], [sW + gW, 5], [550, 5], [550, sR - 50]]

    ## Find the ball position
    bX = rng.uniform(ballXPos[0], ballXPos[1])
    bY = sL * (sW - bX) / sW + sR
    bpos = [
        bX, bY +
        rng.uniform(ballYPos[0], np.min([ballYPos[1], vt.VT_SCALE - bY - 15]))
    ]

    ## Make the world getting into the container
    slopeVerts.reverse()

    if flip_lr:
        slopeVerts = vt.flip_left_right(slopeVerts)
        bpos = vt.flip_left_right(bpos)

    slope = C.add_convex_polygon(vt.convert_phyre_tools_vertices(slopeVerts),
                                 False)
    container, _ = vt.add_container(C,
                                    goalSegs,
                                    10,
                                    False,
                                    True,
                                    flip_lr=flip_lr)
    ball = C.add('dynamic ball',
                 30. / vt.VT_SCALE,
                 center_x=bpos[0] * C.scene.width / vt.VT_SCALE,
                 center_y=bpos[1] * C.scene.height / vt.VT_SCALE)
    strutBox = vt.add_box(C, strutSupportBox, False, flip_lr=flip_lr)
    support = vt.add_box(C, supportBox, True, flip_lr=flip_lr)

    C.update_task(body1=ball,
                  body2=container,
                  relationships=[C.SpatialRelationship.TOUCHING])

    C.set_meta(C.SolutionTier.VIRTUAL_TOOLS)


text:
import ignite.engine as igne
import ignite.metrics as ignm
import numpy as np
import torch
import torch.nn as nn
import torch.utils.data

MAX_HERO_ID = 121

class Hero2Vec(nn.Module):
    # we input 4 heroes and the network should output the fifth one
    INPUT_HEROES = 4

    def __init__(self, emb_dim: int = 128, hero_count: int = MAX_HERO_ID + 1):
        super().__init__()

        self.emb_dim = emb_dim
        self.hero_count = hero_count

        self.hero_embedding = nn.Linear(hero_count, emb_dim, bias=False)
        self.missing_hero_logits = nn.Linear(emb_dim, hero_count)

    def forward(self, inputs):
        embedded_heroes = self.hero_embedding(inputs)
        team_embedding = torch.mean(embedded_heroes, dim=1)
        missing_hero = self.missing_hero_logits(team_embedding)
        return missing_hero

class HeroDataset:
    def __init__(self):
        self.inputs = np.load('data/processed/inputs.npy').astype('float32')
        self.labels = np.load('data/processed/labels.npy').astype('long')

        assert len(self.inputs) == len(self.labels)

    def __getitem__(self, item):
        return self.inputs[item], self.labels[item]

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

def get_train_loader(batch_size=256):
    dataset = HeroDataset()

    train_loader = torch.utils.data.DataLoader(
        dataset,
        batch_size=batch_size,
        shuffle=True,
        num_workers=1,
    )
    return train_loader

def train():
    model = Hero2Vec()
    optimizer = torch.optim.SGD(model.parameters(), lr=0.1, momentum=0.99)
    loss = nn.CrossEntropyLoss()
    trainer = igne.create_supervised_trainer(model, optimizer, loss)
    evaluator = igne.create_supervised_evaluator(
        model,
        metrics={
            'accuracy': ignm.CategoricalAccuracy(),
            'loss': ignm.Loss(loss),
        },
    )

    train_loader = get_train_loader(batch_size=1024)

    @trainer.on(igne.Events.ITERATION_COMPLETED)
    def log_training_loss(trainer):
        if trainer.state.iteration % 100 == 0:
            print(f'[Epoch {trainer.state.epoch}] Loss: {trainer.state.output:.2f}')

    @trainer.on(igne.Events.EPOCH_COMPLETED)
    def log_training_evaluation(trainer):
        evaluator.run(train_loader)
        metrics = evaluator.state.metrics
        print(f'[Epoch {trainer.state.epoch}] Accuracy: {metrics["accuracy"]}, Loss: {metrics["loss"]}')

    evaluator.run(train_loader)
    metrics = evaluator.state.metrics
    print(f'[Before training] Accuracy: {metrics["accuracy"]}, Loss: {metrics["loss"]}')

    trainer.run(train_loader, max_epochs=10)

    import pudb; pu.db
    pass

if __name__ == '__main__':
    train()


text:
# -*- coding: utf-8 -*-
# @Author: patrick
# @Date:   2016-09-01 17:04:53
# @Last Modified by:   Patrick Bos
# @Last Modified time: 2016-10-06 07:30:16

# as per tensorflow styleguide
# https://www.tensorflow.org/versions/r0.11/how_tos/style_guide.html
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
from timeit import default_timer as timer
import time

def apply_constraint(var, constraints):
    var_name = var.name[:var.name.find(':')]
    # low = tf.constant(constraints[var_name][0], dtype=tf.float64)
    # high = tf.constant(constraints[var_name][1], dtype=tf.float64)
    low = constraints[var_name][0]
    high = constraints[var_name][1]
    return tf.assign(var, tf.clip_by_value(var, low, high),
                     name="assign_to_" + var_name)
    # return tf.Variable(tf.clip_by_value(var, low, high), name=var_name + '_clipped')

project_dn = "/home/patrick/projects/apcocsm/"
# project_dn = "/home/pbos/apcocsm/"

m0_num = 5.291
argpar_num = -20.0

constraint = {}

constraint['sigmean'] = (5.20, 5.30)
constraint['sigwidth'] = (0.001, 1.)
constraint['argpar'] = (-100., -1.)
constraint['nsig'] = (0., 10000)
constraint['nbkg'] = (0., 10000)
constraint['mes'] = (5.20, 5.30)

pi = tf.constant(np.pi, dtype=tf.float64, name="pi")
sqrt2pi = tf.constant(np.sqrt(2 * np.pi), dtype=tf.float64, name="sqrt2pi")
two = tf.constant(2, dtype=tf.float64, name="two")
one = tf.constant(1, dtype=tf.float64, name="one")
zero = tf.constant(0, dtype=tf.float64, name="zero")

def gaussian_pdf(x, mean, std):
    val = tf.div(tf.exp(-tf.pow((x - mean) / std, 2) / two), (sqrt2pi * std),
                 name="gaussian_pdf")
    return val

def argus_pdf(m, m0, c, p=0.5):
    t = m / m0
    u = 1 - t * t
    return tf.select(tf.greater_equal(t, one),
                     tf.zeros_like(m),
                     m * tf.pow(u, p) * tf.exp(c * u),
                     name="argus_pdf")
    # return tf.cond(tf.greater_equal(t, one),
    #                lambda: zero,
    #                lambda: m * tf.pow(u, p) * tf.exp(c * u),
    #                name="argus_pdf")
    # N.B.: bij cond moeten de argumenten functies zijn (zonder argumenten)
    #       zodat tf ze pas hoeft te callen / uit te rekenen als ze nodig zijn.
    #       Dat is dus bij select niet mogelijk, daar krijg je meteen beide hele
    #       tensors.

def argus_integral_phalf(m_low, m_high, m0, c):
    """
    Only valid for argus_pdf with p=0.5! Otherwise need to do numerical
    integral.
    """
    def F(x):
        return -0.5 * m0 * m0 * (tf.exp(c * x) * tf.sqrt(x) / c + 0.5 / tf.pow(-c, 1.5) * tf.sqrt(pi) * tf.erf(tf.sqrt(-c * x)))

    a = tf.minimum(m_low, m0)
    b = tf.minimum(m_high, m0)

    x1 = 1 - tf.pow(a / m0, 2)
    x2 = 1 - tf.pow(b / m0, 2)

    area = tf.sub(F(x2), F(x1), name="argus_integral_phalf")
    return area

def argus_integral_phalf_numpy(m_low, m_high, m0, c):
    """
    Only valid for argus_pdf with p=0.5! Otherwise need to do numerical
    integral.
    """
    import scipy.special

    def F(x):
        return -0.5 * m0 * m0 * (np.exp(c * x) * np.sqrt(x) / c + 0.5 / (-c)**1.5 * np.sqrt(np.pi) * scipy.special.erf(np.sqrt(-c * x)))

    a = np.min([m_low, m0])
    b = np.min([m_high, m0])

    x1 = 1 - (a / m0)**2
    x2 = 1 - (b / m0)**2

    area = F(x2) - F(x1)
    return area

argus_numerical_norm = tf.constant(argus_integral_phalf_numpy(constraint['mes'][0],
                                                              constraint['mes'][1],
                                                              m0_num, argpar_num),
                                   dtype=tf.float64, name="argus_numerical_norm")

def argus_pdf_phalf_WN(m, m0, c, m_low, m_high, tf_norm=tf.constant(False)):
    """
    WN: with normalization
    tf_norm: use the tensorflow integral function (True) or the numpy one (False)
    """
    norm = tf.cond(tf_norm,
                   lambda: argus_integral_phalf(m_low, m_high, m0, c),
                   lambda: argus_numerical_norm, name="argus_norm")
    # norm = argus_numerical_norm
    # norm = argus_integral_phalf(m_low, m_high, m0, c)
    return argus_pdf(m, m0, c) / norm

# // --- Observable ---
# RooRealVar mes("mes","m_{ES} (GeV)",5.20,5.30) ;

# // --- Build Gaussian signal PDF ---
# RooRealVar sigmean("sigmean","B^{#pm} mass",5.28,5.20,5.30) ;
# RooRealVar sigwidth("sigwidth","B^{#pm} width",0.0027,0.001,1.) ;

sigmean = tf.Variable(np.float64(5.28), name="sigmean")
sigwidth = tf.Variable(np.float64(0.0027), name="sigwidth")

# RooGaussian gauss("gauss","gaussian PDF",mes,sigmean,sigwidth) ;

# // --- Build Argus background PDF ---
# RooRealVar argpar("argpar","argus shape parameter",-20.0,-100.,-1.) ;
# RooConstVar m0("m0", "resonant mass", 5.291);

argpar = tf.Variable(np.float64(argpar_num), name="argpar")
m0 = tf.constant(np.float64(m0_num), name="m0")

# RooArgusBG argus("argus","Argus PDF",mes,m0,argpar) ;

# // --- Construct signal+background PDF ---
# RooRealVar nsig("nsig","#signal events",200,0.,10000) ;
# RooRealVar nbkg("nbkg","#background events",800,0.,10000) ;

nsig = tf.Variable(np.float64(200), name="nsig")
nbkg = tf.Variable(np.float64(800), name="nbkg")

# RooAddPdf sum("sum","g+a",RooArgList(gauss,argus),RooArgList(nsig,nbkg)) ;

# // --- Generate a toyMC sample from composite PDF ---
# RooDataSet *data = sum.generate(mes,2000) ;

def sum_pdf(mes, nsig, sigmean, sigwidth, nbkg, m0, argpar, mes_low, mes_high):
    add = tf.add(nsig * gaussian_pdf(mes, sigmean, sigwidth),
                 nbkg * argus_pdf_phalf_WN(mes, m0, argpar, mes_low, mes_high),
                 name="sum_pdf")
    return tf.div(add, nsig + nbkg, name="sum_pdf_normalized")

# data in RooFit genereren en importeren
# draai dit in ROOT:
# data.write("roofit_demo_random_data_values.dat");
data_raw = np.loadtxt(project_dn + "roofit_demo_random_data_values.dat",
                      dtype=np.float64)
data = tf.constant(data_raw, name='event_data')

# // --- Perform extended ML fit of composite PDF to toy data ---
# sum.fitTo(*data,"Extended") ;

# convert to tf constants, otherwise you'll get complaints about float32s...
for key in constraint.keys():
    low = constraint[key][0]
    high = constraint[key][1]
    constraint[key] = (tf.constant(low, dtype=tf.float64),
                       tf.constant(high, dtype=tf.float64))

# nll = tf.neg(tf.reduce_sum(tf.log(tf.map_fn(lambda mes: sum_pdf(mes, nsig, sigmean, sigwidth, nbkg, m0, argpar, constraint['mes'][0], constraint['mes'][1]), data))), name="nll")

print("N.B.: using direct data entry")
nll = tf.neg(tf.reduce_sum(tf.log(sum_pdf(data, nsig, sigmean, sigwidth, nbkg, m0, argpar, constraint['mes'][0], constraint['mes'][1]))), name="nll")

# print("N.B.: using unsummed version of nll! This appears to be the way people minimize cost functions in tf...")
# nll = tf.neg(tf.log(sum_pdf(data, nsig, sigmean, sigwidth, nbkg, m0, argpar, constraint['mes'][0], constraint['mes'][1])), name="nll")

# grad = tf.gradients(nll, [mu, sigma])

max_steps = 1000
status_every = 100

# sigmean_c = apply_constraint(sigmean, constraint)
# sigwidth_c = apply_constraint(sigwidth, constraint)
# argpar_c = apply_constraint(argpar, constraint)
# nsig_c = apply_constraint(nsig, constraint)
# nbkg_c = apply_constraint(nbkg, constraint)

# update_vars = [sigmean_c, sigwidth_c, argpar_c, nsig_c, nbkg_c]

variables = tf.all_variables()

# Create an optimizer with the desired parameters.
# opt = tf.train.GradientDescentOptimizer(learning_rate=0.001)
# opt = tf.train.AdagradOptimizer(learning_rate=0.1)
opt = tf.train.AdamOptimizer()
opt_op = opt.minimize(nll)

tf.scalar_summary('nll', nll)

init_op = tf.initialize_all_variables()
# check_op = tf.report_uninitialized_variables()

# start session
with tf.Session() as sess:
    # Merge all the summaries and write them out to /tmp/mnist_logs (by default)
    summarize_merged = tf.merge_all_summaries()
    summary_writer = tf.train.SummaryWriter('./train_%i' % int(time.time()), sess.graph)
    # Run the init operation.
    sess.run(init_op)
    # print sess.run(init_op)
    # print sess.run(check_op)

    true_vars = {}
    for v in variables:
        key = v.name[:v.name.find(':')]
        true_vars[key] = v.eval()

    true_vars['m0'] = m0.eval()

    print("name\t" + "\t".join([v.name.ljust(10) for v in variables]) + "\t | nll\t\t\t | step")
    print("init\t" + "\t".join(["%6.4e" % v for v in sess.run(variables)]) + "\t | %f" % sess.run(nll))
    print

    start = timer()

    for step in xrange(max_steps):
        # print "variables 3:", sess.run(variables)
        summary, _ = sess.run([summarize_merged, opt_op])
        # sess.run([opt_op])
        summary_writer.add_summary(summary, step)

        if step % status_every == 0:
            var_values_opt = sess.run(variables)
            nll_value_opt = sess.run(nll)
            # sess.run(update_vars)
            # var_values_clip = np.array(sess.run(variables))
            # nll_value_clip = np.array(sess.run(nll))
            print("opt\t" + "\t".join(["%6.4e" % v for v in var_values_opt]) + "\t | %f\t | %i" % (nll_value_opt, step))
        # clipped = np.where(var_values_opt == var_values_clip, [" "*10] * len(variables), ["%6.4e" % v for v in var_values_clip])
        # print "clip\t" + "\t".join(clipped) + "\t | %f" % nll_value_clip

    end = timer()

    print("Loop took %f seconds" % (end - start))
    raise Exception

    fit_vars = {}
    for v in variables:
        key = v.name[:v.name.find(':')]
        fit_vars[key] = v.eval()

    fit_vars['m0'] = m0.eval()

    counts, bins = np.histogram(data.eval(), bins=100)
    x_bins = (bins[:-1] + bins[1:]) / 2

    y_fit = [sum_pdf(x, mes_low=constraint['mes'][0], mes_high=constraint['mes'][1], **fit_vars).eval() for x in x_bins]
    argus_fit = [fit_vars['nbkg'] * argus_pdf_phalf_WN(x, fit_vars['m0'], fit_vars['argpar'], m_low=constraint['mes'][0], m_high=constraint['mes'][1]).eval() for x in x_bins]

    y_true = [sum_pdf(x, mes_low=constraint['mes'][0], mes_high=constraint['mes'][1], **true_vars).eval() for x in x_bins]

    # normalize fit values to data counts
    y_fit_norm = np.sum(counts) / np.sum(y_fit)
    y_fit = [y * y_fit_norm for y in y_fit]
    # argus_fit_norm = np.sum(counts) / np.sum(argus_fit)
    argus_fit = [a * y_fit_norm for a in argus_fit]
    y_true_norm = np.sum(counts) / np.sum(y_true)
    y_true = [y * y_true_norm for y in y_true]

    plt.errorbar(x_bins, counts, yerr=np.sqrt(counts), fmt='.g')
    plt.plot(x_bins, y_fit, '-b')
    plt.plot(x_bins, argus_fit, '--b')
    plt.plot(x_bins, y_true, ':k')
    plt.show()

# tf.InteractiveSession()

# sess = tf.Session()

# sess.run(init_op)

# opt = tf.train.GradientDescentOptimizer(learning_rate=1)
# opt_op = opt.minimize(nll, var_list=[sigmean, sigwidth, argpar, nsig, nbkg])
# for step in xrange(10):
#     out = sess.run([opt_op, nll, sigmean, sigwidth, argpar, nsig, nbkg])
#     print out[1:]

# sess.close()

# // --- Plot toy data and composite PDF overlaid ---
# RooPlot* mesframe = mes.frame() ;
# data->plotOn(mesframe) ;
# sum.plotOn(mesframe) ;
# sum.plotOn(mesframe,Components(argus),LineStyle(kDashed)) ;
# mesframe->Draw();

max_stars_repo_name: NLeSC/parallel-roofit-scripts

max_issues_repo_name: NLeSC/parallel-roofit-scripts

max_forks_repo_name: roofit-dev/parallel-roofit-scripts

reason: import numpy,import scipy


text:
from fastapi import FastAPI
from starlette.responses import JSONResponse
import numpy as np
#import torch
#import torch.nn as nn

app = FastAPI()

#redefine model class
#class PytorchMultiClass(nn.Module):
    #def __init__(self, num_features):
        #super(PytorchMultiClass, self).__init__()

        #self.layer_1 = nn.Linear(num_features, 32)
        #self.layer_out = nn.Linear(32, 104)  # from 4, 104 number of classes
        #self.softmax = nn.Softmax(dim=1)

    #def forward(self, x):
        #x = F.dropout(F.relu(self.layer_1(x)), training=self.training)
        #x = self.layer_out(x)
        #return self.softmax(x)

#num_features = 4
#model = PytorchMultiClass(num_features) #*args, **kwargs
#model.load_state_dict(torch.load('../models/NN_4feat.pt'))
#model.eval()

# Inside the main.py file, create a function called read_root() that will describe the project
@app.get("/")
def read_root():
    return "This project was designed to take a few key features from a beer " \
           "styles data set and then build a model that would take those inputs " \
           "and then predict the given beer style. It takes four inputs (review_aroma, " \
           "review_appearance, review_palate, review_taste) and each one must " \
           "be given an integer score from 0-5 and then it will predict the beer " \
           "style. Here is the github repo: https://github.com/freescania/advdsi_at2"

# Inside the main.py file, create a function called healthcheck() that will
# will state 'Neural Network model is all ready to go!'. Add a decorator to it in order
# to add a GET endpoint to app on /health with status code 200
@app.get('/health', status_code=200)
def healthcheck():
    return 'Neural Network model is all ready to go!'

#Inside the main.py file, create a function called format_features() with
# review_aroma, review_appearance, review_palate and review_taste as input parameters that will return a
# dictionary with the names of the features as keys and the inpot parameters as lists
def format_features(review_aroma: int,	review_appearance: int, review_palate: int, review_taste: int):
  return {
        'review aroma (0-5)': ([review_aroma]/5),
        'review appearance (0-5)': ([review_appearance]/5),
        'review palate (0-5)': ([review_palate]/5),
        'review taste (0-5)': ([review_taste]/5)
    }

#Inside the main.py file, Define a function called predict with the following logics:
#input parameters: review_aroma, review_appearance, review_palate, review_taste
#logics: format the input parameters as dict, convert it to a dataframe and make prediction with gmm_pipe
#output: prediction as json
#Add a decorator to it in order to add a GET endpoint to app on /beer/type
@app.get("/beer/type/")
def predict(review_aroma: int,	review_appearance: int, review_palate: int, review_taste: int):
    features = format_features(review_aroma, review_appearance, review_palate, review_taste)
    feats = np.array(list(features.values())).T
    feats = (feats / 5)
    zeroes = np.zeros((32, 4))
    zeroes[0] = feats
    #obs = torch.Tensor(zeroes)
    #change for NN
    output = ['American IPA']
    return JSONResponse(output.tolist())

@app.get("/model/architecture/")
def architecture():
    return 'The model architecture is four layers: layer_1,  Linear, input=4, output=32,' \
           'layer_2,  Linear, input=32, output=104, layer_3,  Linear, input=104, output=32,' \
           'layer_out, Linear, input=32, output=104,  Softmax activation'


text:
import argparse
import pandas as pd
import json
import pickle
import numpy as np
from pathlib import Path
from scipy.special import softmax
import shutil

parser = argparse.ArgumentParser()
parser.add_argument('--config', type=str)
parser.add_argument('--raw_result_dir', type=str)
parser.add_argument('--output_dir', type=str)
args = parser.parse_args()

with open(args.config, 'r') as f:
    config_json = json.load(f)

for config in config_json:
    class_bias_result = []
    for no_try in range(16):
        exp_result_path = Path(
            args.raw_result_dir,
            "{0}_{1}_{2}_{3}/{4}".format(config['network'],
                                        config['training_type'],
                                        config['dataset'],
                                        config['random_seed'],
                                        str(no_try)))

        if not ('adv-conv' in config['training_type']):
            if config['training_type'] == 'no_gender':
                log_dir = 'origin_0'
            else:
                log_dir = config['training_type'].replace('-','_')

            ckpt_path = Path(exp_result_path, 'models', log_dir, 'checkpoint.pth.tar')
        else:
            if config['dataset'] == 'coco':
                if config['training_type'] == 'adv-conv4':
                    log_dir = 'conv4_300_2.0_0.2_0_linear'
                elif config['training_type'] == 'adv-conv5':
                    log_dir = 'conv5_300_2.0_0.2_0_linear'
            elif config['dataset'] == 'imSitu':
                if config['training_type'] == 'adv-conv4':
                    log_dir = 'conv4_300_1.0_0.2_linear_0'
                elif config['training_type'] == 'adv-conv5':
                    log_dir = 'conv5_300_1.0_0.2_linear_0'

            ckpt_path = Path(exp_result_path, 'adv/models', log_dir, 'checkpoint.pth.tar')

        if config['dataset'] == 'coco':
            tech = 'C-'
        elif config['dataset'] == 'imSitu':
            tech = 'I-'

        if config['training_type'] == 'no_gender':
            tech += 'Base'
        elif config['training_type'].startswith('ratio-'):
            tech += 'R' + config['training_type'][-1]
        elif config['training_type'].startswith('adv-conv'):
            tech += 'A' + config['training_type'][-1]

        copy_path = Path(args.output_dir, tech, 'run_' + str(no_try).zfill(2) + '.pth')
        copy_path.parent.mkdir(parents=True, exist_ok=True)
        shutil.copy(ckpt_path, copy_path)

max_stars_repo_name: lin-tan/fairness-variance

max_issues_repo_name: lin-tan/fairness-variance

max_forks_repo_name: lin-tan/fairness-variance


text:
"""
A trainer class.
"""

import torch,math
import torch.nn as nn
import torch.nn.functional as F
from torch.autograd import Variable
import numpy as np
from sklearn import metrics
from model.aggcn import GCNClassifier
from utils import torch_utils

class Trainer(object):
	def __init__(self, opt, emb_matrix=None):
		raise NotImplementedError

	def update(self, batch):
		raise NotImplementedError

	def predict(self, batch):
		raise NotImplementedError

	def update_lr(self, new_lr):
		torch_utils.change_lr(self.optimizer, new_lr)

	def load(self, filename,new_opt = None):
		try:
			checkpoint = torch.load(filename)
		except BaseException:
			print("Cannot load model from {}".format(filename))
			exit()
		self.model.load_state_dict(checkpoint['model'])
		self.opt = checkpoint['config']
		#try to load with new opt
		# if new_opt is not None and new_opt['cuda']:self.model = self.model.cuda()
		# self.opt = checkpoint['config'] if new_opt is None else new_opt

	def save(self, filename, epoch,best_dev_score):
		self.opt['best_dev_score']=best_dev_score
		params = {
				'model': self.model.state_dict(),
				'config': self.opt,
				}
		try:
			torch.save(params, filename)
			print("model saved to {}".format(filename))
		except BaseException:
			print("[Warning: Saving failed... continuing anyway.]")

def unpack_batch(batch, cuda):
	batch = [torch.from_numpy( data_goal ).float() for data_goal in batch ]

	if cuda:
		inputs = Variable(batch[0].cuda())
		goals = Variable( batch[-1].cuda() )
	else:
		inputs = Variable(batch[0])
		goals = Variable(batch[-1])

	return inputs,goals

class GCNTrainer(Trainer):
	def __init__(self, opt, emb_matrix=None):
		self.opt = opt
		self.emb_matrix = emb_matrix
		self.model = GCNClassifier(opt, emb_matrix=emb_matrix)
		self.criterion = nn.CrossEntropyLoss()
		self.parameters = [p for p in self.model.parameters() if p.requires_grad]
		if opt['cuda']:
			self.model.cuda()
			self.criterion.cuda()
		self.optimizer = torch_utils.get_optimizer(opt['optim'], self.parameters, opt['lr'])

	def update(self, batch):
		inputs,goals = unpack_batch(batch, self.opt['cuda'])

		# step forward
		self.model.train()
		self.optimizer.zero_grad()
		predict_goals = self.model(inputs).squeeze()

		Directions = predict_goals * goals
		UnCoDirections = Directions < 0

		loss_fn_RM = torch.nn.MSELoss()
		# loss_fn_RD = torch.nn.BCELoss()
		lossRM = loss_fn_RM(predict_goals,goals)
		loss = lossRM

		if UnCoDirections.sum() > 0:
			loss_UC = (predict_goals-goals)**2 *UnCoDirections.float()
			loss_UC = loss_UC.sum()/UnCoDirections.sum()
			loss_UC = loss_UC**0.5
			loss =lossRM + loss_UC
		print(predict_goals, goals,(predict_goals-goals)**2, UnCoDirections,UnCoDirections.sum(),lossRM,loss,sep='\n')

		loss.backward()
		torch.nn.utils.clip_grad_norm_(self.model.parameters(), self.opt['max_grad_norm'])
		self.optimizer.step()

		return loss.item()

	def predict(self, batch,result_needing=False):
		self.model.eval()
		inputs,goals = unpack_batch(batch, self.opt['cuda'])
		predict_goals = self.model(inputs).squeeze()
		Directions = predict_goals * goals
		CoDirections = Directions >0
		CoDirections_Accurate = sum(CoDirections).float()/CoDirections.flatten().shape[0]

		actions = (predict_goals>0).float()
		profit_loss_s = ( torch.prod( actions*goals/100+1 ) - 1 ) *100

		cover=RD_result=0

		res = {'CoDirections_Accurate': CoDirections_Accurate, 'profit_loss_s': profit_loss_s,
				 'cover': cover,'RD_result':RD_result }

		return res

if __name__ == '__main__':
	from utils import constant, torch_utils
	gcnTrainer = GCNTrainer(constant.default_opt)