text:
import numpy as np
from scipy.special import expit as sigmoid, logsumexp
from scipy.signal import convolve2d
from scipy.stats import norm as norm_dist
import sys, os
import copy
import matplotlib.pyplot as plt

from manipulator_learning.learning.data.tf.img_replay_buffer import ImgReplayBufferDisk
from manipulator_learning.learning.eval.plot.utils import get_roc_stats_qual, get_roc_stats

class FailurePredictorStatistical:
    MODE_OPTIONS = ['qual_per_step', 'cum_qual_full_traj', 'multistep_zero_qual', 'cum_qual_multistep',
                    'gda']

    def __init__(self, discriminator, replay_buffer, expert_replay_buffer,
                 d_model_step, ep_len, num_suc_fail_trajs_avg=25, d_model_step_delay=3000,
                 save_avg_d_fig_dir=None, fp_ignore_quantity=5, percent_change_qual=.05, num_prev_for_prior=25,
                 prev_ts_mean_filter=5, min_tpr_thresh_set=.8, mode='gda'):
        """
        A failure predictor that grounds its predictions based on the average predictions of previous successful
        and failed trajectories, given ground truth labels of success and failure.

        :param discriminator: Discriminator object.
        :param d_model_step: The d_model_step object.
        :param replay_buffer: object containing the replay_buffer
        :param expert_replay_buffer: object containing the expert replay buffer.
        :param num_suc_fail_trajs_avg: Number of previous data to use per timestep for both successful
          and failure data (e.g. 25 here would mean 50 total (s,a) pairs per timestep, or 50 total "trajs".
        :param d_model_step_delay: d_model_step must be higher than this before any predictions are made.
        """
        self.discriminator = discriminator
        self.d_model_step = d_model_step
        self.d_model_step_delay = d_model_step_delay
        self.rb = replay_buffer
        self.erb = expert_replay_buffer
        assert mode in FailurePredictorStatistical.MODE_OPTIONS, 'Invalid mode %s' % mode
        self.mode = mode
        if type(self.rb) == ImgReplayBufferDisk:
            self.obs_is_dict = True
        else:
            raise NotImplementedError('Implement if needed')
        self.ep_len = ep_len

        # indices marked with -1 in these buffers indicate no data yet. Otherwise, these arrays contain the equivalent
        # of num_suc_fail_trajs_avg trajectories of data, but each "trajectory" might not be from one contiguous traj
        # (since trajectories of autonomous data, in intervention learning, do not necessarily all come in at same time)
        #
        # the first index of the final dimension contains the actual data index, the second index contains a
        # 1 if the data comes from the expert buffer, and a 0 if it comes from the non-expert buffer
        self.suc_traj_indices = (np.ones([num_suc_fail_trajs_avg, ep_len, 2]) * -1).astype('int')
        self.fail_traj_indices = (np.ones([num_suc_fail_trajs_avg, ep_len, 2]) * -1).astype('int')

        # variables for keeping track of historical stat data
        self.d_suc_means = np.ones(ep_len) * -1000
        self.d_suc_stds = np.ones(ep_len) * -1000
        self.d_fail_means = np.ones(ep_len) * -1000
        self.d_fail_stds = np.ones(ep_len) * -1000
        self.hell_dist = np.ones(ep_len) * -1000
        self.indicator = np.zeros(ep_len)

        # variables for doing failure prediction
        self.qual_thresh = 0
        self._cur_traj_cum_qual = 0
        self._cur_consec_cond = 0
        self.fp_ignore_quantity = fp_ignore_quantity
        self._ignore_count = 0  # ignore a series of predictions after a false positive
        self._qual_change_amount = [0, 0]  # amount to change qual thresh by after incorrect predictions for [fp, fn]
        if self.mode == 'multistep_zero_qual':
            self._qual_change_amount = [-1, -3]
        self.percent_change_qual = percent_change_qual  # percent of the full range of possible values

        self.save_avg_d_fig_dir = save_avg_d_fig_dir

        self.user_feedback_given = False
        self.print_freq = 10
        self._next_print_count = self.print_freq
        self.latest_raw_d_value = None

        # gda (gaussian discriminant analysis) variables
        self.suc_dists = [None] * ep_len
        self.fail_dists = [None] * ep_len
        self.log_suc_prior = np.log(0.5)
        self.log_fail_prior = np.log(0.5)
        self.prev_results = np.zeros(num_prev_for_prior)
        self._prev_pred_outs = np.ones(prev_ts_mean_filter)
        self.min_tpr_thresh_set = min_tpr_thresh_set

        self.gda_ready = False

        self._always_predict_false = False

    def append_traj(self, successful, expert_rb, traj_data_indices, timesteps, add_to_prev_suc_results=False):
        """
        Append a traj of data indices to the historical data.

        :param successful: Whether the trajectory was successful or not. If only adding a partial trajectory
          if it's partial because an expert intervened, counts as an unsuccessful traj.
        :param: Whether the trajectory is stored in the expert replay buffer. if false, stored in non-expert rb.
        :param traj_data_indices: The trajectory of data indices to add.
        :param timesteps: The t values for each piece of new data.
        :param add_to_prev_suc_results: whether to count this trajectory towards the prev suc results. If this is
          being called for only a small amount of data at the end of a traj, or as a result of an expert demo, this
          should be False.
        :return:
        """
        if len(traj_data_indices) != len(timesteps):
            print("traj data indices: ", traj_data_indices)
            print("timesteps: ", timesteps)
            import ipdb;
            ipdb.set_trace()

        # assert len(traj_data_indices) == len(timesteps), "Num new data indices %d does not match num of " \
        #                                                  "new timesteps %d" % (len(traj_data_indices), len(timesteps))
        if len(traj_data_indices) == 0:
            return

        if successful:
            # only roll the timesteps for which we have new data
            self.suc_traj_indices[:, timesteps, :] = np.roll(self.suc_traj_indices[:, timesteps], -1, axis=0)
            self.suc_traj_indices[-1, timesteps, 0] = traj_data_indices
            self.suc_traj_indices[-1, timesteps, 1] = expert_rb
        else:
            self.fail_traj_indices[:, timesteps, :] = np.roll(self.fail_traj_indices[:, timesteps], -1, axis=0)
            self.fail_traj_indices[-1, timesteps, 0] = traj_data_indices
            self.fail_traj_indices[-1, timesteps, 1] = expert_rb

        if add_to_prev_suc_results:
            self.prev_results = np.roll(self.prev_results, -1)
            self.prev_results[-1] = successful
            self.update_priors()

    def update_priors(self):
        suc_prop = self.prev_results.sum() / len(self.prev_results)
        fail_prop = 1 - suc_prop
        suc_prior = np.clip(suc_prop, .05, .95)
        fail_prior = np.clip(fail_prop, .05, .95)
        self.log_suc_prior = np.log(suc_prior)
        self.log_fail_prior = np.log(fail_prior)

    def discrim_as_batches(self, rb_data, rb, batch_size=512):
        num_data = rb_data[0].shape[0]
        num_batches = num_data // batch_size + 1
        rb_d_values = np.zeros([num_data, 1])
        for i in range(num_batches):
            sl = slice(i * batch_size, min((i + 1) * batch_size, num_data))
            batch_data = tuple([rb_data[d_ind][sl] for d_ind in range(len(rb_data))])
            if len(rb_d_values[sl]) > 0:
                if self.discriminator.num_ensemble == 1:
                    rb_d_values[sl] = self.discriminator(rb.tuple_batch_to_d_input(batch_data)).numpy()
                else:
                    rb_d_values[sl] = np.expand_dims(
                        self.discriminator.get_mean_discrim(rb.tuple_batch_to_d_input(batch_data)).numpy(), -1)

        return rb_d_values

    def calc_d_traj_avgs(self, total_numsteps):
        """ Calculate new statistics. Should be called everytime the Discriminator is updated before trying to
            predict failuures. """
        for traj_indices in [self.suc_traj_indices, self.fail_traj_indices]:
            rb_valid_indices = traj_indices[:, :, 1] == 0
            rb_indices = traj_indices[rb_valid_indices][:, 0]
            rb_data = self.rb.dataset.get_data(rb_indices.astype('int'))

            # instead of running all data thru discrim at the same time, run it through in batches to not overuse gpu mem
            # if batch_size is too high, uses too much mem, too low, too slow
            rb_d_values = self.discrim_as_batches(rb_data, self.rb)

            erb_valid_indices = traj_indices[:, :, 1] == 1
            erb_indices = traj_indices[erb_valid_indices][:, 0]
            if len(erb_indices) > 0:
                erb_data = self.erb.dataset.get_data(erb_indices.astype('int'))
                erb_d_values = self.discrim_as_batches(erb_data, self.erb)
            else:
                erb_d_values = np.array([])

            d_values = np.ones(traj_indices.shape[:2]) * -1
            d_values[rb_valid_indices] = rb_d_values.flatten()
            d_values[erb_valid_indices] = erb_d_values.flatten()
            d_values_sig = sigmoid(d_values)
            if traj_indices is self.suc_traj_indices:
                d_suc_values = copy.deepcopy(d_values)
                d_suc_values_sig = copy.deepcopy(d_values_sig)
                self.d_suc_means = d_values_sig.mean(axis=0)
                self.d_suc_stds = d_values_sig.std(axis=0)
                d_suc_means_sig = d_values_sig.mean(axis=0)
                d_suc_stds_sig = d_values_sig.std(axis=0)
            elif traj_indices is self.fail_traj_indices:
                d_fail_values = copy.deepcopy(d_values)
                d_fail_values_sig = copy.deepcopy(d_values_sig)
                self.d_fail_means = d_values_sig.mean(axis=0)
                self.d_fail_stds = d_values_sig.std(axis=0)
                d_fail_means_sig = d_values_sig.mean(axis=0)
                d_fail_stds_sig = d_values_sig.std(axis=0)

        if self.save_avg_d_fig_dir is not None:
            avg_d_fig = plt.figure()
            avg_d_fig_ax = avg_d_fig.add_subplot(111)
            ep_len = self.suc_traj_indices.shape[1]
            avg_d_fig_ax.axhline(.5, linestyle='--', color='black')
            avg_d_fig_ax.plot(range(ep_len), d_suc_means_sig, label='Successes')
            avg_d_fig_ax.plot(range(ep_len), d_fail_means_sig, label='Failures')
            avg_d_fig_ax.fill_between(range(ep_len), d_suc_means_sig - d_suc_stds_sig, d_suc_means_sig + d_suc_stds_sig,
                                      alpha=.1)
            avg_d_fig_ax.fill_between(range(ep_len), d_fail_means_sig - d_fail_stds_sig,
                                      d_fail_means_sig + d_fail_stds_sig,
                                      alpha=.1)
            avg_d_fig_ax.legend(fontsize=14, loc='lower right')
            # plt.show()
            os.makedirs(self.save_avg_d_fig_dir, exist_ok=True)
            avg_d_fig.savefig(self.save_avg_d_fig_dir + '/d_avg_ckpt_' + str(total_numsteps) + '.png')
            plt.close(avg_d_fig)

        if self.mode == 'gda':
            self.suc_dists = [norm_dist(m, s) for m, s in zip(self.d_suc_means, self.d_suc_stds)]
            self.fail_dists = [norm_dist(m, s) for m, s in zip(self.d_fail_means, self.d_fail_stds)]
            self.update_priors()

            # vectorized dists for getting optimal thresh
            vec_suc_dists = norm_dist(self.d_suc_means, self.d_suc_stds)
            vec_fail_dists = norm_dist(self.d_fail_means, self.d_fail_stds)
            all_d = np.concatenate([d_suc_values_sig, d_fail_values_sig])
            suc_log_probs = np.log(vec_suc_dists.pdf(all_d))
            fail_log_probs = np.log(vec_fail_dists.pdf(all_d))
            log_ds = logsumexp([suc_log_probs + self.log_suc_prior, fail_log_probs + self.log_fail_prior], axis=0)
            per_step_pred_suc = np.exp(suc_log_probs + self.log_suc_prior - log_ds)
            conv_op = np.ones([1, len(self._prev_pred_outs)]) * 1 / len(self._prev_pred_outs)
            conv_out = convolve2d(per_step_pred_suc, conv_op, fillvalue=1.0)[:, :self.ep_len]
            conv_out_min = conv_out.min(axis=1)
            num_suc_trajs = d_suc_values_sig.shape[0]
            num_fail_trajs = d_fail_values_sig.shape[0]
            test_values = np.linspace(conv_out_min.min() - .01, conv_out_min.max() + .01,
                                      (num_suc_trajs + num_fail_trajs) * 2)
            successes = np.concatenate([np.ones(num_suc_trajs), np.zeros(num_fail_trajs)], axis=0).astype('bool')
            roc_points = get_roc_stats_qual(successes, conv_out_min, test_values)
            acceptable_indices = np.argwhere(roc_points[:, 1] > self.min_tpr_thresh_set)
            acceptable_best_thresh = test_values[acceptable_indices][np.linalg.norm(
                roc_points - np.array([0, 1]), axis=1)[acceptable_indices].argmin()]
            self._qual_change_amount = [self.percent_change_qual * (conv_out_min.max() - conv_out_min.min())] * 2
            self.qual_thresh = acceptable_best_thresh

        # Hellinger distance -- see https://en.wikipedia.org/wiki/Hellinger_distance
        elif self.mode in ['cum_qual_full_traj', 'qual_per_step', 'multistep_zero_qual']:
            self.hell_dist = np.sqrt(1 - (np.sqrt(2 * self.d_suc_stds * self.d_fail_stds /
                                                  (self.d_suc_stds ** 2 + self.d_fail_stds ** 2)) *
                                          np.exp(-.25 * (self.d_suc_means - self.d_fail_means) ** 2 /
                                                 (self.d_suc_stds ** 2 + self.d_fail_stds ** 2))))
            self.indicator = (self.d_suc_means > self.d_fail_means)

            # set the threshold optimally based on best ROC point
            if not self.user_feedback_given:
                qual_per_ts_suc, qual_cum_suc, qual_per_ts_min_suc = self.get_quality(d_suc_values_sig)
                qual_per_ts_fail, qual_cum_fail, qual_per_ts_min_fail = self.get_quality(d_fail_values_sig)

                num_suc_trajs = d_suc_values_sig.shape[0]
                num_fail_trajs = d_fail_values_sig.shape[0]
                if self.mode == 'cum_qual_full_traj':
                    qual = np.concatenate([qual_cum_suc, qual_cum_fail], axis=0)
                elif self.mode == 'qual_per_step':
                    qual = np.concatenate([qual_per_ts_min_suc, qual_per_ts_min_fail], axis=0)
                elif self.mode == 'multistep_zero_qual':
                    qual_per_ts = np.concatenate([qual_per_ts_suc, qual_per_ts_fail], axis=0)
                    qual_lt_zero = (qual_per_ts < 0)
                    qual = self.get_max_consecutive_cond(qual_lt_zero)

                qual_test_values = np.linspace(qual.min() - .1, qual.max() + .1,
                                               (num_suc_trajs + num_fail_trajs) * 2)
                successes = np.concatenate([np.ones(num_suc_trajs), np.zeros(num_fail_trajs)], axis=0).astype('bool')

                if self.mode == 'multistep_zero_qual':
                    qual_roc_points = get_roc_stats_qual(successes, qual, qual_test_values,
                                                         qual_should_be_lt_param=False)
                    self._qual_change_amount = [-1, -3]  # since we use - and + for fp and fn respectively
                    self.qual_thresh = int(
                        qual_test_values[np.linalg.norm(qual_roc_points - np.array([0, 1]), axis=1).argmin()])
                elif self.mode in ['cum_qual_full_traj', 'qual_per_step']:
                    qual_roc_points = get_roc_stats_qual(successes, qual, qual_test_values)
                    self._qual_change_amount = [self.percent_change_qual * (qual.max() - qual.min()),
                                                self.percent_change_qual * (qual.max() - qual.min())]
                    self.qual_thresh = qual_test_values[
                        np.linalg.norm(qual_roc_points - np.array([0, 1]), axis=1).argmin()]

        self.gda_ready = True

    def get_max_consecutive_cond(self, cond):
        """ Get the maximum consecutive condition for a set of trajectories. cond should be a (N, ep_len) array
            of booleans indicating the condition."""
        # this is implemented with an outer for loop, but it's such a simple task that it shouldn't significantly slow
        # anything down
        max_consecs = []
        for traj in cond:
            diff = np.diff(traj)
            idx, = diff.nonzero()
            idx += 1
            if traj[0]:
                idx = np.r_[0, idx]
            if traj[-1]:
                idx = np.r_[idx, traj.size]
            idx = idx.reshape(-1, 2)
            lens = idx[:, 1] - idx[:, 0]
            if len(lens) == 0:
                max_consecs.append(0.)
            else:
                max_consecs.append(lens.max())
        return np.array(max_consecs)

    def get_quality(self, d_values):
        """ d_values should be a (N, ep_len) array"""
        qual_per_ts = (self.indicator * self.hell_dist * ((d_values - self.d_suc_means) / self.d_suc_stds +
                                                          (d_values - self.d_fail_means) / self.d_fail_stds))
        qual_cum = qual_per_ts.sum(axis=1)
        qual_per_ts_min = qual_per_ts.min(axis=1)

        return qual_per_ts, qual_cum, qual_per_ts_min

    def get_quality_single_step(self, d_value, t):
        """ Quality for a single d value. needs to have associated timestep. """
        return (self.indicator[t] * self.hell_dist[t] * ((d_value - self.d_suc_means[t]) / self.d_suc_stds[t] +
                                                         (d_value - self.d_fail_means[t]) / self.d_fail_stds[t]))

    def get_gda_pred_single_step(self, d_value, t):
        """ GDA prediction for a single d value and timestep. """
        if not self.gda_ready:
            return 0.5
        suc_log_prob = np.log(self.suc_dists[t].pdf(d_value))
        fail_log_prob = np.log(self.fail_dists[t].pdf(d_value))
        log_d = logsumexp([suc_log_prob + self.log_suc_prior, fail_log_prob + self.log_fail_prior], axis=0)
        pred_suc = np.exp(suc_log_prob + self.log_suc_prior - log_d)
        return pred_suc

    def get_discrim_output(self, obs, action):
        raw_discrim = self.discriminator(self.rb.obs_act_to_d_input(obs, action)).numpy().flatten()
        self.latest_raw_d_value = raw_discrim
        if self.discriminator.num_ensemble > 1:
            return np.mean(raw_discrim)
        else:
            return raw_discrim

    def get_failure_prediction(self, obs, action, t):

        d_out = self.get_discrim_output(obs, action)
        d_out_sig = sigmoid(d_out)
        if self.mode in ['cum_qual_full_traj', 'qual_per_step', 'multistep_zero_qual']:
            qual = self.get_quality_single_step(d_out_sig, t)
        elif self.mode == 'gda':
            qual = self.get_gda_pred_single_step(d_out_sig, t)
        if len(self._prev_pred_outs) > 1:
            self._prev_pred_outs = np.roll(self._prev_pred_outs, -1)
            self._prev_pred_outs[-1] = qual
            qual = np.mean(self._prev_pred_outs)
        if self.mode == 'multistep_zero_qual':
            if qual < 0:
                self._cur_consec_cond += 1  # still update consec values even if ignore count is above 0
            else:
                self._cur_consec_cond = 0
        if self._ignore_count > 0:
            self._ignore_count -= 1
            return False

        # print outputs
        if self._next_print_count == 0:
            if self.mode == 'multistep_zero_qual':
                print('Qual t: %.3f, cur consec: %d, qual thresh: %d' % (qual, self._cur_consec_cond, self.qual_thresh))
            elif self.mode == 'gda':
                print('pred suc t: %.4f, thresh: %.4f' % (qual, self.qual_thresh))
            self._next_print_count = self.print_freq
        else:
            self._next_print_count -= 1

        if self._always_predict_false:
            return False

        if self.mode in ['cum_qual_full_traj', 'qual_per_step', 'gda']:
            if self.mode == 'cum_qual_full_traj':
                self._cur_traj_cum_qual += qual
                qual_test = self._cur_traj_cum_qual
                # print('Qual t: %.3f, cum qual: %.3f, qual thresh: %.3f' % (qual, self._cur_traj_cum_qual, self.qual_thresh))
            elif self.mode in ['qual_per_step', 'gda']:
                qual_test = qual
            return (True if qual_test < self.qual_thresh else False)
        if self.mode == 'multistep_zero_qual':
            return (True if self._cur_consec_cond >= self.qual_thresh else False)

    def false_positive(self):
        """ Called by user after a false positive. """
        self._ignore_count = self.fp_ignore_quantity
        self.qual_thresh -= self._qual_change_amount[0]
        self.user_feedback_given = True

    def false_negative(self):
        """ Called by the user after a false negative """
        self.qual_thresh += self._qual_change_amount[1]
        self.user_feedback_given = True

    def new_ep_reset(self):
        self._ignore_count = 0
        self._cur_traj_cum_qual = 0
        self._cur_consec_cond = 0
        self._next_print_count = self.print_freq
        self._prev_pred_outs.fill(1)

    def set_always_predict_false(self, new_setting):
        self._always_predict_false = new_setting

reason: import numpy,from scipy


text:
import tensorflow as tf
import tensorflow.contrib.eager as tfe
import numpy as np
tf.compat.v1.enable_eager_execution()

# 可以直接计算
x = [[2.]]
m = tf.matmul(x, x)
print(m)

# 也不用担心 lazy loading
x = tf.random.uniform([2, 2])
for i in range(x.shape[0]):
    for j in range(x.shape[1]):
        print(x[i, j])

# 可以像 Numpy Array 一样用
x = tf.constant([1.0, 2.0, 3.0])
assert type(x.numpy()) == np.ndarray
squared = np.square(x)
# Tensors 可以遍历
for i in x:
    print(i)

# 也可以直接进行梯度计算
def square(x):
    return x**2

grad = tfe.gradients_function(square)
print(square(3.))
print(grad(3.))

x = tf.Variable(2.0)
def loss(y):
    return (y - x ** 2) ** 2

grad = tfe.implicit_gradients(loss)
print(loss(7.))
print(grad(7.))

ext: py

max_stars_repo_name: wdxtub/deep-learning-note

max_issues_repo_name: wdxtub/deep-learning-note

max_forks_repo_name: wdxtub/deep-learning-note

reason: import numpy


text:
import numpy as np
import sys
import pandas as pd
from pathlib import Path
import matplotlib as mpl
from matplotlib import pyplot as plt
import stat_tools as st
from datetime import datetime
from scipy import ndimage
from scipy.optimize import minimize
from scipy.optimize import least_squares
import ephem
import configparser as cfg
import yaml
import camcoord
# globals
pi2=np.pi/2.

# Read observed Moon position output from find_moon.py and optimize
# camera parameters by minimizing the mean square angular distance to
# the predicted location from ephem.Moon

# Initial parameters are read from camera_cal_bnl.yaml
#####params: nx0,cy,cx,rotation,beta,azm,c1,c2,c3
# ny0=nx0 are the y and x size of the region of interest 'roi', assumed square
# cy,cx is the central pixel of 'roi'
# rotation is the deviation from North in radian, positive to East.

def meansqdist(x):
    # function to be minimized mean square of angular distance between measured
    # and predicted moon positions
    # vars = string array of parameters to vary in x
    # start with just varying cam.rot
    cam.pr0=x[0]
    cam.cy=x[1]
    cam.cx=x[2]
    cam.rot=x[3]
    cam.beta=x[4]
    cam.azm=x[5]
    cam.c1=x[6]
    cam.c2=x[7]
    cam.c3=(0.5-(cam.c1*pi2+cam.c2*pi2**3))/pi2**5
    cam.azizen()
    ix=(dfref.xmoon+0.5).astype(int)
    iy=(dfref.ymoon+0.5).astype(int)
    otheta=cam.theta0[iy,ix]
    ophi=cam.phi0[iy,ix]
    ptheta=[]
    pphi=[]
    for edate in dfref.ephemDate:
        obs.date=edate
        moon.compute(obs)
        ptheta.append(np.pi/2.-moon.alt)
        pphi.append((moon.az-np.pi)%(2*np.pi))
    ptheta=np.array(ptheta)
    pphi=np.array(pphi)
    return np.mean(camcoord.great_circle_distance(otheta,ophi,ptheta,pphi)**2)

def dist(x):
    # function returns residuals to be
    # minimized by least-squares angular distance between measured
    # and predicted moon positions
    # vars = string array of parameters to vary in x
    # start with just varying cam.rot
    cam.pr0=x[0]
    cam.cy=x[1]
    cam.cx=x[2]
    cam.rot=x[3]
    cam.beta=x[4]
    cam.azm=x[5]
    cam.c1=x[6]
    cam.c2=x[7]
    if constrained_c3:
        cam.c3=(0.5-(cam.c1*pi2+cam.c2*pi2**3))/pi2**5
    else:
        cam.c3=x[8]
    cam.azizen()
    ix=(dfref.xmoon+0.5).astype(int)
    iy=(dfref.ymoon+0.5).astype(int)
    otheta=cam.theta0[iy,ix]
    ophi=cam.phi0[iy,ix]
    try:
        exist=(len(ptheta) == len(otheta))
    except:
        exist=False
    if not exist:
        ptheta=[]
        pphi=[]
        for edate in dfref.ephemDate:
            obs.date=edate
            moon.compute(obs)
            ptheta.append(np.pi/2.-moon.alt)
            pphi.append((moon.az-np.pi)%(2*np.pi))
        ptheta=np.array(ptheta)
        pphi=np.array(pphi)
    return camcoord.great_circle_distance(otheta,ophi,ptheta,pphi)

if __name__ == "__main__":
    ######load the configuration file
    config_file = sys.argv[1] if len(sys.argv) >= 2 else 'camera_calibration.conf'
    config = cfg.ConfigParser()
    config.read(config_file)

    if len(sys.argv) >= 3:
        cameraIDs = sys.argv[2:]
    else:
        cameraIDs = eval(config['camera']['cameraIDs'])
    imagepath = config['path']['imagepath']
    outpath = config['path']['outpath']
    moon_obs_ext = config['path']['moon_obs_ext']
    camera_cal_file_list = config['path']['camera_cal_file_list']
    camera_cal_file_optimized = config['path']['camera_cal_file_optimized']
    begin_date=ephem.date(config['optimization']['begin_date'])
    end_date=ephem.date(config['optimization']['end_date'])
    constrained_c3=eval(config['optimization']['constrained_c3'])
    sitelat = float(config['geolocation']['lat'])
    sitelon = float(config['geolocation']['lon'])
    # big outer loop
    for cameraID in cameraIDs:
        # initial camera parameters from camera_cal_file_optimized
        # which should be a copy of camera_cal_file
        cam=camcoord.camera(cameraID,camera_cal_file=camera_cal_file_optimized)
    # calculate azi and zen for each pixel, to optimize vary camera parameters
    # rot,cx,cy,nr0,beta,azm,c1,c2, and optionally c3
        # (if constrained_c3 is False) and recalculate.
        cam.azizen()
        nx0=ny0=cam.nx0
        # pr0 is nx0/2, i.e. initial radius estimate.
        pr0=cam.pr0
        # obs is an ephem observer object
        obs = ephem.Observer();
        # lat, lon are the only parameters in config file specified in deg.
        # ephem and numpy use radian
        obs.lat = np.deg2rad(cam.lat)
        obs.lon = np.deg2rad(cam.lon)
        # moon is an ephem moon object
        moon=ephem.Moon()
    #
    # read moon position data from find_moon_all output .csv files between
    # begin_date and end_date
        moonobsdir=Path(outpath,cam.camID)
        moonobsfiles=sorted(moonobsdir.glob("*"+moon_obs_ext))
        if len(moonobsfiles) < 1:
            print("no moon obs found in {}, run find_moon_all.py first".format(moonobsdir))
            exit(1)
        for moonobsfile in moonobsfiles:
            datestr=moonobsfile.name.split('_')
            d1,d2=datestr[1].split('-')
            d1=ephem.date(datetime.strptime(d1,'%Y%m%d').strftime('%Y-%m-%d'))
            d2=ephem.date(datetime.strptime(d2,'%Y%m%d').strftime('%Y-%m-%d'))
            if d1>end_date or d2<begin_date:
                continue
            # read into pandas DataFrame
            dfrefi=pd.read_csv(moonobsfile,index_col=0)
            # only use good points
            # i.e. roundness check and minimum radius 3 pixels.
            good=((dfrefi.flag & 512) == 0) & (dfrefi.rmoon > 3.)
            dfrefi=dfrefi[good]
            if len(dfrefi) < 1:
                print("no good data points, skipping "+moonobsfile.name)
                continue
            # now concatenate
            try:
                print('{}: Adding {} records to {}'.format(moonobsfile.name,len(dfrefi),len(dfref)))
                dfref=dfref.append(dfrefi)
            except:
                dfref=dfrefi
    # now minimize mean square distance, second by varying the 3 camera orientation
    # angles
        # constraint optimization pr0,cy,cx should not vary by more than 1%
        ptol=0.99
        # rot should be small -0.5,0.5 (about +-30 deg)
        # beta is the tilt and should be small positive 0,0.2 (11 deg)
        # azm is -pi,pi
        # c1, c2, c3 I am not sure off, but they cannot be 0.
        if constrained_c3:
            x0=np.array([cam.pr0,cam.cy,cam.cx,cam.rot,cam.beta,cam.azm,cam.c1,cam.c2])
            bounds=([cam.pr0*ptol,cam.cy*ptol,cam.cx*ptol,-np.pi,0.,-np.pi,-1.,-1.],
                    [cam.pr0/ptol,cam.cy/ptol,cam.cx/ptol,np.pi,0.2,np.pi,1.,1.])
        else:
            x0=np.array([cam.pr0,cam.cy,cam.cx,cam.rot,cam.beta,cam.azm,cam.c1,cam.c2,cam.c3])
            bounds=([cam.pr0*ptol,cam.cy*ptol,cam.cx*ptol,-np.pi,0.,-np.pi,-1.,-1.,-1.],
                    [cam.pr0/ptol,cam.cy/ptol,cam.cx/ptol,np.pi,0.2,np.pi,1.,1.,1.])
#
        print('before',x0)
        print('with bounds',bounds)
#        res=minimize(meansqdist,x0,method='BFGS',options={'disp': True})
        res=least_squares(dist,x0,verbose=2,bounds=bounds)
        print('after',res.x)
        if constrained_c3:
        # append the constrained c3
            c1=res.x[6]
            c2=res.x[7]
            res.x=np.append(res.x,(0.5-(c1*pi2+c2*pi2**3))/pi2**5)
        # Note: res.x is type float64, which is not one of the canonical
        # guessable data-types and hence yaml.dump adds additional "!!" tags
        # to fully describe the dumped object, which is ugly.
        # converting to float appears to resolve this make sure that's done in
        # these "*_to_yaml" methods.
        cam.save_dict_to_yaml(res.x[:],camera_cal_file_optimized)
        cam.save_list_to_yaml(res.x[:],camera_cal_file_list)


text:
import numpy as np

def P_y_x(x, mean_1, mean_2, stddev_1, stddev_2, P_1=0.5, P_2=0.5):
    """
    Assume P(Y) is uniform i.e. P(Y=c1) = 0.5 and P(Y=c2) = 0.5.
    P(Y=c1|X=x) = P(X=x|Y=c1)P(Y=c1) / P(X=x), P(Y=c2|X=x) = P(X=x|Y=c2)P(Y=c2) / P(X=x)
    :param x: the sample taken form normal
    :param mean_1: mean 1
    :param mean_2: mean 2
    :param stddev_1: standard deviation 1
    :param stddev_2: standard deviation 2
    :param P_1: P(Y=c1)
    :param P_2: P(Y=c2)
    :return: (P(Y=c1|X=x), P(Y=c2|X=x))
    """
    P_X_1 = 1 / (stddev_1 * np.sqrt(2 * np.pi)) * np.exp(- (x - mean_1) ** 2 / (2 * stddev_1 ** 2))
    P_X_2 = 1 / (stddev_2 * np.sqrt(2 * np.pi)) * np.exp(- (x - mean_2) ** 2 / (2 * stddev_2 ** 2))
    return P_X_1 * P_1 / (P_X_1 * P_1 + P_X_2 * P_2), P_X_2 * P_2 / (P_X_1 * P_1 + P_X_2 * P_1)

### NEW 2D ###
def P_y_x_2d(N_1, N_2):
    P_1 = []
    P_2 = []
    for i in range(len(N_1)):
        n_1 = N_1[i]
        n_2 = N_2[i]
        P_1.append(n_1 / (n_1 + n_2))
        P_2.append(n_2 / (n_1 + n_2))
    return P_1, P_2

max_stars_repo_name: adamismaili/MissClassificationError

max_issues_repo_name: adamismaili/MissClassificationError

max_forks_repo_name: adamismaili/MissClassificationError


text:
import copy
import os
from functools import partial

import joblib
import numpy as np
import optuna
import pandas as pd
import xgboost as xgb
from loguru import logger
from sklearn import metrics

from .enums import ProblemType
from .params import get_params

def save_valid_predictions(final_valid_predictions, model_config, target_encoder, output_file_name):
    final_valid_predictions = pd.DataFrame.from_dict(final_valid_predictions, orient="index").reset_index()
    if target_encoder is None:
        final_valid_predictions.columns = [model_config.id_column] + model_config.target_cols
    else:
        final_valid_predictions.columns = [model_config.id_column] + list(target_encoder.classes_)

    final_valid_predictions.to_csv(
        os.path.join(model_config.output_dir, output_file_name),
        index=False,
    )

def save_test_predictions(final_test_predictions, model_config, target_encoder, test_ids, output_file_name):
    final_test_predictions = np.mean(final_test_predictions, axis=0)
    if target_encoder is None:
        final_test_predictions = pd.DataFrame(final_test_predictions, columns=model_config.target_cols)
    else:
        final_test_predictions = pd.DataFrame(final_test_predictions, columns=list(target_encoder.classes_))
    final_test_predictions.insert(loc=0, column=model_config.id_column, value=test_ids)
    final_test_predictions.to_csv(
        os.path.join(model_config.output_dir, output_file_name),
        index=False,
    )

def fetch_xgb_model_params(model_config):
    if model_config.problem_type == ProblemType.binary_classification:
        metric = metrics.log_loss
        xgb_model = xgb.XGBClassifier
        use_predict_proba = True
        direction = "minimize"
        eval_metric = "logloss"
    elif model_config.problem_type == ProblemType.multi_class_classification:
        metric = metrics.log_loss
        xgb_model = xgb.XGBClassifier
        use_predict_proba = True
        direction = "minimize"
        eval_metric = "mlogloss"
    elif model_config.problem_type == ProblemType.multi_label_classification:
        metric = metrics.log_loss
        xgb_model = xgb.XGBClassifier
        use_predict_proba = True
        direction = "minimize"
        eval_metric = "logloss"
    elif model_config.problem_type == ProblemType.single_column_regression:
        metric = partial(metrics.mean_squared_error, squared=False)
        xgb_model = xgb.XGBRegressor
        use_predict_proba = False
        direction = "minimize"
        eval_metric = "rmse"
    elif model_config.problem_type == ProblemType.multi_column_regression:
        metric = partial(metrics.mean_squared_error, squared=False)
        xgb_model = xgb.XGBRegressor
        use_predict_proba = False
        direction = "minimize"
        eval_metric = "rmse"
    else:
        raise NotImplementedError

    return xgb_model, use_predict_proba, eval_metric, metric, direction

def optimize(
    trial,
    xgb_model,
    metric,
    use_predict_proba,
    eval_metric,
    model_config,
):
    params = get_params(trial, model_config)
    early_stopping_rounds = params["early_stopping_rounds"]
    del params["early_stopping_rounds"]

    scores = []

    for fold in range(model_config.num_folds):
        train_feather = pd.read_feather(os.path.join(model_config.output_dir, f"train_fold_{fold}.feather"))
        valid_feather = pd.read_feather(os.path.join(model_config.output_dir, f"valid_fold_{fold}.feather"))
        xtrain = train_feather[model_config.features]
        xvalid = valid_feather[model_config.features]

        ytrain = train_feather[model_config.target_cols].values
        yvalid = valid_feather[model_config.target_cols].values

        # train model
        model = xgb_model(
            random_state=model_config.seed,
            eval_metric=eval_metric,
            use_label_encoder=False,
            **params,
        )

        if model_config.problem_type in (ProblemType.multi_column_regression, ProblemType.multi_label_classification):
            ypred = []
            models = [model] * len(model_config.target_cols)
            for idx, _m in enumerate(models):
                _m.fit(
                    xtrain,
                    ytrain[:, idx],
                    early_stopping_rounds=early_stopping_rounds,
                    eval_set=[(xvalid, yvalid[:, idx])],
                    verbose=False,
                )
                if model_config.problem_type == ProblemType.multi_column_regression:
                    ypred_temp = _m.predict(xvalid)
                else:
                    ypred_temp = _m.predict_proba(xvalid)[:, 1]
                ypred.append(ypred_temp)
            ypred = np.column_stack(ypred)

        else:
            model.fit(
                xtrain,
                ytrain,
                early_stopping_rounds=early_stopping_rounds,
                eval_set=[(xvalid, yvalid)],
                verbose=False,
            )

            if use_predict_proba:
                ypred = model.predict_proba(xvalid)
            else:
                ypred = model.predict(xvalid)

        # calculate metric
        metric_value = metric(yvalid, ypred)
        scores.append(metric_value)

    return sum(scores) / len(scores)

def train_model(model_config):
    xgb_model, use_predict_proba, eval_metric, metric, direction = fetch_xgb_model_params(model_config)

    optimize_func = partial(
        optimize,
        xgb_model=xgb_model,
        metric=metric,
        use_predict_proba=use_predict_proba,
        eval_metric=eval_metric,
        model_config=model_config,
    )
    db_path = os.path.join(model_config.output_dir, "params.db")
    study = optuna.create_study(
        direction=direction,
        study_name="autoxgb",
        storage=f"sqlite:///{db_path}",
        load_if_exists=True,
    )
    study.optimize(optimize_func, n_trials=model_config.num_trials, timeout=model_config.time_limit)
    return study.best_params

def predict_model(model_config, best_params):

    early_stopping_rounds = best_params["early_stopping_rounds"]
    del best_params["early_stopping_rounds"]

    xgb_model, use_predict_proba, eval_metric, metric, _ = fetch_xgb_model_params(model_config)
    scores = []

    final_test_predictions = []
    final_valid_predictions = {}

    target_encoder = joblib.load(f"{model_config.output_dir}/axgb.target_encoder")

    for fold in range(model_config.num_folds):
        train_feather = pd.read_feather(os.path.join(model_config.output_dir, f"train_fold_{fold}.feather"))
        valid_feather = pd.read_feather(os.path.join(model_config.output_dir, f"valid_fold_{fold}.feather"))

        xtrain = train_feather[model_config.features]
        xvalid = valid_feather[model_config.features]

        valid_ids = valid_feather[model_config.id_column].values

        if model_config.test_filename is not None:
            test_feather = pd.read_feather(os.path.join(model_config.output_dir, f"test_fold_{fold}.feather"))
            xtest = test_feather[model_config.features]
            test_ids = test_feather[model_config.id_column].values

        ytrain = train_feather[model_config.target_cols].values
        yvalid = valid_feather[model_config.target_cols].values

        model = xgb_model(
            random_state=model_config.seed,
            eval_metric=eval_metric,
            use_label_encoder=False,
            **best_params,
        )

        if model_config.problem_type in (ProblemType.multi_column_regression, ProblemType.multi_label_classification):
            ypred = []
            test_pred = []
            trained_models = []
            for idx in range(len(model_config.target_cols)):
                _m = copy.deepcopy(model)
                _m.fit(
                    xtrain,
                    ytrain[:, idx],
                    early_stopping_rounds=early_stopping_rounds,
                    eval_set=[(xvalid, yvalid[:, idx])],
                    verbose=False,
                )
                trained_models.append(_m)
                if model_config.problem_type == ProblemType.multi_column_regression:
                    ypred_temp = _m.predict(xvalid)
                    if model_config.test_filename is not None:
                        test_pred_temp = _m.predict(xtest)
                else:
                    ypred_temp = _m.predict_proba(xvalid)[:, 1]
                    test_pred_temp = _m.predict_proba(xtest)[:, 1]

                ypred.append(ypred_temp)
                test_pred.append(test_pred_temp)

            ypred = np.column_stack(ypred)
            test_pred = np.column_stack(test_pred)
            joblib.dump(
                trained_models,
                os.path.join(
                    model_config.output_dir,
                    f"axgb_model.{fold}",
                ),
            )

        else:
            model.fit(
                xtrain,
                ytrain,
                early_stopping_rounds=early_stopping_rounds,
                eval_set=[(xvalid, yvalid)],
                verbose=False,
            )

            joblib.dump(
                model,
                os.path.join(
                    model_config.output_dir,
                    f"axgb_model.{fold}",
                ),
            )

            if use_predict_proba:
                ypred = model.predict_proba(xvalid)
                if model_config.test_filename is not None:
                    test_pred = model.predict_proba(xtest)
            else:
                ypred = model.predict(xvalid)
                if model_config.test_filename is not None:
                    test_pred = model.predict(xtest)

        final_valid_predictions.update(dict(zip(valid_ids, ypred)))
        if model_config.test_filename is not None:
            final_test_predictions.append(test_pred)

        # calculate metric
        metric_value = metric(yvalid, ypred)
        scores.append(metric_value)

    save_valid_predictions(final_valid_predictions, model_config, target_encoder, "oof_predictions.csv")

    if model_config.test_filename is not None:
        save_test_predictions(final_test_predictions, model_config, target_encoder, test_ids, "test_predictions.csv")
    else:
        logger.info("No test data supplied. Only OOF predictions were generated")

max_stars_repo_name: toshihikoyanase/autoxgb

max_issues_repo_name: toshihikoyanase/autoxgb

max_forks_repo_name: toshihikoyanase/autoxgb


text:
"""
This file contains source code from another GitHub project. The comments made there apply. The source code
was licensed under the GNU General Public License v3.0. The license text and a detailed reference can be
found in the license subfolder. Many thanks to the author of the code.

For reasons of clarity unneeded parts of the original code were not taken over. The original project can
be found on the https://github.com/kurapan/EAST page.

The remaining code has been left in its original state.
"""

import keras.backend as K
import numpy as np
import tensorflow as tf
from keras import regularizers
from keras.applications.resnet50 import ResNet50
from keras.layers import Conv2D, concatenate, BatchNormalization, Lambda, Input, Activation
from keras.models import Model

RESIZE_FACTOR = 2

def resize_bilinear(x):
    return tf.image.resize_bilinear(x, size=[K.shape(x)[1] * RESIZE_FACTOR, K.shape(x)[2] * RESIZE_FACTOR])

class EAST_model:
    def __init__(self, input_size=512):
        input_image = Input(shape=(None, None, 3), name='input_image')
        overly_small_text_region_training_mask = Input(shape=(None, None, 1),
                                                       name='overly_small_text_region_training_mask')
        text_region_boundary_training_mask = Input(shape=(None, None, 1), name='text_region_boundary_training_mask')
        target_score_map = Input(shape=(None, None, 1), name='target_score_map')
        resnet = ResNet50(input_tensor=input_image, weights='imagenet', include_top=False, pooling=None)
        x = resnet.get_layer('activation_49').output

        x = Lambda(resize_bilinear, name='resize_1')(x)
        x = concatenate([x, resnet.get_layer('activation_40').output], axis=3)
        x = Conv2D(128, (1, 1), padding='same', kernel_regularizer=regularizers.l2(1e-5))(x)
        x = BatchNormalization(momentum=0.997, epsilon=1e-5, scale=True)(x)
        x = Activation('relu')(x)
        x = Conv2D(128, (3, 3), padding='same', kernel_regularizer=regularizers.l2(1e-5))(x)
        x = BatchNormalization(momentum=0.997, epsilon=1e-5, scale=True)(x)
        x = Activation('relu')(x)

        x = Lambda(resize_bilinear, name='resize_2')(x)
        x = concatenate([x, resnet.get_layer('activation_22').output], axis=3)
        x = Conv2D(64, (1, 1), padding='same', kernel_regularizer=regularizers.l2(1e-5))(x)
        x = BatchNormalization(momentum=0.997, epsilon=1e-5, scale=True)(x)
        x = Activation('relu')(x)
        x = Conv2D(64, (3, 3), padding='same', kernel_regularizer=regularizers.l2(1e-5))(x)
        x = BatchNormalization(momentum=0.997, epsilon=1e-5, scale=True)(x)
        x = Activation('relu')(x)

        x = Lambda(resize_bilinear, name='resize_3')(x)
        # x = concatenate([x, ZeroPadding2D(((1, 0),(1, 0)))(resnet.get_layer('activation_10').output)], axis=3)
        x = concatenate([x, resnet.get_layer('activation_10').output], axis=3)
        x = Conv2D(32, (1, 1), padding='same', kernel_regularizer=regularizers.l2(1e-5))(x)
        x = BatchNormalization(momentum=0.997, epsilon=1e-5, scale=True)(x)
        x = Activation('relu')(x)
        x = Conv2D(32, (3, 3), padding='same', kernel_regularizer=regularizers.l2(1e-5))(x)
        x = BatchNormalization(momentum=0.997, epsilon=1e-5, scale=True)(x)
        x = Activation('relu')(x)

        x = Conv2D(32, (3, 3), padding='same', kernel_regularizer=regularizers.l2(1e-5))(x)
        x = BatchNormalization(momentum=0.997, epsilon=1e-5, scale=True)(x)
        x = Activation('relu')(x)

        pred_score_map = Conv2D(1, (1, 1), activation=tf.nn.sigmoid, name='pred_score_map')(x)
        rbox_geo_map = Conv2D(4, (1, 1), activation=tf.nn.sigmoid, name='rbox_geo_map')(x)
        rbox_geo_map = Lambda(lambda x: x * input_size)(rbox_geo_map)
        angle_map = Conv2D(1, (1, 1), activation=tf.nn.sigmoid, name='rbox_angle_map')(x)
        angle_map = Lambda(lambda x: (x - 0.5) * np.pi / 2)(angle_map)
        pred_geo_map = concatenate([rbox_geo_map, angle_map], axis=3, name='pred_geo_map')

        model = Model(inputs=[input_image, overly_small_text_region_training_mask, text_region_boundary_training_mask,
                              target_score_map], outputs=[pred_score_map, pred_geo_map])

        self.model = model
        self.input_image = input_image
        self.overly_small_text_region_training_mask = overly_small_text_region_training_mask
        self.text_region_boundary_training_mask = text_region_boundary_training_mask
        self.target_score_map = target_score_map
        self.pred_score_map = pred_score_map
        self.pred_geo_map = pred_geo_map


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

# returns date as String YYYYMMDDHHMM and rainfall as float32
def get_data(line, dateIndex=1, rainfallIndex=3):
    a = line.split(";")
    return (a[dateIndex], np.float32(a[rainfallIndex]))

# change representation of String from .txt file
def convert_date_pretty(date):
    return "{}.{}.{} {}:{}".format(date[6:8], date[4:6], date[0:4], date[8:10], date[10:12])

# generate some statistics about selected DWD .txt File
def gen_statistic(path, bins=300, show=False):
    assert PATH.__contains__("02712") # check ID is Valid!
    rainsum_month = np.zeros(12, dtype=np.float32)
    rainsum_reference = np.array([100,20,50,35,50,50,40,40,45,28,0,0])
    all_data = []
    sum = 0
    n_rain = 0
    max_value = 0
    date = "NEVER"
    # open file
    fp = open(path, "r")
    headers = fp.readline()
    for line in fp:
        datum, value = get_data(line)
        all_data.append(value)
        sum += value
        if value > 0.0:
            n_rain += 1
            month = int(datum[4:6])-1
            rainsum_month[month] += value
            if value > max_value:
                max_value = value
                date = datum
        if value > 2.50:
            print(convert_date_pretty(datum), value)

    fp.close()
    if show:
        plt.hist(all_data, bins=bins, log=True)
        plt.show()
    np_data = np.array(all_data, dtype=np.float32)
    print("Ausgewertetes File :", path)
    # ToDo: 0.99 Quantil?
    print("maximalwert:                              {:1.2f} erreicht am {}".format(max_value, date))
    print("Regen/gesamt:                             {:1.2f}%".format(int(n_rain * 100 / np_data.size)))
    print("gesamt Regenmenge:                        {:1.2f}".format(np_data.sum()))
    print("Durchschnittliche regenmenge:             {:1.2f}".format(np_data.sum() / np_data.size))
    print("Durchschnittliche Regenmenge (bei Regen): {:1.2f}".format(np_data.sum() / n_rain))
    print("Monat\tReferenz\tmessung\tdiff")
    for i in range(12):
        rel_err = -1
        if rainsum_reference[i]>0:
            rel_err = abs((rainsum_reference[i]-rainsum_month[i])/rainsum_reference[i])
        print("{}:\t{}\t{:1.2f}\t{:1.2f}".format(i, rainsum_reference[i], rainsum_month[i], rel_err))
    return

if __name__ == '__main__':
    # Path to local copy
    PATH = "produkt_ein_min_rr_20180101_20181125_02712"
    gen_statistic(PATH + ".txt")  # resultat zwei werte > 2.50

max_stars_repo_name: thgnaedi/DeepRain

max_issues_repo_name: thgnaedi/DeepRain

max_forks_repo_name: thgnaedi/DeepRain


text:
#!python3

"""
Truthful auction of Heterogeneous Divisible Goods with an approximation algorithm.

References:

    Yonatan Aumann, Yair Dombb, Avinatan Hassidim (2015):
    "Auctioning time: Truthful auctions of heterogeneous divisible goods"
    ACM Transactions on Economics and Computation, 4(1).

Programmers: Naama Berman and Yonatan Lifshitz
Since: 2019-12
"""

from fairpy.agents import *
from fairpy import Allocation

import random, logging
from typing import *
from networkx import *
from math import *

logger = logging.getLogger(__name__)

def stringify_agent_piece_graph(g: Graph):
    """ Convert an agent-piece graph into a string, for display and testing """
    return str([(agent.name(), piece, data) for (agent,piece,data) in g.edges(data=True)])

def stringify_edge_set(s: set):
    """ Convert an agent-piece graph into a string, for display and testing """
    return str(sorted([(agent.name(), piece) for (agent,piece) in s]))

def equally_sized_pieces(agents: List[Agent], piece_size: float) -> Allocation:
    """
    Algorithm 1.
    Approximation algorithm of the optimal auction for uniform-size pieces.

    Complexity and approximation:
    - Requires only 2 / l values from each agent.
    - Runs in time polynomial in n + 1 / l.
    - Approximates the optimal welfare by a factor of 2.

    :param agents: A list of Agent objects.
    :param piece_size: Size of an equally sized piece (in the paper: l).
    :return: A cake-allocation, not necessarily all the cake will be allocated.

    The doctest will work when the set of edges will return according lexicographic order
    >>> Alice = PiecewiseConstantAgent([100, 1], "Alice")
    >>> Bob = PiecewiseConstantAgent([2, 90], "Bob")
    >>> equally_sized_pieces([Alice, Bob], 0.5)
    Alice gets {(0, 1)} with value 100.
    Bob gets {(1, 2)} with value 90.
    <BLANKLINE>

    The doctest will work when the set of edges will return according lexicographic order
    >>> Alice = PiecewiseConstantAgent([1, 1, 1, 1, 1], "Alice")
    >>> Bob = PiecewiseConstantAgent([3, 3, 3, 1, 1], "Bob")
    >>> equally_sized_pieces([Alice, Bob], 3 / 5)
    Alice gets {(2, 5)} with value 3.
    Bob gets {(0, 3)} with value 9.
    <BLANKLINE>
    """
    # > Bob gets {(0, 3)} with value 9.00

    # Initializing variables and asserting conditions
    num_of_agents = len(agents)
    if num_of_agents == 0:
        raise ValueError("There must be at least one agent")
    if not 0 < piece_size <= 1:
        raise ValueError("Piece size must be between 0 and 1")

    logger.info("Piece size (l) = %f", piece_size)
    delta = 1 - int(1 / piece_size) * piece_size
    logger.info("Delta := 1 - floor(1 / l) * l = %f", delta)

    logger.info("Create the partitions P_0_l and P_d_l")
    # Creating the partition of the pieces that start from 0
    partition_0_l = create_partition(piece_size)
    logger.info("  The partition P_0_l (l-sized pieces starting at 0) = %s", partition_0_l)
    # Creating the partition of the pieces that start from delta
    partition_delta_l = create_partition(piece_size, start=delta)
    logger.info("  The partition P_d_l (l-sized pieces starting at delta) = %s", partition_delta_l)
    # Merging the partitions to one partition
    all_partitions = partition_0_l + partition_delta_l

    length = max([a.cake_length() for a in agents])
    # Normalizing the partitions to match the form of the pieces allocation of the Agents
    normalize_partitions = [(int(p[0] * length), int(p[1] * length)) for p in all_partitions]
    normalize_partitions_0_l = [(int(p[0] * length), int(p[1] * length)) for p in partition_0_l]
    normalize_partitions_delta_l = [(int(p[0] * length), int(p[1] * length)) for p in partition_delta_l]

    # Evaluating the pieces of the partition for every agent there is
    logger.info("For each piece (in both partitions) and agent: compute the agent's value of the piece.")
    evaluations = {}
    # Get evaluation for every piece
    for piece in normalize_partitions:
        # For every piece get evaluation for every agent
        for agent in agents:
            evaluations[(agent, piece)] = agent.eval(start=piece[0], end=piece[1])
    # Create the matching graph
    # One side is the agents, the other side is the partitions and the weights are the evaluations
    logger.info("Create the partition graphs G_0_l and G_d_l")
    g_0_l = create_matching_graph(agents, normalize_partitions_0_l, evaluations)
    logger.info("  The graph G_0_l = %s", stringify_agent_piece_graph(g_0_l))
    g_delta_l = create_matching_graph(agents, normalize_partitions_delta_l, evaluations)
    logger.info("  The graph G_d_l = %s", stringify_agent_piece_graph(g_delta_l))

    # Set the edges to be in order, (Agent, partition)
    logger.info("Compute maximum weight matchings for each graph respectively")
    edges_set_0_l = fix_edges(max_weight_matching(g_0_l))
    logger.info("  The edges in G_0_l = %s", stringify_edge_set(edges_set_0_l))
    edges_set_delta_l = fix_edges(max_weight_matching(g_delta_l))
    logger.info("  The edges in G_d_l = %s", stringify_edge_set(edges_set_delta_l))

    logger.info("Choose the heavier among the matchings")
    # Check which matching is heavier and choose it
    if calculate_weight(g_delta_l, edges_set_delta_l) > calculate_weight(g_0_l, edges_set_0_l):
        edges_set = edges_set_delta_l
    else:
        edges_set = edges_set_0_l

    # Find the agents that are in the allocation that was chosen
    chosen_agents = [agent for (agent,piece) in edges_set]
    chosen_agents.sort(key=lambda agent:agent.name())
    # Create allocation

    pieces = len(chosen_agents)*[None]
    # Add the edges to the allocation
    for edge in edges_set:
        pieces[chosen_agents.index(edge[0])] = [edge[1]]
    return Allocation(chosen_agents, pieces)

def discrete_setting(agents: List[Agent], pieces: List[Tuple[float, float]]) -> Allocation:
    """
    Algorithm 2.
    Approximation algorithm of the optimal auction for a discrete cake with known piece sizes.

    Complexity and approximation:
    - Requires at most 2m values from each agent.
    - Runs in time polynomial in n + log m.
    - Approximates the optimal welfare by a factor of log m + 1.

    :param agents: A list of Agent objects.
    :param pieces: List of sized pieces.
    :return: A cake-allocation.

    The doctest will work when the set of edges will return according lexicographic order
    >>> Alice = PiecewiseConstantAgent([100, 1], "Alice")
    >>> Bob = PiecewiseConstantAgent([2, 90], "Bob")
    >>> discrete_setting([Alice, Bob], [(0, 1), (1, 2)])
    Alice gets {(0, 1)} with value 100.
    Bob gets {(1, 2)} with value 90.
    <BLANKLINE>

    """
    # Set m to be the number of pieces in the given partition
    m = len(pieces)
    # Set r to be log of the number of pieces
    r = int(log(m, 2))

    max_weight = 0
    max_match = None

    logger.info("For every t = 0,...,r create the 2 ^ t-partition, partition sequence of 2 ^ t items.")
    logger.info("Denote the t-th partition by Pt.")
    # Go over the partition by powers of 2
    for t in range(0, r + 1):
        logger.info("Iteration t = %d", t)
        # Change the partition to be a partition with 2^t size of every piece
        partition_i = change_partition(pieces, t)

        logger.info("For each piece and agent: compute the agent's value of the piece.")
        # Evaluate every piece in the new partition
        evaluations = {}
        # Go over every piece in the partition
        for piece in partition_i:
            # Go over each Agent
            for agent in agents:
                # Evaluate the piece according to the Agent
                evaluations[(agent, piece)] = agent.eval(start=piece[0], end=piece[1])

        logger.info("create the partition graph G - Pt=%d", t)
        # Create the matching graph according to the new partition
        g_i = create_matching_graph(agents, partition_i, evaluations)
        logger.info("Compute a maximum weight matching Mt in the graph GPt")
        # Find the max weight matching of the graph and get the set of edges of the matching
        edges_set = max_weight_matching(g_i)
        # Set the edges to be in order, (Agent, partition)
        edges_set = fix_edges(edges_set)
        # Calculate the sum of the weights in the edges set
        weight = calculate_weight(g_i, edges_set)
        # Check for the max weight
        if weight > max_weight:
            max_weight = weight
            # Keep the edges set of the max weight
            max_match = edges_set

    # Get the agents that are part of the edges of the max weight
    chosen_agents = [edge[0] for edge in max_match]
    chosen_agents.sort(key=lambda agent:agent.name())
    # Create the allocation
    pieces = len(chosen_agents)*[None]
    # Add the edges to the allocation
    for edge in max_match:
        pieces[chosen_agents.index(edge[0])] = [edge[1]]
    return Allocation(chosen_agents, pieces)

def continuous_setting(agents: List[Agent]) -> Allocation:
    """
    Algorithm 3.
    Approximation algorithm of the optimal auction for a continuous cake.

    Complexity and approximation:
    - Requires at most 2n2 values from each agent.
    - Runs in time polynomial in n.
    - Approximates the optimal welfare by a factor of O(log n).

    :param agents: A list of Agent objects.
    :return: A cake-allocation.

    >>> Alice1 = PiecewiseConstantAgent([100, 1], "Alice")
    >>> Alice2 = PiecewiseConstantAgent([100, 1], "Alice")
    >>> continuous_setting([Alice1, Alice2])
    Alice gets {(0, 2.0)} with value 101.
    <BLANKLINE>
    """
    # set n to be the number of agents
    n = len(agents)
    logger.info("Choose n/2 agents at random. Denote this set by S.")
    # Choose randomly half of the agents
    s = random.choices(agents, k=n//2)
    # Create a new partition
    partitions = set()
    # Add the start to the partition
    partitions.add(0)

    logger.info("For every agent i in S, ask i to divide [0, 1] into 2n intervals of equal worth")
    # Go over all the agents that were chosen
    for a in s:
        start = 0
        # Get pieces with value of 2n
        for i in range(0,2*n):
            end = a.mark(start, a.total_value() / (2 * n))
            # if the piece is out of boundaries we don't add it to the partition
            if end is None:
                break
            end = float("%.4f" % end)
            # Add the piece to the partition
            partitions.add(end)
            start = end

    partitions = list(partitions)
    # Sort the pieces
    partitions = sorted(partitions)

    logger.info("Generate a partition J by taking the union of all boundary points reported by the agents of S.")
    # Turn the list of the pieces into one partition
    start = partitions[0]
    pieces = []
    # Go over all the parts of the partitions and turn it to one partition
    for part in partitions[1:]:
        p = (start, part)
        pieces.append(p)
        start = part

    logger.info("Invoke Algorithm 2 on the rest of the agents and on the sequence of items in J")
    # Get the agents that were nor chosen
    agents = list(set(agents) - set(s))
    # Find the best allocation for those agents with the partition we generated and use Algo 2 to do that
    res = discrete_setting(agents, pieces)
    # Return the allocation
    return res

def create_partition(size: float, start: float=0) -> List[Tuple[float, float]]:
    """
    Used in algorithm 1.
    Creating a partition of [0, 1] with equally sized pieces of given size starting from a given start.
    :param size: The size of the pieces.
    :param start: The location the pieces will start from.
    :return: A partition as described.

    >>> create_partition(0.5, 0)
    [(0, 0.5), (0.5, 1.0)]

    """
    #print("create_partition: size- ", size, "start- ", start)
    res = []
    end = start + size
    # Iterate until we divide all the cake into pieces
    while end <= 1:
        # add the piece to the list
        res.append((start, end))
        start = end
        end = start + size
    #print("create_partition: return value- ", res)
    return res

def fix_edges(edges_set: Set[Tuple[Agent, Tuple[float, float]]]) -> Set[Tuple[Agent, Tuple[float, float]]]:
    """
    Used in algorithm 1 and 2.
    Fix the edge format, sometimes the edges are written backwards
    since the matching algorithm does not care about the edge direction.
    Each edge contains agent and a piece, this function will make sure the agent comes first in the edge.
    :param edges_set: A set of edges to fix.
    :return: A copy of the fixed set of edges.

    >>> Alice = PiecewiseConstantAgent([100, 1], "Alice")
    >>> Bob = PiecewiseConstantAgent([2, 90], "Bob")
    >>> partitions = [(0, 1), (1, 2)]
    >>> edges_set = {(Alice, (0, 1.4889))}
    >>> fix_edges(edges_set)
    [(Alice is an agent with a Piecewise-constant valuation with values [100   1] and total value=101, (0, 1.4889))]
    """
    ret = []
    # Go over all the edges and check if they are in the right order
    for edge in edges_set:
        # If the partition is first we swap the sides of the edge
        if not isinstance(edge[0], Agent):
            ret.append((edge[1], edge[0]))
        else:
            # The Agent is first and we leave it like that
            ret.append((edge[0], edge[1]))
    # we return the set of edges when all the edges are in the right order of (Agent, partition)

    return ret

def change_partition(partition: List[tuple], t: int) -> List[tuple]:
    """
    Used in algorithm 2.
    Create a partition from original partition where each 2 ^ t pieces are united.
    :param partition: The original partition.
    :param t: Defines the size of the new partition.
    :return: A partition with pieces with 2 ^ t size.

    >>> change_partition([(0.0, 1.0), (1.0, 2.0)], 1)
    [(0.0, 2.0)]

    """

    ret = []
    # Go over all the original partitions with 2^t jumps
    for start in range(0, len(partition) - 2 ** t + 1, 2 ** t):
        end = start + 2 ** t - 1
        # Add the new joined partition to the list
        ret.append((partition[start][0], partition[end][1]))
    return ret

def calculate_weight(g: Graph, edges_set: Set[Tuple[Agent, Tuple[float, float]]]) -> float:
    """
    Used in algorithm 2.
    Calculates the weight of a match over a graph.
    :param g: The graph with all the weights.
    :param edges_set: The edges of the matching - for which we will sum the weight.
    :return: A single number - the total weight.

    >>> Alice = PiecewiseConstantAgent([100, 1], "Alice")
    >>> Bob = PiecewiseConstantAgent([2, 90], "Bob")
    >>> partitions = [(0, 1), (1, 2)]
    >>> g = create_matching_graph([Alice, Bob], partitions, {(Alice, (0, 1)): 100.0, (Bob, (0, 1)): 2.0, (Alice, (1, 2)): 1.0, (Bob, (1, 2)): 90.0})
    >>> calculate_weight(g, {(Bob,(1,2)), (Alice, (0,1))})
    190.0

    """
    ret = 0
    # Go over all the weights of the edges and sum the weights
    for edge in edges_set:
        ret += g.get_edge_data(edge[0], edge[1])['weight']
    return ret

def create_matching_graph(left: List[Agent], right: List[Tuple[float, float]],
                          weights: Dict[Tuple[Agent, Tuple[float, float]], float])-> Graph:
    """
    Used in algorithm 2 and 3.
    Creating a weighted bi-partition graph that represents agents, cake pieces and values.
    :param left: List of agents.
    :param right: List of cake pieces.
    :param weights: A dictionary from agents to pieces - represents the value of each agent to each piece.
    :return: A graph object from the given parameters.

    >>> Alice = PiecewiseConstantAgent([100, 1], "Alice")
    >>> Bob = PiecewiseConstantAgent([2, 90], "Bob")
    >>> partitions = [(0, 1), (1, 2)]
    >>> g = create_matching_graph([Alice, Bob], partitions, {(Alice, (0, 1)): 100.0, (Bob, (0, 1)): 2.0, (Alice, (1, 2)): 1.0, (Bob, (1, 2)): 90.0})
    >>> list(g.edges(data=True))
    [(Alice is an agent with a Piecewise-constant valuation with values [100   1] and total value=101, (0, 1), {'weight': 100.0}), (Alice is an agent with a Piecewise-constant valuation with values [100   1] and total value=101, (1, 2), {'weight': 1.0}), (Bob is an agent with a Piecewise-constant valuation with values [ 2 90] and total value=92, (0, 1), {'weight': 2.0}), (Bob is an agent with a Piecewise-constant valuation with values [ 2 90] and total value=92, (1, 2), {'weight': 90.0})]

    """

    # Create the graph
    g = Graph()
    # g.edges(data = True)
    # Set the left side of the graph to be the Agents
    g.add_nodes_from(left, bipartite=0)
    # Set the right side of the graph to be the partitions
    g.add_nodes_from(right, bipartite=1)
    # Set the edges of the graph with their weights
    for key, value in weights.items():
        g.add_edge(key[0], key[1], weight=value)
    return g

if __name__ == "__main__":
    import doctest
    (failures,tests) = doctest.testmod(report=True)
    print ("{} failures, {} tests".format(failures,tests))

max_stars_repo_name: ShmuelLa/research-algorithms

max_issues_repo_name: ShmuelLa/research-algorithms

max_forks_repo_name: ShmuelLa/research-algorithms

reason: from networkx


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Cesare De Cal
Data Fitting
Exercise 3

"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
from scipy.interpolate import lagrange
from numpy.polynomial.polynomial import Polynomial

# Points provided by the chemistry experiment
# We store the x and y values in numpy arrays
x = np.array([-1, -0.96, -0.86, -0.79, 0.22, 0.50, 0.93])
y = np.array([-1.000, -0.151, 0.894, 0.986, 0.895, 0.500, -0.306])

# Create array of 1000 equi-distant points
xnew = np.linspace(min(x), max(x), 1000)

# Increase plot size
plt.figure(dpi=130)

# Set axis labels
plt.xlabel('x')
plt.ylabel('f(x)')

# Let's first plot the points
plt.plot(x, y, 'bo', label='Data')

# Compute Interpolating Polynomial with Lagrange method
poly = lagrange(x, y)
print('Lagrange Polynomial:\n', poly)

plt.plot(xnew, poly(xnew), '-', label='Interpolating Polynomial')

# Get max precision possibile
np.set_printoptions(precision=100)

# Print coefficients
print('Coefficients of interpolating polynomial:', Polynomial(poly).coef)

# Compute Natural Cubic Spline
ncspline = CubicSpline(x, y, bc_type='natural')
print('Coefficients of cubic spline:', ncspline.c)

# Make y axis a  little taller to show legend
plt.ylim(-1.3, 1.7)
plt.plot(xnew, ncspline(xnew), 'r', label='Natural Cubic Spline')

# Final tweaks to plot
plt.legend(loc=1, prop={'size': 8})
plt.show()

max_stars_repo_path: Data Fitting/CesareDeCal/exercise3.py

max_stars_repo_name: csr/MATLAB-Scientific-Programming

max_issues_repo_path: Data Fitting/CesareDeCal/exercise3.py

max_issues_repo_name: csr/MATLAB-Scientific-Programming

max_forks_repo_path: Data Fitting/CesareDeCal/exercise3.py

max_forks_repo_name: csr/MATLAB-Scientific-Programming

reason: import numpy,from numpy,from scipy


text:
from graphene import ObjectType, Int, List, Field, Argument
from graphql import GraphQLError

from app.services import numpy_service

import numpy as np

# This variable is defined here so that it's not exposed as a field on the GraphQL type.
matrix = List(List(Int))

class MatrixOperations(ObjectType):
	# Return type followed by args
	product = Field(
		matrix, first=Argument(matrix, required=True), second=Argument(matrix, required=True))

	@staticmethod
	def resolve_product(
			parent, info, first: type(matrix), second: type(matrix)) -> type(matrix):
		"""
		Somehow Graphene.List works as an argument to np.asarray and we can also automatically
		convert from a python list back to a Graphene.List. It also seems like exception messages
		are automatically passed to the GraphQL client in the response.
		"""
		first_matrix = np.asarray(first, dtype=np.int32)
		second_matrix = np.asarray(second, dtype=np.int32)

		try:
			product = numpy_service.matrix_product(first_matrix, second_matrix)
		except ValueError:
			raise GraphQLError('Invalid arguments')

		return product.tolist()

max_stars_repo_name: zengspr/DiSPy-API

max_issues_repo_name: zengspr/DiSPy-API

max_forks_repo_name: zengspr/DiSPy-API


text:
# usage: file

import sys
import os
import seaborn as sns
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import numpy as np
from keyname import keyname as kn
from fileshash import fileshash as fsh
import h5py
import colorsys
from tqdm import tqdm
from joblib import delayed, Parallel
from natsort import natsorted

matplotlib.rcParams['pdf.fonttype'] = 42

filename = sys.argv[1]

def RenderAndSave(upd_key, filename):

    file = h5py.File(filename, 'r')
    nlev = int(file.attrs.get('NLEV'))

    live = np.array(file['Live'][upd_key])
    traffic = np.array(file['InboxTraffic']['dir_4'][upd_key])

    image = np.array([
        [
            # dead
            (0.0, 0.0, 0.0)
            if not val_live else
            (lambda n: (1.0, 1.0, 0.0) if n > 4 else {
                0 : (1.0, 1.0, 1.0), # white
                1 : (0.0, 1.0, 0.0), # green
                2 : (0.0, 0.0, 1.0), # blue
                3 : (1.0, 0.0, 1.0), # purple
                4 : (1.0, 0.0, 0.0), # red
            }[n])(val_traffic)
            for val_traffic, val_live in zip(row_traffic, row_live)
        ]
    for row_traffic, row_live in zip(traffic, live)])

    plt.figure(figsize=(18,18))

    plt.imshow(
        image,
        extent = (0, image.shape[1], image.shape[0], 0)
    )

    plt.axis('off')
    plt.grid(b=None)

    lines = LineCollection([
        ((x,y), dest)
        for x in range(image.shape[0])
        for y in range(image.shape[1])
        for dest in ((x+1,y), (x,y+1))
    ], linestyle='solid', colors='black')
    plt.gca().add_collection(lines)

    plt.savefig(
        kn.pack({
            'title' : 'messaging_spiker_viz',
            'update' : ''.join(c for c in upd_key if c.isdigit()),
            'seed' : kn.unpack(filename)['seed'],
            'treat' : kn.unpack(filename)['treat'],
            '_data_hathash_hash' : fsh.FilesHash().hash_files([filename]),
            '_script_fullcat_hash' : fsh.FilesHash(
                                                file_parcel="full_parcel",
                                                files_join="cat_join"
                                            ).hash_files([sys.argv[0]]),
            '_source_hash' :kn.unpack(filename)['_source_hash'],
            'ext' : '.png'
        }),
        transparent=True,
        bbox_inches='tight',
        pad_inches=0
    )

    plt.clf()
    plt.close(plt.gcf())

file = h5py.File(filename, 'r')
upd_keys = natsorted(
    [key for key in file['InboxTraffic']['dir_4']]
)[-16:]

Parallel(n_jobs=-1)(
    delayed(RenderAndSave)(upd_key, filename) for upd_key in tqdm(upd_keys)
)

max_stars_repo_name: schregardusc/dishtiny

max_issues_repo_name: schregardusc/dishtiny

max_forks_repo_name: schregardusc/dishtiny


text:
from matplotlib import pyplot
import matplotlib.pyplot as plt
import numpy

#import Pearson

#templateSpectre = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 42, 40, 26, 26, 24, 24, 24, 26, 26, 10, 12, 10, 8, 6, 6, 6, 6, 6, 6, 8, 6, 8, 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]
#testSpectre =     [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6, 42, 40, 26, 26, 24, 24, 24, 26, 26, 12, 10, 10, 8, 6, 6, 6, 6, 8, 6, 6, 6, 8, 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]

# def remove_outLiers(points):
# 	model = LinearLeastSquaresModel()
# 	minimumNumber = int(len(points)*0.3)
# 	iterations = 1
# 	Y = [p.y for p in points]
# 	mY = sum(Y)/len(Y)
# 	fitDist = mY*0.5
# 	numberOfcloses = int(len(points)*0.5)
# 	inliers = ransac(points, model, minimumNumber, 1, fitDist, numberOfcloses, debug=False,return_all=True, random=False)

# 	return inliers

def leastSquares(points):
	x_d = []
	y_d = []
	for point in points:
		(x, y) = point

		x_d.append(x)
		y_d.append(y)

	n=len(x_d)

	B=numpy.array(y_d)
	A=numpy.array(([[x_d[j], 1] for j in range(n)]))
	X=numpy.linalg.lstsq(A,B)[0]
	a=X[0]; b=X[1]
	#print("Line is: y=",a,"x+",b)

	Xplot = []
	Yplot = []

	x = min(x_d)
	y = a*x + b
	Xplot.append(x)
	Yplot.append(y)

	x = max(x_d)
	y = a*x + b
	Xplot.append(x)
	Yplot.append(y)
	return (Xplot, Yplot, a, b)

fig = plt.figure()

ax1 = fig.add_subplot(1,1,1)
#ax1.set_title('before adjust', fontname='serif', fontsize=16)

templateSpectre = [2, 3, 4, 5, 6, 7, 0, 1]
print("len(templateSpectre) = ", len(templateSpectre))
#print("len(testSpectre) = ", len(testSpectre))

#UFPE
correctPairsIndices =  [(9, 0), (10, 1), (14, 5), (15, 7), (16, 8), (12, 2), (17, 9), (20, 13), (21, 14), (22, 15), (23, 16), (0, 17), (2, 19), (4, 20), (5, 22), (8, 23)]
#correctPairsIndices = [(4, 23), (5, 22), (6, 20), (7, 17), (8, 15), (0, 9), (1, 8), (2, 1), (3, 0)]
#correctPairsIndices = [(14, 0), (21, 3), (25, 6), (33, 5), (2, 2), (4, 9), (8, 1), (11, 7), (38, 8)]

#correctPairsIndices = [(11, 15), (14, 11), (23, 5), (24, 3), (35, 16), (16, 17), (25, 6), (37, 14), (5, 13), (7, 12), (9, 0)]

correctPairsIndices =  [(11, 15), (14, 11), (23, 5), (24, 3), (35, 16), (16, 17), (25, 6), (37, 14), (5, 13), (7, 12), (9, 0)]
Y = [15, 11, 17, 5, 3, 6, 16, 14, 13, 12, 0]

#[(0, 17), (2, 19), (4, 20), (5, 22), (8, 23), (9, 0), (11, 1), (14, 5), (15, 7), (16, 8), (12, 2), (17, 9), (20, 13), (21, 14), (22, 15), (23, 16)]

# tau:  [0, 1, 2, 3, 4, 5, 6, 7, 8]
# Q:   [23, 22, 20, 17, 15, 9, 8, 1, 0]

# UFPE
#Y = [0, 1, 2, 5, 7, 8, 9, 13, 14, 15, 16, 17, 19, 20, 22, 23]#[y for (x,y) in correctPairsIndices]
#Y = [23, 22, 20, 17, 15, 9, 8, 1, 0]
#Y = [0, 3, 6, 5, 8, 2, 9, 1, 7]

#tau:  [5, 7, 11, 14, 25, 37]
#Y =  [15, 11, 6, 14, 13, 12]

X = [x for (x, y) in correctPairsIndices if y in Y]

X.sort()

new_pairs = list(zip(X, Y))
print("new_pairs: ", new_pairs)

T = X#[0, 1, 2, 3, 4, 5, 6, 7]
Q = Y#[1, 3, 4, 5, 6, 7, 1, 8]

points = list(zip(T, Q))

(Xplot, Yplot, a, b) = leastSquares(points)

#plt.plot(T, Q, marker='o', linestyle=':', color='b')
plt.plot(T, Q, marker='o', markersize='5.0', linestyle='None', color='b')
plt.plot(Xplot, Yplot, linestyle='-', linewidth=2.0, color='c')
plt.plot(T, Q, linestyle='--', linewidth=1.5, color='m')
#ax1.text(0, 6.5, r'Template vertex: $T_i$', fontsize=11)
#ax1.text(0, 6.2, r'Query vertex: $Q_j$', fontsize=11)

#for (t, q) in zip(T, Q):
	#text = r"$(T_" + str(t) + ", Q_" + str(q) + ")$"
	#plt.text(t-0.5, q+0.5, r"$(\tau_{" + str(t) + "}, Q_{" + str(q) + "})$", weight='semibold', fontname='Courier', style='italic', fontsize=21, color='black')

# for (t, q) in zip(T, Q):
# 	#text = "$(T_" + str(t) + ", Q_" + str(q) + ")$"
# 	plt.text(t+0.2, q, r"$(T_" + str(t) + ", Q_" + str(q) + ")$", fontname='serif', style='italic', fontsize=14, color='black')

# T = [6, 7]
# Q = [0, 1]
# plt.plot(T, Q, marker='x', linestyle=':', color='r')

plt.show()