text:
#!/usr/bin/env python
#
# reorder the BEAST results to be in ra/dec bins
#
# optimized BEAST runs are done by sets of stars in source density bins
#   sorted by flux and subdivided into smaller files
#   this allow for the BEAST grid to be cut and speeds up the fitting
#
# but this is non-ideal for almost any analysis of the results
#   especially the MegaBEAST!
#
# History: based on work by Heddy Arab prior to 2016
#          reordered code and flow to optimize for speed by Karl Gordon (09/16)

import os
import glob
import math

import h5py
from tqdm import trange, tqdm

import argparse
import numpy as np

from astropy import wcs
from astropy.io import fits
from astropy.table import Table

def reorder_beast_results_spatial(bricknum=None,
                                      stats_filename=None,
                                      region_filebase=None,
                                      output_filebase=None,
                                      reg_size=10.):
    """
    Do the spatial reordering of BEAST results.

    Two options for inputs: either the PHAT brick number (which will
    automatically generate file bases) or specify all of the file bases

    Parameters
    ----------
    bricknum : int or string
        PHAT brick num shortcut (supersedes other inputs)

    stats_filename : string
        Filename of the full run stats, from merge_stats_file.py

    region_filebase : string
        Filebase of the BEAST run sub-files

    output_filebase : string
        Path+prefix to use for output folders/files

    reg_size : float (default=10)
        spatial region size [arcsec]

    """

    if bricknum is not None:
        brick = str(bricknum)
        cat_filename = '/astro/dust_kg2/harab/toothpick_results/v1_1/b' + \
            brick + '_stats_v1_1.fits'
        reg_filebase = '/astro/dust_kg2/kgordon/BEAST_production/b' + \
            brick + '/b' + brick
        out_dir = '/astro/dust_kg2/kgordon/BEAST_production/b' + \
            brick + '/spatial'
        out_filebase = out_dir + '/b' + brick
        if not os.path.exists(out_dir):
            os.makedirs(out_dir)
    elif ((stats_filename is not None) and (region_filebase is not None) and (output_filebase is not None)):
        cat_filename = stats_filename
        reg_filebase = region_filebase
        out_filebase = output_filebase
    else:
        raise ValueError("Must specify either all file information or a PHAT brick number")

    # read in the full brick catalog and setup the spatial subdivided regions
    wcs_info, n_x, n_y = setup_spatial_regions(cat_filename,
                                               pix_size=reg_size)

    # setup array to store number of stars per pixel
    wcs_nstars = np.zeros((n_y, n_x), dtype=int)

    # find all the subdivided BEAST files for this brick
    sub_files = glob.glob(reg_filebase + '*_stats.fits')

    # loop over the files and output to the appropriate spatial region files
    # in loop:
    #      read in the locations of each star and calculated spatial region
    #      append the 1D pdfs
    #      append the sparse nD pdfs
    #      append the completeness function (??)

    #sub_files = sub_files[0:2]
    for cur_file in tqdm(sub_files, desc="orig sub files"):
        # read in the stats file
        cur_cat = Table.read(cur_file)
        n_objs = len(cur_cat)

        # read in the pdf1d file
        #   a number of extensions, one for each entry in the stats file
        hdulist = fits.open(cur_file.replace('_stats.fits','_pdf1d.fits'))
        cur_pdf1d_vals = []
        cur_pdf1d_name = []
        n_qnames = len(hdulist) - 1
        for k in range(n_qnames):
            cur_pdf1d_name.append(hdulist[k+1].header['EXTNAME'])
            cur_pdf1d_vals.append(hdulist[k+1].data)
        hdulist.close()

        # open the lnp file for reading
        cur_lnpfile = h5py.File(cur_file.replace('_stats.fits',
                                                 '_lnp.hd5'), 'r')

        # get the source density and subregion tag
        # allows for unique filenames for the spatial regions for output
        orig_reg_tag = cur_file[cur_file.find('_sd'):cur_file.find('_stats')]

        # determine the subregions for all the objects
        xy_vals = regions_for_objects(cur_cat['RA'],
                                      cur_cat['DEC'],
                                      wcs_info)

        # get the unique xy regions
        xy_names = np.squeeze(xy_vals['name'])
        uniq_xy_names, rindxs = np.unique(xy_names,
                                          return_inverse=True)

        #uniq_xy_names = uniq_xy_names[0:5]

        # loop over the unique xy regions
        for k in trange(len(uniq_xy_names), desc="outputing " + orig_reg_tag,
                        leave=False):
            uxy_name = uniq_xy_names[k]

            # get the indexes for the objects in this region
            indxs, = np.where(rindxs == k)

            # add the number of stars found to summary array
            wcs_nstars[xy_vals['y'][indxs[0]],
                       xy_vals['x'][indxs[0]]] += len(indxs)

            # create region directory if it does not exist
            reg_dir = out_filebase + '_' + uxy_name
            if not os.path.exists(reg_dir):
                os.makedirs(reg_dir)

            # base filename for output
            reg_filebase = out_filebase + '_' + uxy_name + '/' + \
                uxy_name + orig_reg_tag

            # write the stats info
            reg_stats_file = reg_filebase + '_stats.fits'
            cur_cat[indxs].write(reg_stats_file, overwrite=True)

            # write the pdf1d info
            reg_pdf1d_file = reg_filebase + '_pdf1d.fits'

            # setup the primary header and hdulist
            hdulist = fits.HDUList([fits.PrimaryHDU()])

            # generate the extensions
            for k, qname in enumerate(cur_pdf1d_name):
                # get the 1D PDFs for the cur objects
                #   plus the last column giving the values of the bins
                cur_reg_pdf1d = cur_pdf1d_vals[k][np.append(indxs,n_objs),:]

                chdu = fits.PrimaryHDU(cur_reg_pdf1d)
                chdu.header.set('XTENSION','IMAGE')
                chdu.header.set('EXTNAME',qname)

                hdulist.append(chdu)

            # write the 1D PDFs
            hdulist.writeto(reg_pdf1d_file, overwrite=True)

            # write the nD sparse likelihood info
            reg_lnp_file = reg_filebase + '_lnp.hd5'

            # open the file (overwrites an existing file)
            reg_lnpfile = h5py.File(reg_lnp_file, 'w')

            for i, k in enumerate(indxs):
                star_group = reg_lnpfile.create_group('star_%d' % i)

                # transfer the information
                for cp_name, cp_value in cur_lnpfile['star_%d' % k].items():
                    star_group.create_dataset(cp_name,data=cp_value.value)

            reg_lnpfile.close()

        cur_lnpfile.close()

        # Now, write out the WCS info and number of stars per pixel to file
        #   do every subregion file to have an on-the-fly check
        header = wcs_info.to_header()
        hdu = fits.PrimaryHDU(wcs_nstars, header=header)

        # Save to FITS file
        hdu.writeto(out_filebase+'_nstars.fits', overwrite=True)

def setup_spatial_regions(cat_filename,
                          pix_size=10.0):
    """
    The spatial regions are setup via a WCS object

    Parameters
    ----------
    cat_filename : string
       filename of catalog

    pix_size : float
       size of pixels/regions in arcsec

    Returns
    -------
    wcs_info: astropy WCS object
    """

    # read in the catalog file
    cat = Table.read(cat_filename)

    # min/max ra
    min_ra = cat['RA'].min()
    max_ra = cat['RA'].max()
    min_dec = cat['DEC'].min()
    max_dec = cat['DEC'].max()

    # ra/dec delta values
    dec_delt = pix_size/3600.
    ra_delt = dec_delt

    # compute the number of pixels and
    n_y = int(np.rint((max_dec - min_dec)/dec_delt) + 1)
    n_x = int(np.rint(math.cos(0.5*(max_dec+min_dec)*math.pi/180.)*
                      (max_ra-min_ra)/ra_delt) + 1)

    # ra delta should be negative
    ra_delt *= -1.

    print('# of x & y pixels = ', n_x, n_y)

    w = wcs.WCS(naxis = 2)
    w.wcs.crpix = np.asarray([n_x, n_y], dtype = float) / 2.
    w.wcs.cdelt = [ra_delt, dec_delt]
    w.wcs.crval = np.asarray([(min_ra+max_ra), (min_dec+max_dec)]) / 2.
    w.wcs.ctype = ["RA---TAN", "DEC--TAN"]

    return (w, n_x, n_y)

def regions_for_objects(ra,
                        dec,
                        wcs_info):
    """
    Generate the x,y coordinates for each object based on the input
    ra/dec and already created WCS information.

    Parameters
    ----------
    ra : array of float
       right ascension of the objects

    dec : array of float
       declination of the objects

    wcs_info: astropy WCS object
       previously generated WCS object based on the full catalog

    Returns
    -------
    dictonary of:

    x : int array
      x values of regions

    y : int array
      y values of regions

    name : str array
      string array composed of x_y
    """

    # generate the array needed for fast conversion
    world = np.empty((len(ra),2),float)
    world[:,0] = ra
    world[:,1] = dec

    # convert
    pixcrd = wcs_info.wcs_world2pix(world, 1)

    # get the arrays to return
    x = pixcrd[:,0].astype(int)
    y = pixcrd[:,1].astype(int)
    xy_name = [None]*len(ra)

    x_str = x.astype(np.string_)
    y_str = y.astype(np.string_)

    for k in range(len(x)):
        xy_name[k] = str(x[k]) + '_' + str(y[k])

    # return the results as a dictonary
    #   values are truncated to provide the ids for the subregions
    return {'x': x, 'y': y, 'name': xy_name}

if __name__ == '__main__':

    # command line params to specify the run directory
    #   and any other needed parameters

    # commandline parser
    parser = argparse.ArgumentParser()
    parser.add_argument("-b","--bricknum", default=None,
                        help="PHAT brick num shortcut" + \
                        " (supersedes other inputs)")
    parser.add_argument("-s","--stats_filename", default=None,
                        help="Filename of the full run stats")
    parser.add_argument("-r","--region_filebase", default=None,
                        help="Filebase of the run regions")
    parser.add_argument("-o","--output_filebase", default=None,
                        help="Filebase to use for output")
    parser.add_argument("-p","--reg_size", default=10., type=float,
                        help="spatial region size [arcsec]")
    args = parser.parse_args()

    reorder_beast_results_spatial(bricknum=args.bricknum,
                                      stats_filename=args.stats_filename,
                                      region_filebase=args.region_filebase,
                                      output_filebase=args.output_filebase,
                                      reg_size=args.reg_size)

max_stars_repo_name: marthaboyer/beast

max_issues_repo_name: marthaboyer/beast

max_forks_repo_name: marthaboyer/beast

reason: import numpy,from astropy


text:
# coding=utf-8
from __future__ import absolute_import, division, print_function
import pytest
from astwro.phot.phot_error import PhotError

from astwro.phot import phot_error
import numpy as np

@pytest.fixture
def mag_err():
    import os
    import pickle
    pickelfile = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'mag_err.npz')
    dic = np.load(pickelfile)
    return dic

def test_ph_error_tofit(mag_err):
    mag, err = mag_err['mag'], mag_err['mag_err']
    pe = PhotError(mag, err)
    assert pe.err_fitted.shape[0] > 0
    assert pe.err_means.shape[0] > 0
    assert pe.mag_means.shape[0] > 0
    assert len(pe.outlayers_mask) == len(mag)

def test_ph_error_custom_fn(mag_err):
    mag, err = mag_err['mag'], mag_err['mag_err']
    zzz = lambda x: np.percentile(x, 20) if len(x) > 0 else np.nan
    pe = PhotError(mag, err, meanfn=zzz, fitting_order=6, weighted_fit=True)
    assert pe.err_fitted.shape[0] > 0
    assert pe.err_means.shape[0] > 0
    assert pe.mag_means.shape[0] > 0
    assert len(pe.outlayers_mask) == len(mag)

ext: py

max_stars_repo_name: majkelx/astwro

max_issues_repo_name: majkelx/astwro

max_forks_repo_name: majkelx/astwro

reason: import numpy


text:
# pylint: disable=E0401

"""
implements FlexCloud which facilitates point cloud classification by tracking associations between
points, features and labels.

stylish future upgrade for FlexCloud: subclass MutableMapping to index with keys like a dictionary.
make each asset a key. use bracket notation with the facet name (__getitem__) to return an asset
array.
"""

import numpy as np

class FlexCloud(object):
    """
    given a 3d point cloud as a 2d numpy array, shift its points close to the origin and track its
    features. any supplemental information must be added separately, after the object has been
    instantiated with its points. supplemental information is stored as "assets" which can be
    1d or 2d arrays of floats or integers.

    the asset index dictionary looks like the following:
    {
        intensity: {
            index: index_array,
            asset: 1d array of floats,
            meta: "this intensity data is not calibrated"
        },
        geometry_mso_1: {
            index: index_array,
            asset: 2d array of floats,
            meta: {
                voxel: 0.05,
                scales: [0.15, 0.2, 0.25]
            }
        },
        known_label: {
            index: index_array,
            asset: 1d array of ints,
            meta: None
        },
        multilabel_predicted: {
            index: index_array,
            asset: 2d array of ints,
            meta: "this is why i think this is a good idea..."
        }
    }
    """

    #==================================

    def __init__(self, input_cloud):

        if input_cloud.ndim != 2:
            raise ValueError("input point cloud must be a 2D array")
        if input_cloud.shape[1] != 3:
            raise ValueError("must be initialized with a 3D point cloud")
        # now bring the point cloud in to the origin
        self.corner = input_cloud[0]
        self.points = input_cloud - self.corner
        # count how many points we have in the original point cloud
        self.num_points = input_cloud.shape[0]
        self.id_index = np.arange(self.num_points)
        # initialize the asset dictionary
        self.assets = {}

    #==================================

    def add_asset(self, asset_array, index_array, asset_name, meta=None):
        """
        add a new asset array to the cloud's asset index. input index array does not need to be
        sorted or unique, but it will be stored sorted and unique to simplify the set ops later.
        """

        # first make sure this is a good idea
        if asset_name in self.assets:
            raise ValueError("asset {} already exists in asset dictionary".format(asset_name))
        asset_array, index_array = self._validate_asset(asset_array, index_array)

        # assemble the asset
        asset = {
            "asset": asset_array,
            "index": index_array,
            "meta": meta
        }

        self.assets[asset_name] = asset

    #==================================

    def _validate_asset(self, asset_array, index_array):
        """
        unique and sort the index array, and align the asset array to match it.
        raise an exception if the asset to be added to the asset dictionary won't fit in the cloud
        """

        # make sure the asset has no more than 2 dimensions-- we don't support 2d assets per
        # point (yet?)
        if asset_array.ndim > 2:
            raise ValueError("asset array has too many dimensions")
        # make sure the asset and its index array are aligned
        if asset_array.shape[0] != index_array.size:
            raise ValueError("asset and index arrays misaligned")
        # make sure all indices are unique and sorted
        unique_indices, index_to_unique = np.unique(index_array, return_index=True)
        # make sure the index array will index into the cloud
        if index_array.min() < 0 or index_array.max() >= self.num_points:
            raise ValueError("index array addresses outside the extant cloud")

        # now return (assets, indices)
        return asset_array.take(index_to_unique, axis=0), unique_indices

    #==================================

    def intersection(self, asset_names):
        """
        given a list of names of assets, compute the intersection of their index sets and return
        that, along with the horizontal concatenation of all the corresponding assets.
        """

        # fold down the sequence of asset names, accumulating to the identity index set using the
        # intersection operator.
        index_accumulator = self.id_index
        for name in asset_names:
            this_index = self.assets[name]["index"]
            index_accumulator = np.intersect1d(index_accumulator, this_index, assume_unique=True)

        # how many points are there?
        num_points = index_accumulator.size

        # we can now use in1d to find each asset that is present in the intersection
        asset_accumulator = []
        for name in asset_names:
            this_index = self.assets[name]["index"]
            this_asset = self.assets[name]["asset"]
            # find which assets are present in the output index set
            mask = np.in1d(this_index, index_accumulator, assume_unique=True)
            # put them on the accumulator
            asset_accumulator.append(np.compress(mask, this_asset, axis=0).reshape(num_points, -1))

        return_assets = np.concatenate(asset_accumulator, axis=1)

        return index_accumulator, return_assets

    #==================================

    def take(self, index_array=None, original_coordinates=True):
        """
        equivalent to ndarray.take(). return a subset of the FlexCloud's points addressed by an
        index array, in the original coordinates if desired. if no index given, return all.
        """
        if original_coordinates:
            return_points = self.points + self.corner
        else:
            return_points = self.points
        if index_array is not None:
            return return_points.take(index_array, axis=0)
        else:
            return return_points

    #==================================

    #==================================

max_stars_repo_name: grayhem/nimrud

max_issues_repo_name: grayhem/nimrud

max_forks_repo_name: grayhem/nimrud


text:
import heterocl as hcl
import numpy as np
import numpy.testing as tst
import os
import urllib.request
from hlib.op.extern import (
    create_extern_module, register_extern_ip,
    register_tensors, include_dependency)

@register_extern_ip(type="vhls")
def toynn_vhls_ip(input_1, output_1, name=None):
    if name is None: name = "myproject"
    # Function behavior definition
    with hcl.Stage("ExternModule.toyNN") as Module:
        register_tensors([input_1, output_1])

    Module.ext_ip_name = name
    Module.inputs = [input_1, output_1]

    # Include cpp/hpp files
    if not os.path.exists("firmware"):
       urllib.request.urlretrieve("https://raw.githubusercontent.com/Hecmay/debug.trace/main/toynn.tar.gz", filename="toynn.tar.gz")
       os.system("tar -zxvf toynn.tar.gz")

    source = [
        "firmware/myproject.cpp",
        "firmware/nnet_utils/",
        "firmware/weights/"
    ]
    Module.source = include_dependency(source)
    create_extern_module(Module, ip_type="HLS")

def test_toy_nn():
    dtype = hcl.Float(32)
    hcl.init(dtype)

    input_1 = hcl.placeholder((16,), dtype=dtype, name="input_1")
    output_1 = hcl.placeholder((5,), dtype=dtype, name="output_1")
    def math_func(input_1, output_1):
        toynn_vhls_ip(input_1, output_1)

    target = hcl.Platform.aws_f1
    s = hcl.create_schedule([input_1, output_1], math_func)
    s.to(input_1, target.xcel)
    s.to(output_1, target.host)

    target.config(compiler="vitis", mode="debug")
    code = hcl.build(s, target)
    assert "nnet::softmax" in code, code
    os.system("rm -rf firmware toynn.tar.gz")

if __name__ == "__main__":
    test_toy_nn()

max_stars_repo_name: zzzDavid/heterocl

max_issues_repo_name: zzzDavid/heterocl


text:
import numpy as np

# Compute a moving mean along the x-axis.
def moving_mean(x, y, yerr=None, window=3.):

    # Set the weights.
    if yerr is None:
        weights = np.ones_like(y)
    else:
        weights = 1/yerr**2

    # Sums for computing the mean.
    sum1 = np.append(0, np.cumsum(weights*y))
    sum2 = np.append(0, np.cumsum(weights))

    # Indices at the start and end of the window.
    i = np.searchsorted(x, x - window/2.)
    j = np.searchsorted(x, x + window/2.)

    # Compute the mean.
    mean = (sum1[j] - sum1[i])/(sum2[j] - sum2[i])

    return mean

# Fit a curve for local linear
def linfit(lstidx, x, y, sky, mag, emag):

    sort = np.argsort(lstidx)
    invsort = np.argsort(sort)

    lstidx = lstidx[sort]
    x = x[sort]
    y = y[sort]
    sky = sky[sort]
    mag = mag[sort]
    emag = emag[sort]

    _, idx = np.unique(lstidx, return_inverse=True)

    nobs = np.bincount(idx)
    strides = np.append(0, np.cumsum(nobs))

    xbar = np.bincount(idx, x)/np.bincount(idx)
    ybar = np.bincount(idx, y)/np.bincount(idx)

    mat = np.column_stack([np.ones(len(lstidx)), x-xbar[idx], y-ybar[idx], sky])

    pars = np.zeros((len(nobs), 4))
    pars[:,0] = np.bincount(idx, mag/emag**2)/np.bincount(idx, 1/emag**2)

    for i in range(len(nobs)):

        if nobs[i] < 5:
             continue

        i1 = strides[i]
        i2 = strides[i+1]

        pars[i] = np.linalg.lstsq(mat[i1:i2]/emag[i1:i2,None], mag[i1:i2]/emag[i1:i2], rcond=None)[0]

    fit = np.sum(pars[idx]*mat, axis=1)

    fit1 = pars[idx,0]
    fit2 = fit - pars[idx,0]

    return fit1[invsort], fit2[invsort], (nobs > 4)[idx][invsort]

def local_lin(jd, lstseq, x, y, sky, mag, emag, window=5., maxiter=50, dtol=1e-3):

    lstidx = (lstseq % 270)
    fit0 = np.zeros(len(jd))
    fit1 = np.zeros(len(jd))

    fit = np.zeros_like(mag)
    for niter in range(maxiter):
        fit1, fit2, mask = linfit(lstidx, x, y, sky, mag - fit0, emag)
        fit0 = moving_mean(jd, mag - fit1 - fit2, emag, window)

        if niter > 0:

            if np.all(np.abs(fit - fit0 - fit1 - fit2) < dtol):
                break

        fit = fit0 + fit1 + fit2

    return fit#, fit0, fit1, fit2, mask


text:
from __future__ import print_function
from six.moves import range
from PIL import Image

import torch.backends.cudnn as cudnn
import torch
import torch.nn as nn
from torch.autograd import Variable
import torch.optim as optim
import os
import time

import numpy as np
import torchfile
from sklearn.metrics import confusion_matrix

from miscc.config import cfg
from miscc.utils import mkdir_p
from miscc.utils import weights_init
from miscc.utils import save_img_results, save_model
from miscc.utils import KL_loss
from miscc.utils import compute_discriminator_loss, compute_generator_loss

# from tensorboard import summary
import tensorflow as tf
# from tensorflow.summary import FileWriter

from torchsummary import summary

class GANTrainer(object):
    def __init__(self, output_dir):
        if cfg.TRAIN.FLAG:
            self.model_dir = os.path.join(output_dir, 'Model')
            self.image_dir = os.path.join(output_dir, 'Image')
            self.eval_dir = os.path.join(self.image_dir, 'Eval')
            self.log_dir = os.path.join(output_dir, 'Log')
            mkdir_p(self.model_dir)
            mkdir_p(self.image_dir)
            mkdir_p(self.eval_dir)
            mkdir_p(self.log_dir)
            self.summary_writer = tf.summary.FileWriter(self.log_dir)

        self.max_epoch = cfg.TRAIN.MAX_EPOCH
        self.snapshot_interval = cfg.TRAIN.SNAPSHOT_INTERVAL
        s_gpus = cfg.GPU_ID.split(',')
        self.gpus = [int(ix) for ix in s_gpus]
        self.num_gpus = len(self.gpus)
        self.n_output = cfg.GAN.N_OUTPUT
        if cfg.CUDA:
            self.batch_size = cfg.TRAIN.BATCH_SIZE * self.num_gpus
            torch.cuda.set_device(self.gpus[0])
            cudnn.benchmark = True
        else:
            self.batch_size = cfg.TRAIN.BATCH_SIZE
            self.num_gpus = 0

    def load_network_embedding(self):
        # from model import EmbeddingNet
        # netE = EmbeddingNet(cfg.AUDIO.FEATURE_DIM, cfg.AUDIO.DIMENSION)
        # netE.apply(weights_init)
        # if cfg.EMB_NET != '':
        #     state_dict = torch.load(cfg.EMB_NET,
        #                             map_location=lambda storage, loc: storage)
        #     netE.load_state_dict(state_dict)
        #     removed = list(list(netE.children())[0].children())[:-1]
        #     netE = nn.Sequential(*removed)
        print('Load from: ', cfg.EMB_NET)
        netE = torch.load(cfg.EMB_NET)
        for p in netE.parameters():
            p.requires_grad=False
        if cfg.CUDA:
            netE.cuda()
        return netE

    def load_network_stageI(self):
        from model import STAGE1_G, STAGE1_D
        netG = STAGE1_G()
        netG.apply(weights_init)
        netD = STAGE1_D()
        netD.apply(weights_init)

        if cfg.NET_G != '':
            state_dict = \
                torch.load(cfg.NET_G,
                           map_location=lambda storage, loc: storage)
            netG.load_state_dict(state_dict)
            print('Load from: ', cfg.NET_G)
        if cfg.NET_D != '':
            state_dict = \
                torch.load(cfg.NET_D,
                           map_location=lambda storage, loc: storage)
            netD.load_state_dict(state_dict)
            print('Load from: ', cfg.NET_D)
        if cfg.CUDA:
            netG.cuda()
            netD.cuda()
        return netG, netD

    # ############# For training stageII GAN  #############
    def load_network_stageII(self):
        from model import STAGE1_G, STAGE2_G_twostream, STAGE2_D

        Stage1_G = STAGE1_G()
        netG = STAGE2_G_twostream(Stage1_G)
        netG.apply(weights_init)
        # print(netG)
        if cfg.NET_G != '':
            state_dict = \
                torch.load(cfg.NET_G,
                           map_location=lambda storage, loc: storage)
            netG.load_state_dict(state_dict)
            print('Load from: ', cfg.NET_G)
        elif cfg.STAGE1_G != '':
            state_dict = \
                torch.load(cfg.STAGE1_G,
                           map_location=lambda storage, loc: storage)
            netG.STAGE1_G.load_state_dict(state_dict)
            print('Load from: ', cfg.STAGE1_G)
        else:
            print("Please give the Stage1_G path")
            return

        netD = STAGE2_D()
        netD.apply(weights_init)
        if cfg.NET_D != '':
            state_dict = \
                torch.load(cfg.NET_D,
                           map_location=lambda storage, loc: storage)
            netD.load_state_dict(state_dict)
            print('Load from: ', cfg.NET_D)
        # print(netD)

        if cfg.CUDA:
            netG.cuda()
            netD.cuda()
        return netG, netD

    def train(self, data_loader, stage=1, eval_dataset=None):
        if stage == 1:
            netG, netD = self.load_network_stageI()
        else:
            netG, netD = self.load_network_stageII()

        if cfg.DATASET_NAME == 'audioset':
            # print('building')
            netE = self.load_network_embedding()
            # print('done building')
            # summary(netE, input_size=(128, 430))
        else:
            netE = None

        nz = cfg.Z_DIM
        batch_size = self.batch_size
        fixed_noise = torch.FloatTensor(batch_size, nz).normal_(0, 1)
        if self.n_output == 1:
            fake_labels = Variable(torch.FloatTensor(batch_size).fill_(0))
        elif self.n_output > 1:
            fake_labels = Variable(torch.FloatTensor(batch_size).fill_(self.n_output))
        else:
            raise Exception('n_output should be an integer larger than 0.')

        if cfg.CUDA:
            fixed_noise = fixed_noise.cuda()
            fake_labels = fake_labels.cuda()

        generator_lr = cfg.TRAIN.GENERATOR_LR
        discriminator_lr = cfg.TRAIN.DISCRIMINATOR_LR
        lr_decay_step = cfg.TRAIN.LR_DECAY_EPOCH
        optimizerD = \
            optim.Adam(netD.parameters(),
                       lr=cfg.TRAIN.DISCRIMINATOR_LR, betas=(0.9, 0.999))
        netG_para = []
        for p in netG.parameters():
            if p.requires_grad:
                netG_para.append(p)
        optimizerG = optim.Adam(netG_para,
                                lr=cfg.TRAIN.GENERATOR_LR,
                                betas=(0.9, 0.999))
        count = 0
        for epoch in range(self.max_epoch):
            start_t = time.time()
            print("running epoch {}".format(epoch))
            if epoch % lr_decay_step == 0 and epoch > 0:
                generator_lr *= 0.5
                for param_group in optimizerG.param_groups:
                    param_group['lr'] = generator_lr
                discriminator_lr *= 0.5
                for param_group in optimizerD.param_groups:
                    param_group['lr'] = discriminator_lr

            for i, data in enumerate(data_loader, 0):
                embedding, real_labels, real_imgs = self.get_embedding(data, netE)
                lr_fake, fake_imgs, mu, logvar = self.gen_image(embedding, netG)

                ############################
                # (3) Update D network
                ###########################
                if self.n_output == 1:
                    loss_func = nn.BCELoss
                else:
                    loss_func = nn.CrossEntropyLoss

                netD.zero_grad()
                errD, errD_real, errD_wrong, errD_fake = \
                    compute_discriminator_loss(netD, real_imgs, fake_imgs,
                                               real_labels, fake_labels,
                                               mu, loss_func, self.gpus)
                errD.backward()
                optimizerD.step()
                ############################
                # (2) Update G network
                ###########################
                netG.zero_grad()
                errG = compute_generator_loss(netD, fake_imgs,
                                              real_labels, mu, loss_func, self.gpus)
                kl_loss = KL_loss(mu, logvar)
                errG_total = errG + kl_loss * cfg.TRAIN.COEFF.KL
                errG_total.backward()
                optimizerG.step()

                count = count + 1
                # if (epoch+1) % 100 == 0:
                if i % 50 == 0:
                    # summary_D = tf.summary.scalar('D_loss', errD.data[0])
                    # summary_D_r = tf.summary.scalar('D_loss_real', errD_real)
                    # summary_D_w = tf.summary.scalar('D_loss_wrong', errD_wrong)
                    # summary_D_f = tf.summary.scalar('D_loss_fake', errD_fake)
                    # summary_G = tf.summary.scalar('G_loss', errG.data[0])
                    # summary_KL = tf.summary.scalar('KL_loss', kl_loss.data[0])
                    #
                    # self.summary_writer.add_summary(summary_D, count)
                    # self.summary_writer.add_summary(summary_D_r, count)
                    # self.summary_writer.add_summary(summary_D_w, count)
                    # self.summary_writer.add_summary(summary_D_f, count)
                    # self.summary_writer.add_summary(summary_G, count)
                    # self.summary_writer.add_summary(summary_KL, count)

                    # save the image result for each epoch
                    save_img_results(real_imgs, fake_imgs, epoch, self.image_dir)
                    if lr_fake is not None:
                        save_img_results(None, lr_fake, epoch, self.image_dir)
            end_t = time.time()
            print('''[%d/%d][%d/%d] Loss_D: %.4f Loss_G: %.4f Loss_KL: %.4f
                     Loss_real: %.4f Loss_wrong:%.4f Loss_fake %.4f
                     Total Time: %.2fsec
                  '''
                  % (epoch+1, self.max_epoch, i, len(data_loader),
                     errD.data[0], errG.data[0], kl_loss.data[0],
                     errD_real, errD_wrong, errD_fake, (end_t - start_t)))

            if eval_dataset is not None:
                self.evaluate(netG, netE, eval_dataset, epoch, stage=stage, random=False)
                self.evaluate(netG, netE, eval_dataset, epoch, stage=stage, random=True)

            if epoch % self.snapshot_interval == 0:
                save_model(netG, netD, epoch, self.model_dir)
        #
        save_model(netG, netD, self.max_epoch, self.model_dir)
        #
        self.summary_writer.close()

    def evaluate(self, netG, netE, dataset, epoch, stage=1, random=False):
        dataloader = torch.utils.data.DataLoader(
            dataset, batch_size=cfg.TRAIN.BATCH_SIZE,
            shuffle=random, num_workers=1)

        ev = 2 if random else 1
        for i, data in enumerate(dataloader, 0):
            embedding, real_labels, real_imgs = self.get_embedding(data, netE)
            lr_fake, fake_imgs, mu, logvar = self.gen_image(embedding, netG)
            save_img_results(real_imgs, fake_imgs, epoch, self.eval_dir, ev)
            break
        return

    def get_embedding(self, data, netE=None):
        ######################################################
        # (1) Prepare training data
        ######################################################
        if cfg.DATASET_NAME == 'audioset':
            if cfg.CUDA:
                audio_feat = data[0].to(device=torch.device('cuda'), dtype=torch.float)
                real_imgs = data[1].to(device=torch.device('cuda'), dtype=torch.float)
                if self.n_output > 1:
                    real_labels = data[2].to(device=torch.device('cuda'), dtype=torch.float)
                else:
                    real_labels = Variable(torch.FloatTensor(self.batch_size).fill_(1)).cuda()
            else:
                audio_feat = data[0].type(torch.FloatTensor)
                real_imgs = data[1].type(torch.FloatTensor)
                if self.n_output > 1:
                    real_labels = data[2].type(torch.FloatTensor)
                else:
                    real_labels = Variable(torch.FloatTensor(self.batch_size).fill_(1))
            embedding = netE(audio_feat)
            # with torch.no_grad():
            #     m = list(netE._modules.values())[0]
            #     ft = audio_feat
            #     ftlist = []
            #     for j, module in enumerate(m):
            #         nft = module(ft)
            #         ftlist.append(nft)
            #         ft = nft
            #     embedding = ftlist[-2]
        else:
            real_img_cpu, embedding = data
            real_img_cpu = real_img_cpu.type(torch.FloatTensor)
            real_imgs = Variable(real_img_cpu)
            embedding = embedding.type(torch.FloatTensor)
            embedding = Variable(embedding)
        if cfg.CUDA:
            real_imgs = real_imgs.cuda()
            embedding = embedding.cuda()

        return embedding, real_labels, real_imgs

    def gen_image(self, embedding, netG):
        #######################################################
        # (2) Generate fake images
        ######################################################
        noise = Variable(torch.FloatTensor(self.batch_size, cfg.Z_DIM))
        noise.data.normal_(0, 1)
        if cfg.CUDA:
            noise = noise.cuda()
        inputs = (embedding, noise)
        if cfg.CPU:
            lr_fake, fake_imgs, mu, logvar = netG(*inputs)
        else:
            lr_fake, fake_imgs, mu, logvar = \
                nn.parallel.data_parallel(netG, inputs, self.gpus)

        return lr_fake, fake_imgs, mu, logvar

    def sample(self, datapath, stage=1):
        #if stage == 1:
        netG, _ = self.load_network_stageI()
        #else:
        #    netG, _ = self.load_network_stageII()
        #netG.eval()

        # Load text embeddings generated from the encoder
        if cfg.DATASET_NAME == 'coco':
            t_file = torchfile.load(datapath)
            captions_list = t_file.raw_txt
            embeddings = np.concatenate(t_file.fea_txt, axis=0)
            num_embeddings = len(captions_list)
            print('Successfully load sentences from: ', datapath)
            print('Total number of sentences:', num_embeddings)
            print('num_embeddings:', num_embeddings, embeddings.shape)

        # Load Audio files
        elif cfg.DATASET_NAME == 'audioset':
            netE = self.load_network_embedding()

        # path to save generated samples
        save_dir = cfg.NET_G[:cfg.NET_G.find('.pth')]
        mkdir_p(save_dir)

        batch_size = np.minimum(num_embeddings, self.batch_size)
        nz = cfg.Z_DIM
        noise = Variable(torch.FloatTensor(batch_size, nz))
        if cfg.CUDA:
            noise = noise.cuda()
        count = 0
        while count < num_embeddings:
            if count > 3000:
                break
            iend = count + batch_size
            if iend > num_embeddings:
                iend = num_embeddings
                count = num_embeddings - batch_size
            embeddings_batch = embeddings[count:iend]
            # captions_batch = captions_list[count:iend]
            embedding = Variable(torch.FloatTensor(embeddings_batch))
            if cfg.CUDA:
                embedding = embedding.cuda()

            #######################################################
            # (2) Generate fake images
            ######################################################
            noise.data.normal_(0, 1)
            inputs = (embedding, noise)
            if cfg.CPU:
                netG(embedding, noise)
            else:
                _, fake_imgs, mu, logvar = \
                    nn.parallel.data_parallel(netG, inputs, self.gpus)
            for i in range(batch_size):
                save_name = '%s/%d.png' % (save_dir, count + i)
                im = fake_imgs[i].data.cpu().numpy()
                im = (im + 1.0) * 127.5
                im = im.astype(np.uint8)
                # print('im', im.shape)
                im = np.transpose(im, (1, 2, 0))
                # print('im', im.shape)
                im = Image.fromarray(im)
                im.save(save_name)
            count += batch_size

class EmbeddingNetTrainer(object):
    def __init__(self, cfg, output_dir=None, model=None):
        self.output_dir = output_dir
        self.feature_dim = cfg.AUDIO.FEATURE_DIM
        self.output_dim = cfg.AUDIO.DIMENSION
        self.use_gpu = cfg.CUDA
        self.epochs = cfg.TRAIN.MAX_EPOCH
        self.n_workers = int(cfg.WORKERS)
        self.batch_size = cfg.TRAIN.BATCH_SIZE
        self.learning_rate = 8e-4
        self.num_classes = cfg.NUM_CLASSES
        if model is not None:
            self.model = model
        else:
            self.build_model()
        if self.use_gpu:
            self.model = self.model.to(device=torch.device('cuda'))

    def build_model(self):
        from model import EmbeddingNet
        self.embnet = EmbeddingNet(self.feature_dim, self.output_dim)
        self.model = nn.Sequential( self.embnet,
                                    nn.Linear(self.output_dim, self.output_dim*2),
                                    nn.ReLU(),
                                    nn.Linear(self.output_dim*2, self.num_classes),
                                    )
    def train(self, train_set, eval_set=None):
        dataloader = torch.utils.data.DataLoader(
            train_set, batch_size=self.batch_size,
            shuffle=True, num_workers=self.n_workers)
        optimizer = optim.Adam(self.model.parameters(), lr=self.learning_rate)
        best_eval_acc = -1.0
        for e in range(self.epochs):
            print('Epoch %d / %d' % (e, self.epochs))
            st = time.time()
            # Train
            num_correct = 0
            num_samples = 0
            for i, data in enumerate(dataloader):
                self.model = self.model.train()  # put model to training mode
                x = data['audio']
                y = data['label']
                if self.use_gpu:
                    x = x.to(device=torch.device('cuda'), dtype=torch.float)
                    y = y.to(device=torch.device('cuda'), dtype=torch.long)

                scores = self.model(x)
                scores = scores.view(-1, self.num_classes)
                loss = torch.nn.functional.cross_entropy(scores, y)
                optimizer.zero_grad()
                loss.backward()
                optimizer.step()

                with torch.no_grad():
                    _, preds = scores.max(1)
                    num_correct += (preds == y).sum()
                    num_samples += preds.size(0)

                if i % 200 == 0:
                    print('Epoch %d, Iteration %d, loss = %.4f' % (e, i, loss.item()))

            acc = float(num_correct) / num_samples if num_samples > 0 else 0.0
            print('Epoch %d took %.2f mins' % (e, (time.time()-st)/60.0))
            print('Train: Got %d / %d correct (%.2f)' % (num_correct, num_samples, 100 * acc))

            # Evaluate
            if eval_set is not None:
                acc = self.evaluate(eval_set)
                if acc >= best_eval_acc and self.output_dir is not None:
                    save_dir = os.path.join(self.output_dir, 'embnet')
                    mkdir_p(save_dir)
                    torch.save( self.embnet.state_dict(),
                                os.path.join(save_dir, 'embnet.pth')
                                )
                    print('Save Embedding Net model')
                    best_eval_acc = acc

    def evaluate(self, eval_set):
        dataloader = torch.utils.data.DataLoader(
            eval_set, batch_size=self.batch_size, num_workers=self.n_workers)
        num_correct = 0
        num_samples = 0
        model = self.model.eval()  # set model to evaluation mode
        with torch.no_grad():
            for i, data in enumerate(dataloader):
                x = data['audio']
                y = data['label']
                if self.use_gpu:
                    x = x.to(device=torch.device('cuda'), dtype=torch.float)  # move to device, e.g. GPU
                    y = y.to(device=torch.device('cuda'), dtype=torch.long)
                scores = model(x)
                scores = scores.view(-1, self.num_classes)
                loss = torch.nn.functional.cross_entropy(scores, y)
                _, preds = scores.max(1)
                num_correct += (preds == y).sum()
                num_samples += preds.size(0)
            acc = float(num_correct) / num_samples if num_samples > 0 else 0.0
            print('Eval:  Got %d / %d correct (%.2f), loss %.4f' % (num_correct, num_samples, 100 * acc, loss.item()))
        return acc

class EmbeddingNetLSTMTrainer(EmbeddingNetTrainer):
    def build_model(self):
        from model import EmbeddingNetLSTM
        self.model = nn.Sequential( EmbeddingNetLSTM(self.feature_dim, self.output_dim*2),
                                    nn.ReLU(),
                                    nn.Linear(self.output_dim*2, self.num_classes),
                                    )

def print_cm(cm, labels, hide_zeroes=False, hide_diagonal=False, hide_threshold=None):
    """pretty print for confusion matrixes"""
    columnwidth = max([len(x) for x in labels] + [5])  # 5 is value length
    empty_cell = " " * columnwidth
    # Print header
    print("    " + empty_cell, end=" ")
    for label in labels:
        print("%{0}s".format(columnwidth) % label, end=" ")
    print()
    # Print rows
    for i, label1 in enumerate(labels):
        print("    %{0}s".format(columnwidth) % label1, end=" ")
        for j in range(len(labels)):
            cell = "%{0}.1f".format(columnwidth) % cm[i, j]
            if hide_zeroes:
                cell = cell if float(cm[i, j]) != 0 else empty_cell
            if hide_diagonal:
                cell = cell if i != j else empty_cell
            if hide_threshold:
                cell = cell if cm[i, j] > hide_threshold else empty_cell
            print(cell, end=" ")
        print()


text:
"""
@author: Maziar Raissi
"""

from Multistep_NN import Multistep_NN

import numpy as np
from scipy.integrate import odeint

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from plotting import newfig, savefig
import matplotlib.gridspec as gridspec

import scipy.io

def colorline3d(ax, x, y, z, cmap):
    N = len(x)
    skip = int(0.01*N);
    for i in range(0,N,skip):
        ax.plot(x[i:i+skip], y[i:i+skip], z[i:i+skip], color=cmap(int(255*i/N)))

if __name__ == "__main__":

    # Load Data
    data = scipy.io.loadmat('./Cylinder.mat')

    X_star = data['X_star']

    # time points
    dt = 0.02
    t_star = np.arange(0,X_star.shape[0])*dt

    # initial condition
    x0 = X_star[0,:]

    skip = 1
    dt = t_star[skip] - t_star[0]
    X_train = X_star[0::skip,:]

    X_train = np.reshape(X_train, (1,X_train.shape[0],X_train.shape[1]))

    layers = [3, 256, 3]

    M = 1
    scheme = 'AM'
    model = Multistep_NN(dt, X_train, layers, M, scheme)

    N_Iter = 50000
    model.train(N_Iter)

    def learned_f(x,t):
        f = model.predict_f(x[None,:])
        return f.flatten()

    learned_X_star = odeint(learned_f, x0, t_star)

    ####### Plotting ##################
    fig, ax = newfig(1.0, 0.8)
    ax.axis('off')

    gs0 = gridspec.GridSpec(1, 2)
    gs0.update(top=0.95, bottom=0.1, left=0.0, right=0.90, wspace=0.15)

    ax = plt.subplot(gs0[:, 0:1], projection='3d')
    ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    colorline3d(ax, X_star[:,0], X_star[:,1], X_star[:,2], cmap = plt.cm.jet)
    ax.grid(False)
    ax.set_xlim([-200,200])
    ax.set_ylim([-200,200])
    ax.set_zlim([-150,0])
    ax.set_xticks([-200,0,200])
    ax.set_yticks([-150,0,150])
    ax.set_zticks([-150,-75,0])
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_zlabel('$z$')
    ax.set_title('Exact Dynamics', fontsize = 10)

    ax = plt.subplot(gs0[:, 1:2], projection='3d')
    ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
    colorline3d(ax, learned_X_star[:,0], learned_X_star[:,1], learned_X_star[:,2], cmap = plt.cm.jet)
    ax.grid(False)
    ax.set_xlim([-200,200])
    ax.set_ylim([-200,200])
    ax.set_zlim([-150,0])
    ax.set_xticks([-200,0,200])
    ax.set_yticks([-150,0,150])
    ax.set_zticks([-150,-75,0])
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_zlabel('$z$')
    ax.set_title('Learned Dynamics', fontsize = 10)

    # savefig('./figures/Cylinder', crop = False)

    ####### Plotting Vorticity ##################

    data = scipy.io.loadmat('./cylinder_vorticity.mat')
    XX = data['XX']
    YY = data['YY']
    WW = data['WW']

    fig, ax = newfig(1.0, 0.65)
    ax.axis('off')

    gs0 = gridspec.GridSpec(1, 1)
    gs0.update(top=0.85, bottom=0.2, left=0.25, right=0.8, wspace=0.15)

    ax = plt.subplot(gs0[0:1, 0:1])
    h = ax.pcolormesh(XX, YY, WW, cmap='seismic',shading='gouraud', vmin=-5, vmax=5)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_title('Vorticity', fontsize = 10)
    fig.colorbar(h)

    # savefig('./figures/Cylinder_vorticity', crop = False)

max_stars_repo_name: SeoulTechPSE/MultistepNNs

max_issues_repo_name: SeoulTechPSE/MultistepNNs

max_forks_repo_name: SeoulTechPSE/MultistepNNs

reason: import numpy,import scipy,from scipy


text:
# Copyright (C) IBM Corporation 2018
#
# 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.

# This is the main python file. Start the experiment by running this.

import os
import argparse
import random
import string
from multiprocessing.pool import ThreadPool
from time import sleep
import yaml

import numpy as np

import ga
from cloud_wrapper import Cloud

os.environ['REDIS_PW'] = "".join([random.choice(string.ascii_letters + string.digits) for _ in range(20)])
print(os.environ['REDIS_PW'])

import rq_worker  # It will read the password files

# Load settings yaml
with open('settings.yaml', 'rt') as f:
    settings = yaml.load(f)

def main():
    parser = argparse.ArgumentParser()
    required_args = parser.add_argument_group('required arguments')
    required_args.add_argument('-n', '--hostname-prefix', dest='hostname', type=str, required=True,
                               help="Choose a hostname prefix for the machines. "
                                    "If those already exist, they will be used without requesting new ones")
    parser.add_argument('--keep-vms', dest='keep_vms', action='store_true',
                        help="Do not delete the instantiated VMs at the end of the experiment.")
    parser.add_argument('--datacenter', dest='datacenter', type=str, required=True,
                        help="Datacenter to choose.")
    args, _ = parser.parse_known_args()

    # Ask for hostname prefix for the VMs
    hostname = args.hostname
    for i in range(settings['experiment']['count']):
        os.makedirs("out_"+str(i), exist_ok=False)

    cloud_mgr = Cloud(cpus=settings['vm']['cpus'],
                      mem=settings['vm']['mem'],
                      hostname=hostname,
                      count=settings['vm']['count'],
                      datacenter=args.datacenter)

    with ThreadPool(processes=settings['experiment']['count']) as pool:
        pool.map(run_experiment, range(settings['experiment']['count']))

    # Delete all the VMs
    if not args.keep_vms:
        cloud_mgr.cancel_all()

def run_experiment(id_: int):
    # Starting genetic algorithm
    gene_pool = ga.GA(popu_size=settings['experiment']['population'],
                      n_survivors=settings['experiment']['survivors'],
                      n_champions=settings['experiment']['champions'])

    for i in range(settings['experiment']['generations']):
        worker_async_results = [rq_worker.run.delay(seeds) for seeds in gene_pool.genomes]

        scores = wait_for_results(worker_async_results)

        survivors_scores, survivors_genomes = gene_pool.step(scores)

        np.savez("out_" + str(id_) + "/" + str(i) + ".npz", scores=survivors_scores, genomes=survivors_genomes)

        print("=== ID", id_, "Generation", i, "===\n",
              "Best scores:", list(reversed(survivors_scores)), "\n")

def wait_for_results(async_results: [], polling_sleep=2):
    # Get results (non-blocking)
    # If a result is not yet available, we get `None`
    results = [r.result for r in async_results]

    # While `results` contains `None` elements
    while any(e is None for e in results):
        sleep(polling_sleep)
        # Try to get the result for each `None` element
        for j, (r, ar) in enumerate(zip(results, async_results)):
            if r is None:
                results[j] = ar.result

    return results

if __name__ == '__main__':
    main()

max_stars_repo_name: aasseman/distributed-evolutionary-ml

max_issues_repo_name: aasseman/distributed-evolutionary-ml

max_forks_repo_name: aasseman/distributed-evolutionary-ml


text:
# -*- coding: utf-8 -*-
"""
Run the posterior sampling RL algorithm (PSRL) in the RiverSwim environment.

The PSRL algorithm is described in "(More) Efficient Reinforcement Learning via
Posterior Sampling," by I. Osband, B. Van Roy, and D. Russo (2013). It learns
from numerical rewards, rather than preferences.
"""

import scipy.io as io
import os

from Envs.RiverSwim import RiverSwimEnv
from Learning_algorithms.PSRL_numerical_rewards import PSRL

# Define constants:
time_horizon = 50
num_iter = 400     # Number of iterations of the learning algorithm. This is
                   # twice the number of iterations of the preference-based
                   # algorithms, since PSRL rolls out one trajectory/episode
                   # per learning iteration, while the preference-based
                   # algorithms roll out two; thus, the number of episodes is
                   # kept consistent.

run_nums = 100    # Number of times to run the algorithm

# Folder for saving results:
output_folder = 'PSRL/'

if not os.path.isdir(output_folder):
    os.mkdir(output_folder)

# Hyperparameters to use for the PSRL algorithm:
NG_params = [1, 1, 1, 1]      # Prior parameters for normal-gamma reward model

# Instantiate the environment:
env = RiverSwimEnv()

# Run PSRL algorithm:
for run_num in range(run_nums):

    # String to use in status updates:
    run_str = '%s, run %i' % (NG_params, run_num)

    # Run algorithm:
    rewards = PSRL(time_horizon, NG_params, env, num_iter, run_str = run_str)

    # Save results from this algorithm run:
    output_filename = output_folder + 'Iter_' + str(num_iter) + '_params_' + \
                        str(NG_params[0]) + '_' + str(NG_params[1]) \
                         + '_' + str(NG_params[2]) + '_' + str(NG_params[3]) + \
                        '_run_' + str(run_num) + '.mat'

    io.savemat(output_filename, {'rewards': rewards, 'num_iter': num_iter,
              'NG_params': NG_params})

max_stars_repo_name: ernovoseller/DuelingPosteriorSampling

max_issues_repo_name: ernovoseller/DuelingPosteriorSampling

max_forks_repo_name: ernovoseller/DuelingPosteriorSampling

reason: import scipy


text:
"""
Continuous to discrete.
"""

import gym
import numpy as np

class C2DEnv(object):
    """
    Wrapper environment for converting continuous action space to multi discrete action space.

    Parameters
    ----------
    env : gym.Env
    n_bins : int
        Number of bins for converting continuous to discrete.
        e.g. continuous action space is 0 ~ 1 and n_bins=5,
        action space is converted to [0, 0.25, 0.5, 0.75, 1]
    """

    def __init__(self, env, n_bins=30):
        assert isinstance(env.action_space, gym.spaces.Box)
        assert len(env.action_space.shape) == 1
        self.env = env
        self.n_bins = n_bins
        self.action_space = gym.spaces.MultiDiscrete(
            env.action_space.shape[0] * [n_bins])
        self.observation_space = self.env.observation_space
        if hasattr(env, 'original_env'):
            self.original_env = env.original_env
        else:
            self.original_env = env

    @property
    def horizon(self):
        if hasattr(self.env, 'horizon'):
            return self.env._horizon

    def reset(self):
        return self.env.reset()

    def step(self, action):
        continuous_action = []
        for a, low, high in zip(action, self.env.action_space.low, self.env.action_space.high):
            continuous_action.append(np.linspace(low, high, self.n_bins)[a])
        action = np.array(continuous_action)
        next_obs, reward, done, info = self.env.step(action)
        return next_obs, reward, done, info

    def render(self):
        self.env.render()

    def terminate(self):
        self.env.terminate()

max_stars_repo_name: AswinRetnakumar/Machina

max_issues_repo_name: AswinRetnakumar/Machina

max_forks_repo_name: AswinRetnakumar/Machina


text:
from itertools import product
from collections import defaultdict

import pandas
import numpy as np
from libsvmdata import fetch_libsvm
from joblib import Parallel, delayed, parallel_backend

from andersoncd.data.real import load_openml
from andersoncd.lasso import solver_enet, apcg_enet

# to generate the exact fig of the paper:
# dataset_names = [
#     "leukemia", "gina_agnostic", "hiva_agnostic", 'rcv1.binary']
# div_alphas = [10, 100, 1000, 5000]

# if you want to run the file quickly choose instead:
# dataset_names = ["leukemia"]
dataset_names = ["rcv1.binary"]
div_alphas = [1000]

algos = [
    ['apcg', False, 5],
    ['pgd', False, 5],
    ['pgd', True, 5],
    ['fista', False, 5],
    ['cd', False, 5],
    ['rcd', False, 5],
    ['cd', True, 5]
]

dict_maxiter = defaultdict(lambda: 1_000_000, key=None)

dict_maxiter = {}
dict_maxiter["leukemia", 10] = 550
dict_maxiter["mushroom", 10] = 10_000
dict_maxiter["gina_agnostic", 10] = 1_000
dict_maxiter["hiva_agnostic", 10] = 10_000
dict_maxiter["upselling", 10] = 10_000
dict_maxiter["rcv1.binary", 10] = 10_000
dict_maxiter["news20.binary", 10] = 100_000
dict_maxiter["kdda_train", 10] = 10_000
dict_maxiter["finance", 10] = 5_000

dict_maxiter["leukemia", 100] = 5_000
dict_maxiter["mushroom", 100] = 10_000
dict_maxiter["gina_agnostic", 100] = 1_000
dict_maxiter["hiva_agnostic", 100] = 5000
dict_maxiter["upselling", 100] = 10_000
dict_maxiter["rcv1.binary", 100] = 10_000
dict_maxiter["news20.binary", 100] = 10_000
dict_maxiter["kdda_train", 100] = 1_000
dict_maxiter["finance", 100] = 50_000

dict_maxiter["leukemia", 1000] = 100_000
dict_maxiter["mushroom", 1000] = 100_000
dict_maxiter["gina_agnostic", 1000] = 100_000
dict_maxiter["hiva_agnostic", 1000] = 100_000
dict_maxiter["upselling", 1000] = 100_000
dict_maxiter["rcv1.binary", 1000] = 100_000
dict_maxiter["news20.binary", 1000] = 1_000_000
dict_maxiter["kdda_train", 1000] = 1_000
dict_maxiter["finance", 1000] = 50_000

dict_maxiter["leukemia", 5_000] = 300_000
dict_maxiter["mushroom", 5_000] = 300_000
dict_maxiter["gina_agnostic", 5_000] = 300_000
dict_maxiter["hiva_agnostic", 5_000] = 300_000
dict_maxiter["upselling", 5000] = 100_000
dict_maxiter["rcv1.binary", 5000] = 500_000
dict_maxiter["news20.binary", 5000] = 1_000_000
dict_maxiter["kdda_train", 5000] = 1_000

dict_f_gap = {}
dict_f_gap["leukemia"] = 10
dict_f_gap["mushroom"] = 10
dict_f_gap["gina_agnostic"] = 10
dict_f_gap["hiva_agnostic"] = 10
dict_f_gap["upselling"] = 10
dict_f_gap["rcv1.binary"] = 10
dict_f_gap["news20.binary"] = 10
dict_f_gap["kdda_train"] = 10
dict_f_gap["finance"] = 50

dict_tmax = {}
dict_tmax["leukemia", 10] = 10
dict_tmax["leukemia", 100] = 100
dict_tmax["leukemia", 1000] = 1_000
dict_tmax["leukemia", 5000] = 3_600

dict_tmax["hiva_agnostic", 10] = 5
dict_tmax["hiva_agnostic", 100] = 100
dict_tmax["hiva_agnostic", 1000] = 3_600
dict_tmax["hiva_agnostic", 5000] = 7_200

dict_tmax["gina_agnostic", 10] = 5
dict_tmax["gina_agnostic", 100] = 100
dict_tmax["gina_agnostic", 1000] = 100
dict_tmax["gina_agnostic", 5000] = 100

dict_tmax["rcv1.binary", 10] = 5
dict_tmax["rcv1.binary", 100] = 10
dict_tmax["rcv1.binary", 1000] = 300
dict_tmax["rcv1.binary", 5000] = 1500

dict_tmax["news20.binary", 10] = 60
dict_tmax["news20.binary", 100] = 10
dict_tmax["news20.binary", 1000] = 120
dict_tmax["news20.binary", 5000] = 600

def parallel_function(dataset_name, algo, div_alpha):
    algo_name, use_acc, K = algo
    if dataset_name.startswith((
            'rcv1.binary', 'news20.binary', 'kdda_train', 'finance')):
        X, y = fetch_libsvm(dataset_name, normalize=True)
        y -= y.mean()
        y /= np.linalg.norm(y)
    else:
        X, y = load_openml(dataset_name)

    alpha_max = np.max(np.abs(X.T @ y))
    alpha = alpha_max / div_alpha
    tol = 1e-14
    f_gap = dict_f_gap[dataset_name]

    max_iter = dict_maxiter[dataset_name, div_alpha]

    tmax = dict_tmax[dataset_name, div_alpha]

    if algo_name == 'apcg':
        w, E, gaps, times = apcg_enet(
            X, y, alpha, max_iter=max_iter, tol=tol, f_gap=f_gap,
            compute_time=True, tmax=5, verbose=True)
    else:
        w, E, gaps, times = solver_enet(
            X, y, alpha, f_gap=f_gap, max_iter=max_iter, tol=tol,
            use_acc=use_acc, K=K, algo=algo_name, compute_time=True,
            tmax=5, verbose=True, seed=42)

    if algo_name == 'apcg':
        w, E, gaps, times = apcg_enet(
            X, y, alpha, max_iter=max_iter, tol=tol, f_gap=f_gap,
            compute_time=True, tmax=tmax, verbose=True)
    else:
        w, E, gaps, times = solver_enet(
            X, y, alpha, f_gap=f_gap, max_iter=max_iter, tol=tol,
            use_acc=use_acc, K=K, algo=algo_name, compute_time=True,
            tmax=tmax, verbose=True, seed=42)

    my_results = (
        dataset_name, algo_name, use_acc, K, div_alpha, w,
        E, gaps, f_gap, times)
    df = pandas.DataFrame(my_results).transpose()
    df.columns = [
        'dataset', 'algo_name', 'use_acc', 'K', 'div_alpha', "optimum", "E",
        "gaps", "f_gaps", "times"]
    str_results = "results/%s_%s_%s%i.pkl" % (
        dataset_name, algo_name, str(use_acc), div_alpha)
    df.to_pickle(str_results)

print("enter parallel")
n_jobs = len(dataset_names) * len(div_alphas) * len(algos)

# n_jobs = 1
n_jobs = min(n_jobs, 15)

with parallel_backend("loky", inner_max_num_threads=1):
    Parallel(
        n_jobs=n_jobs, verbose=100)(
        delayed(parallel_function)(
            dataset_name, algo, div_alpha)
        for dataset_name, algo, div_alpha in product(
            dataset_names, algos, div_alphas))

print('OK finished parallel')


text:
"""
MIT License

Copyright (c) 2018 Simon Olofsson

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

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

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

import numpy as np

"""
Model super class
"""
class VTHR2014Model:
	"""
	Case study from
	J. Vanlier, C. A. Tiemann, P. A. J. Hilbers, and N. A. W. van Riel
	"Optimal experiment design for model selection in biochemical networks"
	BMC Systems Biology, 8(20), 2014.
	"""
	@property
	def name (self):
		return 'M' + str(self.num_p)
	@property
	def n_outputs (self):
		return 1
	@property
	def x_bounds (self):
		return np.array([[-2.5, 2.5]])
	@property
	def p_bounds (self):
		return np.array([[-10,10]] * self.num_p)

	def __call__ (self, x, P, grad=False):
		F = [
			lambda t: t,
			lambda t: t**2,
			lambda t: np.sin(0.2 * t**3),
		 	lambda t: np.sin(2 * t) * t
		 	]
		F  = F[:self.num_p]

		f  = np.sum([ p * f(x) for p,f in zip(P,F)], axis=0)
		df = np.array([ f(x) for f in F ]).T
		return f if not grad else [f, df]

"""
Models
"""
class M1 (VTHR2014Model):
	@property
	def num_p (self):
		return 1

class M2 (VTHR2014Model):
	@property
	def num_p (self):
		return 2

class M3 (VTHR2014Model):
	@property
	def num_p (self):
		return 3

class M4 (VTHR2014Model):
	@property
	def num_p (self):
		return 4

"""
Data generator
"""
class DataGen (M3):
	@property
	def truemodel (self):
		return 2
	@property
	def measvar (self):
		return np.array([0.5])
	@property
	def p (self):
		return [-1, 0.5, -6]

	def __call__ (self,x):
		state = super(DataGen, self).__call__(x,self.p)
		noise = np.sqrt(self.measvar) * np.random.randn()
		return state + noise

"""
Get model functions
"""
def get ():
	return DataGen(), [M1(),M2(),M3(),M4()]

max_stars_repo_name: scwolof/GPdoemd

max_issues_repo_name: scwolof/GPdoemd

max_forks_repo_name: scwolof/GPdoemd


text:
#-*- coding: utf8

import _trend

from sklearn.base import BaseEstimator
from sklearn.base import ClassifierMixin

import numpy as np

class TrendLearner(BaseEstimator, ClassifierMixin):

    def __init__(self, num_steps, gamma=1):
        self.num_steps = num_steps
        self.gamma = gamma
        self.num_labels = 0
        self.R = None
        self.labels = None

    def fit(self, X, y):

        self.R = np.asanyarray(X, dtype=np.float64, order='C')

        y = np.asanyarray(y)
        unique, labels_flat = np.unique(y, return_inverse=True)
        self.labels = labels_flat.reshape(y.shape)
        self.num_labels = unique.shape[0]

    def predict(self, X):

        X = np.asanyarray(X)[:, :self.num_steps]
        X = np.asanyarray(X, dtype=np.float64, order='C')

        P = _trend.predict(X, self.R, self.labels, self.num_labels, self.gamma)

        return P.argmax(axis=1)

    def predict_proba(self, X):

        X = np.asanyarray(X)[:, :self.num_steps]
        X = np.asanyarray(X, dtype=np.float64, order='C')

        P = _trend.predict(X, self.R, self.labels, self.num_labels, self.gamma)
        P = ((P.T / P.sum(axis=1)).T)

        return P

max_stars_repo_name: flaviovdf/pyksc

max_issues_repo_name: flaviovdf/pyksc

max_forks_repo_name: flaviovdf/pyksc


text:
import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import MinMaxScaler
#X=2次元train_data,出力時に整形して1次元にする
class KMeansFeaturizer:
    """
    数値データをk-meansのクラスタIDに変換します。

    この変換器は入力データに対してk-meansを実行し、各データポイントを最も近い
    クラスタのIDに変換します。ターゲット変数yが存在する場合、クラス分類の境界を
    より重視したクラスタリングの結果を得るために、ターゲット変数をスケーリングして
    入力データに含めてk-meansに渡します。
    """
    #k=出力するクラスタの数

    def __init__(self, k=100, target_scale=5.0, random_state=None):
        self.k = k
        self.target_scale = target_scale
        self.random_state = random_state
        # np.array(range(k)).reshape(-1,1) ->2次元配列にそれぞれ一個ずつ要素が入っている。([[0],[1]...])
        #要素が入ってきたとき、one-hotvectorに変換する?([0]->[1,0,0...])
        self.cluster_encoder = OneHotEncoder(categories='auto').fit(
            np.array(range(k)).reshape(-1, 1))

    def fit(self, X, y=None):
        """
        入力データに対しk-meansを実行し、各クラスタの中心を見つけます。
        """
        if y is None:
            # ターゲット変数がない場合、ふつうのk-meansを実行します。
            km_model = KMeans(n_clusters=self.k, n_init=20,
                              random_state=self.random_state)
            km_model.fit(X)
            self.km_model = km_model
            self.cluster_centers_ = km_model.cluster_centers_
            return self

        # ターゲット変数がある場合。スケーリングして入力データに含めます。
        data_with_target = np.hstack((X, y[:, np.newaxis] * self.target_scale))

        # ターゲットを組み入れたデータで事前学習するためのk-meansモデルを構築します。
        km_model_pretrain = KMeans(
            n_clusters=self.k, n_init=20, random_state=self.random_state)
        km_model_pretrain.fit(data_with_target)

        # ターゲット変数の情報を除いて元の空間におけるクラスタを得るために
        # k-meansを再度実行します。事前学習で見つけたクラスタの中心を
        # 使って初期化し、クラスタの割り当てと中心の再計算を1回にします。
        km_model = KMeans(n_clusters=self.k, init=km_model_pretrain.cluster_centers_[
                          :, :2], n_init=1, max_iter=1)
        km_model.fit(X)

        self.km_model = km_model
        self.cluster_centers_ = km_model.cluster_centers_
        return self

    def transform(self, X, y=None):
        """
        入力データポイントに最も近いクラスタのIDを返します。
        """
        clusters = self.km_model.predict(X)
        #clustersが[[0],[8]...]みたいに代入されてそれがワンホットベクトルとして変換されて出力される?
        return self.cluster_encoder.transform(clusters.reshape(-1, 1))

    def fit_transform(self, X, y=None):
        self.fit(X, y)
        return self.transform(X, y)

    if __name__ == "__main__":
        import feather as ftr
        train = ftr.read_dataframe('../data/input/train.feather')

        train.drop(columns=['target', 'ID_code'], inplace=True)
        train = train.iloc[:10000, :]
        kmf = KMeansFeaturizer()
        k_train = kmf.fit_transform(train)

max_stars_repo_name: zerebom/santander

max_issues_repo_name: zerebom/santander

max_forks_repo_name: zerebom/santander


text:
# This file purify the clip_raw

import pandas as pd
import numpy as np
from math import isnan

def clip_purifier(srcPath, dstPath):
    clipRawDf = pd.read_csv(srcPath)
    clipRawArr = np.array(clipRawDf)
    clipRawList = clipRawArr.tolist()
    rawLen = len(clipRawList)

    clipPureList = []
    for i in range(int(rawLen/2)):
        oddLine = clipRawList[i*2]
        evenLine = clipRawList[i*2 + 1]

        if isnan(oddLine[3]) or (isnan(oddLine[5]) and isnan(oddLine[11])) or (isnan(oddLine[17]) and isnan(oddLine[23])) or \
                (isnan(oddLine[19]) and isnan(oddLine[25])):
            continue
        if isnan(evenLine[3]) or (isnan(evenLine[5]) and isnan(evenLine[11])) or (isnan(evenLine[17]) and isnan(evenLine[23])) or \
                (isnan(evenLine[19]) and isnan(evenLine[25])):
            continue
        if isnan(oddLine[37]) or isnan(evenLine[37]):
            continue

        if isnan(oddLine[5]):
            oddLine[5] = oddLine[11]
            oddLine[6] = oddLine[12]
        if isnan(oddLine[11]):
            oddLine[11] = oddLine[5]
            oddLine[12] = oddLine[6]
        if isnan(oddLine[17]):
            oddLine[17] = oddLine[23]
            oddLine[18] = oddLine[24]
        if isnan(oddLine[23]):
            oddLine[23] = oddLine[17]
            oddLine[24] = oddLine[18]
        if isnan(oddLine[19]):
            oddLine[19] = oddLine[25]
            oddLine[20] = oddLine[26]
        if isnan(oddLine[25]):
            oddLine[25] = oddLine[19]
            oddLine[26] = oddLine[20]

        if isnan(evenLine[5]):
            evenLine[5] = evenLine[11]
            evenLine[6] = evenLine[12]
        if isnan(evenLine[11]):
            evenLine[11] = evenLine[5]
            evenLine[12] = evenLine[6]
        if isnan(evenLine[17]):
            evenLine[17] = evenLine[23]
            evenLine[18] = evenLine[24]
        if isnan(evenLine[23]):
            evenLine[23] = evenLine[17]
            evenLine[24] = evenLine[18]
        if isnan(evenLine[19]):
            evenLine[19] = evenLine[25]
            evenLine[20] = evenLine[26]
        if isnan(evenLine[25]):
            evenLine[25] = evenLine[19]
            evenLine[26] = evenLine[20]

        pureLine = [evenLine[3] - oddLine[3], evenLine[4] - oddLine[4]]
        pureLine.append(evenLine[5] - oddLine[5])
        pureLine.append(evenLine[6] - oddLine[6])
        pureLine.append(evenLine[11] - evenLine[11])
        pureLine.append(evenLine[12] - oddLine[12])
        pureLine.append(evenLine[17] - oddLine[17])
        pureLine.append(evenLine[18] - oddLine[18])
        pureLine.append(evenLine[23] - oddLine[23])
        pureLine.append(evenLine[24] - oddLine[24])
        pureLine.append(evenLine[19] - oddLine[19])
        pureLine.append(evenLine[20] - oddLine[20])
        pureLine.append(evenLine[25] - oddLine[25])
        pureLine.append(evenLine[26] - oddLine[26])

        pureLine.append(oddLine[37])
        pureLine.append(oddLine[0])
        clipPureList.append(pureLine)

    clipPureDf = pd.DataFrame(clipPureList, columns=['nx', 'ny', 'rsx', 'rsy', 'lsx', 'lsy', 'rhx', 'rhy', 'lhx', 'lhy',
                                                     'rkx', 'rky', 'lkx', 'lky', 'label', 'original_index'])
    clipPureDf.dropna(inplace=True)
    clipPureDf.to_csv(dstPath)

clip_purifier('data/clipTrainRaw.csv', 'data/clipTrainPure.csv')
clip_purifier('data/clipTestRaw.csv', 'data/clipTestPure.csv')


text:
import numpy as np
import tensorflow as tf
import math
import os

tf.enable_eager_execution()

import waymo_open_dataset
from waymo_open_dataset.utils import range_image_utils
from waymo_open_dataset.utils import transform_utils
from waymo_open_dataset.utils import  frame_utils
from waymo_open_dataset import dataset_pb2 as open_dataset
from Adapter_utilities import *

def label_ext_with_occlusion(i, filename, Label_all, Label):

    FileName = filename
    dataset = tf.data.TFRecordDataset(FileName, compression_type='')
    # __lidar_list = ['_FRONT', '_FRONT_LEFT', '_FRONT_RIGHT', '_SIDE_LEFT', '_SIDE_RIGHT']
    __lidar_list = ['_FRONT', '_FRONT_RIGHT', '_FRONT_LEFT', '_SIDE_RIGHT', '_SIDE_LEFT']
    __type_list = ['unknown', 'Vehicle', 'Pedestrian', 'Sign', 'Cyclist']

    for data in dataset:

        i_str = '{0:06}'.format(i)

        frame = open_dataset.Frame()
        frame.ParseFromString(bytearray(data.numpy()))

        file_path = os.path.join(Label_all, "%s.txt"%(i_str))

        f_open = open(file_path, 'w+')

        id_to_bbox = dict()
        id_to_name = dict()
        for labels in frame.projected_lidar_labels:
            name = labels.name
            for label in labels.labels:
                bbox = [label.box.center_x - label.box.length / 2, label.box.center_y - label.box.width / 2,
                        label.box.center_x + label.box.length / 2, label.box.center_y + label.box.width / 2]
                id_to_bbox[label.id] = bbox
                id_to_name[label.id] = name - 1

        Tr_velo_to_cam = []
        waymo_cam_RT = np.array([0,-1,0,0,  0,0,-1,0,   1,0,0,0,    0 ,0 ,0 ,1]).reshape(4,4)

        for camera in frame.context.camera_calibrations:
            tmp=np.array(camera.extrinsic.transform).reshape(4,4)
            tmp=np.linalg.inv(tmp)
            tmp = np.matmul(waymo_cam_RT, tmp)
            Tr_velo_to_cam.append(tmp)

        laser = open_dataset.LaserName.TOP
        laser_calib = [obj for obj in frame.context.laser_calibrations if obj.name == laser]
        laser_calib = laser_calib[0]

        RI, CP, RITP = frame_utils.parse_range_image_and_camera_projection(frame)

        pcl, pcl_attr = frame_utils.convert_range_image_to_point_cloud(frame, RI, CP, RITP)
        pcl = np.concatenate(pcl, axis=0)

        vehicle_to_labels = [np.linalg.inv(get_box_transformation_matrix(label.box)) for label in frame.laser_labels]
        vehicle_to_labels = np.stack(vehicle_to_labels)

        pcl1 = np.concatenate((pcl,np.ones_like(pcl[:,0:1])),axis=1)
        proj_pcl = np.einsum('lij,bj->lbi', vehicle_to_labels, pcl1)
        mask = np.logical_and.reduce(np.logical_and(proj_pcl >= -1, proj_pcl <= 1),axis=2)

        counts = mask.sum(1)
        visibility = counts > 10

        for obj, visib in zip(frame.laser_labels, visibility):

            # caculate bounding box
            bounding_box = None
            name = None
            id = obj.id
            for lidar in __lidar_list:
                if id + lidar in id_to_bbox:
                    bounding_box = id_to_bbox.get(id + lidar)
                    name = str(id_to_name.get(id + lidar))
                    break
            if bounding_box == None or name == None:
                continue

            my_type = __type_list[obj.type]
            if visib:
                occluded = 0
            else:
                occluded = 1
            truncated = 0.00
            height = obj.box.height
            width = obj.box.width
            length = obj.box.length
            x = obj.box.center_x
            y = obj.box.center_y
            z = obj.box.center_z
            rotation_y = obj.box.heading

            z -= height/2

            transform_box_to_cam = Tr_velo_to_cam[int(name)] @ get_box_transformation_matrix(obj.box)
            pt1 = np.array([-0.5, 0.5, 0 , 1.])
            pt2 = np.array([0.5, 0.5, 0 , 1.])
            pt1 = np.matmul(transform_box_to_cam, pt1)
            pt2 = np.matmul(transform_box_to_cam, pt2)
            new_ry = math.atan2(pt2[2]-pt1[2],pt2[0]-pt1[0])
            rotation_y = -new_ry

            new_loc = np.matmul(Tr_velo_to_cam[int(name)], np.array([x,y,z,1]).T)
            x, y, z = new_loc[:3]

            beta = math.atan2(x, z)
            alpha = (rotation_y + beta - math.pi / 2) % (2 * math.pi)

            # save the labels
            line = my_type + ' {} {} {} {} {} {} {} {} {} {} {} {} {} {}\n'.format(round(truncated, 2),
                                                                                round(occluded, 2),
                                                                                round(alpha, 2),
                                                                                round(bounding_box[0], 2),
                                                                                round(bounding_box[1], 2),
                                                                                round(bounding_box[2], 2),
                                                                                round(bounding_box[3], 2),
                                                                                round(height, 2),
                                                                                round(width, 2),
                                                                                round(length, 2),
                                                                                round(x, 2),
                                                                                round(y, 2),
                                                                                round(z, 2),
                                                                                round(rotation_y, 2))
            line_all = line[:-1] + ' ' + '\n'
            # # store the label

            name_of_file = ("%s.txt"%(i_str))
            fp_label = open(Label + name + '/' + name_of_file, 'a')
            fp_label.write(line)
            fp_label.close()

            f_open.write(line_all)
        f_open.close()

        i = i + 1
    return i