text:
#!/usr/bin/env python2.7

'''A tool for determining the UVIT FUV and NUV filters.

   Copyright 2017 Prajwel Joseph

   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.

   Changes; when, what
   -------------------
   Dec 22, 2017: bug fixes.
   Dec 23, 2017: deals with cases where FUV is not present.

   The author would like to acknowledge inputs from Dr. Koshy George
   which greatly helped the developement of this script.

'''

import os
import sys
import string
import urllib
import numpy as np

from astropy.io import fits
from requests import Session
from bs4 import BeautifulSoup
from astropy import units as u
from astropy.table import Table, hstack
from astropy.coordinates import SkyCoord

# To get the user input.
instrument = str(sys.argv[1])
RA = str(sys.argv[2])
DEC = str(sys.argv[3])

#instrument = 'uvit'
#RA = "7:36:51.396"
#DEC = "65:36:9.170"

# To check the RA, DEC user input.
if (DEC.count(':') == RA.count(':') == 2):
     pass
else:
   sys.exit(1)

# To do all the stuff in a specific directory.
working_arena = '.'
os.chdir(working_arena)

# instrument and radius of search in arsec.
field_radius = {'uvit'  : 1200,
                'sxt'   : 1500,
                'czti'  : 1680,
                'laxpc' : 1680}

# Functions to convert magnitude to UVIT count rates.
def countfuv(mg):
    caf2 = 1.0
    baf2 = 0.85
    sapphire = 0.63
    silica = 0.22
    mg1 = 18.22

    if mg <= 15.0:
        mg_c = 5.371 + (20.0 * mg - 210.2) ** 0.5
    else:
        mg_c = mg

    cr1 =  caf2 * 10.0 ** ((mg1 - mg_c) * 0.4)
    cr2 =  baf2 * 10.0 ** ((mg1 - mg_c) * 0.4)
    cr3 =  sapphire *10.0 ** ((mg1 - mg_c) * 0.4)
    cr4 =  silica * 10.0 ** ((mg1 - mg_c) * 0.4)
    return mg, mg_c, cr1, cr2, cr3, cr4

countfuv = np.vectorize(countfuv)

def countfuv_abs(mg): # for cases where FUV is absent.
    caf2 = 1.0
    baf2 = 0.85
    sapphire = 0.63
    silica = 0.22
    mg1 = 18.22

    if mg <= 15.0:
        mg_c = 2.634 + (26.316 * mg - 245.329) ** 0.5
    else:
        mg_c = mg

    mg_c = mg_c - 1.65
    cr1 =  caf2 * 10.0 ** ((mg1 - mg_c) * 0.4)
    cr2 =  baf2 * 10.0 ** ((mg1 - mg_c) * 0.4)
    cr3 =  sapphire *10.0 ** ((mg1 - mg_c) * 0.4)
    cr4 =  silica * 10.0 ** ((mg1 - mg_c) * 0.4)
    return mg, mg_c, cr1, cr2, cr3, cr4

countfuv_abs = np.vectorize(countfuv_abs)

def countnuv(mg):
    silica = 1.0
    b4 = 0.22
    b13 = 0.27
    b15 = 0.074
    n2 = 0.055
    mg1 = 20.0
    if mg <= 15.0:
        mg_c = 2.634 + (26.316 * mg - 245.329) ** 0.5
    else:
        mg_c = mg

    cr1 =  silica * 10.0 ** ((mg1 - mg_c) * 0.4)
    cr2 =  b4 * 10 ** ((mg1-mg_c) * 0.4)
    cr3 =  b13 * 10 ** ((mg1 - mg_c) * 0.4)
    cr4 =  b15 * 10 ** ((mg1 - mg_c) * 0.4)
    cr5 =  n2 * 10 ** ((mg1 - mg_c) * 0.4)
    return mg, mg_c, cr1, cr2, cr3, cr4, cr5

countnuv = np.vectorize(countnuv)

# Function to find seperation in celestial coordinates.
cc = SkyCoord(RA, DEC, unit = (u.hourangle, u.deg))
def cel_separation(a, b):
    coo = SkyCoord(a, b, frame = 'icrs', unit = 'deg')
    return coo.separation(cc)

# To get on with MAST website queries.
ra = string.replace(RA, ':', '+')
dec = string.replace(DEC, ':', '+')

# Mast website form data.
mastdata = {'__EVENTTARGET': '""',
            '__EVENTARGUMENT' : '""',
            '__VIEWSTATE' : '[encoded data, 13600 characters]',
           '__VIEWSTATEGENERATOR': 'C84C2718',
           '__EVENTVALIDATION': '/wEdAAue+6xrb6xgp2ityzurA/pfWsTF2CBs9ziYHlDmus7EnHXVqisK/ch+FuYDN4RJj9bNygAwoalISibjyjYgoB7/Pb1PMsXU2LG7o+i6/zoft2ZmqVWZEJyWTGlJer/5/ymk9SeG9Y8RLbkbyiuf4BcRXP2SoyGCMZyu6LfyUjL5ZgAB13huDNxtBirRDFLR6zW3raPnQUy5sK21W/3eiEs/KUQOVtp9GallVy/IsFMIp4yMEruOYx0KrV7GUndYi0m5y40+',
           '_ctl10:txtTargetName': '',
           '_ctl10:resolverDropList': 'SIMBAD',
           '_ctl10:txtRadius': '0.4',
           '_ctl10:txtRA': ra,
           '_ctl10:txtDec': dec,
           '_ctl10:btnSearch': 'Search'}

# Header for the site.
header = {'Host': 'galex.stsci.edu',
          'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64; rv:45.0) Gecko/20100101 Firefox/45.0',
          'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8',
          'Accept-Language': 'en-US,en;q=0.5',
          'Accept-Encoding': 'gzip, deflate',
          'Referer': 'http://galex.stsci.edu/GR6/?page=tilelist&survey=allsurveys',
          'Connection': 'keep-alive'}

# Post request.
session = Session()
response = session.post(url = 'http://galex.stsci.edu/GR6/?page=tilelist&survey=allsurveys',
                        data = mastdata,
                        headers = header)

# To make sense of the mess that is MAST.
soup = BeautifulSoup(response.text, 'html.parser')
ss = soup.find(id = '_ctl10_TileGrid_imgLink_0').get('href')
parent = 'http://galex.stsci.edu/GR6'
son = parent + ss[1:]

# To make sense of the download page.
response2 = session.get(son)
soup2 = BeautifulSoup(response2.text, 'html.parser')
pp = soup2.find_all("a")

# To download the mcat (galex catalogue).
catalogue_link = []
for link in pp:
    somel = link.get('href')
    try:
        if somel[-12:] == 'mcat.fits.gz':
            catalogue_link.append(somel)
    except:
        pass

if len(catalogue_link) != 0:
    catalogue = catalogue_link[0].split('/')[-1]
    urllib.urlretrieve(catalogue_link[0], catalogue)
else:
    sys.exit(1)

# Reading coordinates from catalogue.
hdu = fits.open(catalogue)
alpha = hdu[1].data['alpha_j2000_merged']
delta = hdu[1].data['delta_j2000_merged']

# NUV
nuv_mag = hdu[1].data['nuv_mag']
refined_set = [(al, de, nm) for al, de, nm
                            in zip(alpha, delta, nuv_mag)
                            if int(nm) != -999 and nm <= 22.]

nalpha, ndelta, nuv_mag = zip(*refined_set)

confined_set = [(nm, al, de) for al, de, nm
                             in zip(nalpha, ndelta, nuv_mag)
                             if cel_separation(al, de) <= field_radius[instrument] * u.arcsec]

nd = np.array(sorted(confined_set))[0:5]
ma, ma_c, ta, tb, tc, td, te = countnuv(nd[:,0])
nuv_res = Table([ma, ma_c, ta, tb, tc, td, te],
               names = ('Mag',
                        'Mag_corrected',
                        'silica',
                        'b4',
                        'b13',
                        'b15',
                        'n2'),
               meta = {'name': 'NUV counts'})

nuv_res['Mag'].format = '4.2f'
nuv_res['Mag_corrected'].format = '4.2f'
nuv_res['silica'].format = '4.2f'
nuv_res['b4'].format = '4.2f'
nuv_res['b13'].format = '4.2f'
nuv_res['b15'].format = '4.2f'
nuv_res['n2'].format = '4.2f'

# To convert ra_deg and dec_deg to ra_hms and dec_dms.
coord = SkyCoord(zip(nd[:,1], nd[:,2]),
                 frame = 'icrs',
                 unit = 'deg')

RAhms_DECdms = coord.to_string('hmsdms', sep = ':')
ra_hms, dec_dms = zip(*[hmdm.split(' ') for hmdm in RAhms_DECdms])
xy_tab = Table([ra_hms, dec_dms], names = ('ra_hms', 'dec_dms'))
nuv_res = hstack([xy_tab, nuv_res])

print('\n\n### NUV\n\n{}\n'.format(nuv_res))

# FUV
fuv_mag = hdu[1].data['fuv_mag']

fuv_absent = 'no'
if len(np.unique(fuv_mag)) == 1:  # when FUV data is absent.
    fd = nd
    print('\nFUV observations seem to be absent! Using M_fuv = M_nuv - 1.65.')
    fuv_absent = 'yes'
else:
    refined_set = [(al, de, fm) for al, de, fm
                                in zip(alpha, delta, fuv_mag)
                                if int(fm) != -999 and fm <= 22.]

    falpha, fdelta, fuv_mag = zip(*refined_set)

    confined_set = [(fm, al, de) for al, de, fm
                                 in zip(falpha, fdelta, fuv_mag)
                                 if cel_separation(al, de) <= field_radius[instrument] * u.arcsec]

    fd = np.array(sorted(confined_set))[0:5]

if fuv_absent == 'no':
    ma, ma_c, ta, tb, tc, td = countfuv(fd[:,0])
else:
    ma, ma_c, ta, tb, tc, td = countfuv_abs(fd[:,0])

fuv_res = Table([ma, ma_c, ta, tb, tc, td],
               names = ('Mag',
                        'Mag_corrected',
                        'caf2',
                        'baf2',
                        'sapphire',
                        'silica'),
               meta = {'name': 'FUV counts'})

fuv_res['Mag'].format = '4.2f'
fuv_res['Mag_corrected'].format = '4.2f'
fuv_res['caf2'].format = '4.2f'
fuv_res['baf2'].format = '4.2f'
fuv_res['sapphire'].format = '4.2f'
fuv_res['silica'].format = '4.2f'

# To convert ra_deg and dec_deg to ra_hms and dec_dms.
coord = SkyCoord(zip(fd[:,1], fd[:,2]),
                 frame = 'icrs',
                 unit = 'deg')

RAhms_DECdms = coord.to_string('hmsdms', sep = ':')
ra_hms, dec_dms = zip(*[hmdm.split(' ') for hmdm in RAhms_DECdms])
xy_tab = Table([ra_hms, dec_dms], names = ('ra_hms', 'dec_dms'))
fuv_res = hstack([xy_tab, fuv_res])

print('\n### FUV \n\n{}\n\nDone!\n'.format(fuv_res))

max_stars_repo_path: Gaia (UV filter checking tool)/fast_mcat_uvcheck_V.1.1.py

max_stars_repo_name: prajwel/UVIT-POC

max_issues_repo_name: prajwel/UVIT-POC

max_forks_repo_name: prajwel/UVIT-POC


text:
# -*- coding: utf-8 -*-
"""
Routines for multivariate pressure prediction

Created on Sep 20 2018
"""
from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

__author__ = "yuhao"

import numpy as np

def multivariate_virgin(sigma, phi, vsh, a_0, a_1, a_2, a_3, B):
    """
    Calculate velocity using multivariate virgin curve

    Parameters
    ----------
    sigma : 1-d ndarray
        effective pressure
    phi : 1-d ndarray
        effective porosity
    vsh : 1-d ndarray
        shale volume
    a_0, a_1, a_2, a_3 : float
        coefficients of equation
    B : float
        effective pressure exponential

    Returns
    -------
    out : 1-d ndarray
        velocity array

    Notes
    -----
    .. math:: V = a_0 + a_1\\phi + a_2{V}_{sh} + a_3 {\\sigma}^{B}

    [5]_

    .. [5] Sayers, C., Smit, T., van Eden, C., Wervelman, R., Bachmann, B.,
       Fitts, T., et al. (2003). Use of reflection tomography to predict
       pore pressure in overpressured reservoir sands. In submitted for
       presentation at the SEG 2003 annual meeting.
    """
    return a_0 - a_1 * phi - a_2 * vsh + a_3 * sigma**B

def invert_multivariate_virgin(vel, phi, vsh, a_0, a_1, a_2, a_3, B):
    """
    Calculate effective stress using multivariate virgin curve

    Parameters
    ----------
    vel : 1-d ndarray
        velocity array whose unit is m/s.
    phi : 1-d ndarray
        porosity array
    vsh : 1-d ndarray
        shale volume
    a_0, a_1, a_2, a_3 : scalar
        coefficients

    Returns
    -------
    sigma: 1-d ndarray
    """
    return ((vel - a_0 + a_1 * phi + a_2 * vsh) / a_3)**(1 / B)

def multivariate_unloading(sigma, phi, vsh, a_0, a_1, a_2, a_3, B, U, vmax):
    """
    Calculate velocity using multivariate unloading curve
    """
    sigma_max = invert_multivariate_virgin(
        vmax, phi, vsh, a_0, a_1, a_2, a_3, B)

    sigma_vc = sigma_max*(sigma/sigma_max)**(1/U)

    return multivariate_virgin(sigma_vc, phi, vsh, a_0, a_1, a_2, a_3, B)

def invert_multivariate_unloading(vel, phi, vsh, a_0, a_1, a_2, a_3,
                                  B, U, vmax):
    """
    Calculate effective stress using multivariate unloading curve
    """
    sigma_max = invert_multivariate_virgin(
        vmax, phi, vsh, a_0, a_1, a_2, a_3, B)

    sigma_vc = invert_multivariate_virgin(
        vel, phi, vsh, a_0, a_1, a_2, a_3, B)

    return sigma_max * (sigma_vc/sigma_max)**U

def effective_stress_multivariate(vel, phi, vsh, a_0, a_1, a_2, a_3,
                                  B, U, vmax, start_idx, end_idx=None):
    ves = invert_multivariate_virgin(vel, phi, vsh, a_0, a_1, a_2, a_3, B)
    unloading = invert_multivariate_unloading(
        vel, phi, vsh, a_0, a_1, a_2, a_3, B, U, vmax)
    ves[start_idx: end_idx] = unloading[start_idx: end_idx]
    return ves

def pressure_multivariate(obp, vel, phi, vsh, a_0, a_1, a_2, a_3,
                          B, U, vmax, start_idx, end_idx=None):
    """
    Pressure Prediction using multivariate model
    """
    ves = effective_stress_multivariate(
        vel, phi, vsh, a_0, a_1, a_2, a_3, B, U, vmax, start_idx, end_idx)
    return obp - ves

def pressure_multivariate_varu(obp, vel, phi, vsh, a_0, a_1, a_2, a_3,
                               B, U, vmax, start_idx, buf=20,
                               end_idx=None, end_buffer=10):
    """
    Pressure Prediction using multivariate model
    """
    ves = effective_stress_multivariate_varu(
        vel, phi, vsh, a_0, a_1, a_2, a_3,
        B, U, vmax, start_idx, buf, end_idx, end_buffer)
    return obp - ves

def effective_stress_multivariate_varu(vel, phi, vsh, a_0, a_1, a_2, a_3,
                                       B, U, vmax, start_idx, buf=20,
                                       end_idx=None, end_buffer=10):
    u_array = np.ones(vel.shape)
    u_array[start_idx: end_idx] = U
    # start buffer
    u_buffer = np.linspace(1, U, buf)
    u_array[start_idx-buf+1: start_idx + 1] = u_buffer
    # end buffer
    if end_idx is not None:
        u_array[end_idx: end_idx + end_buffer] = np.linspace(U, 1, end_buffer)
    ves = invert_multivariate_unloading(vel, phi, vsh, a_0, a_1, a_2, a_3,
                                        B, u_array, vmax)
    return ves

ext: py

reason: import numpy


text:
import numpy as SP
import subprocess, sys, os.path

class SnpsFromFile(object): # implements ISnpSet
    '''
     See the Bed class's 'read' method of examples of its use.
     See __init__.py for specification of interface it implements.
     '''
    def __init__(self, filename):
        self.filename = filename

    def addbed(self, bed):
        return SnpsFromFilePlusBed(self.filename, bed)

class SnpsFromFilePlusBed(object): # implements ISnpSetPlusBed
    def __init__(self, filename, bed):
        self.filename = filename
        snplist0 = ReadTokens(filename)
        bed.run_once()
        self.snplist = [val for val in snplist0 if val in bed.snp_to_index] #use the intersection
        if len(snplist0) != len(self.snplist): raise Exception("The file '{0}' contains SNPs not in bed '{1}'".format(filename,bed.basefilename))
        self.bed = bed

    def __str__(self):
        return self.filename

    def __len__(self):
        return len(self.snplist)

    def __iter__(self):
        for snp in self.snplist:
            index = self.bed.snp_to_index[snp]
            yield index

    def read(self):
        return self.bed.read_with_specification(self)

    @property
    def pos(self):
        """
        Returns:
            pos:    position of the SNPs in the specification
        """
        return self.bed.pos[self.to_index]

    @property
    def to_index(self):
        iter = self.__iter__()
        return [i for i in iter]

def ReadTokens(filename):
    token_list = open(filename).read().split() # could be written to 'yield' the tokens one at a time
    return token_list;

max_stars_repo_name: HealthML/FaST-LMM

max_issues_repo_name: epiproject/FaST-LMM-HPC

max_forks_repo_name: epiproject/FaST-LMM-HPC


text:
r'''
Model:

    Z_row,i ~ N(B_row C_{row,i}, Sigma_row)
    Z_col,j ~ N(B_col C_{col,i}, Sigma_col)
    C_ij = <Z_row,i,Z_col,j>
    X_ij ~ p(C_ij;theta_j)

where p can be one of
- bernoulli
- normal
- negativebinomial

'''

import dataclasses
from . import lowlevel
import tensorflow as tf
import numpy.random as npr
import numpy as np
import scipy as sp
import scipy.sparse.linalg
from . import helpers
from . import sparsematrix

import logging
logger = logging.getLogger(__name__)

class NoneArray:
    def __getitem__(self,p):
        return self

@dataclasses.dataclass
class PosteriorGaussian:
    muhat: tf.Variable
    Sighat: tf.Variable
    lr: float=1.0

    @property
    def dtype(self):
        return self.muhat.dtype

    @property
    def double_precision(self):
        return self.dtype==tf.float64

    def copy(self):
        return PosteriorGaussian(
            tf.Variable(self.muhat),
            tf.Variable(self.Sighat),
            self.lr
        )

    def validate(self):
        with tf.device('cpu'):
            return dict(
                symmetric = tf.reduce_all(self.Sighat == helpers.batch_transpose(self.Sighat)).numpy(),
                # pd = (tf.reduce_min(tf.linalg.eigvalsh(self.Sighat)[0]))==1).numpy()
                pd = tf.reduce_min(tf.linalg.eigvalsh(self.Sighat)).numpy()
            )

    @property
    def Nk(self):
        return self.muhat.shape[1]

    def improve(self,omega1,omega2,omega1_prior,omega2_prior):
        lowlevel.assign_stats2natural(omega1,omega2,self.muhat,self.Sighat,self.lr,omega1_prior,omega2_prior)

    def snapshot(self):
        return dict(muhat=self.muhat.numpy(),Sighat=self.Sighat.numpy())

    def sample(self):
        noise=tf.random.normal(self.muhat.shape,dtype=self.muhat.dtype)
        return tf.einsum('ikl,il->ik',tf.linalg.cholesky(self.Sighat),noise)+self.muhat

@dataclasses.dataclass
class PriorGaussian:
    B: tf.Variable  # Nk x Ncov
    C: tf.Variable  # Nsamp x Ncov
    Sig: tf.Variable
    sigmask: tf.Variable = None

    def __post_init__(self):
        if self.sigmask is None:
            self.sigmask=tf.Variable(tf.ones(self.Sig.shape,dtype=self.dtype))

    @property
    def dtype(self):
        return self.B.dtype

    @property
    def double_precision(self):
        return self.dtype==tf.float64

    def copy(self):
        return PriorGaussian(
            tf.Variable(self.B),
            tf.Variable(self.C),
            tf.Variable(self.Sig),
            tf.Variable(self.sigmask),
        )

    def validate(self):
        return dict(
            symmetric = tf.reduce_all(self.Sig == helpers.batch_transpose(self.Sig)).numpy(),
            # pd = (tf.linalg.slogdet(self.Sig)[0]==1).numpy()
            pd = tf.reduce_min(tf.linalg.eigvalsh(self.Sig)).numpy()
        )

    def snapshot(self):
        return dict(B=self.B.numpy(),C=self.C.numpy(),Sig=self.Sig.numpy(),sigmask=self.sigmask.numpy())

    @property
    def mu(self):
        return self.C@tf.transpose(self.B)

@dataclasses.dataclass
class GaussianMatrixVI:
    post: PosteriorGaussian
    prior: PriorGaussian

    def copy(self):
        return GaussianMatrixVI(self.post.copy(),self.prior.copy())

    def validate(self):
        return dict(
            post=self.post.validate(),
            prior=self.prior.validate()
        )

    def update_post_from_suffstats(self,omega1,omega2):
        omega2_prior=tf.linalg.inv(self.prior.Sig)
        omega1_prior = tf.einsum('kl,ik->il',omega2_prior,self.prior.mu)
        self.post.improve(omega1,omega2,omega1_prior,omega2_prior[None,:,:])

    def mult(self,k):
        k2=k[:,None]*k[None,:]
        self.prior.Sig.assign(helpers.symmetrize(self.prior.Sig*k2))
        self.prior.B.assign(self.prior.B*k[:,None])
        self.post.muhat.assign(self.post.muhat*k[None,:])
        self.post.Sighat.assign(helpers.symmetrize(self.post.Sighat*k2))

    def reinitialize_covariates(self,C):
        self.prior.C=tf.Variable(tf.convert_to_tensor(C,dtype=self.dtype))
        self.prior.B=tf.Variable(tf.zeros((self.Nk,self.prior.C.shape[1]),dtype=self.dtype))
        self.update_prior()

    def update_prior(self,numlo=1e-8):
        # get best B
        self.prior.B.assign(tf.transpose(tf.linalg.lstsq(self.prior.C,self.post.muhat)))

        # get sig
        df = self.post.muhat - self.prior.mu
        newsig = tf.reduce_mean(tf.einsum('ij,ik->ijk',df,df)+self.post.Sighat,axis=0)
        newsig = helpers.clip_low_eigenvalues(newsig,numlo)

        # zero out some bits of the sigma
        if self.prior.sigmask is not None:
            newsig=newsig * self.prior.sigmask

        self.prior.Sig.assign(newsig)

    def kl(self):
        klval,badness=lowlevel.prior_KL_ind(self.post.muhat,self.post.Sighat,self.prior.mu,self.prior.Sig)
        return klval

    def kl_subset(self,indices):
        muhat=tf.gather(self.post.muhat,indices)
        Sighat=tf.gather(self.post.Sighat,indices)
        mu=tf.gather(self.prior.mu,indices)
        klval,badness=lowlevel.prior_KL_ind(muhat,Sighat,mu,self.prior.Sig)
        return klval

    def snapshot(self):
        return dict(post=self.post.snapshot(),prior=self.prior.snapshot())

    @property
    def shape(self):
        return self.post.muhat.shape

    def shapecheck(self):
        nsamp,Nk = self.shape
        ncov= self.prior.C.shape[1]
        assert self.post.Sighat.shape==(nsamp,Nk,Nk)
        assert self.prior.C.shape==(nsamp,ncov)
        assert self.prior.B.shape==(Nk,ncov)
        assert self.prior.Sig.shape==(Nk,Nk)
        assert self.prior.sigmask.shape==(Nk,Nk)

    @property
    def Nk(self):
        return self.post.muhat.shape[1]

    @classmethod
    def load(cls,snap,dtype=tf.float64):
        return GaussianMatrixVI(
            PosteriorGaussian(**helpers.tf_variabilify(snap['post'],dtype)),
            PriorGaussian(**helpers.tf_variabilify(snap['prior'],dtype)),
        )

    @classmethod
    def load_from_muhat_and_C(cls,muhat,C,dtype=tf.float64,numlo=1e-8,
                        Sighat_ratio=.001,Sighat_adder=.001,
                        sigmask=None,diagsig=None):
        if sigmask is not None:
            assert diagsig is None
            sigmask=tf.Variable(tf.convert_to_tensor(sigmask,dtye=dtype))
        elif diagsig:
            sigmask=tf.Variable(tf.convert_to_tensor(np.eye(muhat.shape[1]),dtype=dtype))

        Sighat=(np.diag(np.var(muhat,axis=0)+Sighat_adder)[None,:,:])*np.ones(len(muhat))[:,None,None]*Sighat_ratio
        post=PosteriorGaussian(**helpers.tf_variabilify(dict(muhat=muhat,Sighat=Sighat),dtype))
        prior=PriorGaussian(
            C=tf.Variable(tf.convert_to_tensor(C,dtype=dtype)),
            B=tf.Variable(tf.zeros((post.muhat.shape[1],C.shape[1]),dtype=dtype)),
            Sig=tf.Variable(tf.zeros((post.muhat.shape[1],post.muhat.shape[1]),dtype=dtype)),
            sigmask=sigmask
        )
        gmvi=GaussianMatrixVI(post,prior)
        gmvi.shapecheck()
        gmvi.update_prior()
        return gmvi

    @property
    def dtype(self):
        return self.post.dtype

    @property
    def double_precision(self):
        return self.dtype==tf.float64

@dataclasses.dataclass
class Model:
    rowinfo: GaussianMatrixVI
    colinfo: GaussianMatrixVI
    kind: str # what is the datamodel like
    thetas: tf.Tensor # a value for each column
    row_batchsize: int = 100
    col_batchsize: int = 100

    def copy(self):
        return Model(
            self.rowinfo.copy(),
            self.colinfo.copy(),
            self.kind,
            self.thetas,
            self.row_batchsize,
            self.col_batchsize)

    @property
    def double_precision(self):
        return self.rowinfo.double_precision

    @property
    def dtype(self):
        return self.rowinfo.dtype

    def snapshot(self):
        return dict(
            rowinfo=self.rowinfo.snapshot(),
            colinfo=self.colinfo.snapshot(),
            thetas=self._snap_thetas(),
            kind=self.kind,
            row_batchsize=self.row_batchsize,
            col_batchsize=self.col_batchsize,
            double_precision=self.double_precision
        )

    def __hash__(self):
        return id(self)

    def set_lr(self,row_lr=None,col_lr=None):
        if row_lr is not None:
            self.rowinfo.post.lr=row_lr
        if col_lr is not None:
            self.colinfo.post.lr=col_lr

    @classmethod
    def _load_thetas(cls,thetas,dtype):
        return helpers.tf_variabilify_one(thetas,dtype=dtype)

    def _snap_thetas(self):
        if self.thetas is None:
            return None
        else:
            return self.thetas.numpy()

    @classmethod
    def load(cls,snap,double_precision=None):
        kwargs={}
        if 'row_batchsize' in snap:
            kwargs['row_batchsize']=snap['row_batchsize']
        if 'col_batchsize' in snap:
            kwargs['col_batchsize']=snap['col_batchsize']

        # dtype
        if double_precision is None:
            double_precision=snap['double_precision']

        if double_precision:
            dtype=tf.float64
        else:
            dtype=tf.float32

        rez=cls(
            GaussianMatrixVI.load(snap['rowinfo'],dtype=dtype),
            GaussianMatrixVI.load(snap['colinfo'],dtype=dtype),
            snap['kind'],
            cls._load_thetas(snap['thetas'],dtype),
            **kwargs
        )

        rez.rowinfo.validate()
        rez.colinfo.validate()

        return rez

    @property
    def shape(self):
        return self.rowinfo.shape[0],self.colinfo.shape[0]

    @property
    def Nk(self):
        return self.rowinfo.Nk

    @property
    def row_loadings(self):
        return self.rowinfo.post.muhat.numpy()
    @property
    def col_loadings(self):
        return self.colinfo.post.muhat.numpy()

    @property
    def nobs(self):
        shape=self.shape
        return shape[0]*shape[1]

    ############################################################
    # STAT COLLECTION

    def _mom_row(self,st=0,en=None):
        return lowlevel.get_moments(self.rowinfo.post.muhat[st:en],self.rowinfo.post.Sighat[st:en])

    def _mom_col(self,st=0,en=None):
        return lowlevel.get_moments(self.colinfo.post.muhat[st:en],self.colinfo.post.Sighat[st:en])

    def _thetas_for_rows(self):
        if lowlevel.has_theta[self.kind]:
            return self.thetas[None,:]
        else:
            return None

    def _thetas_for_cols(self,st=0,en=None):
        if lowlevel.has_theta[self.kind]:
            return self.thetas[st:en,None]
        else:
            return None

    def _get_row_data_stats(self,X):
        # collect moments about the columns
        mom_col = self._mom_col()

        # accumulate omegas, in a batched way
        def go(st,en):
            mom_row = self._mom_row(st,en)
            o1,o2=lowlevel.accumulate_omega_for_rows(
                X[st:en],mom_row,mom_col,self.kind,self._thetas_for_rows())
            return o1,o2
        return helpers.map_accum(go,self.shape[0],self.row_batchsize,'cc')

    def _get_col_data_stats(self,X):
        # collect moments about the columns
        mom_row = self._mom_row()

        # accumulate omegas, in a batched way
        def go(st,en):
            mom_col = self._mom_col(st,en)
            o1,o2=lowlevel.accumulate_omega_for_rows(
                X.T[st:en],mom_col,mom_row,self.kind,self._thetas_for_cols(st,en))
            return o1,o2
        return helpers.map_accum(go,self.shape[1],self.col_batchsize,'cc')

    def _get_zeta(self,X):
        # collect some info
        if not lowlevel.has_theta[self.kind]:
            return None
        else:
            mom_row =self._mom_row()

            def go(st,en):
                mom_col = self._mom_col(st,en)
                mn,vr=lowlevel.mom2mnvr(mom_col,mom_row)
                zeta = lowlevel.get_zeta(X.T[st:en],mn,vr,self.kind,self._thetas_for_cols(st,en))
                return tf.reduce_sum(zeta,axis=1),
            return helpers.map_accum(go,self.shape[1],self.col_batchsize,'c')[0]

    ############################################################
    # UPDATES
    things_to_update=['rows','cols','prior_rows','prior_cols','thetas','rebalance']

    def update_rebalance(self,X):
        # A=np.diag(self.rowinfo.prior.Sig.numpy())
        # B=np.diag(self.colinfo.prior.Sig.numpy())
        A=tf.reduce_mean(self.rowinfo.post.muhat**2,axis=0)
        B=tf.reduce_mean(self.colinfo.post.muhat**2,axis=0)

        k=tf.convert_to_tensor((A/B)**(.25),dtype=self.dtype)
        self.rowinfo.mult(1/k)
        self.colinfo.mult(k)

    def update_rows(self,X):
        o1r,o2r = self._get_row_data_stats(X)
        self.rowinfo.update_post_from_suffstats(o1r,o2r)

    def update_cols(self,X):
        o1c,o2c = self._get_col_data_stats(X)
        self.colinfo.update_post_from_suffstats(o1c,o2c)

    def update_thetas(self,X,minth=1e-12):
        if not lowlevel.has_theta[self.kind]:
            return
        zetamean = self._get_zeta(X) / self.shape[0]
        newthetas=lowlevel.solve_theta(zetamean,self.kind)
        newthetas=tf.clip_by_value(newthetas,minth,np.inf)
        self.thetas.assign(newthetas)

    def update_prior_rows(self,X):
        self.rowinfo.update_prior()

    def update_prior_cols(self,X):
        self.colinfo.update_prior()

    ###############################################################
    # LOSS

    def _dataterm(self,X):
        mom_col=self._mom_col()
        dataterm=0
        for (st,en) in helpers.batchup(self.shape[0],self.row_batchsize):
            logger.debug(f'dataterm: {st}/{self.shape[0]}')
            mom_row=self._mom_row(st,en)
            mn,vr=lowlevel.mom2mnvr(mom_row,mom_col)
            dataterm=dataterm+tf.reduce_sum(lowlevel.ELBO_dataterm(X[st:en],mn,vr,self.kind,self.thetas))
        return dataterm

    def _loss(self,X):
        dataterm = self._dataterm(X)
        kl_row = self.rowinfo.kl()
        kl_col = self.colinfo.kl()

        nats = -(dataterm -kl_row - kl_col) / (self.nobs)
        rez=dict(loss=nats,kl_row=kl_row,kl_col=kl_col,dataterm=dataterm)
        return rez

    def loss(self,X,calc_raw_nats=True):
        rez = self._loss(X)
        return {x:rez[x].numpy() for x in rez}

    ###############################################################
    # SAMPS

    def posterior_samples(self):
        return self.rowinfo.post.sample().numpy(),self.colinfo.post.sample().numpy()

    def posterior_predictive_sample(self,memlimit=4,sample_loadings=True,sparse=None):
        if sparse is None:
            sparse = lowlevel.counts_data[self.kind]

        if not sparse:
            if 8*self.nobs > 1e9*memlimit:
                raise Exception(f"refusing to create matrix bigger than {memlimit}GB; increase keyword memlimit or use sparse=True")

        if sample_loadings:
            Zrow,Zcol = self.rowinfo.post.sample(),self.colinfo.post.sample()
        else:
            Zrow,Zcol=self.rowinfo.post.muhat,self.colinfo.post.muhat

        # accumulate samples, in a batched way
        if sparse:
            def go(st,en):
                C=Zrow[st:en] @ tf.transpose(Zcol)
                return lowlevel.sample(C,self.kind,self._thetas_for_rows())
            samp = helpers.map_cat_sparse(go,self.shape[0],self.row_batchsize) # <-- returns CSRMatrix
            return samp.to_scipy()
        else:
            def go(st,en):
                C=Zrow[st:en] @ tf.transpose(Zcol)
                return lowlevel.sample(C,self.kind,self._thetas_for_rows()),
            samp, = helpers.map_accum(go,self.shape[0],self.row_batchsize,'c') # <-- returns Tensor
            return samp.numpy()

r'''
 _       _ _   _       _ _          _   _
(_)_ __ (_) |_(_) __ _| (_)______ _| |_(_) ___  _ __
| | '_ \| | __| |/ _` | | |_  / _` | __| |/ _ \| '_ \
| | | | | | |_| | (_| | | |/ / (_| | |_| | (_) | | | |
|_|_| |_|_|\__|_|\__,_|_|_/___\__,_|\__|_|\___/|_| |_|

                                                      '''

def initialize_half(U,C=None,dtype=tf.float64,sigmultiplier=.1,diagsig=False):
    Sig = np.diag(np.var(U,axis=0))

    if C is None:
        C=tf.ones((U.shape[0],1),dtype=dtype)

    Nk=U.shape[1]

    if diagsig:
        sigmask= tf.Variable(tf.eye(Nk,dtype=dtype))
    else:
        sigmask = None

    gmvi=GaussianMatrixVI(
        post = PosteriorGaussian(
            muhat = tf.Variable(tf.convert_to_tensor(U,dtype=dtype)),
            Sighat = tf.Variable(tf.convert_to_tensor((sigmultiplier**2)*np.tile(Sig,[U.shape[0],1,1])*np.std(U),dtype=dtype))
        ),
        prior = PriorGaussian(
            C = tf.Variable(C),
            B = tf.Variable(tf.zeros((U.shape[1],C.shape[1]),dtype=dtype)),
            Sig = tf.Variable(tf.convert_to_tensor(np.cov(U.T).reshape((Nk,Nk)),dtype=dtype)),
            sigmask=sigmask,
        )
    )
    gmvi.update_prior()
    return gmvi

def make_dmhalf_op(data):
    if sp.sparse.issparse(data):
        ASL=sp.sparse.linalg.aslinearoperator
        assert sp.sparse.issparse(data)
        Nr,Nc = data.shape
        one_Nr=ASL(np.ones((Nr,1)))
        one_Nc=ASL(np.ones((Nc,1)))
        return 4*(ASL(data) - .5*(one_Nr@one_Nc.H))
    else:
        return 4*(data-.5)

def make_log_op(data):
    if sp.sparse.issparse(data):
        data=data.tocsr()
        data=sp.sparse.csr_matrix((np.log(1+data.data),data.indices,data.indptr))
        return data
    else:
        return np.log(1+data)

def var_data(data):
    if sp.sparse.issparse(data):
        c = data.copy()
        c.data **= 2
        m2 = np.asarray(np.mean(c,axis=0)).ravel()
        m1 = np.asarray(np.mean(data,axis=0)).ravel()
        return m2-m1**2
    else:
        return np.var(data,axis=0)

def initialize_with_new_rows(n_rows,snap,covariates=None):
    model=Model.load(snap)
    Nk=model.Nk

    if covariates is None:
        model.rowinfo.prior.C=tf.Variable(tf.ones((n_rows,1),dtype=model.dtype))
        model.rowinfo.prior.B=tf.Variable(tf.zeros((model.Nk,1),dtype=model.dtype))
    else:
        model.rowinfo.prior.C=tf.Variable(tf.convert_to_tensor(covariates,dtype=model.dtype))
        model.rowinfo.prior.B=tf.Variable(tf.zeros((model.Nk,covariates.shape[1]),dtype=model.dtype))

    model.rowinfo.post.Sighat=tf.Variable(tf.zeros((n_rows,Nk,Nk),dtype=model.dtype))
    model.rowinfo.post.muhat=tf.Variable(tf.zeros((n_rows,Nk),dtype=model.dtype))
    model.rowinfo.post.Sighat.assign(
        model.rowinfo.prior.Sig[None,:,:]*tf.ones(n_rows,dtype=model.dtype)[:,None,None])

    return model

def initialize(data,Nk,kind,dtype=tf.float64,diagsig=False,row_covariates=None,col_covariates=None,
                                minvar=1e-8):
    if sparsematrix.is_in_tensorflow_land(data):
        data=data._source.astype(np.float64)
    else:
        data=data.astype(np.float64)

    if kind=='normal':
        thetas=var_data(data)
        thetas[thetas<minvar]=minvar
        thetas=tf.convert_to_tensor(thetas,dtype=dtype)
        tot=data
    elif kind=='bernoulli':
        tot= make_dmhalf_op(data)
        thetas= tf.zeros(data.shape[1],dtype=dtype)
    elif kind=='negativebinomial':
        thetas=tf.ones(data.shape[1],dtype=dtype)
        tot = make_log_op(data)
    else:
        raise Exception("NYI")

    U,e,V=sp.sparse.linalg.svds(tot,Nk)
    V=V.T
    U=U@np.diag(np.sqrt(e))
    V=V@np.diag(np.sqrt(e))

    model= Model(
        initialize_half(U,dtype=dtype,diagsig=diagsig),
        initialize_half(V,dtype=dtype,diagsig=diagsig),
        thetas=tf.Variable(thetas),
        kind=kind
    )

    if row_covariates is not None:
        model.rowinfo.reinitialize_covariates(row_covariates)
    if col_covariates is not None:
        model.colinfo.reinitialize_covariates(col_covariates)

    return model

def example_model(N_rows,N_cols,Nk,kind,mag=1,sigmultiplier=.1,dtype=tf.float64):
    return Model(
        initialize_half(npr.randn(N_rows,Nk)*np.sqrt(mag),sigmultiplier=sigmultiplier,dtype=dtype),
        initialize_half(npr.randn(N_cols,Nk)*np.sqrt(mag),sigmultiplier=sigmultiplier,dtype=dtype),
        thetas=tf.Variable(tf.convert_to_tensor(1+npr.rand(N_cols)*.05,dtype=dtype)),
        kind=kind,
    )

max_stars_repo_name: jacksonloper/signedfac

max_issues_repo_name: jacksonloper/signedfac

max_forks_repo_name: jacksonloper/signedfac

reason: import numpy,import scipy


text:
# Copyright 2019 Google Inc. All Rights Reserved.
#
# 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.
# ==============================================================================

import importlib

from mtl.meta.optim import optimizer
import numpy as np

class MetaOptimizer(optimizer.MetaOptimizer):

  def __init__(self, opt):
    self.cfg = importlib.import_module('mtl.config.' + opt.config)
    self.num_exps = len(self.cfg.exp_list)
    self.exps_done = np.zeros((opt.num_samples, self.num_exps))
    self.exp_scores = np.zeros_like(self.exps_done)
    super().__init__(opt)

  def checkpoint_ref_setup(self):
    super().checkpoint_ref_setup()
    self.checkpoint_ref['extra'] = ['score', 'exps_done', 'exp_scores']

  def run(self, cmd_queue, result_queue):
    opt = self.opt
    self.base_cmd = self.cfg.base_cmd

    # Set up any additional arguments to tack on to experiments
    extra_args = []
    unparsed = opt.unparsed
    if '--' in unparsed:
      extra_args = unparsed[unparsed.index('--') + 1:]
    self.extra_child_args = extra_args

    print('Number of experiments:', self.num_exps)
    print('Number of trials:', opt.num_samples)
    print('Extra args:', extra_args)

    # Submit all experiments that haven't been run
    for trial_idx in range(opt.num_samples):
      for exp_idx, e in enumerate(self.cfg.exp_list):
        if not self.exps_done[trial_idx, exp_idx]:
          exp_id = 'trial_%d/%s' % (trial_idx, self.cfg.exp_names[exp_idx])
          if '--' in e:
            extra_child_args = e[e.index('--'):]
            e = e[:e.index('--')]
          else:
            extra_child_args = []

          self.submit_cmd(
              cmd_queue,
              exp_id,
              extra_args=e,
              extra_child_args=extra_child_args)

    # Collect results
    while self.exps_done.sum() != opt.num_samples * self.num_exps:
      result = result_queue.get()
      self.results += [result]

      exp_id = result[0].split('/')
      for tmp_idx, val in enumerate(exp_id):
        if 'trial_' in val:
          trial_idx = int(val.split('_')[-1])
          ref_idx = tmp_idx + 1
      exp_idx = self.cfg.exp_names.index('/'.join(exp_id[ref_idx:]))

      score = result[1]['score']
      self.exp_scores[trial_idx, exp_idx] = score
      self.exps_done[trial_idx, exp_idx] = 1

      if score > self.best_score:
        self.best_score = score
        self.best_exp = result[0]

      self.exp_count += 1

      print('Collected %s with score %.2f' % (result[0], score))

      if opt.num_samples < 100 or self.exp_count % 20 == 0:
        self.save(self.exp_dir + '/snapshot')

    self.save(self.exp_dir + '/snapshot')

max_stars_repo_name: google/multi-task-architecture-search


text:
import numpy as np

import chainer
from chainer import cuda
import chainer.functions as F
import chainer.links as L

def _augmentation(x):
    xp = cuda.get_array_module(x)
    MAX_SHIFT = 2
    batchsize, ch, h, w = x.shape
    h_shift, w_shift = xp.random.randint(-MAX_SHIFT, MAX_SHIFT + 1, size=2)
    a_h_sl = slice(max(0, h_shift), h_shift + h)
    a_w_sl = slice(max(0, w_shift), w_shift + w)
    x_h_sl = slice(max(0, - h_shift), - h_shift + h)
    x_w_sl = slice(max(0, - w_shift), - w_shift + w)
    a = xp.zeros(x.shape)
    a[:, :, a_h_sl, a_w_sl] = x[:, :, x_h_sl, x_w_sl]
    return a.astype(x.dtype)

def _count_params(m, n_grids=6):
    print('# of params', sum(param.size for param in m.params()))
    # The number of parameters in the paper (11.36M) might be
    # of the model with unshared matrices over primary capsules in a same grid
    # when input data are 36x36 images of MultiMNIST (n_grids = 10).
    # Our model with n_grids=10 has 11.349008M parameters.
    # (In the Sec. 4, the paper says "each capsule in the [6, 6] grid
    # is sharing their weights with each other.")
    print('# of params if unshared',
          sum(param.size for param in m.params()) +
          sum(param.size for param in m.Ws.params()) *
          (n_grids * n_grids - 1))

def squash(ss):
    ss_norm2 = F.sum(ss ** 2, axis=1, keepdims=True)
    """
    # ss_norm2 = F.broadcast_to(ss_norm2, ss.shape)
    # vs = ss_norm2 / (1. + ss_norm2) * ss / F.sqrt(ss_norm2): naive
    """
    norm_div_1pnorm2 = F.sqrt(ss_norm2) / (1. + ss_norm2)
    norm_div_1pnorm2 = F.broadcast_to(norm_div_1pnorm2, ss.shape)
    vs = norm_div_1pnorm2 * ss  # :efficient
    # (batchsize, 16, 10)
    return vs

def get_norm(vs):
    return F.sqrt(F.sum(vs ** 2, axis=1))

init = chainer.initializers.Uniform(scale=0.05)

class CapsNet(chainer.Chain):

    def __init__(self, use_reconstruction=False):
        super(CapsNet, self).__init__()
        self.n_iterations = 3  # dynamic routing
        self.n_grids = 6  # grid width of primary capsules layer
        self.n_raw_grids = self.n_grids
        self.use_reconstruction = use_reconstruction
        with self.init_scope():
            self.conv1 = L.Convolution2D(1, 256, ksize=9, stride=1,
                                         initialW=init)
            self.conv2 = L.Convolution2D(256, 32 * 8, ksize=9, stride=2,
                                         initialW=init)
            self.Ws = chainer.ChainList(
                *[L.Convolution2D(8, 16 * 10, ksize=1, stride=1, initialW=init)
                  for i in range(32)])

            self.fc1 = L.Linear(16 * 10, 512, initialW=init)
            self.fc2 = L.Linear(512, 1024, initialW=init)
            self.fc3 = L.Linear(1024, 784, initialW=init)

        _count_params(self, n_grids=self.n_grids)
        self.results = {'N': 0., 'loss': [], 'correct': [],
                        'cls_loss': [], 'rcn_loss': []}

    def pop_results(self):
        merge = dict()
        merge['mean_loss'] = sum(self.results['loss']) / self.results['N']
        merge['cls_loss'] = sum(self.results['cls_loss']) / self.results['N']
        merge['rcn_loss'] = sum(self.results['rcn_loss']) / self.results['N']
        merge['accuracy'] = sum(self.results['correct']) / self.results['N']
        self.results = {'N': 0., 'loss': [], 'correct': [],
                        'cls_loss': [], 'rcn_loss': []}
        return merge

    def __call__(self, x, t):
        if chainer.config.train:
            x = _augmentation(x)
        vs_norm, vs = self.output(x)
        self.loss = self.calculate_loss(vs_norm, t, vs, x)

        self.results['loss'].append(self.loss.data * t.shape[0])
        self.results['correct'].append(self.calculate_correct(vs_norm, t))
        self.results['N'] += t.shape[0]
        return self.loss

    def output(self, x):
        batchsize = x.shape[0]
        n_iters = self.n_iterations
        gg = self.n_grids * self.n_grids

        # h1 = F.relu(self.conv1(x))
        h1 = F.leaky_relu(self.conv1(x), 0.05)
        pr_caps = F.split_axis(self.conv2(h1), 32, axis=1)
        # shapes if MNIST. -> if MultiMNIST
        # x (batchsize, 1, 28, 28) -> (:, :, 36, 36)
        # h1 (batchsize, 256, 20, 20) -> (:, :, 28, 28)
        # pr_cap (batchsize, 8, 6, 6) -> (:, :, 10, 10)

        Preds = []
        for i in range(32):
            pred = self.Ws[i](pr_caps[i])
            Pred = pred.reshape((batchsize, 16, 10, gg))
            Preds.append(Pred)
        Preds = F.stack(Preds, axis=3)
        assert(Preds.shape == (batchsize, 16, 10, 32, gg))

        bs = self.xp.zeros((batchsize, 10, 32, gg), dtype='f')
        for i_iter in range(n_iters):
            cs = F.softmax(bs, axis=1)
            Cs = F.broadcast_to(cs[:, None], Preds.shape)
            assert(Cs.shape == (batchsize, 16, 10, 32, gg))
            ss = F.sum(Cs * Preds, axis=(3, 4))
            vs = squash(ss)
            assert(vs.shape == (batchsize, 16, 10))

            if i_iter != n_iters - 1:
                Vs = F.broadcast_to(vs[:, :, :, None, None], Preds.shape)
                assert(Vs.shape == (batchsize, 16, 10, 32, gg))
                bs = bs + F.sum(Vs * Preds, axis=1)
                assert(bs.shape == (batchsize, 10, 32, gg))

        vs_norm = get_norm(vs)
        return vs_norm, vs

    def reconstruct(self, vs, t):
        xp = self.xp
        batchsize = t.shape[0]
        I = xp.arange(batchsize)
        mask = xp.zeros(vs.shape, dtype='f')
        mask[I, :, t] = 1.
        masked_vs = mask * vs

        x_recon = F.sigmoid(
            self.fc3(F.relu(
                self.fc2(F.relu(
                    self.fc1(masked_vs)))))).reshape((batchsize, 1, 28, 28))
        return x_recon

    def calculate_loss(self, vs_norm, t, vs, x):
        class_loss = self.calculate_classification_loss(vs_norm, t)
        self.results['cls_loss'].append(class_loss.data * t.shape[0])
        if self.use_reconstruction:
            recon_loss = self.calculate_reconstruction_loss(vs, t, x)
            self.results['rcn_loss'].append(recon_loss.data * t.shape[0])
            return class_loss + 0.0005 * recon_loss
        else:
            return class_loss

    def calculate_classification_loss(self, vs_norm, t):
        xp = self.xp
        batchsize = t.shape[0]
        I = xp.arange(batchsize)
        T = xp.zeros(vs_norm.shape, dtype='f')
        T[I, t] = 1.
        m = xp.full(vs_norm.shape, 0.1, dtype='f')
        m[I, t] = 0.9

        loss = T * F.relu(m - vs_norm) ** 2 + \
            0.5 * (1. - T) * F.relu(vs_norm - m) ** 2
        return F.sum(loss) / batchsize

    def calculate_reconstruction_loss(self, vs, t, x):
        batchsize = t.shape[0]
        x_recon = self.reconstruct(vs, t)
        loss = (x_recon - x) ** 2
        return F.sum(loss) / batchsize

    def calculate_correct(self, v, t):
        return (self.xp.argmax(v.data, axis=1) == t).sum()

max_stars_repo_name: soskek/chainer-capsnet


text:
from rlbench.environment import Environment
from rlbench.action_modes import ArmActionMode, ActionMode
from rlbench.observation_config import ObservationConfig

import numpy as np
from romi.trajectory import NamedTrajectory
from romi.movement_primitives import MovementPrimitive, Group, ClassicSpace, \
    LearnTrajectory
from mppca.multiprocess import MultiProcess

from pyrep.backend import sim

from mppca.mixture_ppca import MPPCA
from core.model import RLModel
from core.task_interface import TaskInterface

import matplotlib.pyplot as plt

class ImitationLearning:

    def __init__(self, task_class, n_features, load, n_movements):
        """
        Learn the Movement from demonstration.

        :param task_class: Task that we aim to learn
        :param n_features: Number of RBF in n_features
        :param load: Load from data
        :param n_movements: how many movements do we want to learn
        """

        frequency = 200

        # To use 'saved' demos, set the path below, and set live_demos=False
        live_demos = not load
        DATASET = '' if live_demos else 'datasets'

        obs_config = ObservationConfig()
        obs_config.set_all_low_dim(True)
        obs_config.set_all_high_dim(False)
        self._task_class = task_class
        action_mode = ActionMode(ArmActionMode.ABS_JOINT_POSITION)

        group = Group("rlbench", ["d%d" % i for i in range(7)] + ["gripper"])

        env = Environment(
            action_mode, DATASET, obs_config, headless=True)
        env.launch()

        task = env.get_task(task_class)

        trajectories = []
        states = []

        lengths = []

        print("Start Demo")
        demos = task.get_demos(n_movements, live_demos=live_demos)
        print("End Demo")

        init = True
        for demo in demos:
            trajectory = NamedTrajectory(*group.refs)
            t = 0
            for ob in demo:
                if t == 0:
                    if init:
                        print("State dim: %d" % ob.task_low_dim_state.shape[0])
                        init = False
                    states.append(ob.task_low_dim_state)
                kwargs = {"d%d" % i: ob.joint_positions[i] for i in range(ob.joint_positions.shape[0])}
                kwargs["gripper"] = ob.gripper_open
                trajectory.notify(duration=1/frequency,
                                  **kwargs)
                t += 1
            lengths.append(t/200.)
            trajectories.append(trajectory)

        space = ClassicSpace(group, n_features=n_features, regularization=1E-15)
        z = np.linspace(-0.2, 1.2, 1000)
        Phi = space.get_phi(z)
        for i in range(n_features):
            plt.plot(z, Phi[:, i])
        plt.show()
        parameters = np.array([np.concatenate([s, np.array([l]), LearnTrajectory(space, traj).get_block_params()])
                               for s, l, traj in zip(states, lengths, trajectories)])
        np.save("parameters/%s_%d.npy" % (task.get_name(), n_features), parameters)
        env.shutdown()

class PPCAImitation:

    def __init__(self, task_class, state_dim, n_features, n_cluster, n_latent, parameters=None, headless=False,
                 cov_reg=1E-8, n_samples=50, dense_reward=False, imitation_noise=0.03):

        obs_config = ObservationConfig()
        obs_config.set_all_low_dim(True)
        obs_config.set_all_high_dim(False)
        self._obs_config = obs_config

        self._state_dim = state_dim
        self._headless = headless

        action_mode = ActionMode(ArmActionMode.ABS_JOINT_POSITION)
        self._task_class = task_class
        self._action_mode = action_mode
        self.env = Environment(
            action_mode, "", obs_config, headless=headless)
        self.env.launch()

        self.task = self.env.get_task(task_class)
        if parameters is None:
            self.parameters = np.load("parameters/%s_%d.npy" % (self.task.get_name(), n_features))[:n_samples]
        # parameters = np.concatenate([parameters for _ in range(20)])
        self.imitation_noise = imitation_noise
        self.parameters[:, :3] += imitation_noise * np.random.normal(size=self.parameters[:, :3].shape)
        self.mppca = MPPCA(n_cluster, n_latent, n_iterations=30, cov_reg=cov_reg, n_init=500)
        self.mppca.fit(self.parameters)

        self.rlmppca = None
        self.dense_reward = dense_reward

        print(np.exp(self.mppca.log_pi))

        group = Group("rlbench", ["d%d" % i for i in range(7)] + ["gripper"])
        self.space = ClassicSpace(group, n_features=n_features)

        print("mpcca learned")

    def stop(self, joint_gripper, previous_reward):
        if previous_reward == 0.:
            success = 0.
            for _ in range(50):
                obs, reward, terminate = self.task.step(joint_gripper)
                if reward == 1.0:
                    success = 1.
                    break
            return self.task._task.get_dense_reward(), success
        return self.task._task.get_dense_reward(), 1.

    def collect_samples(self, n_episodes=50, noise=True, isomorphic_noise=False):

        success_list = []
        reward_list = []
        latent = []
        cluster = []
        observations = []
        parameters = []
        ob = self.task.reset()
        for i in range(n_episodes):
            context = ob[1].task_low_dim_state
            observations.append(context)
            w, z, k = self.rlmppca.generate_full(np.expand_dims(context, 0), noise=noise, isomorphic_noise=isomorphic_noise)
            # w = self.parameters[i, 3:]
            parameters.append(w)
            latent.append(z)
            cluster.append(k)
            mp = MovementPrimitive(self.space, MovementPrimitive.get_params_from_block(self.space, w[1:]))
            duration = 1 if w[0] < 0 else w[0]
            print(duration)
            if self._headless:
                trajectory = mp.get_full_trajectory(duration=duration, frequency=200)
            else:
                trajectory = mp.get_full_trajectory(duration=5*duration, frequency=200)
            tot_reward = 0.
            success = 0
            for a1 in trajectory.values:  # , a2 in zip(trajectory.values[:-1, :], trajectory.values[1:, :]):
                joint = a1  # (a2-a1)*20.
                joint_gripper = joint
                obs, reward, terminate = self.task.step(joint_gripper)
                if reward == 1. or terminate==1.:
                    if reward == 1.:
                        success = 1.
                    break
            tot_reward, success = self.stop(joint_gripper, success)
            # tot_reward =  -np.mean(np.abs(context - ob[1].gripper_pose[:3]))
            # tot_reward = -(w[0]-0.2)**2
            print(tot_reward)
            # print("my reward", -np.mean(np.abs(context - ob[1].gripper_pose[:3])))
            # print("parameters", parameters[-1])
            success_list.append(success)

            if self.dense_reward:
                reward_list.append(tot_reward)
            else:
                reward_list.append(success)

            ob = self.task.reset()

        print("-"*50)
        print("Total reward", np.mean(reward_list))
        print("-"*50)
        return np.array(success_list), np.array(reward_list), np.array(parameters), np.array(latent), np.array(cluster), np.array(observations)

    def run_episode(self):
        self.env = Environment(
            self._action_mode, "", self._obs_config, headless=self._headless)
        self.env.launch()

        self.task = self.env.get_task(self._task_class)

    def shutdown(self):
        self.env.shutdown()

class RunModel:

    def __init__(self, task: TaskInterface, rl_model: RLModel, dense_reward=False):
        self.task = task
        self._rl_model = rl_model
        self.dense_reward = dense_reward

    def collect_samples(self, n_episodes=50, noise=True, isomorphic_noise=False):

        success_list = []
        reward_list = []
        latent = []
        cluster = []
        observations = []
        parameters = []

        for i in range(n_episodes):
            self.task.reset()
            context = self.task.read_context()
            observations.append(context)

            # print("RECEIVED CONTEXT", np.expand_dims(context, 0))
            w, z, k = self._rl_model.generate_full(np.expand_dims(context, 0),
                                                   noise=noise, isomorphic_noise=isomorphic_noise)
            print("parameters ::::", w)
            parameters.append(w)
            latent.append(z)
            cluster.append(k)

            success, tot_reward = self.task.send_movement(w[1:], w[0])
            print(success, tot_reward, k, z)
            success_list.append(success)

            if self.dense_reward:
                reward_list.append(tot_reward)
            else:
                reward_list.append(success)

        print("-"*50)
        print("Total reward", np.mean(reward_list))
        print("-"*50)
        return np.array(success_list), np.array(reward_list), np.array(parameters), np.array(latent),\
               np.array(cluster), np.array(observations)


text:
import os
import csv
import numpy as np

from tqdm import tqdm

from sklearn.utils import shuffle
from sklearn.model_selection import train_test_split

seed = 3535999445

def _rocstories(path):
    with open(path, encoding='utf_8') as f:
        f = csv.reader(f)
        st = []
        ct1 = []
        ct2 = []
        y = []
        for i, line in enumerate(tqdm(list(f), ncols=80, leave=False)):
            if i > 0:
                s = ' '.join(line[1:5])
                c1 = line[5]
                c2 = line[6]
                st.append(s)
                ct1.append(c1)
                ct2.append(c2)
                y.append(int(line[-1])-1)
        return st, ct1, ct2, y

def rocstories(data_dir, n_train=1497, n_valid=374):
    storys, comps1, comps2, ys = _rocstories(os.path.join(data_dir, 'cloze_test_val__spring2016.csv'))
    teX1, teX2, teX3, _ = _rocstories(os.path.join(data_dir, 'cloze_test_test__spring2016.csv'))
    tr_storys, va_storys, tr_comps1, va_comps1, tr_comps2, va_comps2, tr_ys, va_ys = train_test_split(storys, comps1, comps2, ys, test_size=n_valid, random_state=seed)
    trX1, trX2, trX3 = [], [], []
    trY = []
    for s, c1, c2, y in zip(tr_storys, tr_comps1, tr_comps2, tr_ys):
        trX1.append(s)
        trX2.append(c1)
        trX3.append(c2)
        trY.append(y)

    vaX1, vaX2, vaX3 = [], [], []
    vaY = []
    for s, c1, c2, y in zip(va_storys, va_comps1, va_comps2, va_ys):
        vaX1.append(s)
        vaX2.append(c1)
        vaX3.append(c2)
        vaY.append(y)
    trY = np.asarray(trY, dtype=np.int32)
    vaY = np.asarray(vaY, dtype=np.int32)
    return (trX1, trX2, trX3, trY), (vaX1, vaX2, vaX3, vaY), (teX1, teX2, teX3)


text:
#Libraries--------------------------------------------------------------------------------
import numpy as np
import pandas as pd
import pytest
from sklearn.datasets import make_blobs
from sklearn.preprocessing import KBinsDiscretizer
import random

import unittest

from quickclus.QuickClus import QuickClus

#Adding some random nulls--------------------------------------------------------------------
def add_random_nulls(df, null_prop = 0.1):
    ix = [(row, col) for row in range(df.shape[0]) for col in range(df.shape[1])]
    for row, col in random.sample(ix, int(round(null_prop * len(ix)))):
        df.iat[row, col] = np.nan

    return df

#Unit test------------------------------------------------------------------------------------

class QuickClusTestCase(unittest.TestCase):

    def setUp(self):
        #Numerical + categorical
        n_clusters = 3

        #Include random null values
        null_prop = 0.1
        add_nulls = True

        X, y = make_blobs(n_samples = 1000, n_features = 8, random_state = 1)

        numerical_features = X[:, :6]
        categorical_features = KBinsDiscretizer(n_bins = 3, encode = "ordinal").fit_transform(X[:, 6:])
        categorical_features = np.where(
            categorical_features == 1.0,
            "M",
            np.where(categorical_features == 2.0, "H", "L"),
        ).astype(str)

        numerical_columns = [f"num_{i}" for i in range(numerical_features.shape[1])]
        self.df = pd.DataFrame(numerical_features, columns = numerical_columns)

        categorical_columns = [f"cat_{i}" for i in range(categorical_features.shape[1])]
        for idx, c in enumerate(categorical_columns):
            self.df[c] = categorical_features[:, idx]

        self.df_numerical = pd.DataFrame(numerical_features, columns = numerical_columns)
        self.df_categorical = pd.DataFrame(categorical_features, columns = categorical_columns)

        if add_nulls:
            self.df = add_random_nulls(df = self.df, null_prop = null_prop)
            self.df_numerical = add_random_nulls(df = self.df_numerical, null_prop = null_prop)
            self.df_categorical = add_random_nulls(df = self.df_categorical, null_prop = null_prop)

        #Numerical + categorical
        self.clf = QuickClus(
            n_components = 2,
            random_state = 42,
            n_neighbors = 10,
            umap_combine_method = "intersection_union_mapper",
        )

    #Tests
    def test_fit_categorical(self):
        self.clf.fit(self.df)
        self.assertEqual(self.clf.umap_categorical_.embedding_.shape, (len(self.df), self.clf.n_components))

    def test_fit_numerical(self):
        self.clf.fit(self.df)
        self.assertEqual(self.clf.umap_numerical_.embedding_.shape, (len(self.df), self.clf.n_components))

    def test_fit_only_categorical(self):
        self.clf.fit(self.df_categorical)
        self.assertEqual(self.clf.umap_categorical_.embedding_.shape, (len(self.df_categorical), self.clf.n_components))

    def test_fit_only_numerical(self):
        self.clf.fit(self.df_numerical)
        self.assertEqual(self.clf.umap_numerical_.embedding_.shape, (len(self.df_numerical), self.clf.n_components))

    def test_umap_embeddings(self):
        self.clf.fit(self.df)
        self.assertEqual(self.clf.umap_combined.embedding_.shape, (len(self.df), self.clf.n_components))

    def test_hdbscan_labels(self):
        self.clf.fit(self.df)
        self.assertEqual(self.clf.hdbscan_.labels_.shape[0], self.df.shape[0])

    def test_hdbscan_optimization(self):
        self.clf.fit(self.df)
        self.clf.tune_model(n_trials = 5, min_cluster_start = 0.01, min_cluster_end = 0.05,
                    min_samples_start = 0.01, min_samples_end = 0.05, max_epsilon = None)
        self.assertEqual(self.clf.hdbscan_.labels_.shape[0], self.df.shape[0])

    def test_quickclus_fit_is_df(self):
        with pytest.raises(TypeError):
            self.clf.fit([1, 2, 3])

    def test_quickclus_method(self):
        with pytest.raises(KeyError):
            _ = QuickClus(umap_combine_method = "notamethod").fit(self.df)

    def test_repr(self):
        self.assertEqual(str(type(self.clf.__repr__)), "<class 'method'>")

    def test_condensed_tree(self):
        self.clf.fit(self.df)
        self.clf.plot_condensed_tree()

    def test_embedding_labels(self):
        self.clf.fit(self.df)
        self.clf.plot_embedding_labels()

    def test_2d_labels(self):
        self.clf.fit(self.df)
        self.clf.plot_2d_labels()

    def test_assing_results(self):
        self.clf.fit(self.df)
        self.results = self.clf.assing_results(self.df)
        self.assertEqual(self.df.shape[0], self.results.shape[0])
        self.assertEqual(self.df.shape[1] + 1, self.results.shape[1])

if __name__ == '__main__':
    unittest.main()

max_stars_repo_name: aitorporcel/QuickClus

max_issues_repo_name: aitorporcel/QuickClus

max_forks_repo_name: aitorporcel/QuickClus


text:
"""
Version 1,
Tester utils, for assessing convergenge accuracy
@author: enzo
"""

from dedicate_code.custom_funcs import benchmark
from dedicate_code.config import reportings_dir
from dedicate_code.setup_testSet import setup_dict

from dedicate_code.data_etl.etl_utils import gen_nonuniform_segment_delta, gen_nonuniform_chebyshev_delta, gen_x_nnunif_mesh, gen_x_mesh, gen_uniform_delta
from dedicate_code.feature_eng.solver_util import de_f_time

#%% Set style
import matplotlib as mpl
axtickfsize = 16
labelfsize = 20
legfsize = labelfsize - 2
txtfsize = labelfsize - 2
lwidth = 3
markersize = 10
markeredgewidth = 0.1
mpl.rcParams['axes.titlesize'] = 24
mpl.rcParams['axes.labelsize'] = labelfsize
mpl.rcParams['xtick.labelsize'] = axtickfsize
mpl.rcParams['ytick.labelsize'] = axtickfsize
mpl.rcParams['font.size'] = txtfsize
mpl.rcParams['legend.fontsize'] = legfsize
mpl.rcParams['lines.linewidth'] = lwidth
mpl.rcParams['lines.markersize'] = markersize
mpl.rcParams['lines.markeredgewidth'] = markeredgewidth

#%% Functions
def plot_final_funcs_x(num_of_cells, find_anlyt_error_func, timeDuration):
    # higher-order function
    import matplotlib.pyplot as plt
    fig, axs = plt.subplots(3, 3, figsize=(16, 9))#
    #1920x1080 (or 1080p)  figsize=(19.20,10.80)  (16, 9)

    formula=setup_dict['eq_formula']
    specifc=setup_dict['text_long']
    title_txt = f'{formula} \n {specifc}'  #
    fig.suptitle(title_txt)

    delta_space = gen_uniform_delta(num_of_cells)
    mesh, x =gen_x_mesh(num_of_cells, delta_space)
    setup_dict['grid_spacing'] = delta_space
    setup_dict['center_spacing'] = delta_space
    three_time_phi = de_f_time(1, setup_dict, mesh, None, timeDuration)  #False  True

    error_u, analyt_unif = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)
    ax=axs[2][0]
    ax.plot(x[0], analyt_unif, label="analytical")
    ax.plot(x[0], three_time_phi[1][-1], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f'T={timeDuration}, err={error_u:.5f}')

    error_u, analyt_unif = find_anlyt_error_func(mesh, three_time_phi[1][-2], timeDuration/2)
    ax=axs[1][0]
    ax.plot(x[0], analyt_unif, label="analytical")
    ax.plot(x[0], three_time_phi[1][-2], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f'T={timeDuration/2}, err={error_u:.5f}')

    error_u, analyt_unif = find_anlyt_error_func(mesh, three_time_phi[1][-3], timeDuration/10)
    ax=axs[0][0]
    ax.plot(x[0], analyt_unif, label="analytical")
    ax.plot(x[0], three_time_phi[1][-2], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f"Uniform mesh {num_of_cells}\n"+f'T={timeDuration/10}, err={error_u:.5f}')

    delta_vect, delta_space = gen_nonuniform_segment_delta(num_of_cells)
    mesh, x =gen_x_nnunif_mesh(delta_vect)
    setup_dict['grid_spacing'] = delta_space
    setup_dict['center_spacing']  = delta_space
    three_time_phi = de_f_time(1, setup_dict, mesh, None, timeDuration)  #False  True

    errors_nu_seg, analyt_nonun = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)
    ax=axs[2][1]
    ax.plot(x[0], analyt_nonun, label="analytical")
    ax.plot(x[0], three_time_phi[1][-1], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f'T={timeDuration}, err={errors_nu_seg:.5f}')

    errors_nu_seg, analyt_nonun = find_anlyt_error_func(mesh, three_time_phi[1][-2], timeDuration/2)
    ax=axs[1][1]
    ax.plot(x[0], analyt_nonun, label="analytical")
    ax.plot(x[0], three_time_phi[1][-2], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f'T={timeDuration/2}, err={errors_nu_seg:.5f}')

    errors_nu_seg, analyt_nonun = find_anlyt_error_func(mesh, three_time_phi[1][-3], timeDuration/10)
    ax=axs[0][1]
    ax.plot(x[0], analyt_nonun, label="analytical")
    ax.plot(x[0], three_time_phi[1][-2], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f"Non-unif (segments) {num_of_cells} mesh\n"+f'T={timeDuration/10}, err={errors_nu_seg:.5f}')

    delta_vect, delta_space = gen_nonuniform_chebyshev_delta(num_of_cells)
    mesh, x =gen_x_nnunif_mesh(delta_vect)
    setup_dict['grid_spacing'] = max(delta_vect)
    setup_dict['center_spacing']  = max(delta_vect)
    three_time_phi = de_f_time(1, setup_dict, mesh, None, timeDuration)  #False  True

    errors_nu_chb, analyt_cheb = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)
    ax=axs[2][2]
    ax.plot(x[0], analyt_cheb, label="analytical")
    ax.plot(x[0], three_time_phi[1][-1], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f'T={timeDuration}, err={errors_nu_chb:.5f}')

    errors_nu_chb, analyt_cheb = find_anlyt_error_func(mesh, three_time_phi[1][-2], timeDuration/2)
    ax=axs[1][2]
    ax.plot(x[0], analyt_cheb, label="analytical")
    ax.plot(x[0], three_time_phi[1][-2], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f'T={timeDuration/2}, err={errors_nu_chb:.5f}')

    errors_nu_chb, analyt_cheb = find_anlyt_error_func(mesh, three_time_phi[1][-3], timeDuration/10)
    ax=axs[0][2]
    ax.plot(x[0], analyt_cheb, label="analytical")
    ax.plot(x[0], three_time_phi[1][-2], linestyle="--", label="u")
    ax.legend(loc="best")
    ax.set_xlabel("x")
    ax.set_ylabel("u")
    ax.set_title(f"Non-unif (Chebyshev) {num_of_cells} mesh\n"+f'T={timeDuration/10}, err={errors_nu_chb:.5f}')

    fig.tight_layout()

    #%% SAVE
    file_name = setup_dict['compactxt'] + f'_{num_of_cells}'
    folder_name = setup_dict['id']
    report_folder = reportings_dir/folder_name
    try:
        report_folder.mkdir(parents=True, exist_ok=False)
    except FileExistsError:
        print(f"Folder {report_folder} is already there")
    else:
        print(f"Folder {report_folder} was created")

    # for pubblication: Nature: Postscript, Vector EPS or PDF format
    #maximum 300 p.p.i.; and min 300dpi
    fig.savefig(report_folder/(file_name+'.pdf'), bbox_inches="tight")
    fig.savefig(report_folder/(file_name+'.jpg'), bbox_inches="tight")
    fig.savefig(report_folder/(file_name+'.png'), bbox_inches="tight")
    fig.savefig(report_folder/(file_name+'.eps'), format='eps', dpi=300)
    fig.savefig(report_folder/(file_name+'.svg'), format='svg', dpi=300)

    return fig

def initializer(nxvec):
    import numpy as np
    errors_u = np.zeros(len(nxvec))
    errors_nu_seg = np.zeros(len(nxvec))
    errors_nu_chb = np.zeros(len(nxvec))
    return errors_u, errors_nu_seg, errors_nu_chb

def init_pool(err_u, err_nu_seg, err_nu_chb):
    global errors_u, errors_nu_seg, errors_nu_chb
    errors_u = err_u
    errors_nu_seg = err_nu_seg
    errors_nu_chb = err_nu_chb

def gen_errs(indx, num_of_cells, find_anlyt_error_func, timeDuration):
    delta_space = gen_uniform_delta(num_of_cells)
    mesh, x =gen_x_mesh(num_of_cells, delta_space)
    setup_dict['grid_spacing'] = delta_space
    setup_dict['center_spacing'] = delta_space
    three_time_phi = de_f_time(indx, setup_dict, mesh, None, timeDuration)  #False  True
    errors_u, _ = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)

    delta_vect, delta_space = gen_nonuniform_segment_delta(num_of_cells)
    mesh, x =gen_x_nnunif_mesh(delta_vect)
    setup_dict['grid_spacing'] = delta_space
    setup_dict['center_spacing']  = delta_space
    three_time_phi = de_f_time(indx, setup_dict, mesh, None, timeDuration)  #False  True
    errors_nu_seg, _ = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)

    delta_vect, delta_space = gen_nonuniform_chebyshev_delta(num_of_cells)
    mesh, x =gen_x_nnunif_mesh(delta_vect)
    setup_dict['grid_spacing'] =  max(delta_vect)
    setup_dict['center_spacing']  = max(delta_vect)
    three_time_phi = de_f_time(indx, setup_dict, mesh, None, timeDuration)  #False  True
    errors_nu_chb, _ = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)
    return indx, errors_u, errors_nu_seg, errors_nu_chb

def merge_result_err(result):
    indx = result[0]
    global errors_u, errors_nu_seg, errors_nu_chb
    errors_u[indx] = result[1]
    errors_nu_seg[indx] = result[2]
    errors_nu_chb[indx] = result[3]

@benchmark
def error_multicore(nxvec, find_anlyt_error_func, timeDuration):
    import multiprocessing as mp
    from functools import partial
    global errors_u, errors_nu_seg, errors_nu_chb
    errors_u, errors_nu_seg, errors_nu_chb = initializer(nxvec)
    pool = mp.Pool(mp.cpu_count()-1, init_pool(errors_u, errors_nu_seg, errors_nu_chb))
    for indx, num_of_cells in enumerate(nxvec):
        pool.apply_async(gen_errs,
        args=(indx, num_of_cells, find_anlyt_error_func, timeDuration),
        callback=merge_result_err)
    pool.close()
    pool.join()
    return errors_u, errors_nu_seg, errors_nu_chb

@benchmark
def error_onecore(nxvec, find_anlyt_error_func, timeDuration):
    global errors_u, errors_nu_seg, errors_nu_chb
    errors_u, errors_nu_seg, errors_nu_chb = initializer(nxvec)
    for indx, num_of_cells in enumerate(nxvec):
        merge_result_err(gen_errs(indx, num_of_cells, find_anlyt_error_func, timeDuration))
    return errors_u, errors_nu_seg, errors_nu_chb

def calulate_all_errors(nxvec, find_anlyt_error_func, timeDuration):
    return error_multicore(nxvec, find_anlyt_error_func, timeDuration)

def  calulate_all_errors_bkup(nxvec, find_anlyt_error_func, timeDuration):
    # higher-order function
    import numpy as np

    errors_u = np.zeros(len(nxvec))
    errors_nu_seg = np.zeros(len(nxvec))
    errors_nu_chb = np.zeros(len(nxvec))
    for indx, num_of_cells in enumerate(nxvec):
        delta_space = gen_uniform_delta(num_of_cells)
        mesh, x =gen_x_mesh(num_of_cells, delta_space)
        setup_dict['grid_spacing'] = delta_space
        setup_dict['center_spacing'] = delta_space
        three_time_phi = de_f_time(1, setup_dict, mesh, None, timeDuration)  #False  True
        errors_u[indx], _ = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)

        delta_vect, delta_space = gen_nonuniform_segment_delta(num_of_cells)
        mesh, x =gen_x_nnunif_mesh(delta_vect)
        setup_dict['grid_spacing'] = delta_space
        setup_dict['center_spacing']  = delta_space
        three_time_phi = de_f_time(1, setup_dict, mesh, None, timeDuration)  #False  True
        errors_nu_seg[indx], _ = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)

        delta_vect, delta_space = gen_nonuniform_chebyshev_delta(num_of_cells)
        mesh, x =gen_x_nnunif_mesh(delta_vect)
        setup_dict['grid_spacing'] =  max(delta_vect)
        setup_dict['center_spacing']  = max(delta_vect)
        three_time_phi = de_f_time(1, setup_dict, mesh, None, timeDuration)  #False  True
        errors_nu_chb[indx], _ = find_anlyt_error_func(mesh, three_time_phi[1][-1], timeDuration)
    return  errors_u, errors_nu_seg, errors_nu_chb

def plot_conv_acc(timeDuration, nxvec, errors_u, errors_nu_seg, errors_nu_chb):
    import numpy as np
    import matplotlib.pyplot as plt

    #%%  Show convergence accuracy
    log_nxvec = np.log(nxvec)
    log_errors_u = np.log(errors_u)
    log_errors_nu_seg = np.log(errors_nu_seg)
    log_errors_nu_chb = np.log(errors_nu_chb)
    fit_u = np.polyfit(log_nxvec, log_errors_u, 1)
    fit_nu_seg = np.polyfit(log_nxvec, log_errors_nu_seg, 1)
    fit_nu_chb = np.polyfit(log_nxvec, log_errors_nu_chb, 1)

    fig, axs = plt.subplots(1, 3, figsize=(16, 9))#
    #1920x1080 (or 1080p)  figsize=(19.20,10.80)  (16, 9)

    formula=setup_dict['eq_formula']
    specifc=setup_dict['text_long']
    title_txt = f'{formula} \n {specifc}; T={timeDuration}'
    fig.suptitle(title_txt)

    ax = axs[0]
    ax.plot(log_nxvec, log_errors_u, label="simulation errors")
    ax.plot(log_nxvec, np.poly1d(fit_u)(log_nxvec), linestyle="--",
            label="Line: {m:1.2f}x+{b:1.2f}".format(m=fit_u[0], b=fit_u[1]))
    ax.legend(loc="best")
    ax.set_xlabel("log(nx)")
    ax.set_ylabel("log(error)")
    ax.set_title("Uniform mesh")

    ax = axs[1]
    ax.plot(log_nxvec, log_errors_nu_seg, label="simulation errors")
    ax.plot(log_nxvec, np.poly1d(fit_nu_seg)(log_nxvec), linestyle="--",
            label="Line: {m:1.2f}x+{b:1.2f}".format(m=fit_nu_seg[0], b=fit_nu_seg[1]))
    ax.legend(loc="best")
    ax.set_xlabel("log(nx)")
    ax.set_ylabel("log(error)")
    ax.set_title("Non-uniform \n (segments) mesh")

    ax = axs[2]
    ax.plot(log_nxvec, log_errors_nu_chb, label="simulation errors")
    ax.plot(log_nxvec, np.poly1d(fit_nu_chb)(log_nxvec), linestyle="--",
            label="Line: {m:1.2f}x+{b:1.2f}".format(m=fit_nu_chb[0], b=fit_nu_chb[1]))
    ax.legend(loc="best")
    ax.set_xlabel("log(nx)")
    ax.set_ylabel("log(error)")
    ax.set_title("Non-uniform \n (Chebyshev) mesh")

    fig.tight_layout()

    #%% SAVE
    file_name = setup_dict['compactxt']
    folder_name = setup_dict['id']
    report_folder = reportings_dir/folder_name
    try:
        report_folder.mkdir(parents=True, exist_ok=False)
    except FileExistsError:
        print(f"Folder {report_folder} is already there")
    else:
        print(f"Folder {report_folder} was created")

    # for pubblication: Nature: Postscript, Vector EPS or PDF format
    #maximum 300 p.p.i.; and min 300dpi
    fig.savefig(report_folder/(file_name+'.pdf'), bbox_inches="tight")
    fig.savefig(report_folder/(file_name+'.jpg'), bbox_inches="tight")
    fig.savefig(report_folder/(file_name+'.png'), bbox_inches="tight")
    fig.savefig(report_folder/(file_name+'.eps'), format='eps', dpi=300)
    fig.savefig(report_folder/(file_name+'.svg'), format='svg', dpi=300)

    return fig

max_stars_repo_name: MthBr/DEDICATE

max_issues_repo_name: MthBr/DEDICATE

max_forks_repo_name: MthBr/DEDICATE


text:
# pylint: disable=no-self-use,invalid-name,protected-access
import numpy
import pytest
import torch
from torch.autograd import Variable

from allennlp.modules.span_extractors import SpanExtractor, BidirectionalEndpointSpanExtractor
from allennlp.common.checks import ConfigurationError
from allennlp.common.params import Params
from allennlp.nn.util import batched_index_select

class TestBidirectonalEndpointSpanExtractor:
    def test_bidirectional_endpoint_span_extractor_can_build_from_params(self):
        params = Params({
                "type": "bidirectional_endpoint",
                "input_dim": 4,
                "num_width_embeddings": 5,
                "span_width_embedding_dim": 3
                })
        extractor = SpanExtractor.from_params(params)
        assert isinstance(extractor, BidirectionalEndpointSpanExtractor)
        assert extractor.get_output_dim() == 2 + 2 + 3

    def test_raises_on_odd_input_dimension(self):
        with pytest.raises(ConfigurationError):
            _ = BidirectionalEndpointSpanExtractor(7)

    def test_correct_sequence_elements_are_embedded(self):
        sequence_tensor = Variable(torch.randn([2, 5, 8]))
        # concatentate start and end points together to form our representation
        # for both the forward and backward directions.
        extractor = BidirectionalEndpointSpanExtractor(input_dim=8,
                                                       forward_combination="x,y",
                                                       backward_combination="x,y")
        indices = Variable(torch.LongTensor([[[1, 3],
                                              [2, 4]],
                                             [[0, 2],
                                              [3, 4]]]))

        span_representations = extractor(sequence_tensor, indices)

        assert list(span_representations.size()) == [2, 2, 16]
        assert extractor.get_output_dim() == 16
        assert extractor.get_input_dim() == 8

        # We just concatenated the start and end embeddings together, so
        # we can check they match the original indices if we split them apart.
        (forward_start_embeddings, forward_end_embeddings,
         backward_start_embeddings, backward_end_embeddings) = span_representations.split(4, -1)

        forward_sequence_tensor, backward_sequence_tensor = sequence_tensor.split(4, -1)

        # Forward direction => subtract 1 from start indices to make them exlusive.
        correct_forward_start_indices = Variable(torch.LongTensor([[0, 1],
                                                                   [-1, 2]]))
        # This index should be -1, so it will be replaced with a sentinel. Here,
        # we'll set it to a value other than -1 so we can index select the indices and
        # replace it later.
        correct_forward_start_indices[1, 0] = 1

        # Forward direction => end indices are the same.
        correct_forward_end_indices = Variable(torch.LongTensor([[3, 4], [2, 4]]))

        # Backward direction => start indices are exclusive, so add 1 to the end indices.
        correct_backward_start_indices = Variable(torch.LongTensor([[4, 5], [3, 5]]))
        # These exclusive end indices are outside the tensor, so will be replaced with the end sentinel.
        # Here we replace them with ones so we can index select using these indices without torch
        # complaining.
        correct_backward_start_indices[0, 1] = 1
        correct_backward_start_indices[1, 1] = 1
        # Backward direction => end indices are inclusive and equal to the forward start indices.
        correct_backward_end_indices = Variable(torch.LongTensor([[1, 2], [0, 3]]))

        correct_forward_start_embeddings = batched_index_select(forward_sequence_tensor.contiguous(),
                                                                correct_forward_start_indices)
        # This element had sequence_tensor index of 0, so it's exclusive index is the start sentinel.
        correct_forward_start_embeddings[1, 0] = extractor._start_sentinel.data
        numpy.testing.assert_array_equal(forward_start_embeddings.data.numpy(),
                                         correct_forward_start_embeddings.data.numpy())

        correct_forward_end_embeddings = batched_index_select(forward_sequence_tensor.contiguous(),
                                                              correct_forward_end_indices)
        numpy.testing.assert_array_equal(forward_end_embeddings.data.numpy(),
                                         correct_forward_end_embeddings.data.numpy())

        correct_backward_end_embeddings = batched_index_select(backward_sequence_tensor.contiguous(),
                                                               correct_backward_end_indices)
        numpy.testing.assert_array_equal(backward_end_embeddings.data.numpy(),
                                         correct_backward_end_embeddings.data.numpy())

        correct_backward_start_embeddings = batched_index_select(backward_sequence_tensor.contiguous(),
                                                                 correct_backward_start_indices)
        # This element had sequence_tensor index == sequence_tensor.size(1),
        # so it's exclusive index is the end sentinel.
        correct_backward_start_embeddings[0, 1] = extractor._end_sentinel.data
        correct_backward_start_embeddings[1, 1] = extractor._end_sentinel.data
        numpy.testing.assert_array_equal(backward_start_embeddings.data.numpy(),
                                         correct_backward_start_embeddings.data.numpy())

    def test_correct_sequence_elements_are_embedded_with_a_masked_sequence(self):
        sequence_tensor = Variable(torch.randn([2, 5, 8]))
        # concatentate start and end points together to form our representation
        # for both the forward and backward directions.
        extractor = BidirectionalEndpointSpanExtractor(input_dim=8,
                                                       forward_combination="x,y",
                                                       backward_combination="x,y")
        indices = Variable(torch.LongTensor([[[1, 3],
                                              [2, 4]],
                                             # This span has an end index at the
                                             # end of the padded sequence.
                                             [[0, 2],
                                              [0, 1]]]))
        sequence_mask = Variable(torch.LongTensor([[1, 1, 1, 1, 1],
                                                   [1, 1, 1, 0, 0]]))

        span_representations = extractor(sequence_tensor, indices, sequence_mask=sequence_mask)

        # We just concatenated the start and end embeddings together, so
        # we can check they match the original indices if we split them apart.
        (forward_start_embeddings, forward_end_embeddings,
         backward_start_embeddings, backward_end_embeddings) = span_representations.split(4, -1)

        forward_sequence_tensor, backward_sequence_tensor = sequence_tensor.split(4, -1)

        # Forward direction => subtract 1 from start indices to make them exlusive.
        correct_forward_start_indices = Variable(torch.LongTensor([[0, 1],
                                                                   [-1, -1]]))
        # These indices should be -1, so they'll be replaced with a sentinel. Here,
        # we'll set them to a value other than -1 so we can index select the indices and
        # replace them later.
        correct_forward_start_indices[1, 0] = 1
        correct_forward_start_indices[1, 1] = 1

        # Forward direction => end indices are the same.
        correct_forward_end_indices = Variable(torch.LongTensor([[3, 4], [2, 1]]))

        # Backward direction => start indices are exclusive, so add 1 to the end indices.
        correct_backward_start_indices = Variable(torch.LongTensor([[4, 5], [3, 2]]))
        # These exclusive backward start indices are outside the tensor, so will be replaced
        # with the end sentinel. Here we replace them with ones so we can index select using
        # these indices without torch complaining.
        correct_backward_start_indices[0, 1] = 1

        # Backward direction => end indices are inclusive and equal to the forward start indices.
        correct_backward_end_indices = Variable(torch.LongTensor([[1, 2], [0, 0]]))

        correct_forward_start_embeddings = batched_index_select(forward_sequence_tensor.contiguous(),
                                                                correct_forward_start_indices)
        # This element had sequence_tensor index of 0, so it's exclusive index is the start sentinel.
        correct_forward_start_embeddings[1, 0] = extractor._start_sentinel.data
        correct_forward_start_embeddings[1, 1] = extractor._start_sentinel.data
        numpy.testing.assert_array_equal(forward_start_embeddings.data.numpy(),
                                         correct_forward_start_embeddings.data.numpy())

        correct_forward_end_embeddings = batched_index_select(forward_sequence_tensor.contiguous(),
                                                              correct_forward_end_indices)
        numpy.testing.assert_array_equal(forward_end_embeddings.data.numpy(),
                                         correct_forward_end_embeddings.data.numpy())

        correct_backward_end_embeddings = batched_index_select(backward_sequence_tensor.contiguous(),
                                                               correct_backward_end_indices)
        numpy.testing.assert_array_equal(backward_end_embeddings.data.numpy(),
                                         correct_backward_end_embeddings.data.numpy())

        correct_backward_start_embeddings = batched_index_select(backward_sequence_tensor.contiguous(),
                                                                 correct_backward_start_indices)
        # This element had sequence_tensor index == sequence_tensor.size(1),
        # so it's exclusive index is the end sentinel.
        correct_backward_start_embeddings[0, 1] = extractor._end_sentinel.data
        # This element has sequence_tensor index == the masked length of the batch element,
        # so it should be the end_sentinel even though it isn't greater than sequence_tensor.size(1).
        correct_backward_start_embeddings[1, 0] = extractor._end_sentinel.data

        numpy.testing.assert_array_equal(backward_start_embeddings.data.numpy(),
                                         correct_backward_start_embeddings.data.numpy())

    def test_forward_raises_with_invalid_indices(self):
        sequence_tensor = Variable(torch.randn([2, 5, 8]))
        extractor = BidirectionalEndpointSpanExtractor(input_dim=8)
        indices = Variable(torch.LongTensor([[[-1, 3],
                                              [7, 4]],
                                             [[0, 12],
                                              [0, -1]]]))

        with pytest.raises(ValueError):
            _ = extractor(sequence_tensor, indices)

max_stars_repo_name: unendin/allennlp

max_issues_repo_name: unendin/allennlp

max_forks_repo_name: unendin/allennlp


text:
import numpy as np
import pandas as pd

from flight_fusion import FusionServiceClient
from flight_fusion.asset_key import AssetKey
from flight_fusion.ipc.v1alpha1 import SaveMode

def test_roundtrip(fusion_client: FusionServiceClient):
    np.random.seed(42)
    df = pd.DataFrame(np.random.randn(5, 3), columns=["col1", "col2", "col3"])

    fds = fusion_client.get_dataset_client(AssetKey(["test_roundtrip", "table"]))
    fds.write_into(df, SaveMode.SAVE_MODE_OVERWRITE)

    df_loaded = fds.load().to_pandas()

    assert df.equals(df_loaded)

def test_roundtrip_versioned(fusion_client: FusionServiceClient):
    np.random.seed(42)
    df = pd.DataFrame(np.random.randn(5, 3), columns=["col1", "col2", "col3"])

    fds = fusion_client.get_versioned_dataset_client(
        AssetKey(["test_roundtrip_versioned", "table"])
    )
    fds.write_into(df, SaveMode.SAVE_MODE_OVERWRITE)

    df_loaded = fds.load().to_pandas()

    assert df.equals(df_loaded)

def test_save_mode(fusion_client: FusionServiceClient):
    np.random.seed(42)
    df = pd.DataFrame(np.random.randn(5, 3), columns=["col1", "col2", "col3"])

    fds = fusion_client.get_dataset_client(AssetKey(["test_save_mode", "table"]))

    fds.write_into(df, SaveMode.SAVE_MODE_OVERWRITE)
    df_loaded = fds.load().to_pandas()
    assert df_loaded.shape == (5, 3)

    fds.write_into(df, SaveMode.SAVE_MODE_OVERWRITE)
    df_loaded = fds.load().to_pandas()
    assert df_loaded.shape == (5, 3)

    fds.write_into(df, SaveMode.SAVE_MODE_APPEND)
    df_loaded = fds.load().to_pandas()
    assert df_loaded.shape == (10, 3)

max_stars_repo_name: roeap/flight-fusion

max_issues_repo_name: roeap/flight-fusion

max_forks_repo_name: roeap/flight-fusion


text:
#!/usr/bin/env python
# -*- encoding: utf-8 -*-
import argparse
import multiprocessing
import time

import gc
import numpy as np

__doc__ = """
splice_profile - a scripts for dataset construction
====================================================

**splice_profile** is a Python script for building dataset with multi-view feature representations.

Main Functions
--------------
Here are just a few of the things that **splice_profile** does well:

  - Multi-process building data sets in equal parts.
  - A single process builds a data set.

Main Program Functions
----------------------

"""

def assembledata_sp10(ddir:str, n=10, *, mulit=True, stop_num=1):
    '''Multi-process builds j times ten equal data sets.

    Args:

        ddir: Data path.

        n: The number of equal parts.

        mulit: Whether to multiprocess.

        stop_num: Processing rounds.

    Returns:

    '''
    checkdata = np.load(ddir + "/portion/ncRNApair_data0.npy")
    dsize = len(checkdata) + 1
    genelen = len(checkdata[0])
    width = len(checkdata[0][0])

    def process_data10(n, j):
        data = np.empty((0, genelen, width), dtype=np.float32)
        label = np.empty((0), dtype=np.int32)
        for i in range(dsize):
            indata = ddir + "/portion/ncRNApair_data" + str(i) + ".npy"
            inlabel = ddir + "/portion/ncRNApair_label" + str(i) + ".npy"

            tmpdata = np.load(indata)
            tmplabel = np.load(inlabel)

            tmpdata_10 = np.array_split(tmpdata, n, axis=0)
            tmplabel_10 = np.array_split(tmplabel, n, axis=0)

            data = np.concatenate([data, tmpdata_10[j]], axis=0)
            label = np.concatenate([label, tmplabel_10[j]], axis=0)
            del tmpdata_10, tmplabel_10
            gc.collect()

        np.save(ddir + "/data_10/ncRNApair_data" + str(j) + ".npy", data)
        np.save(ddir + "/data_10/ncRNApair_label" + str(j) + ".npy", label)
        del data, label
        gc.collect()

    if mulit:
        jobs = []
        for j in range(stop_num):
            p = multiprocessing.Process(target=process_data10, args=(n, j))
            jobs.append(p)
            p.start()

        for job in jobs:
            while job.is_alive() == True:
                continue
            else:
                pass

    else:
        for j in range(stop_num):
            process_data10(n, j)

    return

def assembledata(ddir):
    '''Build a data set that contains all the data.

    Args:

        ddir: Data path.

    Returns:

    '''
    checkdata = np.load(ddir + "/portion/ncRNApair_data0.npy")
    dsize = len(checkdata) + 1
    genelen = len(checkdata[0])
    width = len(checkdata[0][0])

    data = np.empty((0, genelen, width), dtype=np.float32)
    label = np.empty((0), dtype=np.int32)

    print("Makeing data...")
    start = time.time()

    for i in range(dsize):
        indata = ddir + "/portion/ncRNApair_data" + str(i) + ".npy"
        inlabel = ddir + "/portion/ncRNApair_label" + str(i) + ".npy"

        tmpdata = np.load(indata)
        tmplabel = np.load(inlabel)

        data = np.concatenate([data, tmpdata], axis=0)
        label = np.concatenate([label, tmplabel], axis=0)

    elapsed_time = time.time() - start
    print("Loding time:{0}".format(elapsed_time) + "[sec]")
    print("Data size :", len(data), "*", len(data[0]), "*", len(data[0][0]))

    np.save(ddir + "/ncRNApair_data.npy", data)
    np.save(ddir + "/ncRNApair_label.npy", label)

def main(datadir):
    print("Makeing data...")
    start = time.time()

    assembledata(datadir)

    print("Complete makeing data.")
    spend_time = time.time() - start
    print("Loding time:{0}".format(spend_time) + "[sec]")

if __name__ == '__main__':

    parser = argparse.ArgumentParser(description='splice_profile:')
    parser.add_argument('--datadir', '-d',
                        default='../features/profile/',
                        help='paths for multi-view features')
    args = parser.parse_args()

    main(args.datadir)

max_stars_repo_name: oAzv/GCFM

max_issues_repo_name: oAzv/GCFM

max_forks_repo_name: oAzv/GCFM


text:
import pytest
from energy_pandas import EnergyDataFrame, EnergySeries
from numpy.testing import assert_almost_equal
from pandas import read_csv

from archetypal import IDF, settings

@pytest.fixture(scope="module")
def idf(config):
    idfname = "tests/input_data/umi_samples/B_Off_0.idf"
    epw = "tests/input_data/CAN_PQ_Montreal.Intl.AP.716270_CWEC.epw"
    _idf = IDF(idfname, epw)
    yield _idf
    _idf.save()

@pytest.fixture(scope="class")
def rd_es(idf):
    try:
        yield idf.meters.OutputMeter.WaterSystems__MainsWater.values()
    except Exception:
        yield idf.simulate().meters.OutputMeter.WaterSystems__MainsWater.values()

@pytest.fixture(scope="class")
def rd_edf(idf):
    try:
        yield idf.variables.OutputVariable.Zone_Air_Temperature.values()
    except Exception:
        yield idf.simulate().variables.OutputVariable.Zone_Air_Temperature.values()

class TestEnergySeries:
    def test_from_csv(self):
        file = "tests/input_data/test_profile.csv"
        df = read_csv(file, index_col=[0], names=["Heat"])
        ep = EnergySeries.with_timeindex(df.Heat, frequency="1H", units="BTU/hour")
        assert ep.units == settings.unit_registry.parse_expression("BTU/hour").units

    def test_from_report_data(self, rd_es):
        assert_almost_equal(rd_es.sum(), 2.1187133811706036, decimal=3)
        assert rd_es.units == settings.unit_registry.m3

    @pytest.mark.parametrize("kind", ["polygon", "surface", "contour"])
    def test_plot_3d(self, rd_es, kind):
        fig, ax = rd_es.plot3d(
            save=False,
            show=True,
            axis_off=False,
            kind=kind,
            cmap="Reds",
            fig_width=4,
            fig_height=4,
            edgecolors="grey",
            linewidths=0.01,
        )

    def test_plot_2d(self, rd_es):
        fig, ax = rd_es.plot2d(
            axis_off=False,
            cmap="Reds",
            fig_height=2,
            fig_width=6,
            show=True,
            save=False,
            filename=rd_es.name + "_heatmap",
        )

    def test_discretize(self, rd_es):
        res = rd_es.discretize_tsam(noTypicalPeriods=1)
        assert_almost_equal(res.sum(), 2.118713381170598, decimal=3)
        rd_es.discretize_tsam(noTypicalPeriods=1, inplace=True)
        assert_almost_equal(res.sum(), 2.118713381170598, decimal=3)
        # check that the type is maintained
        assert type(rd_es) == EnergySeries

class TestEnergyDataFrame:
    def test_from_report_data(self, rd_edf):
        assert_almost_equal(rd_edf.sum().sum(), 397714.38260559476, decimal=3)
        assert rd_edf.units["CORE"] == settings.unit_registry.C

    def test_discretize(self, rd_edf):
        rd_edf = rd_edf.discretize_tsam(noTypicalPeriods=1)
        assert hasattr(rd_edf, "agg")
        rd_edf.discretize_tsam(noTypicalPeriods=1, inplace=True)
        assert hasattr(rd_edf, "agg")
        # check that the type is maintained
        assert type(rd_edf) == EnergyDataFrame

    def test_plot_2d(self, rd_edf):
        fig, ax = rd_edf.plot2d(
            axis_off=False,
            cmap="Reds",
            fig_height=None,
            fig_width=6,
            show=True,
            save=False,
            extent="tight",
        )

max_stars_repo_name: mingzehuang/archetypal

max_issues_repo_name: mingzehuang/archetypal

max_forks_repo_name: mingzehuang/archetypal

reason: from numpy


text:
# -*- coding: utf-8 -*-
import femm
from numpy import pi

from ...Functions.FEMM.comp_FEMM_Jcus import comp_FEMM_Jcus
from ...Functions.FEMM.set_FEMM_circuit_prop import set_FEMM_circuit_prop
from ...Functions.FEMM.set_FEMM_wind_material import set_FEMM_wind_material

def update_FEMM_simulation(
    output, materials, circuits, is_mmfs, is_mmfr, j_t0, is_sliding_band
):
    """Update the simulation by changing the rotor position and
    updating the currents

    Parameters
    ----------
    output :
        Output object
    """
    angle_rotor = output.get_angle_rotor()

    if is_sliding_band:  # No rotation without sliding band.
        # Rotor rotation using sliding band
        if output.simu.machine.rotor.is_internal:
            femm.mi_modifyboundprop("bc_ag2", 10, 180 * angle_rotor[j_t0] / pi)
        else:
            femm.mi_modifyboundprop("bc_ag2", 11, 180 * angle_rotor[j_t0] / pi)
        # Update currents
        for label in circuits:
            if "Circs" in label:  # Stator
                set_FEMM_circuit_prop(
                    circuits,
                    label,
                    output.elec.Is,
                    is_mmfs,
                    output.simu.machine.stator.winding.Npcpp,
                    j_t0,
                )
            if "Circr" in label:  # Rotor
                set_FEMM_circuit_prop(
                    circuits,
                    label,
                    output.elec.Ir,
                    is_mmfr,
                    output.simu.machine.stator.winding.Npcpp,
                    j_t0,
                )
        # Update winding materials
        for mat in materials:
            if "Js" in mat:  # Stator winding
                Jcus = comp_FEMM_Jcus(
                    output.simu.machine.stator, mat, output.elec.Is, j_t0, is_mmfs
                )
                materials = set_FEMM_wind_material(materials, mat, Jcus)
            elif "Jr" in mat:  # Rotor winding
                Jcus = comp_FEMM_Jcus(
                    output.simu.machine.rotor, mat, output.elec.Ir, j_t0, is_mmfr
                )
                materials = set_FEMM_wind_material(materials, mat, Jcus)

max_stars_repo_name: helene-t/pyleecan

max_issues_repo_name: helene-t/pyleecan

max_forks_repo_name: helene-t/pyleecan