text:
# -*- coding: utf-8 -*-

"""
This module provides evaluators that can be used as fitness evaluators in Genens.
"""

from abc import ABC, abstractmethod
from functools import wraps

import logging
import numpy as np

from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.metrics import accuracy_score

from stopit import ThreadingTimeout as Timeout, TimeoutException

import time
import warnings

def timeout(fn):
    @wraps(fn)
    def with_timeout(self, workflow, *args, **kwargs):
        if not hasattr(self, 'timeout') or self.timeout is None:
            return fn(self, workflow, *args, **kwargs)

        try:
            with Timeout(self.timeout, swallow_exc=False):
                res = fn(self, workflow, *args, **kwargs)
        except TimeoutException:
            logger = logging.getLogger("genens")
            logger.debug(f"Timeouted:\n {workflow}")
            print(f"Timeouted:\n {workflow}")

            res = None

        return res
    return with_timeout

def eval_time(fn):
    @wraps(fn)
    def with_time(self, workflow, *args, **kwargs):
        start_time = time.process_time()

        res = fn(self, workflow, *args, **kwargs)
        if res is None or np.isnan(res):
            return None

        elapsed_time = time.process_time() - start_time

        logger = logging.getLogger("genens")
        logger.debug(f"Score: {res}, evaluation time - {elapsed_time}s,\n{workflow}")

        return res, np.log(elapsed_time + np.finfo(float).eps)

    return with_time

def default_score(workflow, test_X, test_y):
    """
    Default scoring method for evaluators.

    :param workflow: Workflow to be evaluated.
    :param test_X: Features.
    :param test_y: Target vector.
    :return: Accuracy score of ``workflow`` on the data.
    """
    res_y = workflow.predict(test_X)
    return accuracy_score(test_y, res_y)

def _simple_eval(workflow, train_X, train_y, test_X, test_y, scorer=None):
    workflow.fit(train_X, train_y)
    if scorer is not None:
        return scorer(workflow, test_X, test_y)

    return default_score(workflow, test_X, test_y)

class EvaluatorBase(ABC):
    """
    Evaluator interface which is used in Genens during individual evaluation.
    """
    def __init__(self, timeout_s=None):
        self.train_X = None
        self.train_y = None

        self.timeout = timeout_s

    def __repr__(self):
        res = "{}".format(self.__class__.__name__)
        return res

    def fit(self, train_X, train_y):
        """
        Sets the training data for evaluation.

        :param train_X: Features.
        :param train_y: Target vector.
        """
        self.train_X = train_X
        self.train_y = train_y

    @abstractmethod
    def evaluate(self, workflow, scorer=None):
        """
        Evaluates a pipeline (workflow) learned on the training data from this evaluator. Evaluation method
        differs between different evaluators.
        :param workflow: Workflow to be evaluated.
        :param scorer: Scorer to be used, default is the accuracy scorer.
        :return: Score of the workflow.
        """
        pass

    @abstractmethod
    def reset(self):
        pass

    def check_is_fitted(self):
        if self.train_X is None or self.train_y is None:
            raise ValueError("Evaluator is not fitted with training data.")  # TODO specific

    @timeout
    @eval_time
    def score(self, workflow, scorer=None):
        """
        Computes the score of a pipeline (workflow).
        :param workflow: Workflow to be evaluated.
        :param scorer: Scorer to be used, default is the accuracy scorer.
        :return: Score of the workflow.
        """
        self.check_is_fitted()

        try:
            with warnings.catch_warnings():
                warnings.simplefilter('ignore')

                return self.evaluate(workflow, scorer=scorer)
        except Exception as e:
            logger = logging.getLogger("genens")
            logger.debug(f"Workflow failed:\n {e}\nWorkflow: {workflow}.")
            return None

class CrossValEvaluator(EvaluatorBase):
    """
    Evaluator which performs a k-fold cross-validation.
    """
    def __init__(self, cv_k=7, timeout_s=None):
        super().__init__(timeout_s)

        if cv_k < 0:
            raise AttributeError("Cross validation k must be greater than 0.")

        self.cv_k = cv_k

    def __repr__(self):
        res = super().__repr__()
        res += ", cv_k: {}".format(self.cv_k)
        return res

    def evaluate(self, workflow, scorer=None):
        scores = cross_val_score(workflow, self.train_X, self.train_y,
                                 cv=self.cv_k, scoring=scorer, error_score='raise')
        return np.mean(scores)

    def reset(self):
        pass

class FixedTrainTestEvaluator(EvaluatorBase):
    """
    Evaluator that uses fixed train and test samples.
    """
    def __init__(self, test_size=None, random_state=None, timeout_s=None):
        super().__init__(timeout_s)

        self.test_X = None
        self.test_y = None

        self.test_size = test_size
        self.random_state = random_state

    def __repr__(self):
        res = super().__repr__()
        res += ", test_size: {}".format(self.test_size)
        return res

    def fit(self, train_X, train_y):
        self.train_X, self.test_X, self.train_y, self.test_y = \
            train_test_split(train_X, train_y, test_size=self.test_size,
                             random_state=self.random_state)

    def check_is_fitted(self):
        super().check_is_fitted()
        if self.test_X is None or self.test_y is None:
            raise ValueError("Test data missing in evaluator.")

    def evaluate(self, workflow, scorer=None):
        return _simple_eval(workflow, self.train_X, self.train_y, self.test_X, self.test_y,
                            scorer=scorer)

    def reset(self):
        pass

class RandomTrainTestEvaluator(EvaluatorBase):
    """
    Evaluator which performs the train-test split for every evaluation.
    """
    def __init__(self, test_size=None, random_state=None, timeout_s=None):
        super().__init__(timeout_s)

        self.test_size = test_size

        self.rng = None
        if random_state is not None:
            if isinstance(random_state, np.random.RandomState):
                self.rng = random_state

            self.rng = np.random.RandomState(random_state)

    def __repr__(self):
        res = super().__repr__()
        res += ", test_size: {}".format(self.test_size)
        return res

    def evaluate(self, workflow, scorer=None):
        # random state is set only in the constructor
        train_X, test_X, train_y, test_y = train_test_split(self.train_X, self.train_y,
                                                            test_size=self.test_size,
                                                            random_state=self.rng)
        return _simple_eval(workflow, train_X, train_y, test_X, test_y,
                            scorer=scorer)

    def reset(self):
        pass

class TrainTestEvaluator(EvaluatorBase):
    """
    Evaluator with a specified validation set.
    """
    def __init__(self, test_X, test_y, timeout_s=None):
        super().__init__(timeout_s)

        # training set is provided in fit, validation set is set on initialization
        self.test_X = test_X
        self.test_y = test_y

    def evaluate(self, workflow, scorer=None):
        return _simple_eval(workflow, self.train_X, self.train_y, self.test_X, self.test_y,
                            scorer=scorer)

    def reset(self):
        pass

class DataSampler:
    """
    Provides random samples of a given dataset.
    """
    def __init__(self, sample_size=0.20, random_state=None, stratified=True):
        self.full_X = None
        self.full_y = None

        self.sample_size = sample_size
        self.stratified = stratified

        self.rng = None
        if random_state is not None:
            if isinstance(random_state, np.random.RandomState):
                self.rng = random_state

            self.rng = np.random.RandomState(random_state)

    def fit(self, full_X, full_y):
        self.full_X = full_X
        self.full_y = full_y

    def generate_sample(self):
        """
        Generates a random sample from full data.
        :return: A random sample.
        """
        stratify = self.full_y if self.stratified else None

        _, sample_X, _, sample_y = train_test_split(self.full_X, self.full_y,
                                                    test_size=self.sample_size,
                                                    random_state=self.rng,
                                                    stratify=stratify)
        return sample_X, sample_y

class SampleCrossValEvaluator(CrossValEvaluator):
    """
    Evaluator which uses a k-fold cross-validation based on sampling.
    """
    def __init__(self, cv_k=7, timeout_s=None, sample_size=0.20, per_gen=False, random_state=None):
        super().__init__(cv_k=cv_k, timeout_s=timeout_s)

        self.per_gen = per_gen
        self.sampler = DataSampler(sample_size=sample_size, random_state=random_state)

    def __repr__(self):
        res = super().__repr__()
        res += ", sample_size: {}".format(self.sampler.sample_size)
        res += ", per_gen: {}".format(self.per_gen)
        return res

    def _reset_data(self):
        self.train_X, self.train_y = self.sampler.generate_sample()

    def fit(self, train_X, train_y):
        self.sampler.fit(train_X, train_y)
        self._reset_data()

    def evaluate(self, workflow, scorer=None):
        # generate sample for every evaluation
        if not self.per_gen:
            self._reset_data()

        return super().evaluate(workflow, scorer=scorer)

    def reset(self):
        # generate sample once per reset
        if self.per_gen:
            self._reset_data()

class SampleTrainTestEvaluator(FixedTrainTestEvaluator):
    """
    Evaluator which uses a train-test evaluation on samples.
    """
    def __init__(self, test_size=None, timeout_s=None, sample_size=0.20,
                 per_gen=False, random_state=None):

        self.per_gen = per_gen
        self.sampler = DataSampler(sample_size=sample_size, random_state=random_state)

        # rng should be set in DataSampler ctor
        super().__init__(test_size=test_size, random_state=self.sampler.rng, timeout_s=timeout_s)

    def __repr__(self):
        res = super().__repr__()
        res += ", sample_size: {}".format(self.sampler.sample_size)
        res += ", per_gen: {}".format(self.per_gen)
        return res

    def _fit_sample(self):
        train_X, train_y = self.sampler.generate_sample()

        super().fit(train_X, train_y)

    def fit(self, train_X, train_y):
        self.sampler.fit(train_X, train_y)

        self._fit_sample()

    def evaluate(self, workflow, scorer=None):
        # one sample per evaluation
        if not self.per_gen:
            self._fit_sample()

        return super().evaluate(workflow, scorer=scorer)

    def reset(self):
        # generate new fixed sample
        if self.per_gen:
            self._fit_sample()

_eval_names = {
    'crossval': CrossValEvaluator,
    'fixed': FixedTrainTestEvaluator,
    'per_ind': RandomTrainTestEvaluator,
    'train_test': TrainTestEvaluator,
    'sample_crossval': SampleCrossValEvaluator,
    'sample_train_test': SampleTrainTestEvaluator
}

def get_evaluator_cls(name):
    return _eval_names[name]

max_stars_repo_name: gabrielasuchopar/genens

max_issues_repo_name: gabrielasuchopar/genens

max_forks_repo_name: gabrielasuchopar/genens


text:
#!/usr/bin/python
#https://www.hackerrank.com/challenges/np-polynomials/problem

import numpy
l1 = list(map(float,input().split()))
x = int(input())
print(numpy.polyval(l1,x))

ext: py

reason: import numpy


text:
import numpy

# The formula is in the 'Tutorial' tab
# Matrix multiplication operator '@' was introduced in Python 3.5.0

if __name__ == '__main__':
    n = int(input().strip().split(' ')[1])

    X, Y = [], []

    for _ in range(n):
        arr = list(map(float, input().strip().split(' ')))
        X.append([1] + arr[0: -1])
        Y.append(arr[-1])

    X, Y = map(numpy.array, (X, Y))

    XT = numpy.transpose(X)

    B = numpy.linalg.inv(XT @ X) @ XT @ Y

    q = int(input())

    X = []
    for _ in range(q):
        X.append([1] + list(map(float, input().strip().split(' '))))

    X = numpy.array(X)

    Y = X @ B

    print('\n'.join(('{:.2f}'.format(value) for value in Y)))

max_stars_repo_path: 10 Days of Statistics/Day 9: Multiple Linear Regression/solution.py

max_stars_repo_name: filipjeremic/hackerrank

max_issues_repo_name: filipjeremic/hackerrank

max_forks_repo_name: filipjeremic/hackerrank


text:
import pickle
import json
import numpy as np

__columns=None
__location=None
__model=None

def get_location_names():
    return __location

def predict_price(bhk,total_sqrt,bath,location):
    try:
        loc=__columns.index(location.lower())
    except:
        loc=__columns.index('Whitefield'.lower())

    p=np.zeros(len(__columns))
    p[0],p[1],p[2]=bhk,total_sqrt,bath
    p[loc]=1
    return round(__model.predict([p])[0],2)

def load_essentials():
    global __columns
    global __location
    global __model
    with open('columns.json','r') as f:
        __columns=json.load(f)['data_columns']
        __location=__columns[3:]
    with open('banglore_home_price_prediction.pickle','rb') as f:
        __model=pickle.load(f)

    print('model,columns and location loaded')
load_essentials()
if __name__=='__main__':
    print(predict_price(4,2825,4,'1st Phase JP Nagar'))
    print(get_location_names())


text:
import sys
import IPython
import pathlib
import machine_learning as ml
import machine_learning.agnostic as agn
import pandas as pd
import numpy as np

data_path = pathlib.Path(__file__).parents[1] / "data"

def test_fit():
    """
    Montgomery, D. C., & Runger, G. C. (2010). Applied statistics and
    probability for engineers. John Wiley & Sons.
    Example 12-1
    """
    df = pd.read_csv(data_path / "WireBondData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2"]
    model = ml.LinearModel()
    model.fit(df[x], df[y])
    res = "y = 2.26379 + 2.74427*X1 + 0.0125278*X2"
    assert model.get_model_string(precision=6) == res

def test_weighted():
    df = pd.read_csv(data_path / "WireBondData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2"]
    model = ml.LinearModel()
    model.fit(df[x], df[y], weights=df["X4"])
    res = "y = 1.08305 + 2.92825*X1 + 0.0112528*X2"
    assert model.get_model_string(precision=6) == res

def test_without_bias():
    df = pd.read_csv(data_path / "WireBondData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2"]
    model = ml.LinearModel()
    model.fit(df[x], df[y], bias=False)
    res = "y = + 2.82327*X1 + 0.0161177*X2"
    assert model.get_model_string(precision=6) == res

def test_predict():
    df = pd.read_csv(data_path / "WireBondData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2"]
    model = ml.LinearModel()
    model.fit(df[x], df[y])
    X = df[x]
    Y_true = [8.38, 25.60, 33.95, 36.60, 27.91, 15.75, 12.45, 8.40, 28.21,
              27.98, 18.40, 37.46, 41.46, 12.26, 15.81, 18.25, 64.67, 12.34,
              36.47, 46.56, 47.06, 52.56, 56.31, 19.98, 21.00]
    Y_predicted = np.round(model.predict(X), decimals=2)
    assert (Y_true == Y_predicted).all()

def test_Ttest():
    """
    Montgomery, D. C., & Runger, G. C. (2010). Applied statistics and
    probability for engineers. John Wiley & Sons.
    Example 12-4
    """
    df = pd.read_csv(data_path / "WireBondData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2"]
    model = ml.LinearModel()
    p, T = model._T_test(df[x], df[y], "X2")
    assert (round(T, 3), round(p, 4)) == (4.477, 0.0002)

def test_Ftest():
    """
    Montgomery, D. C., & Runger, G. C. (2010). Applied statistics and
    probability for engineers. John Wiley & Sons.
    Example 12-6
    """
    df = pd.read_csv(data_path / "WireBondData.csv", sep=";", decimal=",")
    y = "Y"
    x_full = ["X1", "X2", "X3", "X4"]
    x_sub = ["X1", "X2"]
    model = ml.LinearModel()
    p, F = model._partial_F_test(df, df[y], x_full, x_sub)
    assert (round(F, 3), round(p, 3)) == (4.047, 0.033)

def test_forward():
    df = pd.read_csv(data_path / "WineQualityData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2", "X3", "X4", "X5"]
    model = ml.LinearModel()
    model.fit_forward(df[x], df[y], p=0.1, verbose=False)
    model_str = model.get_model_string(precision=5)
    assert model_str == "y = 6.4672 + 1.1997*X4 - 0.60232*X5 + 0.58012*X2"

def test_backward():
    """
    Montgomery, D. C., & Runger, G. C. (2010). Applied statistics and
    probability for engineers. John Wiley & Sons.
    Example 12-15
    """
    df = pd.read_csv(data_path / "WineQualityData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2", "X3", "X4", "X5"]
    model = ml.LinearModel()
    model.fit_backward(df[x], df[y], p=0.15, verbose=False)
    model_str = model.get_model_string(precision=5)
    assert model_str == "y = 6.4672 + 0.58012*X2 + 1.1997*X4 - 0.60232*X5"

def test_stepwise():
    """
    Montgomery, D. C., & Runger, G. C. (2010). Applied statistics and
    probability for engineers. John Wiley & Sons.
    Example 12-15
    """
    df = pd.read_csv(data_path / "WineQualityData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2", "X3", "X4", "X5"]
    model = ml.LinearModel()
    model.fit_stepwise(df[x], df[y], p=0.1, verbose=False)
    model_str = model.get_model_string(precision=5)
    assert model_str == "y = 6.4672 + 1.1997*X4 - 0.60232*X5 + 0.58012*X2"

def test_confidence():
    """
    Montgomery, D. C., & Runger, G. C. (2010). Applied statistics and
    probability for engineers. John Wiley & Sons.
    Example 12-9
    """
    df = pd.read_csv(data_path / "WireBondData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2"]
    model = ml.LinearModel()
    model.fit(df[x], df[y])
    x = pd.DataFrame(data=[[8., 275.]], columns=["X1", "X2"])
    dy = model.uncertainty(x, alpha=0.05)
    assert round(dy[0], 2) == round(27.66-22.81, 2)

def test_dump():
    df = pd.read_csv(data_path / "WineQualityData.csv", sep=";", decimal=",")
    y = "Y"
    x = ["X1", "X2"]
    model = ml.LinearModel()
    model.fit(df[x], df[y])
    loaded = agn.LinearModel(model.dump)
    eval_model = model.predict(df)
    eval_loaded = loaded(df)
    assert np.isclose(eval_model, eval_loaded, atol=0.).all()

if __name__ == "__main__":
    module = sys.modules[__name__]
    for attr in dir(module):
        if "test_" not in attr:
            continue
        func = getattr(module, attr)
        print(attr)
        func()
    IPython.embed()

max_stars_repo_name: BFavier/pygmalion

max_issues_repo_name: BFavier/pygmalion

max_forks_repo_name: BFavier/pygmalion


text:
import streamlit as st
import os
import requests
import altair as alt
import numpy as np
import pandas as pd
from dotenv import find_dotenv, load_dotenv

def fetch_data_from_aq(city='Calgary', api_key=''):
    query = f'https://api.waqi.info/feed/{city}/?token={api_key}'
    data = requests.get(url=query,
                        params={'token': api_key})

    return data.json()

# load env variables
load_dotenv(find_dotenv())
AQ_TOKEN = os.getenv('AQ_API_KEY', None)
N_IMG = int(os.getenv('N_IMG', "10"))

emoji_json = requests.get(url="https://api.github.com/emojis").json()
img_source = pd.DataFrame({"name": list(emoji_json.keys()), "url": list(emoji_json.values()),
                           'x': np.random.uniform(-50, 50, len(emoji_json)),
                           'y': np.random.uniform(-50, 50, len(emoji_json))})
img_subset =img_source.sample(N_IMG)
ch = alt.Chart(img_subset,width=500,height=500).mark_image(width=20,height=20).encode(
    x='x',
    y='y',
    url='url'
)

st.title('Air Quality check webapp')

st.markdown('This app is a skeleton for what my SPE 2021 Data Science mentees will work on.'
            'Specifically, we would like to compare the air quality levels pre/post COVID in a year-over-year plot '
            'as a baseline')
# get the api key
data = fetch_data_from_aq(city='Houston', api_key=AQ_TOKEN)
st.altair_chart(ch)
st.write(data)


text:
import numpy as np
import pandas as pd
from sklearn.externals import joblib
from matplotlib import pyplot as plt
from sklearn.neural_network import MLPClassifier

# Load loss function values thought training
# Load neural network
ann = joblib.load('../model/mlmodel.pkl')
# Plot Loss function values
fig, ax = plt.subplots()
ax.plot(ann.loss_curve_, 'r', label='Training set')
# plt.style.use('ggplot')
legend = ax.legend(loc=0)
plt.xlabel("Iteration")
plt.ylabel("Mean square error")
plt.grid()
# plt.xlim()
plt.title('Loss Function')
plt.savefig('../graphs/loss_curve_2d', bbox_inches='tight')

max_stars_repo_name: miguelggaspar/neuronalFEM

max_issues_repo_name: miguelggaspar/neuronalFEM

max_forks_repo_name: miguelggaspar/neuronalFEM


text:
"""Contains functions regarding 3D rotation matrices. Used in `scenario.py`."""
import numpy as np

def rot_x(alpha):
    """Returns rotation matrix."""
    return np.asarray(
        [
            [1, 0, 0],
            [0, np.cos(alpha), -np.sin(alpha)],
            [0, np.sin(alpha), np.cos(alpha)]
        ]
    )

def rot_y(alpha):
    """Returns rotation matrix."""
    return np.asarray(
        [
            [np.cos(alpha), 0, np.sin(alpha)],
            [0, 1, 0],
            [-np.sin(alpha), 0, np.cos(alpha)]
        ]
    )

def rot_z(alpha):
    """Returns rotation matrix."""
    return np.asarray(
        [
            [np.cos(alpha), -np.sin(alpha), 0],
            [np.sin(alpha), np.cos(alpha), 0],
            [0, 0, 1]
        ]
    )

max_stars_repo_name: JanekEbb/paderbox

max_issues_repo_name: JanekEbb/paderbox

max_forks_repo_name: JanekEbb/paderbox


text:
import numpy as np

class MeanTracker(object):
    """ simple maximum likelihood mean tracker"""

    def __init__(self, mu_prior=25, n_arms=8):

        self.mu_prior = mu_prior
        self.n_arms = n_arms

        self.means = mu_prior * np.ones((0, n_arms))
        self.counts = np.zeros((0, n_arms), dtype=int)
        self.total_reward = np.zeros((0, n_arms), dtype=float)

        self.visited_blocks = dict()

    def update(self, block, arm, reward):
        if block not in self.visited_blocks.keys():
            self.visited_blocks[block] = len(self.visited_blocks)

            self.means = np.concatenate([self.counts, self.mu_prior * np.ones((1, self.n_arms), dtype=int)])
            self.counts = np.concatenate([self.counts, np.zeros((1, self.n_arms), dtype=int)])
            self.total_reward = np.concatenate([self.total_reward, np.zeros((1, self.n_arms), dtype=int)])

        k = self.visited_blocks[block]

        self.counts[k, arm] += 1
        self.total_reward[k, arm] += reward
        self.means[k, arm] = self.total_reward[k, arm] / self.counts[k, arm]

    def get_mean(self, block, arm):
        if block not in self.visited_blocks.keys():
            return self.mu_prior
        return self.means[self.visited_blocks[block], arm]

max_stars_repo_name: nicktfranklin/StructuredBandits

max_issues_repo_name: nicktfranklin/StructuredBandits

max_forks_repo_name: nicktfranklin/StructuredBandits


text:
# -*- coding: utf-8 -*-
# This software is open source software available under the BSD-3 license.
#
# Copyright (c) 2020 Triad National Security, LLC. All rights reserved.
# Copyright (c) 2020 Lawrence Livermore National Security, LLC. All rights
# reserved.
# Copyright (c) 2020 UT-Battelle, LLC. All rights reserved.
#
# Additional copyright and license information can be found in the LICENSE file
# distributed with this code, or at
# https://raw.githubusercontent.com/MPAS-Dev/MPAS-Analysis/master/LICENSE
#

from __future__ import absolute_import, division, print_function, \
    unicode_literals

import xarray as xr
import numpy
import matplotlib.pyplot as plt

from mpas_tools.cime.constants import constants as cime_constants

from mpas_analysis.shared import AnalysisTask

from mpas_analysis.ocean.compute_anomaly_subtask import ComputeAnomalySubtask
from mpas_analysis.ocean.plot_hovmoller_subtask import PlotHovmollerSubtask
from mpas_analysis.ocean.plot_depth_integrated_time_series_subtask import \
    PlotDepthIntegratedTimeSeriesSubtask

from mpas_analysis.shared.constants import constants as mpas_constants

class TimeSeriesOHCAnomaly(AnalysisTask):
    """
    Performs analysis of ocean heat content (OHC) from time-series output.
    """
    # Authors
    # -------
    # Xylar Asay-Davis, Milena Veneziani, Greg Streletz

    def __init__(self, config, mpasTimeSeriesTask, controlConfig=None):
        # {{{
        """
        Construct the analysis task.

        Parameters
        ----------
        config :  ``MpasAnalysisConfigParser``
            Configuration options

        mpasTimeSeriesTask : ``MpasTimeSeriesTask``
            The task that extracts the time series from MPAS monthly output

        controlConfig :  ``MpasAnalysisConfigParser``, optional
            Configuration options for a control run (if any)
        """
        # Authors
        # -------
        # Xylar Asay-Davis

        # first, call the constructor from the base class (AnalysisTask)
        super(TimeSeriesOHCAnomaly, self).__init__(
            config=config,
            taskName='timeSeriesOHCAnomaly',
            componentName='ocean',
            tags=['timeSeries', 'ohc', 'publicObs'])

        sectionName = 'timeSeriesOHCAnomaly'
        regionNames = config.getExpression(sectionName, 'regions')
        movingAveragePoints = config.getint(sectionName, 'movingAveragePoints')

        self.variableDict = {}
        for suffix in ['avgLayerTemperature', 'sumLayerMaskValue',
                       'avgLayerArea', 'avgLayerThickness']:
            key = 'timeMonthly_avg_avgValueWithinOceanLayerRegion_' + suffix
            self.variableDict[key] = suffix

        mpasFieldName = 'ohc'

        timeSeriesFileName = 'regionAveragedOHCAnomaly.nc'

        anomalyTask = ComputeAnomalySubtask(
            parentTask=self,
            mpasTimeSeriesTask=mpasTimeSeriesTask,
            outFileName=timeSeriesFileName,
            variableList=list(self.variableDict.keys()),
            movingAveragePoints=movingAveragePoints,
            alter_dataset=self._compute_ohc)
        self.add_subtask(anomalyTask)

        for regionName in regionNames:
            caption = 'Trend of {} OHC Anomaly vs depth'.format(
                regionName)
            plotTask = PlotHovmollerSubtask(
                parentTask=self,
                regionName=regionName,
                inFileName=timeSeriesFileName,
                outFileLabel='ohcAnomalyZ',
                fieldNameInTitle='OHC Anomaly',
                mpasFieldName=mpasFieldName,
                unitsLabel=r'[$\times 10^{22}$ J]',
                sectionName='hovmollerOHCAnomaly',
                thumbnailSuffix=u'ΔOHC',
                imageCaption=caption,
                galleryGroup='Trends vs Depth',
                groupSubtitle=None,
                groupLink='trendsvsdepth',
                galleryName=None)

            plotTask.run_after(anomalyTask)
            self.add_subtask(plotTask)

            caption = 'Running Mean of the Anomaly in {} Ocean Heat ' \
                'Content'.format(regionName)
            plotTask = PlotOHCAnomaly(
                parentTask=self,
                regionName=regionName,
                inFileName=timeSeriesFileName,
                outFileLabel='ohcAnomaly',
                fieldNameInTitle='OHC Anomaly',
                mpasFieldName=mpasFieldName,
                yAxisLabel=r'$\Delta$OHC [$\times 10^{22}$ J]',
                sectionName='timeSeriesOHCAnomaly',
                thumbnailSuffix=u'ΔOHC',
                imageCaption=caption,
                galleryGroup='Time Series',
                groupSubtitle=None,
                groupLink='timeseries',
                galleryName=None,
                controlConfig=controlConfig)

            plotTask.run_after(anomalyTask)
            self.add_subtask(plotTask)

        # }}}

    def _compute_ohc(self, ds):  # {{{
        '''
        Compute the OHC time series.
        '''

        # regionNames = self.config.getExpression('regions', 'regions')
        # ds['regionNames'] = ('nOceanRegionsTmp', regionNames)

        # for convenience, rename the variables to simpler, shorter names
        ds = ds.rename(self.variableDict)

        # specific heat [J/(kg*degC)]
        cp = self.namelist.getfloat('config_specific_heat_sea_water')
        # [kg/m3]
        rho = self.namelist.getfloat('config_density0')

        unitsScalefactor = 1e-22

        ds['ohc'] = unitsScalefactor * rho * cp * ds['sumLayerMaskValue'] * \
            ds['avgLayerArea'] * ds['avgLayerThickness'] * \
            ds['avgLayerTemperature']
        ds.ohc.attrs['units'] = '$10^{22}$ J'
        ds.ohc.attrs['description'] = 'Ocean heat content in each region'

        # Note: restart file, not a mesh file because we need refBottomDepth,
        # not in a mesh file
        try:
            restartFile = self.runStreams.readpath('restart')[0]
        except ValueError:
            raise IOError('No MPAS-O restart file found: need at least one '
                          'restart file for OHC calculation')

        # Define/read in general variables
        with xr.open_dataset(restartFile) as dsRestart:
            # reference depth [m]
            # add depths as a coordinate to the data set
            ds.coords['depth'] = (('nVertLevels',),
                                  dsRestart.refBottomDepth.values)

        return ds  # }}}

    # }}}

class PlotOHCAnomaly(PlotDepthIntegratedTimeSeriesSubtask):
    def customize_fig(self, fig):
        """
        A function to override to customize the figure.

        fig : matplotlib.pyplot.Figure
            The figure
        """
        def joules_to_watts_m2(joules):
            watts_m2 = joules/factor
            return watts_m2

        def watts_m2_to_joules(watts_m2):
            joules = factor*watts_m2
            return joules

        # add an axis on the right-hand side
        color = 'tab:blue'
        ax = plt.gca()
        xlim = ax.get_xlim()

        earth_surface_area = (4. * numpy.pi *
                              cime_constants['SHR_CONST_REARTH']**2)

        max_time = xlim[-1]*mpas_constants.sec_per_day

        factor = earth_surface_area*max_time/10**22

        secaxy = ax.secondary_yaxis(
            'right', functions=(joules_to_watts_m2, watts_m2_to_joules))
        secaxy.set_ylabel(r'W/m$^2$', color=color)
        secaxy.tick_params(axis='y', colors=color)
        ax.spines['right'].set_color(color)
        plt.draw()
        yticks = secaxy.get_yticks()
        for ytick in yticks:
            plt.plot(xlim, [0, watts_m2_to_joules(ytick)], color=color,
                     linewidth=0.5)

# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python


text:
import math, unittest
import numpy
from pysparse.tools import poisson
from pysparse.sparse.pysparseMatrix import PysparseMatrix, PysparseIdentityMatrix, PysparseSpDiagsMatrix
from pysparse.direct.pysparseSuperLU import PysparseSuperLUSolver

import sys

def macheps():
    "compute machine epsilon"
    eps = 1.0
    while (1.0 + eps > 1.0):
        eps /= 2.0
    return 2.0 * eps

#
# Tests that are commented out work fine but take a long time...
#

class SpDiagsTestCase(unittest.TestCase):
    def setUp(self):
        self.n = 10000
        self.e = numpy.ones(self.n)
        self.A = PysparseSpDiagsMatrix(size=self.n,
                                       vals=(-2*self.e,self.e,2*self.e),
                                       pos=(-1,0,1))
        self.b = self.A * self.e
        self.tol = 100*macheps()
        self.LU = None
        self.relerr = 0.0
        self.descr = ''
        self.fmt = '\t%10s  %8.2e  %8.2e  %8d  %8d  %6.2f  %6.2f\n'

    def tearDown(self):
        self.LU.fetch_lunz()
        sys.stdout.write(self.fmt % (self.descr, self.relerr, self.tol,
                                     self.A.getNnz(),
                                     self.LU.lunz, self.LU.factorizationTime,
                                     self.LU.solutionTime))

    def computeError(self, x):
        self.relerr = numpy.linalg.norm(x - self.e, ord=numpy.inf)
        return self.relerr

    #def testVanilla(self):
    #    self.descr = 'spdgs-dflt'
    #    self.LU = PysparseSuperLUSolver(self.A)
    #    self.LU.solve(self.b)
    #    self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testTrivialThresh(self):
        self.descr = 'spdgs-trsh'
        self.LU = PysparseSuperLUSolver(self.A, diag_pivot_thresh=0.5)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    #def testTrivialRelax(self):
    #    self.descr = 'spdgs-relx'
    #    self.LU = PysparseSuperLUSolver(self.A, relax=20)
    #    self.LU.solve(self.b)
    #    self.failUnless(self.computeError(self.LU.sol) < self.tol)

    #def testTrivialPanel(self):
    #    self.descr = 'spdgs-size'
    #    self.LU = PysparseSuperLUSolver(self.A, panel_size=1)
    #    self.LU.solve(self.b)
    #    self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testTrivialperm0(self):
        self.descr = 'spdgs-prm0'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=0)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testTrivialperm1(self):
        self.descr = 'spdgs-prm1'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=1)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    #def testTrivialperm2(self):
    #    self.descr = 'spdgs-prm2'
    #    self.LU = PysparseSuperLUSolver(self.A, permc_spec=2)
    #    self.LU.solve(self.b)
    #    self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testTrivialperm3(self):
        self.descr = 'spdgs-prm3'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=3)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

class Poisson1dTestCase(unittest.TestCase):

    def setUp(self):
        self.n = 50000
        self.A = PysparseMatrix(matrix=poisson.poisson1d_sym(self.n))

        self.x_exact = numpy.ones(self.n)/math.sqrt(self.n)
        self.normx = 1.0/math.sqrt(self.n)
        self.b = self.A * self.x_exact

        lmbd_min = 4.0 * math.sin(math.pi/2.0/self.n) ** 2
        lmbd_max = 4.0 * math.sin((self.n - 1)*math.pi/2.0/self.n) ** 2
        cond = lmbd_max/lmbd_min
        self.tol = cond * macheps()
        self.LU = None
        self.relerr = 0.0
        self.descr = ''
        self.fmt = '\t%10s  %8.2e  %8.2e  %8d  %8d  %6.2f  %6.2f\n'

    def tearDown(self):
        self.LU.fetch_lunz()
        sys.stdout.write(self.fmt % (self.descr, self.relerr, self.tol,
                                     self.A.getNnz(),
                                     self.LU.lunz, self.LU.factorizationTime,
                                     self.LU.solutionTime))

    def computeError(self, x):
        absErr = numpy.linalg.norm(x-self.x_exact, ord=numpy.inf)
        self.relerr = absErr/(1+self.normx)
        return self.relerr

    def testPoisson1dDefault(self):
        self.descr = 'poi1d-dflt'
        self.LU = PysparseSuperLUSolver(self.A)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson1dThresh(self):
        self.descr = 'poi1d-trsh'
        self.LU = PysparseSuperLUSolver(self.A, diag_pivot_thresh=0.5)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson1dRelax(self):
        self.descr = 'poi1d-relx'
        self.LU = PysparseSuperLUSolver(self.A, relax=20)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson1dPanel(self):
        self.descr = 'poi1d-size'
        self.LU = PysparseSuperLUSolver(self.A, panel_size=1)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson1dperm0(self):
        self.descr = 'poi1d-prm0'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=0)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson1dperm1(self):
        self.descr = 'poi1d-prm1'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=1)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson1dperm2(self):
        self.descr = 'poi1d-prm2'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=2)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson1dperm3(self):
        self.descr = 'poi1d-prm3'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=3)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

class Poisson2dTestCase(unittest.TestCase):
    def setUp(self):
        self.n = 200
        self.A = PysparseMatrix(matrix=poisson.poisson2d_sym_blk(self.n))

        self.x_exact = numpy.ones(self.n*self.n)/self.n
        self.normx = 1.0/self.n
        self.b = self.A * self.x_exact

        h = 1.0 / self.n
        lmbd_min = 4.0/h/h * (math.sin(math.pi*h/2.0) ** 2 +
                              math.sin(math.pi*h/2.0) ** 2)
        lmbd_max = 4.0/h/h * (math.sin((self.n - 1)*math.pi*h/2.0) ** 2 +
                              math.sin((self.n - 1)*math.pi*h/2.0) ** 2)
        cond = lmbd_max/lmbd_min
        self.tol = cond * macheps()
        self.LU = None
        self.relerr = 0.0
        self.descr = ''
        self.fmt = '\t%10s  %8.2e  %8.2e  %8d  %8d  %6.2f  %6.2f\n'

    def tearDown(self):
        self.LU.fetch_lunz()
        sys.stdout.write(self.fmt % (self.descr, self.relerr, self.tol,
                                     self.A.getNnz(),
                                     self.LU.lunz, self.LU.factorizationTime,
                                     self.LU.solutionTime))

    def computeError(self, x):
        absErr = numpy.linalg.norm(x-self.x_exact, ord=numpy.inf)
        self.relerr = absErr/(1+self.normx)
        return self.relerr

    def testPoisson2dDefault(self):
        self.descr = 'poi2d-dftl'
        self.LU = PysparseSuperLUSolver(self.A)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson2dThresh(self):
        self.descr = 'poi2d-trsh'
        self.LU = PysparseSuperLUSolver(self.A, diag_pivot_thresh=0.5)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson2dRelax(self):
        self.descr = 'poi2d-relx'
        self.LU = PysparseSuperLUSolver(self.A, relax=20)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson2dPanel(self):
        self.descr = 'poi2d-size'
        self.LU = PysparseSuperLUSolver(self.A, panel_size=1)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson2dperm0(self):
        self.descr = 'poi2d-prm0'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=0)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson2dperm1(self):
        self.descr = 'poi2d-prm1'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=1)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson2dperm2(self):
        self.descr = 'poi2d-prm2'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=2)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

    def testPoisson2dperm3(self):
        self.descr = 'poi2d-prm3'
        self.LU = PysparseSuperLUSolver(self.A, permc_spec=3)
        self.LU.solve(self.b)
        self.failUnless(self.computeError(self.LU.sol) < self.tol)

if __name__ == '__main__':
    headfmt = '\t%10s  %8s  %8s  %8s  %8s  %6s  %6s'
    header = headfmt % ('Test','RelErr','Tol','nnz(A)',
                        'nnz(L+U)','Fact','Solve')
    lhead = len(header)
    sys.stderr.write(header + '\n')
    sys.stderr.write('\t' + '-' * lhead + '\n')
    unittest.main()

max_stars_repo_name: heeryerate/pysparse

max_issues_repo_name: heeryerate/pysparse

max_forks_repo_name: heeryerate/pysparse


text:
from __future__ import annotations

import logging
from textwrap import dedent
from typing import Union
import warnings

import numpy as np
import param
import xarray as xr

from GSForge._singledispatchmethod import singledispatchmethod
from ._AnnotatedGEM import AnnotatedGEM
from ._GeneSetCollection import GeneSetCollection
from ..utils import transient_log_handler

logger = logging.getLogger("GSForge")

# TODO: Add .keys() and other dict like functionality.
class Interface(param.Parameterized):
    """The Interface provides common API access for interacting with the ``AnnotatedGEM``
    and ``GeneSetCollection`` objects.
    """

    gem = param.ClassSelector(class_=AnnotatedGEM, doc="""\
    An ``AnnotatedGEM`` object.""", default=None, precedence=-1.0)

    gene_set_collection = param.ClassSelector(class_=GeneSetCollection, doc=dedent("""\
    A ``GeneSetCollection`` object."""), default=None, precedence=-1.0)

    selected_gene_sets = param.ListSelector(default=[None], doc=dedent("""\
    A list of keys from the provided GeneSetCollection (stored in gene_set_collection)
    that are to be used for selecting sets of genes from the count matrix."""))

    selected_genes = param.Parameter(default=None, doc=dedent("""\
    A list of genes to use in indexing from the count matrix. This parameter takes
    priority over all other gene selecting methods. That means that selected
    GeneSets (or combinations thereof) will have no effect."""), precedence=-1.0)

    gene_set_mode = param.ObjectSelector(
        default="union",
        objects=["complete", "union", "intersection"],
        doc=dedent("""\
    Controls how any selected gene sets are returned by the interface.
    **complete**
        Returns the entire gene set of the ``AnnotatedGEM``.
    **union**
        Returns the union of the selected gene sets support.
    **intersection**
        Returns the intersection of the selected gene sets support.
    """))

    sample_subset = param.Parameter(default=None, precedence=-1.0, doc=dedent("""\
    A list of samples to use in a given operation. These can be supplied
    directly as a list of genes, or can be drawn from a given GeneSet."""))

    count_variable = param.ObjectSelector(default=None, precedence=1.0,
                                          doc="The name of the count matrix used.",
                                          objects=[None], check_on_set=False)

    annotation_variables = param.List(doc=dedent("""\
    The name of the active annotation variable(s). These are the annotation columns that will
    be control the subset returned by ``y_annotation_data``."""), precedence=-1.0, default=[None])

    count_mask = param.ObjectSelector(doc=dedent("""\
    The type of mask to use for the count matrix.
    **complete**
        Returns the entire count matrix as numbers.
    **masked**
        Returns the entire count matrix with zero or missing as NaN values.
    **dropped**
        Returns the count matrix without genes that have zero or missing values.
    """), default='complete', objects=["complete", "masked", "dropped"], precedence=1.0)

    annotation_mask = param.ObjectSelector(doc=dedent("""\
    The type of mask to use for the target array.
    **complete**
        Returns the entire target array.
    **dropped**
        Returns the target array without samples that have zero or missing values.
    """), default='complete', objects=["complete", "dropped"], precedence=-1.0)

    count_transform = param.Callable(default=None, precedence=-1.0, doc=dedent("""\
    A transform that will be run on the `x_data` that is supplied by this Interface.
    The transform runs on the subset of the matrix that has been selected."""))

    @singledispatchmethod
    def _interface_dispatch(*args, **params):
        raise TypeError(f"Source of type: {type(args[0])} not supported.")

    def __init__(self, *args, **params):

        # If the user passes a string, place it as a single item within a list.
        if isinstance(params.get("annotation_variables"), str):
            params["annotation_variables"] = [params.get("annotation_variables")]

        if isinstance(params.get("selected_gene_sets"), str):
            params["selected_gene_sets"] = [params.get("selected_gene_sets")]

        if args:
            params = self._interface_dispatch(*args, **params)

        super().__init__(**params)

    @_interface_dispatch.register(AnnotatedGEM)
    @staticmethod
    def _parse_annotated_gem(annotated_gem: AnnotatedGEM, *_args, **params) -> dict:
        """
        Parse arguments for creation of a new `Interface` instance from an `AnnotatedGEM`.

        Parameters
        ----------
        annotated_gem : AnnotatedGEM
            A `GSForge.AnnotatedGEM` object.

        _args :
            Not used.

        params :
            Parameters to initialize this `Interface` with.

        Returns
        -------
        params : dict
            A parsed parameter dictionary.
        """
        params = {"gem": annotated_gem,
                  "count_variable": annotated_gem.count_array_name,
                  **params}
        return params

    @_interface_dispatch.register(GeneSetCollection)
    @staticmethod
    def _parse_gene_set_collection(gene_set_collection: GeneSetCollection, *_args, **params) -> dict:
        """
        Parse arguments for creation of a new `Interface` instance from an `GeneSetCollection`.

        Parameters
        ----------
        annotated_gem : AnnotatedGEM
            A `GSForge.AnnotatedGEM` object.

        _args :
            Not used.

        params :
            Parameters to initialize this `Interface` with.

        Returns
        -------
        params : dict
            A parsed parameter dictionary.
        """
        if gene_set_collection.gem is not None:
            params = {"gem": gene_set_collection.gem,
                      "count_variable": gene_set_collection.gem.count_array_name,
                      **params}
        params = {"gene_set_collection": gene_set_collection, **params}

        return params

    @property
    def active_count_variable(self) -> str:
        """Returns the name of the currently active count matrix."""
        if self.count_variable is not None:
            count_variable = self.count_variable
        else:
            count_variable = self.gem.count_array_name
        return count_variable

    @property
    def gene_index_name(self) -> str:
        """Returns the name of the gene index."""
        return self.gem.gene_index_name

    @property
    def sample_index_name(self) -> str:
        """Returns the name of the sample index."""
        return self.gem.sample_index_name

    def get_sample_index(self) -> np.ndarray:
        """
        Get the currently selected sample index as a numpy array.

        Returns
        -------
        np.ndarray
            An array of the currently selected samples.
        """
        logger.info(f'Determining sample index.')

        if self.sample_subset is not None:
            # We need to load the data if a user-list is supplied to prevent some straneg issues
            # with nested numpy arrays being returned.
            self.gem.data.load()
            subset = self.gem.data.sel({self.gem.sample_index_name: self.sample_subset})
            logger.info('No sample subset selected.')
        else:
            subset = self.gem.data

        if self.annotation_mask == "complete":
            pass
        elif self.annotation_mask == "dropped":
            if self.annotation_variables != [None]:
                subset = subset[self.annotation_variables].dropna(dim=self.gem.sample_index_name)

            subset = subset.dropna(dim=self.gem.sample_index_name)
            logger.info(
                f'Dropped samples with missing labels, {subset[self.gem.sample_index_name].shape[0]} samples remain.')

        selected_samples = subset[self.gem.sample_index_name].copy(deep=True).values
        return selected_samples

    @property
    def get_selection_indices(self) -> dict:
        """Returns the currently selected indexes as a dictionary."""
        return {self.gem.gene_index_name: self.get_gene_index(),
                self.gem.sample_index_name: self.get_sample_index()}

    @property
    def x_count_data(self) -> Union[xr.DataArray, None]:
        """
        Returns the currently selected 'x_data'. Usually this will be a subset of the active count array.

        Note: In constructing the a gene index, the count data is constructed first in order to infer
        coordinate selection based on masking.

        Returns
        -------
        xarray.Dataset
            The selection of the currently active count data.
        """

        gene_set_combinations = {
            "union": lambda sel_gs: self.gene_set_collection.union(sel_gs),
            "intersection": lambda sel_gs: self.gene_set_collection.intersection(sel_gs),
            "joint_difference": lambda sel_gs: self.gene_set_collection.joint_difference(sel_gs),
            "complete": lambda sel_gs: self.gem.gene_index,
        }

        mask_modes = {
            "complete": lambda c: c.fillna(0),
            "masked": lambda c: c.where(c > 0.0),
            "dropped": lambda c: c.where(c > 0.0).dropna(dim=self.gem.gene_index_name),
        }

        # Ensure the correct count array is selected.
        count_variable = self.count_variable if self.count_variable is not None else self.gem.count_array_name

        logger.info(f'Preparing count data from the variable {count_variable}.')
        if self.selected_gene_sets == ['all']:
            self.param.set_param(selected_gene_sets=list(self.gene_set_collection.gene_sets.keys()))

        if self.selected_genes is not None:
            support = self.selected_genes
            logger.info(f'Gene selection of: {support.shape[0]} genes provided.')

        elif self.selected_gene_sets == [None]:
            support = self.gem.gene_index

        # Otherwise, use  some combination of GeneSet supports should be used.
        else:
            support = gene_set_combinations[self.gene_set_mode](self.selected_gene_sets)
            logger.info(
                f'Selected {len(self.selected_gene_sets)} GeneSets, '
                f'using mode: {self.gene_set_mode} for a support of size: {support.shape[0]}.')

        # Now the counts can be selected based on the gene support, and the selected samples.
        sample_support = self.get_sample_index()

        # Check the overlap of the selected genes with those in the GEM index.
        # Give the user a warning if any unavailable genes are requested, but still
        # return the available genes.
        avail_support = np.intersect1d(self.gem.gene_index, support)
        if len(avail_support) < len(support):
            diff = len(support) - len(avail_support)
            warnings.warn(
                f'{diff} Unavailable genes of the {len(support)} requested. '
                f'Using the {len(avail_support)} available.',
                UserWarning)

        counts = self.gem.data.sel({self.gem.gene_index_name: avail_support,
                                    self.gem.sample_index_name: sample_support})[count_variable]
        logger.info(f'Selected count array of shape: {counts.shape}')

        logger.info(f'Preparing count data using mask mode {self.count_mask}.')
        counts = mask_modes[self.count_mask](counts)

        # Optional transform.
        if self.count_transform is not None:
            logger.info(f'Applying given transform to counts...')
            counts = self.count_transform(counts.copy(deep=True))

        return counts

    def get_gene_index(self) -> np.array:
        """
        Get the currently selected gene index as a numpy array.

        Returns
        -------
        np.ndarray
            An array of the currently selected genes.
        """
        logger.info(f'Preparing the gene index, this requires determining x_data.')
        return self.x_count_data[self.gene_index_name].values.copy()

    @property
    def y_annotation_data(self) -> Union[xr.Dataset, xr.DataArray, None]:
        """
        Returns the currently selected 'y_data', or None, based on the `selected_annotation_variables` parameter.

        Returns
        -------
            An ``xarray.Dataset`` of the currently selected y_data.
        """
        if (self.annotation_variables is None) or (self.annotation_variables == [None]):
            logger.info('No annotations selected.')
            return None

        logger.info(f'The following annotations where selected: {self.annotation_variables}.')
        sample_index = self.get_sample_index()
        # If only one label has been selected, return this as an xarray.DataArray.
        if len(self.annotation_variables) == 1:
            return self.gem.data[self.annotation_variables].sel(
                {self.gem.sample_index_name: sample_index})[self.annotation_variables[0]].copy(deep=True)

        return self.gem.data[self.annotation_variables].sel(
                {self.gem.sample_index_name: sample_index}).copy(deep=True)

    def get_gem_data(self, single_object=False, output_type='xarray', **params):
        """
        Returns count [and annotation] data based on the current parameters.

        Users should call gsf.get_gem_data
        """

        if params:
            logger.info(f'params {params}')
            self.param.set_param(**params)

        def xarray_single(self_):
            logger.info('Returning data as a single ``xarray.DataArray``.')
            if self_.y_annotation_data is not None:
                data = xr.merge([self_.x_count_data, self_.y_annotation_data])
            else:
                data = self_.x_count_data
            return data

        def xarray_tuple(self_):
            logger.info('Returning counts as an xarray.DataArray and annotations as an xarray.Dataset.')
            return self_.x_count_data, self_.y_annotation_data

        def pandas_single(self_):
            logger.info('Returning counts and annotations as a single pandas.DataFrame.')
            if self_.y_annotation_data is not None:
                data = xr.merge([self_.x_count_data, self_.y_annotation_data])
            else:
                data = self_.x_count_data
            return data.to_dataframe()

        def pandas_tuple(self_):
            # It may be faster to call .values and create a new dataframe instead of unstacking.
            logger.info('Returning a tuple of counts and annotations each as a pandas.DataFrame.')
            if self_.y_annotation_data is not None:
                return self_.x_count_data.to_dataframe().unstack().droplevel(0, axis=1), \
                       self_.y_annotation_data.to_dataframe()
            else:
                return self_.x_count_data.to_dataframe().unstack().droplevel(0, axis=1), self_.y_annotation_data

        def numpy_single(self_):
            return np.dstack(self_.x_count_data.values, self_.y_annotation_data.values)

        def numpy_tuple(self_):
            return self_.x_count_data.values, self_.y_annotation_data.values

        modes = {
            ('xarray', True): xarray_single,
            ('xarray', False): xarray_tuple,
            ('pandas', True): pandas_single,
            ('pandas', False): pandas_tuple,
            ('numpy', True): numpy_single,
            ('numpy', False): numpy_tuple,
        }
        key = (output_type, single_object)
        # TODO: Clarify error message.
        if key not in modes.keys():
            raise ValueError(f'key given: {key} is not one of the available'
                             f'types: {list(modes.keys())}')

        return modes[key](self)

class CallableInterface(Interface, param.ParameterizedFunction):

    @transient_log_handler
    def __new__(cls, *args, **params):
        logger.debug('Creating a new GSForge.CallableInterface instance.')
        if args:
            params = Interface._interface_dispatch(*args, **params)
            # TODO: Convert to a helper function or consider the implementation in operations.core (commented out).
            if isinstance(params.get("annotation_variables"), str):
                params["annotation_variables"] = [params.get("annotation_variables")]

            if isinstance(params.get("selected_gene_sets"), str):
                params["selected_gene_sets"] = [params.get("selected_gene_sets")]
        inst = cls.instance(**params)  # See the param code for more on this `instance` function.
        return inst.__call__()

    def __call__(self):
        raise NotImplementedError

max_stars_repo_name: SystemsGenetics/GEMprospector

max_issues_repo_name: SystemsGenetics/GEMprospector

max_forks_repo_name: SystemsGenetics/GEMprospector


text:
# This script analyzes the effect of fractional utilities and
# mechanisms for selecting a capacity that maximizes a non
# welfare metric. Since as we will show below, the social
# welfare maximizing metric is sub-optimal for decentralization

# Consider agents who each have a capacity. For each agent
# with capacity at least as large as the system capacity, this
# agent derives $capacity * fraction_of_hashpower$ utility
# and zero otherwise. Since the sum of fractions of hashpower
# is 1, the social welfare maximizing capacity is to just select
# the largest possible capacity for a the given set of agents.
# Thus, the system tolerates high throughput, but the lowest
# decentralization.

import math
import numpy as np
import pprint
import agent
import constants
import copy

num_agents = 10
initial_capacity = 0
behavior_type = constants.HONEST_AGENT
report_type = constants.DIRECT_CAPACITY_REPORT

def increasing_node_capacities(n, start):
    agents = []
    caps = np.arange(start, start + n)

    for i in range(n):
        agents.append(agent.Agent(behavior_type,
                                  caps[i],
                                  math.pow(caps[i], 2),
                                  constants.DIRECT_CAPACITY_REPORT))
    return agents

def random_node_capacities(n, low, high):
    agents = []
    caps = np.random.randint(low, (high + 1) * 10, n)
    hashes = np.random.randint(low, (high + 1), n)

    for i in range(n):
        agents.append(agent.Agent(behavior_type,
                                  caps[i],
                                  hashes[i],
                                  constants.DIRECT_CAPACITY_REPORT))
    return agents

def maximize_obj(reports):
    max_cap = {
        "normal": 0,
        "sqrt_capacity": 0,
        "squared_capacity": 0,
        "cubed_capacity": 0,
        "sqrt_hashpower": 0,
        "squared_hashpower": 0,
        "cubed_hashpower": 0,
        "sqrt_both": 0,
        "squared_hashpower": 0,
        "cubed_hashpower": 0,
    }
    max_obj = {
        "normal": 0,
        "sqrt_capacity": 0,
        "squared_capacity": 0,
        "cubed_capacity": 0,
        "sqrt_hashpower": 0,
        "squared_hashpower": 0,
        "cubed_hashpower": 0,
        "sqrt_both": 0,
        "squared_both": 0,
        "cubed_both": 0,
    }
    hashpower_sum = {
        "normal": 0,
        "sqrt_capacity": 0,
        "squared_capacity": 0,
        "cubed_capacity": 0,
        "sqrt_hashpower": 0,
        "squared_hashpower": 0,
        "cubed_hashpower": 0,
        "sqrt_both": 0,
        "squared_both": 0,
        "cubed_both": 0,
    }

    for r in reports:
        (capacity, _) = r
        summed_hashpower = sum([
            reports[j][1] for j in range(len(reports)) if reports[j][0] >= capacity
        ])

        if summed_hashpower * capacity > max_obj["normal"]:
            max_cap["normal"] = capacity
            max_obj["normal"] = summed_hashpower * capacity
            hashpower_sum["normal"] = summed_hashpower
        if summed_hashpower * math.sqrt(capacity) > max_obj["sqrt_capacity"]:
            max_cap["sqrt_capacity"] = capacity
            max_obj["sqrt_capacity"] = summed_hashpower * math.sqrt(capacity)
            hashpower_sum["sqrt_capacity"] = summed_hashpower
        if summed_hashpower * math.pow(capacity, 2) > max_obj["squared_capacity"]:
            max_cap["squared_capacity"] = capacity
            max_obj["squared_capacity"] = summed_hashpower * math.pow(capacity, 2)
            hashpower_sum["squared_capacity"] = summed_hashpower
        if summed_hashpower * math.pow(capacity, 3) > max_obj["cubed_capacity"]:
            max_cap["cubed_capacity"] = capacity
            max_obj["cubed_capacity"] = summed_hashpower * math.pow(capacity, 3)
            hashpower_sum["cubed_capacity"] = summed_hashpower
        if math.sqrt(summed_hashpower) * capacity > max_obj["sqrt_hashpower"]:
            max_cap["sqrt_hashpower"] = capacity
            max_obj["sqrt_hashpower"] = math.sqrt(summed_hashpower) * capacity
            hashpower_sum["sqrt_hashpower"] = summed_hashpower
        if math.pow(summed_hashpower, 2) * capacity > max_obj["squared_hashpower"]:
            max_cap["squared_hashpower"] = capacity
            max_obj["squared_hashpower"] = math.pow(summed_hashpower, 2) * capacity
            hashpower_sum["squared_hashpower"] = summed_hashpower
        if math.pow(summed_hashpower, 3) * capacity > max_obj["cubed_hashpower"]:
            max_cap["cubed_hashpower"] = capacity
            max_obj["cubed_hashpower"] = math.pow(summed_hashpower, 3) * capacity
            hashpower_sum["cubed_hashpower"] = summed_hashpower
        if math.sqrt(summed_hashpower) * math.sqrt(capacity) > max_obj["sqrt_both"]:
            max_cap["sqrt_both"] = capacity
            max_obj["sqrt_both"] = math.sqrt(summed_hashpower) * math.sqrt(capacity)
            hashpower_sum["sqrt_both"] = summed_hashpower
        if math.pow(summed_hashpower, 2) * math.pow(capacity, 2) > max_obj["squared_both"]:
            max_cap["squared_both"] = capacity
            max_obj["squared_both"] = math.pow(summed_hashpower, 2) * math.pow(capacity, 2)
            hashpower_sum["squared_both"] = summed_hashpower
        if math.pow(summed_hashpower, 3) * math.pow(capacity, 3) > max_obj["cubed_both"]:
            max_cap["cubed_both"] = capacity
            max_obj["cubed_both"] = math.pow(summed_hashpower, 3) * math.pow(capacity, 3)
            hashpower_sum["cubed_both"] = summed_hashpower
    return (max_cap, max_obj, hashpower_sum)

def get_agent_reports(agents):
    return list(map(lambda a: (a.capacity, a.hash_power), agents))

def get_true_hp_sum(reports, capacity):
    return sum([
        reports[j][1] for j in range(len(reports)) if reports[j][0] >= capacity
    ])

# run simulation a bunch of times
for t in range(100):
    # sample random node capacities
    agents = random_node_capacities(10, 10, 100)
    max_capacity_of_all_agents = max([ a.capacity for a in agents ])
    # agents = increasing_node_capacities(100, 10)
    # get reports from agents
    reports = get_agent_reports(agents)
    # find maximizing capacity
    (max_cap, _, hash_sum) = maximize_obj(reports)
    # logging
    print("Run {}, reports = {}".format(t, reports))
    print("Max cap\n")
    pprint.pprint(max_cap)
    # test misreporting above and below for all agents below
    for inx, a in enumerate(agents):
        # agents whose capacity is above or below may be
        # able to affect outcome by misreporting below
        for key in max_cap:
            if a.capacity < max_cap[key]:
                utility_of_agent = 0
            else:
                utility_of_agent = max_cap[key] * (a.hash_power * 1.0 / hash_sum[key])

            if a.capacity < max_cap[key] or a.capacity > max_cap[key]:
                agents_copy = copy.deepcopy(agents)
                ctr = 1
                while a.capacity - ctr > 0:
                    new_capacity = a.capacity - ctr
                    agents_copy[inx].capacity = new_capacity
                    new_reports = get_agent_reports(agents_copy)
                    (new_max_cap, new_values, hp_sum) = maximize_obj(new_reports)
                    ctr += 1

                    if a.capacity < new_max_cap[key]:
                        utility_from_deviation = 0
                    else:
                        utility_from_deviation = new_max_cap[key] * (a.hash_power * 1.0 / get_true_hp_sum(reports, new_max_cap[key]))

                    if utility_from_deviation > utility_of_agent and new_max_cap[key] != max_cap[key]:
                        str = "KEY = {}, NEW_MAX_CAP = {}, TRUE_CAP = {}, DEVIATION = {}, OLD_CAP = {}, NEW_U = {}, OLD_U = {}, NEW_SUM = {}, OLD_SUM = {}"
                        print(str.format(key,
                                         new_max_cap[key],
                                         a.capacity,
                                         new_capacity,
                                         max_cap[key],
                                         utility_from_deviation,
                                         utility_of_agent,
                                         hp_sum[key],
                                         hash_sum[key],
                                        ))

                ctr = 1
                while a.capacity + ctr < max_capacity_of_all_agents:
                    new_capacity = a.capacity + ctr
                    agents_copy[inx].capacity = new_capacity
                    new_reports = get_agent_reports(agents_copy)
                    (new_max_cap, new_values, hp_sum) = maximize_obj(new_reports)
                    ctr += 1

                    if a.capacity < new_max_cap[key]:
                        utility_from_deviation = 0
                    else:
                        utility_from_deviation = new_max_cap[key] * (a.hash_power * 1.0 / get_true_hp_sum(reports, new_max_cap[key]))

                    if utility_from_deviation > utility_of_agent and new_max_cap[key] != max_cap[key]:
                        str = "KEY = {}, NEW_MAX_CAP = {}, TRUE_CAP = {}, DEVIATION = {}, OLD_CAP = {}, NEW_U = {}, OLD_U = {}"
                        print(str.format(key,
                                         new_max_cap[key],
                                         a.capacity,
                                         new_capacity,
                                         max_cap[key],
                                         utility_from_deviation,
                                         utility_of_agent
                                        ))
    print("")

max_stars_repo_name: drewstone/dynamic-governanceq

max_issues_repo_name: drewstone/dynamic-governanceq

max_forks_repo_name: drewstone/dynamic-governanceq


text:
# -*- coding: utf-8 -*-
"""
Created on Wed May 20 15:05:52 2015

@author: joseph
"""

import time
import math

import numpy as np

from memory_profiler import memory_usage

import warnings

import Dispatcher
import Trawler

from os.path import expanduser

class Pool(object):
    memory_test_interval_seconds = 10
    target_process_time_minutes = 60

    MAX_CHECK_STATUS_LOOPS = 1000000
    timeout_time_s = 1 * 60 * 60 #1 hours times 60 minutes per hour times 60 seconds per minute
    TARGET_PROCESS_TIMES_MINUTES = 30.0

    MINIMUM_NUMBER_TRAWLERS = 10

    def __init__(self, maximum_number_of_jobs_running = -1, cluster_dispatcher_type = "SLURM", TEST_MODE = False):
        "initialize a ClusterPool.Pool object to work on a specific cluster architecture"
        self.maximum_number_of_jobs_running = maximum_number_of_jobs_running
        self.cluster_dispatcher_type = cluster_dispatcher_type

        self.TEST_MODE = TEST_MODE

    def map(self, function_to_apply, list_of_objects_to_calculate, index_of_slowest_calculation = -1, time_multiplier_worst_case = 6.0, memory_multiplier_worst_case = 1.2):
        "Use the cluster pool to map a function called calculate() on every member of the supplied list"

        self.list_of_objects_to_calculate = list_of_objects_to_calculate

        self.slowest_calculation = self.list_of_objects_to_calculate[index_of_slowest_calculation]

        self.my_function_to_apply = function_to_apply

        print "testing time and memory requirements..."
        self.slowest_calculation_completed = None

        start_time = time.time()

        if self.TEST_MODE:
            self.slowest_calculation_completed =self.slowest_calculation.calculate_safely()
        else:
            self.slowest_calculation_completed =self.slowest_calculation.calculate()
        self.slowest_time_elapsed_s = time.time() - start_time
        process_completion_time_to_request_minutes = ( time_multiplier_worst_case * self.slowest_time_elapsed_s ) / 60.0
        process_memory_to_request_MB = memory_multiplier_worst_case * self.peak_memory_usage_megabytes()

        number_of_processes_per_trawler = int(math.ceil(Pool.TARGET_PROCESS_TIMES_MINUTES / process_completion_time_to_request_minutes))
        number_of_trawlers_to_create = int(math.ceil(float(len(list_of_objects_to_calculate) / number_of_processes_per_trawler)))
        number_of_trawlers_to_create = max([number_of_trawlers_to_create, Pool.MINIMUM_NUMBER_TRAWLERS])

        print "done!  child processes will request %s MB of RAM and %s minutes per core" % (str(process_memory_to_request_MB), str(process_completion_time_to_request_minutes)  )

        print "creating dispatcher"
        my_dispatcher = Dispatcher.return_appropriate_dispatcher_object(self.cluster_dispatcher_type)(process_completion_time_to_request_minutes, process_memory_to_request_MB)

        process_id_to_list_index_dict = {}
        unfinished_process_ids = []

        calculated_objects = []
        print("creating saved objects")
        for i, obj in enumerate(list_of_objects_to_calculate):
            calculated_objects.append(None)
            if i == index_of_slowest_calculation or i == (index_of_slowest_calculation + len(self.list_of_objects_to_calculate)):
                continue
            new_id = my_dispatcher.save_python_object_file_return_id(obj)
            process_id_to_list_index_dict[new_id] = i
            unfinished_process_ids.append(new_id)

        calculated_objects[index_of_slowest_calculation] = self.slowest_calculation_completed

        print("Spawning Trawlers")
        my_dispatcher.start_subprocess_trawlers(number_of_trawlers_to_create)

        myTrawler = Trawler.Trawler(my_dispatcher.file_directory(), number=-1, TEST_MODE = self.TEST_MODE, IGNORE_MUTEX=True)

        total_completed_calculations = 1
        loop_count = 0

        # time.sleep(5)
        print("checking for finished calculations...")
        while len(unfinished_process_ids) > 0:
            time.sleep(5)
            try:
                myTrawler.find_and_complete_calculation()

            except Trawler.noMoreFilesToCalculate:
                pass

            loop_count = loop_count + 1

            #check for completed processes
            for process_id in unfinished_process_ids:
                try:
                    calc_obj = my_dispatcher.return_calculated_object(process_id)
                    calculated_objects[process_id_to_list_index_dict[process_id]] = calc_obj
                    unfinished_process_ids.remove(process_id)
                    my_dispatcher.delete_all_temp_files(process_id)

                    total_completed_calculations += 1
                    print("%i finished!" % total_completed_calculations)
                    continue

                except:

                    continue
            #TODO: add trawler re-start for killed processes

        my_dispatcher.kill_dispatcher()
        return calculated_objects

    def peak_memory_usage_megabytes(self):
        """Memory usage of the current process in kilobytes."""
        status = None
        result = {'peak': 0, 'rss': 0}
        try:
            # This will only work on systems with a /proc file system
            # (like Linux).
            status = open('/proc/self/status')
            for line in status:
                parts = line.split()
                key = parts[0][2:-1].lower()
                if key in result:
                    result[key] = int(parts[1])
        except IOError:
            mem_usage_MB_history = memory_usage(-1, interval=self.slowest_time_elapsed_s/100.0, timeout=self.slowest_time_elapsed_s)
            return np.max(mem_usage_MB_history)
        finally:
            if status is not None:
                status.close()
        return float(result['peak'] ) / 1024.0

    def close(self):
        pass

max_stars_repo_name: jgoodknight/ClusterPool

max_issues_repo_name: jgoodknight/ClusterPool

max_forks_repo_name: jgoodknight/ClusterPool


text:
#  NLP written by GAMS Convert at 04/21/18 13:53:50
#
#  Equation counts
#      Total        E        G        L        N        X        C        B
#        104       56       15       33        0        0        0        0
#
#  Variable counts
#                   x        b        i      s1s      s2s       sc       si
#      Total     cont   binary  integer     sos1     sos2    scont     sint
#         61       61        0        0        0        0        0        0
#  FX      0        0        0        0        0        0        0        0
#
#  Nonzero counts
#      Total    const       NL      DLL
#        311       92      219        0
#
#  Reformulation has removed 1 variable and 1 equation

from pyomo.environ import *

model = m = ConcreteModel()

m.x1 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x2 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x3 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x4 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x5 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x6 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x7 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x8 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x9 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x10 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x11 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x12 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x13 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x14 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x15 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x16 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x17 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x18 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x19 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x20 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x21 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x22 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x23 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x24 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x25 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x26 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x27 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x28 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x29 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x30 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x31 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x32 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x33 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x34 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x35 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x36 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x37 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x38 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x39 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x40 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x41 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x42 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x43 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x44 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x45 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x46 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x47 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x48 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x49 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x50 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x51 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x52 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x53 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x54 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x55 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x56 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x57 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x58 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x59 = Var(within=Reals,bounds=(None,None),initialize=0)
m.x60 = Var(within=Reals,bounds=(None,None),initialize=0)

m.obj = Objective(expr=1100*m.x55**2 + 500*m.x55 + 850*m.x56**2 + 120*m.x56 + 1225*m.x57**2 + 100*m.x57
                        + 1085, sense=minimize)

m.c2 = Constraint(expr=8.53242320819113*m.x42*m.x48 - 8.53242320819113*m.x39*m.x51 + 8.53242320819113*m.x48*m.x42 -
                       8.53242320819113*m.x51*m.x39 + m.x1 == 0)

m.c3 = Constraint(expr=8.53242320819113*m.x39*m.x51 - 8.53242320819113*m.x42*m.x48 - 8.53242320819113*m.x48*m.x42 +
                       8.53242320819113*m.x51*m.x39 + m.x2 == 0)

m.c4 = Constraint(expr=0.808561236623068*m.x43*m.x44 - 1.61712247324614*m.x43**2 - 6.84898929845422*m.x43*m.x53 +
                       0.808561236623068*m.x44*m.x43 + 6.84898929845422*m.x44*m.x52 + 6.84898929845422*m.x52*m.x44 -
                       1.61712247324614*m.x52**2 + 0.808561236623068*m.x52*m.x53 - 6.84898929845422*m.x53*m.x43 +
                       0.808561236623068*m.x53*m.x52 + m.x3 == 0)

m.c5 = Constraint(expr=0.808561236623068*m.x43*m.x44 + 6.84898929845422*m.x43*m.x53 + 0.808561236623068*m.x44*m.x43 -
                       1.61712247324614*m.x44**2 - 6.84898929845422*m.x44*m.x52 - 6.84898929845422*m.x52*m.x44 +
                       0.808561236623068*m.x52*m.x53 + 6.84898929845422*m.x53*m.x43 + 0.808561236623068*m.x53*m.x52 -
                       1.61712247324614*m.x53**2 + m.x4 == 0)

m.c6 = Constraint(expr=0.641004569212057*m.x41*m.x42 - 1.28200913842411*m.x41**2 - 2.79412248118076*m.x41*m.x51 +
                       0.641004569212057*m.x42*m.x41 + 2.79412248118076*m.x42*m.x50 + 2.79412248118076*m.x50*m.x42 -
                       1.28200913842411*m.x50**2 + 0.641004569212057*m.x50*m.x51 - 2.79412248118076*m.x51*m.x41 +
                       0.641004569212057*m.x51*m.x50 + m.x5 == 0)

m.c7 = Constraint(expr=0.641004569212057*m.x41*m.x42 + 2.79412248118076*m.x41*m.x51 + 0.641004569212057*m.x42*m.x41 -
                       1.28200913842411*m.x42**2 - 2.79412248118076*m.x42*m.x50 - 2.79412248118076*m.x50*m.x42 +
                       0.641004569212057*m.x50*m.x51 + 2.79412248118076*m.x51*m.x41 + 0.641004569212057*m.x51*m.x50 -
                       1.28200913842411*m.x51**2 + m.x6 == 0)

m.c8 = Constraint(expr=0.577543740445048*m.x42*m.x43 - 1.1550874808901*m.x42**2 - 4.89213521318159*m.x42*m.x52 +
                       0.577543740445048*m.x43*m.x42 + 4.89213521318159*m.x43*m.x51 + 4.89213521318159*m.x51*m.x43 -
                       1.1550874808901*m.x51**2 + 0.577543740445048*m.x51*m.x52 - 4.89213521318159*m.x52*m.x42 +
                       0.577543740445048*m.x52*m.x51 + m.x7 == 0)

m.c9 = Constraint(expr=0.577543740445048*m.x42*m.x43 + 4.89213521318159*m.x42*m.x52 + 0.577543740445048*m.x43*m.x42 -
                       1.1550874808901*m.x43**2 - 4.89213521318159*m.x43*m.x51 - 4.89213521318159*m.x51*m.x43 +
                       0.577543740445048*m.x51*m.x52 + 4.89213521318159*m.x52*m.x42 + 0.577543740445048*m.x52*m.x51 -
                       1.1550874808901*m.x52**2 + m.x8 == 0)

m.c10 = Constraint(expr=8*m.x38*m.x53 - 8*m.x44*m.x47 - 8*m.x47*m.x44 + 8*m.x53*m.x38 + m.x9 == 0)

m.c11 = Constraint(expr=8*m.x44*m.x47 - 8*m.x38*m.x53 + 8*m.x47*m.x44 - 8*m.x53*m.x38 + m.x10 == 0)

m.c12 = Constraint(expr=0.971095624357363*m.x40*m.x41 - 1.94219124871473*m.x40**2 - 5.25534102593397*m.x40*m.x50 +
                        0.971095624357363*m.x41*m.x40 + 5.25534102593397*m.x41*m.x49 + 5.25534102593397*m.x49*m.x41 -
                        1.94219124871473*m.x49**2 + 0.971095624357363*m.x49*m.x50 - 5.25534102593397*m.x50*m.x40 +
                        0.971095624357363*m.x50*m.x49 + m.x11 == 0)

m.c13 = Constraint(expr=0.971095624357363*m.x40*m.x41 + 5.25534102593397*m.x40*m.x50 + 0.971095624357363*m.x41*m.x40 -
                        1.94219124871473*m.x41**2 - 5.25534102593397*m.x41*m.x49 - 5.25534102593397*m.x49*m.x41 +
                        0.971095624357363*m.x49*m.x50 + 5.25534102593397*m.x50*m.x40 + 0.971095624357363*m.x50*m.x49 -
                        1.94219124871473*m.x50**2 + m.x12 == 0)

m.c14 = Constraint(expr=8.68055555555556*m.x40*m.x46 - 8.68055555555556*m.x37*m.x49 + 8.68055555555556*m.x46*m.x40 -
                        8.68055555555556*m.x49*m.x37 + m.x13 == 0)

m.c15 = Constraint(expr=8.68055555555556*m.x37*m.x49 - 8.68055555555556*m.x40*m.x46 - 8.68055555555556*m.x46*m.x40 +
                        8.68055555555556*m.x49*m.x37 + m.x14 == 0)

m.c16 = Constraint(expr=0.68259385665529*m.x40*m.x45 + 5.80204778156997*m.x40*m.x54 + 0.68259385665529*m.x45*m.x40 -
                        1.36518771331058*m.x45**2 - 5.80204778156997*m.x45*m.x49 - 5.80204778156997*m.x49*m.x45 +
                        0.68259385665529*m.x49*m.x54 + 5.80204778156997*m.x54*m.x40 + 0.68259385665529*m.x54*m.x49 -
                        1.36518771331058*m.x54**2 + m.x15 == 0)

m.c17 = Constraint(expr=0.68259385665529*m.x40*m.x45 - 1.36518771331058*m.x40**2 - 5.80204778156997*m.x40*m.x54 +
                        0.68259385665529*m.x45*m.x40 + 5.80204778156997*m.x45*m.x49 + 5.80204778156997*m.x49*m.x45 -
                        1.36518771331058*m.x49**2 + 0.68259385665529*m.x49*m.x54 - 5.80204778156997*m.x54*m.x40 +
                        0.68259385665529*m.x54*m.x49 + m.x16 == 0)

m.c18 = Constraint(expr=0.593802189645574*m.x44*m.x45 - 1.18760437929115*m.x44**2 - 2.9875672666543*m.x44*m.x54 +
                        0.593802189645574*m.x45*m.x44 + 2.9875672666543*m.x45*m.x53 + 2.9875672666543*m.x53*m.x45 -
                        1.18760437929115*m.x53**2 + 0.593802189645574*m.x53*m.x54 - 2.9875672666543*m.x54*m.x44 +
                        0.593802189645574*m.x54*m.x53 + m.x17 == 0)

m.c19 = Constraint(expr=0.593802189645574*m.x44*m.x45 + 2.9875672666543*m.x44*m.x54 + 0.593802189645574*m.x45*m.x44 -
                        1.18760437929115*m.x45**2 - 2.9875672666543*m.x45*m.x53 - 2.9875672666543*m.x53*m.x45 +
                        0.593802189645574*m.x53*m.x54 + 2.9875672666543*m.x54*m.x44 + 0.593802189645574*m.x54*m.x53 -
                        1.18760437929115*m.x54**2 + m.x18 == 0)

m.c20 = Constraint(expr=8.53242320819113*m.x39*m.x42 - 17.0648464163823*m.x39**2 + 8.53242320819113*m.x42*m.x39 -
                        17.0648464163823*m.x48**2 + 8.53242320819113*m.x48*m.x51 + 8.53242320819113*m.x51*m.x48 + m.x19
                         == 0)

m.c21 = Constraint(expr=8.53242320819113*m.x39*m.x42 + 8.53242320819113*m.x42*m.x39 - 17.0648464163823*m.x42**2 +
                        8.53242320819113*m.x48*m.x51 + 8.53242320819113*m.x51*m.x48 - 17.0648464163823*m.x51**2 + m.x20
                         == 0)

m.c22 = Constraint(expr=6.84898929845422*m.x43*m.x44 - 13.6234785969084*m.x43**2 + 0.808561236623068*m.x43*m.x53 +
                        6.84898929845422*m.x44*m.x43 - 0.808561236623068*m.x44*m.x52 - 0.808561236623068*m.x52*m.x44 -
                        13.6234785969084*m.x52**2 + 6.84898929845422*m.x52*m.x53 + 0.808561236623068*m.x53*m.x43 +
                        6.84898929845422*m.x53*m.x52 + m.x21 == 0)

m.c23 = Constraint(expr=6.84898929845422*m.x43*m.x44 - 0.808561236623068*m.x43*m.x53 + 6.84898929845422*m.x44*m.x43 -
                        13.6234785969084*m.x44**2 + 0.808561236623068*m.x44*m.x52 + 0.808561236623068*m.x52*m.x44 +
                        6.84898929845422*m.x52*m.x53 - 0.808561236623068*m.x53*m.x43 + 6.84898929845422*m.x53*m.x52 -
                        13.6234785969084*m.x53**2 + m.x22 == 0)

m.c24 = Constraint(expr=2.79412248118076*m.x41*m.x42 - 5.40924496236153*m.x41**2 + 0.641004569212057*m.x41*m.x51 +
                        2.79412248118076*m.x42*m.x41 - 0.641004569212057*m.x42*m.x50 - 0.641004569212057*m.x50*m.x42 -
                        5.40924496236153*m.x50**2 + 2.79412248118076*m.x50*m.x51 + 0.641004569212057*m.x51*m.x41 +
                        2.79412248118076*m.x51*m.x50 + m.x23 == 0)

m.c25 = Constraint(expr=2.79412248118076*m.x41*m.x42 - 0.641004569212057*m.x41*m.x51 + 2.79412248118076*m.x42*m.x41 -
                        5.40924496236153*m.x42**2 + 0.641004569212057*m.x42*m.x50 + 0.641004569212057*m.x50*m.x42 +
                        2.79412248118076*m.x50*m.x51 - 0.641004569212057*m.x51*m.x41 + 2.79412248118076*m.x51*m.x50 -
                        5.40924496236153*m.x51**2 + m.x24 == 0)

m.c26 = Constraint(expr=4.89213521318159*m.x42*m.x43 - 9.67977042636317*m.x42**2 + 0.577543740445048*m.x42*m.x52 +
                        4.89213521318159*m.x43*m.x42 - 0.577543740445048*m.x43*m.x51 - 0.577543740445048*m.x51*m.x43 -
                        9.67977042636317*m.x51**2 + 4.89213521318159*m.x51*m.x52 + 0.577543740445048*m.x52*m.x42 +
                        4.89213521318159*m.x52*m.x51 + m.x25 == 0)

m.c27 = Constraint(expr=4.89213521318159*m.x42*m.x43 - 0.577543740445048*m.x42*m.x52 + 4.89213521318159*m.x43*m.x42 -
                        9.67977042636317*m.x43**2 + 0.577543740445048*m.x43*m.x51 + 0.577543740445048*m.x51*m.x43 +
                        4.89213521318159*m.x51*m.x52 - 0.577543740445048*m.x52*m.x42 + 4.89213521318159*m.x52*m.x51 -
                        9.67977042636317*m.x52**2 + m.x26 == 0)

m.c28 = Constraint(expr=8*m.x38*m.x44 + 8*m.x44*m.x38 - 16*m.x44**2 + 8*m.x47*m.x53 + 8*m.x53*m.x47 - 16*m.x53**2
                         + m.x27 == 0)

m.c29 = Constraint(expr=8*m.x38*m.x44 - 16*m.x38**2 + 8*m.x44*m.x38 - 16*m.x47**2 + 8*m.x47*m.x53 + 8*m.x53*m.x47
                         + m.x28 == 0)

m.c30 = Constraint(expr=5.25534102593397*m.x40*m.x41 - 10.4316820518679*m.x40**2 + 0.971095624357363*m.x40*m.x50 +
                        5.25534102593397*m.x41*m.x40 - 0.971095624357363*m.x41*m.x49 - 0.971095624357363*m.x49*m.x41 -
                        10.4316820518679*m.x49**2 + 5.25534102593397*m.x49*m.x50 + 0.971095624357363*m.x50*m.x40 +
                        5.25534102593397*m.x50*m.x49 + m.x29 == 0)

m.c31 = Constraint(expr=5.25534102593397*m.x40*m.x41 - 0.971095624357363*m.x40*m.x50 + 5.25534102593397*m.x41*m.x40 -
                        10.4316820518679*m.x41**2 + 0.971095624357363*m.x41*m.x49 + 0.971095624357363*m.x49*m.x41 +
                        5.25534102593397*m.x49*m.x50 - 0.971095624357363*m.x50*m.x40 + 5.25534102593397*m.x50*m.x49 -
                        10.4316820518679*m.x50**2 + m.x30 == 0)

m.c32 = Constraint(expr=8.68055555555556*m.x37*m.x40 - 17.3611111111111*m.x37**2 + 8.68055555555556*m.x40*m.x37 -
                        17.3611111111111*m.x46**2 + 8.68055555555556*m.x46*m.x49 + 8.68055555555556*m.x49*m.x46 + m.x31
                         == 0)

m.c33 = Constraint(expr=8.68055555555556*m.x37*m.x40 + 8.68055555555556*m.x40*m.x37 - 17.3611111111111*m.x40**2 +
                        8.68055555555556*m.x46*m.x49 + 8.68055555555556*m.x49*m.x46 - 17.3611111111111*m.x49**2 + m.x32
                         == 0)

m.c34 = Constraint(expr=5.80204778156997*m.x40*m.x45 - 0.68259385665529*m.x40*m.x54 + 5.80204778156997*m.x45*m.x40 -
                        11.5160955631399*m.x45**2 + 0.68259385665529*m.x45*m.x49 + 0.68259385665529*m.x49*m.x45 +
                        5.80204778156997*m.x49*m.x54 - 0.68259385665529*m.x54*m.x40 + 5.80204778156997*m.x54*m.x49 -
                        11.5160955631399*m.x54**2 + m.x33 == 0)

m.c35 = Constraint(expr=5.80204778156997*m.x40*m.x45 - 11.5160955631399*m.x40**2 + 0.68259385665529*m.x40*m.x54 +
                        5.80204778156997*m.x45*m.x40 - 0.68259385665529*m.x45*m.x49 - 0.68259385665529*m.x49*m.x45 -
                        11.5160955631399*m.x49**2 + 5.80204778156997*m.x49*m.x54 + 0.68259385665529*m.x54*m.x40 +
                        5.80204778156997*m.x54*m.x49 + m.x34 == 0)

m.c36 = Constraint(expr=2.9875672666543*m.x44*m.x45 - 5.82213453330859*m.x44**2 + 0.593802189645574*m.x44*m.x54 +
                        2.9875672666543*m.x45*m.x44 - 0.593802189645574*m.x45*m.x53 - 0.593802189645574*m.x53*m.x45 -
                        5.82213453330859*m.x53**2 + 2.9875672666543*m.x53*m.x54 + 0.593802189645574*m.x54*m.x44 +
                        2.9875672666543*m.x54*m.x53 + m.x35 == 0)

m.c37 = Constraint(expr=2.9875672666543*m.x44*m.x45 - 0.593802189645574*m.x44*m.x54 + 2.9875672666543*m.x45*m.x44 -
                        5.82213453330859*m.x45**2 + 0.593802189645574*m.x45*m.x53 + 0.593802189645574*m.x53*m.x45 +
                        2.9875672666543*m.x53*m.x54 - 0.593802189645574*m.x54*m.x44 + 2.9875672666543*m.x54*m.x53 -
                        5.82213453330859*m.x54**2 + m.x36 == 0)

m.c38 = Constraint(expr=m.x1**2 + m.x19**2 <= 9)

m.c39 = Constraint(expr=m.x2**2 + m.x20**2 <= 9)

m.c40 = Constraint(expr=m.x3**2 + m.x21**2 <= 6.25)

m.c41 = Constraint(expr=m.x4**2 + m.x22**2 <= 6.25)

m.c42 = Constraint(expr=m.x5**2 + m.x23**2 <= 2.25)

m.c43 = Constraint(expr=m.x6**2 + m.x24**2 <= 2.25)

m.c44 = Constraint(expr=m.x7**2 + m.x25**2 <= 2.25)

m.c45 = Constraint(expr=m.x8**2 + m.x26**2 <= 2.25)

m.c46 = Constraint(expr=m.x9**2 + m.x27**2 <= 6.25)

m.c47 = Constraint(expr=m.x10**2 + m.x28**2 <= 6.25)

m.c48 = Constraint(expr=m.x11**2 + m.x29**2 <= 6.25)

m.c49 = Constraint(expr=m.x12**2 + m.x30**2 <= 6.25)

m.c50 = Constraint(expr=m.x13**2 + m.x31**2 <= 6.25)

m.c51 = Constraint(expr=m.x14**2 + m.x32**2 <= 6.25)

m.c52 = Constraint(expr=m.x15**2 + m.x33**2 <= 6.25)

m.c53 = Constraint(expr=m.x16**2 + m.x34**2 <= 6.25)

m.c54 = Constraint(expr=m.x17**2 + m.x35**2 <= 6.25)

m.c55 = Constraint(expr=m.x18**2 + m.x36**2 <= 6.25)

m.c56 = Constraint(expr=m.x37**2 + m.x46**2 <= 1.21)

m.c57 = Constraint(expr=m.x38**2 + m.x47**2 <= 1.21)

m.c58 = Constraint(expr=m.x39**2 + m.x48**2 <= 1.21)

m.c59 = Constraint(expr=m.x40**2 + m.x49**2 <= 1.21)

m.c60 = Constraint(expr=m.x41**2 + m.x50**2 <= 1.21)

m.c61 = Constraint(expr=m.x42**2 + m.x51**2 <= 1.21)

m.c62 = Constraint(expr=m.x43**2 + m.x52**2 <= 1.21)

m.c63 = Constraint(expr=m.x44**2 + m.x53**2 <= 1.21)

m.c64 = Constraint(expr=m.x45**2 + m.x54**2 <= 1.21)

m.c65 = Constraint(expr=m.x37**2 + m.x46**2 >= 0.81)

m.c66 = Constraint(expr=m.x38**2 + m.x47**2 >= 0.81)

m.c67 = Constraint(expr=m.x39**2 + m.x48**2 >= 0.81)

m.c68 = Constraint(expr=m.x40**2 + m.x49**2 >= 0.81)

m.c69 = Constraint(expr=m.x41**2 + m.x50**2 >= 0.81)

m.c70 = Constraint(expr=m.x42**2 + m.x51**2 >= 0.81)

m.c71 = Constraint(expr=m.x43**2 + m.x52**2 >= 0.81)

m.c72 = Constraint(expr=m.x44**2 + m.x53**2 >= 0.81)

m.c73 = Constraint(expr=m.x45**2 + m.x54**2 >= 0.81)

m.c74 = Constraint(expr=   m.x55 <= 2.5)

m.c75 = Constraint(expr=   m.x56 <= 3)

m.c76 = Constraint(expr=   m.x57 <= 2.7)

m.c77 = Constraint(expr=   m.x55 >= 0.1)

m.c78 = Constraint(expr=   m.x56 >= 0.1)

m.c79 = Constraint(expr=   m.x57 >= 0.1)

m.c80 = Constraint(expr=   m.x58 <= 3)

m.c81 = Constraint(expr=   m.x59 <= 3)

m.c82 = Constraint(expr=   m.x60 <= 3)

m.c83 = Constraint(expr=   m.x58 >= -3)

m.c84 = Constraint(expr=   m.x59 >= -3)

m.c85 = Constraint(expr=   m.x60 >= -3)

m.c86 = Constraint(expr=   m.x46 == 0)

m.c87 = Constraint(expr=   m.x13 - m.x55 == 0)

m.c88 = Constraint(expr=   m.x10 - m.x56 == 0)

m.c89 = Constraint(expr=   m.x1 - m.x57 == 0)

m.c90 = Constraint(expr=   m.x31 - m.x58 == 0)

m.c91 = Constraint(expr=   m.x28 - m.x59 == 0)

m.c92 = Constraint(expr=   m.x19 - m.x60 == 0)

m.c93 = Constraint(expr=   m.x11 + m.x14 + m.x16 == 0)

m.c94 = Constraint(expr=   m.x5 + m.x12 == -0.9)

m.c95 = Constraint(expr=   m.x2 + m.x6 + m.x7 == 0)

m.c96 = Constraint(expr=   m.x3 + m.x8 == -1)

m.c97 = Constraint(expr=   m.x4 + m.x9 + m.x17 == 0)

m.c98 = Constraint(expr=   m.x15 + m.x18 == -1.25)

m.c99 = Constraint(expr=   m.x29 + m.x32 + m.x34 == 0)

m.c100 = Constraint(expr=   m.x23 + m.x30 == -0.3)

m.c101 = Constraint(expr=   m.x20 + m.x24 + m.x25 == 0)

m.c102 = Constraint(expr=   m.x21 + m.x26 == -0.35)

m.c103 = Constraint(expr=   m.x22 + m.x27 + m.x35 == 0)

m.c104 = Constraint(expr=   m.x33 + m.x36 == -0.5)

reason: from pyomo


text:
import numpy as np
from setuptools import setup, Extension
from Cython.Build import cythonize

hybrid = Extension(
    "hybrid_mandel",
    sources=["python_mandel.py"],
    include_dirs=[np.get_include()],
    define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")]
)

cython = Extension(
    "cython_mandel",
    sources=["cython_mandel.pyx"],
    include_dirs=[np.get_include()],
    define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")]
)

extensions = [hybrid, cython]

setup(
    ext_modules = cythonize(extensions, compiler_directives={"language_level": "3"}),
)

max_stars_repo_name: arifmudi/Applying-Math-with-Python

max_issues_repo_name: arifmudi/Applying-Math-with-Python

max_forks_repo_name: arifmudi/Applying-Math-with-Python