text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
df = pd.read_csv('logit-data.txt', header=None)
df.head()

X = df.iloc[:, :-1]
y = df.iloc[:, -1]

#Add a bias column to the X. The value of the bias column is usually one.

X = np.c_[np.ones((X.shape[0], 1)), X]
X[:10]

#Here, our X is a two-dimensional array and y is a one-dimensional array. Let’s make the ‘y’ two-dimensional to match the dimensions.

y = y[:, np.newaxis]
y[:10]

def sigmoid(x, theta):
    z= np.dot(x, theta)
    return 1/(1+np.exp(-z))

def hypothesis(theta, x):
    return sigmoid(x, theta)

def cost_function(theta, x, y):
    m = X.shape[0]
    h = hypothesis(theta, x)
    return -(1/m)*np.sum(y*np.log(h) + (1-y)*np.log(1-h))

def gradient(theta, x, y):
    m = X.shape[0]
    h = hypothesis(theta, x)
    return (1/m) * np.dot(X.T, (h-y))

theta = np.zeros((X.shape[1], 1))
from scipy.optimize import minimize,fmin_tnc
def fit(x, y, theta):
    opt_weights = fmin_tnc(func=cost_function, x0=theta, fprime=gradient, args=(x, y.flatten()))
    return opt_weights[0]
parameters = fit(X, y, theta)

h = hypothesis(parameters, X)

def predict(h):
    h1 = []
    for i in h:
        if i>=0.5:
            h1.append(1)
        else:
            h1.append(0)
    return h1
y_pred = predict(h)

accuracy = 0
for i in range(0, len(y_pred)):
    if y_pred[i] == y[i]:
        accuracy += 1
accuracy/len(y)
print(accuracy)

max_stars_repo_path: Logistic Regression - Gradient Descent.py

max_stars_repo_name: karthiksekaran/mlprojects

max_issues_repo_name: karthiksekaran/mlprojects

max_forks_repo_name: karthiksekaran/mlprojects


text:
#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
# by the software owners: The Regents of the University of California, through
# Lawrence Berkeley National Laboratory,  National Technology & Engineering
# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
# Research Corporation, et al.  All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
# license information.
#################################################################################
"""
Tests for Flash unit model.
Author: Jaffer Ghouse
"""
import pytest
from pyomo.environ import (ConcreteModel,
                           TerminationCondition,
                           SolverStatus,
                           value,
                           units,
                           Var)
from pyomo.util.check_units import (assert_units_consistent,
                                    assert_units_equivalent)

from idaes.core import (FlowsheetBlock, MaterialBalanceType, EnergyBalanceType,
                        MomentumBalanceType)
from idaes.generic_models.unit_models.flash import Flash, EnergySplittingType
from idaes.generic_models.properties.activity_coeff_models.BTX_activity_coeff_VLE \
    import BTXParameterBlock
from idaes.generic_models.properties import iapws95
from idaes.core.util.model_statistics import (degrees_of_freedom,
                                              number_variables,
                                              number_total_constraints,
                                              number_unused_variables)
from idaes.core.util.testing import (PhysicalParameterTestBlock,
                                     initialization_tester)
from idaes.core.util import get_solver
import idaes.core.util.scaling as iscale

# -----------------------------------------------------------------------------
# Get default solver for testing
solver = get_solver()

# -----------------------------------------------------------------------------
@pytest.mark.unit
def test_config():
    m = ConcreteModel()
    m.fs = FlowsheetBlock(default={"dynamic": False})

    m.fs.properties = PhysicalParameterTestBlock()

    m.fs.unit = Flash(default={"property_package": m.fs.properties})

    # Check unit config arguments
    assert len(m.fs.unit.config) == 11

    assert not m.fs.unit.config.dynamic
    assert not m.fs.unit.config.has_holdup
    assert m.fs.unit.config.material_balance_type == \
        MaterialBalanceType.useDefault
    assert m.fs.unit.config.energy_balance_type == \
        EnergyBalanceType.useDefault
    assert m.fs.unit.config.momentum_balance_type == \
        MomentumBalanceType.pressureTotal
    assert m.fs.unit.config.ideal_separation
    assert m.fs.unit.config.has_heat_transfer
    assert m.fs.unit.config.has_pressure_change
    assert m.fs.unit.config.property_package is m.fs.properties
    assert m.fs.unit.config.energy_split_basis == \
        EnergySplittingType.equal_temperature

# -----------------------------------------------------------------------------
@pytest.mark.unit
def test_calc_scale():
    m = ConcreteModel()
    m.fs = FlowsheetBlock(default={"dynamic": False})
    m.fs.properties = PhysicalParameterTestBlock()
    m.fs.unit = Flash(default={"property_package": m.fs.properties})
    iscale.calculate_scaling_factors(m)

# -----------------------------------------------------------------------------
class TestBTXIdeal(object):
    @pytest.fixture(scope="class")
    def btx(self):
        m = ConcreteModel()
        m.fs = FlowsheetBlock(default={"dynamic": False})

        m.fs.properties = BTXParameterBlock(default={"valid_phase":
                                                     ('Liq', 'Vap'),
                                                     "activity_coeff_model":
                                                     "Ideal"})

        m.fs.unit = Flash(default={"property_package": m.fs.properties})

        m.fs.unit.inlet.flow_mol.fix(1)
        m.fs.unit.inlet.temperature.fix(368)
        m.fs.unit.inlet.pressure.fix(101325)
        m.fs.unit.inlet.mole_frac_comp[0, "benzene"].fix(0.5)
        m.fs.unit.inlet.mole_frac_comp[0, "toluene"].fix(0.5)

        m.fs.unit.heat_duty.fix(0)
        m.fs.unit.deltaP.fix(0)

        return m

    @pytest.mark.build
    @pytest.mark.unit
    def test_build(self, btx):
        assert hasattr(btx.fs.unit.inlet, "flow_mol")
        assert hasattr(btx.fs.unit.inlet, "mole_frac_comp")
        assert hasattr(btx.fs.unit.inlet, "temperature")
        assert hasattr(btx.fs.unit.inlet, "pressure")

        assert hasattr(btx.fs.unit, "vap_outlet")
        assert len(btx.fs.unit.vap_outlet.vars) == 4
        assert hasattr(btx.fs.unit.vap_outlet, "flow_mol")
        assert hasattr(btx.fs.unit.vap_outlet, "mole_frac_comp")
        assert hasattr(btx.fs.unit.vap_outlet, "temperature")
        assert hasattr(btx.fs.unit.vap_outlet, "pressure")

        assert hasattr(btx.fs.unit, "liq_outlet")
        assert len(btx.fs.unit.liq_outlet.vars) == 4
        assert hasattr(btx.fs.unit.liq_outlet, "flow_mol")
        assert hasattr(btx.fs.unit.liq_outlet, "mole_frac_comp")
        assert hasattr(btx.fs.unit.liq_outlet, "temperature")
        assert hasattr(btx.fs.unit.liq_outlet, "pressure")

        assert hasattr(btx.fs.unit, "split")
        assert hasattr(btx.fs.unit, "heat_duty")
        assert hasattr(btx.fs.unit, "deltaP")

        assert number_variables(btx) == 48
        assert number_total_constraints(btx) == 41
        assert number_unused_variables(btx) == 0

    @pytest.mark.component
    def test_units(self, btx):
        assert_units_consistent(btx)
        assert_units_equivalent(btx.fs.unit.heat_duty[0], units.W)
        assert_units_equivalent(btx.fs.unit.deltaP[0], units.Pa)

    @pytest.mark.unit
    def test_dof(self, btx):
        assert degrees_of_freedom(btx) == 0

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_initialize(self, btx):
        initialization_tester(btx)

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solve(self, btx):
        results = solver.solve(btx)

        # Check for optimal solution
        assert results.solver.termination_condition == \
            TerminationCondition.optimal
        assert results.solver.status == SolverStatus.ok

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solution(self, btx):
        assert (pytest.approx(0.603, abs=1e-3) ==
                value(btx.fs.unit.liq_outlet.flow_mol[0]))
        assert (pytest.approx(0.396, abs=1e-3) ==
                value(btx.fs.unit.vap_outlet.flow_mol[0]))
        assert (pytest.approx(368, abs=1e-3) ==
                value(btx.fs.unit.liq_outlet.temperature[0]))
        assert (pytest.approx(101325, abs=1e-3) ==
                value(btx.fs.unit.liq_outlet.pressure[0]))

        assert (pytest.approx(0.412, abs=1e-3) ==
                value(btx.fs.unit.liq_outlet.mole_frac_comp[0, "benzene"]))
        assert (pytest.approx(0.588, abs=1e-3) ==
                value(btx.fs.unit.liq_outlet.mole_frac_comp[0, "toluene"]))
        assert (pytest.approx(0.634, abs=1e-3) ==
                value(btx.fs.unit.vap_outlet.mole_frac_comp[0, "benzene"]))
        assert (pytest.approx(0.366, abs=1e-3) ==
                value(btx.fs.unit.vap_outlet.mole_frac_comp[0, "toluene"]))

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_conservation(self, btx):
        assert abs(value(btx.fs.unit.inlet.flow_mol[0] -
                         (btx.fs.unit.vap_outlet.flow_mol[0] +
                          btx.fs.unit.liq_outlet.flow_mol[0]))) <= 1e-6

    @pytest.mark.ui
    @pytest.mark.unit
    def test_report(self, btx):
        btx.fs.unit.report()

# -----------------------------------------------------------------------------
@pytest.mark.iapws
@pytest.mark.skipif(not iapws95.iapws95_available(),
                    reason="IAPWS not available")
class TestIAPWS(object):
    @pytest.fixture(scope="class")
    def iapws(self):
        m = ConcreteModel()
        m.fs = FlowsheetBlock(default={"dynamic": False})

        m.fs.properties = iapws95.Iapws95ParameterBlock(default={
                "phase_presentation": iapws95.PhaseType.LG})

        m.fs.unit = Flash(default={
                "property_package": m.fs.properties,
                "ideal_separation": False,
                "energy_split_basis": EnergySplittingType.enthalpy_split})

        m.fs.unit.inlet.flow_mol.fix(100)
        m.fs.unit.inlet.enth_mol.fix(24000)
        m.fs.unit.inlet.pressure.fix(101325)

        m.fs.unit.heat_duty.fix(0)
        m.fs.unit.deltaP.fix(0)

        return m

    @pytest.mark.build
    @pytest.mark.unit
    def test_build(self, iapws):
        assert len(iapws.fs.unit.inlet.vars) == 3
        assert hasattr(iapws.fs.unit.inlet, "flow_mol")
        assert hasattr(iapws.fs.unit.inlet, "enth_mol")
        assert hasattr(iapws.fs.unit.inlet, "pressure")

        assert hasattr(iapws.fs.unit, "vap_outlet")
        assert len(iapws.fs.unit.vap_outlet.vars) == 3
        assert hasattr(iapws.fs.unit.vap_outlet, "flow_mol")
        assert hasattr(iapws.fs.unit.vap_outlet, "enth_mol")
        assert hasattr(iapws.fs.unit.vap_outlet, "pressure")

        assert hasattr(iapws.fs.unit, "liq_outlet")
        assert len(iapws.fs.unit.liq_outlet.vars) == 3
        assert hasattr(iapws.fs.unit.liq_outlet, "flow_mol")
        assert hasattr(iapws.fs.unit.liq_outlet, "enth_mol")
        assert hasattr(iapws.fs.unit.liq_outlet, "pressure")

        assert hasattr(iapws.fs.unit, "split")
        assert hasattr(iapws.fs.unit, "heat_duty")
        assert hasattr(iapws.fs.unit, "deltaP")

        assert number_variables(iapws) == 18
        assert number_total_constraints(iapws) == 13
        assert number_unused_variables(iapws) == 0

    @pytest.mark.component
    def test_units(self, iapws):
        assert_units_consistent(iapws)
        assert_units_equivalent(iapws.fs.unit.heat_duty[0], units.W)
        assert_units_equivalent(iapws.fs.unit.deltaP[0], units.Pa)

    @pytest.mark.unit
    def test_dof(self, iapws):
        assert degrees_of_freedom(iapws) == 0

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_initialize(self, iapws):
        initialization_tester(iapws)

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solve(self, iapws):
        results = solver.solve(iapws)

        # Check for optimal solution
        assert results.solver.termination_condition == \
            TerminationCondition.optimal
        assert results.solver.status == SolverStatus.ok

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solution(self, iapws):
        t_in = value(iapws.fs.unit.control_volume.properties_in[0].temperature)
        assert (pytest.approx(373.12, abs=1e-2) == t_in)
        assert (pytest.approx(t_in, abs=1e-2) == value(
                iapws.fs.unit.control_volume.properties_out[0].temperature))
        assert (pytest.approx(t_in, abs=1e-2) == value(
                iapws.fs.unit.split.Vap_state[0].temperature))
        assert (pytest.approx(t_in, abs=1e-2) == value(
                iapws.fs.unit.split.Liq_state[0].temperature))

        assert (pytest.approx(101325.0, abs=1e3) ==
                value(iapws.fs.unit.vap_outlet.pressure[0]))
        assert (pytest.approx(48200, abs=1e3) ==
                value(iapws.fs.unit.vap_outlet.enth_mol[0]))
        assert (pytest.approx(40.467, abs=1e-3) ==
                value(iapws.fs.unit.vap_outlet.flow_mol[0]))

        assert (pytest.approx(101325.0, abs=1e3) ==
                value(iapws.fs.unit.liq_outlet.pressure[0]))
        assert (pytest.approx(7549.4, abs=1) ==
                value(iapws.fs.unit.liq_outlet.enth_mol[0]))
        assert (pytest.approx(59.532, abs=1e-3) ==
                value(iapws.fs.unit.liq_outlet.flow_mol[0]))

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_conservation(self, iapws):
        assert abs(value(iapws.fs.unit.inlet.flow_mol[0] -
                         (iapws.fs.unit.vap_outlet.flow_mol[0] +
                          iapws.fs.unit.liq_outlet.flow_mol[0]))) <= 1e-6
        assert abs(value(iapws.fs.unit.inlet.flow_mol[0] *
                         iapws.fs.unit.inlet.enth_mol[0] -
                         (iapws.fs.unit.vap_outlet.flow_mol[0] *
                          iapws.fs.unit.vap_outlet.enth_mol[0] +
                          iapws.fs.unit.liq_outlet.flow_mol[0] *
                          iapws.fs.unit.liq_outlet.enth_mol[0]))) <= 1e-6

    @pytest.mark.ui
    @pytest.mark.component
    def test_report(self, iapws):
        iapws.fs.unit.report()

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_costing(self, iapws):
        iapws.fs.unit.get_costing()
        assert isinstance(iapws.fs.unit.costing.purchase_cost, Var)
        iapws.fs.unit.diameter.fix(2)
        iapws.fs.unit.length.fix(4)
        # initialize unit with costing block
        iapws.fs.unit.initialize()
        # check costing initialized correct
        assert (pytest.approx(86957.195, abs=1e-3) ==
                value(iapws.fs.unit.costing.purchase_cost))

        results = solver.solve(iapws)
        # Check for optimal solution
        assert results.solver.termination_condition == \
            TerminationCondition.optimal
        assert results.solver.status == SolverStatus.ok
        assert (pytest.approx(63787.06525, abs=1e3) ==
                value(iapws.fs.unit.costing.base_cost))
        assert (pytest.approx(97660.6169, abs=1e3) ==
                value(iapws.fs.unit.costing.purchase_cost))

        assert_units_consistent(iapws.fs.unit.costing)

ext: py

max_stars_repo_name: carldlaird/idaes-pse

max_issues_repo_name: carldlaird/idaes-pse

max_forks_repo_name: carldlaird/idaes-pse

reason: from pyomo


text:
# Copyright 2021 DeepMind Technologies Limited
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Losses for TandemDQN."""

from typing import Any, Callable

import chex
import jax
import jax.numpy as jnp
import rlax

from tandem_dqn import networks

# Batch variants of double_q_learning and SARSA.
batch_double_q_learning = jax.vmap(rlax.double_q_learning)
batch_sarsa_learning = jax.vmap(rlax.sarsa)

# Batch variant of quantile_q_learning with fixed tau input across batch.
batch_quantile_q_learning = jax.vmap(
    rlax.quantile_q_learning, in_axes=(0, None, 0, 0, 0, 0, 0, None))

def _mc_learning(
    q_tm1: chex.Array,
    a_tm1: chex.Numeric,
    mc_return_tm1: chex.Array,
) -> chex.Numeric:
  """Calculates the MC return error."""
  chex.assert_rank([q_tm1, a_tm1], [1, 0])
  chex.assert_type([q_tm1, a_tm1], [float, int])
  return mc_return_tm1 - q_tm1[a_tm1]

# Batch variant of MC learning.
batch_mc_learning = jax.vmap(_mc_learning)

def _qr_loss(q_tm1, q_t, q_target_t, transitions, rng_key):
  """Calculates QR-Learning loss from network outputs and transitions."""
  del q_t, rng_key  # Unused.
  # Compute Q value distributions.
  huber_param = 1.
  quantiles = networks.make_quantiles()
  losses = batch_quantile_q_learning(
      q_tm1.q_dist,
      quantiles,
      transitions.a_tm1,
      transitions.r_t,
      transitions.discount_t,
      q_target_t.q_dist,  # No double Q-learning here.
      q_target_t.q_dist,
      huber_param,
  )
  loss = jnp.mean(losses)
  return loss

def _sarsa_loss(q_tm1, q_t, transitions, rng_key):
  """Calculates SARSA loss from network outputs and transitions."""
  del rng_key  # Unused.
  grad_error_bound = 1. / 32
  td_errors = batch_sarsa_learning(
      q_tm1.q_values,
      transitions.a_tm1,
      transitions.r_t,
      transitions.discount_t,
      q_t.q_values,
      transitions.a_t
  )
  td_errors = rlax.clip_gradient(td_errors, -grad_error_bound, grad_error_bound)
  losses = rlax.l2_loss(td_errors)
  loss = jnp.mean(losses)
  return loss

def _mc_loss(q_tm1, transitions, rng_key):
  """Calculates Monte-Carlo return loss, i.e. regression towards MC return."""
  del rng_key  # Unused.
  errors = batch_mc_learning(q_tm1.q_values, transitions.a_tm1,
                             transitions.mc_return_tm1)
  loss = jnp.mean(rlax.l2_loss(errors))
  return loss

def _double_q_loss(q_tm1, q_t, q_target_t, transitions, rng_key):
  """Calculates Double Q-Learning loss from network outputs and transitions."""
  del rng_key  # Unused.
  grad_error_bound = 1. / 32
  td_errors = batch_double_q_learning(
      q_tm1.q_values,
      transitions.a_tm1,
      transitions.r_t,
      transitions.discount_t,
      q_target_t.q_values,
      q_t.q_values,
  )
  td_errors = rlax.clip_gradient(td_errors, -grad_error_bound, grad_error_bound)
  losses = rlax.l2_loss(td_errors)
  loss = jnp.mean(losses)
  return loss

def _q_regression_loss(q_tm1, q_tm1_target):
  """Loss for regression of all action values towards targets."""
  errors = q_tm1.q_values - jax.lax.stop_gradient(q_tm1_target.q_values)
  loss = jnp.mean(rlax.l2_loss(errors))
  return loss

def make_loss_fn(loss_type: str, active: bool) -> Callable[..., Any]:
  """Create active or passive loss function of given type."""

  if active:
    primary = lambda x: x.active
    secondary = lambda x: x.passive
  else:
    primary = lambda x: x.passive
    secondary = lambda x: x.active

  def sarsa_loss_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """SARSA loss using own networks."""
    del q_t  # Unused.
    return _sarsa_loss(primary(q_tm1), primary(q_target_t), transitions,
                       rng_key)

  def mc_loss_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """MonteCarlo loss."""
    del q_t, q_target_t
    return _mc_loss(primary(q_tm1), transitions, rng_key)

  def double_q_loss_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """Regular DoubleQ loss using own networks."""
    return _double_q_loss(primary(q_tm1), primary(q_t), primary(q_target_t),
                          transitions, rng_key)

  def double_q_loss_v_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """DoubleQ loss using other network's (target) value function."""
    return _double_q_loss(primary(q_tm1), primary(q_t), secondary(q_target_t),
                          transitions, rng_key)

  def double_q_loss_p_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """DoubleQ loss using other network's (online) argmax policy."""
    return _double_q_loss(primary(q_tm1), secondary(q_t), primary(q_target_t),
                          transitions, rng_key)

  def double_q_loss_pv_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """DoubleQ loss using other network's argmax policy & target value fn."""
    return _double_q_loss(primary(q_tm1), secondary(q_t), secondary(q_target_t),
                          transitions, rng_key)

  # Pure regression.
  def q_regression_loss_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """Pure regression of q_tm1(self) towards q_tm1(other)."""
    del q_t, q_target_t, transitions, rng_key   # Unused.
    return _q_regression_loss(primary(q_tm1), secondary(q_tm1))

  # QR loss.
  def qr_loss_fn(q_tm1, q_t, q_target_t, transitions, rng_key):
    """QR-Q loss using own networks."""
    return _qr_loss(primary(q_tm1), primary(q_t), primary(q_target_t),
                    transitions, rng_key)

  if loss_type == 'double_q':
    return double_q_loss_fn
  elif loss_type == 'sarsa':
    return sarsa_loss_fn
  elif loss_type == 'mc_return':
    return mc_loss_fn
  elif loss_type == 'double_q_v':
    return double_q_loss_v_fn
  elif loss_type == 'double_q_p':
    return double_q_loss_p_fn
  elif loss_type == 'double_q_pv':
    return double_q_loss_pv_fn
  elif loss_type == 'q_regression':
    return q_regression_loss_fn
  elif loss_type == 'qr':
    return qr_loss_fn
  else:
    raise ValueError('Unknown loss "{}"'.format(loss_type))

max_stars_repo_name: feifzhou/deepmind-research

reason: import jax


text:
from os.path import join, exists
from time import time

import numpy as np
from pygermanet import load_germanet, Synset
from tqdm import tqdm

from constants import LDA_PATH
from evaluate_topics import parse_args
from utils import load, init_logging, log_args

np.set_printoptions(precision=3)
gn = load_germanet()
tqdm.pandas()

def orth(synset):
    return synset.lemmas[0].orthForm

def compare_synset_lists(synset_list1, synset_list2, sim_func, agg_func):
    try:
        return agg_func(sim_func(ss1, ss2) for ss1 in synset_list1 for ss2 in synset_list2)
    except ValueError:
        return np.nan

def similarities(topic, topn, ignore_unknown=True, sim_func=Synset.sim_lch, agg_func=max):
    arr = np.zeros((topn, topn))
    for j, ssl1 in enumerate(topic.values):
        for k, ssl2 in enumerate(topic.values[j+1:], j+1):
            arr[j, k] = compare_synset_lists(ssl1, ssl2, sim_func, agg_func)
    arr = np.add(arr, arr.T)
    if ignore_unknown:
        arr[arr == 0] = np.nan
    return np.nanmean(arr)

def main():
    (
        dataset, version, params, nbtopics, topn, cores, corpus_type,
        use_coherence, use_w2v, rerank, lsi, args
    ) = parse_args()

    logger = init_logging(name=f'Eval_topics_on_germanet_{dataset}', basic=False, to_stdout=True, to_file=True)
    log_args(logger, args)
    logg = logger.info

    purpose = 'rerank' if rerank else 'topics'
    topics = load(purpose, dataset, version, corpus_type, lsi, *params, *nbtopics)
    if topn > 0:
        topics = topics[:topn]
    else:
        topn = topics.shape[1]
    logg(f'Number of topics {topics.shape}')

    logg('Getting SynSets for topic terms')
    sstopics = topics.applymap(gn.synsets)

    topics['lch'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.sim_lch, agg_func=max, topn=topn
    )
    topics['lch_ignr_unkwn'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.sim_lch, agg_func=max, topn=topn, ignore_unknown=False
    )
    topics['res'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.sim_res, agg_func=max, topn=topn
    )
    topics['res_ignr_unkwn'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.sim_res, agg_func=max, topn=topn, ignore_unknown=False
    )
    topics['jcn'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.dist_jcn, agg_func=min, topn=topn
    )
    topics['jcn_ignr_unkwn'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.dist_jcn, agg_func=min, topn=topn, ignore_unknown=False
    )
    topics['lin'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.sim_lin, agg_func=max, topn=topn
    )
    topics['lin_ignr_unkwn'] = sstopics.progress_apply(
        similarities, axis=1, sim_func=Synset.sim_lin, agg_func=max, topn=topn, ignore_unknown=False
    )

    topics = topics.iloc[:, topn:]
    tpx_path = join(LDA_PATH, version, 'bow', 'topics')
    if rerank:
        file = join(tpx_path, f'{dataset}_reranker-eval_germanet.csv')
    else:
        file = join(
            tpx_path,
            f'{dataset}{"_"+lsi if lsi else ""}_{version}_{corpus_type}_topic-scores_germanet.csv'
        )
    if exists(file):
        file = file.replace('.csv', f'_{str(time()).split(".")[0]}.csv')

    logg(f'Writing {file}')
    topics.to_csv(file)
    logg('done')

if __name__ == '__main__':
    main()

max_stars_repo_name: andifunke/topic-labeling

max_issues_repo_name: andifunke/topic-labeling

max_forks_repo_name: andifunke/topic-labeling

reason: import numpy


text:
"""
This function is intended to wrap the rewards returned by the CityLearn RL environment, and is meant to
be modified by the participants of The CityLearn Challenge.
CityLearn returns the energy consumption of each building as a reward.
This reward_function takes all the electrical demands of all the buildings and turns them into one or
 multiple rewards for the agent(s)

The current code computes a virtual (made-up) electricity price proportional to the total demand for
 electricity in the neighborhood, and multiplies every
reward by that price. Then it returns the new rewards, which should be used by the agent. Participants
 of the CityLearn Challenge are encouraged to completely modify this function
in order to minimize the 5 proposed metrics.
"""
import numpy as np
from utils.standardization import STATIC_LOAD_IDX, SOLAR_POWER_IDX, mean_normalize

# hyper-parameters
NSL_IDX = STATIC_LOAD_IDX
SGEN_IDX = SOLAR_POWER_IDX

# TODO: Renew reward procesing function
def reward_function(rewards, states, alpha=1, beta=1,
                    total_energy_window=12, heat_energy_window=12,
                    price_factor=0.01, ramping_window=6):
    """
    :param rewards: np.array (*, Building, Time),  Negative of building_energy_demand
    :param states: np.array (*, Building, Time, Raw_State),  Raw States of Related Rewards
    :param alpha: float, weight for smoothness of total demand reward
    :param beta: float, weight for heat demand reward
    :param total_energy_window: int, Length of slide window mean for total demand
    :param heat_energy_window: int, Length of slide window mean for heat demand
    :param price_factor: float, weight for global price
    :param ramping_window: int, Length of slide window mean for smoothness of total demand
    :return:
        final_rewards - nd.array (*, Building, 1), processed current reward
        reward_logger - list of min-max value of smoothness reward, heat demand and final reward
    """
    """
    compute demand for DHW & cooling
    rewards[i] = - (cooling_demand[i] + DHW_demand[i] - solar_gen[i] + no_shiftable_load[i])
                = - cooling_demand[i] - DHW_demand[i] + solar_gen[i] - no_shiftable_load[i]
    no_shiftable_load[i] = states[i][23], solar_gen[i] = states[i][24]
    """

    # compute the price
    heat_len = heat_energy_window
    tot_len = total_energy_window
    assert states.shape[:-3] == rewards.shape[:-2]
    lead_dims = states.shape[:-2]
    states = states.reshape(-1, *states.shape[-3:])  # (*, building, Time, State)
    rewards = rewards.reshape(-1, *rewards.shape[-2:])  # (*, building, Time)
    price_total = np.maximum(-rewards[:, :, -1].sum(axis=-1) * price_factor, 0.)

    total_rewards = rewards[:, :, -tot_len:]
    # heat_reward = negative_total_consumption - solar_energy + non-shiftable energy
    heat_rew = rewards[:, :, -heat_len:] - states[:, :, -heat_len:, SGEN_IDX] + states[:, :, -heat_len:, NSL_IDX]
    heat_rew = np.minimum(0., heat_rew)
    current_heat_reward = heat_rew[:, :, -1]

    def get_slider_mean(slider):
        # if slider.shape[-1] == 1 and heat_len == 1:
        if heat_len == 1:
            return 0.
        else:
            return np.mean(slider, axis=-1)

    def ramping_reward(reward, ramp_window):
        if reward.shape[-1] > 1:  # TODO
            return np.abs(reward[:, :, -ramp_window:] - reward[:, :, -ramp_window - 1:-1]).sum(-1)
        else:
            return 0.

    mean_heat_reward = get_slider_mean(heat_rew)
    heat_advantage = beta * mean_normalize(current_heat_reward, mean_heat_reward)

    ramp_len = min(total_rewards.shape[-1] - 1, ramping_window)
    total_ramping = - alpha * ramping_reward(total_rewards, ramp_len)
    final_rewards = np.expand_dims(price_total, -1) * heat_advantage + total_ramping
    final_rewards = final_rewards.reshape(*lead_dims)
    reward_logger = [np.max(heat_advantage), np.min(heat_advantage),
                     np.max(total_ramping), np.min(total_ramping),
                     np.max(final_rewards), np.min(final_rewards),
                     np.max(price_total), np.min(price_total)]

    return final_rewards, reward_logger

# Reward function for the multi-agent (decentralized) agents
def reward_function_ma(electricity_demand):
    return electricity_demand

# Reward function for the single-agent (centralized) agent
def reward_function_sa(electricity_demand):
    total_energy_demand = 0
    for e in electricity_demand:
        total_energy_demand += -e

    price = max(total_energy_demand * 0.01, 0)

    for i in range(len(electricity_demand)):
        electricity_demand[i] = min(price * electricity_demand[i], 0)

    return sum(electricity_demand)

def reward_function_sa_all(electricity_demand):
    total_energy_demand = 0
    for e in electricity_demand:
        total_energy_demand += -e

    price = max(total_energy_demand * 0.01, 0)

    for i in range(len(electricity_demand)):
        electricity_demand[i] = min(price * electricity_demand[i], 0)

    demand_sum = sum(electricity_demand)

    for i in range(len(electricity_demand)):
        electricity_demand[i] = demand_sum

    return electricity_demand

def reward_function_sa_mean(electricity_demand):
    total_energy_demand = 0
    for e in electricity_demand:
        total_energy_demand += -e

    price = max(total_energy_demand * 0.01, 0)

    for i in range(len(electricity_demand)):
        electricity_demand[i] = min(price * electricity_demand[i], 0)

    demand_sum = np.mean(np.array(electricity_demand))

    for i in range(len(electricity_demand)):
        electricity_demand[i] = demand_sum

    return electricity_demand


text:
import albumentations as A
import numpy as np

class ObjectDetectionAugmentation:
    def __init__(self, config: dict):
        '''
        Args:
            config: dictionaly of augmentation params
            (example)
            {
                augment_args:
                    vertical_flip: {p: 0.5}
                    horizontal_flip: {p: 0.5}
                    random_crop: {width: 450, height: 450}
                    random_brightness_contrast: {p: 0.2}
                bbox_param: {format: 'coco', min_area: 100, min_visibility: 0.1, label_fields: []}
                compile: 'Compose'
            }
        '''
        # add other augmentations, if you want
        aug_list = []
        for key in config['augment_args']:
            if key == 'random_crop':
                aug = A.RandomCrop(**config['augment_args'][key])
            elif key == 'horizontal_flip':
                aug = A.HorizontalFlip(**config['augment_args'][key])
            elif key == 'vertical_flip':
                aug = A.VerticalFlip(**config['augment_args'][key])
            elif key =='random_brightness_contrast':
                aug = A.RandomBrightnessContrast(**config['augment_args'][key])
            aug_list.append(aug)

        if config['compile'] == 'compose':
            compile_augs = A.Compose
        elif config['compile'] == 'oneof':
            compile_augs = A.OneOf

        self.transformer = compile_augs(
            aug_list,
            bbox_params=A.BboxParams(**config['bbox_param'])
        )
        self.config = config

    def transform(self, img: np.array, bboxes: list):
        '''Transform set of image and bboxes
        Args:
            img(np.array): input image
            bboxes(list): list of bounding boxes
        Return:
            transformed_img(np.array): output image
            transformed_bboxes(list): list of bounding boxes
        '''
        transformed = self.transformer(image=img, bboxes=bboxes)
        transformed_img = transformed['image']
        transformed_bboxes = transformed['bboxes']
        return transformed_img, transformed_bboxes

    # Todo: Delete loop
    def transform_data(self, img_list: list, bboxes_list: list):
        '''Transform all data
       Args:
            img_list(list[np.array]): list of input images
            bboxes_list(list[list]): list of bounding boxes list
        Return:
            transformed_img_list(list[np.array]): list ofoutput image
            transformed_bboxes_list(list[list]): list of bounding boxes list
        '''
        transformed_img_list = []
        transformed_bboxes_list = []
        for _img, _bboxes in zip(img_list, bboxes_list):
            transformed_img, transformed_bboxes = self.transform(_img, _bboxes)
            transformed_img_list.append(transformed_img)
            transformed_bboxes_list.append(transformed_bboxes)
        return transformed_img_list, transformed_bboxes_list


text:
import numpy as np
import matplotlib.pyplot as plt
mu, sigma = 0, 0.1
V0 = 50
alpha = 0.9
T = 100
epsilon = np.random.normal(mu, [alpha**i*V0 for i in range(T)], T)
X0 = 10
X = [X0 + sum(epsilon[:i]) for i in range(T)]
plt.plot(X)
plt.show()


text:
import cv2
from keras.models import load_model
import numpy as np
print("IM IN CNN TESTING")
def predict(filename):
	model = load_model("./alex-cnn.h5")
	# filename = "/home/atharva/a2iot/deeplearning/DDAE/spectrograms/spectrograms_chunks_noise/chunk24.png"
	# filename = "/home/atharva/a2iot/deeplearning/DDAE/spectrograms/spectrograms_test_enhanced/chunk1107.png"

	img = cv2.imread(filename)
	img = cv2.resize(img, (224,224))
	img = np.reshape(img, (1,224,224,3))
	print(img.shape)
	pred = model.predict(img)
	print(pred)
	return pred

max_stars_repo_name: SharanRajani/SoundQX

max_issues_repo_name: SharanRajani/SoundQX

max_forks_repo_name: SharanRajani/SoundQX


text:
#!/usr/bin/python

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# customization of plotting
plt.rcParams['lines.linewidth']            = 4
plt.rcParams['font.size']                  = 18
plt.rcParams['font.weight']                = 'bold'
plt.rcParams['axes.labelsize']             = 'large'
plt.rcParams['figure.facecolor']           = 'white'
plt.rcParams['image.interpolation']        = 'none'
plt.rcParams['image.origin']               = 'lower'
plt.rcParams['contour.negative_linestyle'] = 'solid'
plt.rcParams['mathtext.default'] = 'regular' # match the font used for
                                             # regular text

# Defining exponential
def exponential(x, a, b):
    return a*np.exp(2*b*x)

files = ["k_0_2/5m_bz2.txt"]

# Loop over files
for i, fileName in enumerate(files):
    minN = 100
    maxR_squared = 0
    data = np.loadtxt(fileName)
    for n in np.linspace(minN, data.shape[0]-1, data.shape[0]-minN):
        x = data[0:n,0]
        y = data[0:n,1]
        try:
            # Main fitting
            popt, pcov = curve_fit(exponential, x, y, p0=(1,0.1))
            # Calculation of R^2
            residuals = y-exponential(x, *popt)
            ss_res = np.sum(residuals**2)
            ss_tot = np.sum((y-np.mean(y))**2)
            R_squared = 1-ss_res/ss_tot
            if R_squared > maxR_squared:
                maxR_squared = R_squared
                bestFit = popt
            percentage = float(n-minN)/(data.shape[0]-minN)*100
            print('{} {:5.2f} %: gamma = {:6.4f} ({:6.4f}); R2 = {:6.4f} ({:6.4f})'.format(fileName, percentage, popt[1], bestFit[1], R_squared, maxR_squared))
        except RuntimeError:
            print("Error - curve_fit failed")
    print('{}: gamma = {}'.format(fileName, bestFit[1]))
    plt.figure(i)
    plt.plot(data[:,0], data[:,1], '.')
    plt.plot(data[:,0], exponential(data[:,0], *bestFit))
    plt.title(fileName)
    plt.xlabel('Time [Ut]')
    plt.ylabel('B_z^2')
    plt.ylim([0,np.max(data[:,1])])
    #plt.draw()
plt.show()

max_stars_repo_name: ammarhakim/gkyl-paper-inp

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: ammarhakim/gkyl-paper-inp

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: ammarhakim/gkyl-paper-inp

max_forks_repo_licenses: MIT Unlicense

reason: import numpy,from scipy


text:
# Step 1: Load dataset, split into training and test sets, and scale features
import numpy as np
from sklearn.datasets import load_boston
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
# load boston housing price dataset
boston = load_boston()
x = boston.data
y = boston.target

# split into training and test sets, namely 80 percent of examples goes for the training, 20 percent goes for the test set
N_train = int(0.8 * x.shape[0])
x_train = x[:N_train, :]
y_train = y[:N_train]
x_test = x[N_train:, :]
y_test = y[N_train:]

# scale features by removing mean and dividing by the standard deviation
x_train_scaled = scaler.fit_transform(x_train)
x_test_scaled = scaler.fit_transform(x_test)

print("Xtrain scaled")
print(x_train_scaled.shape)
print("y train")
print(y_train.shape)
print("Xtest scaled")
print(x_test_scaled.shape)
print("y test")
print(y_test.shape)

# # Step 2: Add intercept terms and initialize parameters
interc_train = np.ones((x_train_scaled.shape[0], 1))
interc_test = np.ones((x_test_scaled.shape[0], 1))
# Note: If you run this step again, please run from step 1 because notebook keeps the value from the previous run
x_train_scaled = np.hstack((interc_train, x_train_scaled))

x_test_scaled = np.hstack((interc_test, x_test_scaled))

print("X_train scaled after bias term")
print(x_train_scaled.shape)
print("X_test scaled after bias term")
print(x_test_scaled.shape)
#
# # init parameters using random values
mu, sigma = 0, 0.5
number_of_features = x_train_scaled.shape[1]
theta = np.random.normal(mu, sigma, number_of_features)
print(theta)

#
#
# Step 3: Implement the gradient and the cost function
# In this step, you have to calculate the gradient. You can use the provided formula but the best way is to vectorize
# that formula for efficiency
def compute_gradient(x, y, theta):
    n = x.shape[0]
    return np.matmul(x.T, (np.matmul(x, theta) - y))/n

def compute_cost(x, y, theta):
    x_theta = np.matmul(x, theta)
    x_theta_minus_y = x_theta - y
    n = x.shape[0]
    return (1 / (2 * n)) * np.matmul(x_theta_minus_y.T, x_theta_minus_y)

# Step 4: Verify the gradient value
# In this step, you need to verify that the computed gradient is correct. The difference betweet the gradient and the
# approximate gradient should be very small (~10^-18)
def approximate_gradient(x, y, theta, epsilon):
    n_features = x.shape[1]
    app_grad = np.zeros(n_features)
    for i in range(n_features):
        temp = theta[i]
        theta[i] = temp + epsilon
        epsilon_plus = compute_cost(x, y, theta)
        theta[i] = temp - epsilon
        epsilon_minus = compute_cost(x, y, theta)
        app_grad[i] = (epsilon_plus - epsilon_minus) / (2 * epsilon)
    return app_grad

grad = compute_gradient(x_train_scaled, y_train, theta)
epsilon = 1e-4
app_grad = approximate_gradient(x_train_scaled, y_train, theta, epsilon)
print('Sum of gradient squared error: ', np.sum((grad - app_grad) ** 2))

# Step 5: Try gradient descent algorithm with different learning rates
import matplotlib.pyplot as plt
import copy

# try different values for the learning rate
learning_rates = [0.001,0.003,0.01,0.03,0.1,0.3]

# this matrix keeps the learned parameters
theta_matrix = np.zeros((len(learning_rates),x_train_scaled.shape[1]))

# number of training iterations
N_iterations = 100

# prepare to plot
plt.subplot(111)

# calculate cost value and update theta
for indx,alpha in enumerate(learning_rates):
    # keep the cost value for each training step
    J = np.zeros(N_iterations)

    # initialize new parameters using random distribution
    theta = 0.5 * np.random.randn(x_train_scaled.shape[1])
    for step in range(N_iterations):
        # update theta
        grad = compute_gradient(x_train_scaled, y_train, theta)
        theta = theta - alpha*grad

        # save the value of theta
        theta_matrix[indx,:] = theta

        # calculate the cost on traing set
        J[step] = compute_cost(x_train_scaled, y_train, theta)
    # plot cost function
    plt.plot(J)
plt.xlabel('Training step')
plt.ylabel('Cost')
plt.legend(('0.001','0.003','0.01','0.03','0.1','0.3'), loc='upper right')
plt.show()
#
#
# # Step 6: Predict the price of house
# # You have to select the best theta you found
theta = theta_matrix[4,:]
predict_price = np.matmul(x_test_scaled, theta)
#
# calculate the cost for the test set
test_cost = compute_cost(x_test_scaled, y_test, theta)
print('test cost: ', test_cost)

# plot the ground truth and the predicted
x_axis = np.linspace(1,len(y_test),len(y_test))
plt.plot(x_axis,y_test,'b',x_axis,predict_price,'r')
plt.legend(('Ground truth','Predicted'))
plt.show()


text:
from functools import partial
from matplotlib.pyplot import xcorr

import numpy as np
import jax
import jax.numpy as jnp
import flax
from flax import linen as nn
import distrax

from .jax_utils import batch_to_jax, extend_and_repeat, next_rng

def update_target_network(main_params, target_params, tau):
    return jax.tree_multimap(
        lambda x, y: tau * x + (1.0 - tau) * y,
        main_params, target_params
    )

# def multiple_action_q_function(forward):
#     # Forward the q function with multiple actions on each state, to be used as a decorator
#     def wrapped(self, observations, actions, **kwargs):
#         multiple_actions = False
#         batch_size = observations.shape[0]
#         if actions.ndim == 3 and observations.ndim == 2:
#             multiple_actions = True
#             observations = extend_and_repeat(observations, 1, actions.shape[1]).reshape(-1, observations.shape[-1])
#             actions = actions.reshape(-1, actions.shape[-1])

#         q_values = forward(self, observations, actions, **kwargs)
#         if multiple_actions:
#             q_values = q_values.reshape(batch_size, -1)
#         return q_values
#     return wrapped

def multiple_action_q_function(forward):
    # Forward the q function with multiple actions on each state, to be used as a decorator
    def wrapped(self, observations, actions, **kwargs):
        multiple_actions = False
        batch_size = observations.shape[0]
        if actions.ndim == 3 and observations.ndim == 2:
            multiple_actions = True
            observations = extend_and_repeat(observations, 1, actions.shape[1]).reshape(-1, observations.shape[-1])
            actions = actions.reshape(-1, actions.shape[-1])

        q_last_layer, q_values = forward(self, observations, actions, **kwargs)
        if multiple_actions:
            q_last_layer = q_last_layer.reshape(batch_size, -1)
            q_values = q_values.reshape(batch_size, -1)
        return q_last_layer, q_values
    return wrapped

def multiple_action_encode_function(forward):
    # Forward the rep function with multiple actions on each state, to be used as a decorator
    def wrapped(self, rng, observations, actions, **kwargs):
        multiple_actions = False
        batch_size = observations.shape[0]
        if actions.ndim == 3 and observations.ndim == 2:
            repeat = actions.shape[1]
            multiple_actions = True
            observations = extend_and_repeat(observations, 1, actions.shape[1]).reshape(-1, observations.shape[-1])
            actions = actions.reshape(-1, actions.shape[-1])
        samples, log_probs = forward(self, rng, observations, actions, **kwargs)
        if multiple_actions:
            samples = samples.reshape(batch_size, repeat, -1)
            log_probs = log_probs.reshape(batch_size, repeat, -1)
        return samples, log_probs
    return wrapped

def multiple_action_decode_function(forward):
    # Forward the rep function with multiple actions on each state, to be used as a decorator
    def wrapped(self, observations, actions, **kwargs):
        multiple_actions = False
        batch_size = observations.shape[0]
        if actions.ndim == 3 and observations.ndim == 2:
            repeat = actions.shape[1]
            multiple_actions = True
            observations = extend_and_repeat(observations, 1, actions.shape[1]).reshape(-1, observations.shape[-1])
            actions = actions.reshape(-1, actions.shape[-1])
        q_values = forward(self, observations, actions, **kwargs)
        if multiple_actions:
            q_values = q_values.reshape(batch_size, repeat, -1)
        return q_values
    return wrapped

class Scalar(nn.Module):
    init_value: float

    def setup(self):
        self.value = self.param('value', lambda x:self.init_value)

    def __call__(self):
        return self.value

class FullyConnectedNetwork(nn.Module):
    output_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False
    # batch_norm: bool = False

    @nn.compact
    def __call__(self, input_tensor):
        x = input_tensor
        hidden_sizes = [int(h) for h in self.arch.split('-')]
        for h in hidden_sizes:
            if self.orthogonal_init:
                x = nn.Dense(
                    h,
                    kernel_init=jax.nn.initializers.orthogonal(jnp.sqrt(2.0)),
                    bias_init=jax.nn.initializers.zeros
                )(x)
            else:
                x = nn.Dense(h)(x)
            # if self.batch_norm:
            #     x = nn.BatchNorm(use_running_average=not train_mode, momentum=0.9,
            #             epsilon=1e-5,
            #             dtype=jnp.float32)(x)
            x = nn.relu(x)

        if self.orthogonal_init:
            output = nn.Dense(
                self.output_dim,
                kernel_init=jax.nn.initializers.orthogonal(1e-2),
                bias_init=jax.nn.initializers.zeros
            )(x)
        else:
            output = nn.Dense(
                self.output_dim,
                kernel_init=jax.nn.initializers.variance_scaling(
                    1e-2, 'fan_in', 'uniform'
                ),
                bias_init=jax.nn.initializers.zeros
            )(x)
        return output

class FullyConnectedNetworkWithLastLayer(nn.Module):
    output_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False

    @nn.compact
    def __call__(self, input_tensor):
        x = input_tensor
        batch, _ = jnp.shape(x)
        hidden_sizes = [int(h) for h in self.arch.split('-')]
        for h in hidden_sizes:
            if self.orthogonal_init:
                x = nn.Dense(
                    h,
                    kernel_init=jax.nn.initializers.orthogonal(jnp.sqrt(2.0)),
                    bias_init=jax.nn.initializers.zeros
                )(x)
            else:
                x = nn.Dense(h)(x)
            # x = nn.LayerNorm()(x)
            x = nn.relu(x)
        normalized = jnp.reshape(jnp.sqrt(jnp.sum(x**2, axis=-1) + 1e-6), (batch,1))
        x = x / normalized
        # x = x / (jnp.sqrt(jnp.sum(x**2, axis=-1) + 1e-6))

        if self.orthogonal_init:
            output = nn.Dense(
                self.output_dim,
                kernel_init=jax.nn.initializers.orthogonal(1e-2),
                bias_init=jax.nn.initializers.zeros
            )(x)
        else:
            output = nn.Dense(
                self.output_dim,
                kernel_init=jax.nn.initializers.variance_scaling(
                    1e-2, 'fan_in', 'uniform'
                ),
                bias_init=jax.nn.initializers.zeros
            )(x)
        return x, output

class FullyConnectedQFunction(nn.Module):
    observation_dim: int
    action_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False

    @nn.compact
    @multiple_action_q_function
    def __call__(self, observations, actions):
        x = jnp.concatenate([observations, actions], axis=-1)
        # x = FullyConnectedNetwork(output_dim=1, arch=self.arch, orthogonal_init=self.orthogonal_init)(x)
        last_layer_x, x = FullyConnectedNetworkWithLastLayer(output_dim=1, arch=self.arch, orthogonal_init=self.orthogonal_init)(x)
        return last_layer_x, jnp.squeeze(x, -1)

class FullyConnectedActionQFunction(nn.Module):
    observation_dim: int
    action_dim: int
    output_dim: int = 1
    arch: str = '256-256'
    orthogonal_init: bool = False
    normalize: bool = False

    @nn.compact
    @multiple_action_q_function
    def __call__(self, observations, actions):
        x = jnp.concatenate([observations, actions], axis=-1)
        batch, _ = jnp.shape(x)
        hidden_sizes = [int(h) for h in self.arch.split('-')]
        for h in hidden_sizes:
            if self.orthogonal_init:
                x = nn.Dense(
                    h,
                    kernel_init=jax.nn.initializers.orthogonal(jnp.sqrt(2.0)),
                    bias_init=jax.nn.initializers.zeros
                )(x)
            else:
                x = nn.Dense(h)(jnp.concatenate([x, actions], axis=-1))
            x = nn.relu(x)

        if self.normalize:
            normalized = jnp.reshape(jnp.sqrt(jnp.sum(x**2, axis=-1) + 1e-6), (batch,1))
            x = x / normalized

        if self.orthogonal_init:
            output = nn.Dense(
                self.output_dim,
                kernel_init=jax.nn.initializers.orthogonal(1e-2),
                bias_init=jax.nn.initializers.zeros
            )(jnp.concatenate([x, actions], axis=-1))
        else:
            output = nn.Dense(
                self.output_dim,
                kernel_init=jax.nn.initializers.variance_scaling(
                    1e-2, 'fan_in', 'uniform'
                ),
                bias_init=jax.nn.initializers.zeros
            )(jnp.concatenate([x, actions], axis=-1))
        return x, jnp.squeeze(output, -1)

class TanhGaussianPolicy(nn.Module):
    observation_dim: int
    action_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False
    log_std_multiplier: float = 1.0
    log_std_offset: float = -1.0
    action_scale: float = 1.0

    def setup(self):
        self.base_network = FullyConnectedNetwork(
            output_dim=2 * self.action_dim, arch=self.arch, orthogonal_init=self.orthogonal_init
        )
        self.log_std_multiplier_module = Scalar(self.log_std_multiplier)
        self.log_std_offset_module = Scalar(self.log_std_offset)

    def log_prob(self, observations, actions):
        if actions.ndim == 3:
            observations = extend_and_repeat(observations, 1, actions.shape[1])
        base_network_output = self.base_network(observations)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        action_distribution = distrax.Transformed(
            distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
            distrax.Block(distrax.Tanh(), ndims=1)
        )
        return action_distribution.log_prob(actions / self.action_scale)

    def __call__(self, rng, observations, deterministic=False, repeat=None):
        if repeat is not None:
            observations = extend_and_repeat(observations, 1, repeat)
        base_network_output = self.base_network(observations)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        action_distribution = distrax.Transformed(
            distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
            distrax.Block(distrax.Tanh(), ndims=1)
        )
        if deterministic:
            mean = jnp.clip(mean, -6, 6)
            samples = jnp.tanh(mean)
            log_prob = action_distribution.log_prob(samples)
        else:
            samples, log_prob = action_distribution.sample_and_log_prob(seed=rng)
        samples = samples * self.action_scale
        return samples, log_prob

class ActionRepresentationPolicy(nn.Module):
    observation_dim: int
    action_dim: int
    latent_action_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False
    no_tanh: bool = False
    log_std_multiplier: float = 1.0
    log_std_offset: float = -1.0
    # batch_norm: bool = True

    def setup(self):
        self.base_network = FullyConnectedNetwork(
            output_dim=2 * self.latent_action_dim, arch=self.arch, orthogonal_init=self.orthogonal_init
        )
        self.log_std_multiplier_module = Scalar(self.log_std_multiplier)
        self.log_std_offset_module = Scalar(self.log_std_offset)

    def log_prob(self, observations, actions, latent_actions):
        if actions.ndim == 3:
            observations = extend_and_repeat(observations, 1, actions.shape[1])
        x = jnp.concatenate([observations, actions], axis=-1)
        base_network_output = self.base_network(x)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        if self.no_tanh:
            action_distribution = distrax.MultivariateNormalDiag(mean, jnp.exp(log_std))
        else:
            action_distribution = distrax.Transformed(
                distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
                distrax.Block(distrax.Tanh(), ndims=1)
            )
        return action_distribution.log_prob(latent_actions)

    @nn.compact
    @multiple_action_encode_function
    def __call__(self, rng, observations, actions, deterministic=False, repeat=None):
        if repeat is not None:
            observations = extend_and_repeat(observations, 1, repeat)
        x = jnp.concatenate([observations, actions], axis=-1)
        base_network_output = self.base_network(x)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        if self.no_tanh:
            action_distribution = distrax.MultivariateNormalDiag(mean, jnp.exp(log_std))
        else:
            action_distribution = distrax.Transformed(
                distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
                distrax.Block(distrax.Tanh(), ndims=1)
            )
        if deterministic:
            samples = jnp.tanh(mean)
            log_prob = action_distribution.log_prob(samples)
        else:
            samples, log_prob = action_distribution.sample_and_log_prob(seed=rng)

        return samples, log_prob

    def get_statistics(self, rng, observations, actions, deterministic=False, repeat=None):
        if repeat is not None:
            observations = extend_and_repeat(observations, 1, repeat)
        x = jnp.concatenate([observations, actions], axis=-1)
        base_network_output = self.base_network(x)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        if self.no_tanh:
            action_distribution = distrax.MultivariateNormalDiag(mean, jnp.exp(log_std))
        else:
            action_distribution = distrax.Transformed(
                distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
                distrax.Block(distrax.Tanh(), ndims=1)
            )
        if deterministic:
            samples = jnp.tanh(mean)
            log_prob = action_distribution.log_prob(samples)
        else:
            samples, log_prob = action_distribution.sample_and_log_prob(seed=rng)

        return samples, mean, log_std

class ActionOnlyRepresentationPolicy(nn.Module):
    action_dim: int
    latent_action_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False
    no_tanh: bool = False
    log_std_multiplier: float = 1.0
    log_std_offset: float = -1.0

    def setup(self):
        self.base_network = FullyConnectedNetwork(
            output_dim=2 * self.latent_action_dim, arch=self.arch, orthogonal_init=self.orthogonal_init
        )
        self.log_std_multiplier_module = Scalar(self.log_std_multiplier)
        self.log_std_offset_module = Scalar(self.log_std_offset)

    def log_prob(self, actions, latent_actions):
        if actions.ndim == 3:
            observations = extend_and_repeat(observations, 1, actions.shape[1])
        base_network_output = self.base_network(actions)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        if self.no_tanh:
            action_distribution = distrax.MultivariateNormalDiag(mean, jnp.exp(log_std))
        else:
            action_distribution = distrax.Transformed(
                distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
                distrax.Block(distrax.Tanh(), ndims=1)
            )
        return action_distribution.log_prob(latent_actions)

    @nn.compact
    def __call__(self, rng, actions, deterministic=False, repeat=None):
        if repeat is not None:
            observations = extend_and_repeat(observations, 1, repeat)
        base_network_output = self.base_network(actions)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        if self.no_tanh:
            action_distribution = distrax.MultivariateNormalDiag(mean, jnp.exp(log_std))
        else:
            action_distribution = distrax.Transformed(
                distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
                distrax.Block(distrax.Tanh(), ndims=1)
            )
        if deterministic:
            samples = jnp.tanh(mean)
            log_prob = action_distribution.log_prob(samples)
        else:
            samples, log_prob = action_distribution.sample_and_log_prob(seed=rng)

        return samples, log_prob

    def get_statistics(self, rng, actions, deterministic=False, repeat=None):
        if repeat is not None:
            observations = extend_and_repeat(observations, 1, repeat)
        base_network_output = self.base_network(actions)
        mean, log_std = jnp.split(base_network_output, 2, axis=-1)
        log_std = self.log_std_multiplier_module() * log_std + self.log_std_offset_module()
        log_std = jnp.clip(log_std, -20.0, 2.0)
        if self.no_tanh:
            action_distribution = distrax.MultivariateNormalDiag(mean, jnp.exp(log_std))
        else:
            action_distribution = distrax.Transformed(
                distrax.MultivariateNormalDiag(mean, jnp.exp(log_std)),
                distrax.Block(distrax.Tanh(), ndims=1)
            )
        if deterministic:
            samples = jnp.tanh(mean)
            log_prob = action_distribution.log_prob(samples)
        else:
            samples, log_prob = action_distribution.sample_and_log_prob(seed=rng)

        return samples, mean, log_std

class ActionDecoder(nn.Module):
    observation_dim: int
    latent_action_dim: int
    action_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False

    @nn.compact
    @multiple_action_decode_function
    def __call__(self, observations, latent_actions):
        x = jnp.concatenate([observations, latent_actions], axis=-1)
        hidden_sizes = [int(h) for h in self.arch.split('-')]
        for h in hidden_sizes:
            if self.orthogonal_init:
                x = nn.Dense(
                    h,
                    kernel_init=jax.nn.initializers.orthogonal(jnp.sqrt(2.0)),
                    bias_init=jax.nn.initializers.zeros
                )(x)
            else:
                x = nn.Dense(h)(x)
            x = nn.relu(x)

        if self.orthogonal_init:
            x = nn.Dense(
                self.action_dim,
                kernel_init=jax.nn.initializers.orthogonal(1e-2),
                bias_init=jax.nn.initializers.zeros
            )(x)
        else:
            x = nn.Dense(
                self.action_dim,
                kernel_init=jax.nn.initializers.variance_scaling(
                    1e-2, 'fan_in', 'uniform'
                ),
                bias_init=jax.nn.initializers.zeros
            )(x)

        output = nn.tanh(x)
        return output

class ActionSeperatedDecoder(nn.Module):
    observation_dim: int
    latent_action_dim: int
    action_dim: int
    arch: str = '256-256'
    orthogonal_init: bool = False
    # batch_norm: bool = True

    @nn.compact
    @multiple_action_decode_function
    def __call__(self, observations, latent_actions):
        x = observations
        hidden_sizes = [int(h) for h in self.arch.split('-')]
        for h in hidden_sizes:
            if self.orthogonal_init:
                x = nn.Dense(
                    h,
                    kernel_init=jax.nn.initializers.orthogonal(jnp.sqrt(2.0)),
                    bias_init=jax.nn.initializers.zeros
                )(jnp.concatenate([x, latent_actions], axis=-1))
            else:
                x = nn.Dense(h)(jnp.concatenate([x, latent_actions], axis=-1))
            # if self.batch_norm:
            #     x = nn.BatchNorm(use_running_average=not train_mode, momentum=0.9,
            #             epsilon=1e-5,
            #             dtype=jnp.float32)(x)
            x = nn.relu(x)

        if self.orthogonal_init:
            x = nn.Dense(
                self.action_dim,
                kernel_init=jax.nn.initializers.orthogonal(1e-2),
                bias_init=jax.nn.initializers.zeros
            )(jnp.concatenate([x, latent_actions], axis=-1))
        else:
            x = nn.Dense(
                self.action_dim,
                kernel_init=jax.nn.initializers.variance_scaling(
                    1e-2, 'fan_in', 'uniform'
                ),
                bias_init=jax.nn.initializers.zeros
            )(jnp.concatenate([x, latent_actions], axis=-1))

        output = nn.tanh(x)
        return output

class Discriminator(nn.Module):
    observation_dim: int
    latent_action_dim: int
    arch: str = '512-256'
    dropout: bool = True

    @nn.compact
    def __call__(self, observations, latent_actions, train=False):
        x = jnp.concatenate([observations, latent_actions], axis=-1)
        hidden_sizes = [int(h) for h in self.arch.split('-')]
        for h in hidden_sizes:
            x = nn.Dense(h)(x)
            # dropout
            # layer norm
            x = nn.leaky_relu(x, 0.2)
            if self.dropout:
                x = nn.Dropout(0.1)(x, deterministic=not train)
        output = nn.Dense(1)(x)
        return output

class SamplerPolicy(object):

    def __init__(self, policy, params):
        self.policy = policy
        self.params = params

    def update_params(self, params):
        self.params = params
        return self

    @partial(jax.jit, static_argnames=('self', 'deterministic'))
    def act(self, params, rng, observations, deterministic):
        return self.policy.apply(params, rng, observations, deterministic, repeat=None)

    def __call__(self, observations, deterministic=False):
        actions, _ = self.act(self.params, next_rng(), observations, deterministic=deterministic)
        assert jnp.all(jnp.isfinite(actions))
        return jax.device_get(actions)

class SamplerDecoder(object):

    def __init__(self, decoder, params):
        self.decoder = decoder
        self.params = params

    def update_params(self, params):
        self.params = params
        return self

    @partial(jax.jit, static_argnames=('self'))
    def act(self, params, observations, actions_rep):
        return self.decoder.apply(params, observations, actions_rep)

    def __call__(self, observations, actions_rep):
        actions = self.act(self.params, observations, actions_rep)
        assert jnp.all(jnp.isfinite(actions))
        return jax.device_get(actions)

class SamplerEncoder(object):

    def __init__(self, encoder, params):
        self.encoder = encoder
        self.params = params

    def update_params(self, params):
        self.params = params
        return self

    @partial(jax.jit, static_argnames=('self'))
    def act(self, params, rng, observations, actions):
        return self.encoder.apply(params, rng, observations, actions)[0]

    def __call__(self, rng, observations, actions):
        actions = self.act(self.params, rng, observations, actions)
        assert jnp.all(jnp.isfinite(actions))
        return jax.device_get(actions)

reason: import numpy,import jax


text:
#!/usr/bin/python
import sys
sys.stdout = sys.stderr
from cgi import parse_qsl
import os
os.environ[ 'HOME' ] = '/tmp/'
import matplotlib
matplotlib.interactive(0)
matplotlib.rcParams['path.simplify'] = True

from matplotlib.figure import Figure
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from cStringIO import StringIO
import numpy as np

path = os.path.dirname(__file__)
histones = {}
for ichr in range(1, 6):
    schr = str(ichr)
    f = os.path.join(path, "histone.%s.bin") % schr
    assert os.access(f, os.R_OK)

    histones[schr] = np.memmap(open(os.path.join(path, 'histone.%s.bin') % schr), dtype=np.float32, mode='r')

colors = [(1, 0, 0), (1, 1, 0), (0, 0, 1)]

def application(env, start_response):
    start_response("200 OK", [("Content-type", "image/png")])

    p = dict(parse_qsl(env['QUERY_STRING']))

    xmin = int(p['xmin'])
    xmax = max(int(p['xmax']), 0)
    xmin1 = max(1, xmin)

    seqid = p['seqid']

    hist = np.array(histones[seqid][xmin - 1: xmax])

    height = int(p.get('height', 300))

    f = Figure(frameon=False)
    f.canvas = FigureCanvas(f)
    dpi = 128.
    f.set_size_inches(int(p['width'])/dpi, height/dpi)

    ax = f.add_axes((0, 0, 1, 1), frameon=False, xticks=(), yticks=())
    if xmax != 0:
        ax.set_autoscale_on(0)
        ax.set_xlim(xmin - 1,xmax)

        ax.set_ylim(0.3, 0.9)
        c = 'k'
        ax.fill_between(
            np.arange(xmin1 - 1, xmax),
            np.zeros_like(hist),
            hist, facecolor=c, edgecolor='none', alpha=0.4) #, lw=0.8)

    io = StringIO()
    f.savefig(io, format='png', dpi=dpi)
    io.seek(0)
    return [io.read()]

ext: wsgi


text:
import os
import torch
import numpy as np
from torch import optim, nn
from torch.utils.data import TensorDataset, ConcatDataset, DataLoader

def agent_wapper(config, agent):
    def fn(obs):
        with torch.no_grad():
            obs = obs.astype(np.float32)
            assert len(obs.shape) == 2
            obs = torch.from_numpy(obs).to(config.device)
            action = agent(obs)
        return action.cpu().numpy()
    return fn

def fit_dataset(config, agent, dataset, n_epochs):
    optimizer = optim.Adam(agent.parameters(), lr=config.lr, weight_decay=config.L2)
    loss_fn = nn.MSELoss()
    dataloader = DataLoader(dataset, batch_size=config.batch_size, shuffle=True)

    step = 0
    best_reward = None
    loss_his = []

    for k in range(n_epochs):
        for batch in dataloader:
            obs, gold_actions = batch
            pred_actions = agent(obs)
            loss = loss_fn(pred_actions, gold_actions)

            optimizer.zero_grad()
            loss.backward()
            optimizer.step()

            loss_his.append(loss.item())

            if step % config.eval_steps == 0:
                avrg_mean, avrg_std = Eval(config, agent_wapper(config, agent))
                avrg_loss = np.mean(loss_his)
                loss_his = []
                print('[epoch {}  step {}] loss: {:.4f}  r_mean: {:.2f}  r_std: {:.2f}'.format(
                    k + 1, step, avrg_loss, avrg_mean, avrg_std))

                avrg_reward = avrg_mean - avrg_std
                if best_reward is None or best_reward < avrg_reward:
                    best_reward = avrg_reward
                    save_model(config, agent, config.model_save_path)

            step += 1

    load_model(config, agent, config.model_save_path)

def BehavioralCloning(config, agent, expert):

    # get expert demonstration
    expert_obs, expert_actions, *_ = run_agent(config, expert, config.n_expert_rollouts)
    expert_obs = torch.from_numpy(expert_obs).to(config.device)
    expert_actions = torch.from_numpy(expert_actions).to(config.device)
    dataset = TensorDataset(expert_obs, expert_actions)

    # training agent
    fit_dataset(config, agent, dataset, config.epochs)

    return agent_wapper(config, agent)

def DAgger(config, agent, expert):
    # get expert demonstration
    expert_obs, expert_actions, *_ = run_agent(config, expert, config.n_expert_rollouts)
    expert_obs = torch.from_numpy(expert_obs).to(config.device)
    expert_actions = torch.from_numpy(expert_actions).to(config.device)
    dataset = TensorDataset(expert_obs, expert_actions)

    for k in range(config.n_dagger_iter):
        # training agent
        fit_dataset(config, agent, dataset, config.epochs)

        # run agent to get new on-policy observations
        new_obs, *_ = run_agent(config, agent_wapper(config, agent), config.n_dagger_rollouts)
        expert_actions = expert(new_obs)

        new_obs = torch.from_numpy(new_obs).to(config.device)
        expert_actions = torch.from_numpy(expert_actions).to(config.device)
        new_data = TensorDataset(new_obs, expert_actions)

        # add new data to dataset
        dataset = ConcatDataset([dataset, new_data])

        avrg_mean, avrg_std = Eval(config, agent_wapper(config, agent))
        print('[DAgger iter {}] r_mean: {:.2f}  r_std: {:.2f}'.format(k + 1, avrg_mean, avrg_std))

    return agent_wapper(config, agent)

def run_agent(config, agent, num_rollouts):
    env = config.env
    max_steps = env.spec.timestep_limit

    returns = []
    observations = []
    actions = []
    for _ in range(num_rollouts):
        obs = env.reset()
        done = False
        totalr = 0
        steps = 0
        while not done:
            action = agent(obs[None, :])
            action = action.reshape(-1)
            observations.append(obs)
            actions.append(action)
            obs, r, done, _ = env.step(action)
            totalr += r
            steps += 1
            if steps >= max_steps:
                break
        returns.append(totalr)

    avrg_mean, avrg_std = np.mean(returns), np.std(returns)
    observations = np.array(observations).astype(np.float32)
    actions = np.array(actions).astype(np.float32)

    return observations, actions, avrg_mean, avrg_std

def Eval(config, agent):
    *_, avrg_mean, avrg_std = run_agent(config, agent, config.n_eval_rollouts)

    return avrg_mean, avrg_std

def save_model(config, model, PATH):
    if not os.path.exists(PATH):
        os.makedirs(PATH)
    PATH = PATH + config.envname + '-' + 'parameters.tar'
    torch.save(model.state_dict(), PATH)
    print('model saved.')

def load_model(config, model, PATH):
    PATH = PATH + config.envname + '-' + 'parameters.tar'
    model.load_state_dict(torch.load(PATH))
    print('model loaded.')


text:
"""Unit test for the diffusion_maps module.

"""

import logging
import unittest

import numpy as np
# import matplotlib.pyplot as plt

from diffusion_maps import DiffusionMaps, downsample

from .misc import make_strip

class DiffusionMapsTest(unittest.TestCase):
    def setUp(self):
        logging.basicConfig(level=logging.DEBUG)
        self.xmin = 0.0
        self.ymin = 0.0
        self.width = 1.0
        self.height = 1e-1
        self.num_samples = 50000
        self.data = make_strip(self.xmin, self.ymin,
                               self.width, self.height,
                               self.num_samples)

    @staticmethod
    def _compute_rayleigh_quotients(matrix, eigenvectors):
        """Compute Rayleigh quotients."""
        N = eigenvectors.shape[0]
        rayleigh_quotients = np.zeros(N)
        for n in range(N):
            v = eigenvectors[n, :]
            rayleigh_quotients[n] = np.dot(v, matrix @ v) / np.dot(v, v)
        rayleigh_quotients = np.sort(np.abs(rayleigh_quotients))
        return rayleigh_quotients[::-1]

    def test_accuracy(self):
        num_samples = 10000
        logging.debug('Computing diffusion maps on a matrix of size {}'
                      .format(num_samples))
        num_eigenpairs = 10
        epsilon = 5e-1
        downsampled_data = downsample(self.data, num_samples)

        dm = DiffusionMaps(downsampled_data, epsilon,
                           num_eigenpairs=num_eigenpairs)

        ew = dm.eigenvalues
        rq = self._compute_rayleigh_quotients(dm.kernel_matrix,
                                              dm.eigenvectors)

        logging.debug('Eigenvalues: {}'.format(ew))
        logging.debug('Rayleigh quotients: {}'.format(rq))

        self.assertTrue(np.allclose(np.abs(ew), np.abs(rq)))

        # import scipy.io
        # scipy.io.mmwrite('kernel_matrix.mtx', dm.kernel_matrix)

    def test_multiple_epsilon_values(self):
        num_samples = 5000
        num_maps = 10
        num_eigenpairs = 10
        epsilon_min, epsilon_max = 1e-1, 1e1
        epsilons = np.logspace(np.log10(epsilon_min),
                               np.log10(epsilon_max), num_maps)

        downsampled_data = downsample(self.data, num_samples)

        evs = np.zeros((num_maps, num_eigenpairs, downsampled_data.shape[0]))
        ews = np.zeros((num_maps, num_eigenpairs))

        logging.basicConfig(level=logging.WARNING)

        for i, epsilon in enumerate(reversed(epsilons)):
            dm = DiffusionMaps(downsampled_data, epsilon,
                               num_eigenpairs=num_eigenpairs)

            evs[i, :, :] = dm.eigenvectors
            ews[i, :] = dm.eigenvalues

            ew = dm.eigenvalues
            rq = self._compute_rayleigh_quotients(dm.kernel_matrix,
                                                  dm.eigenvectors)
            self.assertTrue(np.allclose(np.abs(ew), np.abs(rq)))

            # plt.title('$\\epsilon$ = {:.3f}'.format(epsilon))
            # for k in range(1, 10):
            #     plt.subplot(2, 5, k)
            #     plt.scatter(downsampled_data[:, 0], downsampled_data[:, 1],
            #                 c=evs[i, k, :])
            #     plt.xlim([self.xmin, self.xmin + self.width])
            #     plt.ylim([self.ymin, self.ymin + self.height])
            #     plt.tight_layout()
            #     plt.gca().set_title('$\\psi_{}$'.format(k))
            # plt.subplot(2, 5, 10)
            # plt.step(range(ews[i, :].shape[0]), np.abs(ews[i, :]))
            # plt.title('epsilon = {:.2f}'.format(epsilon))
            # plt.show()

if __name__ == '__main__':
    import os
    verbose = os.getenv('VERBOSE')
    if verbose is not None:
        logging.basicConfig(level=logging.DEBUG, format='%(message)s')
    else:
        logging.basicConfig(level=logging.ERROR, format='%(message)s')
    unittest.main()

max_stars_repo_name: vishalbelsare/diffusion-maps

max_issues_repo_name: vishalbelsare/diffusion-maps

max_forks_repo_name: vishalbelsare/diffusion-maps

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 20 10:33:37 2014

@author: schackv
"""

import matplotlib.pyplot as plt
import shapewarp.ASM as ASM
import shapewarp.warp as warp
import shapewarp.plotting as shplt
import example_data
import scipy.misc
from numpy import isnan

def demo():
    """Demonstrate the shape model and warping functionality"""

    # Read shape data
    data = example_data.examples()
    landmarks = data.landmarks()
    print(landmarks.shape)

    # Plot shapes
    [shplt.plot_shape(L,'k') for L in landmarks]
    plt.show()

    # Active shape model (with generalized Procrustes inside)
    asm = ASM.ASM()
    asm.build(landmarks)

    ######## Plotting of shape model ######
    f, ax = plt.subplots(2,2)
    # Plot mean shape on top of aligned input shapes
    plt.sca(ax[0,0])
    [shplt.plot_shape(L,'k') for L in asm.AlignedShapes.T]
    shplt.plot_shape(asm.MeanShape,'r',linewidth=2.0)
    plt.title('Input shapes and mean shape')

    # Show covariance matrix as image
    plt.sca(ax[0,1])
    plt.imshow(asm.Covariance)
    plt.title('Covariance matrix')

    # Plot the first two eigenmodes
    plt.sca(ax[1,0])
    shplt.plot_mode(asm.MeanShape,asm.PCModes[0],3)
    plt.title('Mean shape +/- 3 std of first mode')
    plt.sca(ax[1,1])
    shplt.plot_mode(asm.MeanShape,asm.PCModes[1],3)
    plt.title('Mean shape +/- 3 std of second mode')

    for a in ax.flatten():
        a.axis('image')
        a.invert_yaxis()

    plt.show()
    ######################################

    ######## Warping demo ################
    image = data.images[ [i.has_image() for i in data.images].index(True) ]
    im = image.load()
    landmarks = image.landmarks()
    W = warp.Warper(asm.MeanShape,scale=0.5)        # Warp to 0.3 of the full scale of the mean shape
    warped_image = W.warp_image(im,landmarks)

    f, ax = plt.subplots(1,3)
    titles = ['Original image','Warped image','Mask']
    for idx,im in enumerate( (im,warped_image,W.mask)):
        plt.sca(ax[idx])
        plt.imshow(im)
        if idx==0:
            shplt.plot_shape(landmarks,'r')
        plt.axis('image')
        plt.title(titles[idx])
    plt.show()

    # Save image
    mask = W.mask
    warped_image[isnan(warped_image)]=0
    mask[isnan(mask)] = 0
    scipy.misc.toimage(warped_image,cmin=0,cmax=1).save('warpedImage.png')
    scipy.misc.toimage(mask,cmin=0,cmax=1).save('mask.png')
    ######################################

# when executed, just run asm_demo():
if __name__ == '__main__':
    demo()

max_stars_repo_name: schackv/shapewarp

max_issues_repo_name: schackv/shapewarp

max_forks_repo_name: schackv/shapewarp

reason: from numpy,import scipy


text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pickle

data_path = 'data/hapt/RawData/'
label = pd.read_csv(data_path+'labels.txt', delimiter=' ', header=None,
                    names=['exp_id', 'user_id', 'act_id', 'start', 'end'])
print(label.head())
windows = []
labels = []
prev_exp = 0
df_acc_exp = None
df_gyro_exp = None
for exp_id, user_id, act_id, start, end in label.values:
    # print(exp_id, user_id)
    if exp_id != prev_exp:
        file_acc = data_path + 'acc_exp{:02}_user{:02}.txt'.format(exp_id, user_id)
        file_gyro = data_path + 'gyro_exp{:02}_user{:02}.txt'.format(exp_id, user_id)

        # print(file_acc)
        # print(file_gyro)

        df_acc_exp = pd.read_csv(file_acc, delimiter=' ', header=None, names=['ax', 'ay', 'az'])
        df_gyro_exp = pd.read_csv(file_gyro, delimiter=' ', header=None, names=['gx', 'gy', 'gz'])
        prev_exp = exp_id
        # print(df_acc_exp.shape, df_gyro_exp.shape)

    acc_win = df_acc_exp.iloc[start:end].values
    gyro_win = df_gyro_exp.iloc[start: end].values
    # print(end - start)

    # print(acc_win)

    win = np.concatenate([acc_win, gyro_win], axis=1)
    windows.append(win)
    labels.append(act_id)
    # return windows, labels

def get_features(list_window):
    features = list()
    for win in list_window:
        max = np.amax(win, axis=0)
        min = np.amin(win, axis=0)
        mean = np.mean(win, axis=0)
        std = np.std(win, axis=0)
        mad = np.median(win, axis=0)
        feature = np.concatenate([max, min, mean, std, mad], axis=0)
        features.append(feature)
    return np.array(features)

# features = get_features(windows)
# print(features.shape)
# print(len(labels))
#
# pickle.dump(features, open('data/features.narray.bin' , 'wb'))
# pickle.dump(np.array(labels), open('data/labels.narray.bin', 'wb'))

features = pickle.load(open('data/features.narray.bin', 'rb'))
print(features.shape)
labels = pickle.load((open('data/labels.narray.bin', 'rb')))
print(labels.shape)

num_train = 1000
x_train = features[:num_train]
y_train = labels[:num_train]
x_test = features[num_train:]
y_test = labels[num_train:]

from sklearn.svm import SVC

model = SVC()
model.fit(x_train, y_train)
print(model.score(x_test, y_test))
predict = model.predict(x_test)
from sklearn.metrics import classification_report, confusion_matrix

print(classification_report(y_test, predict))
lbl = {
    1: 'WALKING',
    2: 'WALKING_UPSTAIRS',
    3: 'WALKING_DOWNSTAIRS',
    4: 'SITTING',
    5: 'STANDING',
    6: 'LAYING',
    7: 'STAND_TO_SIT',
    8: 'SIT_TO_STAND',
    9: 'SIT_TO_LIE',
    10: 'LIE_TO_SIT',
    11: 'STAND_TO_LIE',
    12: 'LIE_TO_STAND'
}

conf_mat = confusion_matrix(y_test, predict)

fig = plt.figure(figsize=(12, 12))
width = np.shape(conf_mat)[1]
height = np.shape(conf_mat)[0]

res = plt.imshow(np.array(conf_mat), cmap=plt.cm.GnBu, interpolation='nearest')
for i, row in enumerate(conf_mat):
    for j, c in enumerate(row):
        if c > 0:
            plt.text(j - .2, i + .1, c, fontsize=16)

cb = fig.colorbar(res)
plt.title('Confusion Matrix')
_ = plt.xticks(range(12), [l for l in lbl.values()], rotation=90)
_ = plt.yticks(range(12), [l for l in lbl.values()])
plt.show()


text:
import numpy as np
import pandas as pd
import pytest

from ..choice import *
from ..sampling import seeded_call

def test_binary_choice():
    p = pd.Series([0, 0.1, 0.9])
    expected_result = pd.Series([False, False, True])

    # 1st test with a provided threshold
    t = 0.5
    c1 = binary_choice(p, t)
    assert (c1 == expected_result).all()

    # next test with randoms (the default)
    # note the randoms produced with this seed should be
    # [0.69646919,  0.28613933,  0.22685145]
    seed = 123
    c2 = seeded_call(seed, binary_choice, p)
    assert (c2 == expected_result).all()

def test_rate_based_binary_choice():
    # expected randoms should be
    # [ 0.69646919  0.28613933  0.22685145  0.55131477]
    seed = 123

    rates = pd.DataFrame({
        'grp1': ['a', 'a', 'b'],
        'grp2': [1, 2, 1],
        'rate': [0.1, 0.6, 0.3]
    })

    agents = pd.DataFrame({
        'grp1': ['a', 'a', 'b', 'b'],
        'grp2': [1, 2, 1, 3]
    })

    expected_result = pd.Series([False, True, True, False])
    c = seeded_call(
        seed,
        rate_based_binary_choice,
        rates, 'rate', agents, ['grp1', 'grp2']
    )
    assert (c == expected_result).all()

def test_logit_binary_choice():
    seed = 123
    agents = pd.DataFrame({
        'do_it': [1, 1, 1, 1, 0, 0, 0],
        'var1': [1, 1, 1, 0, 0, 0, 0],
        'var2': np.arange(7, 0, -1),
        'intercept': np.ones(7)
    })
    coeff = pd.Series(
        [-2, 3, .5],
        index=pd.Index(['intercept', 'var1', 'var2'])
    )
    u, p, c = seeded_call(seed, logit_binary_choice, coeff, agents)

    expected_u = pd.Series([
        90.0171313, 54.59815003, 33.11545196, 1.,
        0.60653066, 0.36787944, 0.22313016
    ])
    assert (u.round(5) == expected_u.round(5)).all()

    expected_p = pd.Series([
        0.98901306, 0.98201379, 0.97068777, 0.5,
        0.37754067, 0.26894142, 0.18242552
    ])
    assert (p.round(5) == expected_p.round(5)).all()

    # expected randoms should be
    # [ 0.69646919  0.28613933  0.22685145  0.55131477
    #   0.71946897  0.42310646 0.9807642 ]
    expected_c = pd.Series([
        True, True, True, False, False, False, False
    ])
    assert (c == expected_c).all()

##########################
# WEIGHTED CHOICE
##########################

@pytest.fixture()
def agents():
    return pd.DataFrame({'col': ['a', 'b', 'c']})

@pytest.fixture()
def alts():
    return pd.DataFrame({
        'w': [1, 5, 10],
        'c': [3, 2, 1]
    })

def test_weighted_choice_noW_noCap(agents, alts):
    c = seeded_call(
        123, weighted_choice,
        agents, alts
    )
    expected_c = pd.Series([2, 1, 2])
    assert (c == expected_c).all()

def test_weighted_choice_noCap(agents, alts):
    c = seeded_call(
        123, weighted_choice,
        agents, alts, 'w'
    )
    expected_c = pd.Series([2, 1, 1])
    assert (c == expected_c).all()

def test_weighted_choice_noW(agents, alts):
    c = seeded_call(
        123, weighted_choice,
        agents, alts, cap_col='c'
    )
    expected_c = pd.Series([0, 1, 1])
    assert (c == expected_c).all()

def test_weighted_choice(agents, alts):
    c = seeded_call(
        123, weighted_choice,
        agents, alts, w_col='w', cap_col='c'
    )
    expected_c = pd.Series([2, 1, 1])
    assert (c == expected_c).all()

def test_weighted_choice_notEnoughCap(agents, alts):
    with pytest.raises(ValueError):
        weighted_choice(agents, alts.loc[2], w_col='w', cap_col='c')

#########################################
# CHOICE WITH SAMPLING OF ALTERNATIVES
#########################################

@pytest.fixture()
def choosers():
    return pd.DataFrame(
        {
            'agent_col1': [10, 20, 30]
        },
        index=pd.Index(list('cba'))
    )

@pytest.fixture()
def alternatives():
    return pd.DataFrame(
        {
            'alt_col1': [100, 200, 300, 400, 500],
            'cap': [2, 2, 0, 1, 1]
        },
        index=pd.Index(np.arange(5, 0, -1,))
    )

def test_get_interaction_data(choosers, alternatives):
    seed = 123
    sample_size = 2

    def check_results(data, num_choosers, size, expected_choosers, exepected_alts):
        for cc in choosers.columns:
            assert cc in data.columns
        for ac in alternatives.columns:
            assert ac in data.columns
        assert len(data) == num_choosers * size

        choosers_idx = data.index.get_level_values(0).values
        alts_idx = data.index.get_level_values(1).values
        assert (choosers_idx == expected_choosers).all()
        assert (alts_idx == expected_alts).all()

    # with replacement
    data, new_size = seeded_call(
        seed, get_interaction_data, choosers, alternatives, sample_size)
    expected_choosers = ['c', 'c', 'b', 'b', 'a', 'a']
    expected_alts = [1, 3, 3, 4, 2, 3]
    assert new_size == sample_size
    check_results(data, len(choosers), new_size, expected_choosers, expected_alts)

    # without replacement -- note not enough alternative to satify sample size
    # so sample size will be adjusted to 1 per agent
    data, new_size = seeded_call(
        seed, get_interaction_data, choosers, alternatives, sample_size, False)
    expected_choosers = ['c', 'b', 'a']
    expected_alts = [4, 2, 1]

def prob_call(interaction_data, num_choosers, sample_size, factor):
    # simple probabilities function that just uses the alt column as a weight
    util = factor * interaction_data['alt_col1'].values.reshape(num_choosers, sample_size)
    return util / util.sum(axis=1, keepdims=True)

def test_choice_with_sampling(choosers, alternatives):
    seed = 123
    sample_size = 4

    # 1st test w/out verbosity
    choices = seeded_call(
        seed,
        choice_with_sampling,
        choosers, alternatives, prob_call, factor=1.0, sample_size=sample_size
    )
    assert (choices['alternative_id'] == [1, 4, 1]).all()

    # now test w/ samples as well
    choices, samples = seeded_call(
        seed,
        choice_with_sampling,
        choosers, alternatives, prob_call, factor=1.0, sample_size=sample_size, verbose=True
    )
    assert (choices['alternative_id'] == [1, 4, 1]).all()
    assert len(samples) == len(choosers) * sample_size
    assert 'alternative_id' in samples.columns
    assert 'chooser_id' in samples.columns
    assert 'prob' in samples.columns

    # test without replacement
    choices, samples = seeded_call(
        seed,
        choice_with_sampling,
        choosers.head(2), alternatives, prob_call,
        sample_size=2, verbose=True, factor=2.0, sample_replace=False
    )
    assert (choices['alternative_id'] == [2, 1]).all()
    assert len(samples) == 4

    # test without replacement, without enough alts
    with pytest.raises(ValueError):
        choices = choice_with_sampling(
            choosers,
            alternatives.head(2),
            prob_call, sample_size=sample_size, factor=2.0, sample_replace=False
        )

def test_capacity_choice_with_sampling(choosers, alternatives):
    seed = 123
    sample_size = 4

    choices, capacities = seeded_call(
        seed,
        capacity_choice_with_sampling,
        choosers,
        alternatives,
        'cap',
        prob_call,
        sample_size,
        factor=2.0
    )
    assert (choices == pd.Series([1, 5, 2], index=choosers.index)).all()
    assert (capacities == pd.Series([1, 2, 0, 0, 0], index=alternatives.index)).all()

def test_capacity_choice_with_sampling_notEnoughCap(choosers, alternatives):
    seed = 123
    sample_size = 4

    choices, capacities = seeded_call(
        seed,
        capacity_choice_with_sampling,
        choosers,
        alternatives.head(1),
        'cap',
        prob_call,
        sample_size,
        factor=2.0
    )
    print choices
    assert choices.loc['c'] == 5
    assert choices.loc['b'] == 5
    assert (choices.isnull().values == [False, False, True]).all()
    assert (capacities == pd.Series([0], index=pd.Index([5]))).all()

def test_capacity_choice_with_sampling_finalI(choosers, alternatives):
    seed = 123
    sample_size = 4

    choices, capacities = seeded_call(
        seed,
        capacity_choice_with_sampling,
        choosers,
        alternatives,
        'cap',
        prob_call,
        sample_size,
        max_iterations=1,
        factor=2.0
    )
    assert (choices == pd.Series([5, 1, 4], index=choosers.index)).all()
    assert (capacities == pd.Series([1, 1, 0, 1, 0], index=alternatives.index)).all()

def test_get_mnl_probs(choosers, alternatives):
    # get some interaction data
    # note this chooses alt IDs [1, 3, 3, 4, 2, 3]
    # resulting alt_col1 values: [500, 300, 300, 200, 400, 300]
    # resulting agent_col1 values: [10, 10, 20, 20, 30, 30]
    sample_size = 2
    data, new_size = seeded_call(
        123, get_interaction_data, choosers, alternatives, sample_size)

    # define some coeeficients
    coeff = pd.Series([0.1, .01], index=pd.Index(['agent_col1', 'alt_col1']))

    # expected results
    exp_p = np.round([0.8807, 0.1192, 0.7310, 0.26894, 0.7310, 0.2689], 2)

    # test w/out an expression
    p1 = get_mnl_probs(data, len(choosers), sample_size, coeff)
    assert (np.round(p1.ravel(), 2) == exp_p).all()

    # test w/ patsy expression
    expr = 'agent_col1 + alt_col1'
    p2 = get_mnl_probs(data, len(choosers), sample_size, coeff, expr)
    assert (np.round(p2.ravel(), 2) == exp_p).all()

    # test w/ an intercept in the coefficients
    coeff.loc['Intercept'] = 1
    p3 = get_mnl_probs(data, len(choosers), sample_size, coeff)
    assert (np.round(p3.ravel(), 2) == exp_p).all

def test_mnl_choice_with_sampling(choosers, alternatives):
    seed = 123
    sample_size = 4

    coeff = pd.Series(
        [0.5, .01],
        index=['agent_col1', 'alt_col1']
    )

    # test 1 - no capacity
    choices = seeded_call(
        seed,
        mnl_choice_with_sampling,
        choosers,
        alternatives,
        coeff
    )
    assert (choices == pd.Series([1, 3, 1], index=choosers.index)).all()

    # test 2 - with capacities, use a patsy expression this time
    expr = 'agent_col1 + alt_col1'
    choices = seeded_call(
        seed,
        mnl_choice_with_sampling,
        choosers,
        alternatives,
        coeff,
        expr,
        sample_size=sample_size,
        cap_col='cap'
    )
    assert (choices == pd.Series([4, 2, 1], index=choosers.index)).all()