text:
"""
    Name: Amit Kachroo
    ----- class for AIC test
    email : amitkac@gmail.com
"""
import pandas as pd
from scipy.stats import exponweib
import numpy as np
from scipy.stats import (gamma, rayleigh, norm,  nakagami, expon,
                         lognorm)

class aic2():

    def __init__(self, data):
        self.data = data
        # self.fitDist()

    def fitDist(self):

        n = len(self.data)

        # gamma distribution
        gammaParam = gamma.fit(self.data)
        gammaNumPar = len(gammaParam)
        gammaSum = -1*np.sum(np.log(gamma.pdf(self.data, *gammaParam)))
        aicGamma = 2*gammaNumPar+2*gammaSum + (2*gammaNumPar*(gammaNumPar+1)
                                               / (n-gammaNumPar-1))

        # rayleigh distribution
        rayleighParam = rayleigh.fit(self.data)
        rayleighNumPar = len(rayleighParam)
        rayleighSum = -1*np.sum(np.log(rayleigh.pdf(self.data,
                                                    *rayleighParam)))
        aicRayleigh = 2*rayleighNumPar+2*rayleighSum+(2*rayleighNumPar
                                                      * (rayleighNumPar+1) /
                                                      (n-rayleighNumPar-1))

        # normal distribution
        normParam = norm.fit(self.data)
        normNumPar = len(normParam)
        normSum = -1*np.sum(np.log(norm.pdf(self.data, *normParam)))
        aicNorm = 2*normNumPar+2*normSum+(2*normNumPar*(normNumPar+1)
                                          / (n-normNumPar-1))

        # LogNormal distribution
        logNormParam = lognorm.fit(self.data)
        logNormNumPar = len(logNormParam)
        logNormSum = -1*np.sum(np.log(lognorm.pdf(self.data, *logNormParam)))
        aicLogNorm = 2*logNormNumPar+2*logNormSum+(2*logNormNumPar *
                                                   (logNormNumPar+1)
                                                   / (n-logNormNumPar-1))

        # Nakagami distribution
        nakagamiParam = nakagami.fit(self.data)
        nakagamiNumPar = len(nakagamiParam)
        nakagamiSum = -1*np.sum(np.log(nakagami.pdf(self.data,
                                                    *nakagamiParam)))
        aicNakagami = 2*nakagamiNumPar+2*nakagamiSum+(2*nakagamiNumPar *
                                                      (nakagamiNumPar+1)
                                                      / (n-nakagamiNumPar-1))

        # exponential distribution
        exponParam = expon.fit(self.data)
        exponNumPar = len(exponParam)
        exponSum = -1*np.sum(np.log(expon.pdf(self.data, *exponParam)))
        aicExpon = 2*exponNumPar+2*exponSum + (2*exponNumPar *
                                               (exponNumPar+1)
                                               / (n-exponNumPar-1))

        # weibul distribution
        exponweibParam = exponweib.fit(self.data)
        exponweibNumPar = len(exponweibParam)
        exponweibSum = -1*np.sum(np.log(exponweib.pdf(self.data,
                                                      *exponweibParam)))
        aicExpWeib = 2*exponweibNumPar+2*exponweibSum+(2*exponweibNumPar *
                                                       (exponweibNumPar+1)
                                                       / (n-exponweibNumPar-1))

        return (aicGamma, aicRayleigh, aicNorm, aicLogNorm, aicNakagami,
                aicExpon, aicExpWeib)

    @classmethod
    def returnDistData(cls, self):
        gammaParam = gamma.fit(10**(self.data/10))
        gammaDist = gamma.pdf(self.data, *gammaParam)

        rayleighParam = rayleigh.fit(self.data)
        rayleighDist = rayleigh.pdf(self.data, *rayleighParam)

        normParam = norm.fit(self.data)
        normDist = norm.pdf(self.data, *normParam)

        logNormParam = lognorm.fit(self.data)
        lognormDist = lognorm.pdf(self.data, *logNormParam)

        nakagamiParam = nakagami.fit(self.data)
        nakagamiDist = nakagami.pdf(self.data, *nakagamiParam)

        exponParam = expon.fit(self.data)
        exponDist = expon.pdf(self.data, *exponParam)

        exponweibParam = exponweib.fit(self.data)
        weibDist = exponweib.pdf(self.data, *exponweibParam)

        distDF = pd.DataFrame(np.column_stack([gammaDist, rayleighDist,
                                               normDist, lognormDist,
                                               nakagamiDist, exponDist,
                                               weibDist]),
                              columns=['gammaDist', 'rayleighDist',
                                       'normDist', 'lognormDist',
                                       'nakagamiDist', 'exponDist',
                                       'weibDist'])
        self.distDF = distDF

    def getDistGamma(self):
        gammaParam = gamma.fit(self.data)
        gammaDist = gamma.pdf(self.data, *gammaParam)
        return gammaDist

    def getDistNorm(self):
        normParam = norm.fit(self.data)
        normDist = norm.pdf(self.data, *normParam)
        return normDist

        # , aicExpWeib)

max_stars_repo_name: amitkac/StatisticalModelling

max_issues_repo_name: amitkac/StatisticalModelling

max_forks_repo_name: amitkac/StatisticalModelling


text:
import unittest
import os
from tempfile import TemporaryDirectory
from typing import cast

import numpy as np
from allennlp.common.params import Params
from allennlp.data.vocabulary import Vocabulary
from allennlp.data.dataset_readers.dataset_reader import DatasetReader

from rulm.settings import REMEMBERING_EXAMPLE, ENCODER_ONLY_MODEL_PARAMS, \
    ENCODER_ONLY_SAMPLED_SOFTMAX_MODEL_PARAMS, DEFAULT_VOCAB_DIR, TRAIN_EXAMPLE, \
    TRAIN_VOCAB_EXAMPLE, TEST_EXAMPLE
from rulm.language_model import LanguageModel
from rulm.models.neural_net import NeuralNetLanguageModel

class TestRNNLM(unittest.TestCase):
    @classmethod
    def setUpClass(cls):
        LanguageModel.set_seed(42)
        configs = (ENCODER_ONLY_MODEL_PARAMS, ENCODER_ONLY_SAMPLED_SOFTMAX_MODEL_PARAMS)
        cls.params_sets = [Params.from_file(config) for config in configs]

        cls.vocabularies = []
        for params in cls.params_sets:
            vocabulary_params = params.pop("vocabulary", default=Params({}))
            reader_params = params.duplicate().pop("reader", default=Params({}))
            reader = DatasetReader.from_params(reader_params)
            dataset = reader.read(REMEMBERING_EXAMPLE)
            cls.vocabularies.append(Vocabulary.from_params(vocabulary_params, instances=dataset))

        cls.train_vocabulary = Vocabulary.from_files(TRAIN_VOCAB_EXAMPLE)

        cls.sentences = []
        with open(REMEMBERING_EXAMPLE, "r", encoding="utf-8") as r:
            for line in r:
                cls.sentences.append(line.strip())

    def _test_model_predictions(self, model, reverse=False):
        for sentence in self.sentences:
            indices = model.text_to_indices(sentence)[1:-1]
            if reverse:
                indices = indices[::-1]
            sentence = " ".join([model.vocab.get_token_from_index(i) for i in indices])
            first_token = model.vocab.get_token_from_index(indices[0])
            for context_right_border in range(1, len(indices)-1):
                if context_right_border == 1 and (first_token in "Я!.?"):
                    continue
                context = indices[:context_right_border]
                text = " ".join([model.vocab.get_token_from_index(i) for i in context])
                prediction = model.sample_decoding(text, k=1)
                self.assertEqual(prediction, sentence)

    @staticmethod
    def _test_model_equality(model1: NeuralNetLanguageModel, model2: NeuralNetLanguageModel):
        old_params = list([param.detach().cpu().numpy() for param in model1.model.parameters()])
        new_params = list([param.detach().cpu().numpy() for param in model2.model.parameters()])
        for o, n in zip(old_params, new_params):
            np.testing.assert_array_almost_equal(o, n)

    def test_model_from_file(self):
        for params, vocabulary in zip(self.params_sets, self.vocabularies):
            params = params.duplicate()
            train_params = params.pop("train")
            model = LanguageModel.from_params(params, vocab=vocabulary)
            model.train(REMEMBERING_EXAMPLE, train_params)
            self._test_model_predictions(model)

    def test_reversed_model(self):
        for params, vocabulary in zip(self.params_sets, self.vocabularies):
            params = params.duplicate()
            train_params = params.pop('train')
            params["reader"]["reverse"] = True
            model_reversed = LanguageModel.from_params(params, vocab=vocabulary)
            model_reversed.train(REMEMBERING_EXAMPLE, train_params)
            self._test_model_predictions(model_reversed, reverse=True)

    def test_save_load(self):
        for params, vocabulary in zip(self.params_sets, self.vocabularies):
            with TemporaryDirectory() as dirpath:
                vocab_dir = os.path.join(dirpath, DEFAULT_VOCAB_DIR)
                os.mkdir(vocab_dir)
                vocabulary.save_to_files(vocab_dir)

                params = params.duplicate()
                train_params = params.pop('train')
                model = LanguageModel.from_params(params, vocab=vocabulary)
                model.train(REMEMBERING_EXAMPLE, train_params, dirpath)

                loaded_model = LanguageModel.load(dirpath,
                                                  params_file=ENCODER_ONLY_MODEL_PARAMS,
                                                  vocabulary_dir=vocab_dir)

                self.assertTrue(isinstance(model, NeuralNetLanguageModel))
                self.assertTrue(isinstance(loaded_model, NeuralNetLanguageModel))
                self._test_model_equality(cast(NeuralNetLanguageModel, model),
                                          cast(NeuralNetLanguageModel, loaded_model))

    def test_valid_ppl(self):
        for params in self.params_sets:
            params = params.duplicate()
            train_params = params.pop('train')
            train_params["trainer"]["num_epochs"] = 1
            train_params["iterator"]["batch_size"] = 50
            model = LanguageModel.from_params(params, vocab=self.train_vocabulary)
            metrics = model.train(TEST_EXAMPLE, train_params, valid_file_name=TEST_EXAMPLE)
            val_loss = metrics["validation_loss"]
            ppl_state = model.measure_perplexity(TEST_EXAMPLE)
            self.assertAlmostEqual(np.log(ppl_state.avg_perplexity), val_loss, places=3)

ext: py

max_stars_repo_name: IlyaGusev/rulm

max_issues_repo_name: IlyaGusev/rulm

max_forks_repo_name: IlyaGusev/rulm

reason: import numpy


text:
# -*- coding: utf8 -*-

from __future__ import division, print_function
import sys

from itertools import permutations
import numpy as np

from ..ext.ordereddict import OrderedDict

try:
    # Try to use the C++ triangle library, faster than the
    # pure Python version.
    # The latest stable release only works with Python 2. The GitHub version
    # works on Python 3 though, but the release has yet to be done.
    import triangle
    assert sys.version_info.major == 2
    _TRIANGLE_AVAILABLE = True
except (ImportError, AssertionError):
    _TRIANGLE_AVAILABLE = False

class Triangulation(object):
    """Constrained delaunay triangulation

    Implementation based on:

        * Domiter, V. and Žalik, B. Sweep‐line algorithm for constrained
          Delaunay triangulation

    Parameters
    ----------
    pts : array
        Nx2 array of points.
    edges : array
        Nx2 array of edges (dtype=int).

    Notes
    -----
    * Delaunay legalization is not yet implemented. This produces a proper
      triangulation, but adding legalisation would produce fewer thin
      triangles.
    * The pts and edges arrays may be modified.
    """
    def __init__(self, pts, edges):
        self.pts = pts[:, :2].astype(np.float32)
        self.edges = edges
        if self.pts.ndim != 2 or self.pts.shape[1] != 2:
            raise TypeError('pts argument must be ndarray of shape (N, 2).')
        if self.edges.ndim != 2 or self.edges.shape[1] != 2:
            raise TypeError('edges argument must be ndarray of shape (N, 2).')

        # described in initialize()
        self.front = None
        self.tris = OrderedDict()
        self.edges_lookup = {}

    def normalize(self):
        # Clean up data   (not discussed in original publication)

        # (i) Split intersecting edges. Every edge that intersects another
        #     edge or point is split. This extends self.pts and self.edges.
        self.split_intersecting_edges()

        # (ii) Merge identical points. If any two points are found to be equal,
        #      the second is removed and the edge table is updated accordingly.
        self.merge_duplicate_points()

        # (iii) Remove duplicate edges
        # TODO

    def initialize(self):
        self.normalize()
        ## Initialization (sec. 3.3)

        # sort points by y, then x
        flat_shape = self.pts.shape[0] * self.pts.shape[1]
        pts = self.pts.reshape(flat_shape).view([('x', np.float32),
                                                 ('y', np.float32)])
        order = np.argsort(pts, order=('y', 'x'))
        pts = pts[order]
        # update edges to match new point order
        invorder = np.argsort(order)
        self.edges = invorder[self.edges]
        self.pts = pts.view(np.float32).reshape(len(pts), 2)

        # make artificial points P-1 and P-2
        xmax = self.pts[:, 0].max()
        xmin = self.pts[:, 0].min()
        ymax = self.pts[:, 1].max()
        ymin = self.pts[:, 1].min()
        xa = (xmax-xmin) * 0.3
        ya = (ymax-ymin) * 0.3
        p1 = (xmin - xa, ymin - ya)
        p2 = (xmax + xa, ymin - ya)

        # prepend artificial points to point list
        newpts = np.empty((self.pts.shape[0]+2, 2), dtype=float)
        newpts[0] = p1
        newpts[1] = p2
        newpts[2:] = self.pts
        self.pts = newpts
        self.edges += 2

        # find topmost point in each edge
        self.tops = self.edges.max(axis=1)
        self.bottoms = self.edges.min(axis=1)

        # inintialize sweep front
        # values in this list are indexes into self.pts
        self.front = [0, 2, 1]

        # empty triangle list.
        # This will contain [(a, b, c), ...] where a,b,c are indexes into
        # self.pts
        self.tris = OrderedDict()

        # For each triangle, maps (a, b): c
        # This is used to look up the thrid point in a triangle, given any
        # edge. Since each edge has two triangles, they are independently
        # stored as (a, b): c and (b, a): d
        self.edges_lookup = {}

    def triangulate(self):
        self.initialize()

        pts = self.pts
        front = self.front

        ## Begin sweep (sec. 3.4)
        for i in range(3, pts.shape[0]):
            pi = pts[i]
            #debug("========== New point %d: %s ==========" % (i, pi))

            # First, triangulate from front to new point
            # This applies to both "point events" (3.4.1)
            # and "edge events" (3.4.2).

            # get index along front that intersects pts[i]
            l = 0
            while pts[front[l+1], 0] <= pi[0]:
                l += 1
            pl = pts[front[l]]

            # "(i) middle case"
            if pi[0] > pl[0]:
                #debug("  mid case")
                # Add a single triangle connecting pi,pl,pr
                self.add_tri(front[l], front[l+1], i)
                front.insert(l+1, i)
            # "(ii) left case"
            else:
                #debug("  left case")
                # Add triangles connecting pi,pl,ps and pi,pl,pr
                self.add_tri(front[l], front[l+1], i)
                self.add_tri(front[l-1], front[l], i)
                front[l] = i

            #debug(front)

            # Continue adding triangles to smooth out front
            # (heuristics shown in figs. 9, 10)
            #debug("Smoothing front...")
            for direction in -1, 1:
                while True:
                    # Find point connected to pi
                    ind0 = front.index(i)
                    ind1 = ind0 + direction
                    ind2 = ind1 + direction
                    if ind2 < 0 or ind2 >= len(front):
                        break

                    # measure angle made with front
                    p1 = pts[front[ind1]]
                    p2 = pts[front[ind2]]
                    err = np.geterr()
                    np.seterr(invalid='ignore')
                    try:
                        angle = np.arccos(self.cosine(pi, p1, p2))
                    finally:
                        np.seterr(**err)

                    # if angle is < pi/2, make new triangle
                    #debug("Smooth angle:", pi, p1, p2, angle)
                    if angle > np.pi/2. or np.isnan(angle):
                        break

                    assert (i != front[ind1] and
                            front[ind1] != front[ind2] and
                            front[ind2] != i)
                    self.add_tri(i, front[ind1], front[ind2], source='smooth1')
                    front.pop(ind1)
            #debug("Finished smoothing front.")

            # "edge event" (sec. 3.4.2)
            # remove any triangles cut by completed edges and re-fill
            # the holes.
            if i in self.tops:
                for j in self.bottoms[self.tops == i]:
                    # Make sure edge (j, i) is present in mesh
                    # because edge event may have created a new front list
                    self.edge_event(i, j)
                    front = self.front

        self.finalize()

        self.tris = np.array(list(self.tris.keys()), dtype=int)

        #debug("Finished with %d tris:" % self.tris.shape[0])
        #debug(str(self.tris))

    def finalize(self):
        ## Finalize (sec. 3.5)

        # (i) Add bordering triangles to fill hull
        #debug("== Fill hull")
        front = list(OrderedDict.fromkeys(self.front))

        l = len(front) - 2
        k = 1
        while k < l-1:
            # if edges lie in counterclockwise direction, then signed area
            # is positive
            if self.iscounterclockwise(front[k], front[k+1], front[k+2]):
                self.add_tri(front[k], front[k+1], front[k+2], legal=False,
                             source='fill_hull')
                front.pop(k+1)
                l -= 1
                continue
            k += 1

        # (ii) Remove all triangles not inside the hull
        #      (not described in article)
        #debug("== Remove triangles outside hull")

        tris = []  # triangles to check
        tri_state = {}  # 0 for outside, 1 for inside

        # find a starting triangle
        for t in self.tris:
            if 0 in t or 1 in t:
                tri_state[t] = 0
                tris.append(t)
                break

        while tris:
            #debug("iterate:", tris)
            next_tris = []
            for t in tris:
                v = tri_state[t]
                for i in (0, 1, 2):
                    edge = (t[i], t[(i + 1) % 3])
                    pt = t[(i + 2) % 3]
                    t2 = self.adjacent_tri(edge, pt)
                    if t2 is None:
                        continue
                    t2a = t2[1:3] + t2[0:1]
                    t2b = t2[2:3] + t2[0:2]
                    if t2 in tri_state or t2a in tri_state or t2b in tri_state:
                        continue
                    if self.is_constraining_edge(edge):
                        tri_state[t2] = 1 - v
                    else:
                        tri_state[t2] = v
                    next_tris.append(t2)
            tris = next_tris

        for t, v in tri_state.items():
            if v == 0:
                self.remove_tri(*t)

    def edge_event(self, i, j):
        """
        Force edge (i, j) to be present in mesh.
        This works by removing intersected triangles and filling holes up to
        the cutting edge.
        """
        front_index = self.front.index(i)

        #debug("  == edge event ==")
        front = self.front

        # First just see whether this edge is already present
        # (this is not in the published algorithm)
        if (i, j) in self.edges_lookup or (j, i) in self.edges_lookup:
            #debug("    already added.")
            return
        #debug("    Edge (%d,%d) not added yet. Do edge event. (%s - %s)" %
        #      (i, j, pts[i], pts[j]))

        # traverse in two different modes:
        #  1. If cutting edge is below front, traverse through triangles. These
        #     must be removed and the resulting hole re-filled. (fig. 12)
        #  2. If cutting edge is above the front, then follow the front until
        #     crossing under again. (fig. 13)
        # We must be able to switch back and forth between these
        # modes (fig. 14)

        # Collect points that draw the open polygons on either side of the
        # cutting edge. Note that our use of 'upper' and 'lower' is not strict;
        # in some cases the two may be swapped.
        upper_polygon = [i]
        lower_polygon = [i]

        # Keep track of which section of the front must be replaced
        # and with what it should be replaced
        front_holes = []  # contains indexes for sections of front to remove

        next_tri = None   # next triangle to cut (already set if in mode 1)
        last_edge = None  # or last triangle edge crossed (if in mode 1)

        # Which direction to traverse front
        front_dir = 1 if self.pts[j][0] > self.pts[i][0] else -1

        # Initialize search state
        if self.edge_below_front((i, j), front_index):
            mode = 1  # follow triangles
            tri = self.find_cut_triangle((i, j))
            last_edge = self.edge_opposite_point(tri, i)
            next_tri = self.adjacent_tri(last_edge, i)
            assert next_tri is not None
            self.remove_tri(*tri)
            # todo: does this work? can we count on last_edge to be clockwise
            # around point i?
            lower_polygon.append(last_edge[1])
            upper_polygon.append(last_edge[0])
        else:
            mode = 2  # follow front

        # Loop until we reach point j
        while True:
            #debug("  == edge_event loop: mode %d ==" % mode)
            #debug("      front_holes:", front_holes, front)
            #debug("      front_index:", front_index)
            #debug("      next_tri:", next_tri)
            #debug("      last_edge:", last_edge)
            #debug("      upper_polygon:", upper_polygon)
            #debug("      lower_polygon:", lower_polygon)
            #debug("      =====")
            if mode == 1:
                # crossing from one triangle into another
                if j in next_tri:
                    #debug("    -> hit endpoint!")
                    # reached endpoint!
                    # update front / polygons
                    upper_polygon.append(j)
                    lower_polygon.append(j)
                    #debug("    Appended to upper_polygon:", upper_polygon)
                    #debug("    Appended to lower_polygon:", lower_polygon)
                    self.remove_tri(*next_tri)
                    break
                else:
                    # next triangle does not contain the end point; we will
                    # cut one of the two far edges.
                    tri_edges = self.edges_in_tri_except(next_tri, last_edge)

                    # select the edge that is cut
                    last_edge = self.intersected_edge(tri_edges, (i, j))
                    #debug("    set last_edge to intersected edge:", last_edge)
                    last_tri = next_tri
                    next_tri = self.adjacent_tri(last_edge, last_tri)
                    #debug("    set next_tri:", next_tri)
                    self.remove_tri(*last_tri)

                    # Crossing an edge adds one point to one of the polygons
                    if lower_polygon[-1] == last_edge[0]:
                        upper_polygon.append(last_edge[1])
                        #debug("    Appended to upper_polygon:", upper_polygon)
                    elif lower_polygon[-1] == last_edge[1]:
                        upper_polygon.append(last_edge[0])
                        #debug("    Appended to upper_polygon:", upper_polygon)
                    elif upper_polygon[-1] == last_edge[0]:
                        lower_polygon.append(last_edge[1])
                        #debug("    Appended to lower_polygon:", lower_polygon)
                    elif upper_polygon[-1] == last_edge[1]:
                        lower_polygon.append(last_edge[0])
                        #debug("    Appended to lower_polygon:", lower_polygon)
                    else:
                        raise RuntimeError("Something went wrong..")

                    # If we crossed the front, go to mode 2
                    x = self.edge_in_front(last_edge)
                    if x >= 0:  # crossing over front
                        #debug("    -> crossed over front, prepare for mode 2")
                        mode = 2
                        next_tri = None
                        #debug("    set next_tri: None")

                        # where did we cross the front?
                        # nearest to new point
                        front_index = x + (1 if front_dir == -1 else 0)
                        #debug("    set front_index:", front_index)

                        # Select the correct polygon to be lower_polygon
                        # (because mode 2 requires this).
                        # We know that last_edge is in the front, and
                        # front[front_index] is the point _above_ the front.
                        # So if this point is currently the last element in
                        # lower_polygon, then the polys must be swapped.
                        if lower_polygon[-1] == front[front_index]:
                            tmp = lower_polygon, upper_polygon
                            upper_polygon, lower_polygon = tmp
                            #debug('    Swap upper/lower polygons')
                        else:
                            assert upper_polygon[-1] == front[front_index]

                    else:
                        assert next_tri is not None

            else:  # mode == 2
                # At each iteration, we require:
                #   * front_index is the starting index of the edge _preceding_
                #     the edge that will be handled in this iteration
                #   * lower_polygon is the polygon to which points should be
                #     added while traversing the front

                front_index += front_dir
                #debug("    Increment front_index: %d" % front_index)
                next_edge = (front[front_index], front[front_index+front_dir])
                #debug("    Set next_edge: %s" % repr(next_edge))

                assert front_index >= 0
                if front[front_index] == j:
                    # found endpoint!
                    #debug("    -> hit endpoint!")
                    lower_polygon.append(j)
                    upper_polygon.append(j)
                    #debug("    Appended to upper_polygon:", upper_polygon)
                    #debug("    Appended to lower_polygon:", lower_polygon)
                    break

                # Add point to lower_polygon.
                # The conditional is because there are cases where the
                # point was already added if we just crossed from mode 1.
                if lower_polygon[-1] != front[front_index]:
                    lower_polygon.append(front[front_index])
                    #debug("    Appended to lower_polygon:", lower_polygon)

                front_holes.append(front_index)
                #debug("    Append to front_holes:", front_holes)

                if self.edges_intersect((i, j), next_edge):
                    # crossing over front into triangle
                    #debug("    -> crossed over front, prepare for mode 1")
                    mode = 1

                    last_edge = next_edge
                    #debug("    Set last_edge:", last_edge)

                    # we are crossing the front, so this edge only has one
                    # triangle.
                    next_tri = self.tri_from_edge(last_edge)
                    #debug("    Set next_tri:", next_tri)

                    upper_polygon.append(front[front_index+front_dir])
                    #debug("    Appended to upper_polygon:", upper_polygon)
                #else:
                    #debug("    -> did not cross front..")

        #debug("Finished edge_event:")
        #debug("  front_holes:", front_holes)
        #debug("  upper_polygon:", upper_polygon)
        #debug("  lower_polygon:", lower_polygon)

        # (iii) triangluate empty areas

        #debug("Filling edge_event polygons...")
        for polygon in [lower_polygon, upper_polygon]:
            dist = self.distances_from_line((i, j), polygon)
            #debug("Distances:", dist)
            while len(polygon) > 2:
                ind = np.argmax(dist)
                #debug("Next index: %d" % ind)
                self.add_tri(polygon[ind], polygon[ind-1],
                             polygon[ind+1], legal=False,
                             source='edge_event')
                polygon.pop(ind)
                dist.pop(ind)

        #debug("Finished filling edge_event polygons.")

        # update front by removing points in the holes (places where front
        # passes below the cut edge)
        front_holes.sort(reverse=True)
        for i in front_holes:
            front.pop(i)

        #debug("Finished updating front after edge_event.")

    def find_cut_triangle(self, edge):
        """
        Return the triangle that has edge[0] as one of its vertices and is
        bisected by edge.

        Return None if no triangle is found.
        """
        edges = []  # opposite edge for each triangle attached to edge[0]
        for tri in self.tris:
            if edge[0] in tri:
                edges.append(self.edge_opposite_point(tri, edge[0]))

        for oedge in edges:
            o1 = self.orientation(edge, oedge[0])
            o2 = self.orientation(edge, oedge[1])
            #debug(edge, oedge, o1, o2)
            #debug(self.pts[np.array(edge)])
            #debug(self.pts[np.array(oedge)])
            if o1 != o2:
                return (edge[0], oedge[0], oedge[1])

        return None

    def edge_in_front(self, edge):
        """ Return the index where *edge* appears in the current front.
        If the edge is not in the front, return -1
        """
        e = (list(edge), list(edge)[::-1])
        for i in range(len(self.front)-1):
            if self.front[i:i+2] in e:
                return i
        return -1

    def edge_opposite_point(self, tri, i):
        """ Given a triangle, return the edge that is opposite point i.
        Vertexes are returned in the same orientation as in tri.
        """
        ind = tri.index(i)
        return (tri[(ind+1) % 3], tri[(ind+2) % 3])

    def adjacent_tri(self, edge, i):
        """
        Given a triangle formed by edge and i, return the triangle that shares
        edge. *i* may be either a point or the entire triangle.
        """
        if not np.isscalar(i):
            i = [x for x in i if x not in edge][0]

        try:
            pt1 = self.edges_lookup[edge]
            pt2 = self.edges_lookup[(edge[1], edge[0])]
        except KeyError:
            return None

        if pt1 == i:
            return (edge[1], edge[0], pt2)
        elif pt2 == i:
            return (edge[1], edge[0], pt1)
        else:
            raise RuntimeError("Edge %s and point %d do not form a triangle "
                               "in this mesh." % (edge, i))

    def tri_from_edge(self, edge):
        """Return the only tri that contains *edge*. If two tris share this
        edge, raise an exception.
        """
        edge = tuple(edge)
        p1 = self.edges_lookup.get(edge, None)
        p2 = self.edges_lookup.get(edge[::-1], None)
        if p1 is None:
            if p2 is None:
                raise RuntimeError("No tris connected to edge %r" % (edge,))
            return edge + (p2,)
        elif p2 is None:
            return edge + (p1,)
        else:
            raise RuntimeError("Two triangles connected to edge %r" % (edge,))

    def edges_in_tri_except(self, tri, edge):
        """Return the edges in *tri*, excluding *edge*.
        """
        edges = [(tri[i], tri[(i+1) % 3]) for i in range(3)]
        try:
            edges.remove(tuple(edge))
        except ValueError:
            edges.remove(tuple(edge[::-1]))
        return edges

    def edge_below_front(self, edge, front_index):
        """Return True if *edge* is below the current front.

        One of the points in *edge* must be _on_ the front, at *front_index*.
        """
        f0 = self.front[front_index-1]
        f1 = self.front[front_index+1]
        return (self.orientation(edge, f0) > 0 and
                self.orientation(edge, f1) < 0)

    def is_constraining_edge(self, edge):
        mask1 = self.edges == edge[0]
        mask2 = self.edges == edge[1]
        return (np.any(mask1[:, 0] & mask2[:, 1]) or
                np.any(mask2[:, 0] & mask1[:, 1]))

    def intersected_edge(self, edges, cut_edge):
        """ Given a list of *edges*, return the first that is intersected by
        *cut_edge*.
        """
        for edge in edges:
            if self.edges_intersect(edge, cut_edge):
                return edge

    def find_edge_intersections(self):
        """
        Return a dictionary containing, for each edge in self.edges, a list
        of the positions at which the edge should be split.
        """
        edges = self.pts[self.edges]
        cuts = {}  # { edge: [(intercept, point), ...], ... }
        for i in range(edges.shape[0]-1):
            # intersection of edge i onto all others
            int1 = self.intersect_edge_arrays(edges[i:i+1], edges[i+1:])
            # intersection of all edges onto edge i
            int2 = self.intersect_edge_arrays(edges[i+1:], edges[i:i+1])

            # select for pairs that intersect
            err = np.geterr()
            np.seterr(divide='ignore', invalid='ignore')
            try:
                mask1 = (int1 >= 0) & (int1 <= 1)
                mask2 = (int2 >= 0) & (int2 <= 1)
                mask3 = mask1 & mask2  # all intersections
            finally:
                np.seterr(**err)

            # compute points of intersection
            inds = np.argwhere(mask3)[:, 0]
            if len(inds) == 0:
                continue
            h = int2[inds][:, np.newaxis]
            pts = (edges[i, 0][np.newaxis, :] * (1.0 - h) +
                   edges[i, 1][np.newaxis, :] * h)

            # record for all edges the location of cut points
            edge_cuts = cuts.setdefault(i, [])
            for j, ind in enumerate(inds):
                if 0 < int2[ind] < 1:
                    edge_cuts.append((int2[ind], pts[j]))
                if 0 < int1[ind] < 1:
                    other_cuts = cuts.setdefault(ind+i+1, [])
                    other_cuts.append((int1[ind], pts[j]))

        # sort all cut lists by intercept, remove duplicates
        for k, v in cuts.items():
            v.sort(key=lambda x: x[0])
            for i in range(len(v)-2, -1, -1):
                if v[i][0] == v[i+1][0]:
                    v.pop(i+1)
        return cuts

    def split_intersecting_edges(self):
        # we can do all intersections at once, but this has excessive memory
        # overhead.
        #int1 = self.intersection_matrix(edges)
        #int2 = int1.T

        # measure intersection point between all pairs of edges
        all_cuts = self.find_edge_intersections()

        # cut edges at each intersection
        add_pts = []
        add_edges = []
        for edge, cuts in all_cuts.items():
            if len(cuts) == 0:
                continue

            #debug("Edge intersections:", edge, self.edges[edge],
            #      self.pts[self.edges[edge]], cuts)

            # add new points
            pt_offset = self.pts.shape[0] + len(add_pts)
            new_pts = [x[1] for x in cuts]
            add_pts.extend(new_pts)
            #debug("Add new points:", new_pts)

            # list of point indexes for all new edges
            pt_indexes = list(range(pt_offset, pt_offset + len(cuts)))
            pt_indexes.append(self.edges[edge, 1])

            # modify original edge
            self.edges[edge, 1] = pt_indexes[0]

            # add new edges
            new_edges = [[pt_indexes[i-1], pt_indexes[i]]
                         for i in range(1, len(pt_indexes))]
            add_edges.extend(new_edges)

        #debug("Adding %d points and %d edges to remove intersections." %
        #      (len(add_pts), len(add_edges)))
        if add_pts:
            add_pts = np.array(add_pts, dtype=self.pts.dtype)
            self.pts = np.append(self.pts, add_pts, axis=0)
        if add_edges:
            add_edges = np.array(add_edges, dtype=self.edges.dtype)
            self.edges = np.append(self.edges, add_edges, axis=0)

    def merge_duplicate_points(self):
        # generate a list of all pairs (i,j) of identical points
        dups = []
        for i in range(self.pts.shape[0]-1):
            test_pt = self.pts[i:i+1]
            comp_pts = self.pts[i+1:]
            eq = test_pt == comp_pts
            eq = eq[:, 0] & eq[:, 1]
            for j in np.argwhere(eq)[:, 0]:
                dups.append((i, i+1+j))

        dups_arr = np.array(dups)
        # remove duplicate points
        pt_mask = np.ones(self.pts.shape[0], dtype=bool)
        for i, inds in enumerate(dups_arr):
            # remove j from points
            # (note we pull the index from the original dups instead of
            # dups_arr because the indexes in pt_mask do not change)
            pt_mask[dups[i][1]] = False

            i, j = inds

            # rewrite edges to use i instead of j
            self.edges[self.edges == j] = i

            #assert not np.any(self.edges[:,0] == self.edges[:,1])

            # decrement all point indexes > j
            self.edges[self.edges > j] -= 1
            dups_arr[dups_arr > j] -= 1
            #assert not np.any(self.edges[:,0] == self.edges[:,1])

        self.pts = self.pts[pt_mask]

        # remove zero-length edges
        mask = self.edges[:, 0] != self.edges[:, 1]
        self.edges = self.edges[mask]

    def distance(self, A, B):
        # Distance between points A and B
        n = len(A)
        assert len(B) == n
        return np.linalg.norm(np.array(list(A)) - np.array(list(B)))

    def distances_from_line(self, edge, points):
        # Distance of a set of points from a given line
        #debug("distance from %r to %r" % (points, edge))
        e1 = self.pts[edge[0]]
        e2 = self.pts[edge[1]]
        distances = []
        for i in points:
            p = self.pts[i]
            proj = self.projection(e1, p, e2)
            distances.append(((p - proj)**2).sum()**0.5)
        assert distances[0] == 0 and distances[-1] == 0
        return distances

    def projection(self, a, b, c):
        """Return projection of (a,b) onto (a,c)
        Arguments are point locations, not indexes.
        """
        ab = b - a
        ac = c - a
        return a + ((ab*ac).sum() / (ac*ac).sum()) * ac

    def cosine(self, A, B, C):
        # Cosine of angle ABC
        a = ((C - B)**2).sum()
        b = ((C - A)**2).sum()
        c = ((B - A)**2).sum()
        d = (a + c - b) / ((4 * a * c)**0.5)
        return d

    #def barycentric(self, A, B, C, p, q, r):
        ## Cartesian coordinates of the point whose barycentric coordinates
        ## with respect to the triangle ABC are [p,q,r]
        #n = len(A)
        #assert len(B) == len(C) == n
        #s = p+q+r
        #p, q, r = p/s, q/s, r/s
        #return tuple([p*A[i]+q*B[i]+r*C[i] for i in range(n)])

    #def trilinear(self, A, B, C, alpha, beta, gamma):
        ## Cartesian coordinates of the point whose trilinear coordinates
        ## with respect to the triangle ABC are [alpha,beta,gamma]
        #a = distance(B, C)
        #b = distance(A, C)
        #c = distance(A, B)
        #return barycentric(A, B, C, a*alpha, b*beta, c*gamma)

    #def circuminfo(self, A, B, C):
        ## Cartesian coordinates of the circumcenter of triangle ABC
        #cosA = cosine(C, A, B)
        #cosB = cosine(A, B, C)
        #cosC = cosine(B, C, A)
        #cc = trilinear(A, B, C, cosA, cosB, cosC)
        ## returns circumcenter and circumradius
        #return cc, distance(cc, A)

    def iscounterclockwise(self, a, b, c):
        # Check if the points lie in counter-clockwise order or not
        A = self.pts[a]
        B = self.pts[b]
        C = self.pts[c]
        return np.cross(B-A, C-B) > 0

    def edges_intersect(self, edge1, edge2):
        """
        Return 1 if edges intersect completely (endpoints excluded)
        """
        h12 = self.intersect_edge_arrays(self.pts[np.array(edge1)],
                                         self.pts[np.array(edge2)])
        h21 = self.intersect_edge_arrays(self.pts[np.array(edge2)],
                                         self.pts[np.array(edge1)])
        err = np.geterr()
        np.seterr(divide='ignore', invalid='ignore')
        try:
            out = (0 < h12 < 1) and (0 < h21 < 1)
        finally:
            np.seterr(**err)
        return out

    def intersection_matrix(self, lines):
        """
        Return a 2D array of intercepts such that
        intercepts[i, j] is the intercept of lines[i] onto lines[j].

        *lines* must be an array of point locations with shape (N, 2, 2), where
        the axes are (lines, points_per_line, xy_per_point).

        The intercept is described in intersect_edge_arrays().
        """
        return self.intersect_edge_arrays(lines[:, np.newaxis, ...],
                                          lines[np.newaxis, ...])

    def intersect_edge_arrays(self, lines1, lines2):
        """Return the intercepts of all lines defined in *lines1* as they
        intersect all lines in *lines2*.

        Arguments are of shape (..., 2, 2), where axes are:

        0: number of lines
        1: two points per line
        2: x,y pair per point

        Lines are compared elementwise across the arrays (lines1[i] is compared
        against lines2[i]). If one of the arrays has N=1, then that line is
        compared against all lines in the other array.

        Returns an array of shape (N,) where each value indicates the intercept
        relative to the defined line segment. A value of 0 indicates
        intersection at the first endpoint, and a value of 1 indicates
        intersection at the second endpoint. Values between 1 and 0 are on the
        segment, whereas values outside 1 and 0 are off of the segment.
        """
        # vector for each line in lines1
        l1 = lines1[..., 1, :] - lines1[..., 0, :]
        # vector for each line in lines2
        l2 = lines2[..., 1, :] - lines2[..., 0, :]
        # vector between first point of each line
        diff = lines1[..., 0, :] - lines2[..., 0, :]

        p = l1.copy()[..., ::-1]  # vectors perpendicular to l1
        p[..., 0] *= -1

        f = (l2 * p).sum(axis=-1)  # l2 dot p
        # tempting, but bad idea!
        #f = np.where(f==0, 1, f)
        err = np.geterr()
        np.seterr(divide='ignore', invalid='ignore')
        try:
            h = (diff * p).sum(axis=-1) / f  # diff dot p / f
        finally:
            np.seterr(**err)

        return h

    def orientation(self, edge, point):
        """ Returns +1 if edge[0]->point is clockwise from edge[0]->edge[1],
        -1 if counterclockwise, and 0 if parallel.
        """
        v1 = self.pts[point] - self.pts[edge[0]]
        v2 = self.pts[edge[1]] - self.pts[edge[0]]
        c = np.cross(v1, v2)  # positive if v1 is CW from v2
        return 1 if c > 0 else (-1 if c < 0 else 0)

    #def legalize(self, p):
        ### Legalize recursively - incomplete
        #return p  # disabled for now

        #f00, f11, p = p

        #debug("Legalizing points = {}, {}, {}".format(f00, f11, p))
        #a = pts[f00]
        #b = pts[f11]
        #c = pts[p]
        #cc, cr = circuminfo(a, b, c)
        #for point in pts:
        #    if np.all(point == a) or np.all(point == b) or np.all(point == c):
        #        continue
        #    elif distance(cc, point) < cr:
        #        #debug("Illegal point")
        #        #debug(point)
        #        pass

        #return (f00, f11, p)

    def add_tri(self, a, b, c, legal=True, source=None):
        # source is just used for #debugging
        #debug("Add triangle [%s]:" % source, (a, b, c))

        # sanity check
        assert a != b and b != c and c != a

        # ignore flat tris
        pa = self.pts[a]
        pb = self.pts[b]
        pc = self.pts[c]
        if np.all(pa == pb) or np.all(pb == pc) or np.all(pc == pa):
            #debug("   Triangle is flat; refusing to add.")
            return

        # check this tri is unique
        for t in permutations((a, b, c)):
            if t in self.tris:
                raise Exception("Cannot add %s; already have %s" %
                                ((a, b, c), t))

        # TODO: should add to edges_lookup after legalization??
        if self.iscounterclockwise(a, b, c):
            #debug("    ", (a, b), (b, c), (c, a))
            assert (a, b) not in self.edges_lookup
            assert (b, c) not in self.edges_lookup
            assert (c, a) not in self.edges_lookup
            self.edges_lookup[(a, b)] = c
            self.edges_lookup[(b, c)] = a
            self.edges_lookup[(c, a)] = b
        else:
            #debug("    ", (b, a), (c, b), (a, c))
            assert (b, a) not in self.edges_lookup
            assert (c, b) not in self.edges_lookup
            assert (a, c) not in self.edges_lookup
            self.edges_lookup[(b, a)] = c
            self.edges_lookup[(c, b)] = a
            self.edges_lookup[(a, c)] = b

        #if legal:
            #tri = self.legalize((a, b, c))
        #else:
        tri = (a, b, c)

        self.tris[tri] = None

    def remove_tri(self, a, b, c):
        #debug("Remove triangle:", (a, b, c))

        for k in permutations((a, b, c)):
            if k in self.tris:
                break
        del self.tris[k]
        (a, b, c) = k

        if self.edges_lookup.get((a, b), -1) == c:
            #debug("    ", (a,b), (b,c), (c,a))
            del self.edges_lookup[(a, b)]
            del self.edges_lookup[(b, c)]
            del self.edges_lookup[(c, a)]
        elif self.edges_lookup.get((b, a), -1) == c:
            #debug("    ", (b,a), (c,b), (a,c))
            del self.edges_lookup[(b, a)]
            del self.edges_lookup[(a, c)]
            del self.edges_lookup[(c, b)]
        else:
            raise RuntimeError("Lost edges_lookup for tri (%d, %d, %d)" %
                               (a, b, c))

        return k

def _triangulate_python(vertices_2d, segments):
    segments = segments.reshape(len(segments) / 2, 2)
    T = Triangulation(vertices_2d, segments)
    T.triangulate()
    vertices_2d = T.pts
    triangles = T.tris.ravel()
    return vertices_2d, triangles

def _triangulate_cpp(vertices_2d, segments):
    T = triangle.triangulate({'vertices': vertices_2d,
                              'segments': segments}, "p")
    vertices_2d = T["vertices"]
    triangles = T["triangles"]
    return vertices_2d, triangles

def triangulate(vertices):
    """Triangulate a set of vertices. Returns a pair (vertices, triangles)."""
    n = len(vertices)
    vertices = np.asarray(vertices)
    zmean = vertices[:, 2].mean()
    vertices_2d = vertices[:, :2]
    segments = np.repeat(np.arange(n + 1), 2)[1:-1]
    segments[-2:] = n - 1, 0

    if _TRIANGLE_AVAILABLE:
        vertices_2d, triangles = _triangulate_cpp(vertices_2d, segments)
    else:
        vertices_2d, triangles = _triangulate_python(vertices_2d, segments)

    vertices = np.empty((len(vertices_2d), 3))
    vertices[:, :2] = vertices_2d
    vertices[:, 2] = zmean
    return vertices, triangles

# Note: using custom #debug instead of logging because
# there are MANY messages and logger might be too expensive.
# After this becomes stable, we might just remove them altogether.
def debug(*args):
    print(*args)

max_stars_repo_name: izaid/vispy

max_issues_repo_name: izaid/vispy

max_forks_repo_name: izaid/vispy


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

import os
import cv2
import numpy as np
from common import config
import glob
from utils import from_img_to_sub_pixel, from_sub_pixel_to_img
from PIL import Image

mat1 = np.array([ 65.481, 128.553, 24.966 ])
mat2 = np.array([-37.797, -74.203, 112.0  ])
mat3 = np.array([ 112.0, -93.786, -18.214])

def rgb2ycbcr(img):
    y = np.round(np.dot(img, mat1) + 16.0)
    cb = np.round(np.dot(img, mat2) + 128.0)
    cr = np.round(np.dot(img, mat3) + 128.0)
    return np.concatenate([y[...,np.newaxis], cb[...,np.newaxis], cr[...,np.newaxis]], axis = 2)

class Dataset():

    def __init__(self, dataset_name):

        test_dataset_path = ['../../dataset/Set14/image_SRF_3']
        train_dataset_path = '../../dataset/Train96'
        self.minibatch_size = config.minibatch_size
        self.ds_name = dataset_name
        self.rng = np.random
        self.ratio = config.ratio
        self.edge = config.edge
        self.stride = config.stride
        self.patch_size = config.patch_size
        train_list = glob.glob(train_dataset_path + '/*.jpg')
        test_list = []
        for single_test_dataset_path in test_dataset_path:
            test_list = test_list + glob.glob(single_test_dataset_path + '/*LR.png')
        self.dataset_meta = {
            'train': (train_list, 4992),
            'test': (test_list, len(test_list)),
        }

    def load(self):

        self.instances = self.dataset_meta[self.ds_name][1]
        self.images = self.dataset_meta[self.ds_name][0]
        self.data_number = len(self.images)
        return self

    @property
    def instances_per_epoch(self):
        return self.instances

    @property
    def minibatchs_per_epoch(self):
        return self.instances // config.minibatch_size

    @property
    def testing_minibatchs_per_epoch(self):
        return self.instances

    def instance_generator(self):
        for i in range(self.instances):
            idx = i % self.data_number
            img = cv2.imread(self.images[idx])
            if self.ds_name == 'test':
                img = rgb2ycbcr(img[:, :, ::-1]/255.0)
                hr_img = rgb2ycbcr(cv2.imread(self.images[idx].split('LR')[0] + 'HR.png')[:, :, ::-1]/255.0)
                yield img.astype(np.float32)/255.0, hr_img.astype(np.float32)/255.0
            elif self.ds_name == 'train':
                h, w = img.shape[:2]
                new_h = h - h % self.ratio
                new_w = w - w % self.ratio
                img = img[0: new_h, 0: new_w, :]
                blurred_img = cv2.GaussianBlur(img, (1,1), 0)
                blurred_img = Image.fromarray(blurred_img.astype(np.uint8))
                lr_img = blurred_img.resize((new_w//config.ratio, new_h//config.ratio),Image.BICUBIC)
                lr_img = np.array(lr_img)
                lr_h, lr_w = lr_img.shape[:2]
                x_start = self.rng.randint(0, (lr_w - self.patch_size)//self.stride + 1 ) * config.stride
                y_start = self.rng.randint(0, (lr_h - self.patch_size)//self.stride + 1 ) * config.stride
                patch = lr_img[y_start:y_start+config.patch_size,x_start:x_start+config.patch_size,:].astype(np.float32)/255.0

                hr_x_start = self.ratio * (x_start + self.edge//2)
                hr_x_end = self.ratio * (x_start + self.patch_size - self.edge//2)
                hr_y_start = self.ratio * (y_start + self.edge//2)
                hr_y_end = self.ratio * (y_start + self.patch_size - self.edge//2)

                hr_img = rgb2ycbcr(img[hr_y_start:hr_y_end, hr_x_start:hr_x_end, ::-1]/255.0)
                hr_patch = from_img_to_sub_pixel(hr_img, self.ratio)
                patch = rgb2ycbcr(patch[:, :, ::-1])

                yield patch.astype(np.float32)/255.0, hr_patch.astype(np.float32)/255.0

if __name__ == "__main__":
    ds = Dataset('train')
    gen = ds.load().instance_generator()

    imggrid = []
    while True:
        for i in range(8):
            img, hr_patch = next(gen)
            img = cv2.resize(img, ((config.patch_size) * config.ratio, (config.patch_size) * config.ratio))
            start_idx = int(config.edge /2 * config.ratio)
            end_idx = int((config.patch_size  - config.edge/2) * config.ratio)
            imggrid.append(img[start_idx:end_idx, start_idx:end_idx, :])
            imggrid.append(from_sub_pixel_to_img(hr_patch, config.ratio))

        imggrid = np.array(imggrid).reshape((4, 4, hr_patch.shape[0] * config.ratio, hr_patch.shape[1] * config.ratio, 3))
        imggrid = imggrid.transpose((0, 2, 1, 3, 4)).reshape((4*hr_patch.shape[0] * config.ratio, 4*hr_patch.shape[1] * config.ratio, 3))
        cv2.imshow('', imggrid)
        c = chr(cv2.waitKey(0) & 0xff)
        if c == 'q':
            exit()
        imggrid = []

max_stars_repo_name: PeiqinSun/tf-tutorials

max_issues_repo_name: megvii-research/tf-tutorials

max_forks_repo_name: PeiqinSun/tf-tutorials


text:
# Copyright 2018 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     https://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.

"""Evaluation functions.
"""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from lib import data
import numpy as np
import scipy.spatial

def closest_line(query_lines, metric='cosine'):
    """Compute the distance to, and parameters for, the closest line to each
    line in query_lines.

    Args:
        - query_lines: Array of lines to compute closest matches for, shape
            (n_lines, width, height, 1)
        - metric: String to pass to scipy.spatial.distance.cdist to choose
            which distance metric to use
    Returns:
        - min_dist, starts, ends: Arrays of shape (n_lines,) denoting the
            distance to the nearest ``true'' line and the start and end points.
    """
    h, w = query_lines.shape[1:-1]
    # Construct 10000 lines with these dimensions
    angles = np.linspace(0, 2*np.pi - 2*np.pi/10000, 10000)
    all_lines = np.array(
        [(data.draw_line(angle, h, w)) for angle in angles])
    # Produce vectorized versions of both for use with scipy.spatial
    flat_query = query_lines.reshape(query_lines.shape[0], -1)
    flat_all = all_lines.reshape(all_lines.shape[0], -1)
    # Compute pairwise distance matrix of query lines with all valid lines
    distances = scipy.spatial.distance.cdist(flat_query, flat_all, metric)
    min_dist_idx = np.argmin(distances, axis=-1)
    min_dist = distances[np.arange(distances.shape[0]), min_dist_idx]
    angles = np.array([angles[n] for n in min_dist_idx])
    return min_dist, angles

def smoothness_score(angles):
    """Computes the smoothness score of a line interpolation according to the
    angles of each line.

    Args:
        - angles: Array of shape (n_interpolations, n_lines_per_interpolation)
            giving the angle of each line in each interpolation.

    Returns:
        - smoothness_scores: Array of shape (n_interpolations,) giving the
            average smoothness score for all of the provided interpolations.
    """
    angles = np.atleast_2d(angles)
    # Remove discontinuities larger than np.pi
    angles = np.unwrap(angles)
    diffs = np.abs(np.diff(angles, axis=-1))
    # Compute the angle difference from the first and last point
    total_diff = np.abs(angles[:, :1] - angles[:, -1:])
    # When total_diff is zero, there's no way to compute this score
    zero_diff = (total_diff < 1e-4).flatten()
    normalized_diffs = diffs/total_diff
    deviation = np.max(normalized_diffs, axis=-1) - 1./(angles.shape[1] - 1)
    # Set score to NaN when we aren't able to compute it
    deviation[zero_diff] = np.nan
    return deviation

def line_eval(interpolated_lines):
    """Given a group of line interpolations, compute mean nearest line distance
    and mean smoothness score for all of the interpolations.

    This version of this metric is meant for vertical lines only.

    Args:
        - interpolated_lines: Collection of line interpolation images, shape
            (n_interpolations, n_lines_per_interpolation, height, width, 1)

    Returns:
        - mean_distance: Average distance to closest ``real'' line.
        - mean_smoothness: Average interpolation smoothness

    """
    original_shape = interpolated_lines.shape
    min_dist, angles = closest_line(
        interpolated_lines.reshape((-1,) + original_shape[2:]))
    mean_distance = np.mean(min_dist)
    smoothness_scores = smoothness_score(
        angles.reshape(original_shape[0], original_shape[1]))
    nan_scores = np.isnan(smoothness_scores)
    # If all scores were NaN, set the mean score to NaN
    if np.all(nan_scores):
        mean_smoothness = np.nan
    # Otherwise only compute mean for non-NaN scores
    else:
        sum_smoothness = np.sum(smoothness_scores[np.logical_not(nan_scores)])
        mean_smoothness = sum_smoothness/float(len(nan_scores))
    return np.float32(mean_distance), np.float32(mean_smoothness)

max_stars_repo_name: hsukyle/acai

max_issues_repo_name: hsukyle/acai

max_forks_repo_name: hsukyle/acai

reason: import numpy,import scipy


text:
#!/usr/bin/env /usr/bin/python3
import numpy as np
import argparse
from matplotlib import pyplot as plt
from pathlib import Path
from scipy import optimize

################################################################################
#===============================================================================
# plot_output.py
#===============================================================================
################################################################################

def smooth_runaway (x):
	return np.where(x>0, np.exp(-1/x), 0)

def smooth_step (x):
	return smooth_runaway(x)/(smooth_runaway(x)+smooth_runaway(1-x))

def smooth_k_step (x,k): # k should be >=1
	return np.where(x>0, np.where(x<1,
		0.5*(np.tanh(k*(2*x-1)/2/np.sqrt(x*(1-x)))+1), 1), 0)

################################################################################

def F (x, x0, alpha, m, b, c):
#	return smooth_step((x0-x)/alpha) * G(x, m, b, c)
	return ((np.tanh((x0-x)/alpha) + 1)/2)**2 * G(x, m, b, c)

def G (x, m, b, c):
	return -m*x+b + np.sqrt( (-m*x+b)**2 + c**2 )

################################################################################

def plot_output (result_file = Path('../flat_output/results.csv'),
				 output_file = Path('../plots/flat_output_plot.svg'),
				 return_ax = False):
	fig, ax = plt.subplots(1)
	x_values = np.linspace(2.8,4.2,11)
	ax.plot(x_values, np.zeros_like(x_values),
			color = 'darkgrey', linestyle = '-', marker = '')
	y_values = np.linspace(0.,1.7,11)
	ax.plot(np.ones_like(y_values)*3.813, y_values,
			color = 'black', linestyle = '--', marker = '')
		#	label = '3.813')
	x = np.linspace(3,4,1000)
	energy_data = np.loadtxt(result_file, comments = '#', delimiter = ',')
	r_param = energy_data[:,0]
	p0_param = energy_data[:,1]
	energy = energy_data[:,2]
	error = energy_data[:,3]
	r_params = np.unique(r_param)[::-1]
	markers = ['*', 'p', 's', 'v', 'o']
	colors = ['purple', 'blue', 'green', 'red', 'black']
	p0 = [3.8, 0.08, 1.0, 3.6, 0.1]
	test = np.array([])
	for index, current_r in enumerate(r_params):
		ax.plot(p0_param[r_param == current_r],
				energy[r_param == current_r],
				color = colors[index], linestyle = '',
				marker = markers[index], markersize = '4.',
				label = f'r = {current_r:1.1f}')
#		ax.errorbar(p0_param[r_param == current_r],
#					energy[r_param == current_r],
#					error[r_param == current_r],
#					color = colors[index], linestyle = '',
#					marker = markers[index], markersize = '4.',
#					label = f'r = {current_r:1.1f}')
		params, covar = optimize.curve_fit(F, p0_param[r_param == current_r],
											  energy[r_param == current_r],
											  p0 = p0,
											  bounds = (0,24))
	#	if current_r == 0.5:
	#		print(params)
		p0 = params
		test = np.append(test,p0[0])
		y = F (x,*params)
		ax.plot(x,y, color = colors[index], linestyle = '-')
	print(np.mean(test))
	print(np.std(test))
	ax.grid(True)
	ax.set_xlim((2.95,4.05))
	ax.set_ylim((-0.05,1.65))
	ax.set_xticks(np.linspace(3.0,4.0,11))
#	ax.set_title('Energy Barrier as a Function of Shape Index')
	ax.set_xlabel(r'$p_0$')
	ax.set_ylabel(r'$\Delta \varepsilon$')
	ax.legend()
	plt.savefig(output_file)
	plt.rc('pgf', texsystem='pdflatex')
	plt.savefig(output_file.with_suffix('.pgf'))
	if return_ax:
		return ax
#	plt.show()

################################################################################

if __name__ == '__main__':
	# Use argparse to get arguements from commandline call
	parser = argparse.ArgumentParser(
							description = 'Plot results of simulation.')
	parser.add_argument('-o', '--outfile',
						nargs = 1,
						default = ['../plots/flat_output_plot.svg'],
						type = str,
						required = False,
						help = 'file to put plot in')
	parser.add_argument('result_file',
						nargs = '?',
						default = '../flat_output/results.csv',
						type = str,
						help = 'file with simulation result data')
	args = parser.parse_args()
	plot_output( result_file = Path(args.result_file),
				 output_file = Path(args.outfile[0]) )

################################################################################
# EOF

reason: import numpy,from scipy


text:
import sys
sys.path.append("/usr/local/cython/")

import numpy as np
import time
import pandas as pd

from leopy.optim.nav2d_factors import *
# from leopy.optim.push2d_factors import *

from leopy.utils import tf_utils, dir_utils, vis_utils
from leopy.utils.logger import Logger

import logging
log = logging.getLogger(__name__)
class GraphOpt():
    def __init__(self):

        self.graph = gtsam.NonlinearFactorGraph()
        self.init_vals = gtsam.Values()

        self.optimizer = self.init_isam2()
        self.est_vals = gtsam.Values()

        self.logger = Logger()
        self.dataframe = pd.DataFrame()
        # self.graph_full = gtsam.NonlinearFactorGraph()

    def init_isam2(self):
        params_isam2 = gtsam.ISAM2Params()
        params_isam2.setRelinearizeThreshold(0.01)
        params_isam2.setRelinearizeSkip(10)

        return gtsam.ISAM2(params_isam2)

def init_vars_step(graphopt, tstep, dataset_type):

    if (dataset_type == "push2d"):
        key_tm1 = gtsam.symbol(ord('o'), tstep-1)
        key_t = gtsam.symbol(ord('o'), tstep)
        graphopt.init_vals.insert(key_t, graphopt.est_vals.atPose2(key_tm1))

        key_tm1 = gtsam.symbol(ord('e'), tstep-1)
        key_t = gtsam.symbol(ord('e'), tstep)
        graphopt.init_vals.insert(key_t, graphopt.est_vals.atPose2(key_tm1))
    else:
        key_tm1 = gtsam.symbol(ord('x'), tstep-1)
        key_t = gtsam.symbol(ord('x'), tstep)
        graphopt.init_vals.insert(key_t, graphopt.est_vals.atPose2(key_tm1))

    return graphopt

def reset_graph(graphopt):
    # graphopt.graph_full.push_back(graphopt.graph)

    graphopt.graph.resize(0)
    graphopt.init_vals.clear()

    return graphopt

def optimizer_update(graphopt):
    graphopt.optimizer.update(graphopt.graph, graphopt.init_vals)
    graphopt.est_vals = graphopt.optimizer.calculateEstimate()

    return graphopt

def print_step(tstep, data, isam2_estimate, dataset_type):

    if (dataset_type == "push2d"):
        key = gtsam.symbol(ord('o'), tstep)
        print('Estimated, Grountruth pose: \n {0} \n {1} '.format(
            tf_utils.pose2_to_vec3(isam2_estimate.atPose2(key)), data['obj_poses_gt'][tstep]))
    else:
        key = gtsam.symbol(ord('x'), tstep)
        print('Estimated, Grountruth pose: \n {0} \n {1} '.format(
            tf_utils.pose2_to_vec3(isam2_estimate.atPose2(key)), data['poses_gt'][tstep]))

def log_step(tstep, logger, data, optimizer, dataset_type):

    if (dataset_type == 'nav2d'):
        logger = log_step_nav2d(tstep, logger, data, optimizer)
    elif (dataset_type == 'push2d'):
        logger = log_step_push2d(tstep, logger, data, optimizer)

    return logger

def log_step_df(tstep, dataframe, data, optimizer, dataset_type):

    if (dataset_type == 'nav2d'):
        logger = log_step_df_nav2d(tstep, dataframe, data, optimizer)
    elif (dataset_type == 'push2d'):
        logger = log_step_df_push2d(tstep, dataframe, data, optimizer)

    return logger

def add_first_pose_priors(graphopt, data, dataset_type):

    prior_cov = np.array([1e-6, 1e-6, 1e-6])

    if(dataset_type == 'nav2d'):
        keys = [gtsam.symbol(ord('x'), 0)]
        graphopt.init_vals.insert(
            keys[0], tf_utils.vec3_to_pose2(data['poses_gt'][0]))
        graphopt.graph = add_unary_factor(
            graphopt.graph, keys, prior_cov, data['poses_gt'][0])
    elif (dataset_type == 'push2d'):
        keys = [gtsam.symbol(ord('o'), 0)]
        graphopt.init_vals.insert(
            keys[0], tf_utils.vec3_to_pose2(data['obj_poses_gt'][0]))
        graphopt.graph = add_unary_factor(
            graphopt.graph, keys, prior_cov, data['obj_poses_gt'][0])
        keys = [gtsam.symbol(ord('e'), 0)]
        graphopt.init_vals.insert(
            keys[0], tf_utils.vec3_to_pose2(data['ee_poses_gt'][0]))
        graphopt.graph = add_unary_factor(
            graphopt.graph, keys, prior_cov, data['ee_poses_gt'][0])

    graphopt = optimizer_update(graphopt)
    graphopt = reset_graph(graphopt)

    return graphopt

def get_optimizer_soln(tstep, graphopt, dataset_type, sampler=False):

    poses_graph = graphopt.optimizer.calculateEstimate()
    pose_vec_graph = np.zeros((poses_graph.size(), 3))

    num_poses = tstep + 1
    if (dataset_type == 'nav2d'):
        keys = [[gtsam.symbol(ord('x'), i)] for i in range(0, num_poses)]
    elif (dataset_type == 'push2d'):
        keys = [[gtsam.symbol(ord('o'), i), gtsam.symbol(ord('e'), i)] for i in range(0, num_poses)]

    key_vec = gtsam.gtsam.KeyVector()
    for key_idx in range(0, len(keys[0])):
        for pose_idx in range(0, num_poses):
            key = keys[pose_idx][key_idx]
            pose2d = poses_graph.atPose2(key)
            pose_vec_graph[key_idx * num_poses + pose_idx, :] = [pose2d.x(), pose2d.y(), pose2d.theta()]

            key_vec.push_back(key)

    mean = pose_vec_graph
    cov = None

    if sampler:
        marginals = gtsam.Marginals(graphopt.optimizer.getFactorsUnsafe(), poses_graph)
        cov = marginals.jointMarginalCovariance(key_vec).fullMatrix()

    return (mean, cov)

def run_optimizer(cost, params=None):
    ''' construct graph from data and optimize '''

    dataset_type = params.dataio.dataset_type

    graphopt = GraphOpt()
    graphopt = add_first_pose_priors(graphopt, cost.data, dataset_type)

    if (dataset_type == "push2d"):
        nsteps = np.minimum(params.optim.nsteps, len(cost.data['obj_poses_gt']))

        planar_sdf = PlanarSDF("{0}/{1}/sdf/{2}.json".format(
            params.BASE_PATH, params.dataio.srcdir_dataset, params.dataio.obj_shape))
        tactile_model = torch.jit.load("{0}/local/models/{1}.pt".format(
            params.BASE_PATH, params.tactile_model.name))
    else:
        nsteps = np.minimum(params.optim.nsteps, len(cost.data['poses_gt']))

    for tstep in range(1, nsteps):

        # filter out curr step factors
        factor_idxs = [idxs for idxs, keys in enumerate(
            cost.data['factor_keyids']) if (max(keys) == tstep)]

        factor_keysyms = [cost.data['factor_keysyms'][idx] for idx in factor_idxs]
        factor_keyids = [cost.data['factor_keyids'][idx] for idx in factor_idxs]
        factor_meas = [cost.data['factor_meas'][idx] for idx in factor_idxs]
        factor_names = [cost.data['factor_names'][idx] for idx in factor_idxs]

        # compute factor costs
        for idx in range(0, len(factor_idxs)):

            key_syms, key_ids = factor_keysyms[idx], factor_keyids[idx]

            keys = [gtsam.symbol(ord(key_syms[i]), key_id)
                    for i, key_id in enumerate(key_ids)]
            factor_name = factor_names[idx]

            sigma_inv_val = cost.theta.get_sigma_inv(factor_name, factor_meas[idx])
            sigma_inv_val = sigma_inv_val.detach().cpu().numpy()
            # factor_cov = np.reciprocal(sigma_inv_val) # factor_cov: sigma format
            factor_cov = np.reciprocal(np.sqrt(sigma_inv_val)) # factor_cov: sigma_sq format

            if (factor_name == 'gps'):
                graphopt.graph = add_unary_factor(
                    graphopt.graph, keys, factor_cov, factor_meas[idx][0:3])
            elif (factor_name == 'odom'):
                graphopt.graph = add_binary_odom_factor(
                    graphopt.graph, keys, factor_cov, factor_meas[idx][0:3])
            elif (factor_name == 'ee_pose_prior'):
                graphopt.graph = add_unary_factor(
                    graphopt.graph, keys, factor_cov, factor_meas[idx])
            elif (factor_name == 'qs_push_motion'):
                graphopt.graph = add_qs_motion_factor(
                    graphopt.graph, keys, factor_cov, params=params)
            elif (factor_name == 'sdf_intersection'):
                graphopt.graph = add_sdf_intersection_factor(
                    graphopt.graph, keys, factor_cov, planar_sdf.sdf)
            elif (factor_name == 'tactile_rel'):
                if (params.tactile_model.oracle == True):
                    tf_pred_net = tactile_oracle_output(cost.data, key_ids)
                else:
                    if (params.dataio.model_type == "fixed_cov_varying_mean"):
                        tf_pred_net = tactile_model_output(
                            cost, img_feats=factor_meas[idx], params=params.tactile_model) # learnable model weights
                        tf_pred_net[3] = np.clip(tf_pred_net[3], -1., 1.)
                    else:
                        tf_pred_net = tactile_model_output_fixed(
                            tactile_model, img_feats=factor_meas[idx], params=params.tactile_model) # fixed model weights
                graphopt.graph = add_tactile_rel_meas_factor(
                    graphopt.graph, keys, factor_cov, tf_pred_net, params=params.tactile_model)
            elif (factor_name == 'binary_interseq_obj'):
                graphopt.graph = add_binary_odom_factor(
                    graphopt.graph, keys, factor_cov, factor_meas[idx])

        # init variables
        graphopt = init_vars_step(graphopt, tstep, dataset_type)

        # optimize
        graphopt = optimizer_update(graphopt)

        # print step
        # print_step(tstep, cost.data, graphopt.est_vals, dataset_type)

        # log step
        # todo: phase out log_step and only use log_step_df. Currently log_step needed for plotting.
        graphopt.logger = log_step(
            tstep, graphopt.logger, cost.data, graphopt.optimizer, dataset_type)

        if params.logger.enable:
            graphopt.dataframe = log_step_df(
                tstep, graphopt.dataframe, cost.data, graphopt.optimizer, dataset_type)

        # vis step
        if (params.optim.vis_step):
            vis_utils.vis_step(tstep, graphopt.logger, params=params)

        # reset graph
        graphopt = reset_graph(graphopt)

    mean, cov = get_optimizer_soln(tstep, graphopt, dataset_type, params.leo.sampler)

    # log final solution
    if params.logger.enable:
        cov_log = cov.tolist() if cov is not None else cov
        data_row = pd.DataFrame({'opt/mean': [mean.tolist()], 'opt/covariance': [cov_log]}, index=[tstep], dtype=object)
        data_row.index.name = 'tstep'
        graphopt.dataframe = pd.concat([graphopt.dataframe, data_row], axis=1)

    if (params.optim.save_fig) | (params.optim.show_fig):
        vis_utils.vis_step(tstep, graphopt.logger, params=params)

    # if params.optim.verbose:
    #     graph_full = graphopt.optimizer.getFactorsUnsafe()
    #     graph_cost = graph_full.error(graphopt.optimizer.calculateEstimate())
    #     print("gtsam cost: {0}".format(graph_cost))
    #     # print("gtsam per factor cost for {0} factors over {1} time steps".format(
    #     # graph_full.size(), nsteps))
    #     # for i in range(graph_full.size()):
    #     #     print(graph_full.at(i).error(graphopt.optimizer.calculateEstimate()))

    return (mean, cov, graphopt.dataframe)

max_stars_repo_name: rpl-cmu/leo

max_issues_repo_name: psodhi/logo

max_forks_repo_name: psodhi/logo


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

# This file is part of the Kramers-Kronig Calculator software package.
#
# Copyright (c) 2013 Benjamin Watts, Daniel J. Lauk
#
# The software is licensed under the terms of the zlib/libpng license.
# For details see LICENSE.txt

"""This module implements the Kramers-Kronig transformation."""

import logging, sys
logger = logging.getLogger(__name__)
if __name__ == '__main__':
	logging.basicConfig(level=logging.DEBUG)
	logging.StreamHandler(stream=sys.stdout)

import math
import numpy
import os
import data

def calc_relativistic_correction(stoichiometry):
	"""Calculate the relativistic correction to the Kramers-Kronig transform.

	Parameters:
	-----------
	stoichiometry : array of integer/float pairs
		Each pair in the list consists of an atomic number and the relative proportion of that element.

	Returns
	-------
	This function returns a ``float`` holding the relativistic
	corection to the Kramers-Kronig transform.

	"""
	correction = 0
	for z, n in stoichiometry:
		correction += (z - (z/82.5)**2.37) * n
	return correction

def KK_General_PP(Eval_Energy, Energy, imaginary_spectrum, orders, relativistic_correction):
	"""Calculate Kramers-Kronig transform with "Piecewise Polynomial"
	algorithm plus the Biggs and Lighthill extended data.

	Parameters
	----------
	Eval_Energy : numpy vector of `float`
		Set of photon energies describing points at which to evaluate the real spectrum
	Energy : numpy vector of `float`
		Set of photon energies describing intervals for which each row of `imaginary_spectrum` is valid
	imaginary_spectrum : two-dimensional `numpy.array` of `float`
		The array consists of columns of polynomial coefficients belonging to the power terms indicated by 'order'
	orders : numpy vector of integers
		The vector represents the polynomial indices corresponding to the columns of imaginary_spectrum
	relativistic_correction : float
		The relativistic correction to the Kramers-Kronig transform.
		You can calculate the value using the `calc_relativistic_correction` function.

	Returns
	-------
	This function returns the real part of the scattering factors evaluated at photon energies specified by Eval_Energy.

	"""
	logger = logging.getLogger(__name__)
	logger.info("Calculate Kramers-Kronig transform using general piecewise-polynomial algorithm")
	# Need to build x-E-n arrays
	X = numpy.tile(Energy[:,numpy.newaxis,numpy.newaxis],(1,len(Eval_Energy),len(orders)))
	E = numpy.tile(Eval_Energy[numpy.newaxis,:,numpy.newaxis],(len(Energy)-1,1,len(orders)))
	C = numpy.tile(imaginary_spectrum[:,numpy.newaxis,:],(1,len(Eval_Energy),1))
	N = numpy.tile(orders[numpy.newaxis,numpy.newaxis,:],(len(Energy)-1,len(Eval_Energy),1))
	poles = numpy.equal(X,numpy.tile(Eval_Energy[numpy.newaxis,:,numpy.newaxis],(len(Energy),1,len(orders))))

	# all N, ln(x+E) and ln(x-E) terms and poles
	Integral = numpy.sum(-C*(-E)**N*numpy.log(numpy.absolute((X[1:,:,:]+E)/(X[:-1,:,:]+E)))-C*E**N*(1-poles[1:,:,:])*numpy.log(numpy.absolute((X[1:,:,:]-E+poles[1:,:,:])/((1-poles[:-1,:,:])*X[:-1,:,:]+poles[:-1,:,:]*X[[0]+list(range(len(Energy)-2)),:,:]-E))),axis=(0,2))

	if numpy.any(orders<=-2): # N<=-2, ln(x) terms
		i = [slice(None,None,None),slice(None,None,None),orders<=-2]
		Integral += numpy.sum(C[i]*((-E[i])**N[i]+E[i]**N[i])*numpy.log(numpy.absolute((X[1:,:,orders<=-2])/(X[:-1,:,orders<=-2]))),axis=(0,2))

	if numpy.any(orders>=0): # N>=0,  x^k terms
		for ni in numpy.where(orders>=0)[0]:
			i = [slice(None,None,None),slice(None,None,None),ni]
			n = orders[ni]
			for k in range(n,0,-2):
				Integral += numpy.sum(C[i]/float(-k)*2*E[i]**(n-k)*(X[1:,:,ni]**k-X[:-1,:,ni]**k),axis=0)

	if numpy.any(orders <=-3): # N<=-3, x^k terms
		for ni in numpy.where(orders<=-3)[0]:
			i = [slice(None,None,None),slice(None,None,None),ni]
			n = orders[ni]
			for k in range(n+2,0,2):
				Integral += numpy.sum(C[i]/float(k)*((-1)**(n-k)+1)*E[i]**(n-k)*(X[1:,:,ni]**k-X[:-1,:,ni]**k),axis=0)

	logger.debug("Done!")
	return Integral / math.pi + relativistic_correction

def KK_PP(Eval_Energy, Energy, imaginary_spectrum, relativistic_correction):
	"""Calculate Kramers-Kronig transform with "Piecewise Polynomial"
	algorithm plus the Biggs and Lighthill extended data.

	Parameters
	----------
	Eval_Energy : numpy vector of `float`
		Set of photon energies describing points at which to evaluate the real spectrum
	Energy : numpy vector of `float`
		Set of photon energies describing intervals for which each row of `imaginary_spectrum` is valid
	imaginary_spectrum : two-dimensional `numpy.array` of `float`
		The array consists of five columns of polynomial coefficients: A_1, A_0, A_-1, A_-2, A_-3
	relativistic_correction : float
		The relativistic correction to the Kramers-Kronig transform.
		You can calculate the value using the `calc_relativistic_correction` function.

	Returns
	-------
	This function returns the real part of the scattering factors evaluated at photon energies specified by Eval_Energy.

	"""
	logger = logging.getLogger(__name__)
	logger.info("Calculate Kramers-Kronig transform using (n from 1 to -3) piecewise-polynomial algorithm")
	X1 = Energy[0:-1]
	X2 = Energy[1:]
	E = numpy.tile(Eval_Energy, (len(Energy)-1, 1)).T
	Full_coeffs = imaginary_spectrum.T
	Symb_1 = (( Full_coeffs[0, :]*E+Full_coeffs[1, :])*(X2-X1)+0.5*Full_coeffs[0, :]*(X2**2-X1**2)-(Full_coeffs[3, :]/E+Full_coeffs[4, :]*E**-2)*numpy.log(numpy.absolute(X2/X1))+Full_coeffs[4, :]/E*(X2**-1-X1**-1))
	Symb_2 = ((-Full_coeffs[0, :]*E+Full_coeffs[1, :])*(X2-X1)+0.5*Full_coeffs[0, :]*(X2**2-X1**2)+(Full_coeffs[3, :]/E-Full_coeffs[4, :]*E**-2)*numpy.log(numpy.absolute(X2/X1))-Full_coeffs[4, :]/E*(X2**-1-X1**-1))+(Full_coeffs[0, :]*E**2-Full_coeffs[1, :]*E+Full_coeffs[2, :]-Full_coeffs[3, :]*E**-1+Full_coeffs[4, :]*E**-2)*numpy.log(numpy.absolute((X2+E)/(X1+E)))
	Symb_3 = (1-1*((X2==E)|(X1==E)))*(Full_coeffs[0, :]*E**2+Full_coeffs[1, :]*E+Full_coeffs[2, :]+Full_coeffs[3, :]*E**-1+Full_coeffs[4, :]*E**-2)*numpy.log(numpy.absolute((X2-E+1*(X2==E))/(X1-E+1*(X1==E))))
	Symb_B = numpy.sum(Symb_2 - Symb_1 - Symb_3, axis=1)  # Sum areas for approximate integral
	# Patch singularities
	hits = Energy[1:-1]==E[:,0:-1]
	E_hits = numpy.append(numpy.insert(numpy.any(hits, axis=0),[0,0],False),[False,False])
	Eval_hits = numpy.any(hits, axis=1)
	X1 = Energy[E_hits[2:]]
	XE = Energy[E_hits[1:-1]]
	X2 = Energy[E_hits[:-2]]
	C1 = Full_coeffs[:, E_hits[2:-1]]
	C2 = Full_coeffs[:, E_hits[1:-2]]
	Symb_singularities = numpy.zeros(len(Eval_Energy))
	Symb_singularities[Eval_hits] = (C2[0, :]*XE**2+C2[1, :]*XE+C2[2, :]+C2[3, :]*XE**-1+C2[4, :]*XE**-2)*numpy.log(numpy.absolute((X2-XE)/(X1-XE)))
	# Finish things off
	KK_Re = (Symb_B-Symb_singularities) / (math.pi*Eval_Energy) + relativistic_correction
	logger.debug("Done!")
	return KK_Re

def improve_accuracy(Full_E, Real_Spectrum, Imaginary_Spectrum, relativistic_correction, tolerance, recursion=50):
	"""Calculate extra data points so that a linear interpolation is more accurate.

	Parameters
	----------
	Full_E : numpy vector of `float`
		Set of photon energies describing intervals for which each row of `imaginary_spectrum` is valid
	Real_Spectrum : numpy vector of `float`
		The real part of the spectrum corresponding to magnitudes at photon energies in Full_E
	Imaginary_Spectrum : two-dimensional `numpy.array` of `float`
		The array consists of five columns of polynomial coefficients: A_1, A_0, A_-1, A_-2, A_-3
	relativistic_correction : float
		The relativistic correction to the Kramers-Kronig transform.
		(You can calculate the value using the `calc_relativistic_correction` function.)
	tolerance : float
		Level of error in linear extrapolation of data values to be allowed.
	recursion : integer
		Number of times an energy interval can be halved before giving up.

	Returns
	-------
	This function returns a numpy array with three columns respectively representing photon energy, the real spectrum and the imaginary spectrum.
	"""
	logger.debug("Improve data accuracy")
	new_points = numpy.cumsum(numpy.ones((len(Full_E)-2,1),dtype=numpy.int8))+1
	Im_values = data.coeffs_to_ASF(Full_E, numpy.vstack((Imaginary_Spectrum,Imaginary_Spectrum[-1])))
	#plot_Im_values = Im_values
	Re_values = Real_Spectrum
	E_values = Full_E
	temp_Im_spectrum = Imaginary_Spectrum[1:]
	count = 0
	improved = 1
	total_improved_points = 0
	while count<recursion and numpy.sum(improved)>0:
		#get E_midpoints
		midpoints = (E_values[new_points-1]+E_values[new_points])/2.
		#evaluate at new points
		Im_midpoints = data.coeffs_to_ASF(midpoints, temp_Im_spectrum)
		Re_midpoints = KK_PP(midpoints, Full_E, Imaginary_Spectrum, relativistic_correction)
		#evaluate error levels
		Im_error = abs((Im_values[new_points-1]+Im_values[new_points])/2. - Im_midpoints)
		Re_error = abs((Re_values[new_points-1]+Re_values[new_points])/2. - Re_midpoints)
		improved = (Im_error>tolerance) | (Re_error>tolerance)
		logger.debug(str(numpy.sum(improved))+" points (out of "+str(len(improved))+") can be improved in pass number "+str(count+1)+".")
		total_improved_points += numpy.sum(improved)
		#insert new points and values
		Im_values = numpy.insert(Im_values,new_points[improved],Im_midpoints[improved])
		Re_values = numpy.insert(Re_values,new_points[improved],Re_midpoints[improved])
		E_values = numpy.insert(E_values,new_points[improved],midpoints[improved])
		#prepare for next loop
		temp_Im_spectrum =numpy.repeat(temp_Im_spectrum[improved],2,axis=0)
		new_points = numpy.where(numpy.insert(numpy.zeros(Im_values.shape, dtype=numpy.bool),new_points[improved],True))[0]
		new_points = numpy.vstack((new_points, new_points+1)).T.flatten()
		count += 1

	#import matplotlib
	#matplotlib.use('WXAgg')
	#import pylab
	#pylab.figure()
	#pylab.plot(Full_E,plot_Im_values,'ok')
	#pylab.plot(Full_E,Real_Spectrum,'og')
	#pylab.plot(midpoints,Im_midpoints,'+b')
	#pylab.plot(midpoints,Re_midpoints,'+r')
	#pylab.plot(E_values,Im_values,'b-')
	#pylab.plot(E_values,Re_values,'r-')
	#pylab.plot(midpoints,Im_error,'b-')
	#pylab.plot(midpoints,Re_error,'r-')
	#pylab.xscale('log')
	#pylab.show()
	logger.info("Improved data accuracy by inserting "+str(total_improved_points)+" extra points.")
	return numpy.vstack((E_values,Re_values,Im_values)).T

def kk_calculate_real(NearEdgeDataFile, ChemicalFormula, load_options=None, input_data_type=None, merge_points=None, add_background=False, fix_distortions=False, curve_tolerance=None, curve_recursion=50):
	"""Do all data loading and processing and then calculate the kramers-Kronig transform.
	Parameters
	----------
	NearEdgeDataFile : string
		Path to file containg near-edge data
	ChemicalFormula : string
		A standard chemical formula string consisting of element symbols, numbers and parentheses.
	merge_points : list or tuple pair of `float` values, or None
		The photon energy values (low, high) at which the near-edge and scattering factor data values
		are set equal so as to ensure continuity of the merged data set.

	Returns
	-------
	This function returns a numpy array with columns consisting of the photon energy, the real and the imaginary parts of the scattering factors.
	"""
	Stoichiometry = data.ParseChemicalFormula(ChemicalFormula)
	Relativistic_Correction = calc_relativistic_correction(Stoichiometry)
	Full_E, Imaginary_Spectrum = data.calculate_asf(Stoichiometry)
	if NearEdgeDataFile is not None:
		NearEdge_Data = data.convert_data(data.load_data(NearEdgeDataFile, load_options),FromType=input_data_type,ToType='asf')
		Full_E, Imaginary_Spectrum = data.merge_spectra(NearEdge_Data, Full_E, Imaginary_Spectrum, merge_points=merge_points, add_background=add_background, fix_distortions=fix_distortions)
	Real_Spectrum = KK_PP(Full_E, Full_E, Imaginary_Spectrum, Relativistic_Correction)
	if curve_tolerance is not None:
		output_data = improve_accuracy(Full_E,Real_Spectrum,Imaginary_Spectrum, Relativistic_Correction, curve_tolerance, curve_recursion)
	else:
		Imaginary_Spectrum_Values = data.coeffs_to_ASF(Full_E, numpy.vstack((Imaginary_Spectrum,Imaginary_Spectrum[-1])))
		output_data = numpy.vstack((Full_E,Real_Spectrum,Imaginary_Spectrum_Values)).T
	return output_data

if __name__ == '__main__':
	#use argparse here to get command line arguments
	#process arguments and pass to a pythonic function

	#I will abuse this section of code for initial testing
	#Output = kk_calculate_real('../../data/Xy_norm_bgsub.txt', 'C10SH14', input_data_type='NEXAFS')
	Output = kk_calculate_real('../../data/LaAlO3/LaAlO3_Exp.csv', 'LaAlO3', input_data_type='NEXAFS', fix_distortions=True, curve_tolerance=0.05)
	#Output = kk_calculate_real('../../data/GaAs/As.xmu.csv', 'GaAs', input_data_type='NEXAFS', fix_distortions=True, curve_tolerance=0.05)

	Stoichiometry = data.ParseChemicalFormula('LaAlO3')
	#Stoichiometry = data.ParseChemicalFormula('GaAs')
	Relativistic_Correction = calc_relativistic_correction(Stoichiometry)
	ASF_E, ASF_Data = data.calculate_asf(Stoichiometry)
	ASF_Data3 = data.coeffs_to_linear(ASF_E, ASF_Data, 0.1)
	ASF_Data2 = data.coeffs_to_ASF(ASF_E, numpy.vstack((ASF_Data,ASF_Data[-1])))

	#Test_E = (Output[1:,0]+Output[0:-1,0])*0.5
	#Test_E = numpy.linspace(41257.87,41259.87,num=21)
	#Real_Spectrum2 = KK_PP(Test_E, Output[:,0], Im, Relativistic_Correction)

	import matplotlib
	matplotlib.use('WXAgg')
	import pylab

	pylab.figure()
	pylab.plot(Output[:,0],Output[:,1],'xg-',Output[:,0],Output[:,2],'xb-')
	pylab.plot(ASF_E,ASF_Data2,'+r')
	#pylab.plot(ASF_E,ASF_Data22,'xr')
	pylab.plot(ASF_Data3[0],ASF_Data3[1],'r-')
	#pylab.plot(Test_E,Real_Spectrum2,'*y')
	pylab.xscale('log')
	pylab.show()

max_stars_repo_name: benajamin/kkcalc

max_issues_repo_name: benajamin/kkcalc

max_forks_repo_name: benajamin/kkcalc


text:
import numpy as np
import matplotlib.pyplot as plt
from utils.show_array import show_array

def specshow(x, spacing=256, sr=44100, max_frames=1960, skip=1, gamma=6, use_mag=True, cmap='inferno'):
    frames, bins = x.shape
    seconds = (frames * spacing) / sr
    minutes = int(seconds // 60)
    seconds = seconds - (minutes * 60)
    sx = x[:max_frames*skip:skip]
    if not np.iscomplexobj(sx):
        spec = np.copy(sx)
    else:
        mag = np.abs(sx)
        spec = mag if use_mag else np.log(1+mag**2)
    spec -= spec.min()
    spec = spec.T / spec.max()
    spec **= 1 / gamma
    cm = plt.get_cmap(cmap)
    spec = cm(spec)[:,:,:3]
    show_array(255 * np.flipud(spec), retina=True)
    print(f'{minutes}:{seconds:04.2f} @ {sr}Hz, {frames} frames x {bins} bins @ {spacing} spacing')

max_stars_repo_name: CreativeInquiry/TeenieHarrisProject

max_issues_repo_name: CreativeInquiry/TeenieHarrisProject

max_forks_repo_name: CreativeInquiry/TeenieHarrisProject


text:
#!/usr/bin/python3
#   produces a bar graph chart for sync failure tests
#   uses graphs described by variable GRAPH
#   uses database described by variable DATABASE
#
# Usage
# python printFailureTest.py #normprob #max_iter #num_trial
# normprob = normalized probability
# #max iter = maximum allowed iteration for LP to converge (typically 1000)
# num_trial = will try for num_trial time and estimate failure rate for each (typically 100)
# example: python3  printFailureTest.py 10 1000 10

import sys
import numpy as np
import pandas as pd
import subprocess
import matplotlib
matplotlib.use('Agg')
matplotlib.rcParams['backend'] = "Qt4Agg"
import matplotlib.pyplot as plt
import os
import sqlite3
from matplotlib import cm

# assumes the data is saved in the database

# modify following list of graph to choose which ones do you want
GRAPHS = ["kron_g500-simple-logn18", "rgg_n_2_18_s0",\
     "astro-ph","cond-mat","cond-mat-2005" , "caidaRouterLevel"]

# get the experiment parameters
fault_rate = int(sys.argv[1])
max_iter = int(sys.argv[2])
num_trials = int(sys.argv[3])

#opening the  data base
DATABASE = 'failureExp.db'
fdb = sqlite3.connect(DATABASE)

# get the data from sqlite database
sqlCmd= "select * from fTable where graphName in ({seq}) \
            and algmType=? and normFaultRate=? and numTrial=? and maxIter=? \
            ".format(seq=','.join(['?']*len(GRAPHS)))
sqlArg=  GRAPHS + ["sync", fault_rate, num_trials, max_iter]

print(sqlCmd);
print(sqlArg);
qResults = fdb.execute(sqlCmd, sqlArg);

qResTuple = qResults.fetchall()

graphName = [ row[1] for row in qResTuple ];
FiSV =      [ int(row[6]) for row in qResTuple ];
SsSV =      [ int(row[7]) for row in qResTuple ];
SshSV =     [ int(row[8]) for row in qResTuple ];
TmrSV =     [ int(row[9]) for row in qResTuple ];

N = len(graphName);
ind = 1+np.arange(N)  # the x locations for the groups
width = 0.1       # the width of the bars

plt.style.use('ggplot')
# plt.style.use('fivethirtyeight')
# plt.style.use('bmh')
fig, ax = plt.subplots()

color_list= ['#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00']
color_list=['#66c2a5','#fc8d62','#8da0cb','#e78ac3','#a6d854']
color_list = plt.cm.Set3(np.linspace(0, 1, 5))
color_list = plt.cm.Accent(np.linspace(0, 1, 5))
color_list = ['#1b9e77','#d95f02','#7570b3','#e7298a']
rects1 = ax.bar(-0.5+ind, 100-np.array(FiSV), width, color=color_list[0], edgecolor='white')
rects2 = ax.bar(-0.5+ind+width, 100-np.array(SsSV), width, color=color_list[1],edgecolor='white')
rects3 = ax.bar(-0.5+ind+2*width, 100-np.array(SshSV), width, color=color_list[2],edgecolor='white')
rects4 = ax.bar(-0.5+ind+3*width, 100-np.array(TmrSV), width, color=color_list[3],edgecolor='white')

# setting x and y axis
ax.set_ylim(-5, 110)

#setting X ticks
plt.xticks(1+np.array(list(range(N))))
ax.set_xticklabels(graphName, rotation=60, ha='center')

plt.savefig('output.pdf')
plt.show()

# gname = os.path.basename(fname).split(".")[0]
# print("Running for graphs: " + gname)
# # command
# algo_types = ['FI', 'SS', 'SSH', 'TM']
# fault_rates = range(9, 16)

# npdata = np.zeros((len(fault_rates), len(algo_types) + 1))

# i = 0
# for fault_rate in fault_rates:
#     cmd = "PRINT_GRAPH=0 MAX_ITER=%s NUM_TRIAL=%d NORM_PROB=%s ../../failureTestSync %s"%(max_iter, num_trials, fault_rate, fname)

#     p = subprocess.Popen(cmd,shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
#     op,err = p.communicate()

#     op = op.decode("utf-8")
#     op = op.split("\n")
#     op = [o.strip() for o in op]
#     print op
#     op2 = op[1].split()

#     npdata[i,0] = fault_rate
#     npdata[i,-1] = int(op2[-1])
#     npdata[i,-2] = int(op2[-2])
#     npdata[i,-3] = int(op2[-3])
#     npdata[i,-4] = int(op2[-4])

#     i+= 1

# output = pd.DataFrame(npdata)

# output.columns = ['fault_rate'] + algo_types
# output.to_csv(gname+'_fault_rate.csv')

# fig, ax = plt.subplots()
# ind = np.arange(len(fault_rates))
# width = 0.2

# rects1 = ax.bar(ind, 100*output['FI']/num_trials, width, color='r')
# rects2 = ax.bar([i + width for i in ind], 100*output['SS']/num_trials, width, color='g')
# rects3 = ax.bar([i + 2*width for i in ind], 100*output['SSH']/num_trials, width, color='b')
# rects4 = ax.bar([i + 3*width for i in ind], 100*output['TM']/num_trials, width, color='purple')

# ax.set_ylabel('Failures %')
# ax.set_ylim(0, 110)
# ax.set_title(gname + 'Failure Rate')
# ax.set_xticks(ind + 2*width)
# ax.set_xticklabels(['2^-'+str(fr) for fr in fault_rates])

# ax.legend((rects1[0], rects2[0], rects3[0], rects4[0]), ('FI', 'SS', 'SSH', 'TM'))

# plt.savefig('')

max_stars_repo_name: ogreen/approxGraphComp

max_issues_repo_name: ogreen/approxGraphComp

max_forks_repo_name: ogreen/approxGraphComp


text:
''' Eingenvector for an eingenvalue of 1

Author: Andrei-Claudiu Roibu, 2019

This code has been created to compare the speed difference between a numpy dot product and a slow dot product. This code was written in support of my learning.

The original code comes from these two sources:

    # https://deeplearningcourses.com/c/deep-learning-prerequisites-the-numpy-stack-in-python
    # https://www.udemy.com/deep-learning-prerequisites-the-numpy-stack-in-python

Description:

Defines a matrix A and a vector V which is uniform and sums to one.

This code multiplies A and V 25 times, obtaining a new v every time.

The code then plots the euclidian distance between each new v, and the original v, ie. |v'-v| . This converges to 0.

This means that if |v'-v| = 0, we have found the eigenvector of A, for which the corresponding eingenvalue is 1.
'''

import numpy as np
import matplotlib.pyplot as plt

def data_generator():
    '''
    This function generates the A and v matrices.

    Args:
        none
    Returns:
        A (array): Starting Matrix
        v (vector): Original value of vector
    '''

    A = np.array(
        [[0.3, 0.6, 0.1],
        [0.5, 0.2, 0.3],
        [0.4, 0.1, 0.5]]
    )

    v = np.ones(3) / 3.0

    return A,v

def eigenvectors(iterations=25):
    '''
    This function calculates the eingenvector

    Args:
        iterations (int): The number of iterations for running the function

    Returns:
        distances (list): A list containing the different euclidian distances
    '''

    A, v = data_generator()

    distances = []

    for _ in range(iterations):
        v_new = v.dot(A)
        distance = np.linalg.norm(v_new - v)
        distances.append(distance)
        v = v_new

    return distances

if __name__ == '__main__':
    distances = eigenvectors()
    plt.plot(distances)
    plt.savefig('Eigenvalue.png')