text:
import math
import numpy as np

from . import BaseAnnotationDefinition

class Cube(BaseAnnotationDefinition):
    """
        Cube annotation object
    """
    type = "cube"

    def __init__(self, label, front_tl, front_tr, front_br, front_bl,
                 back_tl, back_tr, back_br, back_bl, angle=None,
                 attributes=None, description=None):
        super().__init__(description=description, attributes=attributes)
        self.front_bl = front_bl
        self.front_br = front_br
        self.front_tr = front_tr
        self.front_tl = front_tl
        self.back_bl = back_bl
        self.back_br = back_br
        self.back_tr = back_tr
        self.back_tl = back_tl

        self._angle = angle
        self.label = label

        self.keys = ["front_tl", "front_tr", "front_br", "front_bl",
                     "back_tl", "back_tr", "back_br", "back_bl"]

    @staticmethod
    def calculate_angle(b, c):
        a = [b[0] + 200, b[1]]
        if b in (a, c):
            return 0

        ang = math.degrees(
            math.atan2(c[1] - b[1], c[0] - b[0]) - math.atan2(a[1] - b[1], a[0] - b[0]))
        return ang + 360 if ang < 0 else ang

    @property
    def angle(self):
        if self._angle is None:
            self._angle = Cube.calculate_angle(self.front_tl, self.front_tr)
        return self._angle

    @property
    def x(self):
        return self.geo[:, 0]

    @property
    def y(self):
        return self.geo[:, 1]

    @property
    def geo(self):
        return np.asarray([self.front_tl,
                           self.front_tr,
                           self.front_br,
                           self.front_bl,
                           self.back_tl,
                           self.back_tr,
                           self.back_br,
                           self.back_bl
                           ])

    def show(self, image, thickness, with_text, height, width, annotation_format, color):
        """
        Show annotation as ndarray
        :param image: empty or image to draw on
        :param thickness:
        :param with_text: not required
        :param height: item height
        :param width: item width
        :param annotation_format: options: list(dl.ViewAnnotationOptions)
        :param color: color
        :return: ndarray
        """
        try:
            import cv2
        except (ImportError, ModuleNotFoundError):
            self.logger.error(
                'Import Error! Cant import cv2. Annotations operations will be limited. import manually and fix errors')
            raise

        image = cv2.polylines(image,
                              pts=[np.asarray(
                                  [self.front_bl, self.front_br, self.front_tr, self.front_tl]).round().astype(int)],
                              isClosed=True,
                              color=color,
                              thickness=thickness)
        image = cv2.polylines(image,
                              pts=[np.asarray([self.back_bl, self.back_br, self.back_tr, self.back_tl]).round().astype(
                                  int)],
                              isClosed=True,
                              color=color,
                              thickness=thickness)
        image = cv2.line(image,
                         pt1=tuple(np.asarray(self.front_bl).round().astype(int)),
                         pt2=tuple(np.asarray(self.back_bl).round().astype(int)),
                         color=color,
                         thickness=thickness)
        image = cv2.line(image,
                         pt1=tuple(np.asarray(self.front_br).round().astype(int)),
                         pt2=tuple(np.asarray(self.back_br).round().astype(int)),
                         color=color,
                         thickness=thickness)
        image = cv2.line(image,
                         pt1=tuple(np.asarray(self.front_tl).round().astype(int)),
                         pt2=tuple(np.asarray(self.back_tl).round().astype(int)),
                         color=color,
                         thickness=thickness)
        image = cv2.line(image,
                         pt1=tuple(np.asarray(self.front_tr).round().astype(int)),
                         pt2=tuple(np.asarray(self.back_tr).round().astype(int)),
                         color=color,
                         thickness=thickness)
        return image

    def to_coordinates(self, color):
        coordinates = {self.keys[idx]: {"x": float(x), "y": float(y), "z": 0}
                       for idx, [x, y] in enumerate(self.geo)}
        coordinates['angle'] = self.angle
        return coordinates

    @staticmethod
    def from_coordinates(coordinates):
        geo = list()
        for key, pt in enumerate(coordinates):
            geo.append([pt["x"], pt["y"]])
        return np.asarray(geo)

    @classmethod
    def from_json(cls, _json):
        if "coordinates" in _json:
            key = "coordinates"
        elif "data" in _json:
            key = "data"
        else:
            raise ValueError('can not find "coordinates" or "data" in annotation. id: {}'.format(_json["id"]))

        return cls(
            front_bl=np.asarray([_json[key]["front_bl"]['x'], _json[key]["front_bl"]['y']]),
            front_br=np.asarray([_json[key]["front_br"]['x'], _json[key]["front_br"]['y']]),
            front_tl=np.asarray([_json[key]["front_tl"]['x'], _json[key]["front_tl"]['y']]),
            front_tr=np.asarray([_json[key]["front_tr"]['x'], _json[key]["front_tr"]['y']]),
            back_bl=np.asarray([_json[key]["back_bl"]['x'], _json[key]["back_bl"]['y']]),
            back_br=np.asarray([_json[key]["back_br"]['x'], _json[key]["back_br"]['y']]),
            back_tl=np.asarray([_json[key]["back_tl"]['x'], _json[key]["back_tl"]['y']]),
            back_tr=np.asarray([_json[key]["back_tr"]['x'], _json[key]["back_tr"]['y']]),
            label=_json["label"],
            angle=_json[key]["angle"],
            attributes=_json.get("attributes", None)
        )

    @staticmethod
    def rotate(center, point, angle):
        angle = math.radians(angle)
        cx, cy = center
        px, py = point

        qx = cx + math.cos(angle) * (px - cx) - math.sin(angle) * (py - cy)
        qy = cy + math.sin(angle) * (px - cx) + math.cos(angle) * (py - cy)
        return [qx, qy]

    @classmethod
    def from_boxes_and_angle(cls,
                             front_left, front_top, front_right, front_bottom,
                             back_left, back_top, back_right, back_bottom,
                             label, angle=0, attributes=None):
        """
        Create cuboid by given front and back boxes with angle
        the angle calculate fom the center of each box
        """
        if angle != 0:
            front_center = [front_left + (front_right - front_left) / 2, front_top + (front_bottom - front_top) / 2]
            back_center = [back_left + (back_right - back_left) / 2, back_top + (back_bottom - back_top) / 2]
            front_tl = Cube.rotate(center=front_center, point=[front_left, front_top], angle=angle)
            front_tr = Cube.rotate(center=front_center, point=[front_right, front_top], angle=angle)
            front_br = Cube.rotate(center=front_center, point=[front_right, front_bottom], angle=angle)
            front_bl = Cube.rotate(center=front_center, point=[front_left, front_bottom], angle=angle)
            back_tl = Cube.rotate(center=back_center, point=[back_left, back_top], angle=angle)
            back_tr = Cube.rotate(center=back_center, point=[back_right, back_top], angle=angle)
            back_br = Cube.rotate(center=back_center, point=[back_right, back_bottom], angle=angle)
            back_bl = Cube.rotate(center=back_center, point=[back_left, back_bottom], angle=angle)
        else:
            front_tl = [front_left, front_top]
            front_tr = [front_right, front_top]
            front_br = [front_right, front_bottom]
            front_bl = [front_left, front_bottom]
            back_tl = [back_left, back_top]
            back_tr = [back_right, back_top]
            back_br = [back_right, back_bottom]
            back_bl = [back_left, back_bottom]

        return cls(
            front_tl=front_tl, front_tr=front_tr, front_br=front_br, front_bl=front_bl,
            back_tl=back_tl, back_tr=back_tr, back_br=back_br, back_bl=back_bl,
            label=label,
            angle=angle,
            attributes=attributes
        )

max_stars_repo_name: dataloop-ai/dtlpy

max_issues_repo_name: dataloop-ai/dtlpy

max_forks_repo_name: dataloop-ai/dtlpy


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from absl import app
from absl import flags
from absl import logging

import math
from math import ceil
from math import floor
from math import exp
from math import log

import numpy as np
from numpy import dot
from numpy import arange
from numpy import zeros
from numpy import ones
from numpy import multiply
from numpy import diagflat as diag
from numpy.random import rand
from numpy.linalg import svd
from numpy.linalg import cholesky as chol

import pandas as pd
from pandas import concat

import scipy
from scipy.special import gamma
from scipy.special import gammaincc

def getallobj(p,n,R,opt_type):

  # calculate all objective values for different shifts
  allobj=zeros(n)
  l=len(p)
  if opt_type == 1:
    for k in range(n):
      o1=np.sum(multiply((p[1:l-k-1] - p[1+k:l-1]), np.log(p[1:l-k-1] / p[1+k:l-1])))
      q_plus= p[l-1]*np.power(R, arange(0,k+1,1))
      q_minus=p[0]*np.power(R, arange(k,-1,- 1))
      o2=sum(multiply((p[l-k-1:l] - q_plus), np.log(p[l-k-1:l] / q_plus))) + sum(multiply((p[1:2+k] - q_minus), np.log(p[1:2+k] / q_minus)))
      o3=(p[0] + p[l-1])*k*(1 - R ** (k+1)) / (1 - R)*(- log(R))

      allobj[k]=o1 + o2 + o3
  else:
    for k in arange(n):
      allobj[k]=sum(multiply((p[k:l] - p[:l-k]),np.log(p[:l-k] / p[:l-k])))

  return allobj

def getgrad(p,k,R,opt_type):
  # Calculate gradient for a shift of k
  l=len(p)
  if opt_type == 1:
    rat1=np.concatenate((p[2 + k:],p[l-1]*np.power(R, arange(1,k+1)))) / p[1:l-1]
    rat2=np.concatenate((p[0]*np.power(R,arange(k,0,-1)),p[0:l-k-2])) / p[1:l-1]
    gradmid=(- np.log(rat1) + 1 - rat1) + (- np.log(rat2) + 1 - rat2)

    Nfactor=(1+k)*(1 - R ** (1+k)) / (1 - R)*(- log(R))
    gradN=sum(multiply(np.power(R,arange(0,k+1)), np.log(rat1[-k-1:]) + 1 - 1.0 / rat1[-k-1:])) + Nfactor
    gradmN=sum(multiply(np.power(R,arange(k,-1,-1)),(np.log(rat2[:k+1]) + 1 - 1.0 / rat2[:k+1]))) + Nfactor

    grad=np.concatenate(([gradmN],gradmid,[gradN]))

  else:
    rat=p[-l+k:] / p[:l-k]
    rati=p[:l-k] / p[-l+k:]
    logr=np.log(rat)
    grad=np.concatenate((zeros((k+1,1)),[logr + 1 - rati], [- logr + 1 - rat],zeros((k+1,1))))

  return grad

def getH(p,n,evals,R,opt_type):

  # Calculate the Hessian matrix
  l=len(p)

  if opt_type == 1:

    d=zeros(l)
    for k in range(n):
      dN = sum(np.power(R, arange(0,k+1)))/p[l-1] + sum(p[-k-1:]) / p[l-1] ** 2
      dmN=sum(np.power(R,arange(k-1,-1,-1))) / p[0] + sum(p[1:1+k]) / p[0] ** 2
      dmid=2.0/p[1:l-1] + (np.concatenate((p[2+k:], np.power(R, arange(1,k+1))*p[l-1])) + np.concatenate((np.power(R, arange(k,0,-1))*p[0], p[:l-2-k] )))/ p[1:l-1]**2
      d=d + evals[k]*np.concatenate(([dmN],dmid,[dN]))

    H=np.diagflat(d)

    q=zeros((n,1))
    for k in range(n):
      j = 1 + k
      q[k]=sum(multiply(evals[k:n],(- np.power(R,arange(0,n-k))) / p[j] - 1 / p[0]))
    H[0,1:n+1]=q.reshape(-1)
    H[1:n+1,0]=q.reshape(-1)

    # deal with i=l
    q=zeros((n,1))
    for k in range(n):
      i = l-k-1
      q[k]=sum(multiply(evals[k:n], -np.power(R,arange(0,n-k))/p[i] - 1.0/p[l-1]))
    H[l-1,range(l-2,l-n-2,-1) ]=q.reshape(-1)
    H[range(l-2,l-n-2,-1),l-1]=q.reshape(-1)

    for k in range(1,l-2):
      kmax = min(n,l-k-1)
      q=multiply(evals[:kmax], (-1/p[k]-1.0/p[k+1:k+1+kmax]))
      H[k, k+1:k+1+kmax]=q.reshape(-1)
      H[k+1:k+1+kmax,k]=q.reshape(-1)
  else:
    d=zeros(l,1)
    for a in arange(1,n).reshape(-1):
      d=d + dot(evals(a),(concat([[1.0 / p(arange(1,l - a)) + p(arange(1 + a,l)) / p(arange(1,l - a)) ** 2],[zeros(a,1)]]) + concat([[zeros(a,1)],[1.0 / p(arange(1 + a,l)) + p(arange(1,l - a)) / p(arange(1 + a,l)) ** 2]])))
    H=diag(d)
    for i in arange(1,l - 1).reshape(-1):
      amax=min(n,l - i)
      q=multiply(evals(arange(1,amax)),(- 1 / p(i) - 1.0 / p(arange(i + 1,i + amax))))
      H[i,arange(i + 1,i + amax)]=q
      H[arange(i + 1,i + amax),i]=q

  return H

def sample(x,CDF,r,n):

  # produce one sample from the distribution

  l=len(x)
  i=find(rand(1) < CDF,1)
  s=x(i)
  if i == 1 or i == l:
    # sample from a geometric with parameter R
    j=floor(log(rand(1)) / log(r))
    if i == 1:
      s=s - j / n
    else:
      s=s + j / n

  s=s + (rand(1) - 0.5) / n
  return s

def main(n=200, xmax=8, C=0.25, c_exp=2, opt_type=1, r=0.9, tol=1e-08, verbose=2):

  N = math.ceil(n*xmax)
  xmax = N / n

  x=arange(-N,N+1) / n
  l=len(x)

  # The two constraints (probability normalization, and the variance constraint) are handled by A*p=b as follows.
  # Note that even though the variance constraint is an inequality, it will always be tight at optimal, so we treat it as an equality constraint.

  if opt_type == 1:
    c=zeros(l)
    for i in range(l):
      if x[i]>0:
        c[i] = dot(((x[i] + 1 / 2 / n) ** (c_exp + 1) - (x[i] - 1 / 2 / n) ** (c_exp + 1)),(n / (c_exp + 1)))
      elif x[i]==0:
        c[i] = (1 / 2 / n) ** c_exp / (1 + c_exp)
      else:
        c[i]=dot(((abs(x[i]) + 1 / 2 / n) ** (c_exp + 1) - (abs(x[i]) - 1 / 2 / n) ** (c_exp + 1)),(n / (c_exp + 1)))

    # slight over-estimate, but getting the exact value requires a nasty infinite sum. This upper bounds the sum by an integral.
    # This still gives a true bound, without much loss
    cN = n*r**(-N-1/2)*(-n*np.log(r))**(-c_exp-1)*gamma(c_exp+1)*gammaincc(c_exp+1,-np.log(r)*(N-1/2))
    c[0]=cN
    c[l-1]=cN
    d = np.ones(l)
    d[0] = d[l-1] = 1/(1-r)
    A = np.stack((c,d),axis = 0)
    # A = [c,d]
  else:
    # If we're doing the under-estimating model, let c be the minimal value of the cost function within each quantization bin.
    c=max(0,abs(x)-1/(2*n)) ** c_exp
    A=np.stack((c,ones(l)),axis = 0)

  b = np.stack(([C],[1]),axis = 0)
  # We start with a guess for the distribution p.
  # This is a distribution designed to have farily heavy tails (which seems to make the solver work better), and to roughly satisfy the cost constraint.
  q = C ** (- 1 - 2 / c_exp)
  p = 1/(1+q*abs(x)**(c_exp+2))
  # make sure that p is normalized (it actually doesn't need to satisfy the cost constraint exactly initially -- the algorithm with ensure that)
  p = p/dot(A[1],p)
  iter=0

  t=1

  isfeasible = False

  allobj = getallobj(p,n,r,opt_type)
  primobj=np.amax(allobj)
  fval=np.log(sum(np.exp(t*(allobj - primobj)))) / t + primobj

  while 1:

    iter = iter + 1
    eta = np.exp(t*(allobj - primobj))
    eta=eta / sum(eta)
    allgrads=zeros((l,n))

    for k in range(n):
      allgrads[:,k]=getgrad(p,k,r,opt_type)
    grad=np.matmul(allgrads,eta.reshape(200,1))

    rpri=np.matmul(A,p.reshape(len(p),1))-b

    # The next few lines are typically the majority of the run-time.
    # Compute Hessian matrix
    H=getH(p,n,eta,r,opt_type)
    # strictly positive definite. This doesn't much change overall performance
    fullH = H + np.matmul(np.matmul(allgrads,np.diagflat(t*eta)),allgrads.T) - (1 - 1e-05)*t*np.matmul(grad,grad.T)

    try:
      R=scipy.linalg.cholesky(fullH)

      gtilde=np.linalg.solve(R.T,grad)

      AR=np.matmul(A, np.linalg.inv(R))
      U,S,V = svd(AR,full_matrices=False)
      #U,S,V=svd(AR,'econ',nargout=3)
      s=S
      vtilde = - gtilde - np.matmul(V.T,multiply(1.0/s.reshape(2,1), np.matmul(U.T,rpri))) + np.matmul(V.T,(np.matmul(V,gtilde)))
      v=np.linalg.solve(R,vtilde)

    except:
      if verbose >=1:
          print('Chlesky failed')
      vw = np.linalg.solve(np.stack((np.concatenate((fullH, A.T)), np.concatenate((A, zeros(2)))),axis = 0), np.stack((-grad, -rpri),axis = 0))
      v = vw[:]

    nwt_dec=np.matmul(- grad.T,v) / 2

    # implementing a line search. We move to p+dst*v
    dst=1
    while 1:
      pnew = p + dst*v.reshape(-1)
      if min(pnew) > 0:

        allobj=getallobj(pnew,n,r,opt_type)
        primobj=max(allobj)
        newfval=np.log(sum(np.exp(t*(allobj - primobj)))) / t + primobj
        if not isfeasible:
          break
        if newfval < fval + 0.1*dst*(np.matmul(grad.T,v)):
          break
        if dst < 1e-08:
          break
      dst=dst / 2

    if  not isfeasible and dst == 1:
      isfeasible=True
      if verbose >= 1:
          print('Feasible!')
    p=pnew
    fval=newfval
    if verbose >= 1:
      print('iter=%d  dst=%1.1e  primobj=%1.4f  nwt_dec=%1.2e  t=%1.1e \n' %(iter,dst,primobj,nwt_dec,t));
    # if verbose >=2:
    #     drawnow(draw_fig)
    if isfeasible:
      # n/e/t is an estimate for the duality gap, if we are at exact optimality for the soft-max problem.
      # Thus nwt_dec+n/e/1 is an estimate for the total duality gap
      if nwt_dec + n / exp(1) / t < tol:
          break
      # if we take a full Newton step or get extremely close to optimal, then increase t.
      if dst == 1 or nwt_dec < tol / 2:
          # This is a fairly modest increase in t.
          # This seems to make the solver the most robust, if not necessarily the most efficient.
          t=t*1.25

  cdf = np.cumsum(multiply(p,A[1].T))

  filename=('test_data/cactus_cdf_v%.2f.csv' %C)
  np.savetxt(filename, cdf, delimiter=",")

  filename=('test_data/cactus_x_v%.2f.csv' %C)
  np.savetxt(filename, x, delimiter=",")

  filename=('test_data/cactus_pdf_v%.2f.csv' %C)
  np.savetxt(filename, p, delimiter=",")

#   return x, p

if __name__ == '__main__':
  app.run(main(n=200, xmax=8, C=0.25, c_exp=2, opt_type=1, r=0.9, tol=1e-05, verbose=2))

ext: py

max_stars_repo_name: FredericW/privacy

max_issues_repo_name: FredericW/privacy

max_forks_repo_name: FredericW/privacy

reason: import numpy,from numpy,import scipy,from scipy


text:
import pandas as pd
import numpy as np
from algorithms.Algorithm import Algorithm
from datetime import datetime, timedelta

class Gini(Algorithm):

    def __init__(self, risk_score, model_name):
        super().__init__( risk_score, model_name )

    def get_optimal_portfolio(self):
        selected_prices_value = self.prices_df[self.selected_assets].dropna()
        num_portfolios = 1000
        years = len(selected_prices_value) / 253
        starting_value = selected_prices_value.iloc[0, :]
        ending_value = selected_prices_value.iloc[len(selected_prices_value) - 1, :]
        total_period_return = ending_value / starting_value
        annual_returns = (total_period_return ** (1 / years)) - 1
        port_returns = []
        port_gini_coefficient = []
        sharpe_ratio = []
        stock_weights = []
        num_assets = len(self.selected_assets)
        np.random.seed(101)

        for single_portfolio in range(num_portfolios):
            weights = np.random.random(num_assets)
            weights /= np.sum(weights)
            returns = np.dot(weights, annual_returns)
            gini_coefficient = (selected_prices_value.mad()).mad()
            sharpe = returns / gini_coefficient
            sharpe_ratio.append(sharpe)
            port_returns.append(returns * 100)
            port_gini_coefficient.append(gini_coefficient * 100)
            stock_weights.append(weights)
        portfolio = {'Returns': port_returns,
                     'Volatility': port_gini_coefficient,
                     'Sharpe Ratio': sharpe_ratio}
        for i, symbol in enumerate(self.selected_assets):
            portfolio[symbol] = [Weight[i] for Weight in stock_weights]
        df = pd.DataFrame(portfolio)
        columns = ['Returns', 'Volatility', 'Sharpe Ratio'] + [stock for stock in self.selected_assets]
        df = df[columns]
        best_sharpe_portfolio = df.loc[df['Sharpe Ratio'] == df['Sharpe Ratio'].max()]
        sharpe_portfolio = pd.DataFrame(columns=['Ticker', 'Weight'])
        for i in range(len(self.selected_assets)):
            ticker = self.selected_assets[i]
            weight = best_sharpe_portfolio.loc[:, ticker].iloc[0]
            sharpe_portfolio = sharpe_portfolio.append({'Ticker': ticker, 'Weight': weight}, ignore_index=True)
        sharpe_portfolio = sharpe_portfolio.set_index('Ticker')
        return sharpe_portfolio

max_stars_repo_name: avielfedida/RoboAdvisor

max_issues_repo_name: avielfedida/RoboAdvisor

max_forks_repo_name: avielfedida/RoboAdvisor

reason: import numpy


text:
from firedrake import *
from firedrake.petsc import PETSc
from firedrake import COMM_WORLD
import numpy as np
from scipy.stats import linregress
import os
import porousdrake.DPP.convergence.exact_solution as sol

try:
    import matplotlib.pyplot as plt

    plt.rcParams["contour.corner_mask"] = False
    plt.close("all")
except:
    warning("Matplotlib not imported")

def compute_error(computed_sol, analytical_sol, var_name, norm_type="L2"):
    # Now we compute the various metrics. First we
    # simply compute the L2 error between the analytic
    # solutions and the computed ones.
    error = errornorm(analytical_sol, computed_sol, norm_type=norm_type)

    # We keep track of all metrics using a Python dictionary
    error_dictionary = {var_name: error}

    return error_dictionary

def convergence_hp(
    solver,
    min_degree=1,
    max_degree=4,
    numel_xy=[2, 4, 8, 16, 32, 64, 128, 256],
    norm_type="L2",
    quadrilateral=True,
    name="",
    **kwargs
):
    if name:
        name += "_"
    for degree in range(min_degree, max_degree):
        p1_errors = np.array([])
        p2_errors = np.array([])
        v1_errors = np.array([])
        v2_errors = np.array([])
        num_cells = np.array([])
        mesh_size = np.array([])
        for n in numel_xy:
            nel_x = nel_y = n
            mesh = UnitSquareMesh(nel_x, nel_y, quadrilateral=quadrilateral)
            num_cells = np.append(num_cells, mesh.num_cells())
            mesh_size = np.append(mesh_size, 1.0 / n)

            p1_sol, v1_sol, p2_sol, v2_sol, p_e_1, v_e_1, p_e_2, v_e_2 = solver(
                mesh=mesh, degree=degree, **kwargs
            )
            error_dictionary = {}
            error_dictionary.update(compute_error(p1_sol, p_e_1, "p1_error", norm_type=norm_type))
            error_dictionary.update(compute_error(p2_sol, p_e_2, "p2_error", norm_type=norm_type))
            p1_errors = np.append(p1_errors, error_dictionary["p1_error"])
            p2_errors = np.append(p2_errors, error_dictionary["p2_error"])
            error_dictionary.update(compute_error(v1_sol, v_e_1, "v1_error"), norm_type=norm_type)
            error_dictionary.update(compute_error(v2_sol, v_e_2, "v2_error"), norm_type=norm_type)
            v1_errors = np.append(v1_errors, error_dictionary["v1_error"])
            v2_errors = np.append(v2_errors, error_dictionary["v2_error"])
        p1_errors_log2 = np.log2(p1_errors)
        p2_errors_log2 = np.log2(p2_errors)
        v1_errors_log2 = np.log2(v1_errors)
        v2_errors_log2 = np.log2(v2_errors)
        num_cells_log2 = np.log2(num_cells)
        mesh_size_log2 = np.log2(mesh_size)
        p1_slope, intercept1, r_value1, p_value1, stderr1 = linregress(
            mesh_size_log2, p1_errors_log2
        )
        p2_slope, intercept2, r_value2, p_value2, stderr2 = linregress(
            mesh_size_log2, p2_errors_log2
        )
        PETSc.Sys.Print(
            "\n--------------------------------------\nDegree %d: p1 slope error %f"
            % (degree, np.abs(p1_slope)),
            "\nDegree %d: p2 slope error %f" % (degree, np.abs(p2_slope)),
        )
        v1_slope, intercept_v1, r_value_v1, p_value_v1, stderr_v1 = linregress(
            mesh_size_log2, v1_errors_log2
        )
        v2_slope, intercept_v2, r_value_v2, p_value_v2, stderr_v2 = linregress(
            mesh_size_log2, v2_errors_log2
        )
        PETSc.Sys.Print(
            "\n--------------------------------------\nDegree %d: v1 slope error %f"
            % (degree, np.abs(v1_slope)),
            "\nDegree %d: v2 slope error %f" % (degree, np.abs(v2_slope)),
        )
        # _plot_errors(mesh_size, p1_errors, p1_slope, degree, name='p1_errors')
        # _plot_errors(mesh_size, p2_errors, p2_slope, degree, name='p2_errors')
        # _plot_errors(mesh_size, v1_errors, v1_slope, degree, name='v1_errors')
        # _plot_errors(mesh_size, v2_errors, v2_slope, degree, name='v2_errors')
        os.makedirs("results_%s" % name, exist_ok=True)
        np.savetxt(
            ("results_%s/errors_degree%d.dat" % (name, degree)),
            np.transpose(
                [-mesh_size_log2, p1_errors_log2, p2_errors_log2, v1_errors_log2, v2_errors_log2]
            ),
        )

    return

def _plot_errors(mesh_size, errors, slope, degree, name="Error"):
    plt.figure()
    plt.loglog(mesh_size, errors, "-o", label=(r"k = %d; slope = %f" % (degree, np.abs(slope))))
    plt.legend(loc="best")
    plt.grid(True)
    plt.savefig("%s_deg%i.png" % (name, degree))
    return

max_stars_repo_name: volpatto/porousdrake

max_issues_repo_name: volpatto/porousdrake

max_forks_repo_name: volpatto/porousdrake

reason: import numpy,from scipy


text:
import numpy as np
from numba import jit,njit
from scipy.stats import norm
from itertools import combinations as comb

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

def interaction_matrix(N,K,shape="roll"):
    """Creates an interaction matrix for a given K

    Args:
        N (int): Number of bits
        K (int): Level of interactions
        shape (str): Shape of interactions. Takes values 'roll' (default), 'random', 'diag'.

    Returns:
        numpy.ndarray: an NxN numpy array with diagonal values equal to 1 and rowSums=colSums=K+1
    """

    output = None
    if K == 0:
        output = np.eye(N,dtype=int)
    elif shape=="diag":
        tmp = [np.diag(np.ones(N-abs(z)),z) for z in range((-K),(K+1))]
        tmp = np.array(tmp)
        output = np.sum(tmp,0)
    elif shape=="updiag":
        tmp = [np.diag(np.ones(N-abs(z)),z) for z in range(0,(K+1))]
        tmp = np.array(tmp)
        output = np.sum(tmp,0)
    elif shape=="downdiag":
        tmp = [np.diag(np.ones(N-abs(z)),z) for z in range((-K),1)]
        tmp = np.array(tmp)
        output = np.sum(tmp,0)
    elif shape=="sqdiag":
        tmp = np.eye(N,dtype=int)
        tmp = tmp.repeat(K+1,axis=0)
        tmp = tmp.repeat(K+1,axis=1)
        output = tmp[0:N,0:N]
    elif shape=="roll":
        tmp = [1]*(K+1) + [0]*(N-K-1)
        tmp = [np.roll(tmp,z) for z in range(N)]
        tmp = np.array(tmp)
        output = tmp.transpose()
    elif shape=="random":
        output = random_binary_matrix(N,K+1,1)
    elif shape=="chess":
        print(f"Uncrecognized interaction type '{type}'")
    # print(f"Interaction shape '{shape}' selected")
    return output

def binary_combinations(N,R):
    """Generates all binary vectors with sum R

    Args:
        N (int): Size of a binary vector
        R (int): Sum of elements of a vector

    Returns:
        numpy.ndarray: an (N R)xN numpy array, where rows are vectors of size N and sum of elements R
    """

    tmp = []
    idx = comb(range(N),R)
    for i in idx:
        A = [0]*N
        for j in i:
            A[j] = 1
        tmp.append(A)
    output = np.reshape(tmp,(-1,N))
    return(output)

def random_binary_matrix(N,R,diag=None):
    """Generates a random binary square matrix with a given row/col sum

    Args:
        N (int): Number of rows/cols
        R (int): Sum of rows/cols
        diag (int or None): Fixed value for diagonal. Takes values None (default), 0, 1.

    Returns:
        numpy.ndarray: an NxN numpy array
    """

    if N==R:
        tmp = np.ones((N,R),dtype=int)
        return tmp
    elif N<R:
        print("Incorrect binary matrix. Check parameters.")
        return
    # create a minimal 2d matrix of zeros for easier indexing
    tmp = np.zeros(2*N,dtype=int).reshape(2,N)
    cl = binary_combinations(N,R)
    for i in range(N):
        colsums = np.sum(tmp,0)
        # remove excess ones
        idx = np.empty(0,dtype=int)
        for j in np.where(colsums>=R)[0]:
            k = np.where(cl[:,j]>0)[0]
            idx = np.union1d(idx,k)
        cl = np.delete(cl,idx,0)
        # remove excess zeros
        inx = np.empty(0,dtype=int)
        for j in np.where(colsums+N-i == R)[0]:
            k = np.where(cl[:,j]==0)[0]
            inx = np.union1d(inx,k)
        cl = np.delete(cl,inx,0)
        # temporarily ignore diagonal 1s or 0s
        cli = cl.copy()
        if diag is not None:
            ivx = np.where(cl[:,i]==diag)[0]
            cli = cl[ivx,]
            clk = (cli + colsums)[:,i+1:]
            tp = N-i-1 if diag==0 else 0 # tuning parameter
            ikx = np.where(clk+tp==R)[0]
            cli = np.delete(cli,ikx,0)

        ncl = cl.shape[0]
        ncli = cli.shape[0]
        if ncli > 0:
            ind = np.random.choice(ncli)
            tmp = np.vstack((tmp,cli[ind,:]))
        elif ncli==0 and ncl>0:
            print('Error creating non-zero diagonals. Rerun the function')
            return 0
        else:
            print('Incorrect binary matrix. Check the dimensions.')
            return
    output = np.delete(tmp,[0,1],0) # remove first 2 empty rows created above
    return(output)

###############################################################################
def binx(x,size=4,out=None):
    """Converts values to binary and back

    Args:
        x: Input value (can be any type)
        size (int): Desired output vector size (adds leading zeros if the output size is less than desired, ignores otherwise)
        out (str or None): Specifies output type. Is ignored at the moment.

    Returns:
        list: A list of 0s and 1s if x is int.
        int: A decimal equivalent if x is str, numpy.ndarray, list.
    """

    tmp = x
    if type(tmp) is int:
        tmp = np.binary_repr(tmp,size)
        tmp = [int(z) for z in tmp]
    elif type(tmp) is str:
        tmp = int(tmp,2)
    elif type(tmp) is np.ndarray or type(tmp) is list:
        tmp = np.sum(np.flip(tmp) * 2 ** (np.arange(len(tmp))))
    else:
        print("incorrect input for function binx")
    return tmp

def contrib_define(p,n,k,c,s,rho):
    """Defines a contribution matrix for given parameters

    Args:
        p (int): Number of landscapes (population size)
        n (int): Number of tasks per landscape
        k (int): Number of internal bits interacting with each bit
        c (int): Number of external bits interacting with each bit
        s (int): Number of landscapes considered for external bits
        rho (float): Correlation coefficient between landscapes

    Returns:
        numpy.ndarray: An (N*P)x2**(1+K+C*S) matrix of P contribution matrices with correlation RHO
    """

    #output = np.random.uniform(0,1,(2**(1+k+c*s),n*p))
    corrmat = np.repeat(rho,p*p).reshape(p,p) + (1-rho) * np.eye(p)
    corrmat = 2*np.sin((np.pi / 6 ) * corrmat)
    tmp = np.random.multivariate_normal(mean=[0]*p,cov=corrmat,size=(n*2**(1+k+c*s)))
    tmp = norm.cdf(tmp)
    tmp = np.reshape(tmp.T,(n*p,(2**(1+k+c*s)))).T
    output = tmp
    return output

@njit
def contrib_solve(x,imat,cmat,n,p):
    """Computes a performance for vector x from given contribution matrix and interaction matrix

    Notes:
        Uses Numba's njit compiler.

    Args:
        x: An input vector
        imat (numpy.ndarray): Interaction matrix
        cmat (numpy.ndarray): Contribution matrix
        n (int): Number of tasks per landscape
        p (int): Number of landscapes (population size)

    Returns:
        float: A mean performance of an input vector x given cmat and imat.
    """

    n_p = n*p
    phi = np.zeros(n_p)
    for i in range(n_p):
        tmp = x[np.where(imat[:,i]>0)] # coupled bits
        tmp_loc = np.sum(np.flip(tmp) * 2 ** (np.arange(len(tmp)))) # convert to integer
        phi[i] = cmat[tmp_loc,i]

    output = [0.0]*p
    for i in range(p):
        output[i] = phi[i*n : (i+1)*n].mean()
    return output

def contrib_full(imat,cmat,n,p):
    """Computes performances for all binary vectors of size N*P

    Notes:
        The most processing-heavy part of any simulation. get_global_max() is a 'lazy' alternative

    Args:
        imat (numpy.ndarray): Interaction matrix
        cmat (numpy.ndarray): Contribution matrix
        n (int): Number of tasks per landscape
        p (int): Number of landscapes (population size)

    Returns:
        numpy.ndarray: Array of performances for all vectors of size N*P, normalized by the global maximum
        float: the global maximum value
    """

    n_p = n*p
    perfmat = np.empty((2**n_p,p),dtype=float)
    for i in range(2**n_p):
        bstring = np.array(binx(i,n_p))
        bval = contrib_solve(bstring,imat,cmat,n,p) # !!! processing heavy !!!
        perfmat[i,] = bval
    idxmax = np.argmax(np.mean(perfmat,1))
    perfmax = perfmat[idxmax]

    perfmean = np.mean(perfmat,axis=1)
    perfglobalmax = perfmean.max()
    #perfargmax = perfmean.argmax()

    output1 = perfmat / perfmax
    #output2 = perfargmax
    #output3 = perfglobalmax
    output3 = perfmax
    return output1, output3#, output2, output3

def get_globalmax(imat,cmat,n,p,brute=True,t0=20,t1=0,alpha=0.1,kk=1):
    """Estimates a global maximum for a landscape using Simulated Annealing algorithm

    Args:
        imat (numpy.ndarray): Interaction matrix
        cmat (numpy.ndarray): Contribution matrix
        n (int): Number of tasks per landscape
        p (int): Number of landscapes (population size)
        t0 (float): Initial temperature. Suggested values are 20 (default) or 1.
        t1 (float): Final temperature. Suggested value is 0 (default)
        alpha (float): Step for temperature. Suggested value are 0.1 (default) or 0.001
        kk (float): Adjustment parameter. Not used at the moment.

    Returns:
        numpy.ndarray: The estimates of the global maximum for each of the P landscapes
    """

    n_p = n*p
    output = None

    if brute:
        bstrings = map(lambda x: np.array(binx(x,n_p)), range(2**n_p))

        perfmax = [0.0]*p#np.zeros(p,dtype=float)
        for i in bstrings:#range(2**n_p):
            bstring = i
            bval = contrib_solve(bstring,imat,cmat,n,p) # !!! processing heavy !!!
            if np.mean(bval)>np.mean(perfmax):
                perfmax = bval

        output = np.array(perfmax)
    else:
        state = np.array(binx(0,n_p))
        value = contrib_solve(state,imat,cmat,n,p)

        t = t0
        while t>t1:
            state_ = random_neighbour(state,0,n_p)
            value_ = contrib_solve(state_,imat,cmat,n,p)
            value_ = np.array(value_)

            d_mean = np.mean(value_) - np.mean(value)
            #d_sep = value_ - value

            if d_mean > 0 or np.exp(d_mean/t) > np.random.uniform():
                state = state_
                value = value_
            t -= alpha

        output = value
    return output
###############################################################################

def assign_tasks(N,POP,shape="solo"):
    """Assigns N tasks to POP agents

    Args:
        N (int): Number of tasks
        POP (int): Number of agents (population size)
        shape (str): Type of allocation. Takes values 'solo' (default) and 'overlap' (not used at the moment)

    Returns:
        numpy.ndarray: Returns a POPxN matrix where rows represent agents and cols represent tasks
    """

    output = None
    if shape=="solo":
        perCapita = N / POP
        tmp = np.eye(POP,dtype=int)
        tmp = tmp.repeat(perCapita,axis=1)
        output = tmp
    else:
        print("Task assignment shape unrecognized")
    # print(f"Assignment shape {shape} selected")
    return output

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

def generate_network(POP,S=2,pcom=1.0,shape="random",absval=False):
    """Generates a unidirectional network topology

    Args:
        POP (int): Number of agents (population size)
        S (int): Network degree
        pcom (float): Probability of communicating through the channel
        shape (str): A network topology. Takes values 'random' (default), 'cycle' (ring topology), 'star' (not used at the moment)
        absval (bool): Indexing convention. If True, converts negative indices to positive.

    Returns:
        numpy.ndarray: A POPxPOP matrix with probabilities of connecting to other agents
        list: A list of S-sized vectors for every agent, if shape is 'cycle' (to be fixed)
    """

    if S>=POP:
        print("Error: wrong network degree")
        return 0
    if pcom>1 or pcom<0:
        print("Error: wrong probability of communication")
        return 0
    output = None
    if S == 0:
        output = []
    elif shape=="cycle":
        tmp = [[z-1] + [_%POP for _ in range(z+1,z+S)] for z in range(POP)]
        if absval==True:
            tmp = [[(z-1)%POP] + [_%POP for _ in range(z+1,z+S)] for z in range(POP)]
        output = tmp
    elif shape == "random":
        tmp = random_binary_matrix(POP,S,0)
        output = tmp * pcom
    else:
        print(f"Unrecognized network shape '{shape}'")

    return(output)

def generate_couples(POP,S=2,shape="cycle"):
    """Generates couplings between landscapes (external interaction)

    Args:
        POP (int): Number of landscapes (population size)
        S (int): Number of landscapes considered for external bits
        shape (str): A network topology. Takes values 'cycle' (default)

    Returns:
        list: A list of S-sized vectors with couples for every landscape.
    """

    if S>=POP:
        print("Error: wrong network degree")
        return 0
    output = None
    if S == 0:
        output = []
    elif shape=="cycle":
        tmp = [[(z-1)%POP] + [_%POP for _ in range(z+1,z+S)] for z in range(POP)]
        output = tmp
    else:
        print(f"Unrecognized network shape '{shape}'")
    return(output)
###############################################################################

def get_neighbours(vec,count):
    """Generates binary vectors that are 1-bit away (a unit Hamming distance)

    Args:
        vec (list or numpy.ndarray): An input vector
        count (int): Number of neighbours to generate

    Returns:
        list: A list of 1-bit neighbours of vec
        list: A list of decimal equivalents of the above
    """

    tmpv = []
    tmpi = []
    subbset = np.random.choice(np.arange(len(vec)),count,replace=False)
    for i in subbset:
        y = vec.copy()
        y[i] = 1 - y[i]
        tmpv.append(y)
        tmpi.append(binx(y))
    return(tmpv, tmpi)

def random_neighbour(vec,myid,n):
    """Generates a random binary vector that is 1-bit away (a unit Hamming distance)

    Args:
        vec (list or numpy.ndarray): An input vector
        myid (int): An id of an agent of interest
        n (int): Number of tasks allocated to a single agent

    Returns:
        list: A vector with one bit flipped for agent myid
    """

    tmp = vec.copy()
    rnd = np.random.choice(range(myid*n,(myid+1)*n))
    tmp[rnd] = 1- tmp[rnd]
    output = tmp
    return output

def get_index(vec,myid,n):
    """Gets a decimal equivalent for the bitstring for agent myid

    Args:
        vec (list or numpy.ndarray): An input vector
        myid (int): An id of an agent of interest
        n (int): Number of tasks allocated to a single agent

    Returns:
        int: A decimal value of vec for myid
    """

    tmp = vec.copy()
    tmp = binx(tmp[myid*n:(myid+1)*n])
    output = tmp
    return output

def with_noise(vec,pb=0):
    """Adds a noise to a vector

    Args:
        vec (list or numpy.ndarray): An input vector
        pb (float): Probability of error/noise. Takes values between 0 (default) and 1.

    Returns:
        numpy.ndarray: A vector of same size as vec but with errors
    """

    if vec.size==0:
        return vec
    tmp = vec
    nvec = len(vec)
    noise = np.random.choice(2,p=[1-pb, pb])
    if noise:
        tmp = random_neighbour(vec,0,nvec)
    output = tmp
    return output
###############################################################################

def beta_mean(x,y):
    """Calculates the mean of the Beta(x,y) distribution

    Args:
        x (float): alpha (shape) parameter
        y (float): beta (scale) parameter

    Returns:
        float: A float that returns x/(x+y)
    """

    output = x / (x+y)
    return output

def flatten(x):
    """Converts a list of lists into a flat list

    Args:
        x (list): An input list

    Returns:
        list: A flat list
    """

    output = [_ for z in x for _ in z]
    return output

#def calculate_freq(x,vec):
#    tmp = flatten(vec)
#    if len(tmp)==0:
#        tmp = 0
#    else:
#        tmp = tmp.count(x)/len(tmp)
#    output = tmp
#    return output

def calculate_freq(x,vec):
    """Calclates frequency of vec in x

    Args:
        x (numpy.ndarray): A 2d matrix
        vec (numpy.ndarray): A 1d array

    Returns:
        float: Frequency of vec in x
    """

    if x is None or vec.size==0 or vec[0,0]==-1:
        return 0.0
    freq = np.mean(vec,axis=0)
    tmp1 = np.multiply(x,freq)
    tmp2 = np.multiply(1-x,1-freq)
    tmp = tmp1 + tmp2
    output = np.mean(tmp)
    return output

def extract_soc(x,myid,n,nsoc):
    """Extracts social bits from a bitstring

    Args:
        x (numpy.ndarray): An input vector
        myid (int): An id of an agent of interest
        n (int): Number of tasks allocated to a single agent
        nsoc (int): Number of social tasks (exogeneous)

    Returns:
        numpy.ndarray: A vector of size nsoc
    """

    tmp = x[(myid+1)*n-nsoc:(myid+1)*n]
    #tmp = str(tmp).replace(" ","").replace("[","").replace("]","")
    output = tmp
    return output

#def extract_pub(x,myid,n,p,npub):
#    tmp = np.reshape(x,(p,n))
#    tmp = tmp[:,-npub:n]
#    output = tmp
#    return output

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

def pick(vec,count):
    """Picks random elements form a vector

    Args:
        vec (numpy.ndarray): An input vector
        count (int): Number of elements to pick

    Returns:
        numpy.ndarray: A vector of size count
    """

    tmp = np.random.choice(range(vec),count,replace=False)
    output = np.sort(tmp)
    return output

def artify(n,p,r):
    """depreciated"""

    tmp = np.arange(n*p)
    tmp = tmp.reshape(p,n)
    fnc = lambda z: np.random.choice(z,r,replace=False)
    tmp = np.apply_along_axis(fnc,1,tmp)
    output = tmp
    return output

def calculate_match(x,art):
    """depreciated"""

    if art==[]:
        return 0.0
    tmp = [x[z[0]]==z[1] for z in art]
    output = sum(tmp) / len(tmp)
    return output

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

def cobb_douglas(weights,vec):
    """A Cobb-Douglas utility function with given weights

    Args:
        weights (list): Weights
        vec (numpy.ndarray): An input vector

    Returns:
        float: A Cobb-Douglas value
    """

    x = [z+1 for z in vec]
    w = weights
    tmp = np.power(x,w).prod()
    output = tmp
    return output

def satisfice(x,y,u):
    """A satisficing utility

    Args:
        x (float): first value
        y (float): second value
        u (float): first goal

    Returns:
        float: Utility value
    """

    cond = (x >= u)
    tmp = cond * (y + u) + (1-cond) * x
    output = tmp
    return output

def weighted(x,y,p1,p2):
    """A weighted sum

    Args:
        x (float): first value
        y (float): second value
        p1 (float): first weight
        p2 (float): second weight

    Returns:
        float: A weighted sum
    """

    tmp = p1*x + p2*y
    output = tmp
    return output

def goal_prog(perf1,perf2,u,p1,p2):
    """A goal programming

    Args:
        perf1 (float): first performance
        perf2 (float): second performance
        u (list): goals
        p1 (float): first weigh
        p2 (float): second weight

    Returns:
        float: A GP output
    """

    d1 = np.max((u[0]-perf1,0))
    d2 = np.max((u[1]-perf2,0))
    tmp = p1*d1 + p2*d2
    output = -tmp
    return output

def schism(perf1,perf2,social):
    """depreciated"""

    tmp = None
    if social is True:
        tmp = perf2
    else:
        tmp = perf1
    output = tmp
    return output

max_stars_repo_name: ravshansk/NKPackage

max_issues_repo_name: ravshansk/NKPackage

max_forks_repo_name: ravshansk/NKPackage

reason: import numpy,from scipy,from numba


text:
from numba import jit,vectorize,njit

import numpy as np
import functools as fc
import os
import time
from tqdm.notebook import tqdm_notebook as tqdm
# from tqdm import tqdm

# tqdm  = lambda x : x

# def co_variance(X,bias=0):
#     nx,i = X.shape ## X is input data matrix
#     # ans = np.zeros((i,i))#,dtype=float)
#     # for a in X:
#     #     ans += np.outer(a,a)
#     ans = X.T.dot(X)
#     return ans/(nx-bias)

def co_variance(x,bias=0):
    return x.T.dot(x)/(x.shape[0]-bias)
co_variance.name = "co_variance" ## Abuse of function notations

def ex_variance(cm):
    i,_ = cm.shape
    fv = np.outer(cm,cm)
    return fv.reshape(i,i**3)

def ra_kurtosis(X,bias=0):
    nx,i = X.shape
    ans = np.zeros((i**3,i))
    looper = tqdm(X)
    for a in looper:
        ans += np.outer(np.outer(np.outer(a,a),a),a)
    return ans.T/(nx-bias)
ra_kurtosis.name="Raw_kurtosis"

def co_kurtosis(rand_mat,bias=0):
    ck = ra_kurtosis(rand_mat,bias)
    cm = co_variance(rand_mat,1) #UNBIASED est
    ev = ex_variance(cm)
    return ck- 3*ev
co_kurtosis.name = "co_kurtosis"

@jit
def val_substraction(CK,CV):
    nvar = CV.shape[-1]
    for i in range(nvar):
        for j in range(nvar):
            for k in range(nvar):
                for l in range(nvar):
                    CK[i,j,k,l] = (CK[i,j,k,l] - CV[i,j]*CV[k,l] - CV[i,k]*CV[j,l] - CV[i,l]*CV[j,k])
    return CK

# def outer_Variance()
def val_kurtosis(xscaled):
    n,nvar = xscaled.shape

    CK = ra_kurtosis(xscaled).reshape(nvar,nvar,nvar,nvar)
    CV= co_variance(xscaled)

    CK = val_substraction(CK,CV)
    CK_m = CK.reshape(nvar, nvar*nvar*nvar)
    return CK_m
val_kurtosis.name="val_kurtosis"

## -------------------------------------------------------------------------

### Following function isued for testing accuracy with index definitio of kurtosis
### Warning-~very slow execution without numba jit

@jit
def ex_kurtosis(u):
    nx,nv = u.shape
    mom = np.zeros((nv, nv))#, dtype=float, order='F')
    # compute covariance matrix
    for j in range(nv):
        for i in range(nv):
            for n in range(nx):
                mom[i,j] = mom[i,j] + u[n,i] * u[n,j]
    mom2 = mom/nx

    tmp = np.zeros((nv,nv,nv,nv))#, dtype=float, order='F')
    # compute cokurtosis matrix
    for l in range(nv):
        for k in range(nv):
            for j in range(nv):
                for i in range(nv):
                    for n in range(nx):
                        tmp[i,j,k,l] = tmp[i,j,k,l] + u[n,i] * u[n,j] * u[n,k] * u[n,l]

    tmp=tmp/nx

    for l in range(nv):
        for k in range(nv):
            for j in range(nv):
                for i in range(nv):
                    tmp[i,j,k,l] = tmp[i,j,k,l] - mom2[i,j]*mom2[k,l] - mom2[i,k]*mom2[j,l] - mom2[i,l]*mom2[j,k]

    return tmp.reshape(nv,nv**3)
ex_kurtosis.name="For loop kurtosis"

max_stars_repo_name: kulkarnisp/dim-red

max_issues_repo_name: kulkarnisp/dim-red

max_forks_repo_name: kulkarnisp/dim-red

reason: import numpy,from numba


text:
import wave
import numpy as np
import pylab as plt
f=wave.open("./train/1.wav")
params=f.getparams()
nchannels, sampwidth, framerate, nframes = params[:4]
str_data=f.readframes(nframes)
f.close()
wave_data=np.fromstring(str_data,dtype=np.short)
wave_data.shape=-1,2
wave_data=wave_data.T
time=np.arange(0,nframes)/framerate
print(params)
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(time,wave_data[0])
plt.subplot(2,1,2)
plt.plot(time,wave_data[1],c="r")
plt.xlabel("time")
plt.show()


text:
import os
import os.path
import shutil
import numpy

from amuse.lab import *
"""
from amuse.units import units, constants
from amuse.datamodel import Particles, Particle
from amuse.io import write_set_to_file
"""
from amuse.community.mesa.interface import MESA as stellar_evolution_code
from amuse.ext.star_to_sph import convert_stellar_model_to_SPH
from amuse.ext.relax_sph import relax

def evolve_star(mass, age):
    stellar = MESA()
    stellar.parameters.metallicity = 0.02
    star = stellar.particles.add_particle(Particle(mass=mass))
    stellar.evolve_model(age)
    return stellar

def evolve_star_and_convert_to_sph(stellar, omega, Nsph):
    new_age = stellar.model_time
    star = stellar.particles[0]
    if star.core_mass>zero:
        with_core_particle = True
        target_core_mass  = coremass
    else:
        with_core_particle = False
        target_core_mass  = 0|units.MSun

    m_sph = (star.mass-target_core_mass)/float(Nsph)
    print("N_sph=", Nsph, m_sph.in_(units.MSun))
    print("Target core mass:", target_core_mass)
    model = convert_stellar_model_to_SPH(
        star,
        Nsph,
        with_core_particle = with_core_particle,
        target_core_mass  = target_core_mass,
        do_store_composition = False,
        base_grid_options=dict(type="fcc")
    )
    print("Final star:", star)
    age = stellar.model_time
    M_scale = star.mass
    R_scale = star.radius
    stellar.stop()

    core = model.core_particle
    print("core", core)
    if core is None:
        print("Make zero mass core")
        core = Particle(mass=0|units.MSun, radius=1|units.RSun)
        core.position = (0,0,0) | units.AU
        core.velocity = (0,0,0) | units.kms
    sph_star = model.gas_particles
    print("Add Spin to star:", omega)
    sph_star.add_spin(omega)

    relaxed_sph_star, core_particles = relax_sph_realization(sph_star, core)
    return relaxed_sph_star, core_particles, age

def relax_sph_realization(sph_star, core):

    dynamical_timescale = sph_star.dynamical_timescale()
    converter = nbody_system.nbody_to_si(dynamical_timescale, 1|units.RSun)
    hydro = Gadget2(converter, number_of_workers=2)
    hydro.parameters.time_max = 3 * dynamical_timescale
    hydro.parameters.max_size_timestep = dynamical_timescale / 100
    hydro.parameters.time_limit_cpu = 1.0 | units.Gyr

    t_end_in_t_dyn = 2.5 # Relax for this many dynamical timescales
    t_end = t_end_in_t_dyn * sph_star.dynamical_timescale(mass_fraction=0.9)
    n_steps = 250
    velocity_damp_factor = 1.0 - (2.0*numpy.pi*t_end_in_t_dyn)/n_steps # Critical damping

    hydro.gas_particles.add_particles(sph_star)
    to_framework = hydro.gas_particles.new_channel_to(sph_star)

    if core.mass>0|units.MSun:
        hydro.dm_particles.add_particles(core)
    for i_step, time in enumerate(t_end * numpy.linspace(1.0/n_steps, 1.0, n_steps)):
        hydro.evolve_model(time)
        to_framework.copy()
        hydro.gas_particles.velocity = velocity_damp_factor * hydro.gas_particles.velocity

    return hydro.gas_particles, hydro.dm_particles

def new_option_parser():
    from amuse.units.optparse import OptionParser
    result = OptionParser()
    result.add_option("-N",
                      dest="Nsph", type="int",
                      default = 1000,
                      help="number of SPH particles[%default]")
    result.add_option("-M", unit=units.MSun,
                      dest="mass", type="float",
                      default = 0.6 | units.MSun,
                      help="stellar mass [%default]")
    result.add_option("-o", unit=units.day**-1,
                      dest="omega", type="float",
                      default = 24 | units.day**-1,
                      help="stellar rotation [%default]")
    result.add_option("-t", unit=units.Myr,
                      dest="age", type="float",
                      default = 10 | units.Myr,
                      help="stellar age [%default]")
    return result

if __name__ == "__main__":
    o, arguments = new_option_parser().parse_args()
    stellar = evolve_star(o.mass, o.age)
    star, core, age = evolve_star_and_convert_to_sph(stellar, o.omega, o.Nsph)
    print("age=", age.in_(units.Myr))
    print(star)
    print(star.mass.sum().in_(units.MSun), core.mass)
    filename = 'Hydro_M%2.2dMSun.h5'%o.mass.value_in(units.MSun)
    print(filename)

    write_set_to_file(star, filename, format='hdf5', append_to_file=False)

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import wobble
from time import time
import h5py
import os

if __name__ == "__main__":
    starname = '51peg'
    K_star = 0
    K_t = 3
    niter = 150 # for optimization
    plots = True
    epochs = [0, 50] # to plot
    movies = False

    star_reg_file = '../wobble/regularization/{0}_star_K{1}.hdf5'.format(starname, K_star)
    tellurics_reg_file = '../wobble/regularization/{0}_t_K{1}.hdf5'.format(starname, K_t)
    plot_dir = '../results/plots_{0}_Kstar{1}_Kt{2}/'.format(starname, K_star, K_t)

    if False:
        # quick test on single order
        data = wobble.Data(filename='../data/'+starname+'_e2ds.hdf5', orders=[30, 56])
        results = wobble.Results(data=data)
        for r in range(data.R):
            model = wobble.Model(data, results, r)
            model.add_star('star', variable_bases=K_star,
                            regularization_par_file=star_reg_file,
                            learning_rate_template=0.01, learning_rate_rvs=1.)
            model.add_telluric('tellurics', rvs_fixed=True, variable_bases=K_t,
                                regularization_par_file=tellurics_reg_file,
                                learning_rate_template=0.01)
            wobble.optimize_order(model, niter=niter, save_history=False, rv_uncertainties=False,
                                  template_uncertainties=False, basename='../results/test',
                                  epochs=epochs, movies=movies)
        results.write('../results/test_{0}_Kstar{1}_Kt{2}.hdf5'.format(starname, K_star, K_t))
        assert False

    print("running wobble on star {0} with K_star = {1}, K_t = {2}".format(starname, K_star, K_t))
    start_time = time()
    orders = np.arange(72)
    data = wobble.Data(filename='../data/'+starname+'_e2ds.hdf5', orders=orders)
    results = wobble.Results(data=data)

    print("data loaded")
    print("time elapsed: {0:.2f} min".format((time() - start_time)/60.0))
    elapsed_time = time() - start_time

    if plots:
        print("plots will be saved under directory: {0}".format(plot_dir))
        if not os.path.exists(plot_dir):
            os.makedirs(plot_dir)
    star_learning_rate = 0.1
    telluric_learning_rate = 0.01
    for r,o in enumerate(orders):
        model = wobble.Model(data, results, r)
        model.add_star('star', variable_bases=K_star,
                        regularization_par_file=star_reg_file,
                        learning_rate_template=star_learning_rate)
        model.add_telluric('tellurics', rvs_fixed=True, variable_bases=K_t,
                            regularization_par_file=tellurics_reg_file,
                            learning_rate_template=telluric_learning_rate)
        print("--- ORDER {0} ---".format(o))
        if plots:
            wobble.optimize_order(model, niter=niter, save_history=True,
                                  basename=plot_dir+'history', movies=movies, epochs_to_plot=epochs)
            fig, ax = plt.subplots(1, 1, figsize=(8,5))
            ax.plot(data.dates, results.star_rvs[r] + data.bervs - data.drifts - np.mean(results.star_rvs[r] + data.bervs),
                    'k.', alpha=0.8, ms=4)
            ax.plot(data.dates, data.pipeline_rvs + data.bervs - np.mean(data.pipeline_rvs + data.bervs),
                    'r.', alpha=0.5, ms=4)
            ax.set_ylabel('RV (m/s)', fontsize=14)
            ax.set_xlabel('BJD', fontsize=14)
            plt.savefig(plot_dir+'results_rvs_o{0}.png'.format(o))
            plt.close(fig)
            for e in epochs:
                results.plot_spectrum(r, e, data, plot_dir+'results_synth_o{0}_e{1}.png'.format(o, e))
        else:
            wobble.optimize_order(model, niter=niter)
        del model # not sure if this does anything
        print("order {1} optimization finished. time elapsed: {0:.2f} min".format((time() - start_time)/60.0, o))
        print("this order took {0:.2f} min".format((time() - start_time - elapsed_time)/60.0))
        elapsed_time = time() - start_time

    print("all orders optimized.")
    print("time elapsed: {0:.2f} minutes".format((time() - start_time)/60.0))

    # do post-processing:
    results.combine_orders('star')
    results.apply_drifts('star')
    results.apply_bervs('star')

    if plots:
        fig, (ax, ax2) = plt.subplots(2, 1, gridspec_kw = {'height_ratios':[3, 1]})
        ax.scatter(data.dates, data.pipeline_rvs - np.mean(data.pipeline_rvs),
                    c='r', label='DRS', alpha=0.7, s=12)
        ax.scatter(data.dates, results.star_time_rvs - np.mean(results.star_time_rvs),
                    c='k', label='wobble', alpha=0.7, s=12)
        ax.legend()
        ax.set_xticklabels([])
        ax2.scatter(data.dates, results.star_time_rvs - data.pipeline_rvs, c='k', s=12)
        ax2.set_ylabel('JD')
        fig.tight_layout()
        fig.subplots_adjust(hspace=0.05)
        plt.savefig(plot_dir+'results_rvs.png')
        plt.close(fig)

        fig, (ax, ax2) = plt.subplots(2, 1, gridspec_kw = {'height_ratios':[3, 1]})
        ax.scatter(data.dates % 4.2308, data.pipeline_rvs - np.mean(data.pipeline_rvs),
                    c='r', label='DRS', alpha=0.7, s=12)
        ax.scatter(data.dates % 4.2308, results.star_time_rvs - np.mean(results.star_time_rvs),
                    c='k', label='wobble', alpha=0.7, s=12)
        ax.legend()
        ax.set_xticklabels([])
        ax2.scatter(data.dates, results.star_time_rvs - data.pipeline_rvs, c='k', s=12)
        ax2.set_ylabel('JD')
        fig.tight_layout()
        fig.subplots_adjust(hspace=0.05)
        plt.savefig(plot_dir+'results_rvs_phased.png')
        plt.close(fig)

    print("final RVs calculated.")
    print("time elapsed: {0:.2f} minutes".format((time() - start_time)/60.0))

    # save output:
    results_file = '../results/results_{0}_Kstar{1}_Kt{2}.hdf5'.format(starname, K_star, K_t)
    results.write(results_file)
    star_rvs_file = '../results/rvs_{0}_Kstar{1}_Kt{2}.txt'.format(starname, K_star, K_t)
    results.write_rvs('star', star_rvs_file, all_orders=True)

    print("results saved as: {0} & {1}".format(results_file, star_rvs_file))
    print("-----------------------------")
    print("HARPS pipeline std = {0:.3f} m/s".format(np.std(data.pipeline_rvs)))
    print("wobble std = {0:.3f} m/s".format(np.std(results.star_time_rvs)))
    print("total runtime:{0:.2f} minutes".format((time() - start_time)/60.0))

max_stars_repo_name: exord/wobble

max_issues_repo_name: exord/wobble

max_forks_repo_name: exord/wobble


text:
import pyccl as ccl
import numpy as np
import scipy.integrate
import math
from functools import partial

cosmoin=ccl.Cosmology(
    Omega_c=0.25,
    Omega_b=0.05,
    h=0.7,
    sigma8=0.8,
    n_s=0.96,
    Neff=0,
    m_nu=0.0,
    w0=-1.,
    wa=0.,
    T_CMB=2.7,
    m_nu_type='normal',
    Omega_g=0,
    Omega_k=0,
    transfer_function='bbks',
    matter_power_spectrum='linear')

k_lin=np.logspace(-4., 3., 10000)
a_arr = np.logspace(-1., 0., 10)
knl_arr = np.zeros((len(a_arr),))
for i in range(len(a_arr)):
    a = a_arr[i]
    #I will need linear power at z
    pk_lin_z=ccl.linear_matter_power(cosmoin, k_lin, a)
    interp_pk_lin_z=scipy.interpolate.interp1d(k_lin,pk_lin_z)

    #---kNL prediction----
    [intknl,errintknl]=scipy.integrate.quad(interp_pk_lin_z,min(k_lin),max(k_lin),epsabs=0,epsrel=1e-6)
    print("Check error is small:",errintknl/intknl)
    knlemin2=1./(6*math.pi**2)*intknl
    knl=1./math.sqrt(knlemin2)
    print('knl[1/Mpc]=',knl)
    knl_arr[i] = knl

header_knl="[0] a, [1] k_NL Mpc^-1"

np.savetxt("kNL.txt", np.transpose(np.vstack((a_arr ,knl_arr))), header=header_knl)

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL

reason: import numpy,import scipy


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

# Lint as: python3
"""Graph construction utility functions.

Functions for graph manipulation and creation.
TODO(tsitsulin): add headers, tests, and improve style.
"""

import pickle
import sys

import networkx as nx
import numpy as np
import scipy.sparse
from scipy.sparse import csr_matrix
from scipy.sparse import lil_matrix
from scipy.sparse.base import spmatrix
from sklearn.neighbors import kneighbors_graph
import tensorflow as tf

def construct_knn_graph(data, k = 15, symmetrize = True):
  graph = kneighbors_graph(data, k)
  if symmetrize:
    graph = graph + graph.T
    graph.data = np.ones(graph.data.shape)
  return graph

def normalize_graph(graph,  # pylint: disable=missing-function-docstring
                    normalized,
                    add_self_loops = True):
  if add_self_loops:  # Bröther may i have some self-lööps
    graph = graph + scipy.sparse.identity(graph.shape[0])
  degree = np.squeeze(np.asarray(graph.sum(axis=1)))
  if normalized:
    with np.errstate(divide='ignore'):
      degree = 1. / np.sqrt(degree)
    degree[degree == np.inf] = 0
    degree = scipy.sparse.diags(degree)
    return degree @ graph @ degree
  else:
    with np.errstate(divide='ignore'):
      degree = 1. / degree
    degree[degree == np.inf] = 0
    degree = scipy.sparse.diags(degree)
    return degree @ graph

def scipy_to_tf(matrix):
  matrix = matrix.tocoo()
  return tf.sparse.SparseTensor(
      np.vstack([matrix.row, matrix.col]).T, matrix.data.astype(np.float32),
      matrix.shape)

def load_npz_to_sparse_graph(file_name):  # pylint: disable=missing-function-docstring
  with np.load(open(file_name, 'rb'), allow_pickle=True) as loader:
    loader = dict(loader)
    adj_matrix = csr_matrix(
        (loader['adj_data'], loader['adj_indices'], loader['adj_indptr']),
        shape=loader['adj_shape'])

    if 'attr_data' in loader:
      # Attributes are stored as a sparse CSR matrix
      attr_matrix = csr_matrix(
          (loader['attr_data'], loader['attr_indices'], loader['attr_indptr']),
          shape=loader['attr_shape']).todense()
    elif 'attr_matrix' in loader:
      # Attributes are stored as a (dense) np.ndarray
      attr_matrix = loader['attr_matrix']
    else:
      raise Exception('No attributes in the data file', file_name)

    if 'labels_data' in loader:
      # Labels are stored as a CSR matrix
      labels = csr_matrix((loader['labels_data'], loader['labels_indices'],
                           loader['labels_indptr']),
                          shape=loader['labels_shape'])
      label_mask = labels.nonzero()[0]
      labels = labels.nonzero()[1]
    elif 'labels' in loader:
      # Labels are stored as a numpy array
      labels = loader['labels']
      label_mask = np.ones(labels.shape, dtype=np.bool)
    else:
      raise Exception('No labels in the data file', file_name)

  return adj_matrix, attr_matrix, labels, label_mask

def _parse_index_file(filename):
  index = []
  for line in open(filename):
    index.append(int(line.strip()))
  return index

def _sample_mask(idx, l):
  """Create mask."""
  mask = np.zeros(l)
  mask[idx] = 1
  return np.array(mask, dtype=np.bool)

def load_kipf_data(path_str, dataset_str):  # pylint: disable=missing-function-docstring
  names = ['x', 'y', 'tx', 'ty', 'allx', 'ally', 'graph']
  objects = []
  for i in range(len(names)):
    with open('{}/ind.{}.{}'.format(path_str, dataset_str, names[i]),
              'rb') as f:
      if sys.version_info > (3, 0):
        objects.append(pickle.load(f, encoding='latin1'))
      else:
        objects.append(pickle.load(f))

  x, y, tx, ty, allx, ally, graph = tuple(objects)  # pylint: disable=unbalanced-tuple-unpacking
  test_idx_reorder = _parse_index_file('{}/ind.{}.test.index'.format(
      path_str, dataset_str))
  test_idx_range = np.sort(test_idx_reorder)

  if dataset_str == 'citeseer':
    # Fix citeseer dataset (there are some isolated nodes in the graph)
    # Find isolated nodes, add them as zero-vecs into the right position
    test_idx_range_full = range(
        min(test_idx_reorder),
        max(test_idx_reorder) + 1)
    tx_extended = lil_matrix((len(test_idx_range_full), x.shape[1]))
    tx_extended[test_idx_range - min(test_idx_range), :] = tx
    tx = tx_extended
    ty_extended = np.zeros((len(test_idx_range_full), y.shape[1]))
    ty_extended[test_idx_range - min(test_idx_range), :] = ty
    ty = ty_extended

  features = scipy.sparse.vstack((allx, tx)).tolil()
  features[test_idx_reorder, :] = features[test_idx_range, :]
  adj = nx.adjacency_matrix(nx.from_dict_of_lists(graph))

  labels = np.vstack((ally, ty))
  labels[test_idx_reorder, :] = labels[test_idx_range, :]

  idx_test = test_idx_range.tolist()
  idx_train = range(len(y))
  idx_val = range(len(y), len(y) + 500)

  train_mask = _sample_mask(idx_train, labels.shape[0])
  val_mask = _sample_mask(idx_val, labels.shape[0])
  test_mask = _sample_mask(idx_test, labels.shape[0])

  y_train = np.zeros(labels.shape)
  y_val = np.zeros(labels.shape)
  y_test = np.zeros(labels.shape)
  y_train[train_mask, :] = labels[train_mask, :]
  y_val[val_mask, :] = labels[val_mask, :]
  y_test[test_mask, :] = labels[test_mask, :]

  labels = (y_train + y_val + y_test).nonzero()[1]
  label_mask = (y_train + y_val + y_test).nonzero()[0]

  return adj, features.todense(), labels, label_mask

reason: import numpy,import scipy,from scipy,import networkx


text:
"""Flexible method to load confounds generated by fMRIprep.

Authors: Hanad Sharmarke, Dr. Pierre Bellec, Francois Paugam
"""
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
import warnings
import os

# Global variables listing the admissible types of noise components
all_confounds = ["motion", "high_pass", "wm_csf", "global", "compcor"]

def _add_suffix(params, model):
    """
    Add suffixes to a list of parameters.
    Suffixes includes derivatives, power2 and full
    """
    params_full = params.copy()
    suffix = {
        "basic": {},
        "derivatives": {"derivative1"},
        "power2": {"power2"},
        "full": {"derivative1", "power2", "derivative1_power2"},
    }
    for par in params:
        for suff in suffix[model]:
            params_full.append(f"{par}_{suff}")
    return params_full

def _check_params(confounds_raw, params):
    """Check that specified parameters can be found in the confounds."""
    for par in params:
        if not par in confounds_raw.columns:
            raise ValueError(
                f"The parameter {par} cannot be found in the available confounds. You may want to use a different denoising strategy'"
            )

    return None

def _find_confounds(confounds_raw, keywords):
    """Find confounds that contain certain keywords."""
    list_confounds = []
    for key in keywords:
        key_found = False
        for col in confounds_raw.columns:
            if key in col:
                list_confounds.append(col)
                key_found = True
        if not key_found:
            raise ValueError(f"could not find any confound with the key {key}")
    return list_confounds

def _load_global(confounds_raw, global_signal):
    """Load the regressors derived from the global signal."""
    global_params = _add_suffix(["global_signal"], global_signal)
    _check_params(confounds_raw, global_params)
    return confounds_raw[global_params]

def _load_wm_csf(confounds_raw, wm_csf):
    """Load the regressors derived from the white matter and CSF masks."""
    wm_csf_params = _add_suffix(["csf", "white_matter"], wm_csf)
    _check_params(confounds_raw, wm_csf_params)
    return confounds_raw[wm_csf_params]

def _load_high_pass(confounds_raw):
    """Load the high pass filter regressors."""
    high_pass_params = _find_confounds(confounds_raw, ["cosine"])
    return confounds_raw[high_pass_params]

def _label_compcor(confounds_raw, compcor_suffix, n_compcor):
    """Builds list for the number of compcor components."""
    compcor_cols = []
    for nn in range(n_compcor + 1):
        nn_str = str(nn).zfill(2)
        compcor_col = compcor_suffix + "_comp_cor_" + nn_str
        if compcor_col not in confounds_raw.columns:
            warnings.warn(f"could not find any confound with the key {compcor_col}")
        else:
            compcor_cols.append(compcor_col)

    return compcor_cols

def _load_compcor(confounds_raw, compcor, n_compcor):
    """Load compcor regressors."""
    if compcor == "anat":
        compcor_cols = _label_compcor(confounds_raw, "a", n_compcor)

    if compcor == "temp":
        compcor_cols = _label_compcor(confounds_raw, "t", n_compcor)

    if compcor == "full":
        compcor_cols = _label_compcor(confounds_raw, "a", n_compcor)
        compcor_cols.extend(_label_compcor(confounds_raw, "t", n_compcor))

    compcor_cols.sort()
    _check_params(confounds_raw, compcor_cols)
    return confounds_raw[compcor_cols]

def _load_motion(confounds_raw, motion, n_motion):
    """Load the motion regressors."""
    motion_params = _add_suffix(
        ["trans_x", "trans_y", "trans_z", "rot_x", "rot_y", "rot_z"], motion
    )
    _check_params(confounds_raw, motion_params)
    confounds_motion = confounds_raw[motion_params]

    # Optionally apply PCA reduction
    if n_motion > 0:
        confounds_motion = _pca_motion(confounds_motion, n_components=n_motion)

    return confounds_motion

def _pca_motion(confounds_motion, n_components):
    """Reduce the motion paramaters using PCA."""
    n_available = confounds_motion.shape[1]
    if n_components > n_available:
        raise ValueError(
            f"User requested n_motion={n_components} motion components, but found only {n_available}."
        )
    confounds_motion = confounds_motion.dropna()
    confounds_motion_std = scale(
        confounds_motion, axis=0, with_mean=True, with_std=True
    )
    pca = PCA(n_components=n_components)
    motion_pca = pd.DataFrame(pca.fit_transform(confounds_motion_std))
    motion_pca.columns = ["motion_pca_" + str(col + 1) for col in motion_pca.columns]
    return motion_pca

def _sanitize_strategy(strategy):
    """Defines the supported denoising strategies."""
    if isinstance(strategy, list):
        for conf in strategy:
            if not conf in all_confounds:
                raise ValueError(f"{conf} is not a supported type of confounds.")
    else:
        raise ValueError("strategy needs to be a list of strings")
    return strategy

def _confounds_to_df(confounds_raw):
    """Load raw confounds as a pandas DataFrame."""
    if not isinstance(confounds_raw, pd.DataFrame):
        if "nii" in confounds_raw[-6:]:
            suffix = "_space-" + confounds_raw.split("space-")[1]
            confounds_raw = confounds_raw.replace(
                suffix, "_desc-confounds_timeseries.tsv",
            )
            # fmriprep has changed the file suffix between v20.1.1 and v20.2.0 with respect to BEP 012.
            # cf. https://neurostars.org/t/naming-change-confounds-regressors-to-confounds-timeseries/17637
            # Check file with new naming scheme exists or replace, for backward compatibility.
            if not os.path.exists(confounds_raw):
                confounds_raw = confounds_raw.replace(
                    "_desc-confounds_timeseries.tsv", "_desc-confounds_regressors.tsv",
                )

        confounds_raw = pd.read_csv(confounds_raw, delimiter="\t", encoding="utf-8")
    return confounds_raw

def _sanitize_confounds(confounds_raw):
    """Make sure the inputs are in the correct format."""
    # we want to support loading a single set of confounds, instead of a list
    # so we hack it
    flag_single = isinstance(confounds_raw, str) or isinstance(
        confounds_raw, pd.DataFrame
    )
    if flag_single:
        confounds_raw = [confounds_raw]

    return confounds_raw, flag_single

def _confounds_to_ndarray(confounds, demean):
    """Convert confounds from a pandas dataframe to a numpy array."""
    # Convert from DataFrame to numpy ndarray
    labels = confounds.columns
    confounds = confounds.values

    # Derivatives have NaN on the first row
    # Replace them by estimates at second time point,
    # otherwise nilearn will crash.
    mask_nan = np.isnan(confounds[0, :])
    confounds[0, mask_nan] = confounds[1, mask_nan]

    # Optionally demean confounds
    if demean:
        confounds = scale(confounds, axis=0, with_std=False)

    return confounds, labels

class Confounds:
    """
    Confounds from fmriprep

    Parameters
    ----------
    strategy : list of strings
        The type of noise confounds to include.
        "motion" head motion estimates.
        "high_pass" discrete cosines covering low frequencies.
        "wm_csf" confounds derived from white matter and cerebrospinal fluid.
        "global" confounds derived from the global signal.

    motion : string, optional
        Type of confounds extracted from head motion estimates.
        "basic" translation/rotation (6 parameters)
        "power2" translation/rotation + quadratic terms (12 parameters)
        "derivatives" translation/rotation + derivatives (12 parameters)
        "full" translation/rotation + derivatives + quadratic terms + power2d derivatives (24 parameters)

    n_motion : float
        Number of pca components to keep from head motion estimates.
        If the parameters is strictly comprised between 0 and 1, a principal component
        analysis is applied to the motion parameters, and the number of extracted
        components is set to exceed `n_motion` percent of the parameters variance.
        If the n_components = 0, then no PCA is performed.

    wm_csf : string, optional
        Type of confounds extracted from masks of white matter and cerebrospinal fluids.
        "basic" the averages in each mask (2 parameters)
        "power2" averages and quadratic terms (4 parameters)
        "derivatives" averages and derivatives (4 parameters)
        "full" averages + derivatives + quadratic terms + power2d derivatives (8 parameters)

    global_signal : string, optional
        Type of confounds extracted from the global signal.
        "basic" just the global signal (1 parameter)
        "power2" global signal and quadratic term (2 parameters)
        "derivatives" global signal and derivative (2 parameters)
        "full" global signal + derivatives + quadratic terms + power2d derivatives (4 parameters)

    compcor : string,optional
        Type of confounds extracted from a component based noise correction method
        "anat" noise components calculated using anatomical compcor
        "temp" noise components calculated using temporal compcor
        "full" noise components calculated using both temporal and anatomical

    n_compcor : int, optional
        The number of noise components to be extracted.

    demean : boolean, optional
        If True, the confounds are standardized to a zero mean (over time).
        This step is critical if the confounds are regressed out of time series
        using nilearn with no or zscore standardization, but should be turned off
        with "spc" normalization.

    Attributes
    ----------
    `confounds_` : ndarray
        The confounds loaded using the specified model

    `columns_`: list of str
        The labels of the different confounds

    Notes
    -----
    The predefined strategies implemented in this class are from
    adapted from (Ciric et al. 2017). Band-pass filter is replaced
    by high-pass filter, as high frequencies have been shown to carry
    meaningful signal for connectivity analysis.

    References
    ----------
    Ciric et al., 2017 "Benchmarking of participant-level confound regression
    strategies for the control of motion artifact in studies of functional
    connectivity" Neuroimage 154: 174-87
    https://doi.org/10.1016/j.neuroimage.2017.03.020
    """

    def __init__(
        self,
        strategy=["motion", "high_pass", "wm_csf"],
        motion="full",
        n_motion=0,
        wm_csf="basic",
        global_signal="basic",
        compcor="anat",
        n_compcor=10,
        demean=True,
    ):
        """Default parameters."""
        self.strategy = _sanitize_strategy(strategy)
        self.motion = motion
        self.n_motion = n_motion
        self.wm_csf = wm_csf
        self.global_signal = global_signal
        self.compcor = compcor
        self.n_compcor = n_compcor
        self.demean = demean

    def load(self, confounds_raw):
        """
        Load fMRIprep confounds

        Parameters
        ----------
        confounds_raw : Pandas Dataframe or path to tsv file(s), optionally as a list.
            Raw confounds from fmriprep

        Returns
        -------
        confounds :  ndarray or list of ndarray
            A reduced version of fMRIprep confounds based on selected strategy and flags.
            An intercept is automatically added to the list of confounds.
        """
        confounds_raw, flag_single = _sanitize_confounds(confounds_raw)
        confounds_out = []
        columns_out = []
        for file in confounds_raw:
            conf, col = self._load_single(file)
            confounds_out.append(conf)
            columns_out.append(col)

        # If a single input was provided,
        # send back a single output instead of a list
        if flag_single:
            confounds_out = confounds_out[0]
            columns_out = columns_out[0]

        self.confounds_ = confounds_out
        self.columns_ = columns_out
        return confounds_out

    def _load_single(self, confounds_raw):
        """Load a single confounds file from fmriprep."""
        # Convert tsv file to pandas dataframe
        confounds_raw = _confounds_to_df(confounds_raw)

        confounds = pd.DataFrame()

        if "motion" in self.strategy:
            confounds_motion = _load_motion(confounds_raw, self.motion, self.n_motion)
            confounds = pd.concat([confounds, confounds_motion], axis=1)

        if "high_pass" in self.strategy:
            confounds_high_pass = _load_high_pass(confounds_raw)
            confounds = pd.concat([confounds, confounds_high_pass], axis=1)

        if "wm_csf" in self.strategy:
            confounds_wm_csf = _load_wm_csf(confounds_raw, self.wm_csf)
            confounds = pd.concat([confounds, confounds_wm_csf], axis=1)

        if "global" in self.strategy:
            confounds_global_signal = _load_global(confounds_raw, self.global_signal)
            confounds = pd.concat([confounds, confounds_global_signal], axis=1)

        if "compcor" in self.strategy:
            confounds_compcor = _load_compcor(
                confounds_raw, self.compcor, self.n_compcor
            )
            confounds = pd.concat([confounds, confounds_compcor], axis=1)

        confounds, labels = _confounds_to_ndarray(confounds, self.demean)

        return confounds, labels


text:
"""

 Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
 SPDX-License-Identifier: MIT-0

 Permission is hereby granted, free of charge, to any person obtaining a copy of this
 software and associated documentation files (the "Software"), to deal in the Software
 without restriction, including without limitation the rights to use, copy, modify,
 merge, publish, distribute, sublicense, and/or sell copies of the Software, and to
 permit persons to whom the Software is furnished to do so.

 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
 INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
 OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

"""

import tensorflow as tf
import argparse
import os
import numpy as np
import json

# for data processing
def _load_training_data(base_dir):
    """ load training data """
    df_train = np.load(os.path.join(base_dir, 'train.npy'))
    user_train, item_train, y_train = np.split(np.transpose(df_train).flatten(), 3)
    return user_train, item_train, y_train

def batch_generator(x, y, batch_size, n_batch, shuffle, user_dim, item_dim):
    """ batch generator to supply data for training and testing """

    user_df, item_df = x

    counter = 0
    training_index = np.arange(user_df.shape[0])

    if shuffle:
        np.random.shuffle(training_index)

    while True:
        batch_index = training_index[batch_size*counter:batch_size*(counter+1)]
        user_batch = tf.one_hot(user_df[batch_index], depth=user_dim)
        item_batch = tf.one_hot(item_df[batch_index], depth=item_dim)
        y_batch = y[batch_index]
        counter += 1
        yield [user_batch, item_batch], y_batch

        if counter == n_batch:
            if shuffle:
                np.random.shuffle(training_index)
            counter = 0

# network
def _get_user_embedding_layers(inputs, emb_dim):
    """ create user embeddings """
    user_gmf_emb = tf.keras.layers.Dense(emb_dim, activation='relu')(inputs)

    user_mlp_emb = tf.keras.layers.Dense(emb_dim, activation='relu')(inputs)

    return user_gmf_emb, user_mlp_emb

def _get_item_embedding_layers(inputs, emb_dim):
    """ create item embeddings """
    item_gmf_emb = tf.keras.layers.Dense(emb_dim, activation='relu')(inputs)

    item_mlp_emb = tf.keras.layers.Dense(emb_dim, activation='relu')(inputs)

    return item_gmf_emb, item_mlp_emb

def _gmf(user_emb, item_emb):
    """ general matrix factorization branch """
    gmf_mat = tf.keras.layers.Multiply()([user_emb, item_emb])

    return gmf_mat

def _mlp(user_emb, item_emb, dropout_rate):
    """ multi-layer perceptron branch """
    def add_layer(dim, input_layer, dropout_rate):
        hidden_layer = tf.keras.layers.Dense(dim, activation='relu')(input_layer)

        if dropout_rate:
            dropout_layer = tf.keras.layers.Dropout(dropout_rate)(hidden_layer)
            return dropout_layer

        return hidden_layer

    concat_layer = tf.keras.layers.Concatenate()([user_emb, item_emb])

    dropout_l1 = tf.keras.layers.Dropout(dropout_rate)(concat_layer)

    dense_layer_1 = add_layer(64, dropout_l1, dropout_rate)

    dense_layer_2 = add_layer(32, dense_layer_1, dropout_rate)

    dense_layer_3 = add_layer(16, dense_layer_2, None)

    dense_layer_4 = add_layer(8, dense_layer_3, None)

    return dense_layer_4

def _neuCF(gmf, mlp, dropout_rate):
    concat_layer = tf.keras.layers.Concatenate()([gmf, mlp])

    output_layer = tf.keras.layers.Dense(1, activation='sigmoid')(concat_layer)

    return output_layer

def build_graph(user_dim, item_dim, dropout_rate=0.25):
    """ neural collaborative filtering model """

    user_input = tf.keras.Input(shape=(user_dim))
    item_input = tf.keras.Input(shape=(item_dim))

    # create embedding layers
    user_gmf_emb, user_mlp_emb = _get_user_embedding_layers(user_input, 32)
    item_gmf_emb, item_mlp_emb = _get_item_embedding_layers(item_input, 32)

    # general matrix factorization
    gmf = _gmf(user_gmf_emb, item_gmf_emb)

    # multi layer perceptron
    mlp = _mlp(user_mlp_emb, item_mlp_emb, dropout_rate)

    # output
    output = _neuCF(gmf, mlp, dropout_rate)

    # create the model
    model = tf.keras.Model(inputs=[user_input, item_input], outputs=output)

    return model

def model(x_train, y_train, n_user, n_item, num_epoch, batch_size):

    num_batch = np.ceil(x_train[0].shape[0]/batch_size)

    # build graph
    model = build_graph(n_user, n_item)

    # compile and train
    optimizer = tf.keras.optimizers.Adam(learning_rate=1e-3)

    model.compile(optimizer=optimizer,
                  loss=tf.keras.losses.BinaryCrossentropy(),
                  metrics=['accuracy'])

    model.fit_generator(
        generator=batch_generator(
            x=x_train, y=y_train,
            batch_size=batch_size, n_batch=num_batch,
            shuffle=True, user_dim=n_user, item_dim=n_item),
        epochs=num_epoch,
        steps_per_epoch=num_batch,
        verbose=2
    )

    return model

def _parse_args():
    parser = argparse.ArgumentParser()

    parser.add_argument('--model_dir', type=str)
    parser.add_argument('--sm-model-dir', type=str, default=os.environ.get('SM_MODEL_DIR'))
    parser.add_argument('--train', type=str, default=os.environ.get('SM_CHANNEL_TRAINING'))
    parser.add_argument('--hosts', type=list, default=json.loads(os.environ.get('SM_HOSTS')))
    parser.add_argument('--current-host', type=str, default=os.environ.get('SM_CURRENT_HOST'))
    parser.add_argument('--epochs', type=int, default=3)
    parser.add_argument('--batch_size', type=int, default=256)
    parser.add_argument('--n_user', type=int)
    parser.add_argument('--n_item', type=int)

    return parser.parse_known_args()

if __name__ == "__main__":
    args, unknown = _parse_args()

    # load data
    user_train, item_train, train_labels = _load_training_data(args.train)

    # build model
    ncf_model = model(
        x_train=[user_train, item_train],
        y_train=train_labels,
        n_user=args.n_user,
        n_item=args.n_item,
        num_epoch=args.epochs,
        batch_size=args.batch_size
    )

    if args.current_host == args.hosts[0]:
        # save model to an S3 directory with version number '00000001'
        ncf_model.save(os.path.join(args.sm_model_dir, '000000001'), 'neural_collaborative_filtering.h5')


text:
# !/usr/bin/env python3
import sys
import pandas as pd
import numpy as np
import csv
import sklearn
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

with open('train.csv', 'w') as file:
    for line in sys.stdin:
        file.write(line)

#train = pd.read_csv('train.csv', encoding='unicode_escape')
#test = pd.read_csv('Test.csv', encoding='unicode_escape')

#X_train = train.iloc[:, :48]
#Y_train = train.iloc[:, 48:]

#X_train, X_test, Y_train, Y_test = train_test_split(X_train, Y_train, test_size=0.2)

train = pd.read_csv('train.csv',encoding='unicode_escape')
test = pd.read_csv('test_norm.csv',encoding='unicode_escape')
#train = pd.read_csv('Train.csv',encoding='unicode_escape')
train1 = train.head(40941)
#train1 = train.head(20)
#test = train.tail(15)
X_train = train1.iloc[:, 1:49]
Y_train = train1.iloc[:, 49:]
#print(Y_train)
Y_train = Y_train.astype({'48':'int32'})
#print(Y_train)
test = test.iloc[:, 1:49]
#print(Y_train.typ)
#X_test = test.iloc[:, 1:49]
#Y_test = test.iloc[:, 49:]
#Y_test = Y_test.astype({'48':'int32'})
#Y_test.to_csv('y_pred.csv',encoding='utf-8')

# x_train = []
# y_train = []
# x_test = []
# for line in sys.stdin:
#    line = line.rstrip()
#    list_of_elements = line.split(',')
#    list_of_elements = list(map(float,list_of_elements))
#    if len(list_of_elements) == 49:
#        x_train.append(list_of_elements[:47])
#        y_train.append(list_of_elements[48])
#        if len(list_of_elements) == 48:
#            x_test.append(list_of_elements[:47])
#        X_train = np.array(x_train)
#        Y_train = np.array(y_train)
#        X_test = np.array(x_test)
#        Y_train = Y_train.reshape(Y_train.shape[0],1)

X_train = X_train.to_numpy()
Y_train = Y_train.to_numpy()

# X_test = test.to_numpy()

# function that calculates euclidean distance using tile
def euclidean(X_test, center):
    z = pd.np.tile(center, (len(X_test), 1))
    ans = np.subtract(X_test, z)
    ans = pow(ans, 2)
    ans = np.sum(ans, axis=1)
    ans = np.sqrt(ans)
    return ans

def normalize(X):
    scaler = MinMaxScaler()
    scaler.fit(X)
    norm_X = scaler.transform(X)
    return norm_X

# K-NEAREST NEIGHBOR FUNCTION
def KNN(X_train, X_test, Y_train):
    """
  :type X_train: numpy.ndarray
  :type X_test: numpy.ndarray
  :type Y_train: numpy.ndarray

  :rtype: numpy.ndarray
  """
    #X_train = normalize(X_train)
    #X_test = normalize(X_test)
    i = 0
    out = []
    labels = []
    for train_sample in X_train:
        ans = euclidean(X_test, train_sample)
        out.append(ans)
        labels.append(Y_train[i][0])
        i += 1
    out1 = np.array(out).transpose()
    #print(out1.shape)
    #print(len(labels))
    #print(len(out1))
    #print(len(out1[0]))
    for i in range(len(out1)):
        for j in range(len(out1[0])):
            #print(i,j)
            #print("%s,%s,%s"%(str(i), out1[i,j],str(labels[j])))
            print(str(i), out1[i,j], str(labels[j]))

        #t = " ".join(str(a) for a in ans)
        #print(t + " " + str(Y_train[i][0]))

KNN(X_train, test, Y_train)


text:
'''
    Model constructor of an affine model
'''

from sympy import sympify  # Concert numbers in symbpy objects
from sympy import printing  # print sympy expressions in c++
from sympy import diff  # differentiation
from sympy import Poly  # Polynomials
from sympy import polys  # Polynomials
from .modelconstructor import modelConstructorSolver

class modelConstructorAffine(modelConstructorSolver):
    def __init__(self, mn, n, m, _path='./'):
        modelConstructorSolver.__init__(self, mn, n, m, _path)
        self.aus = []
        self.gs = []
        self.hs = []
        self.g0 = sympify(0)
        self.h0s = []
        self.adim = 0

    def init(self, f, rc, fc):
        modelConstructorSolver.init(self, f, rc, fc)
        res = []

        # Here we compute in which components of the
        # control u the vector field and the running cost
        # are linear, the vector res will have
        # the value of one for each affine control
        for u in self.control_symbols_:
            res.append(1)
            c0 = Poly(self.f0, u).all_coeffs()
            for fi in self.funcion_symbols_:
                try:
                    cs = Poly(fi, u).all_coeffs()
                except polys.polyerrors.PolynomialError:
                    cs = [0, 0, 0, 0]
                n = max(len(cs), len(c0))
                if (n > 2):
                    res[-1] = 0

        for i in res:
            if (i == 1):
                self.adim = self.adim + 1

        print('affine vector: ' + str(res))

        self.gs = self.funcion_symbols_[:]
        for i, gi in zip(range(0, self.state_dim_), self.gs):
            for j, uj in zip(
                    range(0, self.control_dim_), self.control_symbols_):
                if (res[j]):
                    cs = Poly(gi, uj).all_coeffs()
                    gi = cs[-1]
                    if (len(cs) == 2):
                        self.hs.append(cs[0])
                    else:
                        self.hs.append(sympify(0))
            self.gs[i] = gi

        self.g0 = self.f0
        for j, uj in zip(range(0, self.control_dim_), self.control_symbols_):
            if (res[j]):
                cs = Poly(self.g0, uj).all_coeffs()
                self.g0 = cs[-1]
                if (len(cs) == 2):
                    self.h0s.append(cs[0])
                else:
                    self.h0s.append(sympify(0))

    def writeAffineFuncs(self):
        self.output_file_.write("""
    void dynSysAffinePart
        (double t,const double x[],const double u2[],double g[]){\n""")
        i = 0
        for gi in self.gs:
            self.output_file_.write("\t\tg[%(n)i]= %(c)s;\n" % {
                'n': i,
                'c': printing.ccode(gi)
            })
            i = i + 1

        self.output_file_.write("\t}\n")

        self.output_file_.write("""
        void dynSysLinearPart(double t,const double x[],double h[]){\n""")
        counter = 0
        for ex in self.hs:
            self.output_file_.write("\t\th[%(n)i]= %(c)s;" % {
                'n': counter,
                'c': printing.ccode(ex.evalf())
            })
            if (counter % self.control_dim_):
                self.output_file_.write("// ---- end of row\n")
            else:
                self.output_file_.write("\n")

            counter = counter + 1
        self.output_file_.write("\t}\n")

        self.output_file_.write("""
        void dynSysLinearPart_1G(double t,const double *x,double dhdx[]){\n""")
        counter = 0
        for x in self.state_symbols_:
            for ex in self.hs:
                self.output_file_.write(
                    "\t\tdhdx[%(n)i]= %(c)s;" % {
                        'n': counter,
                        'c': printing.ccode(diff(ex, x).evalf())
                    })
                if (counter % self.adim):
                    self.output_file_.write("// ----- end of row\n")
                else:
                    self.output_file_.write("\n")
                counter = counter + 1
            self.output_file_.write("//--end of matrix\n\n")
        self.output_file_.write("\t}\n")

        self.output_file_.write("""
        void dynSysLinearPart_dt(double t,const double *x,double dhdt[]){\n""")
        counter = 0
        for ex in self.hs:
            self.output_file_.write(
                "\t\tdhdt[%(n)i]= %(c)s;\n" % {
                    'n': counter,
                    'c': printing.ccode(diff(ex, self.time_symbol_).evalf())
                })
            counter = counter + 1
        self.output_file_.write("\t}\n")

        self.output_file_.write("""
        double runningCostAffinePart(double t,const double *x,const double u2[]){
               return  %s;
        }\n""" % printing.ccode(self.g0.evalf()))

        self.output_file_.write("""
    void runningCostLinearPart(double t,const double x[],double h0[]){\n""")
        counter = 0
        for ex in self.h0s:
            self.output_file_.write("\t\th0[%(n)i]= %(c)s;\n" % {
                'n': counter,
                'c': printing.ccode(ex.evalf())
            })
            counter = counter + 1
        self.output_file_.write("\t}\n")

        self.output_file_.write("""
    void runningCostLinearPart_1G(double t,const double x[],double dh0dx[]){\n
""")
        counter = 0
        for x in self.state_symbols_:
            for ex in self.h0s:
                self.output_file_.write(
                    "\t\tdh0dx[%(n)i]= %(c)s;\n" % {
                        'n': counter,
                        'c': printing.ccode(diff(ex, x).evalf())
                    })
                counter = counter + 1
        self.output_file_.write("\t}\n")

        self.output_file_.write("""
    void runningCostLinearPart_dt(double t,const double *x,double h0t[]){\n""")
        counter = 0
        for ex in self.h0s:
            self.output_file_.write(
                "\t\th0t[%(n)i]= %(c)s;\n" % {
                    'n': counter,
                    'c': printing.ccode(diff(ex, self.time_symbol_).evalf())
                })
            counter = counter + 1
        self.output_file_.write("\t}\n")

    def create(self):
        self.writeClassHeader()
        self.writeSolverFuncs()
        self.write_max_hamiltonian()
        self.writeAuxFuncs()
        self.writeAffineFuncs()
        self.writeClassClosure()

    def write_max_hamiltonian(self):
        if self.control_dim_ == 1:
            self.output_file_.write("""
    void maxHamiltonian(double t,const double x[],
        const double lambda[],double u[]){

            compute_1Control(t,x,umin,umax,u,1.0e-6);
    }\n""")
        else:
            self.write_void_max_hamiltonian()

    def writeClassHeader(self):
        self.output_file_.write("""
#ifndef """ + self.class_name_.upper() + """_H
#define """ + self.class_name_.upper() + """_H

#include<affinesys.h>
#include<math.h>
class c""" + self.class_name_ + """:public affinesys{
    public:
    double x0[%(dim)i];
    double u0[%(adim)i];
    double umin[%(adim)i];
    double umax[%(adim)i];\n""" % {
            'dim': self.state_dim_,
            'adim': self.adim
        })

        for i in self.aux_vars_symbols_:
            self.output_file_.write("\tdouble %s;\n" % i.name)
        self.output_file_.write(
            """
    c%(cn)s():\n\taffinesys(%(dim)i,%(cdim)i,%(adim)i,1.0e-3,20.0,1.0e-7,x0)"""
            % {
                'cn': self.class_name_,
                'dim': self.state_dim_,
                'cdim': self.control_dim_,
                'adim': self.adim
            })
        for i, ex in zip(self.aux_vars_symbols_, self.aux_vars_values_):
            self.output_file_.write(",\n\t\t{:s} ( {:s} ) ".format(
                i.name, printing.ccode(ex)))

        self.output_file_.write("\n\t{\n")
        for i in range(self.state_dim_):
            self.output_file_.write("\t\tx0[%i]= 0.0;\n" % i)
        for i in range(self.control_dim_):
            self.output_file_.write("\t\tu0[{:d}]= {:s};\n".format(
                i, printing.ccode(self.control0_values_[i])))
            self.output_file_.write("\t\tumax[{:d}]= {:s};\n".format(
                i, printing.ccode(self.control_max_[i])))
            self.output_file_.write("\t\tumin[{:d}]= {:s};\n".format(
                i, printing.ccode(self.control_min_[i])))
        self.output_file_.write("\n\t}\n")

reason: from sympy


text:
from selection10 import genInterval
import pandas as pd
import sys
import numpy as np
import matplotlib.pyplot as plt
import os

genInt = genInterval(os.getcwd() + '/')
genInt.plotGenInt()

print 'Created plot GenInt_PLOT.png'