text:
#!usr/bin/python
"""
Script to download Open Street Map data of land use. Converts to GeoJSON
Lukas Arnold
WB-DIME
Jan 31 2017
"""
from country_bounding_boxes import (
      country_subunits_containing_point,
      country_subunits_by_iso_code
    )
from utils.overpass_client import OverpassClient
from numpy import floor
from modules.getFeatures import find_between
import json
import urllib2

# Query string to OSM to provide landuse data, with the query key later
# placed in the middle
query_begin = '''\
way
  ['''
query_end=''']
  ({query_bb_s},{query_bb_w},{query_bb_n},{query_bb_e});
(._;>;);
out;
'''
# Output is OSM format

datatype="GeoJSON"

def script(countryISO='US',query='landuse',outputFolder='data/',partOfData=1,
	outputFile='OSMdata_'):
	"""
	Main function executed by top

	'countryISO': Country for which BBox data should be downloaded. Can also contain custom boundary box
	'query': Tag for OSM query to search for
	'partOfData': Part of total data of a country to be processed
	Returns list with [filename,datatype], where datatype is the
		GDAL_CODE
	"""
	partOfData=float(partOfData)
	subunits=[]
	countryISOlist=countryISO.split()
	if countrISOlist[1]: #if there is more than one entry
		bbox=[]
		bbox.append(float(contryISOlist[0]))
		bbox.append(float(contryISOlist[1]))
		bbox.append(float(contryISOlist[2]))
		bbox.append(float(contryISOlist[3]))

	else:
		#Load country data
		for c in country_subunits_by_iso_code(countryISO):
			subunits.append(c)
		#Chose subunits, if more than 1
		subunit=1
		if len(subunits)>1: #if there are subunits
			cnt = 1
			print 	"Subunits:"
			for c in subunits:
				print cnt,"- ",c.name
				cnt += 1
			subunit=input('Chose subunit: ')
		elif len(subunits)==0: #if nothing found
			print "Error: No country or entry with ISO code",countryISO
			exit()
		#Get BBox data for country
		print "Acquiring data for",subunits[subunit-1].name
		bbox = subunits[subunit-1].bbox #0-w, 1-s, 2-e, 3-n
	w = bbox[0]
	s = bbox[1]
	e = bbox[2]
	n = bbox[3]
	print "Coordinates:",w,s,e,n
	print "Key:",query
	# Country is split into 100 boxes, as (for the us) sample is too big
	# (timeout)
	# Number of Boxes = (samples-1)^2 boxes.
	#calculate number of boxes
	mindiff=min([abs(w-e),abs(n-s)])
	samples=int(floor((mindiff*8)+2))
	print "Number of queries:",(samples-1)**2
	#samples = 11  # 100 boxes
	fullquery = query_begin+query+query_end
	#Get Elements from OSM
	overpass_client = OverpassClient(endpoint='fr')
	d = overpass_client.get_bbox_elements(
	    ql_template=fullquery,
	    bb_s=s, bb_w=w, bb_n=n, bb_e=e,
	    samples=samples)
	lene=len(d)
	print 'Total elements found: %d' % lene
	boundery_index=int(floor(partOfData*lene))
	d=d[0:boundery_index]
	dr=list(reversed(d))
	lene=boundery_index

	fileName=outputFolder+'/'+outputFile+str(subunits[subunit-1].name).replace(" ","_")+".json"

	#Create GeoJSON string
	geojson=[]
	geojson.append('''
		{
	    "type": "FeatureCollection",
		"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::4326" } },
	    "features": [''')
	#create faster node searchs
	print "Create library for faster node searches"
	#get max id
	cnt = 0.
	ids={}
	for e in d:
		print "\t Library creation:",int(cnt*100./lene),"% done.\r",
		cnt+=1
		if e['type']=='node':
			ids[e['id']]=[e['lon'],e['lat']]

	#creade list of nodes with ids.
	#coordlist=[]
	#cnt = 0.
	#for i in range(0,maxid+1):
	#	print "\t Step 2/3:",cnt*100./maxid+1,"% done.\r",
	#	cnt+=1
	#	coordlist.append([0,0])
	#cnt = 0.
	#rint ""
	#for e in d:
	#	print "\t Step 3/3:",cnt*100./lene,"% done.\r",
	#	cnt+=1
	#	if e['type']=='node':
	#		coordlist(e['id'])[0]=e['lon']
	#		coordlist(e['id'])[1]=e['lat']
	#loop through elements and append to GeoJSON string
	cnt = 0.
	cnte = 0
	print ""
	print "Convert to GeoJSON file",fileName
	writetag=[]
	writecoord=[]
	for e in d :
		cnt+=1
		print "\tGet elemenents:",int(cnt*100./lene),"% done.\r",
		if e['type']=='way':
	#		if  e['area']=='yes':
				writetag.append(e['tags'][query])
				writecoord.append([])
				for node in e['nodes']:
					try:
						lon=str(ids[node][0])
						lat=str(ids[node][1])
					except KeyError:
						print ''
						print '\tNode',node,'not found in library. Download informations from openstreetmap.org ...'
						response=urllib2.urlopen('http://api.openstreetmap.org/api/0.6/node/'+str(node))
						fullxml = str(response.read())
						lon=find_between(fullxml,"lon=\"","\"",lastfirst=True)
						lat=find_between(fullxml,"lat=\"","\"",lastfirst=True)
					writecoord[cnte].append([lon,lat])
				cnte+=1
	cnte2=0
	print ''
	for tag in writetag:
		print "\tCreate GeoJSON:",int(cnte2*100./cnte),"% done.\r",
		geojson.append('''
		    {
			"type": "Feature",'''+\
				"\n\t\t\t\"properties\":{\"Descriptio\":\""+\
				tag+"\"},")
		geojson.append('''
			"geometry" : {
			    "type": "MultiPolygon",
			    "coordinates":[[[''')
		for coord in writecoord[cnte2]:
			geojson.append("["+coord[0]+","+coord[1]+"],")
		cnte2+=1
		geojson[-1]=geojson[-1][0:-1]+"]]]}},"
		#geojson=geojson[0:-1]+"]]]}},"
	geojson=''.join(geojson)
	geojson=geojson[0:-1]+"\n]\n}"
	print " "
	with open(fileName, 'w+') as f:
			json.dump(geojson,f) #save as geojson file
	#replace escape characters
	with open(fileName, 'r') as file :
	  filedata = file.read()[1:-1]
	# Replace the target string
	filedata = filedata.replace('\\n', '')
	filedata = filedata.replace('\\t', '')
	filedata = filedata.replace('\\"', '"')
	# Save the result
	with open(fileName, 'w') as file:
	  file.write(filedata)

	print "Written to",fileName
	return [fileName,datatype]

reason: from numpy


text:
"""Sky brightnes approzimation using Zernike polynomials

The form and notation used here follow:

Thibos, L. N., Applegate, R. A., Schwiegerling, J. T., Webb, R. & VSIA
Standards Taskforce Members. Vision science and its
applications. Standards for reporting the optical aberrations of
eyes. J Refract Surg 18, S652-660 (2002).
"""

# imports
from math import factorial
import logging
import os
import warnings
from glob import glob
from functools import lru_cache
import numpy as np
import pandas as pd
from numexpr import NumExpr
from sklearn.linear_model import LinearRegression
import scipy.optimize
from scipy.interpolate import interp1d
import palpy
import healpy

import rubin_sim.utils as utils
from rubin_sim.data import get_data_dir

# constants

logging.basicConfig(format="%(asctime)s %(message)s")
LOGGER = logging.getLogger(__name__)

TELESCOPE = utils.Site("LSST")
SIDEREAL_TIME_SAMPLES_RAD = np.radians(np.arange(361, dtype=float))
BANDS = ("u", "g", "r", "i", "z", "y")

# exception classes

# interface functions

def fit_pre(npy_fname, npz_fname, *args, **kwargs):
    """Fit Zernike coefficients to a pre-computed data set

    Parameters
    ----------
    npy_fname : `str`
        File name of the SkyBrightessPre <MJD>_<MDJ>.npy file
    npz_fname : `str`
        File name of the SkyBrightessPre <MJD>_<MDJ>.npz file

    other arguments are passed to the ZernikeSky constructor.

    Returns
    -------
    zernike_coeffs : `pd.DataFrame`
        A DataFrame with the coefficients, indexed by band and mjd.
    """
    # Load the pre-computed data
    npz = np.load(npz_fname, allow_pickle=True)
    npz_hdr = npz["header"][()]
    npz_data = npz["dict_of_lists"][()]
    pre_sky = np.load(npy_fname, allow_pickle=True)

    mjds = npz_data["mjds"]
    alt = npz_hdr["alt"]
    az = npz_hdr["az"]
    zernike_coeffs_by_band = []
    zernike_sky = ZernikeSky(*args, **kwargs)
    for band in pre_sky.dtype.fields.keys():
        LOGGER.info("Starting %s band", band)
        zernike_coeff_arrays = []
        for mjd_idx, mjd in enumerate(mjds):
            zernike_coeff_arrays.append(
                zernike_sky.fit_coeffs(alt, az, pre_sky[band][mjd_idx], mjd)
            )
            if mjd_idx % 1000 == 0:
                msg = f"Finished {mjd_idx*100.0/float(len(mjds)):.2f}%"
                LOGGER.debug(msg)
        zernike_coeffs_by_band.append(
            pd.DataFrame(
                zernike_coeff_arrays,
                columns=np.arange(len(zernike_coeff_arrays[0])),
                index=pd.MultiIndex.from_arrays(
                    [np.full_like(mjds, band, dtype=type(band)), mjds],
                    names=["band", "mjd"],
                ),
            )
        )

    zernike_coeffs = pd.concat(zernike_coeffs_by_band)
    return zernike_coeffs

def bulk_zernike_fit(data_dir, out_fname, *args, **kwargs):
    """Fit Zernike coeffs to all SkyBrightnessPre files in a directory.

    Parameters
    ----------
    data_dir : `str`
        Name of the directory in which to look for SkyBrightnessPre
        data files.
    out_fname: `str`
        Name of the file in which to save fit coefficients.

    other arguments are passed to the ZernikeSky constructor.

    Returns
    -------
    zernike_coeffs : `pd.DataFrame`
        A DataFrame with the coefficients, indexed by band and mjd.
    """
    zernike_coeff_batches = []
    for npz_fname in glob(os.path.join(data_dir, "?????_?????.npz")):
        LOGGER.info("Processing %s", npz_fname)
        npy_fname = os.path.splitext(npz_fname)[0] + ".npy"
        zernike_coeff_batch = fit_pre(npy_fname, npz_fname, *args, **kwargs)
        zernike_coeff_batches.append(zernike_coeff_batch)

    zernike_coeffs = pd.concat(zernike_coeff_batches)
    zernike_coeffs.sort_index(level="mjd", inplace=True)

    if out_fname is not None:
        zernike_coeffs.to_hdf(out_fname, "zernike_coeffs", complevel=6)

        zernike_sky = ZernikeSky(*args, **kwargs)
        zernike_metadata = pd.Series(
            {"order": zernike_sky.order, "max_zd": zernike_sky.max_zd}
        )
        zernike_metadata.to_hdf(out_fname, "zernike_metadata")

    return zernike_coeffs

    # classes

class ZernikeSky:
    """Zernike sky approximator.

    Parameters
    ----------
    order : `int`, optional
        The order of the Zernike polynomial to use. Default is 6.
    nside : `int`, optional
        The nside of the healpix array to pre-compute Zernike Z terms for.
        Default is 32.
    max_zd : `float`, optional
        The maximum zenith distance, in degrees. This value will correspond
        to rho=1 in the Thibos et al. (2002) notation.
        Default is 67.
    dtype : `type`: optional
        The numpy type to use for all calculations. Default is `np.float64`.
    """

    def __init__(self, order=6, nside=32, max_zd=67, dtype=np.float64):
        self.order = order
        self.dtype = dtype
        self.nside = nside

        # Sets the value of zd where rho (radial coordinate of the
        # unit disk in which Zernike polynomials are orthogonal) = 1
        self.max_zd = max_zd

        # a list of functions to calculate big Z given rho, phi,
        # following eqn 1 of Thibos et al. (2002). The jth element of
        # the list returns the jth Z, following the indexing
        # convertions of Thibos et al. eqn 4.
        #
        # Should switch to using functools.cached_property in python 3.8
        self._z_function = self._build_z_functions()

        # A function that calculates the full Zernike approximation,
        # taking rho and phi as arguments.
        #
        # numexpr can only compile functions with a limited number of
        # arguments. If the order is too high, sum the terms
        # separately
        if order <= 7:
            self._zern_function = self._build_zern_function()
        else:
            self._zern_function = self._compute_sky_by_sum

        # big Z values for all m,n at all rho, phi in the
        # pre-defined healpix coordinate, following eqn 1 of Thibos et
        # al. (2002) The array returned should be indexed with j,
        # Should switch to using functools.cached_property in python 3.8
        self.healpix_z = self._compute_healpix_z()
        self._interpolate_healpix_z = interp1d(
            SIDEREAL_TIME_SAMPLES_RAD, self.healpix_z, axis=0, kind="nearest"
        )

        # A pd.DataFrame of zernike coeffs, indexed by mjd, providing the
        # Zernike polynomial coefficients for the approximation of the
        # sky at that time. That is, self._coeffs[5, 3] is the
        # j=3 coefficient of the approximation of the sky at
        # mjd=self.mjds[5], where j is defined as in Thibos et al. eqn 4.
        self._coeffs = pd.DataFrame()

    def load_coeffs(self, fname, band):
        """Load Zernike coefficients from a file.

        Parameters
        ----------
        fname : `str`
            The file name of the hdf5 file with the Zernike coeffs.
        band : `str`
            The band to load.

        """
        zernike_metadata = pd.read_hdf(fname, "zernike_metadata")
        assert self.order == zernike_metadata["order"]
        assert self.max_zd == zernike_metadata["max_zd"]
        all_zernike_coeffs = pd.read_hdf(fname, "zernike_coeffs")
        self._coeffs = all_zernike_coeffs.loc[band]
        self._coeff_calc_func = interp1d(
            self._coeffs.index.values, self._coeffs.values, axis=0
        )

    def compute_sky(self, alt, az, mjd=None):
        """Estimate sky values

        Parameters
        ----------
        alt : `np.ndarray`, (N)
            An array of altitudes above the horizon, in degrees
        az : `np.ndarray`, (N)
            An array of azimuth coordinates, in degrees
        mjd : `float`
            The time (floating point MJD) at which to estimate the sky.

        Returns
        -------
        `np.ndarray` (N) of sky brightnesses (mags/asec^2)
        """
        rho = self._calc_rho(alt)
        phi = self._calc_phi(az)
        result = self._zern_function(rho, phi, *tuple(self.coeffs(mjd)))
        return result

    def _compute_sky_by_sum(self, rho, phi, *coeffs):
        z = self._compute_z(rho, phi)
        if len(z.shape) == 2:
            result = np.sum(np.array(coeffs) * z, axis=1)
        else:
            result = np.sum(np.array(coeffs) * z)
        return result

    def compute_healpix(self, hpix, mjd=None):
        """Estimate sky values

        Parameters
        ----------
        hpix : `int`, (N)
            Array of healpix indexes of the desired coordinates.
        mjd : `float`
            The time (floating point MJD) at which to estimate the sky.

        Returns
        -------
        `np.ndarray` (N) of sky brightnesses (mags/asec^2)
        """
        interpolate_healpix_z = self._interpolate_healpix_z
        gmst = palpy.gmst(mjd)
        mjd_healpix_z = interpolate_healpix_z(gmst)
        # mjd_healpix_z = self.healpix_z[int(np.degrees(gmst))]
        if hpix is None:
            result = np.sum(self.coeffs(mjd) * mjd_healpix_z, axis=1)
        else:
            result = np.sum(self.coeffs(mjd) * mjd_healpix_z[hpix], axis=1)
        return result

    def coeffs(self, mjd):
        """Zerinke coefficients at a time

        Parameters
        ----------
        mjd : `float`
            The time (floating point MJD) at which to estimate the sky.

        Returns
        -------
        `np.ndarray` of Zernike coefficients following the OSA/ANSI
        indexing convention described in Thibos et al. (2002).
        """
        if len(self._coeffs) == 1:
            these_coeffs = self._coeffs.loc[mjd]
        else:
            calc_these_coeffs = self._coeff_calc_func
            these_coeffs = calc_these_coeffs(mjd)
        return these_coeffs

    def fit_coeffs(self, alt, az, sky, mjd, min_moon_sep=10, maxdiff=False):
        """Fit Zernike coefficients to a set of points

        Parameters
        ----------
        alt : `np.ndarray`, (N)
            An array of altitudes above the horizon, in degrees
        az : `np.ndarray`, (N)
            An array of azimuth coordinates, in degrees
        sky : `np.ndarray`, (N)
            An array of sky brightness values (mags/asec^2)
        mjd : `float`
            The time (floating point MJD) at which to estimate the sky.
        maxdiff : `bool`
            Minimize the maximum difference between the estimate and data,
            rather than the default RMS.

        """

        # Do not fit too close to the moon
        alt_rad, az_rad = np.radians(alt), np.radians(az)
        gmst_rad = palpy.gmst(mjd)
        lst_rad = gmst_rad + TELESCOPE.longitude_rad
        moon_ra_rad, moon_decl_rad, moon_diam = palpy.rdplan(
            mjd, 3, TELESCOPE.longitude_rad, TELESCOPE.latitude_rad
        )
        moon_ha_rad = lst_rad - moon_ra_rad
        moon_az_rad, moon_el_rad = palpy.de2h(
            moon_ha_rad, moon_decl_rad, TELESCOPE.latitude_rad
        )
        moon_sep_rad = palpy.dsepVector(
            np.full_like(az_rad, moon_az_rad),
            np.full_like(alt_rad, moon_el_rad),
            az_rad,
            alt_rad,
        )
        moon_sep = np.degrees(moon_sep_rad)

        rho = self._calc_rho(alt)
        phi = self._calc_phi(az)
        good_points = np.logical_and(rho <= 1.0, moon_sep > min_moon_sep)

        rho = rho[good_points]
        phi = phi[good_points]
        sky = sky[good_points]
        alt = alt[good_points]
        az = az[good_points]
        num_points = len(alt)
        assert len(az) == num_points
        assert len(sky) == num_points

        z = np.zeros((num_points, self._number_of_terms), dtype=self.dtype)
        for j in np.arange(self._number_of_terms):
            compute_z = self._z_function[j]
            z[:, j] = compute_z(rho, phi)

        # If the points being fit were evenly distributed across the sky,
        # we might be able to get away with a multiplication rather than
        # a linear regression, but we might be asked to fit masked data
        zern_fit = LinearRegression(fit_intercept=False).fit(z, sky)

        fit_coeffs = zern_fit.coef_

        if maxdiff:

            def max_abs_diff(test_coeffs):
                max_resid = np.max(
                    np.abs(np.sum(test_coeffs * z, axis=1) - sky)
                )
                return max_resid

            min_fit = scipy.optimize.minimize(max_abs_diff, fit_coeffs)
            fit_coeffs = min_fit.x

        self._coeffs = pd.DataFrame(
            [fit_coeffs],
            columns=np.arange(len(fit_coeffs)),
            index=pd.Index([mjd], name="mjd"),
        )
        return fit_coeffs

    def _compute_healpix_z(self):
        # Compute big Z values for all m,n at all rho, phi in the
        # pre-defined healpix coordinate, following eqn 1 of Thibos et
        # al. (2002) The array returned should be indexed with j,
        # following the conventions of eqn 4.
        sphere_npix = healpy.nside2npix(self.nside)
        sphere_ipix = np.arange(sphere_npix)
        ra, decl = healpy.pix2ang(self.nside, sphere_ipix, lonlat=True)

        num_st = len(SIDEREAL_TIME_SAMPLES_RAD)
        healpix_z = np.full(
            [num_st, sphere_npix, self._number_of_terms], np.nan
        )
        for st_idx, gmst_rad in enumerate(SIDEREAL_TIME_SAMPLES_RAD):
            lst_rad = gmst_rad + TELESCOPE.longitude_rad
            ha_rad = lst_rad - np.radians(ra)
            az_rad, alt_rad = palpy.de2hVector(
                ha_rad, np.radians(decl), TELESCOPE.latitude_rad
            )
            sphere_az, sphere_alt = np.degrees(az_rad), np.degrees(alt_rad)

            # We only need the half sphere above the horizen
            visible_ipix = sphere_ipix[sphere_alt > 0]
            alt, az = sphere_alt[visible_ipix], sphere_az[visible_ipix]
            rho = self._calc_rho(alt)
            phi = self._calc_phi(az)
            healpix_z[st_idx, visible_ipix] = self._compute_z(rho, phi)

        return healpix_z

    def _compute_horizan_healpix_z(self):
        # Compute big Z values for all m,n at all rho, phi in the
        # pre-defined healpix coordinate, following eqn 1 of Thibos et
        # al. (2002) The array returned should be indexed with j,
        # following the conventions of eqn 4.
        sphere_npix = healpy.nside2npix(self.nside)
        sphere_ipix = np.arange(sphere_npix)
        sphere_az, sphere_alt = healpy.pix2ang(
            self.nside, sphere_ipix, lonlat=True
        )

        # We only need the half sphere above the horizen
        ipix = sphere_ipix[sphere_alt > 0]
        alt, phi_deg = sphere_alt[ipix], sphere_az[ipix]
        rho = self._calc_rho(alt)
        rho, phi = (90.0 - alt) / self.max_zd, np.radians(phi_deg)
        healpix_z = self._compute_z(rho, phi)
        return healpix_z

    def _compute_z(self, rho, phi):
        # Compute big Z values for all m,n at rho, phi
        # following eqn 1 of Thibos et al. (2002)
        # The array returned should be indexed with j,
        # following the conventions of eqn 4.
        try:
            npix = len(rho)
            z = np.zeros((npix, self._number_of_terms), dtype=self.dtype)
            for j in np.arange(self._number_of_terms):
                compute_z = self._z_function[j]
                z[:, j] = compute_z(rho, phi)
        except TypeError:
            z = np.zeros(self._number_of_terms, dtype=self.dtype)
            for j in np.arange(self._number_of_terms):
                compute_z = self._z_function[j]
                z[j] = compute_z(rho, phi)
        return z

    def _build_z_functions(self):
        z_functions = []
        for j in np.arange(self._number_of_terms):
            z_functions.append(self._make_z_function(j))
        return z_functions

    def _build_zern_function(self):
        coeffs = [f"c{j}" for j in np.arange(self._number_of_terms)]

        expression = ""
        for j, coeff in enumerate(coeffs):
            zern_z_expr = self._make_z_expression(j)
            if zern_z_expr == "(1)":
                term = f"{coeff}"
            else:
                term = f"{coeff}*({zern_z_expr})"
            if expression == "":
                expression = term
            else:
                expression += f" + {term}"

        arg_types = []
        if expression.find("rho") >= 0:
            arg_types.append(("rho", self.dtype),)
        if expression.find("phi") >= 0:
            arg_types.append(("phi", self.dtype),)
        for coeff in coeffs:
            arg_types.append((coeff, self.dtype),)
        arg_types = tuple(arg_types)

        zern_function = NumExpr(expression, arg_types)
        return zern_function

    @property
    def _number_of_terms(self):
        n_terms = np.sum(np.arange(self.order) + 1)
        return n_terms

    def _make_r_expression(self, m, n):
        if (n - m) % 2 == 1:
            return 0
        assert n >= m
        assert m >= 0
        m = int(m)
        n = int(n)
        num_terms = 1 + (n - m) // 2
        expression = "("
        for k in range(num_terms):
            # From eqn 2 of Thibos et al. (2002)
            coeff = (((-1) ** k) * factorial(n - k)) / (
                factorial(k)
                * factorial(int((n + m) / 2 - k))
                * factorial(int((n - m) / 2 - k))
            )
            assert coeff == int(coeff)
            coeff = int(coeff)
            power = n - 2 * k
            if len(expression) > 1:
                expression += " + "
            if power == 0:
                expression += f"{coeff}"
            elif power == 1:
                expression += f"{coeff}*rho"
            else:
                expression += f"{coeff}*rho**{power}"
        expression += ")"
        return expression

    def _make_z_expression(self, j=None, mprime=None, n=None):
        if j is None:
            assert mprime is not None
            assert n is not None
        else:
            assert mprime is None
            assert n is None
            # From eqn 5 in Thibos et al. (2002)
            n = np.ceil((-3 + np.sqrt(9 + 8 * j)) / 2).astype(int)

            # From eqn 6 in Thibos et al. (2002)
            mprime = 2 * j - n * (n + 2)

        m = np.abs(mprime)
        r = self._make_r_expression(m, n)

        # From eqn. 3 of Thibos et al. 2002, again
        delta = 1 if m == 0 else 0
        big_nsq = 2 * (n + 1) / (1 + delta)
        assert int(big_nsq) == big_nsq
        big_nsq = int(big_nsq)

        if mprime == 0:
            expression = f"sqrt({big_nsq})*{r}"
        elif mprime > 0:
            expression = f"sqrt({big_nsq})*{r}*cos({m}*phi)"
        elif mprime < 0:
            expression = f"sqrt({big_nsq})*{r}*sin({m}*phi)"
        else:
            assert False

        return expression

    def _make_z_function(self, j=None, mprime=None, n=None):
        expression = self._make_z_expression(j, mprime, n)

        arg_types = []
        if expression.find("rho") >= 0:
            arg_types.append(("rho", self.dtype),)
        if expression.find("phi") >= 0:
            arg_types.append(("phi", self.dtype),)
        arg_types = tuple(arg_types)

        raw_z_function = NumExpr(expression, arg_types)

        # Create functions with dummy arguments so that
        # terms that do not require both phi and rho can
        # still accept them, such that all z_functions
        # can be called in the same way.
        if len(arg_types) == 0:

            def z_function(rho=None, phi=None):
                return raw_z_function()

        elif len(arg_types) == 1:

            def z_function(rho, phi=None):
                return raw_z_function(rho)

        else:
            z_function = raw_z_function

        return z_function

    def _calc_rho(self, alt):
        zd = 90.0 - alt
        if np.isscalar(alt) and zd > self.max_zd:
            return np.nan

        rho = zd / self.max_zd

        if not np.isscalar(alt):
            rho[zd > self.max_zd] = np.nan

        return rho

    def _calc_phi(self, az):
        phi = np.radians(az)
        return phi

class SkyBrightnessPreData:
    """Manager for raw pre-computed sky brightness data

    Parameters
    ----------
    base_fname : `str`
        Base name for data files to load.
    bands: `List` [`str`]
        Name of bands to read.
    pre_data_dir : `str`
        Name of source directory for pre-computed sky brightness data.
    max_num_mjds : `int`
        If there are more than this number of MJDs in the requested
        data files, sample this many out of the total.
    """

    def __init__(
        self, fname_base, bands, pre_data_dir=None, max_num_mjds=None
    ):
        if pre_data_dir is None:
            try:
                self.pre_data_dir = os.environ["SIMS_SKYBRIGHTNESS_DATA"]
            except KeyError:
                self.pre_data_dir = "."
        else:
            self.pre_data_dir = pre_data_dir

        self.fname_base = fname_base
        self.max_num_mjds = max_num_mjds
        self.times = None
        self.sky = None
        self.metadata = {}
        self.load(fname_base, bands)

    def load(self, fname_base, bands="ugrizy"):
        """Load pre-computed sky values.

        Parameters
        ----------
        base_fname : `str`
            Base name for data files to load.
        bands: `List` [`str`]
            Name of bands to read.
        """

        npz_fname = os.path.join(self.pre_data_dir, fname_base + "." + "npz")
        npy_fname = os.path.join(self.pre_data_dir, fname_base + "." + "npy")
        npz = np.load(npz_fname, allow_pickle=True)
        npz_hdr = npz["header"][()]
        npz_data = npz["dict_of_lists"][()]
        pre_sky = np.load(npy_fname, allow_pickle=True)

        alt = npz_hdr["alt"]
        az = npz_hdr["az"]
        alt_rad, az_rad = np.radians(alt), np.radians(az)

        self.metadata = npz_hdr

        self.times = pd.DataFrame(
            {
                k: npz_data[k]
                for k in npz_data.keys()
                if npz_data[k].shape == npz_data["mjds"].shape
            }
        )

        read_mjds = len(self.times)
        if self.max_num_mjds is not None:
            read_mjd_idxs = pd.Series(np.arange(read_mjds))
            mjd_idxs = read_mjd_idxs.sample(self.max_num_mjds)
        else:
            mjd_idxs = np.arange(read_mjds)

        skies = []
        for mjd_idx in mjd_idxs:
            mjd = npz_data["mjds"][mjd_idx]
            gmst_rad = palpy.gmst(mjd)
            lst_rad = gmst_rad + TELESCOPE.longitude_rad
            ha_rad, decl_rad = palpy.dh2eVector(
                az_rad, alt_rad, TELESCOPE.latitude_rad
            )
            ra_rad = (lst_rad - ha_rad) % (2 * np.pi)
            moon_ra_rad = npz_data["moonRAs"][mjd_idx]
            moon_decl_rad = npz_data["moonDecs"][mjd_idx]
            moon_ha_rad = lst_rad - moon_ra_rad
            moon_az_rad, moon_el_rad = palpy.de2h(
                moon_ha_rad, moon_decl_rad, TELESCOPE.latitude_rad
            )
            moon_sep = palpy.dsepVector(
                np.full_like(az_rad, moon_az_rad),
                np.full_like(alt_rad, moon_el_rad),
                az_rad,
                alt_rad,
            )
            for band in bands:
                skies.append(
                    pd.DataFrame(
                        {
                            "band": band,
                            "mjd": npz_data["mjds"][mjd_idx],
                            "gmst": np.degrees(gmst_rad),
                            "lst": np.degrees(lst_rad),
                            "alt": alt,
                            "az": az,
                            "ra": np.degrees(ra_rad),
                            "decl": np.degrees(decl_rad),
                            "moon_ra": np.degrees(
                                npz_data["moonRAs"][mjd_idx]
                            ),
                            "moon_decl": np.degrees(
                                npz_data["moonDecs"][mjd_idx]
                            ),
                            "moon_alt": np.degrees(
                                npz_data["moonAlts"][mjd_idx]
                            ),
                            "moon_az": np.degrees(moon_az_rad),
                            "moon_sep": np.degrees(moon_sep),
                            "sun_ra": np.degrees(npz_data["sunRAs"][mjd_idx]),
                            "sun_decl": np.degrees(
                                npz_data["sunDecs"][mjd_idx]
                            ),
                            "sun_alt": np.degrees(
                                npz_data["sunAlts"][mjd_idx]
                            ),
                            "sky": pre_sky[band][mjd_idx],
                        }
                    )
                )

        self.sky = pd.concat(skies).set_index(
            ["band", "mjd", "alt", "az"], drop=False
        )
        self.sky.sort_index(inplace=True)

        if self.max_num_mjds is not None:
            self.times = self.times.iloc[mjd_idxs]

    def __getattr__(self, name):
        return self.metadata[name]

class SkyModelZernike:
    """Interface to zernike sky that is more similar to SkyModelPre

    Parameters
    ----------
    data_file : `str`, optional
        File name from which to load Zernike coefficients. Default None uses default data directory.

    """

    def __init__(self, data_file=None, **kwargs):
        if data_file is None:
            if "SIMS_SKYBRIGHTNESS_DATA" in os.environ:
                data_dir = os.environ["SIMS_SKYBRIGHTNESS_DATA"]
            else:
                data_dir = os.path.join(get_data_dir(),
                                        "sims_skybrightness_pre")

            data_file = os.path.join(data_dir, "zernike", "zernike.h5")

            zernike_metadata = pd.read_hdf(data_file, "zernike_metadata")

            order = int(zernike_metadata["order"])
            if 'order' in kwargs:
                assert order == kwargs['order']
            else:
                kwargs['order'] = order

            max_zd = zernike_metadata["max_zd"]
            if 'max_zd' in kwargs:
                assert max_zd == kwargs['max_zd']
            else:
                kwargs['max_zd'] = max_zd

        self.zernike_model = {}
        for band in BANDS:
            sky = ZernikeSky(**kwargs)
            sky.load_coeffs(data_file, band)
            self.zernike_model[band] = sky
        self.nside = sky.nside

    def returnMags(
        self,
        mjd,
        indx=None,
        badval=healpy.UNSEEN,
        filters=["u", "g", "r", "i", "z", "y"],
        extrapolate=False,
    ):
        """
        Return a full sky map or individual pixels for the input mjd

        Parameters
        ----------
        mjd : float
            Modified Julian Date to interpolate to
        indx : List of int(s) (None)
            indices to interpolate the sky values at. Returns full sky if None. If the class was
            instatiated with opsimFields, indx is the field ID, otherwise it is the healpix ID.
        badval : float (-1.6375e30)
            Mask value. Defaults to the healpy mask value.
        filters : list
            List of strings for the filters that should be returned.
        extrapolate : bool (False)
            In indx is set, extrapolate any masked pixels to be the same as the nearest non-masked
            value from the full sky map.

        Returns
        -------
        sbs : dict
            A dictionary with filter names as keys and np.arrays as values which
            hold the sky brightness maps in mag/sq arcsec.
        """
        sky_brightness = {}

        sun_el = _calc_sun_el(mjd)
        if sun_el > 0:
            warnings.warn('Requested MJD between sunrise and sunset')
            if indx is None:
                nside = self.zernike_model[filters[0]].nside
                npix = healpy.nside2npix(nside)
            else:
                npix = len(indx)

            for band in filters:
                sky_brightness[band] = np.full(npix, badval)

            return sky_brightness

        if extrapolate:
            raise NotImplementedError

        for band in filters:
            band_brightness = self.zernike_model[band].compute_healpix(
                indx, mjd
            )
            badval_idxs = np.where(~np.isfinite(band_brightness))
            band_brightness[badval_idxs] = badval
            sky_brightness[band] = band_brightness

        return sky_brightness

def cut_pre_dataset(
    fname_base="59823_60191",
    num_mjd=3,
    pre_dir="/data/des91.b/data/neilsen/LSST/horizon_skybrightness_pre",
    cut_dir=".",
):
    """Cut a per-computed dataset to specified number of MJDs

    The purpose of this command is to create small input datafiles
    that can be used for testing

    Parameters
    ----------
    fname_base : `str`
        The base from which to construct file names
    num_mjd : `int`
        The number of MJDs to include in the cut file.
    pre_dir : `str`
        The directory from which to read data files
    cut_dur : `str`
        The directory into which to write cut data files
    """
    npy_fname = os.path.join(pre_dir, fname_base + ".npy")
    npz_fname = os.path.join(pre_dir, fname_base + ".npz")

    npz = np.load(npz_fname, allow_pickle=True)
    npz_hdr = npz["header"][()]
    npz_data = npz["dict_of_lists"][()]
    pre_sky = np.load(npy_fname, allow_pickle=True)

    kept_mjds = np.sort(npz_hdr["required_mjds"])[:num_mjd].copy()
    kept_mjd_idxs = np.where(np.in1d(npz_data["mjds"], kept_mjds))

    npz_hdr["required_mjds"] = kept_mjds
    del npz_hdr["version"]
    del npz_hdr["fingerprint"]

    for data_key in npz_data.keys():
        npz_data[data_key] = npz_data[data_key][kept_mjd_idxs].copy()

    pre_sky = pre_sky[kept_mjd_idxs].copy()

    min_mjd = int(np.floor(kept_mjds.min()))
    max_mjd = int(np.floor(kept_mjds.max()))
    out_fname_base = f"{min_mjd}_{max_mjd}"

    cut_npz_fname = os.path.join(cut_dir, out_fname_base + ".npz")
    np.savez(cut_npz_fname, header=npz_hdr, dict_of_lists=npz_data)

    cut_npy_fname = os.path.join(cut_dir, out_fname_base + ".npy")
    np.save(cut_npy_fname, pre_sky)

# internal functions & classes

@lru_cache()
def _calc_moon_az_rad(mjd):
    ra_rad, decl_rad, diam = palpy.rdplan(
        mjd, 3, TELESCOPE.longitude_rad, TELESCOPE.latitude_rad
    )
    ha_rad = palpy.gmst(mjd) + TELESCOPE.longitude_rad - ra_rad
    az_rad, el_rad = palpy.de2h(ha_rad, decl_rad, TELESCOPE.latitude_rad)
    return az_rad

@lru_cache()
def _calc_sun_el(mjd):
    ra_rad, decl_rad, diam = palpy.rdplan(
        mjd, 0, TELESCOPE.longitude_rad, TELESCOPE.latitude_rad
    )
    ha_rad = palpy.gmst(mjd) + TELESCOPE.longitude_rad - ra_rad
    az_rad, el_rad = palpy.de2h(ha_rad, decl_rad, TELESCOPE.latitude_rad)
    el = np.degrees(el_rad)
    return el

ext: py

max_stars_repo_name: RileyWClarke/flarubin

max_issues_repo_name: RileyWClarke/flarubin

max_forks_repo_name: RileyWClarke/flarubin

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


text:
# Copyright (c) 2020 Microsoft
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

from azureml.core import Run
import argparse
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os

from sklearn.model_selection import train_test_split
from sklearn.linear_model import ElasticNet, Ridge
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import joblib

import warnings
warnings.filterwarnings("ignore", message='l1_reg="auto" is deprecated')
warnings.filterwarnings(
    "ignore", message='Using a non-tuple sequence for multidimensional indexing is deprecated')

# Create logger from current context
# If this is run in AML Compute, it will return the current run
# If this runs locally, it'll only print metrics to standard out
run = Run.get_context()

def load_train_test_dataset(dataset_name):
    # Load DF from CSV
    diabetes_df = run.input_datasets[dataset_name].to_pandas_dataframe()

    # Split out X and Y variables
    y = diabetes_df.pop('target').values
    X = diabetes_df

    column_names = X.columns

    # Split training and test datasets
    X_train, X_test, y_train, y_test = train_test_split(
        X.values, y, test_size=0.2, random_state=0)

    print(
        f"Data contains {len(X_train)} training samples and {len(X_test)} test samples")
    return X_train, X_test, y_train, y_test, column_names

def train_elasticnet(X, y, alpha, l1_ratio):
    run.log('model_type', 'ElasticNet')

    # Create a new model object
    model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio,
                       copy_X=True, random_state=40760)

    # Fit the model and return
    model.fit(X, y)

    # Save model to "outputs" folder
    os.makedirs('./outputs', exist_ok=True)
    joblib.dump(value=model, filename='./outputs/model.pkl')

    return model

def train_ridge_regression(X, y, alpha):
    run.log('model_type', "Ridge")

    # Create a new model object
    model = Ridge(alpha=alpha)

    # Fit the model and return
    model.fit(X, y)

    # Save model to "outputs" folder
    os.makedirs('./outputs', exist_ok=True)
    joblib.dump(value=model, filename='./outputs/model.pkl')

    return model

def train_gradient_boosted_regressor(X, y, alpha):
    run.log('model_type', "GradientBoostingRegressor")

    model = GradientBoostingRegressor(loss='huber', alpha=alpha)

    # Fit the model and return
    model.fit(X, y)

    # Save model to "outputs" folder
    os.makedirs('./outputs', exist_ok=True)
    joblib.dump(value=model, filename='./outputs/model.pkl')

    return model

def predict_and_log_performance(model, X_test, y_test):
    # Get the predicted values
    preds = model.predict(X_test)

    # Log the metrics to the AML run
    run.log("rmse", np.sqrt(mean_squared_error(y_test, preds)))
    run.log('mae', mean_absolute_error(y_test, preds))
    run.log('r2', r2_score(y_test, preds))

    return preds

def plot_residuals_v_actuals(y, y_hat):
    """Residuals (y-axis) vs. Actuals (x-axis) - colored green"""
    resids = y - y_hat

    fig = plt.figure()
    sns.regplot(y, resids, color='g')

    plt.title('Residual vs. Actual')
    plt.xlabel("Actual Value")
    plt.ylabel("Residuals")

    plt.close(fig)

    run.log_image(name='residuals-v-actuals', plot=fig)
    return fig

def plot_predictions(y, y_hat):
    """Predictions (y-axis) vs. Actuals (x-axis)"""
    fig = plt.figure()

    sns.regplot(y, y_hat, color='b')

    plt.title("Prediction vs. Actual")
    plt.xlabel("Actual Value")
    plt.ylabel("Predicted Value")

    plt.close(fig)
    run.log_image(name='prediction-v-actual', plot=fig)
    return fig

def plot_resid_histogram(y, y_hat):
    resids = y - y_hat

    fig = plt.figure()
    sns.distplot(resids, color='g')

    plt.title("Residual Histogram")

    plt.close(fig)

    run.log_image(name='residuals-histogram', plot=fig)
    return fig

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument('--dataset-name', '-d', type=str, default='data',
                        help="The path where the data file is located")
    parser.add_argument('--alpha', '-a', type=float,
                        help="The alpha value for training", default=0.03)
    parser.add_argument('--file-name', '-f', type=str, default='diabetes.csv',
                        help="The file name of the diabetes csv dataset")
    parser.add_argument('--l1-ratio', type=float, default=0.05,
                        help='The l1_ratio of the Scikit-Learn ElasticNet model')
    parser.add_argument('--learning-rate', '-l', type=float, default=0.001,
                        help='The learning rate for the GBT model')
    parser.add_argument('--model-name', '-n', type=str, default="ridge",
                        help="The name of the model to try. Supported values are 'ridge', 'elastic', and 'gbt'")
    args = parser.parse_args()

    X_train, X_test, y_train, y_test, column_names = load_train_test_dataset(
        args.dataset_name)

    if args.model_name == "ridge":
        # log the hyperparameters
        run.log('alpha', args.alpha)

        # set the arguments and algorithm choice
        args = {'alpha': args.alpha}
        algo = train_ridge_regression

    elif args.model_name == "elastic":
        # log the hyperparameters
        run.log('alpha', args.alpha)
        run.log('l1_ratio', args.l1_ratio)

        # set the arguments and algorithm choice
        args = {'alpha': args.alpha,
                'l1_ratio': args.l1_ratio}
        algo = train_elasticnet

    elif args.model_name == 'gbt':
        # log the hyperparameters
        run.log('alpha', args.alpha)
        run.log('l1_ratio', args.l1_ratio)

        # set the arguments and algorithm choice
        args = {'alpha': args.alpha}
        algo = train_gradient_boosted_regressor

    # Train the model
    model = algo(X=X_train,
                 y=y_train,
                 **args)

    # Generate predictions and log performance metrics
    preds = predict_and_log_performance(model=model,
                                        X_test=X_test,
                                        y_test=y_test)

    # Plot the residuals
    resid_fig = plot_residuals_v_actuals(y_test, preds)
    resid_hist = plot_resid_histogram(y_test, preds)
    pred_plt = plot_predictions(y_test, preds)

max_stars_repo_name: ezwiefel/aml-workshop

max_issues_repo_name: ezwiefel/aml-workshop

max_forks_repo_name: ezwiefel/aml-workshop

reason: import numpy


text:
import os
import time
import load_data
import torch
import torch.nn.functional as F
from torch.autograd import Variable
import torch.optim as optim
import numpy as np
from models.LSTM import LSTMClassifier

class TextClassificationModel(object):

    def __init__(self, model):
        self._model = model
        pass

    def _clip_gradient(self, clip_value):
        params = list(filter(lambda p: p.grad is not None, self._model.parameters()))
        for p in params:
            p.grad.data.clamp_(-clip_value, clip_value)

    def train(self, train_iter, epoch, loss_fn = F.cross_entropy,
        learning_rate = 2e-5,
            batch_size = 32,
            output_size = 2,
            hidden_size = 256,
            embedding_length = 300):
        total_epoch_loss = 0
        total_epoch_acc = 0
        self._model.cuda()
        optim = torch.optim.Adam(filter(lambda p: p.requires_grad, self._model.parameters()))
        steps = 0
        self._model.train()
        for idx, batch in enumerate(train_iter):
            text = batch.text[0]
            target = batch.label
            target = torch.autograd.Variable(target).long()
            if torch.cuda.is_available():
                text = text.cuda()
                target = target.cuda()
            if (text.size()[0] is not 32):# One of the batch returned by BucketIterator has length different than 32.
                continue
            optim.zero_grad()
            prediction = self._model(text)
            loss = loss_fn(prediction, target)
            num_corrects = (torch.max(prediction, 1)[1].view(target.size()).data == target.data).float().sum()
            acc = 100.0 * num_corrects/len(batch)
            loss.backward()
            self._clip_gradient(1e-1)
            optim.step()
            steps += 1

            if steps % 100 == 0:
                print (f'Epoch: {epoch+1}, Idx: {idx+1}, Training Loss: {loss.item():.4f}, Training Accuracy: {acc.item(): .2f}%')

            total_epoch_loss += loss.item()
            total_epoch_acc += acc.item()

        return total_epoch_loss/len(train_iter), total_epoch_acc/len(train_iter)

    def eval(self, val_iter, score_fn ):
        total_epoch_loss = 0
        total_epoch_acc = 0
        self._model.eval()
        with torch.no_grad():
            for _, batch in enumerate(val_iter):
                text = batch.text[0]
                if (text.size()[0] is not 32):
                    continue
                target = batch.label
                target = torch.autograd.Variable(target).long()
                if torch.cuda.is_available():
                    text = text.cuda()
                    target = target.cuda()
                prediction = self._model(text)
                loss = score_fn(prediction, target)
                num_corrects = (torch.max(prediction, 1)[1].view(target.size()).data == target.data).sum()
                acc = 100.0 * num_corrects/len(batch)
                total_epoch_loss += loss.item()
                total_epoch_acc += acc.item()

        return total_epoch_loss/len(val_iter), total_epoch_acc/len(val_iter)


text:
from __future__ import print_function
from keras.preprocessing.image import ImageDataGenerator
import numpy as np
import os
import glob
import skimage.io as io
import skimage.transform as trans
import itertools

Sky = [128,128,128]
Building = [128,0,0]
Pole = [192,192,128]
Road = [128,64,128]
Pavement = [60,40,222]
Tree = [128,128,0]
SignSymbol = [192,128,128]
Fence = [64,64,128]
Car = [64,0,128]
Pedestrian = [64,64,0]
Bicyclist = [0,128,192]
Unlabelled = [0,0,0]

COLOR_DICT = np.array([Sky, Building, Pole, Road, Pavement,
                          Tree, SignSymbol, Fence, Car, Pedestrian, Bicyclist, Unlabelled])

def adjustData(img,mask,flag_multi_class,num_class):
    if(flag_multi_class):
        # img = img / 255
        img /= 255.0
        mask = mask[:,:,:,0] if(len(mask.shape) == 4) else mask[:,:,0]
        new_mask = np.zeros(mask.shape + (num_class,))
        for i in range(num_class):
            #for one pixel in the image, find the class in mask and convert it into one-hot vector
            #index = np.where(mask == i)
            #index_mask = (index[0],index[1],index[2],np.zeros(len(index[0]),dtype = np.int64) + i) if (len(mask.shape) == 4) else (index[0],index[1],np.zeros(len(index[0]),dtype = np.int64) + i)
            #new_mask[index_mask] = 1
            new_mask[mask == i,i] = 1
        new_mask = np.reshape(new_mask,(new_mask.shape[0],new_mask.shape[1]*new_mask.shape[2],new_mask.shape[3])) if flag_multi_class else np.reshape(new_mask,(new_mask.shape[0]*new_mask.shape[1],new_mask.shape[2]))
        mask = new_mask
    elif(np.max(img) > 1):
        # img = img / 255
        img /= 255.0
        # mask = mask / 255
        mask /= 255.0
        mask[mask > 0.5] = 1
        mask[mask <= 0.5] = 0
    return (img,mask)

def trainGenerator(batch_size,train_path,image_folder,mask_folder,aug_dict,image_color_mode = "grayscale",
                    mask_color_mode = "grayscale",image_save_prefix  = "image",mask_save_prefix  = "mask",
                    flag_multi_class = False,num_class = 2,save_to_dir = None,target_size = (512,512),seed = 1):
    '''
    can generate image and mask at the same time
    use the same seed for image_datagen and mask_datagen to ensure the transformation for image and mask is the same
    if you want to visualize the results of generator, set save_to_dir = "your path"
    '''
    image_datagen = ImageDataGenerator(**aug_dict)
    mask_datagen = ImageDataGenerator(**aug_dict)
    image_generator = image_datagen.flow_from_directory(
        train_path,
        classes = [image_folder],
        class_mode = None,
        color_mode = image_color_mode,
        target_size = target_size,
        batch_size = batch_size,
        save_to_dir = save_to_dir,
        save_prefix  = image_save_prefix,
        seed = seed)
    mask_generator = mask_datagen.flow_from_directory(
        train_path,
        classes = [mask_folder],
        class_mode = None,
        color_mode = mask_color_mode,
        target_size = target_size,
        batch_size = batch_size,
        save_to_dir = save_to_dir,
        save_prefix  = mask_save_prefix,
        seed = seed)
    # train_generator = zip(image_generator, mask_generator)
    train_generator = itertools.izip(image_generator, mask_generator)
    for (img,mask) in train_generator:
        img,mask = adjustData(img,mask,flag_multi_class,num_class)
        yield (img,mask)

def testGenerator(test_path,num_image = 223, target_size = (512,512),flag_multi_class = False,as_gray = True):
    for i in range(num_image):
        img = io.imread(os.path.join(test_path,"%d.jpg"%i),as_gray = True)

        img = img / 255.
        # img /= 255.

        img = trans.resize(img,target_size)
        img = np.reshape(img,img.shape+(1,)) if (not flag_multi_class) else img
        img = np.reshape(img,(1,)+img.shape)
        # img = (img * 255).astype(np.uint8)
        yield img

def geneTrainNpy(image_path,mask_path,flag_multi_class = False,num_class = 2,image_prefix = "image",mask_prefix = "mask",image_as_gray = True,mask_as_gray = True):
    image_name_arr = glob.glob(os.path.join(image_path,"%s*.png"%image_prefix))
    image_arr = []
    mask_arr = []
    for index,item in enumerate(image_name_arr):
        img = io.imread(item,as_gray = image_as_gray)
        img = np.reshape(img,img.shape + (1,)) if image_as_gray else img
        mask = io.imread(item.replace(image_path,mask_path).replace(image_prefix,mask_prefix),as_gray = mask_as_gray)
        mask = np.reshape(mask,mask.shape + (1,)) if mask_as_gray else mask
        img,mask = adjustData(img,mask,flag_multi_class,num_class)
        image_arr.append(img)
        mask_arr.append(mask)
    image_arr = np.array(image_arr)
    mask_arr = np.array(mask_arr)
    return image_arr,mask_arr

def labelVisualize(num_class,color_dict,img):
    img = img[:,:,0] if len(img.shape) == 3 else img
    img_out = np.zeros(img.shape + (3,))
    for i in range(num_class):
        img_out[img == i,:] = color_dict[i]
    return img_out / 255.0

def saveResult(save_path,npyfile,flag_multi_class = False,num_class = 2):
    for i,item in enumerate(npyfile):
        img = labelVisualize(num_class,COLOR_DICT,item) if flag_multi_class else item[:,:,0]
        io.imsave(os.path.join(save_path,"%d.jpg"%i),img)


text:
"""
The *utils* module contains utilities not strictly related to visualization which we often use (eg harmonics computation).
"""

import math
from datetime import datetime
from os.path import dirname

import numpy as np
import pandas as pd
from colour import Color
from scipy import stats
from statsmodels.tsa._stl import STL
from statsmodels.tsa.seasonal import seasonal_decompose

def harmonics(
    dates: pd.DataFrame,
    period: int,
    n: int,
    epoch: datetime = datetime(1900, 1, 1),
) -> pd.DataFrame:
    """
    Computes harmonics for the given dates. Each harmonic is made of a couple of sinusoidal and cosinusoidal waves
    with frequency i/period, i = 1...n. The argument of the functions is the number of hours from the starting epoch.

    :param dates: A pandas series of dates.
    :param period: The base period of the harmonics.
    :param n: The number of harmonics to include.
    :param epoch: The epoch used to compute the argument of the sin.
    :return: A Pandas DataFrame with dates as index and harmonics as columns.
    """
    d = pd.DataFrame(index=dates)
    hours = (dates - epoch) / pd.Timedelta(hours=1)

    for i in range(1, n + 1):
        d[f"Sin_{round(period)}_{i}"] = np.sin(2 * i * np.pi * hours / period)
        d[f"Cos_{round(period)}_{i}"] = np.cos(2 * i * np.pi * hours / period)

    return d

def boxcox(x):
    return stats.boxcox(x)[0]

def yeojohnson(x):
    return stats.yeojohnson(x)[0]

def moving_average(x, w):
    return np.convolve(x, np.ones(w), "same") / w

def set_time_index(df, time_col):
    if time_col is not None:
        df.set_index(time_col, inplace=True)

def get_components(result):
    components = {}
    for c in valid_components:
        if c in dir(result):
            components[c] = getattr(result, c)
    return components

def make_granular_colorscale(col1, col2, n):
    return [c.hex for c in Color(col1).range_to(col2, n)]

def apply_grad_color_to_traces(fig, col1, col2):
    n_traces = len(fig["data"])
    custom_colorscale = make_granular_colorscale(col1, col2, n_traces)
    for i in range(n_traces):
        fig["data"][i]["line"]["color"] = custom_colorscale[i]

root_path = dirname(dirname(__file__))

decomp_methods = {
    "STL": {STL: {}},
    "seasonal_additive": {seasonal_decompose: {"model": "additive"}},
    "seasonal_multiplicative": {
        seasonal_decompose: {"model": "multiplicative"}
    },
}

valid_components = [
    "level",
    "trend",
    "seasonal",
    "freq_seasonal",
    "cycle",
    "autoregressive",
    "resid",
]

valid_seasons = {
    "grouping": {
        "minute": lambda x: x.minute,
        "hour": lambda x: x.hour,
        "day": lambda x: x.isocalendar().day,
        "week": lambda x: x.isocalendar().week,
        "month": lambda x: x.month,
        "quarter": lambda x: (x.month - 1) // 3 + 1,
        "year": lambda x: x.isocalendar().year,
    },
    "granularity": {
        "minute": lambda x: x.second,
        "hour": lambda x: x.minute,
        "day": lambda x: x.hour,
        "week": lambda x: x.isocalendar().day,
        "month": lambda x: x.day,
        "quarter": lambda x: (x - pd.PeriodIndex(x, freq="Q").start_time).days
        + 1,
        "year": lambda x: x.dayofyear,
    },
}

transform_dict = {
    "Box-Cox": boxcox,
    "Yeo-Johnson": yeojohnson,
    "log": np.vectorize(math.log),
    "moving_average": moving_average,
}

max_stars_repo_name: xtreamsrl/pytsviz

max_issues_repo_name: xtreamsrl/pytsviz

max_forks_repo_name: xtreamsrl/pytsviz

reason: import numpy,from scipy,from statsmodels


text:
#importing the required packages
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt

from tensorflow.examples.tutorials.mnist import input_data
mnist_dataset = input_data.read_data_sets('MNIST_data', validation_size=0)

# Plotting one image from the training set.
image = mnist_dataset.train.images[2]
plt.imshow(image.reshape((28, 28)), cmap='Greys_r')

learning_rate = 0.001

# Define the placeholder variable sfor the input and target values
inputs_values = tf.placeholder(tf.float32, (None, 28,28,1), name="inputs_values")
targets_values = tf.placeholder(tf.float32, (None, 28,28,1), name="targets_values")

# Defining the Encoder part of the netowrk
# Defining the first convolution layer in the encoder parrt
# The output tenosor will be in the shape of 28x28x16
conv_layer_1 = tf.layers.conv2d(inputs=inputs_values, filters=16, kernel_size=(3,3), padding='same', activation=tf.nn.relu)

# The output tenosor will be in the shape of 14x14x16
maxpool_layer_1 = tf.layers.max_pooling2d(conv_layer_1, pool_size=(2,2), strides=(2,2), padding='same')

# The output tenosor will be in the shape of 14x14x8
conv_layer_2 = tf.layers.conv2d(inputs=maxpool_layer_1, filters=8, kernel_size=(3,3), padding='same', activation=tf.nn.relu)

# The output tenosor will be in the shape of 7x7x8
maxpool_layer_2 = tf.layers.max_pooling2d(conv_layer_2, pool_size=(2,2), strides=(2,2), padding='same')

# The output tenosor will be in the shape of 7x7x8
conv_layer_3 = tf.layers.conv2d(inputs=maxpool_layer_2, filters=8, kernel_size=(3,3), padding='same', activation=tf.nn.relu)

# The output tenosor will be in the shape of 4x4x8
encoded_layer = tf.layers.max_pooling2d(conv_layer_3, pool_size=(2,2), strides=(2,2), padding='same')

# Defining the Decoder part of the netowrk
# Defining the first upsampling layer in the decoder part
# The output tenosor will be in the shape of 7x7x8
upsample_layer_1 = tf.image.resize_images(encoded_layer, size=(7,7), method=tf.image.ResizeMethod.NEAREST_NEIGHBOR)

# The output tenosor will be in the shape of 7x7x8
conv_layer_4 = tf.layers.conv2d(inputs=upsample_layer_1, filters=8, kernel_size=(3,3), padding='same', activation=tf.nn.relu)

# The output tenosor will be in the shape of 14x14x8
upsample_layer_2 = tf.image.resize_images(conv_layer_4, size=(14,14), method=tf.image.ResizeMethod.NEAREST_NEIGHBOR)

# The output tenosor will be in the shape of 14x14x8
conv_layer_5 = tf.layers.conv2d(inputs=upsample_layer_2, filters=8, kernel_size=(3,3), padding='same', activation=tf.nn.relu)

# The output tenosor will be in the shape of 28x28x8
upsample_layer_3 = tf.image.resize_images(conv_layer_5, size=(28,28), method=tf.image.ResizeMethod.NEAREST_NEIGHBOR)

# The output tenosor will be in the shape of 28x28x16
conv6 = tf.layers.conv2d(inputs=upsample_layer_3, filters=16, kernel_size=(3,3), padding='same', activation=tf.nn.relu)

# The output tenosor will be in the shape of 28x28x1
logits_layer = tf.layers.conv2d(inputs=conv6, filters=1, kernel_size=(3,3), padding='same', activation=None)

# feeding the logits values to the sigmoid activation function to get the reconstructed images
decoded_layer = tf.nn.sigmoid(logits_layer)

# feeding the logits to sigmoid while calculating the cross entropy
model_loss = tf.nn.sigmoid_cross_entropy_with_logits(labels=targets_values, logits=logits_layer)

# Getting the model cost and defining the optimizer to minimize it
model_cost = tf.reduce_mean(model_loss)
model_optimizer = tf.train.AdamOptimizer(learning_rate).minimize(model_cost)

# Starting the training process
sess = tf.Session()

num_epochs = 20
train_batch_size = 200
sess.run(tf.global_variables_initializer())

# Visualizing the results
for e in range(num_epochs):
    for ii in range(mnist_dataset.train.num_examples//train_batch_size):
        input_batch = mnist_dataset.train.next_batch(train_batch_size)
        input_images = input_batch[0].reshape((-1, 28, 28, 1))
        input_batch_cost, _ = sess.run([model_cost, model_optimizer], feed_dict={inputs_values: input_images,
                                                         targets_values: input_images})

        print("Epoch: {}/{}...".format(e+1, num_epochs),
              "Training loss: {:.3f}".format(input_batch_cost))

fig, axes = plt.subplots(nrows=2, ncols=10, sharex=True, sharey=True, figsize=(20,4))
input_images = mnist_dataset.test.images[:10]
reconstructed_images = sess.run(decoded_layer, feed_dict={inputs_values: input_images.reshape((10, 28, 28, 1))})

for imgs, row in zip([input_images, reconstructed_images], axes):
    for img, ax in zip(imgs, row):
        ax.imshow(img.reshape((28, 28)), cmap='Greys_r')
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)

fig.tight_layout(pad=0.1)


text:
import numpy as np
import os
import fnmatch

#center shape by removing the mean of each column
def centerShape(shape):
    shape=shape-shape.mean(axis=0)
    return shape

#scale shape by divinding by frobinius norm
def scaleShape(shape):
    shape=shape/np.linalg.norm(shape)
    return shape

#align one shape to a reference shape, soley by rotation
def alignShape(refShape,shape):
    """
    Align the shape to the mean.
    """
    u,s,v=np.linalg.svd(np.dot(np.transpose(refShape),shape), full_matrices=True)
    rotationMatrix=np.dot(np.transpose(v), np.transpose(u))
    shape=np.dot(shape,rotationMatrix)

    return shape

def procrustesAlign(refShape,shape):
   #center both shapes
    refShape=centerShape(refShape)
    shape=centerShape(shape)

    #scale both shapes
    refShape=scaleShape(refShape)
    shape=scaleShape(shape)

    # rotate shape to match the refshape
    shape=alignShape(refShape,shape)

    return shape

def procrustesAlignNoScale(refShape,shape):
   #center both shapes
    refShape=centerShape(refShape)
    shape=centerShape(shape)
    shape=alignShape(refShape,shape)

    return shape

# <codecell>

def makeTwoDim(monsters):
    i,j,k=monsters.shape
    tmp=np.zeros((i*j,k))
    for x in range(k):
        vec=np.reshape(monsters[:,:,x],(i*j),order='F')
        tmp[:,x]=vec
    return tmp

def calcMean(vec):
    i,j=vec.shape
    meanVec=np.zeros((i))
    for x in range(j):
        meanVec+=vec[:,x]/float(j)
    return meanVec

def calcCov(vec):
    i,j=vec.shape
    meanVec=calcMean(vec)
    covMatrix=np.zeros((i,i))
    for x in range(j):
        t1=((vec[:,x]-meanVec).T).reshape(i,1)
        t2=(vec[:,x]-meanVec).reshape(1,i)
        covMatrix+=np.dot(t1,t2)/float(j)
    return covMatrix

def sortEig(eVal, eVec):
    i,j=eVec.shape
    ePair=list(range(j))
    for y in range(j):
        ePair[y]=[(np.abs(eVal[y]), eVec[:,y])]
        ePair[y].sort()
        ePair[y].reverse()
    return ePair

def makeTransformMatrix(ePair,pcA,pcB):
    tmp=ePair[0][0][1]
    i=tmp.shape
    i=i[0]
    vec1=ePair[pcA][0][1]
    vec2=ePair[pcB][0][1]

    transform=np.zeros((i,2))
    transform[:,0]=vec1.reshape(i).real
    transform[:,1]=vec2.reshape(i).real
    return transform

def plotTanProj(monsters,pcA,pcB):
    twoDim=makeTwoDim(monsters)
    mShape=calcMean(twoDim)
    covMatrix=calcCov(twoDim)
    eigVal, eigVec=np.linalg.eig(covMatrix)
    eigPair=sortEig(eigVal,eigVec)
    transform=makeTransformMatrix(eigPair,pcA,pcB)
    transform=np.transpose(transform)
    coords=np.dot(transform,twoDim)
    tmp=np.column_stack((coords[0,:],coords[1,:]))
    return tmp

#GPA
def meanShape(monsters):
    return monsters.mean(axis=2)

def alignToOne(monsters):
    i,j,k=monsters.shape
    #scale and center the first monster
    monsters[:,:,0]=centerShape(monsters[:,:,0])
    monsters[:,:,0]=scaleShape(monsters[:,:,0])

    #align other monster to it
    for x in range(1,k):
        monsters[:,:,x]=procrustesAlign(monsters[:,:,0],monsters[:,:,x])

    return monsters
def alignToOneNoScale(monsters):
    i,j,k=monsters.shape
    #scale and center the first monster
    monsters[:,:,0]=centerShape(monsters[:,:,0])

    #align other monster to it
    for x in range(1,k):
        monsters[:,:,x]=procrustesAlignNoScale(monsters[:,:,0],monsters[:,:,x])
    return monsters

def alignToMean(monsters,trys):
    mean1=meanShape(monsters)
    i,j,k=monsters.shape
    for x in range(k):
        monsters[:,:,x]=procrustesAlign(mean1,monsters[:,:,x])
    mean2=meanShape(monsters)

    diff=np.linalg.norm(mean1-mean2)
    trys=trys+1
    if diff>.00001 and trys< 50:
        alignToMean(monsters, trys)

    return monsters

def alignToMeanNoScale(monsters,trys):
    mean1=meanShape(monsters)
    i,j,k=monsters.shape
    for x in range(k):
        monsters[:,:,x]=procrustesAlignNoScale(mean1,monsters[:,:,x])
    mean2=meanShape(monsters)

    diff=np.linalg.norm(mean1-mean2)
    trys=trys+1
    if diff>.00001 and trys< 50:
        alignToMeanNoScale(monsters, trys)

    return monsters

def centSize(monsters):
    i,j,k=monsters.shape
    size=np.zeros(k)
    for x in range(k):
        size[x]=np.linalg.norm(monsters[:,:,x])
    return size

def doGPA(monsters):
    monsters=alignToOne(monsters)
    mShape=meanShape(monsters)
    monsters=alignToMean(monsters,1)
    mShape=meanShape(monsters)
    return monsters, mShape

def doGPANoScale(monsters):
    monsters=alignToOneNoScale(monsters)
    mShape=meanShape(monsters)
    monsters=alignToMeanNoScale(monsters,1)
    mShape=meanShape(monsters)
    return monsters, mShape

def procDist(monsters,mshape):
    i,j,k=monsters.shape

    procDists=np.zeros(k)
    for x in range(k):
        tmp=monsters[:,:,x]-mshape
        procDists[x]=np.linalg.norm(tmp,'fro')

    return procDists

def procDistPP(monsters):
    i,j,k=monsters.shape

    procDists=np.zeros((k,k))
    for x in range(k):
        for y in range(k):
            tmp=monsters[:,:,x]-monsters[:,:,y]
            procDists[x,y]=np.linalg.norm(tmp,'fro')

    return procDists

def makeTwoDim2(monsters):
    i,j,k=monsters.shape
    tmp=np.zeros((i*j,k))
    for x in range(k):
        vec=np.reshape(monsters[:,:,x],(i*j),order='C')
        tmp[:,x]=vec
    return tmp

max_stars_repo_name: jamesobutler/SlicerMorph

max_issues_repo_name: jamesobutler/SlicerMorph

max_forks_repo_name: jamesobutler/SlicerMorph


text:
"""Transform a roidb into a trainable roidb by adding a bunch of metadata."""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import pickle

from ..datasets.imdb import imdb
import numpy as np
from ..utils.config import cfg
from ..datasets.factory import get_imdb
import PIL
import pdb

def prepare_roidb(imdb_obj):
  """Enrich the imdb's roidb by adding some derived quantities that
  are useful for training. This function precomputes the maximum
  overlap, taken over ground-truth boxes, between each ROI and
  each ground-truth box. The class with maximum overlap is also
  recorded.
  """

  roidb = imdb_obj.roidb
  if not (imdb_obj.name.startswith('coco')):
    cache_file = os.path.join(imdb_obj.cache_path, imdb_obj.name + '_sizes.pkl')
    if os.path.exists(cache_file):
      print('Image sizes loaded from %s' % cache_file)
      with open(cache_file, 'rb') as f:
        sizes = pickle.load(f)
    else:
      print('Extracting image sizes... (It may take long time)')
      sizes = [PIL.Image.open(imdb_obj.image_path_at(i)).size
               for i in range(imdb_obj.num_images)]
      with open(cache_file, 'wb') as f:
        pickle.dump(sizes, f)
      print('Done!!')

  for i in range(len(imdb_obj.image_index)):
    roidb[i]['img_id'] = imdb_obj.image_id_at(i)
    roidb[i]['image'] = imdb_obj.image_path_at(i)
    if not (imdb_obj.name.startswith('coco')):
      roidb[i]['width'] = sizes[i][0]
      roidb[i]['height'] = sizes[i][1]
    # need gt_overlaps as a dense array for argmax
    gt_overlaps = roidb[i]['gt_overlaps'].toarray()
    # max overlap with gt over classes (columns)
    max_overlaps = gt_overlaps.max(axis=1)
    # gt class that had the max overlap
    max_classes = gt_overlaps.argmax(axis=1)
    roidb[i]['max_classes'] = max_classes
    roidb[i]['max_overlaps'] = max_overlaps
    # sanity checks
    # max overlap of 0 => class should be zero (background)
    zero_inds = np.where(max_overlaps == 0)[0]
    assert all(max_classes[zero_inds] == 0)
    # max overlap > 0 => class should not be zero (must be a fg class)
    nonzero_inds = np.where(max_overlaps > 0)[0]
    assert all(max_classes[nonzero_inds] != 0)

def rank_roidb_ratio(roidb):
    # rank roidb based on the ratio between width and height.
  ratio_large = 2  # largest ratio to preserve.
  ratio_small = 0.5  # smallest ratio to preserve.

  ratio_list = []
  for i in range(len(roidb)):
    width = roidb[i]['width']
    height = roidb[i]['height']
    ratio = width / float(height)

    if ratio > ratio_large:
      roidb[i]['need_crop'] = 1
      ratio = ratio_large
    elif ratio < ratio_small:
      roidb[i]['need_crop'] = 1
      ratio = ratio_small
    else:
      roidb[i]['need_crop'] = 0

    ratio_list.append(ratio)

  ratio_list = np.array(ratio_list)
  ratio_index = np.argsort(ratio_list)
  return ratio_list[ratio_index], ratio_index

def filter_roidb(roidb):
  # filter the image without bounding box.
  print('before filtering, there are %d images...' % (len(roidb)))
  i = 0
  while i < len(roidb):
    if len(roidb[i]['boxes']) == 0:
      del roidb[i]
      i -= 1
    i += 1

  print('after filtering, there are %d images...' % (len(roidb)))
  return roidb

def combined_roidb(imdb_names, training=True):
  """
  Combine multiple roidbs
  """

  def get_training_roidb(imdb_obj):
    """Returns a roidb (Region of Interest database) for use in training."""
    if cfg.TRAIN.USE_FLIPPED:
      print('Appending horizontally-flipped training examples...')
      imdb_obj.append_flipped_images()
      print('done')

    print('Preparing training data...')

    prepare_roidb(imdb_obj)
    # ratio_index = rank_roidb_ratio(imdb)
    print('done')

    return imdb_obj.roidb

  def get_roidb(imdb_name):
    imdb_obj = get_imdb(imdb_name)
    print('Loaded dataset `{:s}`'.format(imdb_obj.name))
    imdb_obj.set_proposal_method(cfg.TRAIN.PROPOSAL_METHOD)
    print('Set proposal method: {:s}'.format(cfg.TRAIN.PROPOSAL_METHOD))
    roidb = get_training_roidb(imdb_obj)
    return roidb

  roidbs = [get_roidb(s) for s in imdb_names.split('+')]
  roidb = roidbs[0]

  if len(roidbs) > 1:
    for r in roidbs[1:]:
      roidb.extend(r)
    tmp = get_imdb(imdb_names.split('+')[1])
    imdb_obj = imdb(imdb_names, tmp.classes)
  else:
    imdb_obj = get_imdb(imdb_names)

  if training:
    roidb = filter_roidb(roidb)

  ratio_list, ratio_index = rank_roidb_ratio(roidb)

  return imdb_obj, roidb, ratio_list, ratio_index

max_stars_repo_name: wangvation/torch-mobilenet

max_issues_repo_name: wangvation/torch-mobilenet

max_forks_repo_name: wangvation/torch-mobilenet


text:
import argparse
import datetime
import gym
import numpy as np

import os
import logging
import sys
import socket
import torch

from ppo_model import Memory, PPO, ActorCritic

'''
parser = argparse.ArgumentParser(description='PyTorch REINFORCE example')
parser.add_argument('--env_name', default="HalfCheetah-v2",
                    help='name of the environment to run')
parser.add_argument('--seed', type=int, default=456, metavar='N',
                    help='random seed (default: 456)')
'''

hostname = socket.gethostname()
if not os.path.exists('./ppo_log/' + hostname):
    os.makedirs('./ppo_log/' + hostname)
output_file = './ppo_log/' + hostname + '/output.log'
cal_file = './ppo_log/' + hostname + '/cal.log'
logger = logging.getLogger('mylogger')
logger.setLevel(logging.INFO)

file_handler = logging.FileHandler(output_file, mode='a')
file_handler.setLevel(logging.INFO)
file_handler.setFormatter(logging.Formatter("%(asctime)s - %(levelname)s - %(message)s"))
stdout_handler = logging.StreamHandler(sys.stdout)
stdout_handler.setLevel(logging.INFO)
logger.addHandler(file_handler)
logger.addHandler(stdout_handler)
logger_cal = logging.getLogger('loggercal')
logger_cal.setLevel(logging.INFO)
cal_f_handler = logging.FileHandler(cal_file, mode='a')
file_handler.setLevel(logging.INFO)
logger_cal.addHandler(cal_f_handler)

#env = gym.make(args.env_name)
'''
env.seed(args.seed)
torch.manual_seed(args.seed)
np.random.seed(args.seed)
'''

def main():
    ############## Hyperparameters ##############
    env_name = "HalfCheetah-v2"
    render = True#False
    #solved_reward = 300         # stop training if avg_reward > solved_reward
    #log_interval = 20           # print avg reward in the interval
    max_episodes = 10000        # max training episodes
    max_timesteps = 4000        # max timesteps in one episode

    update_timestep = 1000      # update policy every n timesteps
    #action_std = 0.5            # constant std for action distribution (Multivariate Normal)
    K_epochs = 4               # update policy for K epochs
    eps_clip = 0.2              # clip parameter for PPO
    gamma = 0.99                # discount factor

    lr = 0.0003                 # parameters for Adam optimizer
    betas = (0.9, 0.999)
    hidden_dim = 256
    random_seed = None
    #############################################

    # creating environment
    env = gym.make(env_name)
    state_dim = env.observation_space.shape[0]
    action_dim = env.action_space.shape[0]

    print('state_dim',state_dim, '   action_dim',action_dim)
    print('load env succeed!*********************************')

    if random_seed:
        print("Random Seed: {}".format(random_seed))
        torch.manual_seed(random_seed)
        env.seed(random_seed)
        np.random.seed(random_seed)

    memory = Memory()
    ppo = PPO(state_dim, action_dim, hidden_dim, lr, betas, gamma, K_epochs, eps_clip)
    ppo.policy.load_state_dict(torch.load('ppo_model/ppo_500.pth'))
    print(lr,betas)

    # logging variables
    episode_reward = 0
    avg_length = 0
    time_step = 0

    # training loop
    for i_episode in range(1, max_episodes+1):
        state = env.reset()
        for t in range(max_timesteps):
            time_step +=1
            # Running policy_old:
            action, mean = ppo.select_action(state, memory)
            next_state, reward, done, _ = env.step(action)
            # Saving reward:
            '''
            memory.rewards.append(reward)
            # update if its time
            if time_step % update_timestep == 0:
                ppo.update(memory)
                print('update')
                memory.clear_memory()
                time_step = 0
            episode_reward += reward
            '''
            state = next_state
            if render:
                env.render()
            if done:
                break

        #avg_length += t
        '''
        logger.info("Episode: {}, total numsteps: {}, reward: {}".format(i_episode, time_step, round(episode_reward, 2)))
        #logger.info('Env %d, Goal (%05.1f, %05.1f), Episode %05d, setp %03d, Reward %-5.1f,%s' % (env.index, env.goal_point[0], env.goal_point[1], id + 1, step, ep_reward, result))
        logger_cal.info(episode_reward)

        # save every 500 episodes
        if i_episode % 300 == 0:
            torch.save(ppo.policy.state_dict(), './ppo_model/ppo_{}.pth'.format(i_episode))
        '''

if __name__ == '__main__':
    main()


text:
import os
import cflearn

import numpy as np

from cflearn.misc.toolkit import get_latest_workplace

base = cflearn.ml.CarefreePipeline
file_folder = os.path.dirname(__file__)
train_file = os.path.join(file_folder, "train.csv")
test_file = os.path.join(file_folder, "test.csv")
data = cflearn.MLData.with_cf_data(
    train_file,
    cf_data_config={"label_name": "Survived"},
)
m = base().fit(data)
assert isinstance(m, base)

idata = cflearn.MLInferenceData(test_file)
results = m.predict(idata, make_loader_kwargs={"contains_labels": False})
predictions = results[cflearn.PREDICTIONS_KEY]

export_folder = "titanic"
m.save(export_folder)
m2 = base.load(export_folder)
results = m2.predict(idata, make_loader_kwargs={"contains_labels": False})
assert np.allclose(predictions, results[cflearn.PREDICTIONS_KEY])

latest = get_latest_workplace("_logs")
assert latest is not None
m3 = base.load(base.pack(latest))
results = m3.predict(idata, make_loader_kwargs={"contains_labels": False})
assert np.allclose(predictions, results[cflearn.PREDICTIONS_KEY])

onnx_folder = "titanic_onnx"
m.to_onnx(onnx_folder)
m4 = cflearn.ml.CarefreePipeline.from_onnx(onnx_folder)
results = m4.predict(idata, make_loader_kwargs={"contains_labels": False})
assert np.allclose(predictions, results[cflearn.PREDICTIONS_KEY], atol=1.0e-4)


text:
import numpy as np
import os
import ctypes as c

CAMERA_ID_MAX = 128

BAYER_RG = 0
BAYER_BG = 1
BAYER_GR = 2
BAYER_GB = 3

IMG_RAW8 = 0
IMG_RGB24 = 1
IMG_RAW16 = 2
IMG_Y8 = 3
IMG_END = -1

GUIDE_NORTH = 0
GUIDE_SOUTH = 1
GUIDE_EAST = 2
GUIDE_WEST = 3

FLIP_NONE = 0
FLIP_HORIZ = 1
FLIP_VERT = 2
FLIP_BOTH = 3

MODE_NORMAL = 0
MODE_TRIG_SOFT_EDGE = 1
MODE_TRIG_RISE_EDGE = 2
MODE_TRIG_FALL_EDGE = 3
MODE_TRIG_SOFT_LEVE = 4
MODE_TRIG_HIGH_LEVEL = 5
MODE_TRIG_LOW_LEVEL = 6
MODE_END = -1

SUCCESS = 0
ERROR_INVALID_INDEX = 1
ERROR_INVALID_ID = 2
ERROR_INVALID_CONTROL_TYPE = 3
ERROR_CAMERA_CLOSED = 4
ERROR_CAMERA_REMOVED = 5
ERROR_INVALID_PATH = 6
ERROR_INVALID_FILEFORMAT = 7
ERROR_INVALID_SIZE = 8
ERROR_INVALID_IMGTYPE = 9
ERROR_OUTOF_BOUNDARY = 10
ERROR_TIMEOUT = 11
ERROR_INVALID_SEQUENCE = 12
ERROR_BUFFER_TOO_SMALL = 13
ERROR_VIDEO_MODE_ACTIVE = 14
ERROR_EXPOSURE_IN_PROGRESS = 15
ERROR_GENERAL_ERROR = 16
ERROR_INVALID_MODE = 17
ERROR_END = 18

FALSE = 0
TRUE = 1

GAIN = 0
EXPOSURE = 1
GAMMA = 2
WB_R = 3
WB_B = 4
OFFSET = 5
BANDWIDTHOVERLOAD = 6
OVERCLOCK = 7
TEMPERATURE = 8
FLIP = 9
AUTO_MAX_GAIN = 10
AUTO_MAX_EXP = 11
AUTO_TARGET_BRIGHTNESS = 12
HARDWARE_BIN = 13
HIGH_SPEED_MODE = 14
COOLER_POWER_PERC = 15
TARGET_TEMP = 16
COOLER_ON = 17
MONO_BIN = 18
FAN_ON = 19
PATTERN_ADJUST = 20
ANTI_DEW_HEATER = 21

BRIGHTNESS = OFFSET
AUTO_MAX_BRIGHTNESS = AUTO_TARGET_BRIGHTNESS

EXP_IDLE = 0
EXP_WORKING = 1
EXP_SUCCESS = 2
EXP_FAILED = 3

class CAMERA_INFO(c.Structure):
    _fields_ = [
        ("Name", c.c_char * 64),
        ("CameraID", c.c_int),
        ("MaxHeight", c.c_long),
        ("MaxWidth", c.c_long),
        ("IsColorCam", c.c_int),
        ("BayerPattern", c.c_int),
        ("SupportedBins", c.c_int * 16),
        ("SupportedVideoFormat", c.c_int * 8),
        ("PixelSize", c.c_double),
        ("MechanicalShutter", c.c_int),
        ("ST4Port", c.c_int),
        ("IsCoolerCam", c.c_int),
        ("IsUSB3Host", c.c_int),
        ("IsUSB3Camera", c.c_int),
        ("ElecPerADU", c.c_float),
        ("BitDepth", c.c_int),
        ("IsTriggerCam", c.c_int),
        ("Unused", c.c_byte * 16),
    ]

class CONTROL_CAPS(c.Structure):
    _fields_ = [
        ("Name", c.c_char * 64),
        ("Description", c.c_char * 128),
        ("MaxValue", c.c_long),
        ("MinValue", c.c_long),
        ("DefaultValue", c.c_long),
        ("IsAutoSupported", c.c_int),
        ("IsWritable", c.c_int),
        ("ControlType", c.c_int),
        ("Unused", c.c_byte * 32),
    ]

class ID(c.Structure):
    _fields_ = [
        ("id", c.c_char * 8),
    ]

class SUPPORTED_MODE(c.Structure):
    _fields_ = [
        ("SupportedCameraMode", c.c_int * 16),
    ]

_im_Bpp = (1, 3, 2, 1)
_im_ds = (np.uint8, np.uint8, np.uint16, np.uint8)
requires = ('C', 'W', 'A')

my_path = os.path.dirname(os.path.abspath(__file__))
lib = np.ctypeslib.load_library('libASICamera2', my_path)

lib.ASIGetNumOfConnectedCameras.restype = c.c_int
lib.ASIGetNumOfConnectedCameras.argtypes = []

lib.ASIGetProductIDs.restype = c.c_int
lib.ASIGetProductIDs.argtypes = [c.POINTER(c.c_int)]

lib.ASIGetCameraProperty.restype = c.c_int
lib.ASIGetCameraProperty.argtypes = [c.POINTER(CAMERA_INFO), c.c_int]

lib.ASIOpenCamera.restype = c.c_int
lib.ASIOpenCamera.argtypes = [c.c_int]

lib.ASIInitCamera.restype = c.c_int
lib.ASIInitCamera.argtypes = [c.c_int]

lib.ASICloseCamera.restype = c.c_int
lib.ASICloseCamera.argtypes = [c.c_int]

lib.ASIGetNumOfControls.restype = c.c_int
lib.ASIGetNumOfControls.argtypes = [c.c_int, c.POINTER(c.c_int)]

lib.ASIGetControlCaps.restype = c.c_int
lib.ASIGetControlCaps.argtypes = [c.c_int, c.c_int, c.POINTER(CONTROL_CAPS)]

lib.ASIGetControlValue.restype = c.c_int
lib.ASIGetControlValue.argtypes = [c.c_int, c.c_int, c.POINTER(c.c_long),
                                   c.POINTER(c.c_int)]

lib.ASISetControlValue.restype = c.c_int
lib.ASISetControlValue.argtypes = [c.c_int, c.c_int, c.c_long, c.c_int]

lib.ASISetROIFormat.restype = c.c_int
lib.ASISetROIFormat.argtypes = [c.c_int, c.c_int, c.c_int, c.c_int, c.c_int]

lib.ASIGetROIFormat.restype = c.c_int
lib.ASIGetROIFormat.argtypes = [c.c_int, c.POINTER(c.c_int),
                                c.POINTER(c.c_int), c.POINTER(c.c_int),
                                c.POINTER(c.c_int)]

lib.ASISetStartPos.restype = c.c_int
lib.ASISetStartPos.argtypes = [c.c_int, c.c_int, c.c_int]

lib.ASIGetStartPos.restype = c.c_int
lib.ASIGetStartPos.argtypes = [c.c_int, c.POINTER(c.c_int), c.POINTER(c.c_int)]

lib.ASIGetDroppedFrames.restype = c.c_int
lib.ASIGetDroppedFrames.argtypes = [c.c_int, c.POINTER(c.c_int)]

lib.ASIEnableDarkSubtract.restype = c.c_int
lib.ASIEnableDarkSubtract.argtypes = [c.c_int, c.c_char_p]

lib.ASIDisableDarkSubtract.restype = c.c_int
lib.ASIDisableDarkSubtract.argtypes = [c.c_int]

lib.ASIStartVideoCapture.restype = c.c_int
lib.ASIStartVideoCapture.argtypes = [c.c_int]

lib.ASIStopVideoCapture.restype = c.c_int
lib.ASIStopVideoCapture.argtypes = [c.c_int]

lib.ASIGetVideoData.restype = c.c_int
lib.ASIGetVideoData.argtypes = [
    c.c_int, np.ctypeslib.ndpointer(dtype=np.uint8, ndim=1,
                                    flags=requires),
    c.c_long, c.c_int]

lib.ASIPulseGuideOn.restype = c.c_int
lib.ASIPulseGuideOn.argtypes = [c.c_int, c.c_int]

lib.ASIPulseGuideOff.restype = c.c_int
lib.ASIPulseGuideOff.argtypes = [c.c_int, c.c_int]

lib.ASIStartExposure.restype = c.c_int
lib.ASIStartExposure.argtypes = [c.c_int, c.c_int]

lib.ASIStopExposure.restype = c.c_int
lib.ASIStopExposure.argtypes = [c.c_int]

lib.ASIGetExpStatus.restype = c.c_int
lib.ASIGetExpStatus.argtypes = [c.c_int, c.POINTER(c.c_int)]

lib.ASIGetDataAfterExp.restype = c.c_int
lib.ASIGetDataAfterExp.argtypes = [
    c.c_int, np.ctypeslib.ndpointer(dtype=np.uint8, ndim=1,
                                    flags=requires),
    c.c_int]

lib.ASIGetID.restype = c.c_int
lib.ASIGetID.argtypes = [c.c_int, c.POINTER(c.c_int)]

lib.ASISetID.restype = c.c_int
lib.ASISetID.argtypes = [c.c_int, c.c_int]

lib.ASIGetGainOffset.restype = c.c_int
lib.ASIGetGainOffset.argtypes = [c.c_int,
                                 c.POINTER(c.c_int), c.POINTER(c.c_int),
                                 c.POINTER(c.c_int), c.POINTER(c.c_int)]

lib.ASIGetSDKVersion.restype = c.c_char_p
lib.ASIGetSDKVersion.argtypes = []

def GetNumOfConnectedCameras():
    return lib.ASIGetNumOfConnectedCameras()

def GetProductIDs():
    pids = (c.c_int * CAMERA_ID_MAX)()
    n = lib.ASIGetProductIDs(pids)
    return pids[:n]

def GetSDKVersion():
    return lib.ASIGetSDKVersion()

class Error(Exception):
    def __init__(self, err):
        super().__init__("ASI Camera error: %d" % err)
        self.code = err

class Camera:

    def __init__(self, i):
        self.i = i
        self.width = -1
        self.height = -1
        self.img_type = -1

    def GetCameraProperty(self):
        info = CAMERA_INFO()
        err = lib.ASIGetCameraProperty(c.byref(info), self.i)
        if err != 0:
            raise Error(err)
        return info

    def OpenCamera(self):
        err = lib.ASIOpenCamera(self.i)
        if err != 0:
            raise Error(err)

    def InitCamera(self):
        err = lib.ASIInitCamera(self.i)
        if err != 0:
            raise Error(err)

    def CloseCamera(self):
        err = lib.ASICloseCamera(self.i)
        if err != 0:
            raise Error(err)

    def GetNumOfControls(self):
        n = c.c_int()
        err = lib.ASIGetNumOfControls(self.i, c.byref(n))
        if err != 0:
            raise Error(err)
        return n.value

    def GetControlCaps(self, idx):
        caps = CONTROL_CAPS()
        err = lib.ASIGetControlCaps(self.i, idx, c.byref(caps))
        if err != 0:
            raise Error(err)
        return caps

    def GetControlValue(self, idx):
        val = c.c_long()
        auto = c.c_int()
        err = lib.ASIGetControlValue(self.i, idx, c.byref(val), c.byref(auto))
        if err != 0:
            raise Error(err)
        return (val.value, bool(auto.value))

    def SetControlValue(self, idx, val, auto):
        err = lib.ASISetControlValue(self.i, idx, c.c_long(val), c.c_int(auto))
        if err != 0:
            raise Error(err)

    def SetROIFormat(self, width, height, binning, img_type):
        err = lib.ASISetROIFormat(self.i, width, height, binning, img_type)
        if err != 0:
            raise Error(err)
        self.width = width
        self.height = height
        self.img_type = img_type

    def GetROIFormat(self):
        width = c.c_int()
        height = c.c_int()
        binning = c.c_int()
        img_type = c.c_int()
        err = lib.ASIGetROIFormat(self.i, c.byref(width), c.byref(height),
                                  c.byref(binning), c.byref(img_type))
        if err != 0:
            raise Error(err)
        self.width = width.value
        self.height = height.value
        self.img_type = img_type.value
        return (width.value, height.value, binning.value, img_type.value)

    def SetStartPos(self, x, y):
        err = lib.ASISetStartPos(self.i, x, y)
        if err != 0:
            raise Error(err)

    def GetStartPos(self):
        x = c.c_int()
        y = c.c_int()
        err = lib.ASIGetStartPos(self.i, c.byref(x), c.byref(y))
        if err != 0:
            raise Error(err)
        return x.value, y.value

    def GetDroppedFrames(self):
        dropped = c.c_int()
        err = lib.ASIGetDroppedFrames(self.i, c.byref(dropped))
        if err != 0:
            raise Error(err)
        return dropped.value

    def EnableDarkSubtract(self, fname):
        err = lib.ASIEnableDarkSubtract(self.i, fname)
        if err != 0:
            raise Error(err)

    def DisableDarkSubtract(self):
        err = lib.ASIDisableDarkSubtract(self.i)
        if err != 0:
            raise Error(err)

    def StartVideoCapture(self):
        err = lib.ASIStartVideoCapture(self.i)
        if err != 0:
            raise Error(err)

    def StopVideoCapture(self):
        err = lib.ASIStopVideoCapture(self.i)
        if err != 0:
            raise Error(err)

    def convert(self, x):
        if self.img_type == -1:
            self.GetROIFormat()
        if self.img_type == 1:
            imr = x.reshape(self.height, self.width, 3)
            return imr
        else:
            return x.view(
                dtype=_im_ds[self.img_type]).reshape(
                    self.height, self.width)

    def bufAlloc(self):
        if self.img_type == -1:
            self.GetROIFormat()
        t = self.width * self.height * _im_Bpp[self.img_type]
        return (np.require(np.zeros(t), np.uint8, requires), t)

    def GetVideoData(self, wait):
        (b, l) = self.bufAlloc()
        err = lib.ASIGetVideoData(self.i, b, l, wait)
        if err != 0:
            raise Error(err)
        return self.convert(b)

    def PulseGuideOn(self, dire):
        err = lib.ASIPulseGuideOn(self.i, dire)
        if err != 0:
            raise Error(err)

    def PulseGuideOff(self, dire):
        err = lib.ASIPulseGuideOff(self.i, dire)
        if err != 0:
            raise Error(err)

    def StartExposure(self, dark):
        err = lib.ASIStartExposure(self.i, c.c_int(dark))
        if err != 0:
            raise Error(err)

    def StopExposure(self):
        err = lib.ASIStopExposure(self.i)
        if err != 0:
            raise Error(err)

    def GetExpStatus(self):
        st = c.c_int()
        err = lib.ASIGetExpStatus(self.i, c.byref(st))
        if err != 0:
            raise Error(err)
        return st.value

    def GetDataAfterExp(self):
        (b, l) = self.bufAlloc()
        err = lib.ASIGetDataAfterExp(self.i, b, l)
        if err != 0:
            raise Error(err)
        return self.convert(b)

    def GetID(self):
        id = c.c_int()
        err = lib.ASIGetID(self.i, c.byref(id))
        if err != 0:
            raise Error(err)
        return id.value

    def SetID(self, id):
        err = lib.ASISetID(self.i, id)
        if err != 0:
            raise Error(err)

    def GetGainOffset(self):
        d1 = c.c_int()
        d2 = c.c_int()
        d3 = c.c_int()
        d4 = c.c_int()
        err = lib.ASIGetGainOffset(self.i, c.byref(d1), c.byref(d2),
                                   c.byref(d3), c.byref(d4))
        if err != 0:
            raise Error(err)
        return d1.value, d2.value, d3.value, d4.value

max_stars_repo_name: chripell/fit-image-helper

max_issues_repo_name: chripell/fit-image-helper

max_forks_repo_name: chripell/fit-image-helper