text:
# -*- coding: utf-8 -*-
import logging
import numpy as np
import pandas as pd
import geopandas as gpd
from geopy.exc import GeopyError
from time import time, sleep
from tqdm import tqdm

from vespid.data.async_requests import async_request_url_list
from vespid.data import extract_named_value_from_lists_of_dicts, replace_special_solr_characters

LOG_FORMAT = '%(asctime)s: %(levelname)s (%(name)s) - %(message)s'
MAX_CONCURRENT_REQUESTS = 100
REQUEST_INTERVAL_SECONDS = 2

logging.basicConfig(
    format=LOG_FORMAT,
    level=logging.INFO,
    datefmt='%m/%d/%Y %H:%M:%S')

logger = logging.getLogger(__name__)

def get_from_openstreetmap(
    df,
    **kwargs
):
    '''
    Given a list of address strings, iterate through them and find the lat/long,
    using as much address info as is possible.

    Note: this currently only has accuracy at the city level, as OSM results get
    pretty wonky/error-prone when you start querying at the street address level
    (especially since a lot of our addresses are weird, being universities).

    Parameters
    ----------
    df: pandas DataFrame with some location data included as columns.

    street: str. Indicates which column in ``df`` has the full address.

    city: str. Indicates which column in ``df`` has the city name.

    state: str. Indicates which column in ``df`` has the state name, if there is one.

    postalcode: str. Indicates which column in ``df`` has the ZIP/postal code, if there is one.

    country: str. Indicates which column in ``df`` has the country name.

    Returns
    -------
    Copy of ``df`` with columns 'address', 'latitude', and 'longitude'.
    '''

    def _build_urls(
        full_address_strings=None,
        location_data=None,
        **kwargs
    ):
        '''
        Given some address data, build OSM URLs to query for them.

        Parameters
        ----------
        full_address_strings: pandas Series of addresses with a format like "Cleveland, OH, USA".
            These will be used to construct free-text queries. All other parameters will be ignored
            if this isn't None.

        location_data: pandas DataFrame. If not using ``full_address_strings``, this DataFrame should
            contain structured data to be used in an address-segment-specific query (e.g. city=X, country=Y),
            with each unique address segment/entity type (e.g. city or state) contained in separate columns.

        street: str. If using ``location_data``, this is the name of the column in which the
            street address strings can be found. If None, this will not be included in the query.
            Format of this should be "<house number> <street name>".

        city: str. If using ``location_data``, this is the name of the column in which the
            city name strings can be found. If None, this will not be included in the query.

        state: str. If using ``location_data``, this is the name of the column in which the
            state name strings can be found. If None, this will not be included in the query.

        postalcode: str. If using ``location_data``, this is the name of the column in which the
            postal code strings can be found. If None, this will not be included in the query.

        country: str. If using ``location_data``, this is the name of the column in which the
            country name strings can be found. If None, this will not be included in the query.

        Returns
        ----------
        pandas Series of query URL strings.
        '''
        base_url = 'https://nominatim.openstreetmap.org/search?'
        url_suffix = '&format=json&limit=1'
        segments = {**kwargs}

        if full_address_strings is not None:
            queries = 'q=' + replace_special_solr_characters(full_address_strings)
            idx = full_address_strings.index

        elif location_data is not None and len(segments) > 0:
            idx = location_data.index
            for i, (field, value) in enumerate(segments.items()):
                if i == 0:
                    queries = field + "=" + replace_special_solr_characters(location_data[value]).fillna('')

                else:
                    queries += "&" + field + "=" + replace_special_solr_characters(location_data[value]).fillna('')

        else:
            raise ValueError("No data passed for ``location_data`` or no columns identified to use from it.")

        output = base_url + queries + url_suffix

        # Make sure our input data index is maintained
        output.index = idx

        return output

    def _geocode(urls):
        '''
        Handles the geocode query buildings, sending, and error catching.
        '''

        # Allows us to re-try the server if we error out every once in a while
        idx_name = 'index'

        def _process_response(response):
            '''
            Assumes that we are only doing limit = 1 for geocode queries
            '''
            try:
                # No guarantees that the JSON is any good...
                data = response.json()
                if len(data) == 0: # empty list? Will return as all-null row
                    return pd.Series(dtype='object')

                elif response.text != '':
                    output = pd.json_normalize(data).replace({None: np.nan})

                    column_name_dict = {'display_name': 'address', 'lat': 'latitude', 'lon': 'longitude'}
                    output.rename(columns=column_name_dict,
                                inplace=True)

                    # Return only the columns we bothered renaming
                    output = output[column_name_dict.values()].iloc[0]

                else:
                    output = pd.Series(dtype='object')

            except Exception:
                if response.status_code == 200:
                    logger.error(f"Something went wrong with our query using URL {response.url}"
                                f"\nHere's the bad response JSON: \n\n{response.json()}",
                                exc_info=False)
                else:
                    logger.error(f"Something went wrong with our query using URL {response.url}"
                        f"\nHere's the bad response text: \n\n{response.text}",
                                exc_info=False)

                output = pd.Series(dtype='object')

            output['url'] = response.url
            return output

        idx_name = 'index'
        headers = {'User-Agent': 'vespid'}

        output = async_request_url_list(
            urls,
            data_processing_fn=_process_response,
            max_workers=1,
            rate_limit_interval_secs=1,
            return_results_ordered=True,
            flatten_result_dict_include_idx=True,
            flattened_idx_field_name=idx_name,
            index=urls.index,
            use_tqdm=True,
            request_verb_kwargs={'headers': headers}
        )

        output = pd.DataFrame(output).set_index(idx_name)
        output.index = output.index.astype(urls.index)
        output.index.name = ''

        # For some reason, they come back as strings
        if 'latitude' in output.columns and 'longitude' in output.columns:
            output[['latitude', 'longitude']] = output[['latitude', 'longitude']].astype(float)

        return output

    # Build out URLs using structured data (e.g. say 'city=CITY' instead of 'query=CITY,STATE,ETC.')
    # Note that OSM is pretty rigid about its tag/address segment values with this type of query
    # so we will likely have a number of null results we'll need to fill in another way (see below)
    location_dict = {**kwargs}

    urls = _build_urls(
        location_data=df,
        city=location_dict['city'],
        state=location_dict['state'],
        postalcode=location_dict['postalcode'],
        country=location_dict['country']
    )

    results = df.copy()
    location_columns = ['address', 'latitude', 'longitude', 'url']

    results[location_columns] = \
        _geocode(urls)

    # Deal with failed responses
    num_null_results = results['address'].isnull().sum()
    address_separator_string = ', '

    # It all errored out!
    if results.empty:
        logger.warn("No geocoding results found! Retrying with the simplest free-text search: city + country...")
        queries = df[location_dict['city']].fillna('') + address_separator_string +\
            df[location_dict['country']].fillna('')

        urls = _build_urls(full_address_strings=queries)
        results[location_columns] = \
            _geocode(urls)

    # Fill in those that couldn't be processed via structured query
    elif num_null_results > 0:
        logger.info(f"{num_null_results} locations didn't get geocoded, retrying as full-text address search...")
        remaining_index = results[results['address'].isnull()].index

        columns = [
            location_dict['city'],
            location_dict['state'],
            location_dict['postalcode'],
            location_dict['country']
        ]
        for i,column in enumerate(columns):
            if i == 0:
                queries = df.loc[remaining_index, column].fillna('').str.strip()

            else:
                queries += address_separator_string + df.loc[remaining_index, column].fillna('').str.strip()

        urls = _build_urls(full_address_strings=queries)
        results.loc[remaining_index, location_columns] = \
            _geocode(urls)

    num_null_results = results['address'].isnull().sum()
    # We have some results, but not fully filled out
    if num_null_results > 0:
        logger.info(f"{num_null_results} locations still not geocoded, "
        "retrying as full-text address search but with only city + country...")
        remaining_index = results[results['address'].isnull()].index

        columns = [
            location_dict['city'],
            location_dict['country']
        ]
        for i,column in enumerate(columns):
            if i == 0:
                queries = df.loc[remaining_index, column].fillna('').str.strip()

            else:
                queries += address_separator_string + df.loc[remaining_index, column].fillna('').str.strip()

        urls = _build_urls(full_address_strings=queries)
        results.loc[remaining_index, location_columns] = \
            _geocode(urls)

    return results


text:
from galry import *
import numpy as np

x = 10 * np.random.rand(100, 50)
imshow(x)

show()

ext: py

max_stars_repo_name: rossant/galry

max_issues_repo_name: rossant/galry

max_forks_repo_name: rossant/galry

reason: import numpy


text:
#-------------------------------------------------------------------------------
#
#  Spherical Harmonic Coefficients.
#
# Author: Martin Paces <martin.paces@eox.at>
#
#-------------------------------------------------------------------------------
# Copyright (C) 2018 EOX IT Services GmbH
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies of this Software or works derived from this Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#-------------------------------------------------------------------------------
# pylint: disable=too-few-public-methods,abstract-method

from numpy import (
    inf, array, zeros, dot, digitize, argsort, abs as aabs, stack,
)
from ..time_util import decimal_year_to_mjd2000, mjd2000_to_decimal_year

class SHCoefficients(object):
    """ Abstract base class for the spherical harmonic coefficients. """

    def __init__(self, is_internal=True, **kwargs):
        self.is_internal = is_internal
        self.validity = self._get_converted_validity(**kwargs)

    @staticmethod
    def _get_converted_validity(validity_start=None, validity_end=None,
                                to_mjd2000=None, **_):
        """ Get validity range converted to MJD2000 using the optionally
        provided conversion function.
        """

        def _get_converted_value_or_default(value, default, conversion_function):
            if value is None:
                value = default
            if conversion_function is not None:
                value = conversion_function(value)
            return value

        return (
            _get_converted_value_or_default(validity_start, -inf, to_mjd2000),
            _get_converted_value_or_default(validity_end, +inf, to_mjd2000),
        )

    def is_valid(self, time):
        """ Check if the time is within the coefficients validity range. """
        validity_start, validity_end = self.validity
        return (validity_start <= time) & (time <= validity_end)

    @property
    def degree(self):
        """ Get [maximum] model degree. """
        raise NotImplementedError

    @property
    def min_degree(self):
        """ Get minimum model degree.
        Below this degree all model coefficients are zero.
        """
        raise NotImplementedError

    def __call__(self, time, **parameters):
        """ Return the matrix of the full model coefficients. """
        raise NotImplementedError

class CombinedSHCoefficients(SHCoefficients):
    """ Model composed of multiple coefficient sets. """

    def __init__(self, *items):
        if len(items) < 1:
            raise ValueError(
                "The composed model must be composed from at least one "
                "coefficient set."
            )

        item = items[0]
        is_internal = item.is_internal
        validity_start, validity_end = item.validity
        degree = item.degree
        min_degree = item.min_degree

        for item in items[1:]:
            if is_internal != item.is_internal:
                raise ValueError(
                    "Mixing of external and internal coefficient sets!"
                )
            new_start, new_end = item.validity
            validity_start = max(validity_start, new_start)
            validity_end = min(validity_end, new_end)
            degree = max(degree, item.degree)
            min_degree = min(min_degree, item.min_degree)

        SHCoefficients. __init__(
            self, is_internal=is_internal, validity_start=validity_start,
            validity_end=validity_end,
        )

        self._degree = degree
        self._min_degree = min_degree
        self._items = items

    @property
    def degree(self):
        return self._degree

    @property
    def min_degree(self):
        return self._min_degree

    def __call__(self, time, **parameters):
        max_degree = parameters.get("max_degree", -1)
        degree = self.degree if max_degree < 0 else min(self.degree, max_degree)
        coeff_full = zeros((coeff_size(degree), 2))
        for item in self._items:
            item_coeff, item_degree = item(time, **parameters)
            coeff_full[:coeff_size(item_degree), :] += item_coeff
        return coeff_full, degree

class SparseSHCoefficients(SHCoefficients):
    """ Base class for sparse spherical harmonic coefficients. """
    def __init__(self, indices, coefficients, **kwargs):
        SHCoefficients.__init__(self, **kwargs)
        n_idx, m_idx = indices[..., 0], indices[..., 1]
        self._degree = n_idx.max()
        self._min_degree = n_idx.min()
        self._index = stack((
            aabs(m_idx) + (n_idx*(n_idx + 1))//2,
            (m_idx < 0).astype('int'),
            n_idx,
        ), 1)
        self._coeff = coefficients

    @property
    def degree(self):
        return self._degree

    @property
    def min_degree(self):
        return self._min_degree

    def _subset(self, min_degree, max_degree):
        """ Get subset of the coefficients for the give min. and max. degrees.
        """
        default_max_degree = self._degree
        default_min_degree = self._min_degree

        #degree = self._degree
        index = self._index
        coeff = self._coeff

        if max_degree < 0:
            max_degree = default_max_degree

        if min_degree < 0:
            min_degree = default_min_degree

        if min_degree > default_min_degree or max_degree < default_max_degree:
            idx, = (
                (index[:, 2] <= max_degree) & (index[:, 2] >= min_degree)
            ).nonzero()
            coeff = coeff[idx]
            index = index[idx]
            degree = default_max_degree

            if index.shape[0] > 0:
                degree = index[:, 2].max()
            else:
                degree = 0
        else:
            degree = default_max_degree

        return degree, coeff, index[:, 0], index[:, 1]

class SparseSHCoefficientsConstant(SparseSHCoefficients):
    """ Time invariant sparse spherical harmonic coefficients. """
    def __init__(self, indices, coefficients, **kwargs):
        SparseSHCoefficients.__init__(self, indices, coefficients, **kwargs)

    def __call__(self, time, **parameters):
        degree, coeff, index, kind = self._subset(
            parameters.get("min_degree", -1), parameters.get("max_degree", -1)
        )
        coeff_full = zeros((coeff_size(degree), 2))
        coeff_full[index, kind] = coeff
        return coeff_full, degree

class SparseSHCoefficientsTimeDependent(SparseSHCoefficients):
    """ Time dependent sparse spherical harmonic coefficients
    evaluated by piecewise linear interpolation of a time series of
    coefficients snapshots.
    """
    def __init__(self, indices, coefficients, times, **kwargs):
        order = argsort(times)
        kwargs['validity_start'] = kwargs.get('validity_start', times[order[0]])
        kwargs['validity_end'] = kwargs.get('validity_end', times[order[-1]])
        self._times = _convert(times[order], kwargs.get("to_mjd2000"))
        SparseSHCoefficients.__init__(
            self, indices, coefficients[:, order], **kwargs
        )

    def __call__(self, time, **parameters):
        degree, coeff, index, kind = self._subset(
            parameters.get("min_degree", -1), parameters.get("max_degree", -1)
        )
        coeff_full = zeros((coeff_size(degree), 2))
        coeff_full[index, kind] = self._interpolate_coefficients(time, coeff)
        return coeff_full, degree

    def _interpolate_coefficients(self, time, coeff):
        """ Return interpolated coefficients. """
        idx, basis = self._interpolation_basis(time)
        return dot(coeff[:, idx], basis).reshape(coeff.shape[0])

    def _interpolation_basis(self, time):
        """ Return indices and corresponding non-zero basis functions. """
        times = self._times
        idx1 = digitize([time], times)[0]
        idx1 = min(times.size - 1, max(1, idx1))
        idx0 = idx1 - 1
        alpha = (time -times[idx0])/(times[idx1] - times[idx0])
        return [idx0, idx1], array([1.0 - alpha, alpha]).reshape((2, 1))

class SparseSHCoefficientsTimeDependentDecimalYear(SparseSHCoefficientsTimeDependent):
    """ Time dependent sparse spherical harmonic coefficients
    evaluated by piecewise linear interpolation of a time series of
    coefficients snapshots with time in decimal year interpolated
    in the decimal years time domain.
    """

    def __init__(self, indices, coefficients, times,
                 to_mjd2000=decimal_year_to_mjd2000,
                 to_decimal_year=mjd2000_to_decimal_year,
                 **kwargs):
        SparseSHCoefficientsTimeDependent.__init__(
            self, indices, coefficients, times, **kwargs
        )
        self._to_decimal_year = to_decimal_year
        # Fix the validity range to be in the expected MJD2000.
        self.validity = self._get_converted_validity(
            *self.validity, to_mjd2000=to_mjd2000
        )

    def __call__(self, time, **parameters):
        return SparseSHCoefficientsTimeDependent.__call__(
            self, self._to_decimal_year(time), **parameters
        )

def coeff_size(degree):
    """ Size of the full coefficient array. """
    return ((degree + 2)*(degree + 1))//2

def _convert(value, conversion_function):
    """ Convert value using the optional conversion function. """
    if conversion_function is not None:
        value = conversion_function(value)
    return value

max_stars_repo_name: CHEN-Zhaohui/geoist

max_issues_repo_name: CHEN-Zhaohui/geoist

max_forks_repo_name: CHEN-Zhaohui/geoist

reason: from numpy


text:
# https://machinelearningmastery.com/multi-class-classification-tutorial-keras-deep-learning-library/
# https://matplotlib.org/users/pyplot_tutorial.html
import matplotlib.pyplot as plt
import numpy
import pandas
import operator
from keras import Sequential
from keras.layers import Dense
from keras.utils import np_utils
from sklearn.preprocessing import LabelEncoder

def find_max(x):
    max_index, max_value = max(enumerate(x), key=operator.itemgetter(1))
    return max_index
    # return max(x)

numpy.random.seed(0)

# Читаем ирисы через пандас
dataframe = pandas.read_csv("iris.csv", header=None)
# Конвертируем dataframe в "контейнеры", связаны
data = dataframe.values
print(type(data))
numpy.random.shuffle(data)
print(data)
# Достаем X и Y, связаны
X = data[:,0:4].astype(float)
y = data[:,4]
print(type(X))
# Создаем энкодер для названий, связаны
encoder = LabelEncoder()
# Классификация, связаны
encoder.fit(y)
# Трансформируемся, связаны
encoded_Y = encoder.transform(y)
# У нас есть вектор числе типа [0,0,0,.....,1,1,1,.....,2,2,2,2] нужно его преобразовать [[ 1.  0.  0.],....,
# [ 0.  1.  0.],....., [ 0.  0.  1.]], связаны
dummy_y = np_utils.to_categorical(encoded_Y)

# Создамка тренеровочные данные и тестовые, чтоб все как у людей
persent = 0.8
len_train = int(len(X)*persent)
trainX = X[0:len_train]
trainY = dummy_y[0:len_train]
testX =  X[0:-len_train]
testY = dummy_y[0:-len_train]

# Создаем модель сети, для кого я это пишу вообще?
model = Sequential()
model.add(Dense(6, input_dim=4, activation='relu'))
model.add(Dense(3, activation='softmax'))
model.compile(loss='mse', optimizer='adam', metrics=['accuracy'])
model.fit(trainX, trainY, batch_size=1, epochs=10,verbose=1)
# Делаем предсказания класса на основе testX
prediction = model.predict(testX)
# Вычисляем ошибку
eval = model.evaluate(testX,testY,verbose=1)
print(len(trainX))
print(len(testX))
print(eval)
# Подготовливаем данные для графика
graph_prediction = []
for i in range(len(prediction)):
    graph_prediction.append(find_max(prediction[i]))
# Просто очень нравится list'ы
graph_real = []
for i in range(len(testY)):
    graph_real.append(find_max(testY[i]))

plt.plot(graph_prediction)
plt.plot(graph_real)
plt.ylabel('Class')
plt.xlabel('Number of flower')
plt.title('Loss value = ' + str(eval[0]))
plt.show()


text:
import cv2
import numpy as np

def cartoonize_image(img, ds_factor=4, sketch_mode=False):
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    img_gray = cv2.medianBlur(img_gray, 7)

    edges = cv2.Laplacian(img_gray, cv2.CV_8U, ksize=5)
    ret, mask = cv2.threshold(edges, 100, 255, cv2.THRESH_BINARY_INV)

    if sketch_mode:
        return cv2.cvtColor(mask, cv2.COLOR_GRAY2BGR)
        #img_sketch = cv2.cvtColor(mask, cv2.COLOR_GRAY2BGR)
        #kernel = np.ones((3,3), np.uint8)
        #img_eroded = cv2.erode(img_sketch, kernel, iterations=1)
        #return cv2.medianBlur(img_eroded, 5)

    img_small = cv2.resize(img, None, fx=1.0/ds_factor, fy=1.0/ds_factor, interpolation=cv2.INTER_AREA)

    num_repetitions = 10
    sigma_color = 5
    sigma_space = 7
    size = 5
    for i in xrange(num_repetitions):
        img_small = cv2.bilateralFilter(img_small, size, sigma_color, sigma_space)

    img_output = cv2.resize(img_small, None, fx=ds_factor, fy=ds_factor, interpolation=cv2.INTER_LINEAR)

    dst = np.zeros(img_gray.shape)
    dst = cv2.bitwise_and(img_output, img_output, mask=mask)

    return dst

if __name__=='__main__':
    cap = cv2.VideoCapture(0)

    cur_char = -1
    prev_char = -1

    while True:
        ret, frame = cap.read()
        frame = cv2.resize(frame, None, fx=0.5, fy=0.5, interpolation=cv2.INTER_AREA)

        c = cv2.waitKey(1)
        if c == 27:
            break

        if c > -1 and c != prev_char:
            cur_char = c

        prev_char = c

        if cur_char == ord('s'):
            cv2.imshow('Cartoonize', cartoonize_image(frame, sketch_mode=True))
        elif cur_char == ord('c'):
            cv2.imshow('Cartoonize', cartoonize_image(frame, sketch_mode=False))
        else:
            cv2.imshow('Cartoonize', frame)

    cap.release()
    cv2.destroyAllWindows()

max_stars_repo_path: Module 2/2/cartoonize.py

max_issues_repo_path: Module 2/2/cartoonize.py

max_forks_repo_path: Module 2/2/cartoonize.py


text:
#!/usr/bin/env python3

"""
program shows turtle following sinuswave
"""

__author__ = "F.Feenstra"

## CODE

#IMPORTS
import numpy as np
import sys
import turtle

#CONSTANTS
MARGE = 10
SCREEN = turtle.Screen()
SAMPLE_RATE = int(SCREEN.window_width()/2 -MARGE) - int(MARGE - SCREEN.window_width()/2)
AMPLITUDE = SCREEN.window_height()/4
PHASE = 0
FREQUENCY = 2

#FUNCTIONS

def comp_y(t):
    """compute the y value of the sin wave at the for a specific sample"""
    y_t = AMPLITUDE * np.sin(2 * np.pi * FREQUENCY * (t/SAMPLE_RATE) + PHASE)
    return y_t

def comp_amplitude(x):
    """compute the y value of the sin wave at the for each sample"""
    y = []
    for t in x:
        y.append(comp_y(t))
    return y

def turtle_wave(start, y):
    """Draw a turtle and move it over the sinewave, starting at _start_"""
    dude = turtle.Turtle(shape="turtle", visible=False)
    dude.penup()
    i = start
    for j in range(0, SAMPLE_RATE):
        pos = dude.pos()
        dude.setheading(np.arctan2(y[j]-pos[1], i-pos[0]) * 180 / np.pi)
        dude.goto(i, y[j])
        if i == start:
            dude.pendown()
            dude.showturtle()
        i +=1
    dude.hideturtle()
    return

def turtle_wave_back(end, y):
    """Draw a turtle and move it backwards over the sinewave, from _end_ to _start_"""
    dude = turtle.Turtle(shape="turtle", visible=False)
    dude.penup()
    i = end
    for j in range(SAMPLE_RATE-1, 0, -1):
        pos = dude.pos()
        dude.setheading(np.arctan2(y[j]-pos[1], i-pos[0]) * 180 / np.pi)
        dude.goto(i, y[j])
        if i == end:
            dude.pendown()
            dude.showturtle()
        i -=1
    return

def main():
    """Main function to calculate sine wave and draw turtle"""
    start = int(MARGE - SCREEN.window_width()/2)
    end = int(SCREEN.window_width()/2 - MARGE)
    x = np.arange(start, end) # the points on the x axis for plotting
    y = comp_amplitude(x) # the points on the y axis for plotting
    turtle_wave(start, y)
    turtle_wave_back(end, y)
    return 0

if __name__ == '__main__':
    exitcode = main()
    sys.exit(exitcode)


text:
from mnist import load_mnist
import numpy as np
# import pickle

# def get_data():
#     (x_train, t_train), (x_test, t_test) = load_mnist(normalize=True, flatten=True, one_hot_label=False)
#     return x_test, t_test

# def init_network():
#     with open("sample_weight.pkl", "rb") as f:
#         network = pickle.load(f)
#     return network

# x, t = get_data()
# network = init_network()

(x_train, t_train), (x_test, t_test) = load_mnist(normalize=True, one_hot_label=True)
# print(x_train.shape)
# print(t_train.shape)
train_size = x_train.shape[0]
batch_size = 10
batch_mask = np.random.choice(train_size, batch_size)

# print(batch_mask)

x_batch = x_train[batch_mask]
t_batch = t_train[batch_mask]

max_stars_repo_name: sdbzzhaotq/DeepLearning

max_issues_repo_name: sdbzzhaotq/DeepLearning

max_forks_repo_name: sdbzzhaotq/DeepLearning


text:
import cv2 as cv
import numpy as np
import sys
import os

def height(img):
    return img.shape[0]

def width(img):
    return img.shape[1]

def draw_staff_lines(img_no_staff, staff):
    img_colored_staff_base = cv.cvtColor(img_no_staff, cv.COLOR_GRAY2RGB)
    for i in range(len(staff.bases)):
        for b in staff.bases[i]:
            for line in b-(np.array(range(5)) * (staff.gap[i] + staff.line_thickness[i])):
                for col in range(i * staff.chunk_size, (i+1) * staff.chunk_size):
                    if (col >= width(img_no_staff) or line >= height(img_no_staff)):
                        break
                    img_colored_staff_base[line][col][0] = 0
                    img_colored_staff_base[line][col][1] = 0
                    img_colored_staff_base[line][col][2] = 255

    return img_colored_staff_base

def show_wait_destroy(winname, img, resize=True):
    if (resize):
        img_resized = cv.resize(img.copy(), (720, 480))
    else:
        img_resized = img
    cv.imshow(winname, img_resized)
    cv.moveWindow(winname, 500, 0)
    cv.waitKey(0)
    cv.destroyWindow(winname)

def read_image_from_args():
    img_file = sys.argv[1]
    img = cv.imread(img_file, 0)
    if img is None:
        raise RuntimeError("Image path not found. Did you forget to specify the images/ folder?")
    return img

def read_image_from_file(filename, color = False):
    if color:
        img = cv.imread(filename, cv.IMREAD_COLOR)
    else:
        img = cv.imread(filename, cv.IMREAD_GRAYSCALE)
    if img is None:
        raise RuntimeError("Image path not found. Did you forget to specify the images/ folder?")
    return img

def convert_to_grayscale(img):
    return cv.cvtColor(img.copy(), cv.COLOR_RGB2GRAY)

max_stars_repo_name: MatthewTurney/SheetMusicRecognition

max_issues_repo_name: MatthewTurney/SheetMusicRecognition

max_forks_repo_name: MatthewTurney/SheetMusicRecognition


text:
# coding:utf-8
"""
Create by Wangmeng Song
July 21,2017
"""

import shapefile as sf
from shapely.geometry import Polygon, Point
import os
import inspect
import numpy as np
import json

filedir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) + "/" + "area"
pingchearea = u'pincheArea.dbf'  # 拼车的区域，已更新
westoutside = u'area26.dbf'   # 已变成area26
eastpick = u'area28.dbf'   # area28
westpick = u'area27.dbf'    # area27
# 成都
chengdujieji = u'zhuanchejieji.dbf'
chengdusongji = u'zhuanchesongji.dbf'
# 宜宾
yibing = ''
# 绵阳
mianyang = ''
# 重庆
chongqing = ''
# 西安
xian = ''
# 在东边的区域
ateastarray = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 27]

class SIDE:
    # 不再用此排班范围
    # def chengduArea(self, advanceGetOnTheCar, rmtspID, rmtspLoc, rmtspSeat):
    #     latfilename = filedir + "/" + schedulearea
    #     polys = sf.Reader(latfilename)
    #     polygon = polys.shapes()
    #     shpfilePoints = []
    #     for shape in polygon:
    #         shpfilePoints = shape.points
    #     polygon = Polygon(shpfilePoints)
    #     delindex = []
    #     for i in range(len(rmtspLoc)):
    #         lng = rmtspLoc[i][1]
    #         lat = rmtspLoc[i][0]
    #         point = Point(lng, lat)
    #         if polygon.contains(point):
    #             continue
    #         else:
    #             advanceGetOnTheCar.append(rmtspID[i])   # 去除显示区域外的订单
    #             delindex.append(i)
    #     # 删除提前上车的订单
    #     for deli in reversed(delindex):
    #         del (rmtspID[deli])
    #         del (rmtspLoc[deli])
    #         del (rmtspSeat[deli])
    # 地区东边和西边的代码，地图东为1，西为2,array
    def ateast(self, orderNum, arealoclist):
        sideNo = np.zeros([orderNum], dtype=int)
        for i in range(len(arealoclist)):
            if arealoclist[i] in ateastarray:
                sideNo[i] = 1
            else:
                sideNo[i] = 2
        return sideNo

    def atwest2out(self, westLoc, orderNo):
        westsideNo = np.zeros([orderNo], dtype=int)
        westoutfilename = filedir + "/" + westoutside
        polys = sf.Reader(westoutfilename)
        polygon = polys.shapes()
        shpfilePoints = []
        for shape in polygon:
            shpfilePoints = shape.points
        polygon = Polygon(shpfilePoints)
        for i in range(len(westLoc)):
            lng = westLoc[i][1]
            lat = westLoc[i][0]
            point = Point(lng, lat)
            if polygon.contains(point):
                westsideNo[i] = 1         # 1表示在2环到2.5环之间
            else:
                westsideNo[i] = 2           # 2表示在西边2环内
        return westsideNo

    # 判断是否在西边2环和2.5环
    # def ateast2out(self, eastLoc, orderNo):
    #     eastsideNo = np.zeros([orderNo], dtype=int)
    #     eastoutfilename = filedir + "/" + eastoutside
    #     polys = sf.Reader(eastoutfilename)
    #     polygon = polys.shapes()
    #     shpfilePoints = []
    #     for shape in polygon:
    #         shpfilePoints = shape.points
    #     polygon = Polygon(shpfilePoints)
    #     for i in range(len(eastLoc)):
    #         lng = eastLoc[i][1]
    #         lat = eastLoc[i][0]
    #         point = Point(lng, lat)
    #         if polygon.contains(point):
    #             eastsideNo[i] = 1         # 1表示在2环到2.5环之间
    #         else:
    #             eastsideNo[i] = 2           # 2表示在东边2环内
    #     return eastsideNo

    def orderinchengdutwofive(self, orderdata):
        latfilename = filedir + "/" + pingchearea
        polys = sf.Reader(latfilename)
        polygon = polys.shapes()
        shpfilePoints = []
        for shape in polygon:
            shpfilePoints = shape.points
        polygon = Polygon(shpfilePoints)
        lng = orderdata['bdlng']
        lat = orderdata['bdlat']
        point = Point(lng, lat)
        if polygon.contains(point):
            orderdata['inside'] = True
        else:
            orderdata['inside'] = False
        jsondatar = json.dumps(orderdata, ensure_ascii=False, separators=(',', ':')).encode('utf-8')
        return jsondatar

    # 判断是否在同一个区域里
    def judgeinarea(self, filepath, lat, lng):
        polys = sf.Reader(filepath)
        polygon = polys.shapes()
        shpfilePoints = []
        for shape in polygon:
            shpfilePoints = shape.points
        polygon = Polygon(shpfilePoints)
        point = Point(lng, lat)
        if polygon.contains(point):
            return True
        else:
            return False

    def specificitywholeChengDu(self, SpecificityOrderdata):
        try:
            city = SpecificityOrderdata['city']
        except:
            # 没有传城市编码，默认为成都
            city = 510100
        tripnum = SpecificityOrderdata['triptype']  # 1为接机，2为送机
        lng = SpecificityOrderdata['bdlng']
        lat = SpecificityOrderdata['bdlat']
        # 宜宾
        if city == 'YBP':
            # 接机和送机类型
            SpecificityOrderdata['inside'] = False
            # latfilename = filedir + "/" + yibing
            # if self.judgeinarea(latfilename, lat, lng):
            #     SpecificityOrderdata['inside'] = True
            # else:
            #     SpecificityOrderdata['inside'] = False
        # 绵阳
        elif city == 510700:
            SpecificityOrderdata['inside'] = True
            # latfilename = filedir + "/" + mianyang
            # if self.judgeinarea(latfilename, lat, lng):
            #     SpecificityOrderdata['inside'] = True
            # else:
            #     SpecificityOrderdata['inside'] = False
        # 重庆
        elif city == 'CKG':
            SpecificityOrderdata['inside'] = False
            # latfilename = filedir + "/" + chongqing
            # if self.judgeinarea(latfilename, lat, lng):
            #     SpecificityOrderdata['inside'] = True
            # else:
            #     SpecificityOrderdata['inside'] = False
        # 西安
        elif city == 610100:
            SpecificityOrderdata['inside'] = True
            # latfilename = filedir + "/" + xian
            # if self.judgeinarea(latfilename, lat, lng):
            #     SpecificityOrderdata['inside'] = True
            # else:
            #     SpecificityOrderdata['inside'] = False
        # 成都
        elif city == 510100:
            if tripnum == 1:  # 1为接机
                latfilename = filedir + "/" + chengdujieji
            else:   # 2为送机
                latfilename = filedir + "/" + chengdusongji
            if self.judgeinarea(latfilename, lat, lng):
                SpecificityOrderdata['inside'] = True
            else:
                SpecificityOrderdata['inside'] = False
        jsondatar = json.dumps(SpecificityOrderdata, ensure_ascii=False, separators=(',', ':')).encode('utf-8')
        return jsondatar

    def eastpick(self, pickpoint):
        eastfile = filedir + "/" + eastpick
        polys = sf.Reader(eastfile)
        polygon = polys.shapes()
        shpfilePoints = []
        for shape in polygon:
            shpfilePoints = shape.points
        polygon = Polygon(shpfilePoints)
        lng = pickpoint[1]
        lat = pickpoint[0]
        point = Point(lng, lat)
        if polygon.contains(point):
            return True
        else:
            return False

    def westpick(self, pickpoint):
        westfile = filedir + "/" + westpick
        polys = sf.Reader(westfile)
        polygon = polys.shapes()
        shpfilePoints = []
        for shape in polygon:
            shpfilePoints = shape.points
        polygon = Polygon(shpfilePoints)
        lng = pickpoint[1]
        lat = pickpoint[0]
        point = Point(lng, lat)
        if polygon.contains(point):
            return True
        else:
            return False


text:
import numpy as np
import cv2
import pickle

from helpers import Helpers
from digit import Digit

class Cells(object):
    '''
    Extracts each cell from the sudoku grid obtained
    from the Extractor
    '''

    def __init__(self, sudoku):
        print 'Extracting cells...',
        self.helpers = Helpers()
        self.cells = self.extractCells(sudoku)
        print 'done.'

    def extractCells(self, sudoku):
        cells = []
        W, H = sudoku.shape
        cell_size = W / 9
        i, j = 0, 0
        for r in range(0, W, cell_size):
            row = []
            j = 0
            for c in range(0, W, cell_size):
                cell = sudoku[r:r + cell_size, c:c + cell_size]
                cell = self.helpers.make_it_square(cell, 28)
                #self.helpers.show(cell, 'Before clean')
                cell = self.clean(cell)
                digit = Digit(cell).digit
                #self.helpers.show(digit, 'After clean')
                digit = self.centerDigit(digit)
                #self.helpers.show(digit, 'After centering')
                row.append(digit // 255)
                j += 1
            cells.append(row)
            i += 1
        return cells

    def clean(self, cell):
        contour = self.helpers.largestContour(cell.copy())
        x, y, w, h = cv2.boundingRect(contour)
        cell = self.helpers.make_it_square(cell[y:y + h, x:x + w], 28)
        kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (2, 2))
        cell = cv2.morphologyEx(cell, cv2.MORPH_CLOSE, kernel)
        cell = 255 * (cell / 130)
        return cell

    def centerDigit(self, digit):
        digit = self.centerX(digit)
        digit = self.centerY(digit)
        return digit

    def centerX(self, digit):
        topLine = self.helpers.getTopLine(digit)
        bottomLine = self.helpers.getBottomLine(digit)
        if topLine is None or bottomLine is None:
            return digit
        centerLine = (topLine + bottomLine) >> 1
        imageCenter = digit.shape[0] >> 1
        digit = self.helpers.rowShift(
            digit, start=topLine, end=bottomLine, length=imageCenter - centerLine)
        return digit

    def centerY(self, digit):
        leftLine = self.helpers.getLeftLine(digit)
        rightLine = self.helpers.getRightLine(digit)
        if leftLine is None or rightLine is None:
            return digit
        centerLine = (leftLine + rightLine) >> 1
        imageCenter = digit.shape[1] >> 1
        digit = self.helpers.colShift(
            digit, start=leftLine, end=rightLine, length=imageCenter - centerLine)
        return digit

max_stars_repo_name: sumanthrao/sudo-sudoku

max_issues_repo_name: sumanthrao/SnapSudoku

max_forks_repo_name: sumanthrao/SnapSudoku


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import time
import sympy as sp
import numpy as np
import pandas as pd
from tqdm import tqdm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.collections import EllipseCollection
from matplotlib import patches
from fractions import Fraction

# ------------------------------------------------------------------------------
# logger
# ------------------------------------------------------------------------------

# logging.basicConfig(format='%(levelname)s %(asctime)s @%(name)s # %(message)s',
#                     datefmt='%m/%d/%Y %I:%M:%S %p',
#                     # filename='example.log',
#                     level=logging.DEBUG)

# ------------------------------------------------------------------------------
# utilities
# ------------------------------------------------------------------------------

def mat_show(mat):
    """Graphical visualization of a 2D matrix.
    """
    fig = plt.figure()
    ax = fig.add_subplot(111)
    cax = ax.matshow(to_np_mat(mat), interpolation='nearest')
    fig.colorbar(cax)

def is_symmetric(a, tol=1e-8):
    """Check if matrix is symmetric.

    """
    return np.allclose(a, a.T, atol=tol)

def mat(array):
    """For a given 2D array return a numpy matrix.
    """
    return np.matrix(array)

def vec(vector):
    """Construct a column vector of type numpy matrix.
    """
    return np.matrix(vector).reshape(-1, 1)

def to_np_mat(sympy_mat):
    """Cast sympy Matrix to numpy matrix of float type.

    Parameters
    ----------
    m: sympy 2D matrix

    Returns
    -------
    a numpy asmatrix

    """
    return np.asmatrix(sympy_mat.tolist(), dtype=np.float)

def to_np_array(sympy_mat):
    """Cast sympy Matrix to numpy matrix of float type. Works for N-D matrices as
    compared to to_np_mat().

    Parameters
    ----------
    m: sympy 2D matrix

    Returns
    -------
    a numpy asmatrix

    """
    return np.asarray(sympy_mat.tolist(), dtype=np.float)

def to_np_vec(sympy_vec):
    """Transforms a 1D sympy vector (e.g. 5 x 1) to numpy array (e.g. (5,)).

    Parameters
    ----------
    v: 1D sympy vector

    Returns
    -------
    a 1D numpy array

    """
    return np.asarray(sp.flatten(sympy_vec), dtype=np.float)

def is_pd(A):
    """Check if matrix is positive definite
    """
    return np.all(np.linalg.eigvals(A) > 0)

def nullspace(A, atol=1e-13, rtol=0):
    """Compute an approximate basis for the nullspace of A.

    The algorithm used by this function is based on the singular value
    decomposition of `A`.

    Parameters
    ----------
    A : ndarray
        A should be at most 2-D.  A 1-D array with length k will be treated
        as a 2-D with shape (1, k)
    atol : float
        The absolute tolerance for a zero singular value.  Singular values
        smaller than `atol` are considered to be zero.
    rtol : float
        The relative tolerance.  Singular values less than rtol*smax are
        considered to be zero, where smax is the largest singular value.

    If both `atol` and `rtol` are positive, the combined tolerance is the
    maximum of the two; that is::
        tol = max(atol, rtol * smax)
    Singular values smaller than `tol` are considered to be zero.

    Return value
    ------------
    ns : ndarray
        If `A` is an array with shape (m, k), then `ns` will be an array
        with shape (k, n), where n is the estimated dimension of the
        nullspace of `A`.  The columns of `ns` are a basis for the
        nullspace; each element in numpy.dot(A, ns) will be approximately
        zero.
    """

    A = np.atleast_2d(A)
    u, s, vh = np.linalg.svd(A)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T
    return ns

def lrs_inequality_vertex_enumeration(A, b):
    """Find the vertices given an inequality system A * x <= b. This function
    depends on lrs library.

    Parameters
    ----------

    A: numpy array [m x n]

    b: numpy array [m]

    Returns
    -------

    v: numpy array [k x n]
        the vertices of the polytope

    """
    # export H-representation
    with open('temp.ine', 'w') as file_handle:
        file_handle.write('Feasible_Set\n')
        file_handle.write('H-representation\n')
        file_handle.write('begin\n')
        file_handle.write(str(A.shape[0]) + ' ' +
                          str(A.shape[1] + 1) + ' rational\n')
        for i in range(0, A.shape[0]):
            file_handle.write(str(Fraction(b[i])))
            for j in range(0, A.shape[1]):
                file_handle.write(' ' + str(Fraction(-A[i, j])))

            file_handle.write('\n')

        file_handle.write('end\n')

    # call lrs
    try:
        os.system('lrs temp.ine > temp.ext')
    except OSError as e:
        raise RuntimeError(e)

    # read the V-representation
    vertices = []
    with open('temp.ext', 'r') as file_handle:
        begin = False
        for line in file_handle:
            if begin:
                if 'end' in line:
                    break

                comp = line.split()
                try:
                    v_type = comp.pop(0)
                except:
                    print('No feasible solution')

                if v_type is '1':
                    v = [float(Fraction(i)) for i in comp]
                    vertices.append(v)

            else:
                if 'begin' in line:
                    begin = True

    # delete temporary files
    try:
        os.system('rm temp.ine temp.ext')
    except OSError as e:
        pass

    return vertices

def ccd_inequality_vertex_enumeration(A, b):
    """Find the vertices given an inequality system A * x <= b. This function
    depends on pycddlib (cdd).

    Parameters
    ----------

    A: numpy array [m x n]

    b: numpy array [m]

    Returns
    -------

    v: numpy array [k x n]
        the vertices of the polytope

    """
    import cdd
    # try floating point, if problem fails try exact arithmetics (slow)
    try:
        M = cdd.Matrix(np.hstack((b.reshape(-1, 1), -A)),
                       number_type='float')
        M.rep_type = cdd.RepType.INEQUALITY
        p = cdd.Polyhedron(M)
    except:
        print('Warning: switch to exact arithmetics')
        M = cdd.Matrix(np.hstack((b.reshape(-1, 1), -A)),
                       number_type='fraction')
        M.rep_type = cdd.RepType.INEQUALITY
        p = cdd.Polyhedron(M)

    G = np.array(p.get_generators())

    if not G.shape[0] == 0:
        return G[np.where(G[:, 0] == 1.0)[0], 1:].tolist()
    else:
        raise ValueError('Infeasible Inequality')

def optimization_based_sampling(A, b, optimziation_samples,
                                closiness_tolerance, max_opt_iterations):
    """Efficient method for sampling the feasible set for a large system of
    inequalities (A x <= b). When the dimension of the set (x) is large,
    deterministic and pure randomized techniques fail to solve this problem
    efficiently.

    This method uses constrained optimization in order to find n solutions that
    satisfy the inequality. Each iteration new randomized objective function is
    assigned so that the optimization will find a different solution.

    Parameters
    ----------
    A: numpy array [m x n]

    b: numpy array [m]

    optimziation_samples: integer
        number of samples to generate

    closiness_tolerance: float
        accept solution if distance is larger than the provided tolerance

    max_opt_iterations: integer
        maximum iteration of the optimization algorithm

    Returns
    -------

    solutions: list

    """
    from scipy.optimize import minimize

    nullity = A.shape[1]
    solutions = []
    j = 0
    while j < optimziation_samples:
        # change objective function randomly (Dirichlet distribution ensures
        # that w sums to 1)
        w = np.random.dirichlet(np.ones(nullity), size=1)

        def objective(x):
            return np.sum((w * x) ** 2)

        # solve the optimization_based_sampling
        def inequality_constraint(x):
            return A.dot(x) - b

        x0 = np.random.uniform(-1, 1, nullity)
        bounds = tuple([(None, None) for i in range(0, nullity)])
        constraints = ({'type': 'ineq', 'fun': inequality_constraint})
        options = {'maxiter': max_opt_iterations}
        sol = minimize(objective, x0, method='SLSQP',
                       bounds=bounds,
                       constraints=constraints,
                       options=options)

        x = sol.x
        # check if solution satisfies the system
        if np.all(A.dot(x) <= b):
            # if solution is not close to the rest then accept
            close = False
            for xs in solutions:
                if np.linalg.norm(xs - x, 2) < closiness_tolerance:
                    close = True
                    break

            if not close:
                solutions.append(x)
                j = j + 1

    return solutions

def convex_bounded_vertex_enumeration(A,
                                      b,
                                      convex_combination_passes=1,
                                      method='lrs'):
    """Sample a convex, bounded inequality system A * x <= b. The vertices of the
    convex polytope are first determined. Then the convexity property is used to
    generate additional solutions within the polytope.

    Parameters
    ----------

    A: numpy array [m x n]

    b: numpy array [m]

    convex_combination_passes: int (default 1)
        recombine vertices to generate additional solutions using the convex
        property

    method: str (lrs or cdd or rnd)

    Returns
    -------

    v: numpy array [k x n]
        solutions within the convex polytope

    """
    # find polytope vertices
    if method == 'lrs':
        solutions = lrs_inequality_vertex_enumeration(A, b)
    elif method == 'cdd':
        solutions = ccd_inequality_vertex_enumeration(A, b)
    elif method == 'rnd':
        solutions = optimization_based_sampling(A, b,
                                                A.shape[0] ** 2,
                                                1e-3, 1000)
    else:
        raise RuntimeError('Unsupported method: choose "lrs" or "cdd" or "rnd"')

    # since the feasible space is a convex set we can find additional solution
    # in the form z = a * x_i + (1-a) x_j
    for g in range(0, convex_combination_passes):
        n = len(solutions)
        for i in range(0, n):
            for j in range(0, n):
                if i == j:
                    continue

                a = 0.5
                x1 = np.array(solutions[i])
                x2 = np.array(solutions[j])
                z = a * x1 + (1 - a) * x2
                solutions.append(z.tolist())

    # remove duplicates from 2D list
    solutions = [list(t) for t in set(tuple(element) for element in solutions)]

    return np.array(solutions, np.float)

def test_inequality_sampling(d):
    A = np.array([[0, 0, -1],
                  [0, -1, 0],
                  [1, 0, 0],
                  [-1, 0, 0],
                  [0, 1, 0],
                  [0, 0, 1]])
    b = 0.5 * np.ones(6)
    print(A, b)
    t1 = time.time()
    # solutions = optimization_based_sampling(A, b, 20, 0.3, 1000)
    solutions = convex_bounded_vertex_enumeration(A, b, d)
    t2 = time.time()
    print('Execution time: ' + str(t2 - t1) + 's')

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(solutions[:, 0], solutions[:, 1], solutions[:, 2])
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    plt.show()
    fig.tight_layout()
    fig.savefig('results/inequality_sampling/inequality_sampling_d' + str(d) + '.png',
                format='png', dpi=300)
    fig.savefig('results/inequality_sampling/inequality_sampling_d' + str(d) + '.pdf',
                format='pdf', dpi=300)
    fig.savefig('results/inequality_sampling/inequality_sampling_d' + str(d) + '.eps',
                format='eps', dpi=300)

# test_inequality_sampling(0)
# test_inequality_sampling(1)
# test_inequality_sampling(2)

def tensor3_vector_product(T, v):
    """Implements a product of a rank-3 tensor (3D array) with a vector using
    tensor product and tensor contraction.

    Parameters
    ----------

    T: sp.Array of dimensions n x m x k

    v: sp.Array of dimensions k x 1

    Returns
    -------

    A: sp.Array of dimensions n x m

    Example
    -------

    >>>T = sp.Array([[[1, 4, 7, 10], [2, 5, 8, 11], [3, 6, 9, 12]],
                     [[13, 16, 19, 22], [14, 17, 20, 23], [15, 18, 21, 24]]])
    ⎡⎡1  4  7  10⎤  ⎡13  16  19  22⎤⎤
    ⎢⎢           ⎥  ⎢              ⎥⎥
    ⎢⎢2  5  8  11⎥  ⎢14  17  20  23⎥⎥
    ⎢⎢           ⎥  ⎢              ⎥⎥
    ⎣⎣3  6  9  12⎦  ⎣15  18  21  24⎦⎦
    >>>v = sp.Array([1, 2, 3, 4]).reshape(4, 1)
    ⎡1⎤
    ⎢ ⎥
    ⎢2⎥
    ⎢ ⎥
    ⎢3⎥
    ⎢ ⎥
    ⎣4⎦
    >>>tensor3_vector_product(T, v)
    ⎡⎡70⎤  ⎡190⎤⎤
    ⎢⎢  ⎥  ⎢   ⎥⎥
    ⎢⎢80⎥  ⎢200⎥⎥
    ⎢⎢  ⎥  ⎢   ⎥⎥
    ⎣⎣90⎦  ⎣210⎦⎦

    """
    import sympy as sp
    assert(T.rank() == 3)
    # reshape v to ensure 1D vector so that contraction do not contain x 1
    # dimension
    v.reshape(v.shape[0], )
    p = sp.tensorproduct(T, v)
    return sp.tensorcontraction(p, (2, 3))

def test_tensor_product():
    T = sp.Array([[[1, 4, 7, 10], [2, 5, 8, 11], [3, 6, 9, 12]],
                  [[13, 16, 19, 22], [14, 17, 20, 23], [15, 18, 21, 24]]])
    v = sp.Array([1, 2, 3, 4]).reshape(4, 1)
    display(T, v)
    display(tensor3_vector_product(T, v))

# test_tensor_product()

def draw_ellipse(ax, xc, A, scale=1.0, show_axis=False):
    """Construct an ellipse representation of a 2x2 matrix.

    Parameters
    ----------
    ax: plot axis

    xc: np.array 2 x 1
        center of the ellipse
    mat: np.array 2 x 2

    scale: float (default=1.0)
        scale factor of the principle axes

    """
    eigen_values, eigen_vectors = np.linalg.eig(A)
    idx = np.abs(eigen_values).argsort()[::-1]
    eigen_values = eigen_values[idx]
    eigen_vectors = eigen_vectors[:, idx]
    phi = np.rad2deg(np.arctan2(eigen_vectors[1, 0], eigen_vectors[0, 0]))

    ellipse = patches.Ellipse(xy=(xc[0, 0], xc[1, 0]),
                              width=2 * scale * eigen_values[0],
                              height=2 * scale * eigen_values[1],
                              angle=phi,
                              linewidth=2, fill=False)
    ax.add_patch(ellipse)

    # axis
    if show_axis:
        x_axis = np.array([[xc[0, 0], xc[1, 0]],
                           [xc[0, 0] + scale * np.abs(eigen_values[0]) * eigen_vectors[0, 0],
                            xc[1, 0] + scale * np.abs(eigen_values[0]) * eigen_vectors[1, 0]]])
        y_axis = np.array([[xc[0, 0], xc[1, 0]],
                           [xc[0, 0] + scale * eigen_values[1] * eigen_vectors[0, 1],
                            xc[1, 0] + scale * eigen_values[1] * eigen_vectors[1, 1]]])
        ax.plot(x_axis[:, 0], x_axis[:, 1], '-r', label='x-axis')
        ax.plot(y_axis[:, 0], y_axis[:, 1], '-g', label='y-axis')

    return phi, eigen_values, eigen_vectors

def test_ellipse():
    fig, ax = plt.subplots()
    xc = vec([0, 0])
    M = mat([[2, 1], [1, 2]])
    # M = mat([[-2.75032375, -11.82938331], [-11.82938331, -53.5627191]])
    print np.linalg.matrix_rank(M)
    phi, l, v = draw_ellipse(ax, xc, M, 1, True)
    print(phi, l, v)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.axis('equal')
    ax.legend()
    fig.show()

# test_ellipse()

def calculate_feasible_muscle_set(feasible_muscle_set_analysis, base_name,
                                  t_start, t_end, dt, speed):
    """ Calculates the feasible muscle space of a simulation.

    Parameters
    ----------

    feasible_muscle_set_analysis: FeasibleMuscleSetAnalysis

    base_name: base name of simulation files

    t_start: t start

    t_end: t end

    dt: time interval for reporting

    speed: speed of animation

    """
    print('Calculating feasible muscle set ...')
    time = np.linspace(t_start, t_end, t_end / dt + 1, endpoint=True)
    for i, t in enumerate(tqdm(time)):
        visualize_feasible_muscle_set(feasible_muscle_set_analysis, t,
                                      base_name + str(i).zfill(6), 'png')

    command = 'convert -delay ' + \
              str(speed * dt) + ' -loop 0 ' + base_name + \
              '*.png ' + base_name + 'anim.gif'

    print(command)
    try:
        os.system(command)
    except:
        print('unable to execute command')

def visualize_feasible_muscle_set(feasible_muscle_set_analysis, t,
                                  fig_name='fig/feasible_muscle_set', format='png'):
    """ Visualization of the feasible muscle space.

    Parameters
    ----------

    feasible_muscle_set_analysis: FeasibleMuscleSetAnalysis
    t: time instance to evaluate the feasible
    fig_name: figure name for saving
    format: format (e.g. .png, .pdf, .eps)

    """
    fig, ax = plt.subplots(1, 3, figsize=(15, 5))
    feasible_muscle_set_analysis.visualize_simple_muscle(t, ax)
    fig.suptitle('t = ' + str(np.around(t, decimals=2)),
                 y=1.00, fontsize=12, fontweight='bold')
    fig.tight_layout()
    fig.savefig(fig_name + '.' + format, format=format, dpi=300)
    fig.savefig(fig_name + '.pdf', format='pdf', dpi=300)
    fig.savefig(fig_name + '.eps', format='eps', dpi=300)

def apply_generalized_force(f):
    """Applies a generalized force (f) in a manner that is consistent with Newton's
    3rd law.

    Parameters
    ----------
    f: generalized force

    """
    n = len(f)
    tau = []
    for i in range(0, n):
        if i == n - 1:
            tau.append(f[i])
        else:
            tau.append(f[i] - f[i + 1])

    return tau

def custom_exponent(q, A, k, q_lim):
    """ Sympy representation of custom exponent function.

    f(q) = A e^(k (q - q_lim)) / (150) ** k

    """
    return A * sp.exp(k * (q - q_lim)) / (148.42) ** k

def coordinate_limiting_force(q, q_low, q_up, a, b):
    """A continuous coordinate limiting force for a rotational joint.

    It applies an exponential force when approximating a limit. The convention
    is that positive force is generated when approaching the lower limit and
    negative when approaching the upper. For a = 1, F ~= 1N at the limits.

    Parameters
    ----------
    q: generalized coordinate
    q_low: lower limit
    q_up: upper limit
    a: force at limits
    b: rate of the slop

    Note: q, q_low, q_up must have the same units (e.g. rad)

    """
    return custom_exponent(q_low + 5, a, b, q) - custom_exponent(q, a, b, q_up - 5)

def test_limiting_force():
    """
    """
    q = np.linspace(0, np.pi / 4, 100, endpoint=True)
    f = [coordinate_limiting_force(qq, 0, np.pi / 4, 1, 50) for qq in q]

    plt.plot(q, np.array(f))
    plt.show()

def gaussian(x, a, m, s):
    """Gaussian function.

    f(x) = a e^(-(x - m)^2 / (2 s ^2))

    Parameters
    ----------
    x: x
    a: peak
    m: mean
    s: standard deviation

    For a good approximation of an impulse at t = 0.3 [x, 1, 0.3, 0.01].

    """
    return a * np.exp(- (x - m) ** 2 / (2 * s ** 2))

def test_gaussian():
    """
    """
    t = np.linspace(0, 2, 200)
    y = [gaussian(tt, 0.4, 0.4, 0.01) for tt in t]
    plt.plot(t, y)
    plt.show()

def rotate(origin, point, angle):
    """Rotate a point counterclockwise by a given angle around a given origin.

    The angle should be given in radians.

    """
    R = np.asmatrix([[np.cos(angle), - np.sin(angle)],
                     [np.sin(angle), np.cos(angle)]])

    q = origin + R * (point - origin)

    return q

def sigmoid(t, t0, A, B):
    """Implementation of smooth sigmoid function.

    Parameters
    ----------
    t: time to be evalutaed
    t0: delay
    A: magnitude
    B: slope

    Returns
    -------
    (y, y', y'')

    """
    return (A * (np.tanh(B * (t - t0)) + 1) / 2,
            A * B * (- np.tanh(B * (t - t0)) ** 2 + 1) / 2,
            - A * B ** 2 * (- np.tanh(B * (t - t0)) ** 2 + 1) * np.tanh(B * (t - t0)))

def test_sigmoid():
    """
    """
    t, A, B, t0 = sp.symbols('t A B t0')
    y = A / 2 * (sp.tanh(B * (t - t0 - 1)) + 1)
    yd = sp.diff(y, t)
    ydd = sp.diff(yd, t)
    print('\n', y, '\n', yd, '\n', ydd)

    tt = np.linspace(-2, 2, 100)
    yy = np.array([sigmoid(x, 0.5, 2, 5) for x in tt])
    plt.plot(tt, yy)
    plt.show()

def plot_corr_ellipses(data, ax=None, **kwargs):
    """For a given correlation matrix "data", plot the correlation matrix in terms
    of ellipses.

    parameters
    ----------
    data: Pandas dataframe containing the correlation of the data (df.corr())
    ax: axis (e.g. fig, ax = plt.subplots(1, 1))
    kwards: keywords arguments (cmap="Greens")

    https://stackoverflow.com/questions/34556180/
    how-can-i-plot-a-correlation-matrix-as-a-set-of-ellipses-similar-to-the-r-open

    """
    M = np.array(data)
    if not M.ndim == 2:
        raise ValueError('data must be a 2D array')
    if ax is None:
        fig, ax = plt.subplots(1, 1, subplot_kw={'aspect': 'equal'})
        ax.set_xlim(-0.5, M.shape[1] - 0.5)
        ax.set_ylim(-0.5, M.shape[0] - 0.5)

    # xy locations of each ellipse center
    xy = np.indices(M.shape)[::-1].reshape(2, -1).T

    # set the relative sizes of the major/minor axes according to the strength of
    # the positive/negative correlation
    w = np.ones_like(M).ravel()
    h = 1 - np.abs(M).ravel()
    a = 45 * np.sign(M).ravel()

    ec = EllipseCollection(widths=w, heights=h, angles=a, units='x', offsets=xy,
                           transOffset=ax.transData, array=M.ravel(), **kwargs)
    ax.add_collection(ec)

    # if data is a DataFrame, use the row/column names as tick labels
    if isinstance(data, pd.DataFrame):
        ax.set_xticks(np.arange(M.shape[1]))
        ax.set_xticklabels(data.columns, rotation=90)
        ax.set_yticks(np.arange(M.shape[0]))
        ax.set_yticklabels(data.index)

    return ec

def get_cmap(n, name='hsv'):
    """Returns a function that maps each index in 0, 1, ..., n-1 to a distinct RGB
    color; the keyword argument name must be a standard mpl colormap name.

    """
    return plt.cm.get_cmap(name, n)

def assert_if_same(A, B):
    """Assert whether two quantities (value, vector, matrix) are the same."""
    assert np.isclose(
        np.array(A).astype(np.float64),
        np.array(B).astype(np.float64)).all() == True, 'quantities not equal'

def christoffel_symbols(M, q, i, j, k):
    """
    M [n x n]: inertia mass matrix
    q [n x 1]: generalized coordinates
    i, j, k  : the indexies to be computed
    """
    return sp.Rational(1, 2) * (sp.diff(M[i, j], q[k]) + sp.diff(M[i, k], q[j]) -
                                sp.diff(M[k, j], q[i]))

def coriolis_matrix(M, q, dq):
    """
    Coriolis matrix C(q, qdot) [n x n]
    Coriolis forces are computed as  C(q, qdot) *  qdot [n x 1]
    """
    n = M.shape[0]
    C = sp.zeros(n, n)
    for i in range(0, n):
        for j in range(0, n):
            for k in range(0, n):
                C[i, j] = C[i, j] + christoffel_symbols(M, q, i, j, k) * dq[k]
    return C

def substitute(symbols, constants):
    """For a given symbolic sequence substitute symbols."""
    return np.array([substitute(exp, constants)
                     if hasattr(exp, '__iter__')
                     else exp.subs(constants) for exp in symbols])

reason: import numpy,from scipy,import sympy


text:
import cv2
#import cv
import pylab
from SIGBTools import RegionProps
from SIGBTools import getLineCoordinates
from SIGBTools import ROISelector
from SIGBTools import getImageSequence
import numpy as np
import sys
from scipy.cluster.vq import *
from scipy.misc import imresize
from matplotlib.pyplot import *
import scipy

inputFile = "Sequences/eye1.m4v"
outputFile = "eyeTrackerResult.mp4"

#--------------------------
#         Global variable
#--------------------------
global imgOrig,leftTemplate,rightTemplate,frameNr
imgOrig = [];
#These are used for template matching
leftTemplate = []
rightTemplate = []
frameNr =0;
img = []

def GetPupil(gray,thr,minS,maxS):
	'''Given a gray level image, gray and threshold value return a list of pupil locations'''
	tempResultImg = cv2.cvtColor(gray,cv2.COLOR_GRAY2BGR) #used to draw temporary results
	#cv2.circle(tempResultImg,(100,200), 2, (0,0,255),4) #draw a circle
	#cv2.imshow("TempResults",tempResultImg)

	props = RegionProps()
	val,binI =cv2.threshold(gray, thr, 255, cv2.THRESH_BINARY_INV)
	cv2.imshow("Threshold",binI)
	#Calculate blobs
	_, contours, hierarchy= cv2.findContours(binI, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
	pupils = [];
	pupilEllipses = [];
	Centroids = [];
	# YOUR IMPLEMENTATION HERE !!!!

	#CalcContourProperties

	for cnt in contours:
		vals = props.CalcContourProperties(cnt,['Area','Length','Centroid','Extend','ConvexHull'])
		if vals['Area']> minS and vals['Area']< maxS:
			if vals['Extend'] < 1:

				pupils.append(cnt)
				#print vals['Centroid']
				#cv2.circle(tempResultImg,(int(vals['Centroid'][0]),int(vals['Centroid'][1])), 2, (255,0,0),-1) #draw a circle
				Centroids.append(vals['Centroid'])
				pupilEllipse = cv2.fitEllipse(cnt)
				pupilEllipses.append(pupilEllipse)
				#cv2.ellipse(tempResultImg,ellipse,(0,255,0),2)
	#cv2.imshow("TempResults",tempResultImg)
	return (pupils,pupilEllipses,Centroids)

def GetGlints(gray,thr,minS,maxS):
	''' Given a gray level image, gray and threshold
	value return a list of glint locations'''
	# YOUR IMPLEMENTATION HERE !!!!
	tempResultImg = cv2.cvtColor(gray,cv2.COLOR_GRAY2BGR) #used to draw temporary results
	#cv2.circle(tempResultImg,(100,200), 2, (0,0,255),4) #draw a circle
	#cv2.imshow("TempResults",tempResultImg)

	props = RegionProps()
	val,binI =cv2.threshold(gray, thr, 255, cv2.THRESH_BINARY)
	cv2.imshow("ThresholdInverse",binI)
	#Calculate blobs
	_, contours, hierarchy= cv2.findContours(binI, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
	glints = [];
	glintEllipses = [];
	Centroids = [];
	# YOUR IMPLEMENTATION HERE !!!!

	#CalcContourProperties

	for cnt in contours:
		vals = props.CalcContourProperties(cnt,['Area','Length','Centroid','Extend','ConvexHull'])
		if vals['Area']> minS and vals['Area']< maxS:
			if vals['Extend'] < 1:

				glints.append(cnt)
				#print vals['Centroid']
				#cv2.circle(tempResultImg,(int(vals['Centroid'][0]),int(vals['Centroid'][1])), 2, (255,0,0),-1) #draw a circle
				Centroids.append(vals['Centroid'])
				glintEllipse = cv2.fitEllipse(cnt)
				glintEllipses.append(glintEllipse)
				#cv2.ellipse(tempResultImg,ellipse,(0,255,0),2)
	#cv2.imshow("TempResultsglint",tempResultImg)
	return (glints,glintEllipses,Centroids)

def GetIrisUsingThreshold(gray,pupil):
	''' Given a gray level image, gray and threshold
	value return a list of iris locations'''
	# YOUR IMPLEMENTATION HERE !!!!
	pass

def circularHough(gray):
	''' Performs a circular hough transform of the image, gray and shows the  detected circles
	The circe with most votes is shown in red and the rest in green colors '''
 #See help for http://opencv.itseez.com/modules/imgproc/doc/feature_detection.html?highlight=houghcircle#cv2.HoughCircles
	blur = cv2.GaussianBlur(gray, (31,31), 11)

	dp = 6; minDist = 30
	highThr = 20 #High threshold for canny
	accThr = 850; #accumulator threshold for the circle centers at the detection stage. The smaller it is, the more false circles may be detected
	maxRadius = 50;
	minRadius = 155;
	circles = cv2.HoughCircles(blur,cv2.cv.CV_HOUGH_GRADIENT, dp,minDist, None, highThr,accThr,maxRadius, minRadius)

	#Make a color image from gray for display purposes
	gColor = cv2.cvtColor(gray, cv2.COLOR_GRAY2BGR)
	if (circles !=None):
		#print circles
		all_circles = circles[0]
		M,N = all_circles.shape
		k=1
		for c in all_circles:
			cv2.circle(gColor, (int(c[0]),int(c[1])),c[2], (int(k*255/M),k*128,0))
			K=k+1
		c=all_circles[0,:]
		cv2.circle(gColor, (int(c[0]),int(c[1])),c[2], (0,0,255),5)
		cv2.imshow("hough",gColor)

def GetIrisUsingNormals(gray,pupil,normalLength):
	''' Given a gray level image, gray and the length of the normals, normalLength
	 return a list of iris locations'''
	# YOUR IMPLEMENTATION HERE !!!!
	pass

def GetIrisUsingSimplifyedHough(gray,pupil):
	''' Given a gray level image, gray
	return a list of iris locations using a simplified Hough transformation'''
	# YOUR IMPLEMENTATION HERE !!!!
	pass

def GetEyeCorners(leftTemplate, rightTemplate,pupilPosition=None):
	pass

def FilterPupilGlint(pupilcenters,glintcenters,minDisG,maxDisG):
	''' Given a list of pupil candidates and glint candidates returns a list of pupil and glints'''
	for pupilcenter1 in pupilcenters:
		for pupilcenter2 in pupilcenters:
			if minDisG< Distance(pupilcenter1, pupilcenter2) < maxDisG:
				return (pupilcenter1, pupilcenter2)

def Distance(u,v):
	#dst= scipy.distance.euclidean(u,v)
	distance=np.linalg.norm(np.array(u)-np.array(v))
	return distance

def update(I):
	'''Calculate the image features and display the result based on the slider values'''
	#global drawImg
	global frameNr,drawImg
	img = I.copy()
	sliderVals = getSliderVals()
	gray = cv2.cvtColor(img,cv2.COLOR_RGB2GRAY)
	# Do the magic
	#(pupils,pupilEllipses,Centroids) = GetPupil(gray,sliderVals['pupilThr'],sliderVals['minSize'],sliderVals['maxSize'])
	#(glints,glintEllipses,glintCentroids) = GetGlints(gray,sliderVals['glintThr'],sliderVals['minSizeG'],sliderVals['maxSizeG'])
	#FilterPupilGlint(Centroids,glintCentroids,1,100000000)

	#Kmeans implementation
	#detectPupilKMeans(gray,20,2,reSize=(40,40))

	#Hough Transform for pupils
	#detectPupilHough(gray,sliderVals['minDist'],sliderVals['highThr'],sliderVals['accThr'],sliderVals['minRadius'],sliderVals['maxRadius'])

	#Hough Transform for Iris
	detectIrisHough(gray,sliderVals['minDistI'],sliderVals['highThrI'],sliderVals['accThrI'],sliderVals['minRadiusI'],sliderVals['maxRadiusI'])

	#Do template matching
	global leftTemplate
	global rightTemplate
	GetEyeCorners(leftTemplate, rightTemplate)
	#Display results
	global frameNr,drawImg
	x,y = 10,10
	#setText(img,(x,y),"Frame:%d" %frameNr)
	sliderVals = getSliderVals()

	# for non-windows machines we print the values of the threshold in the original image
	if sys.platform != 'win32':
		step=18
		cv2.putText(img, "pupilThr :"+str(sliderVals['pupilThr']), (x, y+step), cv2.FONT_HERSHEY_PLAIN, 1.0, (255, 255, 255), lineType=cv2.LINE_AA)
		cv2.putText(img, "glintThr :"+str(sliderVals['glintThr']), (x, y+2*step), cv2.FONT_HERSHEY_PLAIN, 1.0, (255, 255, 255), lineType=cv2.LINE_AA)
	#cv2.imshow('Result',img)

		#Uncomment these lines as your methods start to work to display the result in the
		#original image
		# for pupil in pupils:
	"""if pupilEllipses is not None and Centroids is not None:

		for ellipse in pupilEllipses:

			cv2.ellipse(img,ellipse,(0,255,0),1)
		for centroid in Centroids:

			C = (int(centroid[0]),int(centroid[1]))
			cv2.circle(img,C, 2, (0,0,255),4)
	for glint in glintCentroids:
		C = int(glint[0]),int(glint[1])
		cv2.circle(img,C, 2,(255,0,255),5)
		#     cv2.imshow("Result", img)

		#For Iris detection - Week 2
		#circularHough(gray)

	"""
	cv2.imshow('Result',img)
	#copy the image so that the result image (img) can be saved in the movie
	drawImg = img.copy()

def printUsage():
	print "Q or ESC: Stop"
	print "SPACE: Pause"
	print "r: reload video"
	print 'm: Mark region when the video has paused'
	print 's: toggle video  writing'
	print 'c: close video sequence'

def run(fileName,resultFile='eyeTrackingResults.avi'):

	''' MAIN Method to load the image sequence and handle user inputs'''
	global imgOrig, frameNr,drawImg
	setupWindowSliders()
	props = RegionProps();
	cap,imgOrig,sequenceOK = getImageSequence(fileName)
	videoWriter = 0

	frameNr =0
	if(sequenceOK):
		update(imgOrig)
	printUsage()
	frameNr=0;
	saveFrames = False

	while(sequenceOK):
		sliderVals = getSliderVals();
		frameNr=frameNr+1
		ch = cv2.waitKey(1)
		#Select regions
		if(ch==ord('m')):
			if(not sliderVals['Running']):
				roiSelect=ROISelector(imgOrig)
				ptsL,regionSelectedL= roiSelect.SelectArea('Select left eye corner',(400,200))
				if(regionSelectedL):
					leftTemplate = imgOrig[ptsL[0][1]:ptsL[1][1],ptsL[0][0]:ptsL[1][0]]

				ptsR,regionSelectedR= roiSelect.SelectArea('Select right eye corner',(600,300))
				if(regionSelectedR):
					rightTemplate = imgOrig[ptsR[0][1]:ptsR[1][1],ptsR[0][0]:ptsR[1][0]]

		if ch == 27:
			break
		if (ch==ord('s')):
			if((saveFrames)):
				videoWriter.release()
				saveFrames=False
				print "End recording"
			else:
				imSize = np.shape(imgOrig)
				videoWriter = cv2.VideoWriter(resultFile, cv.CV_FOURCC('D','I','V','3'), 15.0,(imSize[1],imSize[0]),True) #Make a video writer
				saveFrames = True
				print "Recording..."

		if(ch==ord('q')):
			break
		if(ch==32): #Spacebar
			sliderVals = getSliderVals()
			cv2.setTrackbarPos('Stop/Start','Threshold',not sliderVals['Running'])
		if(ch==ord('r')):
			frameNr =0
			sequenceOK=False
			cap,imgOrig,sequenceOK = getImageSequence(fileName)
			update(imgOrig)
			sequenceOK=True

		sliderVals=getSliderVals()
		if(sliderVals['Running']):
			sequenceOK, imgOrig = cap.read()
			if(sequenceOK): #if there is an image
				update(imgOrig)
			if(saveFrames):
				videoWriter.write(drawImg)
	if(videoWriter!=0):
		videoWriter.release()
        print "Closing videofile..."
#------------------------

""""def detectPupilKMeans(gray,K=2,distanceWeight=2,reSize=(40,40)):
	''' Detects the pupil in the image, gray, using k-means
		gray              : grays scale image
		K                 : Number of clusters
		distanceWeight    : Defines the weight of the position parameters
		reSize            : the size of the image to do k-means on
	'''
	#Resize for faster performance
	smallI = cv2.resize(gray, reSize)
	M,N = smallI.shape
	#Generate coordinates in a matrix
	X,Y = np.meshgrid(range(M),range(N))
	#Make coordinates and intensity into one vectors
	z = smallI.flatten()
	x = X.flatten()
	y = Y.flatten()
	O = len(x)
	#make a feature vectors containing (x,y,intensity)
	features = np.zeros((O,3))
	features[:,0] = z;
	features[:,1] = y/distanceWeight; #Divide so that the distance of position weighs less than intensity
	features[:,2] = x/distanceWeight;
	features = np.array(features,'f')
	# cluster data
	centroids,variance = kmeans(features,K)
	#use the found clusters to map
	label,distance = vq(features,centroids)
	# re-create image from
	labelIm = np.array(np.reshape(label,(M,N)))
	f = figure(1)
	imshow(labelIm)
	f.canvas.draw()
	f.show()
"""
#------------------------------------------------
#   Methods for segmentation
#------------------------------------------------
def detectPupilKMeans(gray,K=2,distanceWeight=2,reSize=(40,40)):
	''' Detects the pupil in the image, gray, using k-means
			gray              : grays scale image
			K                 : Number of clusters
			distanceWeight    : Defines the weight of the position parameters
			reSize            : the size of the image to do k-means on
		'''
	#Resize for faster performance
	smallI = cv2.resize(gray, reSize)
	M,N = smallI.shape
	#Generate coordinates in a matrix
	X,Y = np.meshgrid(range(M),range(N))
	#Make coordinates and intensity into one vectors
	z = smallI.flatten()
	x = X.flatten()
	y = Y.flatten()
	O = len(x)
	#make a feature vectors containing (x,y,intensity)
	features = np.zeros((O,3))
	features[:,0] = z;
	features[:,1] = y/distanceWeight; #Divide so that the distance of position weighs less than intensity
	features[:,2] = x/distanceWeight;
	features = np.array(features,'f')
	# cluster data
	centroids,variance = kmeans(features,K)
	#use the found clusters to map
	label,distance = vq(features,centroids)
	# re-create image from
	labelIm = np.array(np.reshape(label,(M,N)))
	f = figure(1)
	imshow(labelIm)
	f.canvas.draw()
	f.show()

def detectPupilHough(gray,minDist,highThr,accThr,minRadius,maxRadius):
	#Using the Hough transform to detect ellipses
	blur = cv2.GaussianBlur(gray, (9,9),3)
	##Pupil parameters
	dp = 6
	#minDist = 10
	#highThr = 10 #High threshold for canny
	#accThr = 200; #accumulator threshold for the circle centers at the detection stage. The smaller it is, the more false circles may be detected
	#maxRadius = 50;
	#minRadius = 20;
	#See help for http://opencv.itseez.com/modules/imgproc/doc/feature_detection.html?highlight=houghcircle#cv2.HoughCirclesIn thus
	circles = cv2.HoughCircles(blur,cv2.HOUGH_GRADIENT, dp,minDist, None, highThr,accThr,minRadius, maxRadius)
	#Print the circles
	gColor = cv2.cvtColor(gray, cv2.COLOR_GRAY2BGR)
	if (circles !=None):
		#print circles
		all_circles = circles[0]
		M,N = all_circles.shape
		k=1
		for c in all_circles:
			cv2.circle(gColor, (int(c[0]),int(c[1])),c[2], (int(k*255/M),k*128,0))
			#setText(gColor, (int(c[0]),int(c[1])), "Radius:" + str(c[2]))
			cv2.putText(gColor,"Radius:" + str(c[2]), (int(c[0]),int(c[1])), cv2.FONT_HERSHEY_SIMPLEX, .5, 255)
			K=k+1
			#Circle with max votes
		c=all_circles[0,:]
		cv2.circle(gColor, (int(c[0]),int(c[1])),c[2], (0,0,255))
	#text labels

	setText(gColor, (10, 10), "minDist:" + str(minDist) + " highThr:" + str(highThr) +" accThr:" + str(accThr) +" minRadius:" + str(minRadius) +" maxRadius:" + str(maxRadius) )

	cv2.imshow("hough",gColor)

def detectIrisHough(gray,minDist,highThr,accThr,minRadius,maxRadius):
	#Using the Hough transform to detect ellipses
	blur = cv2.GaussianBlur(gray, (9,9),3)
	##Pupil parameters
	dp = 6
	#minDist = 10
	#highThr = 10 #High threshold for canny
	#accThr = 200; #accumulator threshold for the circle centers at the detection stage. The smaller it is, the more false circles may be detected
	#maxRadius = 50;
	#minRadius = 20;
	#See help for http://opencv.itseez.com/modules/imgproc/doc/feature_detection.html?highlight=houghcircle#cv2.HoughCirclesIn thus
	circles = cv2.HoughCircles(blur,cv2.HOUGH_GRADIENT, dp,minDist, None, highThr,accThr,minRadius, maxRadius)
	#Print the circles
	gColor = cv2.cvtColor(gray, cv2.COLOR_GRAY2BGR)
	if (circles !=None):
		#print circles
		all_circles = circles[0]
		M,N = all_circles.shape
		k=1
		for c in all_circles:
			cv2.circle(gColor, (int(c[0]),int(c[1])),c[2], (int(k*255/M),k*128,0))
			#setText(gColor, (int(c[0]),int(c[1])), "Radius:" + str(c[2]))
			cv2.putText(gColor,"Radius:" + str(c[2]), (int(c[0]),int(c[1])), cv2.FONT_HERSHEY_SIMPLEX, .5, 255)
			K=k+1
			#Circle with max votes
		c=all_circles[0,:]
		cv2.circle(gColor, (int(c[0]),int(c[1])),c[2], (0,0,255))
	#text labels

	setText(gColor, (10, 10), "minDist:" + str(minDist) + " highThr:" + str(highThr) +" accThr:" + str(accThr) +" minRadius:" + str(minRadius) +" maxRadius:" + str(maxRadius) )

	cv2.imshow("houghI",gColor)
#--------------------------
#         UI related
#--------------------------

def setText(dst, (x, y), s):
	cv2.putText(dst, s, (x+1, y+1), cv2.FONT_HERSHEY_PLAIN, 1.0, (0, 0, 0), thickness = 2)
	cv2.putText(dst, s, (x, y), cv2.FONT_HERSHEY_PLAIN, 1.0, (255, 255, 255))

def setupWindowSliders():
	''' Define windows for displaying the results and create trackbars'''
	cv2.namedWindow("Result")
	cv2.namedWindow('Threshold')
	cv2.namedWindow("TempResults")
	cv2.namedWindow("hough")
	cv2.namedWindow("houghI")
	#Threshold value for the pupil intensity
	cv2.createTrackbar('pupilThr','Threshold', 90, 255, onSlidersChange)
	#Threshold value for the glint intensities
	cv2.createTrackbar('glintThr','Threshold', 240, 255,onSlidersChange)
	#define the minimum and maximum areas of the pupil
	cv2.createTrackbar('minSize','Threshold', 20, 200, onSlidersChange)
	cv2.createTrackbar('maxSize','Threshold', 200,200, onSlidersChange)
	cv2.createTrackbar('minSizeG','Threshold', 20, 200, onSlidersChange)
	cv2.createTrackbar('maxSizeG','Threshold', 30,200, onSlidersChange)
	#Value to indicate whether to run or pause the video
	cv2.createTrackbar('Stop/Start','Threshold', 0,1, onSlidersChange)

	#Hough Transform pupils

	cv2.createTrackbar('minDist','hough', 1,100, onSlidersChange)
	cv2.createTrackbar('highThr','hough', 1,50, onSlidersChange)
	cv2.createTrackbar('accThr','hough', 50,800, onSlidersChange)
	cv2.createTrackbar('minRadius','hough', 10,100, onSlidersChange)
	cv2.createTrackbar('maxRadius','hough', 20,300, onSlidersChange)

	#Hough Transform Iris

	cv2.createTrackbar('minDistI','houghI', 10,200, onSlidersChange)
	cv2.createTrackbar('highThrI','houghI', 10,100, onSlidersChange)
	cv2.createTrackbar('accThrI','houghI', 50,800, onSlidersChange)
	cv2.createTrackbar('minRadiusI','houghI', 100,180, onSlidersChange)
	cv2.createTrackbar('maxRadiusI','houghI', 110,180, onSlidersChange)

def getSliderVals():
	'''Extract the values of the sliders and return these in a dictionary'''
	sliderVals={}
	sliderVals['pupilThr'] = cv2.getTrackbarPos('pupilThr', 'Threshold')
	sliderVals['glintThr'] = cv2.getTrackbarPos('glintThr', 'Threshold')
	sliderVals['minSize'] = 50*cv2.getTrackbarPos('minSize', 'Threshold')
	sliderVals['maxSize'] = 50*cv2.getTrackbarPos('maxSize', 'Threshold')
	sliderVals['minSizeG'] = cv2.getTrackbarPos('minSizeG', 'Threshold')
	sliderVals['maxSizeG'] = 50*cv2.getTrackbarPos('maxSizeG', 'Threshold')
	sliderVals['Running'] = 1==cv2.getTrackbarPos('Stop/Start', 'Threshold')

	#Hough pupils
	sliderVals['minDist'] = cv2.getTrackbarPos('minDist', 'hough')
	sliderVals['highThr'] = cv2.getTrackbarPos('highThr', 'hough')
	sliderVals['accThr'] = cv2.getTrackbarPos('accThr', 'hough')
	sliderVals['minRadius'] = cv2.getTrackbarPos('minRadius', 'hough')
	sliderVals['maxRadius'] = cv2.getTrackbarPos('maxRadius', 'hough')

	#Hough Iriss
	sliderVals['minDistI'] = cv2.getTrackbarPos('minDistI', 'houghI')
	sliderVals['highThrI'] = cv2.getTrackbarPos('highThrI', 'houghI')
	sliderVals['accThrI'] = cv2.getTrackbarPos('accThrI', 'houghI')
	sliderVals['minRadiusI'] = cv2.getTrackbarPos('minRadiusI', 'houghI')
	sliderVals['maxRadiusI'] = cv2.getTrackbarPos('maxRadiusI', 'houghI')

	return sliderVals

def onSlidersChange(dummy=None):
	''' Handle updates when slides have changed.
	 This  function only updates the display when the video is put on pause'''
	global imgOrig;
	sv=getSliderVals()
	if(not sv['Running']): # if pause
		update(imgOrig)

#--------------------------
#         main
#--------------------------
run(inputFile)

max_stars_repo_name: tonybeltramelli/Graphics-And-Vision

max_issues_repo_name: tonybeltramelli/Graphics-And-Vision

max_forks_repo_name: tonybeltramelli/Graphics-And-Vision

reason: import numpy,import scipy,from scipy