text:
import math
import numpy as np
import pandas as pd
from pyomo import environ as pe
from pyomo.opt import SolverStatus, TerminationCondition

from hatchet import config
from hatchet.utils.solve.utils import Random

class ILPSubset:
    def __init__(self, n, cn_max, d, mu, ampdel, copy_numbers, f_a, f_b, w):

        # Each ILPSubset maintains its own data, so make a deep-copy of passed-in DataFrames
        f_a, f_b = f_a.copy(deep=True), f_b.copy(deep=True)

        assert f_a.shape == f_b.shape
        assert np.all(f_a.index == f_b.index)
        assert np.all(f_a.columns == f_b.columns)

        self.m, self.k = f_a.shape
        self.f_a = f_a
        self.f_b = f_b
        self.cluster_ids = f_a.index
        self.sample_ids = f_a.columns

        self.n = n
        self.cn_max = cn_max
        self.d = d
        self.mu = mu
        self.ampdel = ampdel
        self.copy_numbers = copy_numbers
        self.w = w

        self.tol = 0.001

        self.mode = 'FULL'

        # Values we want to optimize for, as dataframes
        self.cA = [[np.nan for _ in range(self.n)] for _ in range(self.m)]
        self.cB = [[np.nan for _ in range(self.n)] for _ in range(self.m)]
        self.u = [[np.nan for _ in range(self.k)] for _ in range(self.n)]

        # Fixed values of cA/cB/u
        self._fixed_cA = [[np.nan for _ in range(n)] for _ in range(self.m)]
        self._fixed_cB = [[np.nan for _ in range(n)] for _ in range(self.m)]
        self._fixed_u = [[np.nan for _ in range(self.k)] for _ in range(self.n)]

        self.warmstart = False  # set on hot_start()
        self.model = None  # initialized on create_model()

    def __copy__(self):
        return ILPSubset(
            n=self.n, cn_max=self.cn_max, d=self.d, mu=self.mu, ampdel=self.ampdel, copy_numbers=self.copy_numbers,
            f_a=self.f_a, f_b=self.f_b, w=self.w
        )

    @property
    def M(self):
        return math.floor(math.log2(self.cn_max)) + 1

    @property
    def base(self):
        return min(2, len(self.copy_numbers))

    @staticmethod
    def symmCoeff(i):
        return math.pow(i + 1, 1)

    @property
    def optimized_cA(self):
        if self.mode == 'UARCH':
            return self._fixed_cA
        else:
            return self.cA

    @property
    def optimized_cB(self):
        if self.mode == 'UARCH':
            return self._fixed_cB
        else:
            return self.cB

    @property
    def optimized_u(self):
        if self.mode == 'CARCH':
            return self._fixed_u
        else:
            return self.u

    def create_model(self):

        m, n, k = self.m, self.n, self.k
        f_a, f_b = self.f_a, self.f_b
        cn_max = self.cn_max
        ampdel = self.ampdel
        copy_numbers = self.copy_numbers
        mode_t = self.mode
        d = self.d
        _M = self.M
        _base = self.base

        model = pe.ConcreteModel()

        fA = {}
        fB = {}
        yA = {}
        yB = {}
        adA = {}
        adB = {}

        for _m in range(m):

            cluster_id = f_a.index[_m]

            # upper bound for solver
            ub = max(sum(copy_numbers.get(cluster_id, (0, 0))), cn_max)

            for _k in range(k):
                yA[(_m, _k)] = pe.Var(bounds=(0, np.inf), domain=pe.Reals)
                model.add_component(f'yA_{_m + 1}_{_k + 1}', yA[(_m, _k)])
                yB[(_m, _k)] = pe.Var(bounds=(0, np.inf), domain=pe.Reals)
                model.add_component(f'yB_{_m + 1}_{_k + 1}', yB[(_m, _k)])
                fA[(_m, _k)] = pe.Var(bounds=(0, ub), domain=pe.Reals)
                model.add_component(f'fA_{_m + 1}_{_k + 1}', fA[(_m, _k)])
                fB[(_m, _k)] = pe.Var(bounds=(0, ub), domain=pe.Reals)
                model.add_component(f'fB_{_m + 1}_{_k + 1}', fB[(_m, _k)])

        if mode_t in ('FULL', 'CARCH'):
            for _m in range(m):

                cluster_id = f_a.index[_m]

                # upper bound for solver
                ub = max(sum(copy_numbers.get(cluster_id, (0, 0))), cn_max)

                for _n in range(n):
                    self.cA[_m][_n] = pe.Var(bounds=(0, ub), domain=pe.Integers)
                    model.add_component(f'cA_{_m + 1}_{_n + 1}', self.cA[_m][_n])
                    self.cB[_m][_n] = pe.Var(bounds=(0, ub), domain=pe.Integers)
                    model.add_component(f'cB_{_m + 1}_{_n + 1}', self.cB[_m][_n])

            if ampdel:
                for _m in range(m):
                    cluster_id = f_a.index[_m]
                    if cluster_id not in copy_numbers:
                        adA[_m] = pe.Var(bounds=(0, 1), domain=pe.Binary)
                        model.add_component(f'adA_{_m + 1}', adA[_m])
                        adB[_m] = pe.Var(bounds=(0, 1), domain=pe.Binary)
                        model.add_component(f'adB_{_m + 1}', adB[_m])

        bitcA = {}
        bitcB = {}
        if (mode_t == 'FULL') or (d > 0 and mode_t == 'CARCH'):
            for _b in range(_M):
                for _m in range(m):
                    for _n in range(n):
                        bitcA[(_b, _m, _n)] = pe.Var(bounds=(0, 1), domain=pe.Binary)
                        model.add_component(f'bitcA_{_b + 1}_{_m + 1}_{_n + 1}', bitcA[(_b, _m, _n)])
                        bitcB[(_b, _m, _n)] = pe.Var(bounds=(0, 1), domain=pe.Binary)
                        model.add_component(f'bitcB_{_b + 1}_{_m + 1}_{_n + 1}', bitcB[(_b, _m, _n)])

        if mode_t in ('FULL', 'UARCH'):
            for _n in range(n):
                for _k in range(k):
                    self.u[_n][_k] = pe.Var(bounds=(0, 1), domain=pe.Reals)
                    model.add_component(f'u_{_n + 1}_{_k + 1}', self.u[_n][_k])

        vA = {}
        vB = {}
        if mode_t == 'FULL':
            for _b in range(_M):
                for _m in range(m):
                    for _n in range(n):
                        for _k in range(k):
                            vA[(_b, _m, _n, _k)] = pe.Var(bounds=(0, 1), domain=pe.Reals)
                            model.add_component(f'vA_{_b + 1}_{_m + 1}_{_n + 1}_{_k + 1}', vA[(_b, _m, _n, _k)])
                            vB[(_b, _m, _n, _k)] = pe.Var(bounds=(0, 1), domain=pe.Reals)
                            model.add_component(f'vB_{_b + 1}_{_m + 1}_{_n + 1}_{_k + 1}', vB[(_b, _m, _n, _k)])

        x = {}
        if (mode_t in ('FULL', 'UARCH')) and (self.mu > 0):
            for _n in range(n):
                for _k in range(k):
                    x[(_n, _k)] = pe.Var(domain=pe.Binary)
                    model.add_component(f'x_{_n + 1}_{_k + 1}', x[(_n, _k)])

        # buildOptionalVariables
        z = {}
        if (mode_t in ('FULL', 'CARCH')) and d > 0:
            for _m in range(self.m):
                for _n in range(1, self.n):
                    for _d in range(d):
                        z[(_m, _n, _d)] = pe.Var(bounds=(0, 1), domain=pe.Binary)
                        model.add_component(f'z_{_m + 1}_{_n + 1}_{_d + 1}', z[(_m, _n, _d)])

        # CONSTRAINTS
        model.constraints = pe.ConstraintList()

        for _m in range(m):

            cluster_id = f_a.index[_m]
            f_a_values = f_a.loc[cluster_id].values
            f_b_values = f_b.loc[cluster_id].values

            for _k in range(k):
                _yA, _yB = yA[(_m, _k)], yB[(_m, _k)]
                _fA, _fB = fA[(_m, _k)], fB[(_m, _k)]

                model.constraints.add(float(f_a_values[_k]) - _fA <= _yA)
                model.constraints.add(_fA - float(f_a_values[_k]) <= _yA)
                model.constraints.add(float(f_b_values[_k]) - _fB <= _yB)
                model.constraints.add(_fB - float(f_b_values[_k]) <= _yB)

        if mode_t == 'FULL':
            for _m in range(m):
                for _k in range(k):

                    sum_a = 0
                    for _n in range(n):
                        for _b in range(_M):
                            sum_a += vA[(_b, _m, _n, _k)] * math.pow(2, _b)
                            model.constraints.add(vA[(_b, _m, _n, _k)] <= bitcA[(_b, _m, _n)])
                            model.constraints.add(vA[(_b, _m, _n, _k)] <= self.u[_n][_k])
                            model.constraints.add(vA[(_b, _m, _n, _k)] >= bitcA[(_b, _m, _n)] + self.u[_n][_k] - 1)

                    model.constraints.add(fA[(_m, _k)] == sum_a)

            for _m in range(m):
                for _k in range(k):

                    sum_b = 0
                    for _n in range(n):
                        for _b in range(_M):
                            sum_b += vB[(_b, _m, _n, _k)] * math.pow(2, _b)
                            model.constraints.add(vB[(_b, _m, _n, _k)] <= bitcB[(_b, _m, _n)])
                            model.constraints.add(vB[(_b, _m, _n, _k)] <= self.u[_n][_k])
                            model.constraints.add(vB[(_b, _m, _n, _k)] >= bitcB[(_b, _m, _n)] + self.u[_n][_k] - 1)

                    model.constraints.add(fB[(_m, _k)] == sum_b)

            for _n in range(n):
                for _k in range(k):
                    _sum = 0
                    for _m in range(m):
                        for _b in range(_M):
                            _sum += bitcA[(_b, _m, _n)] + bitcB[(_b, _m, _n)]
                    model.constraints.add(_sum >= self.u[_n][_k])

        if (mode_t == 'FULL') or (d > 0 and mode_t == 'CARCH'):
            for _m in range(m):
                cluster_id = f_a.index[_m]
                # upper bound for solver
                ub = max(sum(copy_numbers.get(cluster_id, (0, 0))), cn_max)

                for _n in range(n):
                    sum_a = 0
                    sum_b = 0
                    for _b in range(_M):
                        sum_a += bitcA[(_b, _m, _n)] * math.pow(2, _b)
                        sum_b += bitcB[(_b, _m, _n)] * math.pow(2, _b)

                    model.constraints.add(self.cA[_m][_n] == sum_a)
                    model.constraints.add(self.cB[_m][_n] == sum_b)
                    model.constraints.add(self.cA[_m][_n] + self.cB[_m][_n] <= ub)

        if mode_t == 'CARCH':
            # TODO: These loops can be collapsed once validation against C++ is complete
            for _m in range(m):
                for _k in range(k):
                    _sumA = 0
                    _sumB = 0
                    for _n in range(n):
                        if self._fixed_u[_n][_k] >= self.mu - self.tol:
                            _sumA += self.cA[_m][_n] * self._fixed_u[_n][_k]
                            _sumB += self.cB[_m][_n] * self._fixed_u[_n][_k]
                    model.constraints.add(fA[(_m, _k)] == _sumA)
                    model.constraints.add(fB[(_m, _k)] == _sumB)

                cluster_id = f_a.index[_m]
                # upper bound for solver
                ub = max(sum(copy_numbers.get(cluster_id, (0, 0))), cn_max)
                for _n in range(n):
                    model.constraints.add(self.cA[_m][_n] + self.cB[_m][_n] <= ub)

        if mode_t in ('FULL', 'CARCH'):
            for _m in range(m):
                model.constraints.add(self.cA[_m][0] == 1)
                model.constraints.add(self.cB[_m][0] == 1)

            if ampdel:
                for _m in range(m):
                    cluster_id = f_a.index[_m]
                    if cluster_id not in copy_numbers:
                        for _n in range(1, n):
                            model.constraints.add(self.cA[_m][_n] <= cn_max * adA[_m] + _base - _base * adA[_m])
                            model.constraints.add(self.cA[_m][_n] >= _base * adA[_m])
                            model.constraints.add(self.cB[_m][_n] <= cn_max * adB[_m] + _base - _base * adB[_m])
                            model.constraints.add(self.cB[_m][_n] >= _base * adB[_m])

        if mode_t == 'UARCH':
            # TODO: These loops can be collapsed once validation against C++ is complete
            for _m in range(m):
                for _k in range(k):
                    _sumA = 0
                    _sumB = 0
                    for _n in range(n):
                        _sumA += int(self._fixed_cA[_m][_n]) * self.u[_n][_k]
                        _sumB += int(self._fixed_cB[_m][_n]) * self.u[_n][_k]
                    model.constraints.add(fA[(_m, _k)] == _sumA)
                    model.constraints.add(fB[(_m, _k)] == _sumB)

        if mode_t in ('FULL', 'UARCH'):
            for _k in range(k):
                _sum = 0
                for _n in range(n):
                    _sum += self.u[_n][_k]
                model.constraints.add(_sum == 1)

        if (mode_t in ('FULL', 'UARCH')) and self.mu > 0:
            for _k in range(k):
                for _n in range(1, n):
                    model.constraints.add(x[(_n, _k)] >= self.u[_n][_k])
                    model.constraints.add(self.u[_n][_k] >= self.mu * x[(_n, _k)])

        # buildOptionalConstraints
        if (mode_t in ('FULL', 'CARCH')) and d > 0:
            for _m in range(self.m):
                for _n in range(1, self.n):
                    _sum = 0
                    for _d in range(self.d):
                        _sum += z[(_m, _n, _d)]
                    model.constraints.add(_sum == 1)

            for _m in range(self.m):
                for _M in range(self.M):
                    for _d in range(d):
                        for _i in range(1, self.n - 1):
                            for _j in range(1, self.n):
                                model.constraints.add(
                                    bitcA[(_M, _m, _i)] - bitcA[(_M, _m, _j)] <= 2 - z[(_m, _i, _d)] - z[(_m, _j, _d)])
                                model.constraints.add(
                                    bitcA[(_M, _m, _j)] - bitcA[(_M, _m, _i)] <= 2 - z[(_m, _i, _d)] - z[(_m, _j, _d)])
                                model.constraints.add(
                                    bitcB[(_M, _m, _i)] - bitcB[(_M, _m, _j)] <= 2 - z[(_m, _i, _d)] - z[(_m, _j, _d)])
                                model.constraints.add(
                                    bitcB[(_M, _m, _j)] - bitcB[(_M, _m, _i)] <= 2 - z[(_m, _i, _d)] - z[(_m, _j, _d)])

            for _m in range(self.m):
                for _d in range(d - 1):
                    _sum_l = _sum_l1 = 0
                    for _n in range(1, self.n):
                        _sum_l += z[(_m, _n, _d)] * self.symmCoeff(_n)
                        _sum_l1 += z[(_m, _n, _d + 1)] * self.symmCoeff(_n)
                        model.constraints.add(_sum_l <= _sum_l1)

        if mode_t in ('FULL', 'CARCH'):
            self.build_symmetry_breaking(model)
            self.fix_given_cn(model)

        if mode_t == 'FULL':
            self.hot_start()

        obj = 0
        for _m in range(m):
            for _k in range(k):
                cluster_id = self.cluster_ids[_m]
                obj += (yA[(_m, _k)] + yB[(_m, _k)]) * self.w[cluster_id]

        model.obj = pe.Objective(expr=obj, sense=pe.minimize)
        self.model = model

    def build_symmetry_breaking(self, model):
        for i in range(1, self.n - 1):
            _sum1 = 0
            _sum2 = 0
            for _m in range(self.m):
                _symCoeff = self.symmCoeff(_m)
                _sum1 += self.cA[_m][i] * _symCoeff + self.cB[_m][i] * _symCoeff
                _sum2 += self.cA[_m][i] * _symCoeff + self.cB[_m][i] * _symCoeff
            model.constraints.add(_sum1 <= _sum2)

    def fix_given_cn(self, model):
        for _m in range(self.m):
            cluster_id = self.f_a.index[_m]
            if cluster_id in self.copy_numbers:
                _cnA, _cnB = self.copy_numbers[cluster_id]
                for _n in range(1, self.n):
                    model.constraints.add(self.cA[_m][_n] == _cnA)
                    model.constraints.add(self.cB[_m][_n] == _cnB)

    def first_hot_start(self):

        if self.d > 0:
            targetA = np.empty((self.m, self.d))
            targetB = np.empty_like(targetA)
            for _m in range(self.m):
                targetA[_m, :] = np.random.choice(self.f_a.iloc[_m].values, self.d, replace=False)
                targetB[_m, :] = np.random.choice(self.f_b.iloc[_m].values, self.d, replace=False)
            targetA = targetA.round()
            targetB = targetB.round()
        else:
            targetA = np.round(self.f_a.values)
            targetB = np.round(self.f_b.values)

        hcA = np.zeros((self.m, self.n))
        hcB = np.zeros((self.m, self.n))

        for _m in range(self.m):
            cluster_id = self.f_a.index[_m]

            adA = 0
            adB = 0
            for _n in range(self.n):
                if _n == 0:
                    hcA[_m][0] = 1
                    hcB[_m][0] = 1
                else:
                    if cluster_id in self.copy_numbers:
                        hcA[_m][_n] = self.copy_numbers[cluster_id][0]
                        hcB[_m][_n] = self.copy_numbers[cluster_id][1]
                    else:
                        mod = min(self.n, self.k)
                        a = min(targetA[_m][_n % mod], self.cn_max)
                        b = min(targetB[_m][_n % mod], self.cn_max)

                        if self.ampdel:
                            base = self.base
                            if adA == 0 and a > base:
                                adA = 1
                            if adA == 0 and a < base:
                                adA = -1
                            if adB == 0 and b > base:
                                adB = 1
                            if adB == 0 and b < base:
                                adB = -1

                            a = max(a, base) if adA >= 0 else min(a, base)
                            b = max(b, base) if adB >= 0 else min(b, base)

                            if a + b > self.cn_max:
                                a = self.cn_max - base
                                b = self.cn_max - a

                        if a + b <= self.cn_max:
                            hcA[_m][_n] = a
                            hcB[_m][_n] = b
                        else:
                            hcA[_m][_n] = a
                            hcB[_m][_n] = max(self.cn_max - a, 0)

                        assert hcA[_m][_n] + hcB[_m][_n] <= self.cn_max
                        assert hcA[_m][_n] >= 0
                        assert hcB[_m][_n] >= 0

        return hcA, hcB

    def hot_start(self, f_a=None, f_b=None):
        def get_rank(hcA, hcB):
            # m, n = hcA.shape
            m, n = len(hcA), len(hcA[0])
            rank = np.zeros(n)
            rank[0] = -1
            for _m in range(m):
                for _n in range(1, n):
                    rank[_n] += hcA[_m][_n] * self.symmCoeff(_m) + hcB[_m][_n] * self.symmCoeff(_m)

            return rank

        if f_a is None and f_b is None:
            f_a, f_b = self.first_hot_start()
        rank = get_rank(f_a, f_b)
        rank_indices = np.argsort(rank)

        for _m in range(self.m):
            for _n in range(0, self.n):
                self.cA[_m][rank_indices[_n]].value = f_a[_m][_n]
                self.cB[_m][rank_indices[_n]].value = f_b[_m][_n]
        self.warmstart = True

    def fix_u(self, u):
        self._fixed_u[:] = u
        self.mode = 'CARCH'
        # TODO: sanity checks in fixU

    def fix_c(self, cA, cB):
        self._fixed_cA = cA
        self._fixed_cB = cB
        self.mode = 'UARCH'
        # TODO: sanity checks in fixC

    def build_random_u(self, random_seed=None):
        def _calculate_size_bubbles(mu):
            if mu <= 0.1:
                return 10
            elif mu <= 0.15:
                return 6
            elif mu <= 0.2:
                return 5
            else:
                return 3

        def _build_partition_vector(n, n_parts, size_bubbles, mu=0.03):
            """
            Return fractional components of n components as an ndarray of size n
            """
            # initialize all fractions to 0
            result = np.zeros(n)
            # generate n_parts random positions from all available positions
            positions = np.random.choice(np.arange(n), n_parts, replace=False)

            bubbles = np.sort(np.random.choice(np.arange(1, size_bubbles), n_parts - 1, replace=False)) / size_bubbles
            _result = np.diff(bubbles, prepend=0, append=1)

            # set fractions at selected positions to successive difference between bubbles
            result[positions] = _result
            # for fractions that are within tol of mu, clamp them up to mu
            result[(mu - self.tol <= result) & (result < mu)] = mu

            return result

        size_bubbles = _calculate_size_bubbles(self.mu)
        U = np.empty((self.n, self.k))

        with Random(random_seed):
            for _k in range(self.k):
                # Generate 2 random ints between [1, n] and take the max - biasing towards larger numbers
                _n0 = np.random.randint(1, self.n + 1)
                _n1 = np.random.randint(1, self.n + 1)
                n_parts = min(max(_n0, _n1), size_bubbles)  # no. of clones with non-zero proportion in the mix
                v = _build_partition_vector(self.n, n_parts, size_bubbles, mu=self.mu)
                U[:, _k] = v
        return U

    def run(self, solver_type='gurobi', timelimit=None, write_path=None):
        if solver_type == 'gurobipy':
            solver = pe.SolverFactory('gurobi', solver_io='python')
        else:
            solver = pe.SolverFactory(solver_type)

        kwargs = {'report_timing': False}
        if timelimit is not None:
            kwargs['timelimit'] = int(timelimit)
        if solver.warm_start_capable():
            kwargs['warmstart'] = self.warmstart

        results = solver.solve(self.model, **kwargs)
        if (
                (results.solver.status == SolverStatus.ok) and (
                results.solver.termination_condition in (TerminationCondition.optimal,
                                                         TerminationCondition.feasible))
                or
                (results.solver.status == SolverStatus.aborted and
                 results.solver.termination_condition == TerminationCondition.maxTimeLimit)
        ):
            pass
        else:
            return None

        if write_path is not None:
            self.model.write(write_path)

        return (
            self.model.obj(),
            [[int(getattr(x, 'value', x)) for x in row] for row in self.optimized_cA],
            [[int(getattr(x, 'value', x)) for x in row] for row in self.optimized_cB],
            [[getattr(x, 'value', x) for x in row] for row in self.optimized_u],
            self.cluster_ids,
            self.sample_ids
        )

reason: import numpy,from pyomo


text:
import tensorflow as tf
import tensorflow_datasets as tfds
import tensorflow.keras as tfk
from sfh.datasets.tng100 import tng100

# Global variables used as parameter for preprocessing the data
mean_stellar_light = 167915.38
mean_velocity      = 1.2179685
mean_velocity_disp = 80.03517
std_velocity       = 47.085832102830075
std_velocity_disp  = 29.16650153197798
std_stellar_light  = 18212.175615239412
stellar_light_compression = 3.0

def preprocessing(example):
  import numpy as np

  img0 = example['stellar_light']
  img1 = example['velocity_map']
  img2 = example['velocity_dispersion_map']

  # Replace NaNs by zeros
  img0 = tf.where(tf.math.is_nan(img0), tf.zeros_like(img0), img0)
  img1 = tf.where(tf.math.is_nan(img1), tf.zeros_like(img1), img1)
  img2 = tf.where(tf.math.is_nan(img2), tf.zeros_like(img2), img2)
  # Replace InFs by zeros
  img0 = tf.where(tf.math.is_inf(img0), tf.zeros_like(img0), img0)
  img1 = tf.where(tf.math.is_inf(img1), tf.zeros_like(img1), img1)
  img2 = tf.where(tf.math.is_inf(img2), tf.zeros_like(img2), img2)

  # Normalize data
  img0 = tf.math.asinh(img0 / tf.constant(std_stellar_light) * tf.constant(stellar_light_compression) ) / tf.constant(stellar_light_compression)
  img1 = (img1 - tf.constant(mean_velocity)) / tf.constant(std_velocity)
  img2 = (img2 - tf.constant(mean_velocity_disp))/ tf.constant(std_velocity_disp)

  # Reduce images size (drop stellar light)
  img = tf.stack([img1, img2], axis=-1)
  # Resize images
  img = tf.image.resize(img, [128, 128])

  return img, example['last_major_merger']

def input_fn(mode='train', batch_size=64):
  """
  mode: 'train' or 'test'
  """
  # Jean-Zay datasets diretory:
  data_dir='/gpfsscratch/rech/wvb/commun/tensorflow_datasets'
  #data_dir='/Users/benjamin/SCRATCH/sfh/content/data/'

  # 3 cases: train (75%), test (15%) and valid (10%). Full dataset (100%) by default
  if (mode == 'train'):
    dataset = tfds.load('tng100', split='train[:75%]', data_dir=data_dir)
    dataset = dataset.map(preprocessing) # Apply data preprocessing
    dataset = dataset.shuffle(20000)
    dataset = dataset.repeat()
  elif (mode == 'test'):
    dataset = tfds.load('tng100', split='train[75%:90%]', data_dir=data_dir)
    dataset = dataset.map(preprocessing) # Apply data preprocessing
  elif (mode == 'valid'):
    dataset = tfds.load('tng100', split='train[90%:100%]', data_dir=data_dir)
    dataset = dataset.map(preprocessing) # Apply data preprocessing
  else:
    dataset = tfds.load('tng100', split='train', data_dir=data_dir)
    dataset = dataset.map(preprocessing) # Apply data preprocessing

  dataset = dataset.batch(batch_size, drop_remainder=True)
  dataset = dataset.prefetch(-1) # fetch next batches while training current one (-1 for autotune)

  return dataset

ext: py

max_stars_repo_name: mhuertascompany/sfh-inference

max_issues_repo_name: mhuertascompany/sfh-inference

max_forks_repo_name: mhuertascompany/sfh-inference

reason: import numpy


text:
##############################################
# sudo apt-get install -y python3-picamera
# sudo -H pip3 install imutils --upgrade
##############################################

import multiprocessing as mp
import sys
from time import sleep

import argparse
import cv2
import numpy as np
import time

try:
    from armv7l.openvino.inference_engine import IENetwork, IEPlugin
except:
    from openvino.inference_engine import IENetwork, IEPlugin
import heapq
import threading

try:
    from imutils.video.pivideostream import PiVideoStream
    from imutils.video.filevideostream import FileVideoStream
    import imutils
except:
    pass

lastresults = None
threads = []
processes = []
frameBuffer = None
results = None
fps = ""
detectfps = ""
framecount = 0
detectframecount = 0
time1 = 0
time2 = 0
cam = None
vs = None
window_name = ""
elapsedtime = 0.0

g_plugin = None
g_inferred_request = None
g_heap_request = None
g_inferred_cnt = 0
g_number_of_allocated_ncs = 0

LABELS = ["neutral", "happy", "sad", "surprise", "anger"]
COLORS = np.random.uniform(0, 255, size=(len(LABELS), 3))

def camThread(LABELS, resultsEm, frameBuffer, camera_width, camera_height, vidfps, number_of_camera, mode_of_camera):
    global fps
    global detectfps
    global lastresults
    global framecount
    global detectframecount
    global time1
    global time2
    global cam
    global vs
    global window_name

    if mode_of_camera == 0:
        cam = cv2.VideoCapture(number_of_camera)
        if cam.isOpened() != True:
            print("USB Camera Open Error!!!")
            sys.exit(0)
        cam.set(cv2.CAP_PROP_FPS, vidfps)
        cam.set(cv2.CAP_PROP_FRAME_WIDTH, camera_width)
        cam.set(cv2.CAP_PROP_FRAME_HEIGHT, camera_height)
        window_name = "USB Camera"
    else:
        vs = PiVideoStream((camera_width, camera_height), vidfps).start()
        sleep(3)
        window_name = "PiCamera"

    cv2.namedWindow(window_name, cv2.WINDOW_AUTOSIZE)

    while True:
        t1 = time.perf_counter()

        # USB Camera Stream or PiCamera Stream Read
        color_image = None
        if mode_of_camera == 0:
            s, color_image = cam.read()
            if not s:
                continue
        else:
            color_image = vs.read()

        if frameBuffer.full():
            frameBuffer.get()
        frames = color_image

        height = color_image.shape[0]
        width = color_image.shape[1]
        frameBuffer.put(color_image.copy())
        res = None

        if not resultsEm.empty():
            res = resultsEm.get(False)
            # print("[LOG] ".format(type(res)))
            # print(res)
            detectframecount += 1
            imdraw = overlay_on_image(frames, res)
            lastresults = res
        else:
            imdraw = overlay_on_image(frames, lastresults)

        cv2.imshow(window_name, cv2.resize(imdraw, (width, height)))

        if cv2.waitKey(1) & 0xFF == ord('q'):
            sys.exit(0)

        ## Print FPS
        framecount += 1
        if framecount >= 25:
            fps = "(Playback) {:.1f} FPS".format(time1 / 25)
            detectfps = "(Detection) {:.1f} FPS".format(detectframecount / time2)
            framecount = 0
            detectframecount = 0
            time1 = 0
            time2 = 0
        t2 = time.perf_counter()
        elapsedTime = t2 - t1
        time1 += 1 / elapsedTime
        time2 += elapsedTime

# l = Search list
# x = Search target value
def searchlist(l, x, notfoundvalue=-1):
    if x in l:
        return l.index(x)
    else:
        return notfoundvalue

def async_infer(ncsworkerFd, ncsworkerEm):
    while True:
        ncsworkerFd.predict_async()
        ncsworkerEm.predict_async()

class BaseNcsWorker():

    def __init__(self, devid, model_path, number_of_ncs):
        global g_plugin
        global g_inferred_request
        global g_heap_request
        global g_inferred_cnt
        global g_number_of_allocated_ncs

        self.devid = devid
        if number_of_ncs == 0:
            self.num_requests = 4
        elif number_of_ncs == 1:
            self.num_requests = 4
        elif number_of_ncs == 2:
            self.num_requests = 2
        elif number_of_ncs >= 3:
            self.num_requests = 1

        print("g_number_of_allocated_ncs =", g_number_of_allocated_ncs, "number_of_ncs =", number_of_ncs)

        if g_number_of_allocated_ncs < 1:
            self.plugin = IEPlugin(device="MYRIAD")
            self.inferred_request = [0] * self.num_requests
            self.heap_request = []
            self.inferred_cnt = 0
            g_plugin = self.plugin
            g_inferred_request = self.inferred_request
            g_heap_request = self.heap_request
            g_inferred_cnt = self.inferred_cnt
            g_number_of_allocated_ncs += 1
        else:
            self.plugin = g_plugin
            self.inferred_request = g_inferred_request
            self.heap_request = g_heap_request
            self.inferred_cnt = g_inferred_cnt

        self.model_xml = model_path + ".xml"
        self.model_bin = model_path + ".bin"
        self.net = IENetwork(model=self.model_xml, weights=self.model_bin)
        self.input_blob = next(iter(self.net.inputs))
        self.exec_net = self.plugin.load(network=self.net, num_requests=self.num_requests)

class NcsWorkerFd(BaseNcsWorker):

    def __init__(self, devid, frameBuffer, resultsFd, model_path, number_of_ncs):

        super().__init__(devid, model_path, number_of_ncs)
        self.frameBuffer = frameBuffer
        self.resultsFd = resultsFd

    def image_preprocessing(self, color_image):

        prepimg = cv2.resize(color_image, (300, 300))
        prepimg = prepimg[np.newaxis, :, :, :]  # Batch size axis add
        prepimg = prepimg.transpose((0, 3, 1, 2))  # NHWC to NCHW
        return prepimg

    def predict_async(self):
        try:

            if self.frameBuffer.empty():
                return

            color_image = self.frameBuffer.get()
            prepimg = self.image_preprocessing(color_image)
            reqnum = searchlist(self.inferred_request, 0)

            if reqnum > -1:
                self.exec_net.start_async(request_id=reqnum, inputs={self.input_blob: prepimg})
                self.inferred_request[reqnum] = 1
                self.inferred_cnt += 1

                if self.inferred_cnt == sys.maxsize:
                    self.inferred_request = [0] * self.num_requests
                    self.heap_request = []
                    self.inferred_cnt = 0

                self.exec_net.requests[reqnum].wait(-1)
                out = self.exec_net.requests[reqnum].outputs["detection_out"].flatten()

                detection_list = []
                face_image_list = []

                for detection in out.reshape(-1, 7):

                    confidence = float(detection[2])

                    if confidence > 0.3:
                        detection[3] = int(detection[3] * color_image.shape[1])
                        detection[4] = int(detection[4] * color_image.shape[0])
                        detection[5] = int(detection[5] * color_image.shape[1])
                        detection[6] = int(detection[6] * color_image.shape[0])
                        if (detection[6] - detection[4]) > 0 and (detection[5] - detection[3]) > 0:
                            detection_list.extend(detection)
                            face_image_list.extend([color_image[int(detection[4]):int(detection[6]),
                                                    int(detection[3]):int(detection[5]), :]])

                if len(detection_list) > 0:
                    self.resultsFd.put([detection_list, face_image_list])

                self.inferred_request[reqnum] = 0

        except:
            import traceback
            traceback.print_exc()

class NcsWorkerEm(BaseNcsWorker):

    def __init__(self, devid, resultsFd, resultsEm, model_path, number_of_ncs):

        super().__init__(devid, model_path, number_of_ncs)
        self.resultsFd = resultsFd
        self.resultsEm = resultsEm

    def image_preprocessing(self, color_image):

        try:
            prepimg = cv2.resize(color_image, (64, 64))
        except:
            prepimg = np.full((64, 64, 3), 128)
        prepimg = prepimg[np.newaxis, :, :, :]  # Batch size axis add
        prepimg = prepimg.transpose((0, 3, 1, 2))  # NHWC to NCHW
        return prepimg

    def predict_async(self):
        try:

            if self.resultsFd.empty():
                return

            resultFd = self.resultsFd.get()
            detection_list = resultFd[0]
            face_image_list = resultFd[1]
            emotion_list = []
            max_face_image_list_cnt = len(face_image_list)
            image_idx = 0
            end_cnt_processing = 0
            heapflg = False
            cnt = 0
            dev = 0

            if max_face_image_list_cnt <= 0:
                detection_list.extend([""])
                self.resultsEm.put([detection_list])
                return

            while True:
                reqnum = searchlist(self.inferred_request, 0)

                if reqnum > -1 and image_idx <= (max_face_image_list_cnt - 1) and len(face_image_list[image_idx]) > 0:

                    if len(face_image_list[image_idx]) == []:
                        image_idx += 1
                        continue
                    else:
                        prepimg = self.image_preprocessing(face_image_list[image_idx])
                        image_idx += 1

                    self.exec_net.start_async(request_id=reqnum, inputs={self.input_blob: prepimg})
                    self.inferred_request[reqnum] = 1
                    self.inferred_cnt += 1
                    if self.inferred_cnt == sys.maxsize:
                        self.inferred_request = [0] * self.num_requests
                        self.heap_request = []
                        self.inferred_cnt = 0
                    heapq.heappush(self.heap_request, (self.inferred_cnt, reqnum))
                    heapflg = True

                if heapflg:
                    cnt, dev = heapq.heappop(self.heap_request)
                    heapflg = False

                if self.exec_net.requests[dev].wait(0) == 0:
                    self.exec_net.requests[dev].wait(-1)
                    out = self.exec_net.requests[dev].outputs["prob_emotion"].flatten()
                    emotion = LABELS[int(np.argmax(out))]
                    detection_list.extend([emotion])
                    self.resultsEm.put([detection_list])
                    self.inferred_request[dev] = 0
                    end_cnt_processing += 1
                    if end_cnt_processing >= max_face_image_list_cnt:
                        break
                else:
                    heapq.heappush(self.heap_request, (cnt, dev))
                    heapflg = True

        except:
            import traceback
            traceback.print_exc()

def inferencer(resultsFd, resultsEm, frameBuffer, number_of_ncs, fd_model_path, em_model_path):
    # Init infer threads
    threads = []
    for devid in range(number_of_ncs):
        # Face Detection, Emotion Recognition start
        thworker = threading.Thread(target=async_infer,
                                    args=(NcsWorkerFd(devid, frameBuffer, resultsFd, fd_model_path, number_of_ncs),
                                          NcsWorkerEm(devid, resultsFd, resultsEm, em_model_path, 0),))
        thworker.start()
        threads.append(thworker)
        print("Thread-" + str(devid))

    for th in threads:
        th.join()

def overlay_on_image(frames, object_infos):
    try:

        color_image = frames

        if isinstance(object_infos, type(None)):
            return color_image

        # Show images
        height = color_image.shape[0]
        width = color_image.shape[1]
        img_cp = color_image.copy()
        for object_info in object_infos:

            if object_info[2] == 0.0:
                break

            if (not np.isfinite(object_info[0]) or
                    not np.isfinite(object_info[1]) or
                    not np.isfinite(object_info[2]) or
                    not np.isfinite(object_info[3]) or
                    not np.isfinite(object_info[4]) or
                    not np.isfinite(object_info[5]) or
                    not np.isfinite(object_info[6])):
                continue

            min_score_percent = 60

            percentage = int(object_info[2] * 100)

            if (percentage <= min_score_percent):
                continue

            box_left = int(object_info[3])
            box_top = int(object_info[4])
            box_bottom = int(object_info[6])
            emotion = str(object_info[7])

            label_text = emotion + " (" + str(percentage) + "%)"
            label_text_color = (255, 255, 255)

            # info fps
            cv2.putText(img_cp, fps, (width - 170, 15), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (38, 0, 255), 1, cv2.LINE_AA)
            cv2.putText(img_cp, detectfps, (width - 170, 30), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (38, 0, 255), 1,
                        cv2.LINE_AA)

            # background of expression list
            overlay = img_cp.copy()
            opacity = 0.4
            cv2.rectangle(img_cp, (box_left + box_bottom + 10 - 250, box_top - 25),
                          (box_left + box_bottom - 50, box_top + 25),
                          (64, 64, 64), cv2.FILLED)
            cv2.addWeighted(overlay, opacity, img_cp, 1 - opacity, 0, img_cp)

            # connect face and expressions
            cv2.line(img_cp, (int((box_left + box_left + box_bottom - 250) / 2), box_top + 15),
                     (box_left + box_bottom - 250, box_top - 20),
                     (255, 255, 255), 1)
            cv2.line(img_cp, (box_left + box_bottom - 250, box_top - 20),
                     (box_left + box_bottom + 10 - 250, box_top - 20),
                     (255, 255, 255), 1)

            cv2.putText(img_cp, label_text, (int(box_left + box_bottom + 15 - 250), int(box_top - 12 + 20)),
                        cv2.FONT_HERSHEY_SIMPLEX, 0.5, label_text_color,
                        1)

        return img_cp

    except:
        import traceback
        traceback.print_exc()

if __name__ == '__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument('-cm', '--modeofcamera', dest='mode_of_camera', type=int, default=0,
                        help='Camera Mode. 0:=USB Camera, 1:=PiCamera (Default=0)')
    parser.add_argument('-cn', '--numberofcamera', dest='number_of_camera', type=int, default=0,
                        help='USB camera number. (Default=0)')
    parser.add_argument('-wd', '--width', dest='camera_width', type=int, default=640,
                        help='Width of the frames in the video stream. (Default=640)')
    parser.add_argument('-ht', '--height', dest='camera_height', type=int, default=480,
                        help='Height of the frames in the video stream. (Default=480)')
    parser.add_argument('-numncs', '--numberofncs', dest='number_of_ncs', type=int, default=1,
                        help='Number of NCS. (Default=1)')
    parser.add_argument('-vidfps', '--fpsofvideo', dest='fps_of_video', type=int, default=30,
                        help='FPS of Video. (Default=30)')
    parser.add_argument('-fdmp', '--facedetectionmodelpath', dest='fd_model_path',
                        default='./model/face-detection-retail-0004',
                        help='Face Detection model path. (xml and bin. Except extension.)')
    parser.add_argument('-emmp', '--emotionrecognitionmodelpath', dest='em_model_path',
                        default='./model/emotions-recognition-retail-0003',
                        help='Emotion Recognition model path. (xml and bin. Except extension.)')

    args = parser.parse_args()
    mode_of_camera = args.mode_of_camera
    number_of_camera = args.number_of_camera
    camera_width = args.camera_width
    camera_height = args.camera_height
    number_of_ncs = args.number_of_ncs
    vidfps = args.fps_of_video
    fd_model_path = args.fd_model_path
    em_model_path = args.em_model_path

    try:

        mp.set_start_method('forkserver')
        frameBuffer = mp.Queue(10)
        resultsFd = mp.Queue()  # Face Detection Queue
        resultsEm = mp.Queue()  # Emotion Recognition Queue

        # Start streaming
        p = mp.Process(target=camThread,
                       args=(LABELS, resultsEm, frameBuffer, camera_width, camera_height, vidfps, number_of_camera,
                             mode_of_camera),
                       daemon=True)
        p.start()
        processes.append(p)

        # Start detection MultiStick
        # Activation of inferencer
        p = mp.Process(target=inferencer,
                       args=(resultsFd, resultsEm, frameBuffer, number_of_ncs, fd_model_path, em_model_path),
                       daemon=True)
        p.start()
        processes.append(p)

        while True:
            sleep(1)

    except:
        import traceback

        traceback.print_exc()
    finally:
        for p in range(len(processes)):
            processes[p].terminate()

        print("\n\nFinished\n\n")


text:
import tensorflow as tf
from nptyping import NDArray
import cv2 as cv
from pathlib import Path
import numpy as np

FLAGS = tf.app.flags.FLAGS
tf.app.flags.DEFINE_string('template_dir', 'data', """Path to template image, template image has to be name as template.""")
tf.app.flags.DEFINE_string('image_path', 'data/KITTI/testing/image_2/IMG_7178.jpg', """Path to template image, template image has to be name as template.""")

def select_tmp(input_image: NDArray, tmp_dir: str):
    resize_input_image = cv.resize(input_image, (500, 600))
    template_path_list = list(Path(tmp_dir).glob('**/template*'))
    similarity_list = [compute_similarity(resize_input_image, template_path) for template_path in template_path_list]
    return template_path_list[similarity_list.index(min(similarity_list))]

def compute_similarity(image: NDArray, template_path: str):
    template_img = cv.resize(cv.imread(str(template_path)), (500, 600))
    feature_map_result = map_feature(image, template_img, .75)
    cv.imwrite('./unwarpped'+template_path.name,  feature_map_result)

    template_img = cv.medianBlur(template_img, 11)
    feature_map_result = cv.medianBlur(feature_map_result, 11)
    blured_img = cv.GaussianBlur(template_img, (201, 201), 0)
    return np.mean((image - feature_map_result)**2)

def map_feature(image, template, ratio_thresh):
    akaze = cv.AKAZE_create()
    key_points, descriptors = akaze.detectAndCompute(template, None)
    key_points2, descriptors2 = akaze.detectAndCompute(image, None)
    matcher = cv.DescriptorMatcher_create(cv.DescriptorMatcher_BRUTEFORCE_HAMMING)
    knn_matches = matcher.knnMatch(descriptors, descriptors2, 2)
    good_matches = [m for m, n in knn_matches if m.distance < ratio_thresh * n.distance]
    img_matches = np.empty((max(template.shape[0], image.shape[0]), template.shape[1]+image.shape[1], 3), dtype=np.uint8)
    cv.drawMatches(template, key_points, image, key_points2, good_matches, img_matches, flags=cv.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)
    obj = np.empty((len(good_matches), 2), dtype=np.float32)
    scene = np.empty((len(good_matches), 2), dtype=np.float32)
    for i in range(len(good_matches)):
        obj[i, 0] = key_points[good_matches[i].queryIdx].pt[0]
        obj[i, 1] = key_points[good_matches[i].queryIdx].pt[1]
        scene[i, 0] = key_points2[good_matches[i].trainIdx].pt[0]
        scene[i, 1] = key_points2[good_matches[i].trainIdx].pt[1]
    H, _ = cv.findHomography(scene, obj, cv.RANSAC)
    return cv.warpPerspective(src=image, M=H, dsize=template.shape[1::-1])

if __name__ == "__main__":
    input_image = cv.imread(FLAGS.image_path)
    print(select_tmp(input_image=input_image, tmp_dir=FLAGS.template_dir))

max_stars_repo_name: eclipsesource/squeezeDet

max_issues_repo_name: eclipsesource/squeezeDet

max_forks_repo_name: eclipsesource/squeezeDet


text:
# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import csv
import numpy as np
import os
import sys

from observations.util import maybe_download_and_extract

def fertility(path):
  """Fertility

  Fertility measurements for a sample of women

  A dataset with 333 observations on the following 10 variables.

  `Age`

  Age (in years)

  `LowAFC`

  Smallest antral follicle count

  `MeanAFC`

  Average antral follicle count

  `FSH`

  Maximum follicle stimulating hormone level

  `E2`

  Fertility level

  `MaxE2`

  Maximum fertility level

  `MaxDailyGn`

  Maximum daily gonadotropin level

  `TotalGn`

  Total gonadotropin level

  `Oocytes`

  Number of egg cells

  `Embryos`

  Number of embryos

  We thank Dr. Priya Maseelall and her research team for sharing these

  Args:

    path: str.
      Path to directory which either stores file or otherwise file will
      be downloaded and extracted there.
      Filename is `fertility.csv`.

  Returns:

    Tuple of np.ndarray `x_train` with 333 rows and 10 columns and
    dictionary `metadata` of column headers (feature names).
  """
  import pandas as pd
  path = os.path.expanduser(path)
  filename = 'fertility.csv'
  if not os.path.exists(os.path.join(path, filename)):
    url = 'http://dustintran.com/data/r/Stat2Data/Fertility.csv'
    maybe_download_and_extract(path, url,
                               save_file_name='fertility.csv',
                               resume=False)

  data = pd.read_csv(os.path.join(path, filename), index_col=0,
                     parse_dates=True)
  x_train = data.values
  metadata = {'columns': data.columns}
  return x_train, metadata


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

import numpy as np
import os
import sys
import theano
import theano.tensor as T

class Corpus(object):
    def __init__(self, corpus_path):
        self.corpus = {}
        self.length = 0

        for fname in os.listdir(corpus_path):
            fpath = os.path.join(corpus_path, fname)

            with open(fpath, "r") as f:
                self.corpus[fname.replace(".txt", "")] = f.read().decode("utf-8")

        characters = set()

        for work_name, work in self.corpus.iteritems():
            for c in work:
                characters.add(c)
                self.length += 1

        self.characters = sorted(characters)

    def character_encoder(self, char):
        return self.characters.index(char)

    def __iter__(self):
        for work_name, work in self.corpus.iteritems():
            for char in work:
                yield self.character_encoder(char)
            yield self.character_encoder(u"\n")

    def __len__(self):
        return self.length

print >> sys.stderr, "Getting corpus"
sys.stderr.flush()
corpus = Corpus("corpus/borges")

NT = len(corpus)  # Number of examples (timesteps)
n_in = len(corpus.characters)  # Size of the input data (one-hot vector of a character)
n_out = len(corpus.characters)  # Size of the output data (one-hot vector of a character)
n_h = 50  # Size of the hidden layer

print >> sys.stderr, "Declaring Theano variables"
sys.stderr.flush()
# Stateless variables to handle the input
X = T.ivector('X')
y = T.ivector('y')

W_hx = theano.shared(
    value=np.random.uniform(
        low=-1.0,
        high=1.0,
        size=(n_h, n_in)
    ).astype(theano.config.floatX),
    name='W_hx',
    borrow=True
)

b_h = theano.shared(
    value=np.zeros(n_h, dtype=theano.config.floatX),
    name='b_h',
    borrow=True
)

W_hh = theano.shared(
    value=np.random.uniform(
        low=-1.0,
        high=1.0,
        size=(n_h, n_h)
    ).astype(theano.config.floatX),
    name='W_hh',
    borrow=True
)

W_S = theano.shared(
    value=np.random.uniform(
        low=-1.0,
        high=1.0,
        size=(n_out, n_h)
    ).astype(theano.config.floatX),
    name='W_S',
    borrow=True
)

b_S = theano.shared(
    value=np.zeros(n_out, dtype=theano.config.floatX),
    name='b_S',
    borrow=True
)

def forward_propagation_step(x_t, h_t_prev, W_hx, W_hh, b_h, W_S, b_S):
    h_t = T.tanh(W_hx[:, x_t] + T.dot(h_t_prev, W_hh) + b_h)
    y_t = T.nnet.softmax(T.dot(W_S, h_t) + b_S)

    return [h_t, y_t]

print >> sys.stderr, "Declaring Theano Loop and Graph"
sys.stderr.flush()

[h, y_out], _ = theano.scan(
    forward_propagation_step,
    sequences=X,
    outputs_info=[dict(initial=T.zeros(n_h)), None],
    non_sequences=[W_hx, W_hh, b_h, W_S, b_S],
    truncate_gradient=50,
    n_steps=X.shape[0]
)

p_y_given_x = y_out[:, 0, :]

y_pred = T.argmax(p_y_given_x, axis=1)

loss = T.nnet.categorical_crossentropy(p_y_given_x, y).mean()

print >> sys.stderr, "Getting gradients"
sys.stderr.flush()
dWhx = T.grad(loss, wrt=W_hx)
dWhh = T.grad(loss, wrt=W_hh)
dbh = T.grad(loss, wrt=b_h)
dWS = T.grad(loss, wrt=W_S)
dbS = T.grad(loss, wrt=b_S)

print >> sys.stderr, "Compiling functions"
forward_propagation = theano.function([X], y_out)
loss_calculation = theano.function([X, y], loss)
predict = theano.function([X], y_pred)
# bbtt = theano.function([X, y], [dWhx, dWhh, dbh, dWS, dbS])

alpha = T.scalar('alpha')

print >> sys.stderr, "Declaring gradients"
sys.stderr.flush()
updates = [
    (W_hx, W_hx - alpha * dWhx),
    (W_hh, W_hh - alpha * dWhh),
    (b_h, b_h - alpha * dbh),
    (W_S, W_S - alpha * dWS),
    (b_S, b_S - alpha * dbS)
]

gradient_step = theano.function(
    inputs=[X, y, alpha],
    outputs=loss,
    updates=updates
)

print >> sys.stderr, "Getting dataset"
sys.stderr.flush()

X_train = []
y_train = []

for char in corpus:
    X_train.append(char)
    y_train.append(char)

X_train = np.array(X_train[:-1]).astype('int32')
y_train = np.array(y_train[1:]).astype('int32')

print >> sys.stderr, "Begin training"
sys.stderr.flush()

for i in xrange(1, 1000):  # We train for epochs times
    for j in xrange(y_train.shape[0], 128):
        gradient_step(X_train[j:j+128], y_train[j:j+128], 0.01)

    if i % 100 == 0:
        print >> sys.stderr, "Loss for iteration {}: {}".format(
            i, loss_calculation(X_train, y_train)
        )
        sys.stderr.flush()

        # Generate a 2000 characters text

        characters_indexes = [np.random.randint(28, 82)]
        # The first character is alphabetic random

        for j in xrange(2000):
            characters_indexes.append(predict(np.hstack(characters_indexes).astype('int32'))[-1])

        print "".join([corpus.characters[char] for char in characters_indexes])
        print
        print "#" * 100
        print
        sys.stdout.flush()

reason: import numpy,import theano


text:
from __future__ import division

from dolo.symbolic.model import SModel
from dolo.symbolic.symbolic import Equation

import numpy as np
import numpy.linalg

class MirFacCompiler(object):

    def __init__(self,model):
        self.model = model
        self.__transformed_model__ = None
        # we assume model has already been checked

    def read_model(self):

        if self.__transformed_model__:
            return self.__transformed_model__

        dmodel = SModel(**self.model) # copy the model
        dmodel.check_consistency(auto_remove_variables=False)

        def_eqs = [eq for eq in dmodel.equations if eq.tags['eq_type'] in ('def', 'auxiliary')]

        from dolo.symbolic.symbolic import map_function_to_expression
        from dolo.symbolic.symbolic import Variable
        def timeshift(v,n):
            if isinstance(v,Variable):
                return v(n)
            else:
                return v

        import sympy

        #### build substitution dict
        def_dict = {}
        for eq in def_eqs:
            v = eq.lhs
            rhs = sympy.sympify( eq.rhs )
            def_dict[v] = rhs
            def_dict[v(1)] = map_function_to_expression( lambda x: timeshift(x,1), rhs)

        new_equations = []
        tbr = []
        for i,eq in enumerate(dmodel.equations) :
            if not ('def' == eq.tags['eq_type']):
                lhs = sympy.sympify( eq.lhs ).subs(def_dict)
                rhs = sympy.sympify( eq.rhs ).subs(def_dict)
                neq = Equation(lhs,rhs).tag(**eq.tags)
                new_equations.append( neq )

        dmodel['equations'] = new_equations
        dmodel.check_consistency()

        f_eqs = [eq for eq in dmodel.equations if eq.tags['eq_type'] in ('f','arbitrage')]
        g_eqs = [eq for eq in dmodel.equations if eq.tags['eq_type'] in ('g','transition')]
        h_eqs = [eq for eq in dmodel.equations if eq.tags['eq_type'] in ('h','expectation')]

        states_vars = [eq.lhs for eq in g_eqs]
        exp_vars =  [eq.lhs for eq in h_eqs]
        controls = set(dmodel.variables) - set(states_vars + exp_vars)
        controls = list(controls)

        states_vars = [v for v in dmodel.variables if v in states_vars]
        exp_vars = [v for v in dmodel.variables if v in exp_vars]
        controls = [v for v in dmodel.variables if v in controls]

        # now we remove the left side of equations
        f_eqs = [eq.gap for eq in f_eqs]
        g_eqs = [eq.rhs for eq in g_eqs]
        h_eqs = [eq.rhs for eq in h_eqs]

        g_eqs = [map_function_to_expression(lambda x: timeshift(x,1),eq) for eq in g_eqs]
        #h_eqs = [map_function_to_expression(lambda x: timeshift(x,-1),eq) for eq in h_eqs] #no

    #    sub_list[v] = v.name

        # read complementarity conditions
        compcond = {}
        of_eqs = [eq for eq in dmodel.equations if eq.tags['eq_type'] in ('f','arbitrage')]
        locals = {}
        import sympy
        locals['inf'] = sympy.Symbol('inf')
        locals['log'] = sympy.log # this should be more generic
        locals['exp'] = sympy.exp

        for v in dmodel.variables + dmodel.parameters:
            locals[v.name] = v
        import re
        compregex = re.compile('(.*)<=(.*)<=(.*)')
        for eq in of_eqs:
            tg = eq.tags['complementarity']
            [lhs,mhs,rhs] = compregex.match(tg).groups()
            [lhs,mhs,rhs] = [dmodel.eval_string(x) for x in [lhs,mhs,rhs] ]
            compcond[mhs] = (lhs,rhs)

        complementarities = [compcond[v] for v in controls]

        inf_bounds = [c[0] for c in complementarities]
        sup_bounds = [c[1] for c in complementarities]

        data = {
            'f_eqs': f_eqs,
            'g_eqs': g_eqs,
            'h_eqs': h_eqs,
            'controls': controls,
            'states_vars': states_vars,
            'exp_vars': exp_vars,
            'inf_bounds': inf_bounds,
            'sup_bounds': sup_bounds
        }

        self.__transformed_model__ = data # cache computation

        return data

    def perturbation_solution(self,with_bounds=False):
        """
        Returns perturbation solution around the steady-state
        Result is [X,Y,Z] representing :
        s_t = X s_{t-1} + Y e_t
        x_t = Z s_t
        """

        model = self.model
        data = self.read_model()

        controls_i = [model.variables.index(v) for v in data['controls'] ]
        states_i = [model.variables.index(v) for v in data['states_vars'] ]

        from dolo.numeric.perturbations import solve_decision_rule
        dr = solve_decision_rule(model,order=1)

        A = dr['g_a'][states_i,:]
        B = dr['g_a'][controls_i,:]

        [Z,err,rank,junk] = np.linalg.lstsq(A.T,B.T)
        Z = Z.T

        X = A[:,states_i] + np.dot( A[:,controls_i], Z )
        Y = dr.ghu[states_i,:]

        # steady_state values
        s_ss = dr['ys'][states_i,]
        x_ss = dr['ys'][controls_i,]

        # We also compute bounds for distribution at a given probability interval
        if with_bounds:
            M = np.linalg.solve( 1-X, Y)

            Sigma = np.array(model.covariances).astype(np.float64)
            [V,P] = np.linalg.eigh(Sigma)
            # we have Sigma == P * diag(V) * P.T
            # unconditional distribution is ( P*diag(V) ) * N
            # where N is the normal distribution associated
            H = np.dot(Y,np.dot(P,np.diag(V)))
            n_s = Sigma.shape[0]
            I = np.eye(n_s)
            points = np.concatenate([H,-H],axis=1)
            lam = 2 # this coefficient should be more cleverly defined
            p_infs = np.min(points,axis = 1) * lam
            p_max  = np.max(points,axis = 1) * lam

            bounds = np.row_stack([
                s_ss + p_infs,
                s_ss + p_max
            ])
        else:
            bounds = None

        return [[s_ss,x_ss],[X,Y,Z],bounds]

    def process_output_recs(self, solution_order=False, fname=None):

        data = self.read_model()
        dmodel = self.model
        model = dmodel

        f_eqs = data['f_eqs']
        g_eqs = data['g_eqs']
        h_eqs = data['h_eqs']
        states_vars = data['states_vars']
        controls = data['controls']
        exp_vars = data['exp_vars']
        inf_bounds = data['inf_bounds']
        sup_bounds = data['sup_bounds']

        controls_f = [v(1) for v in controls]
        states_f = [v(1) for v in states_vars]

        sub_list = dict()
        for i,v in enumerate(exp_vars):
            sub_list[v] = 'z(:,{0})'.format(i+1)

        for i,v in enumerate(controls):
            sub_list[v] = 'x(:,{0})'.format(i+1)
            sub_list[v(1)] = 'xnext(:,{0})'.format(i+1)

        for i,v in enumerate(states_vars):
            sub_list[v] = 's(:,{0})'.format(i+1)
            sub_list[v(1)] = 'snext(:,{0})'.format(i+1)

        for i,v in enumerate(dmodel.shocks):
            sub_list[v] = 'e(:,{0})'.format(i+1)

        for i,v in enumerate(dmodel.parameters):
            sub_list[v] = 'p({0})'.format(i+1)

        import sympy
        sub_list[sympy.Symbol('inf')] = 'Inf'

        text = '''function [out1,out2,out3,out4,out5] = {mfname}(flag,s,x,z,e,snext,xnext,p,out);

output = struct('F',1,'Js',0,'Jx',0,'Jsn',0,'Jxn',0,'Jz',0,'hmult',0);

if nargin == 9
  output                     = catstruct(output,out);
  voidcell                   = cell(1,5);
  [out1,out2,out3,out4,out5] = voidcell{{:}};
else
  if nargout >= 2, output.Js = 1; end
  if nargout >= 3, output.Jx = 1; end
  if nargout >= 4
    if strcmpi(flag, 'f')
      output.Jz = 1;
    else
      output.Jsn = 1;
    end
  end
  if nargout >= 5, output.Jxn = 1; end
end

switch flag

  case 'b';
    n = size(s,1);
{eq_bounds_block}

  case 'f';
    n = size(s,1);
{eq_fun_block}
  case 'g';
    n = size(s,1);
{state_trans_block}
  case 'h';
    n = size(snext,1);
{exp_fun_block}
  case 'e';
    out1 = [];

  case 'model'; % informations about the model
{model_info}

end
'''

        from dolo.compiler.common import DicPrinter

        dp = DicPrinter(sub_list)

        def write_eqs(eq_l,outname='out1',ntabs=0):
            eq_block = '  ' * ntabs + '{0} = zeros(n,{1});'.format(outname,len(eq_l))
            for i,eq in enumerate(eq_l):
                eq_block += '\n' + '  ' * ntabs + '{0}(:,{1}) = {2};'.format( outname,  i+1,  dp.doprint_matlab(eq,vectorize=True) )
            return eq_block

        def write_der_eqs(eq_l,v_l,lhs,ntabs=0):
            eq_block = '  ' * ntabs + '{lhs} = zeros(n,{0},{1});'.format(len(eq_l),len(v_l),lhs=lhs)
            eq_l_d = eqdiff(eq_l,v_l)
            for i,eqq in enumerate(eq_l_d):
                for j,eq in enumerate(eqq):
                    s = dp.doprint_matlab( eq, vectorize=True )
                    eq_block += '\n' + '  ' * ntabs + '{lhs}(:,{0},{1}) = {2}; % d eq_{eq_n} w.r.t. {vname}'.format(i+1,j+1,s,lhs=lhs,eq_n=i+1,vname=str(v_l[j]) )
            return eq_block

        eq_bounds_block = write_eqs(inf_bounds,ntabs=2)
        eq_bounds_block += '\n'
        eq_bounds_block += write_eqs(sup_bounds,'out2',ntabs=2)

        eq_f_block = '''
    % f
    if output.F
{0}
    end

    % df/ds
    if output.Js
{1}
    end

    % df/dx
    if output.Jx
{2}
    end

    % df/dz
    if output.Jz
{3}
    end
        '''.format( write_eqs(f_eqs,'out1',3),
                    write_der_eqs(f_eqs,states_vars,'out2',3),
                    write_der_eqs(f_eqs,controls,'out3',3),
                    write_der_eqs(f_eqs,exp_vars,'out4',3)
            )

        eq_g_block = '''
    % g
    if output.F
{0}
    end

    if output.Js
{1}
    end

    if output.Jx
{2}
    end
        '''.format( write_eqs(g_eqs,'out1',3),
                    write_der_eqs(g_eqs,states_vars,'out2',3),
                    write_der_eqs(g_eqs,controls,'out3',3)
            )

        eq_h_block = '''
    %h
    if output.F
{0}
    end

    if output.Js
{1}
    end

    if output.Jx
{2}
    end

    if output.Jsn
{3}
    end

    if output.Jxn
{4}
    end
        '''.format(
             write_eqs(h_eqs,'out1',3),
             write_der_eqs(h_eqs,states_vars,'out2',3),
             write_der_eqs(h_eqs,controls,'out3',3),
             write_der_eqs(h_eqs,states_f,'out4',3),
             write_der_eqs(h_eqs,controls_f,'out5',3)
        )

        # if not with_param_names:
        #    eq_h_block = 's=snext;\nx=xnext;\n'+eq_h_block

        # param_def = 'p = [ ' + str.join(',',[p.name for p in dmodel.parameters])  + '];'

        from dolo.misc.matlab import value_to_mat

        # read model informations
        [y,x,params_values] = model.read_calibration()
        #params_values = '[' + str.join(  ',', [ str( p ) for p in params] ) + '];'
        vvs = model.variables
        s_ss = [ y[vvs.index(v)] for v in model['variables_groups']['states'] ]
        x_ss = [ y[vvs.index(v)] for v in model['variables_groups']['controls'] ]

        model_info = '''
    mod = struct;
    mod.s_ss = {s_ss};
    mod.x_ss = {x_ss};
    mod.params = {params_values};
'''.format(
    s_ss = value_to_mat(s_ss),
    x_ss = value_to_mat(x_ss),
    params_values = value_to_mat(params_values)
)
        if solution_order:
            from dolo.numeric.perturbations_to_states import approximate_controls

            ZZ = approximate_controls(self.model,order=solution_order)
            n_c = len(controls)

            ZZ = [np.array(e) for e in ZZ]
            ZZ = [e[:n_c,...] for e in ZZ] # keep only control vars. (x) not expectations (h)

            solution = "    mod.X = cell({0},1);\n".format(len(ZZ))
            for i,zz in enumerate(ZZ):
                solution += "    mod.X{{{0}}} = {1};\n".format(i+1,value_to_mat(zz))
            model_info += solution
        model_info += '    out1 = mod;\n'

        text = text.format(
            eq_bounds_block = eq_bounds_block,
            mfname = fname if fname else 'mf_' + model.fname,
            eq_fun_block=eq_f_block,
            state_trans_block=eq_g_block,
            exp_fun_block=eq_h_block,
#            solution = solution,
            model_info = model_info
        )

        return text

def eqdiff(leq,lvars):
    resp = []
    for eq in leq:
        el = [ eq.diff(v) for v in lvars]
        resp += [el]
    return resp

if __name__ == "__main__":
    from dolo import dynare_import
    model = dynare_import('../../../examples/global_models/optimal_growth.mod')
    model.check()
    mfcp = MirFacCompiler(model)
    print mfcp.process_output_recs()
    #print mfcp.process_output_python()

max_stars_repo_name: christophe-gouel/dolo

max_issues_repo_name: christophe-gouel/dolo

max_forks_repo_name: christophe-gouel/dolo

reason: import numpy,import sympy


text:
r"""
=======================================
Robust vs Empirical covariance estimate
=======================================

The usual covariance maximum likelihood estimate is very sensitive to the
presence of outliers in the data set. In such a case, it would be better to
use a robust estimator of covariance to guarantee that the estimation is
resistant to "erroneous" observations in the data set. [1]_, [2]_

Minimum Covariance Determinant Estimator
----------------------------------------
The Minimum Covariance Determinant estimator is a robust, high-breakdown point
(i.e. it can be used to estimate the covariance matrix of highly contaminated
datasets, up to
:math:`\frac{n_\text{samples} - n_\text{features}-1}{2}` outliers) estimator of
covariance. The idea is to find
:math:`\frac{n_\text{samples} + n_\text{features}+1}{2}`
observations whose empirical covariance has the smallest determinant, yielding
a "pure" subset of observations from which to compute standards estimates of
location and covariance. After a correction step aiming at compensating the
fact that the estimates were learned from only a portion of the initial data,
we end up with robust estimates of the data set location and covariance.

The Minimum Covariance Determinant estimator (MCD) has been introduced by
P.J.Rousseuw in [3]_.

Evaluation
----------
In this example, we compare the estimation errors that are made when using
various types of location and covariance estimates on contaminated Gaussian
distributed data sets:

- The mean and the empirical covariance of the full dataset, which break
  down as soon as there are outliers in the data set
- The robust MCD, that has a low error provided
  :math:`n_\text{samples} > 5n_\text{features}`
- The mean and the empirical covariance of the observations that are known
  to be good ones. This can be considered as a "perfect" MCD estimation,
  so one can trust our implementation by comparing to this case.

References
----------
.. [1] Johanna Hardin, David M Rocke. The distribution of robust distances.
    Journal of Computational and Graphical Statistics. December 1, 2005,
    14(4): 928-946.
.. [2] Zoubir A., Koivunen V., Chakhchoukh Y. and Muma M. (2012). Robust
    estimation in signal processing: A tutorial-style treatment of
    fundamental concepts. IEEE Signal Processing Magazine 29(4), 61-80.
.. [3] P. J. Rousseeuw. Least median of squares regression. Journal of American
    Statistical Ass., 79:871, 1984.

"""
print(__doc__)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.font_manager

from sklearn.covariance import EmpiricalCovariance, MinCovDet

# example settings
n_samples = 80
n_features = 5
repeat = 10

range_n_outliers = np.concatenate(
    (np.linspace(0, n_samples / 8, 5),
     np.linspace(n_samples / 8, n_samples / 2, 5)[1:-1])).astype(int)

# definition of arrays to store results
err_loc_mcd = np.zeros((range_n_outliers.size, repeat))
err_cov_mcd = np.zeros((range_n_outliers.size, repeat))
err_loc_emp_full = np.zeros((range_n_outliers.size, repeat))
err_cov_emp_full = np.zeros((range_n_outliers.size, repeat))
err_loc_emp_pure = np.zeros((range_n_outliers.size, repeat))
err_cov_emp_pure = np.zeros((range_n_outliers.size, repeat))

# computation
for i, n_outliers in enumerate(range_n_outliers):
    for j in range(repeat):

        rng = np.random.RandomState(i * j)

        # generate data
        X = rng.randn(n_samples, n_features)
        # add some outliers
        outliers_index = rng.permutation(n_samples)[:n_outliers]
        outliers_offset = 10. * \
            (np.random.randint(2, size=(n_outliers, n_features)) - 0.5)
        X[outliers_index] += outliers_offset
        inliers_mask = np.ones(n_samples).astype(bool)
        inliers_mask[outliers_index] = False

        # fit a Minimum Covariance Determinant (MCD) robust estimator to data
        mcd = MinCovDet().fit(X)
        # compare raw robust estimates with the true location and covariance
        err_loc_mcd[i, j] = np.sum(mcd.location_ ** 2)
        err_cov_mcd[i, j] = mcd.error_norm(np.eye(n_features))

        # compare estimators learned from the full data set with true
        # parameters
        err_loc_emp_full[i, j] = np.sum(X.mean(0) ** 2)
        err_cov_emp_full[i, j] = EmpiricalCovariance().fit(X).error_norm(
            np.eye(n_features))

        # compare with an empirical covariance learned from a pure data set
        # (i.e. "perfect" mcd)
        pure_X = X[inliers_mask]
        pure_location = pure_X.mean(0)
        pure_emp_cov = EmpiricalCovariance().fit(pure_X)
        err_loc_emp_pure[i, j] = np.sum(pure_location ** 2)
        err_cov_emp_pure[i, j] = pure_emp_cov.error_norm(np.eye(n_features))

# Display results
font_prop = matplotlib.font_manager.FontProperties(size=11)
plt.subplot(2, 1, 1)
lw = 2
plt.errorbar(range_n_outliers, err_loc_mcd.mean(1),
             yerr=err_loc_mcd.std(1) / np.sqrt(repeat),
             label="Robust location", lw=lw, color='m')
plt.errorbar(range_n_outliers, err_loc_emp_full.mean(1),
             yerr=err_loc_emp_full.std(1) / np.sqrt(repeat),
             label="Full data set mean", lw=lw, color='green')
plt.errorbar(range_n_outliers, err_loc_emp_pure.mean(1),
             yerr=err_loc_emp_pure.std(1) / np.sqrt(repeat),
             label="Pure data set mean", lw=lw, color='black')
plt.title("Influence of outliers on the location estimation")
plt.ylabel(r"Error ($||\mu - \hat{\mu}||_2^2$)")
plt.legend(loc="upper left", prop=font_prop)

plt.subplot(2, 1, 2)
x_size = range_n_outliers.size
plt.errorbar(range_n_outliers, err_cov_mcd.mean(1),
             yerr=err_cov_mcd.std(1),
             label="Robust covariance (mcd)", color='m')
plt.errorbar(range_n_outliers[:(x_size // 5 + 1)],
             err_cov_emp_full.mean(1)[:(x_size // 5 + 1)],
             yerr=err_cov_emp_full.std(1)[:(x_size // 5 + 1)],
             label="Full data set empirical covariance", color='green')
plt.plot(range_n_outliers[(x_size // 5):(x_size // 2 - 1)],
         err_cov_emp_full.mean(1)[(x_size // 5):(x_size // 2 - 1)],
         color='green', ls='--')
plt.errorbar(range_n_outliers, err_cov_emp_pure.mean(1),
             yerr=err_cov_emp_pure.std(1),
             label="Pure data set empirical covariance", color='black')
plt.title("Influence of outliers on the covariance estimation")
plt.xlabel("Amount of contamination (%)")
plt.ylabel("RMSE")
plt.legend(loc="upper center", prop=font_prop)

#plt.show()