text:
import os

import pytest
import numpy as np
from urllib.error import HTTPError, URLError

from astropy.time import Time
from astropy import units as u
from astropy.constants import c
from astropy.coordinates.builtin_frames import GCRS, TETE
from astropy.coordinates.earth import EarthLocation
from astropy.coordinates.sky_coordinate import SkyCoord
from astropy.coordinates.representation import CartesianRepresentation, UnitSphericalRepresentation
from astropy.coordinates.solar_system import (get_body, get_moon, BODY_NAME_TO_KERNEL_SPEC,
                                              _get_apparent_body_position, solar_system_ephemeris,
                                              get_body_barycentric, get_body_barycentric_posvel)
from astropy.coordinates.funcs import get_sun
from astropy.tests.helper import assert_quantity_allclose
from astropy.units import allclose as quantity_allclose
from astropy.utils.data import download_file, get_pkg_data_filename
from astropy.utils.compat.optional_deps import (HAS_JPLEPHEM,  # noqa
                                                HAS_SKYFIELD)

if HAS_SKYFIELD:
    from skyfield.api import Loader, Topos

de432s_separation_tolerance_planets = 5*u.arcsec
de432s_separation_tolerance_moon = 5*u.arcsec
de432s_distance_tolerance = 20*u.km

skyfield_angular_separation_tolerance = 1*u.arcsec
skyfield_separation_tolerance = 10*u.km

@pytest.mark.remote_data
@pytest.mark.skipif('not HAS_SKYFIELD')
def test_positions_skyfield(tmpdir):
    """
    Test positions against those generated by skyfield.
    """
    load = Loader(tmpdir)
    t = Time('1980-03-25 00:00')
    location = None

    # skyfield ephemeris
    try:
        planets = load('de421.bsp')
        ts = load.timescale()
    except OSError as e:
        if os.environ.get('CI', False) and 'timed out' in str(e):
            pytest.xfail('Timed out in CI')
        else:
            raise

    mercury, jupiter, moon = planets['mercury'], planets['jupiter barycenter'], planets['moon']
    earth = planets['earth']

    skyfield_t = ts.from_astropy(t)

    if location is not None:
        earth = earth+Topos(latitude_degrees=location.lat.to_value(u.deg),
                            longitude_degrees=location.lon.to_value(u.deg),
                            elevation_m=location.height.to_value(u.m))

    skyfield_mercury = earth.at(skyfield_t).observe(mercury).apparent()
    skyfield_jupiter = earth.at(skyfield_t).observe(jupiter).apparent()
    skyfield_moon = earth.at(skyfield_t).observe(moon).apparent()

    if location is not None:
        frame = TETE(obstime=t, location=location)
    else:
        frame = TETE(obstime=t)

    ra, dec, dist = skyfield_mercury.radec(epoch='date')
    skyfield_mercury = SkyCoord(ra.to(u.deg), dec.to(u.deg), distance=dist.to(u.km),
                                frame=frame)
    ra, dec, dist = skyfield_jupiter.radec(epoch='date')
    skyfield_jupiter = SkyCoord(ra.to(u.deg), dec.to(u.deg), distance=dist.to(u.km),
                                frame=frame)
    ra, dec, dist = skyfield_moon.radec(epoch='date')
    skyfield_moon = SkyCoord(ra.to(u.deg), dec.to(u.deg), distance=dist.to(u.km),
                             frame=frame)

    # planet positions w.r.t true equator and equinox
    moon_astropy = get_moon(t, location, ephemeris='de430').transform_to(frame)
    mercury_astropy = get_body('mercury', t, location, ephemeris='de430').transform_to(frame)
    jupiter_astropy = get_body('jupiter', t, location, ephemeris='de430').transform_to(frame)

    assert (moon_astropy.separation(skyfield_moon) <
            skyfield_angular_separation_tolerance)
    assert (moon_astropy.separation_3d(skyfield_moon) < skyfield_separation_tolerance)

    assert (jupiter_astropy.separation(skyfield_jupiter) <
            skyfield_angular_separation_tolerance)
    assert (jupiter_astropy.separation_3d(skyfield_jupiter) <
            skyfield_separation_tolerance)

    assert (mercury_astropy.separation(skyfield_mercury) <
            skyfield_angular_separation_tolerance)
    assert (mercury_astropy.separation_3d(skyfield_mercury) <
            skyfield_separation_tolerance)
    planets.close()

class TestPositionsGeocentric:
    """
    Test positions against those generated by JPL Horizons accessed on
    2016-03-28, with refraction turned on.
    """

    def setup(self):
        self.t = Time('1980-03-25 00:00')
        self.apparent_frame = TETE(obstime=self.t)
        # Results returned by JPL Horizons web interface
        self.horizons = {
            'mercury': SkyCoord(ra='22h41m47.78s', dec='-08d29m32.0s',
                                distance=c*6.323037*u.min, frame=self.apparent_frame),
            'moon': SkyCoord(ra='07h32m02.62s', dec='+18d34m05.0s',
                             distance=c*0.021921*u.min, frame=self.apparent_frame),
            'jupiter': SkyCoord(ra='10h17m12.82s', dec='+12d02m57.0s',
                                distance=c*37.694557*u.min, frame=self.apparent_frame),
            'sun': SkyCoord(ra='00h16m31.00s', dec='+01d47m16.9s',
                            distance=c*8.294858*u.min, frame=self.apparent_frame)}

    @pytest.mark.parametrize(('body', 'sep_tol', 'dist_tol'),
                             (('mercury', 7.*u.arcsec, 1000*u.km),
                              ('jupiter', 78.*u.arcsec, 76000*u.km),
                              ('moon', 20.*u.arcsec, 80*u.km),
                              ('sun', 5.*u.arcsec, 11.*u.km)))
    def test_erfa_planet(self, body, sep_tol, dist_tol):
        """Test predictions using erfa/plan94.

        Accuracies are maximum deviations listed in erfa/plan94.c, for Jupiter and
        Mercury, and that quoted in Meeus "Astronomical Algorithms" (1998) for the Moon.
        """
        astropy = get_body(body, self.t, ephemeris='builtin')
        horizons = self.horizons[body]

        # convert to true equator and equinox
        astropy = astropy.transform_to(self.apparent_frame)

        # Assert sky coordinates are close.
        assert astropy.separation(horizons) < sep_tol

        # Assert distances are close.
        assert_quantity_allclose(astropy.distance, horizons.distance,
                                 atol=dist_tol)

    @pytest.mark.remote_data
    @pytest.mark.skipif('not HAS_JPLEPHEM')
    @pytest.mark.parametrize('body', ('mercury', 'jupiter', 'sun'))
    def test_de432s_planet(self, body):
        astropy = get_body(body, self.t, ephemeris='de432s')
        horizons = self.horizons[body]

        # convert to true equator and equinox
        astropy = astropy.transform_to(self.apparent_frame)

        # Assert sky coordinates are close.
        assert (astropy.separation(horizons) <
                de432s_separation_tolerance_planets)

        # Assert distances are close.
        assert_quantity_allclose(astropy.distance, horizons.distance,
                                 atol=de432s_distance_tolerance)

    @pytest.mark.remote_data
    @pytest.mark.skipif('not HAS_JPLEPHEM')
    def test_de432s_moon(self):
        astropy = get_moon(self.t, ephemeris='de432s')
        horizons = self.horizons['moon']

        # convert to true equator and equinox
        astropy = astropy.transform_to(self.apparent_frame)

        # Assert sky coordinates are close.
        assert (astropy.separation(horizons) <
                de432s_separation_tolerance_moon)

        # Assert distances are close.
        assert_quantity_allclose(astropy.distance, horizons.distance,
                                 atol=de432s_distance_tolerance)

class TestPositionKittPeak:
    """
    Test positions against those generated by JPL Horizons accessed on
    2016-03-28, with refraction turned on.
    """

    def setup(self):
        kitt_peak = EarthLocation.from_geodetic(lon=-111.6*u.deg,
                                                lat=31.963333333333342*u.deg,
                                                height=2120*u.m)
        self.t = Time('2014-09-25T00:00', location=kitt_peak)
        self.apparent_frame = TETE(obstime=self.t, location=kitt_peak)
        # Results returned by JPL Horizons web interface
        self.horizons = {
            'mercury': SkyCoord(ra='13h38m58.50s', dec='-13d34m42.6s',
                                distance=c*7.699020*u.min, frame=self.apparent_frame),
            'moon': SkyCoord(ra='12h33m12.85s', dec='-05d17m54.4s',
                             distance=c*0.022054*u.min, frame=self.apparent_frame),
            'jupiter': SkyCoord(ra='09h09m55.55s', dec='+16d51m57.8s',
                                distance=c*49.244937*u.min, frame=self.apparent_frame)}

    @pytest.mark.parametrize(('body', 'sep_tol', 'dist_tol'),
                             (('mercury', 7.*u.arcsec, 500*u.km),
                              ('jupiter', 78.*u.arcsec, 82000*u.km)))
    def test_erfa_planet(self, body, sep_tol, dist_tol):
        """Test predictions using erfa/plan94.

        Accuracies are maximum deviations listed in erfa/plan94.c.
        """
        # Add uncertainty in position of Earth
        dist_tol = dist_tol + 1300 * u.km

        astropy = get_body(body, self.t, ephemeris='builtin')
        horizons = self.horizons[body]

        # convert to true equator and equinox
        astropy = astropy.transform_to(self.apparent_frame)

        # Assert sky coordinates are close.
        assert astropy.separation(horizons) < sep_tol

        # Assert distances are close.
        assert_quantity_allclose(astropy.distance, horizons.distance,
                                 atol=dist_tol)

    @pytest.mark.remote_data
    @pytest.mark.skipif('not HAS_JPLEPHEM')
    @pytest.mark.parametrize('body', ('mercury', 'jupiter'))
    def test_de432s_planet(self, body):
        astropy = get_body(body, self.t, ephemeris='de432s')
        horizons = self.horizons[body]

        # convert to true equator and equinox
        astropy = astropy.transform_to(self.apparent_frame)

        # Assert sky coordinates are close.
        assert (astropy.separation(horizons) <
                de432s_separation_tolerance_planets)

        # Assert distances are close.
        assert_quantity_allclose(astropy.distance, horizons.distance,
                                 atol=de432s_distance_tolerance)

    @pytest.mark.remote_data
    @pytest.mark.skipif('not HAS_JPLEPHEM')
    def test_de432s_moon(self):
        astropy = get_moon(self.t, ephemeris='de432s')
        horizons = self.horizons['moon']

        # convert to true equator and equinox
        astropy = astropy.transform_to(self.apparent_frame)

        # Assert sky coordinates are close.
        assert (astropy.separation(horizons) <
                de432s_separation_tolerance_moon)

        # Assert distances are close.
        assert_quantity_allclose(astropy.distance, horizons.distance,
                                 atol=de432s_distance_tolerance)

    @pytest.mark.remote_data
    @pytest.mark.skipif('not HAS_JPLEPHEM')
    @pytest.mark.parametrize('bodyname', ('mercury', 'jupiter'))
    def test_custom_kernel_spec_body(self, bodyname):
        """
        Checks that giving a kernel specifier instead of a body name works
        """
        coord_by_name = get_body(bodyname, self.t, ephemeris='de432s')
        kspec = BODY_NAME_TO_KERNEL_SPEC[bodyname]
        coord_by_kspec = get_body(kspec, self.t, ephemeris='de432s')

        assert_quantity_allclose(coord_by_name.ra, coord_by_kspec.ra)
        assert_quantity_allclose(coord_by_name.dec, coord_by_kspec.dec)
        assert_quantity_allclose(coord_by_name.distance, coord_by_kspec.distance)

@pytest.mark.remote_data
@pytest.mark.skipif('not HAS_JPLEPHEM')
def test_horizons_consistency_with_precision():
    """
    A test to compare at high precision against output of JPL horizons.

    Tests ephemerides, and conversions from ICRS to GCRS to TETE. We are aiming for
    better than 2 milli-arcsecond precision.

    We use the Moon since it is nearby, and moves fast in the sky so we are
    testing for parallax, proper handling of light deflection and aberration.
    """
    # JPL Horizon values for 2020_04_06 00:00 to 23:00 in 1 hour steps
    # JPL Horizons has a known offset (frame bias) of 51.02 mas in RA. We correct that here
    ra_apparent_horizons = [
        170.167332531, 170.560688674, 170.923834838, 171.271663481, 171.620188972, 171.985340827,
        172.381766539, 172.821772139, 173.314502650, 173.865422398, 174.476108551, 175.144332386,
        175.864375310, 176.627519827, 177.422655853, 178.236955730, 179.056584831, 179.867427392,
        180.655815385, 181.409252074, 182.117113814, 182.771311578, 183.366872837, 183.902395443
    ] * u.deg + 51.02376467 * u.mas
    dec_apparent_horizons = [
        10.269112037, 10.058820647, 9.837152044, 9.603724551, 9.358956528, 9.104012390, 8.840674927,
        8.571162442, 8.297917326, 8.023394488, 7.749873882, 7.479312991, 7.213246666, 6.952732614,
        6.698336823, 6.450150213, 6.207828142, 5.970645962, 5.737565957, 5.507313851, 5.278462034,
        5.049521497, 4.819038911, 4.585696512
    ] * u.deg
    with solar_system_ephemeris.set('de430'):
        loc = EarthLocation.from_geodetic(-67.787260*u.deg, -22.959748*u.deg, 5186*u.m)
        times = Time('2020-04-06 00:00') + np.arange(0, 24, 1)*u.hour
        astropy = get_body('moon', times, loc)

        apparent_frame = TETE(obstime=times, location=loc)
        astropy = astropy.transform_to(apparent_frame)
        usrepr = UnitSphericalRepresentation(ra_apparent_horizons, dec_apparent_horizons)
        horizons = apparent_frame.realize_frame(usrepr)
    assert_quantity_allclose(astropy.separation(horizons), 0*u.mas, atol=1.5*u.mas)

@pytest.mark.remote_data
@pytest.mark.skipif('not HAS_JPLEPHEM')
@pytest.mark.parametrize('time', (Time('1960-01-12 00:00'),
                                  Time('1980-03-25 00:00'),
                                  Time('2010-10-13 00:00')))
def test_get_sun_consistency(time):
    """
    Test that the sun from JPL and the builtin get_sun match
    """
    sun_jpl_gcrs = get_body('sun', time, ephemeris='de432s')
    builtin_get_sun = get_sun(time)
    sep = builtin_get_sun.separation(sun_jpl_gcrs)
    assert sep < 0.1*u.arcsec

def test_get_moon_nonscalar_regression():
    """
    Test that the builtin ephemeris works with non-scalar times.

    See Issue #5069.
    """
    times = Time(["2015-08-28 03:30", "2015-09-05 10:30"])
    # the following line will raise an Exception if the bug recurs.
    get_moon(times, ephemeris='builtin')

def test_barycentric_pos_posvel_same():
    # Check that the two routines give identical results.
    ep1 = get_body_barycentric('earth', Time('2016-03-20T12:30:00'))
    ep2, _ = get_body_barycentric_posvel('earth', Time('2016-03-20T12:30:00'))
    assert np.all(ep1.xyz == ep2.xyz)

def test_earth_barycentric_velocity_rough():
    # Check that a time near the equinox gives roughly the right result.
    ep, ev = get_body_barycentric_posvel('earth', Time('2016-03-20T12:30:00'))
    assert_quantity_allclose(ep.xyz, [-1., 0., 0.]*u.AU, atol=0.01*u.AU)
    expected = u.Quantity([0.*u.one,
                           np.cos(23.5*u.deg),
                           np.sin(23.5*u.deg)]) * -30. * u.km / u.s
    assert_quantity_allclose(ev.xyz, expected, atol=1.*u.km/u.s)

def test_earth_barycentric_velocity_multi_d():
    # Might as well test it with a multidimensional array too.
    t = Time('2016-03-20T12:30:00') + np.arange(8.).reshape(2, 2, 2) * u.yr / 2.
    ep, ev = get_body_barycentric_posvel('earth', t)
    # note: assert_quantity_allclose doesn't like the shape mismatch.
    # this is a problem with np.testing.assert_allclose.
    assert quantity_allclose(ep.get_xyz(xyz_axis=-1),
                             [[-1., 0., 0.], [+1., 0., 0.]]*u.AU,
                             atol=0.06*u.AU)
    expected = u.Quantity([0.*u.one,
                           np.cos(23.5*u.deg),
                           np.sin(23.5*u.deg)]) * ([[-30.], [30.]] * u.km / u.s)
    assert quantity_allclose(ev.get_xyz(xyz_axis=-1), expected,
                             atol=2.*u.km/u.s)

@pytest.mark.remote_data
@pytest.mark.skipif('not HAS_JPLEPHEM')
@pytest.mark.parametrize(('body', 'pos_tol', 'vel_tol'),
                         (('mercury', 1000.*u.km, 1.*u.km/u.s),
                          ('jupiter', 100000.*u.km, 2.*u.km/u.s),
                          ('earth', 10*u.km, 10*u.mm/u.s),
                          ('moon', 18*u.km, 50*u.mm/u.s)))
def test_barycentric_velocity_consistency(body, pos_tol, vel_tol):
    # Tolerances are about 1.5 times the rms listed for plan94 and epv00,
    # except for Mercury (which nominally is 334 km rms), and the Moon
    # (which nominally is 6 km rms).
    t = Time('2016-03-20T12:30:00')
    ep, ev = get_body_barycentric_posvel(body, t, ephemeris='builtin')
    dp, dv = get_body_barycentric_posvel(body, t, ephemeris='de432s')
    assert_quantity_allclose(ep.xyz, dp.xyz, atol=pos_tol)
    assert_quantity_allclose(ev.xyz, dv.xyz, atol=vel_tol)
    # Might as well test it with a multidimensional array too.
    t = Time('2016-03-20T12:30:00') + np.arange(8.).reshape(2, 2, 2) * u.yr / 2.
    ep, ev = get_body_barycentric_posvel(body, t, ephemeris='builtin')
    dp, dv = get_body_barycentric_posvel(body, t, ephemeris='de432s')
    assert_quantity_allclose(ep.xyz, dp.xyz, atol=pos_tol)
    assert_quantity_allclose(ev.xyz, dv.xyz, atol=vel_tol)

@pytest.mark.remote_data
@pytest.mark.skipif('not HAS_JPLEPHEM')
@pytest.mark.parametrize('time', (Time('1960-01-12 00:00'),
                                  Time('1980-03-25 00:00'),
                                  Time('2010-10-13 00:00')))
def test_url_or_file_ephemeris(time):
    # URL for ephemeris de432s used for testing:
    url = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de432s.bsp'

    # Pass the ephemeris directly as a URL.
    coord_by_url = get_body('earth', time, ephemeris=url)

    # Translate the URL to the cached location on the filesystem.
    # Since we just used the url above, it should already have been downloaded.
    filepath = download_file(url, cache=True)

    # Get the coordinates using the file path directly:
    coord_by_filepath = get_body('earth', time, ephemeris=filepath)

    # Using the URL or filepath should give exactly the same results:
    assert_quantity_allclose(coord_by_url.ra, coord_by_filepath.ra)
    assert_quantity_allclose(coord_by_url.dec, coord_by_filepath.dec)
    assert_quantity_allclose(coord_by_url.distance, coord_by_filepath.distance)

@pytest.mark.remote_data
@pytest.mark.skipif('not HAS_JPLEPHEM')
def test_url_ephemeris_wrong_input():
    # Try loading a non-existing URL:
    time = Time('1960-01-12 00:00')
    with pytest.raises((HTTPError, URLError)):
        get_body('earth', time, ephemeris=get_pkg_data_filename('path/to/nonexisting/file.bsp'))

@pytest.mark.skipif('not HAS_JPLEPHEM')
def test_file_ephemeris_wrong_input():
    time = Time('1960-01-12 00:00')
    # Try loading a non-existing file:
    with pytest.raises(ValueError):
        get_body('earth', time, ephemeris='/path/to/nonexisting/file.bsp')

    # NOTE: This test currently leaves the file open (ResourceWarning).
    # To fix this issue, an upstream fix is required in jplephem
    # package.
    # Try loading a file that does exist, but is not an ephemeris file:
    with pytest.warns(ResourceWarning), pytest.raises(ValueError):
        get_body('earth', time, ephemeris=__file__)

def test_regression_10271():
    t = Time(58973.534052125986, format='mjd')
    # GCRS position of ALMA at this time
    obs_p = CartesianRepresentation(5724535.74068625, -1311071.58985697, -2492738.93017009, u.m)

    geocentre = CartesianRepresentation(0, 0, 0, u.m)
    icrs_sun_from_alma = _get_apparent_body_position('sun', t, 'builtin', obs_p)
    icrs_sun_from_geocentre = _get_apparent_body_position('sun', t, 'builtin', geocentre)

    difference = (icrs_sun_from_alma - icrs_sun_from_geocentre).norm()
    assert_quantity_allclose(difference, 0.13046941*u.m, atol=1*u.mm)

max_stars_repo_name: emirkmo/astropy

max_issues_repo_name: emirkmo/astropy

max_forks_repo_name: emirkmo/astropy

reason: import numpy,from astropy


text:
from numpy import sum
from numpy import zeros

from gwlfe.Memoization import memoize
from gwlfe.Output.Loading.StreamBankEros_1 import StreamBankEros_1
from gwlfe.Output.Loading.StreamBankEros_1 import StreamBankEros_1_f

def AvStreamBankEros(NYrs, DaysMonth, Temp, InitSnow_0, Prec, NRur, NUrb, Area, CNI_0, AntMoist_0, Grow_0, CNP_0, Imper,
                     ISRR, ISRA, CN, UnsatStor_0, KV, PcntET, DayHrs, MaxWaterCap, SatStor_0, RecessionCoef, SeepCoef,
                     Qretention, PctAreaInfil, n25b, Landuse, TileDrainDensity, PointFlow, StreamWithdrawal,
                     GroundWithdrawal, NumAnimals, AvgAnimalWt, StreamFlowVolAdj, SedAFactor_0, AvKF, AvSlope,
                     SedAAdjust, StreamLength, n42b, n46c, n85d, AgLength, n42, n45, n85, UrbBankStab):
    result = zeros(12)
    stream_bank_eros_f = StreamBankEros_1(NYrs, DaysMonth, Temp, InitSnow_0, Prec, NRur, NUrb, Area,
                                          CNI_0, AntMoist_0, Grow_0, CNP_0, Imper, ISRR, ISRA, CN, UnsatStor_0,
                                          KV, PcntET, DayHrs, MaxWaterCap, SatStor_0, RecessionCoef, SeepCoef
                                          , Qretention, PctAreaInfil, n25b, Landuse, TileDrainDensity, PointFlow,
                                          StreamWithdrawal, GroundWithdrawal, NumAnimals, AvgAnimalWt, StreamFlowVolAdj,
                                          SedAFactor_0, AvKF, AvSlope, SedAAdjust, StreamLength, n42b, n46c, n85d,
                                          AgLength, n42, n45, n85, UrbBankStab)
    for Y in range(NYrs):
        for i in range(12):
            result[i] += stream_bank_eros_f[Y][i] / NYrs
    return result

@memoize
def AvStreamBankEros_f(NYrs, DaysMonth, Temp, InitSnow_0, Prec, NRur, NUrb, Area, CNI_0, AntMoist_0, Grow_0, CNP_0,
                       Imper, ISRR, ISRA, CN, UnsatStor_0, KV, PcntET, DayHrs, MaxWaterCap, SatStor_0, RecessionCoef,
                       SeepCoef, Qretention, PctAreaInfil, n25b, Landuse, TileDrainDensity, PointFlow, StreamWithdrawal,
                       GroundWithdrawal, NumAnimals, AvgAnimalWt, StreamFlowVolAdj, SedAFactor_0, AvKF, AvSlope,
                       SedAAdjust, StreamLength, n42b, n46c, n85d, AgLength, n42, n45, n85, UrbBankStab):
    stream_bank_eros_f = StreamBankEros_1_f(NYrs, DaysMonth, Temp, InitSnow_0, Prec, NRur, NUrb, Area,
                                            CNI_0, AntMoist_0, Grow_0, CNP_0, Imper, ISRR, ISRA, CN, UnsatStor_0,
                                            KV, PcntET, DayHrs, MaxWaterCap, SatStor_0, RecessionCoef, SeepCoef
                                            , Qretention, PctAreaInfil, n25b, Landuse, TileDrainDensity, PointFlow,
                                            StreamWithdrawal, GroundWithdrawal, NumAnimals, AvgAnimalWt,
                                            StreamFlowVolAdj,
                                            SedAFactor_0, AvKF, AvSlope, SedAAdjust, StreamLength, n42b, n46c, n85d,
                                            AgLength, n42, n45, n85, UrbBankStab)
    return sum(stream_bank_eros_f, axis=0) / NYrs

ext: py

max_stars_repo_name: mudkipmaster/gwlf-e

max_issues_repo_name: mudkipmaster/gwlf-e

max_forks_repo_name: mudkipmaster/gwlf-e

reason: from numpy


text:
import os

import numpy as np
from IPython.display import IFrame
from matplotlib import pyplot as plt
from scipy.io import loadmat
from scipy.ndimage.filters import median_filter

from config import BASE_DIR, SHAPENET_IM
from utils import mkdir_p
from uuid import uuid4
with open(os.path.join(BASE_DIR, 'pyntcloud.js'), 'r') as f:
    TEMPLATE_POINTS = f.read()

def array_to_color(array, cmap="Oranges"):
    s_m = plt.cm.ScalarMappable(cmap=cmap)
    return s_m.to_rgba(array, norm=False)[:, :-1]

def plot_points(xyz,
                colors=None,
                size=0.1,
                axis=False,
                title=None,
                html_out=None):
    positions = xyz.reshape(-1).tolist()
    mkdir_p('vis')
    if html_out is None:
        html_out = os.path.join('vis', 'pts{:s}.html'.format(uuid4()))
    if title is None:
        title = "PointCloud"
    camera_position = xyz.max(0) + abs(xyz.max(0))

    look = xyz.mean(0)

    if colors is None:
        colors = [1, 0.5, 0] * len(positions)

    elif len(colors.shape) > 1:
        colors = colors.reshape(-1).tolist()

    if axis:
        axis_size = xyz.ptp() * 1.5
    else:
        axis_size = 0

    with open(html_out, "w") as html:
        html.write(
            TEMPLATE_POINTS.format(
                title=title,
                camera_x=camera_position[0],
                camera_y=camera_position[1],
                camera_z=camera_position[2],
                look_x=look[0],
                look_y=look[1],
                look_z=look[2],
                positions=positions,
                colors=colors,
                points_size=size,
                axis_size=axis_size))

    return IFrame(html_out, width=1024, height=768)

def image_grid(ims, mask=None):
    if mask is not None:
        ims[np.logical_not(mask)] = None
    gh, gw, h, w, ch = ims.shape
    disp_im = np.zeros([gh * h, gw * w, ch])
    for y in range(gh):
        for x in range(gw):
            disp_im[y * h:(y + 1) * h, x * w:(x + 1) * w, :] = ims[y][x]
    return disp_im

def voxel_grid(voxels, thresh=0.4, cmap='viridis'):
    gh, gw, h, w, d, ch = voxels.shape
    all_pts, all_clr = [], []
    for bx in range(gh):
        for ix in range(gw):
            pts, clr = voxels2pts(voxels[bx][ix], cmap=cmap)
            pts[:, 0] += ix * w
            pts[:, 1] += bx * h
            all_pts.append(pts)
            all_clr.append(clr)
        all_pts.append(pts)
        all_clr.append(clr)
    vis_pts, vis_clr = np.concatenate(
        all_pts, axis=0), np.concatenate(
            all_clr, axis=0)
    return vis_pts, vis_clr

def voxels2pts(voxels, thresh=0.4, cmap="Oranges"):
    if voxels.ndim == 4:
        fvox = voxels[..., 0]
    elif voxels.ndim == 3:
        fvox = voxels
    else:
        print('Invalid number of dimension in voxel grid')
        return
    vox = (fvox > thresh).astype(np.int)
    points = np.argwhere(vox > 0)
    colors = array_to_color(fvox[vox > 0], cmap=cmap)
    return points, colors

def voxel2mesh(voxels):
    cube_verts = [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0],
                  [1, 0, 1], [1, 1, 0], [1, 1, 1]]  # 8 points

    cube_faces = [[0, 1, 2], [1, 3, 2], [2, 3, 6], [3, 7, 6], [0, 2, 6],
                  [0, 6, 4], [0, 5, 1], [0, 4, 5], [6, 7, 5], [6, 5, 4],
                  [1, 7, 3], [1, 5, 7]]  # 12 face

    cube_verts = np.array(cube_verts)
    cube_faces = np.array(cube_faces) + 1

    l, m, n = voxels.shape

    scale = 1.0 / voxels.shape[0]
    cube_dist_scale = 1.1
    verts = []
    faces = []
    curr_vert = 0
    for i in range(l):
        for j in range(m):
            for k in range(n):
                # If there is a non-empty voxel
                if voxels[i, j, k] > 0:
                    verts.extend(
                        scale *
                        (cube_verts + cube_dist_scale * np.array([[i, j, k]])))
                    faces.extend(cube_faces + curr_vert)
                    curr_vert += len(cube_verts)

    verts = np.array(verts) - 0.5
    return verts, np.array(faces)

def write_obj(filename, verts, faces):
    """ write the verts and faces on file."""
    with open(filename, 'w') as f:
        # write vertices
        f.write('g\n# %d vertex\n' % len(verts))
        for vert in verts:
            f.write('v %f %f %f\n' % tuple(vert))

        # write faces
        f.write('# %d faces\n' % len(faces))
        for face in faces:
            f.write('f %d %d %d\n' % tuple(face))

def voxel2obj(filename, pred):
    verts, faces = voxel2mesh(pred)
    write_obj(filename, verts, faces)

def unproject_depth(d_im, K, R, im=None, dmin=1.0, dmax=3.0):
    px, py, f = K[0, 2], K[1, 2], K[0, 0]
    size = d_im.shape
    x, y = np.meshgrid(range(size[0]), range(size[1]))
    x, y = (x - px) * d_im / f, (y - py) * d_im / f
    xyz = np.stack([x, y, d_im], axis=0)
    xyz = np.reshape(xyz, [3, -1])
    mask = np.logical_and(xyz[-1, :] < dmax, xyz[-1, :] > dmin)
    xyz = xyz[:, mask]
    clr = None
    if im is not None:
        im = np.transpose(im, [2, 0, 1])
        im = np.reshape(im, [3, -1])
        clr = im[:, mask]

    tr = -np.matmul(R[:3, :3].T, R[:, 3][..., np.newaxis])
    Rt = np.concatenate([R[:3, :3].T, tr], axis=1)
    Xw = np.matmul(Rt,
                   np.concatenate([xyz, np.ones((1, xyz.shape[1]))], axis=0))
    return np.transpose(Xw), np.transpose(clr)

def depth2mesh(classId,
               shapeId,
               dmap,
               im_idx,
               dmin=1.15,
               dmax=2.85,
               discThresh=0.035,
               smooth=True,
               obj=None):
    shapeNetFolder = os.path.join(SHAPENET_IM, classId, shapeId)
    camera_f = os.path.join(shapeNetFolder, 'camera_{:d}.mat'.format(im_idx))
    mat = loadmat(camera_f)
    K = mat['K']
    px, py, f = K[0, 2], K[1, 2], K[0, 0]
    R = mat['extrinsic'][:3, :]
    tr = -np.matmul(R[:3, :3].T, R[:, 3][..., np.newaxis])
    Rt = np.concatenate([R[:3, :3].T, tr], axis=1)
    if obj is None:
        out_f = 'depth/mesh_{}_{}_{}.obj'.format(classId, shapeId, im_idx)
    else:
        out_f = obj

    depthMap = dmap
    if smooth:
        depthMap = median_filter(depthMap, (3, 3))

    h, w = depthMap.shape
    allPoints = np.ndarray(shape=(h, w, 5))

    with open(out_f, 'w') as obj:
        ind = 1
        for y in range(0, h):
            for x in range(0, w):
                d_im = depthMap[y, x]
                allPoints[y, x, 3] = depthMap[y, x]
                if (d_im < dmax and d_im > dmin):
                    x_c, y_c = (x - px) * d_im / f, (y - py) * d_im / f
                    upointLocal = np.array([[x_c], [y_c], [d_im], [1]])
                    upoint = np.matmul(Rt, upointLocal)
                    allPoints[y, x, :3] = upoint[:3, 0]
                    allPoints[y, x, 4] = ind
                    ind = ind + 1
                    obj.write('v {} {} {}\n'.format(upoint[0, 0], upoint[1, 0],
                                                    upoint[2, 0]))

        for y in range(0, h - 1):
            for x in range(0, w - 1):
                v = allPoints[y, x, 4]
                vd = allPoints[y, x, 3]
                vx = allPoints[y, x + 1, 4]
                vxd = allPoints[y, x + 1, 3]
                vy = allPoints[y + 1, x, 4]
                vyd = allPoints[y + 1, x, 3]
                vxy = allPoints[y + 1, x + 1, 4]
                vxyd = allPoints[y + 1, x + 1, 3]

                t1_minD = min(vd, vxd, vyd)
                t1_maxD = max(vd, vxd, vyd)

                if (t1_minD > dmin and t1_maxD < dmax and
                        t1_maxD - t1_minD < discThresh):
                    obj.write(
                        'f {:d} {:d} {:d}\n'.format(int(vy), int(vx), int(v)))

                t2_minD = min(vxyd, vxd, vyd)
                t2_maxD = max(vxyd, vxd, vyd)

                if (t2_minD > dmin and t2_maxD < dmax and
                        t2_maxD - t2_minD < discThresh):
                    try:
                        obj.write('f {:d} {:d} {:d}\n'.format(
                            int(vxy), int(vx), int(vy)))
                    except:
                        print 'Error', vxy, vx, vy
                        return
    return out_f

max_issues_repo_name: Ajithbalakrishnan/lsm

max_forks_repo_name: Ajithbalakrishnan/lsm

reason: import numpy,from scipy


text:
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import sys
import numpy
import os

os.environ["CC"] = "g++"
os.environ["CXX"] = "g++"

ext_modules = [Extension("average_trajectory_cython", ["average_trajectory_cython.pyx"], language='c++')]

setup(
    name='average_trajectory_cython',
    cmdclass={'build_ext': build_ext},
    include_dirs=[numpy.get_include()],
    ext_modules=ext_modules,
)

max_stars_repo_name: amirdel/dispersion-continua

max_issues_repo_name: amirdel/dispersion-continua

max_forks_repo_name: amirdel/dispersion-continua

reason: import numpy


text:
import numpy as np
import matplotlib as plt
import pandas as pd
import copy
from constants import *
from DispData import *
import matplotlib.pyplot as plt
from refl_disp import *
from scipy.optimize import minimize

"""
This class optimizes the spectrum of a multi-layer structure according to a target spectrum.
"""

#%%
class MAoptim:
    def __init__(self, disp_data, list_material_names, material_in, material_sub, num_layers, num_population, r_target_input, \
                 theta_in=0, d_max=DEFAULT_D_MAX, max_iter=DEFAULT_MAX_ITER):
        """
        Constructor for generating a multi-layer optimization instance

        Parameters:

            disp_data: an object from DispData containing the optimization wavelengths and materials of interest
            list_material_names: a list of material names used in the multi-layer
            material_in: material of incidence medium
            material_sub: material of substrate medium
            num_layers: number of layers
            num_population: size of the memetic algorithm population
            r_input_target: array that contains the target spectrum's information
                format: array([[wvlen_min1, wvlen_max1, R, weight], ...])
            theta_in: incidence angle
            d_max: maximum layer thickness
            max_iter: maximum number of memetic iterations

        Attributes:
            disp_data: an object from DispData
            N: size of the memetic algorithm population
            K: number of layers
            d_max: maximum layer thickness

            list_material_names: a list of material names used in the multi-layer
            num_materials: number of materials used
            list_materials: list_material_names represented by integers
            material_in: material of incidence medium
            material_sub: material of substrate medium
            n_in: refractive index incidence material
            n_in_fine: finely sampled refractive index incidence material
            n_sub: refractive index of substrate
            n_sub_fine: finely sampled refractive index of substrate
            theta_in: incidence angle

            r_target_input: array input of target reflection spectrum
                format: array([[wvlen_min1, wvlen_max1, R, weight], ...])
            merit_function: merit function
            r_target: target reflection spectrum

            max_iter: maximum number of memetic iterations
            rate_top: reselection rate - top percentages to keep
            rate_good: reselection rate - percentages selected from the top half of population
            mutation_rate: mutation rate of population
            refine_period: number of iterations between population refinement
            refine_num: number of individuals refined

            fitness_generation: fitness of the current generation (size 2N)
            fitness_best: best fitness of the generation
            fitness_avg: average fitness of the generation
            fitness_ind: sorted fitness indices according to structure
            fitness_best_so_far: best fitness of the algorithm so far

            n_material_best: materials of the best structure (size K)
            d_best: thicknesses of the best structure (size K)
            n_material_generation: materials of each generation (size 2N by K)
            d_generation: thicknesses of each generation (size 2N by K)
            n_material_parents: materials of parents (size N by K)
            d_parents: thicknesses of parents (size N by K)
            n_material_childern: materials of childern (size N by K)
            d_childern: thicknesses of childern (size N by K)

            optim_R_spectrum_fine: spectrum of the optimum structure

        """

        #%% Initialize a memetic algorithm's population
        self.disp_data = disp_data
        self.N = num_population # N = number of population size
        self.K = num_layers # K = number of layers
        self.list_material_names = list_material_names
        self.num_materials = len(self.list_material_names)
        self.list_materials = self.matname2list(self.list_material_names)
        self.d_max = d_max

        #%% Incidence and substrate materials
        self.material_in = material_in
        self.material_sub = material_sub

        self.n_in, self.n_in_fine = self.material_to_ind(self.material_in)
        self.n_sub, self.n_sub_fine = self.material_to_ind(self.material_sub)

        self.theta_in = theta_in;

        #%% Initialize the merit function
        self.r_target_input = r_target_input # [[wvlen_min1, wvlen_min2, R, weight], ...]
        self.merit_function, self.r_target = self.generate_mf()

        #%% Initialize the memetic operator values
        self.max_iter = max_iter
        self.rate_top = DEFAULT_RATE_TOP
        self.rate_good = DEFAULT_RATE_GOOD
        self.mutation_rate = DEFAULT_MUTATION_RATE
        self.refine_period = DEFAULT_REFINE_PERIOD
        self.refine_num = DEFAULT_REFINE_NUM

        #%% Initialize state variables
        self.fitness_generation = []
        self.fitness_best = np.zeros(self.max_iter)
        self.fitness_avg = np.zeros(self.max_iter)
        self.fitness_ind = []
        self.fitness_best_so_far = 1e9;

        self.n_material_best = []
        self.d_best = []

        self.n_material_generation = []
        self.d_generation = []

        self.n_material_parents = []
        self.d_parents = []

        self.n_material_childern = []
        self.d_childern = []

        self.optim_R_spectrum_fine = np.zeros([disp_data.N_wvlen_fine])

    #%% Convert material names to integers
    def matname2list(self, list_material_names):
        """
        Converts material names to integer representations.

        Parameters:
            list_material_names: list of material names

        Output:
            list_materials: list of integers, each representing a material

        """

        list_materials = []

        for i in range(self.num_materials):
            list_materials.append(self.disp_data.material_key.index(list_material_names[i]))

        return list_materials

    #%% Function that generates the population
    def generate_population(self):
        """
        Generates a population for memetic algorithm

        This method updates self.n_material_parents and self.d_parents
        so that it produces N number of K-layer structures randomly

        """

        # num_population, num_layers, list_materials, d_max
        max_index = len(self.list_materials)
        temp_ordering = np.random.randint(0, max_index, [self.N, self.K])

        self.n_material_parents = np.zeros((self.N, self.K)).astype(int)
        for i in range(self.N):
            for j in range(self.K):
                self.n_material_parents[i, j] = self.list_materials[temp_ordering[i, j]]

        self.d_parents = self.d_max/2 * np.ones((self.N, self.K)) \
                      + self.d_max * (np.random.rand(self.N, self.K) - 0.5)

        return

    #%% Function that performs crossover
    def crossover(self):
        """
        Performs crossover amongst parents

        This method randomly pairs up two parent structures, perform crossover
        to produce two new childern structures. This doubles the size of the
        individuals in the generation.

        """

        # Initialize childern arrays
        self.n_material_childern = np.empty([self.N, self.K]).astype(int)
        self.d_childern = np.empty([self.N, self.K])

        # Shuffle the layers
        shuffle_index = np.random.permutation(self.N)

        self.n_material_parents = self.n_material_parents[shuffle_index, :]
        self.d_parents = self.d_parents[shuffle_index, :]

        # Perform crossover for each pair of parents
        for j in range(int(self.N/2)):
            n_material_parent1 = copy.deepcopy(self.n_material_parents[2*j, :])
            n_material_parent2 = copy.deepcopy(self.n_material_parents[2*j+1, :])

            d_parent1 = copy.deepcopy(self.d_parents[2*j, :])
            d_parent2 = copy.deepcopy(self.d_parents[2*j+1, :])

            # Choose a random crossover point
            cross_pt = np.random.randint(self.K-1)+1

            self.n_material_childern[2*j, :] = np.concatenate((n_material_parent1[0:cross_pt], n_material_parent2[cross_pt:self.K]), axis=0)
            self.n_material_childern[2*j+1, :] = np.concatenate((n_material_parent2[0:cross_pt], n_material_parent1[cross_pt:self.K]), axis=0)

            self.d_childern[2*j, :] = np.concatenate((d_parent1[0:cross_pt], d_parent2[cross_pt:self.K]), axis=0)
            self.d_childern[2*j+1, :] = np.concatenate((d_parent2[0:cross_pt], d_parent1[cross_pt:self.K]), axis=0)

        # Generation now double in size
        self.n_material_generation = np.concatenate((self.n_material_parents, self.n_material_childern), axis=0)
        self.d_generation = np.concatenate((self.d_parents, self.d_childern), axis=0)

        return

    #%%
    def mutation(self):
        """
        Performs mutations on a small number of individuals

        This method randomly selects individuals with probability mutation_rate
        to perform mutation. When mutation happens to an individual, it changes
        a random layer's material and thickness

        """

        for n in range(2*self.N):
            if np.random.rand() < self.mutation_rate:
                layer_change = np.random.randint(0, self.K)
                self.d_generation[n, layer_change] = self.d_max/2. + self.d_max * (np.random.rand() - 0.5)

                index_max = len(self.list_materials)
                temp_index = np.random.randint(0, index_max)
                self.n_material_generation[n, layer_change] = self.list_materials[temp_index]
        return

    #%% Evaluate the fitness of the population
    def eval_fitness(self):
        """
        Evaluates the fitness of the individuals in a generation

        """

        if self.fitness_generation == []:
            self.fitness_generation = np.zeros([2*self.N])

        for n in range(2*self.N):
            n_generation = self.disp_data.ref_ind[:, self.n_material_generation[n, :]]
            self.fitness_generation[n] = self.merit_function(n_generation, self.d_generation[n, :])

        return

    #%%
    def sort_population(self):
        """
        Sorts the individuals in a generation according to their merit function.

        """
        self.fitness_ind =sorted(range(2*self.N), key = lambda k : self.fitness_generation[k])
        self.fitness_ind = np.array(self.fitness_ind)

        # Arrange materials according to fitness_ind (input passed by reference)
        self.n_material_generation = self.n_material_generation[self.fitness_ind, :]
        self.d_generation = self.d_generation[self.fitness_ind, :]
        self.fitness_generation = self.fitness_generation[self.fitness_ind]

        return

    #%%
    def reselection(self):
        """
        Performs reselection on a generation's population

        This method selects the parent individuals for the next generation of
        memetic optimization according to the reselection rule.

        Reselection rule:
            Keep top rate_top of the population as parents
            Select the rest of the rate_good of the parents from the top-half of population
            Select the remaining rate_poor of the parents from the bottom half

        """
        #
        N_top = round(self.N * self.rate_top)
        N_good = round(self.N * self.rate_good - self.rate_top)
        N_poor = self.N - N_top - N_good

        N_good_pool = self.N - N_top

        # Keep the N_top best structures
        selected_ind_top = np.array(range(N_top))

        # Select from the good structures
        temp_rand_index = N_top + np.random.permutation(N_good_pool)
#        temp_selected_good = temp_rand_index[0:N_good]
#        selected_ind_good = self.fitness_ind[temp_selected_good]
        selected_ind_good = temp_rand_index[0:N_good]

        # Select from the poor structures
        temp_rand_index = np.random.permutation(self.N) + self.N
#        selected_ind_poor = self.fitness_ind[temp_rand_index[0:N_poor]]
        selected_ind_poor = temp_rand_index[0:N_poor]

        # Concatenate the selected indices
        selected_inds = np.concatenate((selected_ind_top, selected_ind_good, selected_ind_poor)).astype(int)

#        print(selected_inds)
        self.n_material_parents = copy.deepcopy(self.n_material_generation[selected_inds, :])
        self.d_parents = copy.deepcopy(self.d_generation[selected_inds, :])

        return

    #%% Keep and display the best structure
    def keep_best(self, iter):
        """
        Keeps the best structures up to the current optimization step

        """

        self.fitness_best[iter] = self.fitness_generation[0]
        self.fitness_avg[iter] = np.mean(self.fitness_generation)

        if self.fitness_best[iter] < self.fitness_best_so_far:
            print('Found a temporary best structure in iteration ', iter)

            self.n_material_best = self.n_material_generation[0, :]
            self.d_best = self.d_generation[0, :]

            # Display the best structure so far and its fitness
            self.fitness_best_so_far = self.fitness_best[iter]
            print('Best fitness so far is ', self.fitness_best_so_far)
            self.print_struct(self.n_material_best, self.d_best)
            print('\n')

    #%% Perform local optimization
    def refinement(self, iter):
        """
        Performs refinement of a few structures

        This method performs a local optimization on the thicknesses of some
        structures periodically.

        """

        if np.mod(iter, self.refine_period) == self.refine_period-1:
            fitness_parents = np.zeros([self.N])
            # Resort fitness
            for n in range(self.N):
                n_parents = self.disp_data.ref_ind[:, self.n_material_parents[n, :]]
                fitness_parents[n] = self.merit_function(n_parents, self.d_parents[n, :])

            fitness_ind_parents = sorted(range(self.N), key = lambda k : fitness_parents[k])
            self.n_material_parents = self.n_material_parents[fitness_ind_parents, :]
            self.d_parents = self.d_parents[fitness_ind_parents, :]

            # Perform local optimization
            for m in range(self.refine_num):
                n_layers_m = self.disp_data.ref_ind[:, self.n_material_parents[m, :]]

                obj_fun = lambda x: self.merit_function(n_layers_m, x)

                x0 = copy.deepcopy(self.d_parents[m, :])

                x_opt = minimize(obj_fun, x0)
                # print(x_opt.x)

                self.d_parents[m, :] = copy.deepcopy(x_opt.x)
                print('Merit function of individual = ', obj_fun(x_opt.x))

    #%%  Print the material
    def print_struct(self, n_index_layers, d_layers):
        """
        This method prints a structure's material composition and thicknesses.

        Parameters:
            n_index_layers: an array of integers representing materials
            d_layers: thickness of each layer

        """

        struct_name = []

        for i in range(self.K):
            struct_name.append(self.disp_data.material_key[n_index_layers[i]])

        d = {'material' : struct_name, 'thickness' : d_layers}
        df = pd.DataFrame(data=d)

        print(df)

    #%% Set optimization rates

    #%% Public function: optimize the population given
    def ma_optimize(self):
        """
        Performs the memetic algorithm optimization

        This method calls all the memetic algorithm operators to perform
        optimization and obtain the optimum structure

        """

        #%% Convergence curve
        fig_conv = plt.figure()
        ax_conv = fig_conv.add_subplot(111)

        for i in range(self.max_iter):
            if np.mod(i, 5) == 0:
                print('iter = ', i, '\n')

            # Single point crossover
            self.crossover()

            # Perform mutation
            self.mutation()

            # Evaluate fitness and sort the population
            self.eval_fitness()
            self.sort_population()

            # Keep best structure
            self.keep_best(i)

            # Perform reselection
            self.reselection()

            # Plot
            self.plot_convergence(i, fig_conv)

            # Refinement
            self.refinement(i)

        return []

    #%%
    def plot_spectrum(self):
        """
        Plots the spectrum of the most optimum structure

        """

        n_material_final = self.n_material_best
        d_final = self.d_best

        print('The optimum structure found was: ')
        self.print_struct(n_material_final, d_final)

        n_final_fine = self.disp_data.ref_ind_fine[:, n_material_final]

        # Display the best structure (including the fixed materials above the substrate)
#        structure_best = convert2struct(material_key, n_material_final, d_final)

        # Compute the optimized structure's spectrum
        self.optim_R_spectrum_fine = refl_disp(self.n_in_fine, self.n_sub_fine, n_final_fine, d_final, self.disp_data.wvlen_vec_fine, self.theta_in)

        # Plot the optimized spectrum
        fig_spect = plt.figure()
        ax_spect = fig_spect.add_subplot(111)
        plt.plot(self.disp_data.wvlen_vec_fine, self.optim_R_spectrum_fine, 'b', label='Optimized spectrum')
        plt.plot(self.disp_data.wvlen_vec, self.r_target, 'g.', label='Target')

        plt.axis([min(self.disp_data.wvlen_vec_fine), max(self.disp_data.wvlen_vec_fine), -0.01, 1.01])
        plt.xlabel('wavelength ($\mu$m)')
        plt.ylabel('reflectivity')
        plt.legend()

        plt.show()

    #%%
    def plot_r_target(self):
        """
        Plot the target reflection spectrum

        """
        fig_tgt = plt.figure()
        ax_tgt = fig_tgt.add_subplot(111)
        plt.plot(self.disp_data.wvlen_vec, self.r_target, 'g.', label='Target')
        plt.axis([min(self.disp_data.wvlen_vec_fine), max(self.disp_data.wvlen_vec_fine), -0.01, 1.01])
        plt.xlabel('wavelength ($\mu$m)')
        plt.ylabel('reflectivity')
        plt.legend()

        plt.show()

    #%% Plot the convergence
    def plot_convergence(self, i, fig_conv):
        """
        Plots the convergence curve

        This method plots the convergence of the best fitness and average
        fitness of the population over each generation.

        Parameters:
            i: current iteration
            fig_conv: figure that refers to the convergence plot

        """

        if i == 1:
            plt.semilogy([i-1, i], [self.fitness_best[i-1], self.fitness_best[i]], 'b', label='Best merit function')
            plt.semilogy([i-1, i], [self.fitness_avg[i-1], self.fitness_avg[i]], 'r--', label='Average merit function')
            plt.legend(loc='upper right')
            fig_conv.show()
            plt.pause(0.05)

        elif i > 1:
            plt.semilogy([i-1, i], [self.fitness_best[i-1], self.fitness_best[i]], 'b')
            plt.semilogy([i-1, i], [self.fitness_avg[i-1], self.fitness_avg[i]], 'r--')
            fig_conv.show()
            plt.pause(0.05)

        return

    #%% function
    def generate_mf(self):
        """
        Generates the merit function

        This method generates the a function handle that represents the merit
        function of a given structure according to its material composition and
        layer thicknesses.

        Output:
            mf: merit function handle
            r_target: target spectrum

        """

        r_target = np.zeros([self.disp_data.N_wvlen])
        weight = np.zeros([self.disp_data.N_wvlen])

        if len(self.r_target_input.shape) == 1:
            wvlen_min_temp = self.r_target_input[0]
            wvlen_max_temp = self.r_target_input[1]
            r_temp = self.r_target_input[2]
            w_temp = self.r_target_input[3]

            indicator_interval = (self.disp_data.wvlen_vec >= wvlen_min_temp) & (self.disp_data.wvlen_vec < wvlen_max_temp)
            r_target = r_temp * indicator_interval
            weight = w_temp * indicator_interval

        else:
            N_intervals = len(self.r_target_input[:, 0])

            for i in range(N_intervals):
                # Parse input parameters
                wvlen_min_temp = self.r_target_input[i, 0]
                wvlen_max_temp = self.r_target_input[i, 1]
                r_temp = self.r_target_input[i, 2]
                w_temp = self.r_target_input[i, 3]

                # Generate the target spectrum and weight function
                indicator_interval = (self.disp_data.wvlen_vec >= wvlen_min_temp) & (self.disp_data.wvlen_vec < wvlen_max_temp)
                r_target += r_temp * indicator_interval
                weight += w_temp * indicator_interval

        mf = lambda n_layers, d: np.linalg.norm(weight * (refl_disp(self.n_in, self.n_sub, n_layers, abs(d), self.disp_data.wvlen_vec, self.theta_in) - r_target))**2

        return mf, r_target

    #%% This function converts a material name to its dispersive refractive indices
    def material_to_ind(self, material_input):
        """
        Converts a material's name to its refractive indices

        Parameters:
            material_input: a single material input name

        Outputs:
            ref_ind_output: refractive index of the material_input
            ref_ind_output_fine: finely sampled refractive index of the material input

        """

        if type(material_input) == str:
            material_index = self.disp_data.material_key.index(material_input)
            ref_ind_output = self.disp_data.ref_ind[:, material_index]
            ref_ind_output_fine = self.disp_data.ref_ind_fine[:, material_index]

        elif type(material_input) == float or type(material_input) == int:
            ref_ind_output = material_input
            ref_ind_output_fine = material_input

        return ref_ind_output, ref_ind_output_fine

    #%% Set the memetic optimization parameters
    def set_optim_params(self, **kwargs):
        """
        Adjust the optimization parameter to users' input

        Parameters:
            **kwargs: key item that contains {'optim_param', value}
            'optim_params' include:
                'max_iter'
                'rate_top'
                'rate_good'
                'mutation_rate'
                'refine_period'
                'refine_num'

        """

        for key, value in kwargs.items():
            if key == 'max_iter':
                self.max_iter = value
            elif key == 'rate_top':
                self.rate_top = value
            elif key == 'rate_good':
                self.rate_good = value
            elif key == 'mutation_rate':
                self.mutation_rate = value
            elif key == 'refine_period':
                self.refine_period = value
            elif key == 'refine_num':
                self.refine_num = value
            else:
                raise ValueError('"%s" is not a valid optimization parameter' % key)


text:
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import MultinomialNB
from sklearn import svm
from sklearn.feature_extraction import DictVectorizer
import gensim
from gensim.models import KeyedVectors
import pandas as pd
import numpy as np
import sys
import csv

# Train models with features represented as one-hot encodings

feature_to_index = {'token': 0, 'pos': 1, 'chunk': 2, 'previous': 4, 'following': 5, 'capital': 6}

def extract_features_and_labels(trainingfile, selected_features):
    '''
    Extract features and gold labels from a preprocessed file with the training data and return them as lists

    :param trainingfile: path to training file
    :param selected_features: list of features that will be used to train the model

    :type trainingfile: string
    :type selected_features: list

    :return features: features as a list of dictionaries
    :return gold_labels: list of gold labels
    '''
    features = []
    gold_labels = []

    conllinput = open(trainingfile, 'r')
    csvreader = csv.reader(conllinput, delimiter='\t',quotechar='|')

    for row in csvreader:
        feature_value = {}
        # Only extract the selected features
        for feature_name in selected_features:
            row_index = feature_to_index.get(feature_name)
            feature_value[feature_name] = row[row_index]
        features.append(feature_value)

        # Gold is in the third column
        gold_labels.append(row[3])

    return features, gold_labels

def extract_features(testfile, selected_features):
    '''Extract features from a preprocessed file with the test data and return them as a list

    :param trainingfile: path to test file
    :param selected_features: list of features that were selected to train the model

    :type testfile: string
    :type selected_features: list

    :return features: features as a list of dictionaries'''

    features = []
    gold_labels = []

    conllinput = open(testfile, 'r')
    csvreader = csv.reader(conllinput, delimiter='\t',quotechar='|')

    for row in csvreader:
        feature_value = {}
        for feature_name in selected_features:
            row_index = feature_to_index.get(feature_name)
            feature_value[feature_name] = row[row_index]
        features.append(feature_value)

    return features

def create_classifier(train_features, train_labels, modelname):
    '''Create a classifier and train it with vectorized features and corresponding gold labels

    input train_features: features to be transformed into vectors
    input train_labels: gold labels corresponding to features
    input modelname: name of the model that will be trained

    output model: trained classifier
    output vec: DictVectorizer'''

    if modelname ==  'logreg':
        model = LogisticRegression(max_iter=10000)
    elif modelname == 'NB':
        model = MultinomialNB()
    elif modelname == 'SVM':
        model = svm.LinearSVC(max_iter=10000)

    vec = DictVectorizer()

    features_vectorized = vec.fit_transform(train_features)
    model.fit(features_vectorized, train_labels)

    return model, vec

def classify_data(model, vec, inputdata, outputfile, selected_features):
    '''Make predictions on test data and write the results to a file

    input model: classifier that will make predictions
    input vec: DictVectorizer object to transform the features into vectors
    input inputdata: path to input data
    input outputfile: path to output file where the predictions for each feature will be written
    input selected_features: list of features that were selected to train the model

    output: external file with results'''

    # Extract features from inputdata
    features = extract_features(inputdata, selected_features)
    features = vec.transform(features)

    # Make predictions
    predictions = model.predict(features)

    # Write results to an outputfile
    outfile = open(outputfile, 'w')
    counter = 0
    for line in open(inputdata, 'r'):
        if len(line.rstrip('\n').split()) > 0:
            outfile.write(line.rstrip('\n') + '\t' + predictions[counter] + '\n')
            counter += 1
    outfile.close()

# Train model with only word embeddings as features

def extract_embeddings_as_features_and_gold(trainingfile, word_embedding_model):
    '''Extracts features and gold labels from file with training data using word embeddings

    input trainingfile: path to training file
    input word_embedding_model: a pretrained word embedding model

    output features: list of vector representation of tokens
    output gold_labels: list of gold labels'''

    features = []
    gold_labels = []

    conllinput = open(trainingfile, 'r')
    csvreader = csv.reader(conllinput, delimiter='\t',quotechar='|')

    for row in csvreader:
        if row[0] in word_embedding_model:
            vector = word_embedding_model[row[0]]
        else:
            vector = [0]*300
        features.append(vector)

        # Gold is in the third column
        gold_labels.append(row[3])

    return features, gold_labels

def extract_embeddings_as_features(testfile, word_embedding_model):
    '''Extracts features from file with test data using word embeddings

    input testfile: path to conll file
    input word_embedding_model: a pretrained word embedding model

    output features: list of vector representation of tokens'''

    features = []

    conllinput = open(testfile, 'r')
    csvreader = csv.reader(conllinput, delimiter='\t',quotechar='|')

    for row in csvreader:
        if row[0] in word_embedding_model:
            vector = word_embedding_model[row[0]]
        else:
            vector = [0]*300
        features.append(vector)

    return features

def create_classifier_embeddings(train_features, train_labels):
    '''Create an SVM classifier and train it with vectorized features and corresponding gold labels

    input train_features: features to be transformed into vectors
    input train_labels: gold labels corresponding to features

    output model: trained classifier'''

    model = svm.LinearSVC(max_iter=10000)
    model.fit(train_features, train_labels)

    return model

def classify_data_embeddings(model, inputdata, outputfile, word_embedding_model):
    '''Let a classifier make predictions on new data

    input model: classifier that will make predictions
    input inputdata: path to input data
    input outputfile: path to output file, where the predictions for each feature will be written
    input word_embed'''

    # Extract features from inputdata
    features = extract_embeddings_as_features(inputdata, word_embedding_model)

    # Make predictions
    predictions = model.predict(features)

    # Write results to an outputfile
    outfile = open(outputfile, 'w')
    counter = 0
    for line in open(inputdata, 'r'):
        if len(line.rstrip('\n').split()) > 0:
            outfile.write(line.rstrip('\n') + '\t' + predictions[counter] + '\n')
            counter += 1
    outfile.close()

# Train mixed model (with embeddings and features represented as one-hot)

def combine_sparse_and_dense_features(dense_vectors, sparse_features):
    '''Function that takes sparse and dense feature representations and appends their vector representation

    input dense_vectors: list of dense vector representations
    input sparse_features: list of sparse vector representations

    output: list of arrays in which sparse and dense vectors are concatenated '''

    combined_vectors = []
    sparse_vectors = np.array(sparse_features.toarray())

    for index, vector in enumerate(sparse_vectors):
        combined_vector = np.concatenate((vector,dense_vectors[index]))
        combined_vectors.append(combined_vector)

    return combined_vectors

def extract_traditional_features_and_embeddings_plus_gold_labels(conllfile, word_embedding_model, vectorizer=None):
    '''Function that extracts traditional features as well as embeddings and gold labels using word embeddings for current and preceding token

    input conllfile: path to conll file
    input word_embedding_model: a pretrained word embedding model

    output features: list of vector representation of tokens
    output labels: list of gold labels'''

    # Get the word embeddings and a selection of traditional features plus gold labels using the functions from before
    dense_vectors = extract_embeddings_as_features(conllfile, word_embedding_model)
    traditional_features,labels = extract_features_and_labels(conllfile, ['pos','chunk','capital'])

    # Create vector representation of traditional features
    if vectorizer is None:
        vectorizer = DictVectorizer()
        vectorizer.fit(traditional_features)
    sparse_features = vectorizer.transform(traditional_features)
    combined_vectors = combine_sparse_and_dense_features(dense_vectors, sparse_features)

    return combined_vectors, vectorizer, labels

def label_data_with_combined_features(testfile, outputfile, classifier, vectorizer, word_embedding_model):
    '''Labels data with model using both sparse and dense features'''

    # Extract features and gold labels from inputdata
    feature_vectors, vectorizer, goldlabels = extract_traditional_features_and_embeddings_plus_gold_labels(testfile, word_embedding_model,vectorizer)

    # Make predictions
    predictions = classifier.predict(feature_vectors)

    # Write the results to an output file
    outfile = open(outputfile, 'w')
    counter = 0
    for line in open(testfile, 'r'):
        if len(line.rstrip('\n').split()) > 0:
            outfile.write(line.rstrip('\n') + '\t' + predictions[counter] + '\n')
            counter += 1
    outfile.close()

def main():

    argv = sys.argv

    trainingfile = argv[1] #path to training file
    inputfile = argv[2] #path to dev file that will be used to generate output
    outputfile = argv[3] # name of output file that will be created
    type_of_system = argv[4] # 0 == traditional, 1==embedding, 2==mixed
    language_model = '../../models/GoogleNews-vectors-negative300.bin.gz' #if at some point there are others, add as argument argv[5]

    if type_of_system == '0': #traditional_system

        # Select features to train the model (update list for feature ablation)
        selected_features = ['token', 'pos', 'chunk', 'previous', 'following', 'capital']

        # Get the selected training features and gold labels
        training_features, gold_labels = extract_features_and_labels(trainingfile, selected_features)

        # Train three different models with the features, let them classify the data and write the result to new conll files
        for modelname in ['logreg', 'NB', 'SVM']:
            ml_model, vec = create_classifier(training_features, gold_labels, modelname)
            classify_data(ml_model, vec, inputfile, outputfile.replace('.conll','.' + modelname + '.conll'), selected_features)

            df = pd.read_table(outputfile.replace('.conll','.' + modelname + '.conll'))
            df = df.set_axis([*df.columns[:-1], 'NER2'], axis=1, inplace=False)
            df.to_csv(outputfile.replace('.conll','.' + modelname + '.conll'), sep='\t')

    elif type_of_system == '1': #embeddings_system

        # Load word embedding model
        language_model = gensim.models.KeyedVectors.load_word2vec_format(language_model, binary=True)

        # Get training features and gold labels
        training_features, gold_labels = extract_embeddings_as_features_and_gold(trainingfile, language_model)

        # Train an SVM classifier with the features, let it classify the data and write result to conll file
        ml_model = create_classifier_embeddings(training_features, gold_labels)
        classify_data_embeddings(ml_model, inputfile, outputfile.replace('.conll','.embed.conll'), language_model)

        df = pd.read_table(outputfile.replace('.conll','.embed.conll'))
        df = df.set_axis([*df.columns[:-1], 'NER2'], axis=1, inplace=False)
        df.to_csv(outputfile.replace('.conll','.embed.conll'), sep='\t')

    elif type_of_system == '2': #mixed_system

        # Load word embedding model
        language_model = gensim.models.KeyedVectors.load_word2vec_format(language_model, binary=True)

        # Get the combined feature vectors for training, vectorizer and gold labels
        combined_vectors, vectorizer, gold_labels = extract_traditional_features_and_embeddings_plus_gold_labels(inputfile, language_model, vectorizer=None)

        # Train an SVM classifier with the features, let it classify the data and write result to conll file
        ml_model = create_classifier_embeddings(combined_vectors, gold_labels)
        label_data_with_combined_features(inputfile,outputfile.replace('.conll','.mixed.conll'),ml_model,vectorizer,language_model)

        df = pd.read_table(outputfile.replace('.conll','.mixed.conll'))
        df = df.set_axis([*df.columns[:-1], 'NER2'], axis=1, inplace=False)
        df.to_csv(outputfile.replace('.conll','.mixed.conll'), sep='\t')

if __name__ == '__main__':
    main()

#python3 optimized_ner_machine_learning.py ../../data/preprocessed/conll2003.train-preprocessed_with_features.conll ../../data/preprocessed/conll2003.dev-preprocessed_with_features.conll ../../models/output/output_system.conll 0
#python3 optimized_ner_machine_learning.py ../../data/preprocessed/conll2003.train-preprocessed_with_features.conll ../../data/preprocessed/conll2003.dev-preprocessed_with_features.conll ../../models/output/output_system.conll 1
#python3 optimized_ner_machine_learning.py ../../data/preprocessed/conll2003.train-preprocessed_with_features.conll ../../data/preprocessed/conll2003.dev-preprocessed_with_features.conll ../../models/output/output_system.conll 2

#for ablation tests we have 6 features in total. remove one by one from last until we only have first
#eg: remove the last one, AKA 'capital' has the name output_system_rm1.conll
#remove 'capital' and 'following' has the name output_system_rm2.conll


text:
import os,sys
import pandas as pd
import numpy as np
import json,time
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense,Activation, Dropout, BatchNormalization,Conv2D,Conv1D,Flatten,LSTM,MaxPool1D
from sklearn.preprocessing import LabelEncoder,LabelBinarizer
from sklearn.model_selection import train_test_split
from tensorflow.keras.callbacks import TensorBoard

class dataModel:
    def __init__(self,dataName):
        self.xTrain = 0
        self.yTrain = 0
        self.xEval = 0
        self.yEval = 0
        self.xTest = 0
        self.yTest = 0
        self.nClasses = 0

        self.model = 0

        self.loadData(dataName)

    def loadData(self,dataName):
        with open(dataName,'r') as fileIn:
            data = json.load(fileIn)
        self.stratifyData(data)

    def stratifyData(self,data):
        xTrain,xEval,yTrain,yEval =  train_test_split(data["trials"],data['labels'],test_size = 0.3,random_state=42,stratify=data['labels'])
        xEval,xTest,yEval,yTest =  train_test_split(xEval,yEval,test_size = 0.5,random_state=42,stratify=yEval)

        #encoder = LabelEncoder()
        encoder = LabelBinarizer()
        #transformedLabaelencoder.fit_transform(list(set(yTrain)))
        #transfomed_label = encoder.fit_transform(["dog", "cat", "bird"])
        self.xTrain = np.array(xTrain).reshape((len(xTrain),128,3))
        self.xEval = np.array(xEval).reshape((len(xEval),128,3))
        self.xTest = np.array(xTest).reshape((len(xTest),128,3))
        self.yTrain = encoder.fit_transform(yTrain)
        self.yEval = encoder.fit_transform(yEval)
        self.yTest = encoder.fit_transform(yTest)

        #after transforming to one hot label vectors, the length is equal to the amount of classes
        self.nClasses = len(self.yTest[0])

    def buildNet2D(self):
        model = Sequential()
        model.add(Conv2D(input_shape = (128,3,1,),kernel_size = (8,3), padding = "valid", strides = (2,1), filters = 3, activation = "relu"))
        model.add(BatchNormalization())
        #model.add(Dropout(0.5))
        model.add(Conv2D(kernel_size = (11,1), strides = (2,1), padding = "valid", filters = 1, activation = "elu"))
        model.add(BatchNormalization())
        model.add(Flatten())
        #model.add(Dropout(0.5))
        model.add(Dense(self.nClasses, activation = "softmax" ))
        model.compile(loss='categorical_crossentropy', optimizer = 'adam', metrics = ['accuracy'])
        self.model = model
        print(model.summary())

    def buildNets1D(self):

        models = []

        denseLayers = [0,1,2]
        CNNLayers = [1,2,3,4,5]
        filters = [2,4,8,16,32,64]
        for dense in denseLayers:
            for CNNLayer in CNNLayers:
                for filt in filters:
                    nameOfModel = "{}-conv-{}-filter-{}-dense-{}".format(CNNLayer,filt,dense,int(time.time()))
                    model = Sequential()
                    model.add(Conv1D(input_shape = (128,3,),kernel_size = (3), padding = "valid", filters = filt))
                    model.add(Activation("elu"))
                    model.add(BatchNormalization())
                    model.add(MaxPool1D(pool_size=2,padding="valid"))

                    for _ in range(CNNLayer):
                        model.add(Conv1D(kernel_size = (2), padding = "valid", filters = filt))
                        model.add(Activation("elu"))
                        model.add(BatchNormalization())
                        model.add(MaxPool1D(pool_size=2,padding="valid"))

                    model.add(Flatten())

                    for _ in range(dense):
                        model.add(Dense(filt))
                        model.add(Activation("elu"))

                    model.add(Dense(self.nClasses, activation = "softmax" ))
                    model.compile(loss='categorical_crossentropy', optimizer = 'adam', metrics = ['accuracy'])

                    keyVals = {"model":model,
                                "name":nameOfModel}
                    models.append(keyVals)

        self.models = models

    def trainNetwork(self,epochs,batchSize):
        print("training network")
        self.model.fit(self.xTrain,self.yTrain,epochs = epochs, batch_size = batchSize,
        verbose = 1, validation_data = (self.xEval,self.yEval))

    def validateNetwork(self):
        accuracy = self.model.evaluate(self.xEval,self.yEval)
        print("The accurcy on the Eval Data was: " + str(accuracy))

    def testNetwork(self):
        accuracy = self.model.evaluate(self.xTest,self.yTest)
        print("The accurcy on the test Data was: " + str(accuracy))

dataModel1 = dataModel("labeledAccelerometer.json")
#dataModel1.buildNet2D()
dataModel1.buildNets1D()

wFile = open("results.text","w")

for keyVals in dataModel1.models:
    dataModel1.model = keyVals["model"]
    name = keyVals["name"]

    wFile.write("\n")
    wFile.write(name)

    dataModel1.trainNetwork(epochs=100,batchSize =128)
    sys.stdout = wFile

    dataModel1.validateNetwork()
    sys.stdout = sys.__stdout__

wFile.close()


text:
#! /usr/bin/env python

'''
Example vectorize usage.
'''

import math
import sys
from timeit import default_timer as time

import numpy as np

from numba import vectorize

def generate_input(n, dtype):
    A = np.array(np.random.sample(n), dtype=dtype)
    B = np.array(np.random.sample(n) + 10, dtype=dtype)
    C = np.array(np.random.sample(n), dtype=dtype)
    return A, B, C

def check_answer(ans, A, B, C):
    for d, a, b, c in zip(ans, A, B, C):
        gold = discriminant(a, b, c)
        assert np.allclose(d, gold), (d, gold)

def discriminant(a, b, c):
    '''a ufunc kernel to compute the discriminant of quadratic equation
    '''
    return math.sqrt(b ** 2 - 4 * a * c)

def main():

    N = 1e+8 // 2
    N = int(N)
    print('Data size', N)

    targets = ['cpu', 'parallel']

    # run just one target if is specified in the argument
    for t in targets:
        if t in sys.argv[1:]:
            targets = [t]
            break

    for target in targets:
        print('== Target', target)
        vect_discriminant = vectorize(['f4(f4, f4, f4)', 'f8(f8, f8, f8)'],
                                      target=target)(discriminant)

        A, B, C = generate_input(N, dtype=np.float32)
        D = np.empty(A.shape, dtype=A.dtype)

        ts = time()
        D = vect_discriminant(A, B, C)
        te = time()

        total_time = (te - ts)

        print('Execution time %.4f' % total_time)
        print('Throughput %.4f' % (N / total_time))

        if '-verify' in sys.argv[1:]:
            check_answer(D, A, B, C)

if __name__ == '__main__':
    main()

max_stars_repo_name: SPIRV/NUMBA

max_issues_repo_name: SPIRV/NUMBA

max_forks_repo_name: SPIRV/NUMBA

reason: import numpy,from numba


text:
#reading different model files

import numpy as np
from numpy import recfromtxt, genfromtxt
import pandas as pd
from astropy import units as u

import logging
# Adding logging support
logger = logging.getLogger(__name__)

from tardis.util import parse_quantity

class ConfigurationError(Exception):
    pass

def read_density_file(filename, filetype):
    """
    read different density file formats

    Parameters
    ----------

    filename: ~str
        filename or path of the density file

    filetype: ~str
        type of the density file

    Returns
    -------
    time_of_model: ~astropy.units.Quantity
        time at which the model is valid

    velocity: ~np.ndarray
        the array containing the velocities

    unscaled_mean_densities: ~np.ndarray
        the array containing the densities

    """
    file_parsers = {'artis': read_artis_density,
                    'simple_ascii': read_simple_ascii_density}

    (time_of_model, velocity,
     unscaled_mean_densities) = file_parsers[filetype](filename)
    v_inner = velocity[:-1]
    v_outer = velocity[1:]
    invalid_volume_mask = (v_outer - v_inner) <= 0
    if invalid_volume_mask.sum() > 0:
        message = "\n".join(["cell {0:d}: v_inner {1:s}, v_outer "
                             "{2:s}".format(i, v_inner_i, v_outer_i) for i,
                               v_inner_i, v_outer_i in
                             zip(np.arange(len(v_outer))[invalid_volume_mask],
                                 v_inner[invalid_volume_mask],
                                 v_outer[invalid_volume_mask])])
        raise ConfigurationError("Invalid volume of following cell(s):\n"
                                 "{:s}".format(message))

    return time_of_model, velocity, unscaled_mean_densities

def read_abundances_file(abundance_filename, abundance_filetype,
                         inner_boundary_index=None, outer_boundary_index=None):
    """
    read different density file formats

    Parameters
    ----------

    abundance_filename: ~str
        filename or path of the density file

    abundance_filetype: ~str
        type of the density file

    inner_boundary_index: int
        index of the inner shell, default None

    outer_boundary_index: int
        index of the outer shell, default None

    """

    file_parsers = {'simple_ascii': read_simple_ascii_abundances,
                    'artis': read_simple_ascii_abundances}

    index, abundances = file_parsers[abundance_filetype](abundance_filename)
    if outer_boundary_index is not None:
        outer_boundary_index_m1 = outer_boundary_index - 1
    else:
        outer_boundary_index_m1 = None
    index = index[inner_boundary_index:outer_boundary_index]
    abundances = abundances.ix[:, slice(inner_boundary_index, outer_boundary_index_m1)]
    abundances.columns = np.arange(len(abundances.columns))
    return index, abundances

def read_simple_ascii_density(fname):
    """
    Reading a density file of the following structure (example; lines starting with a hash will be ignored):
    The first density describes the mean density in the center of the model and is not used.
    5 s
    #index velocity [km/s] density [g/cm^3]
    0 1.1e4 1.6e8
    1 1.2e4 1.7e8

    Parameters
    ----------

    fname: str
        filename or path with filename

    Returns
    -------

    time_of_model: ~astropy.units.Quantity
        time at which the model is valid

    data: ~pandas.DataFrame
        data frame containing index, velocity (in km/s) and density
    """

    with open(fname) as fh:
        time_of_model_string = fh.readline().strip()
        time_of_model = parse_quantity(time_of_model_string)

    data = recfromtxt(fname, skip_header=1,
                      names=('index', 'velocity', 'density'),
                      dtype=(int, float, float))
    velocity = (data['velocity'] * u.km / u.s).to('cm/s')
    mean_density = (data['density'] * u.Unit('g/cm^3'))[1:]

    return time_of_model, velocity, mean_density

def read_artis_density(fname):
    """
    Reading a density file of the following structure (example; lines starting with a hash will be ignored):
    The first density describes the mean density in the center of the model and is not used.
    5
    #index velocity [km/s] log10(density) [log10(g/cm^3)]
    0 1.1e4 1.6e8
    1 1.2e4 1.7e8

    Parameters
    ----------

    fname: str
        filename or path with filename

    Returns
    -------

    time_of_model: ~astropy.units.Quantity
        time at which the model is valid

    data: ~pandas.DataFrame
        data frame containing index, velocity (in km/s) and density
    """

    with open(fname) as fh:
        for i, line in enumerate(file(fname)):
            if i == 0:
                no_of_shells = np.int64(line.strip())
            elif i == 1:
                time_of_model = u.Quantity(float(line.strip()), 'day').to('s')
            elif i == 2:
                break

    artis_model_columns = ['index', 'velocities', 'mean_densities_0', 'ni56_fraction', 'co56_fraction', 'fe52_fraction',
                           'cr48_fraction']
    artis_model = recfromtxt(fname, skip_header=2, usecols=(0, 1, 2, 4, 5, 6, 7), unpack=True,
                                dtype=[(item, np.float64) for item in artis_model_columns])

    velocity = u.Quantity(artis_model['velocities'], 'km/s').to('cm/s')
    mean_density = u.Quantity(10 ** artis_model['mean_densities_0'], 'g/cm^3')[1:]

    return time_of_model, velocity, mean_density

def read_simple_ascii_abundances(fname):
    """
    Reading an abundance file of the following structure (example; lines starting with hash will be ignored):
    The first line of abundances describe the abundances in the center of the model and are not used.
    #index element1, element2, ..., element30
    0 0.4 0.3, .. 0.2

    Parameters
    ----------

    fname: str
        filename or path with filename

    Returns
    -------

    index: ~np.ndarray
        containing the indices

    abundances: ~pandas.DataFrame
        data frame containing index, element1 - element30 and columns according to the shells
    """
    data = np.loadtxt(fname)

    index = data[1:,0].astype(int)
    abundances = pd.DataFrame(data[1:,1:].transpose(), index=np.arange(1, data.shape[1]))

    return index, abundances

max_stars_repo_name: laudb/tardis

max_issues_repo_name: laudb/tardis

max_forks_repo_name: laudb/tardis

reason: import numpy,from numpy,from astropy


text:
import tensorflow as tf
import numpy as np
from tensorflow import keras

xs = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
ys = np.array([1.0, 1.5, 2.0, 2.5, 3.0, 3.5])

model = keras.Sequential([keras.layers.Dense(units=1, input_shape=[1])])
model.compile(optimizer='sgd', loss='mean_squared_error')
model.fit(xs, ys, epochs=500)

print(model.predict([7.0]))

max_stars_repo_name: cimomo/tf-lib

max_issues_repo_name: cimomo/tf-lib

max_forks_repo_name: cimomo/tf-lib


text:
import re, os
import numpy as np

from .normalized_hermite import hermval_weighted, hermval

def read_hermite_basis(fname='spectral_basis.desc'):
    """
    Keyword Arguments:
    fname -- (default 'spectral_basis.desc')
    """
    pattern = 'H_([0-9]+),\s+fw_(.*)\s+H_([0-9]+),\s+fw_(.*)'
    f = open(fname,'r')

    elements = []
    for l in f.readlines():
        match = re.match(pattern, l)
        if match:
            k1 = int(match.group(1))
            w1 = float(match.group(2))
            k2 = int(match.group(3))
            w2 = float(match.group(4))
            elements.append( {'k1' : k1, 'w1' : w1, 'k2' : k2, 'w2' : w2} )

    return elements

def hermite_evaluatez(C, Z, N, weighted=True):
    """
    Evaluate 2-d Hermite series with weight exp(-|z|^2/2)

    Keyword Arguments:
    C -- coefficients as 1d array
    Z -- evaluation points Z = X + 1j Y
    N -- max degree in one direction
    """

    X = np.real(Z)
    Y = np.imag(Z)

    degx = np.array([]).astype(np.int)
    degy = np.array([]).astype(np.int)

    for i in range(N+1):
        t = np.arange(i+1).astype(np.int)

        degx = np.hstack((degx,t))
        degy = np.hstack((degy,t[::-1]))

    #return degx, degy

    assert(len(degx) == ((N+1)*(N+2))/2)

    if weighted:
        Fx = hermval_weighted(X, np.diag(np.ones(N+1)))
        Fy = hermval_weighted(Y, np.diag(np.ones(N+1)))
    else:
        Fx = hermval(X, np.diag(np.ones(N+1)))
        Fy = hermval(Y, np.diag(np.ones(N+1)))

    return np.einsum(Fx[...,degx]*Fy[...,degy], [1, 2, 3], C, [3])


text:
import json
import glob, os
import dash
import datetime
import dash_core_components as dcc
import dash_html_components as html
from dash.dependencies import Input, Output
import plotly.express as px
from django_plotly_dash import DjangoDash
from django.apps import apps
import plotly.graph_objects as go
from Lab_Misc.General import *
from dbfread import DBF
import pandas as pd
from django.http import HttpResponse
from django.utils import timezone
import numpy as np
import time
from Exp_Main.models import SEL
from Exp_Sub.models import LSP
from plotly.subplots import make_subplots
from Lab_Misc import General

def conv(x):
    return x.replace(',', '.').encode()
def Gen_dash(dash_name):
    class Gen_fig():
        def load_data(self, target_id):
            try:
                entry = SEL.objects.get(id = target_id)
                self.entry = entry
                file = os.path.join( rel_path, entry.Link_XLSX)
                df = pd.read_excel(file, 'Tabelle1')
                new_vals = df[df>1]/1000000#correct for wrong format
                Curr_Dash = self.entry.Dash
                df.update(new_vals)
                self.data = df
                self.data["Time (min.)"] = Curr_Dash.Start_datetime_elli + pd.TimedeltaIndex(self.data["Time (min.)"], unit='m')
                self.data["Time (min.)"] = self.data["Time (min.)"].dt.tz_convert(timezone.get_current_timezone())
                return_str = 'The following data could be loaded: Ellisometry'
                self.sub_data, return_str_sub = self.get_subs()
                return_str += return_str_sub
                os.chdir(cwd)
                return True, return_str
            except:
                return False, 'No data found!'

        def has_sub(self):
            try:
                #Sub_Exp = self.entry.Sub_Exp.all()
                #print(Sub_Exp)
                #self.get_subs()
                return True
            except:
                return False

        def get_subs(self):
            Sub_Exps = self.entry.Sub_Exp.all()
            Sub_Exps_dic = {}
            return_str = ''
            for Sub_Exp in Sub_Exps:
                Device = Sub_Exp.Device
                model = apps.get_model('Exp_Sub', str(Device.Abbrev))
                Exp_in_Model = model.objects.get(id = Sub_Exp.id)
                if Device.Abbrev == 'MFL':
                    Gas = Exp_in_Model.Gas.first()
                    if Gas.Name == 'H2O':
                        MFL_H2O_data = self.get_sub_csv(Exp_in_Model)
                        MFL_H2O_data['Date_Time'] = pd.to_datetime(MFL_H2O_data['Date_Time'], format='%d.%m.%Y %H:%M:%S', errors="coerce")
                        MFL_H2O_data['Date_Time'] = MFL_H2O_data['Date_Time'].dt.tz_localize(timezone.get_current_timezone())
                        Sub_Exps_dic.update(MFL_H2O_data = MFL_H2O_data)
                        return_str += ', massflow of water stream'
                    if Gas.Name == 'N2':
                        MFL_N2_data = self.get_sub_csv(Exp_in_Model)
                        MFL_N2_data['Date_Time'] = pd.to_datetime(MFL_N2_data['Date_Time'], format='%d.%m.%Y %H:%M:%S', errors="coerce")
                        MFL_N2_data['Date_Time'] = MFL_N2_data['Date_Time'].dt.tz_localize(timezone.get_current_timezone())
                        Sub_Exps_dic.update(MFL_N2_data = MFL_N2_data)
                        return_str += ', massflow of nitrogen stream'
                if Device.Abbrev == 'HME':
                    for pos_env in Exp_in_Model.PossibleEnvironments:
                        if pos_env[0] == Exp_in_Model.Environments:
                            if pos_env[1] == 'Cell':
                                Humidity_data = self.get_sub_dbf(Exp_in_Model)
                                Humidity_data['UHRZEIT'] = pd.to_datetime(Humidity_data['DATUM'] + Humidity_data['UHRZEIT'], format='%d.%m.%Y    %H:%M:%S', errors="coerce")
                                Humidity_data['UHRZEIT'] = Humidity_data['UHRZEIT'].dt.tz_localize(timezone.get_current_timezone())
                                Sub_Exps_dic.update(HME_data = Humidity_data)
                                return_str += ', humidity measurements'
            return Sub_Exps_dic, return_str

        def get_sub_dbf(self, model):
            file = os.path.join( rel_path, model.Link)
            table = DBF(file, load=True)
            df = pd.DataFrame(iter(table))
            return df

        def get_sub_csv(self, model):
            file = os.path.join( rel_path, model.Link)
            #file_name = file[get_LastIndex(file, '\\')+1:get_LastIndex(file, '.')]
            df = pd.read_csv(file, sep=';', error_bad_lines=False, decimal = ',', parse_dates=[['Date', 'Time']])#skips bad lines
            return df

        def slice_data(self, data):
            DashTab = self.entry.Dash
            if isinstance(DashTab.CA_high_degree, float):
                slice_CA_high = (data['CA_L']<DashTab.CA_high_degree) & (data['CA_R']<DashTab.CA_high_degree)
                data = data[slice_CA_high]

            if isinstance(DashTab.CA_low_degree, float):
                slice_CA_low = (data['CA_L']>DashTab.CA_low_degree) & (data['CA_R']>DashTab.CA_low_degree)
                data = data[slice_CA_low]

            if isinstance(DashTab.BD_high_mm, float):
                slice_BD = (data['BI_left']<DashTab.BD_high_mm) & (data['BI_right']<DashTab.BD_high_mm)
                data = data[slice_BD]

            if isinstance(DashTab.BD_low_mm, float):
                slice_BD = (data['BI_left']>DashTab.BD_low_mm) & (data['BI_right']>DashTab.BD_low_mm)
                data = data[slice_BD]

            if isinstance(DashTab.Time_high_sec, float):
                slice_time = data['Age']<DashTab.Time_high_sec
                data = data[slice_time]

            if isinstance(DashTab.Time_low_sec, float):
                slice_time = data['Age']>DashTab.Time_low_sec
                data = data[slice_time]
            return data

        def CA_Time(self):
            fig = go.Figure()
            fig.add_trace(go.Scattergl(x=self.data['Time (min.)'], y=self.data['Thickness # 3'],
                        mode='markers',
                        name='CA left')
            )
            return fig

        def CA_BD(self):
            fig = go.Figure()
            fig.add_trace(go.Scattergl(x=self.data['BD'], y=self.data['CA_L'],
                        mode='markers',
                        name='CA left')
            )
            fig.add_trace(go.Scattergl(x=self.data['BD'], y=self.data['CA_R'],
                        mode='markers',
                        name='CA right')
            )
            fig.add_trace(go.Scattergl(x=self.data['BD'], y=self.data['CA_M'],
                        mode='markers',
                        name='CA average')
            )
            fig.update_layout(  xaxis_title='Base diameter [mm]',
                                yaxis_title='Contact angle [°]')
            return fig

        def Plot_3(self):
            Humidity_data = self.sub_data['HME_data']
            start_date = self.data['Time (min.)'][0]
            end_date = self.data['Time (min.)'][len(self.data)-1]
            date = start_date
            time_step = datetime.timedelta(seconds = 15)
            thickness = []
            humidities = []
            times  = []
            while(date<end_date):
                thicknes = np.mean(self.data['Thickness # 3'][(self.data['Time (min.)']<date+time_step)&(self.data['Time (min.)']>date)])
                thickness.append(thicknes)
                humidity = np.mean(Humidity_data['CHN1RH'][(Humidity_data['UHRZEIT']<date+time_step)&(Humidity_data['UHRZEIT']>date)])
                humidities.append(humidity)
                times.append(date)
                date = date + time_step
            self.cal_dat = pd.DataFrame(times, columns = ['times'])
            self.cal_dat["times"] = self.cal_dat["times"].dt.tz_convert('UTC')#time already shifted
            self.cal_dat ['thickness'] = thickness
            self.cal_dat ['humidities'] = humidities

    def get_figure(df, x_col, y_col, selectedpoints, selectedpoints_local):

        if selectedpoints_local and selectedpoints_local['range']:
            ranges = selectedpoints_local['range']
            selection_bounds = {'x0': ranges['x'][0], 'x1': ranges['x'][1],
                                'y0': ranges['y'][0], 'y1': ranges['y'][1]}
        else:
            selection_bounds = {'x0': np.min(df[x_col]), 'x1': np.max(df[x_col]),
                                'y0': np.min(df[y_col]), 'y1': np.max(df[y_col])}

        # set which points are selected with the `selectedpoints` property
        # and style those points with the `selected` and `unselected`
        # attribute. see
        # https://medium.com/@plotlygraphs/notes-from-the-latest-plotly-js-release-b035a5b43e21
        # for an explanation
        fig = px.scatter(df, x=df[x_col], y=df[y_col], text=df.index)

        fig.update_traces(selectedpoints=selectedpoints,
                        customdata=df.index,
                        mode='markers+text', marker={ 'color': 'rgba(0, 116, 217, 0.7)', 'size': 5 }, unselected={'marker': { 'color': 'rgba(200, 116, 0, 0.1)', 'size': 5 }, 'textfont': { 'color': 'rgba(0, 0, 0, 0)' }})

        fig.update_layout(margin={'l': 20, 'r': 0, 'b': 15, 't': 5}, dragmode='select', hovermode=False)

        fig.add_shape(dict({'type': 'rect',
                            'line': { 'width': 1, 'dash': 'dot', 'color': 'darkgrey' }},
                        **selection_bounds))
        return fig

    value = 'temp'

    app = DjangoDash(name=dash_name, id='target_id')
    cwd = os.getcwd()
    rel_path = General.get_BasePath()
    GenFig = Gen_fig()

    fig = {
                'data': [{
                    'y': [1]
                }],
                'layout': {
                    'height': 800
                }
            }

    app.layout = html.Div(children=[
        html.Div([
            html.Div(
                dcc.Graph(id='g1', config={'displayModeBar': True}),
                className='four columns',
                style={'width': '33%', 'display': 'inline-block'},
                ),
            html.Div(
                dcc.Graph(id='g2', config={'displayModeBar': True}),
                className='four columns',
                style={'width': '33%', 'display': 'inline-block'},
                ),
            html.Div(
                dcc.Graph(id='g3', config={'displayModeBar': True}),
                className='four columns',
                style={'width': '33%', 'display': 'inline-block'},
                )
            ], style={"display": "block"}, className='row'),
        dcc.Input(id='target_id', type='hidden', value='1'),
        html.Button('Load data', id='Load_Data'),
        dcc.Loading(
            id="loading",
            children=[html.Div([html.Div(id="loading-output")])],
            type="default",
        ),
    ])

    @app.callback(
        Output("loading-output", "children"),
        [Input('Load_Data', 'n_clicks'),
        Input('target_id', 'value'),])
    def update_output(n_clicks, target_id, *args,**kwargs):
        data_was_loaded, return_str = GenFig.load_data(target_id)
        if data_was_loaded:
            return_str += '\n Select the desired plot at the dropdown.'
        GenFig.Plot_3()
        return return_str
    # this callback defines 3 figures
    # as a function of the intersection of their 3 selections
    @app.callback(
        [Output('g1', 'figure'),
        Output('g2', 'figure'),
        Output('g3', 'figure')],
        [Input('target_id', 'value'),
        Input('g1', 'selectedData'),
        Input('g2', 'selectedData'),
        Input('g3', 'selectedData')]
        )
    def callback(target_id, selection1, selection2, selection3):
        no_data = True
        while no_data:
            try:
                GenFig.cal_dat
                no_data = False
            except:
                time.sleep(1)
        selectedpoints = GenFig.cal_dat.index
        for selected_data in [selection1, selection2, selection3]:
            if selected_data and selected_data['points']:
                selectedpoints = np.intersect1d(selectedpoints,
                    [p['customdata'] for p in selected_data['points']])

        return [get_figure(GenFig.cal_dat, "times", "thickness", selectedpoints, selection1),
                get_figure(GenFig.cal_dat, "times", "humidities", selectedpoints, selection2),
                get_figure(GenFig.cal_dat, "humidities", "thickness", selectedpoints, selection3)]


text:
import numpy as np

# arr1 = np.zeros((10,0)) Creates Empty Array if either dimension is 0
arr1 = np.zeros((10,))
arr2 = np.ones((10,))
arr3 = np.full((5,2),41)
arr4 = np.empty((10)) # (10,) == (10)

print(arr1)
print(arr2)
print(arr3)
print(arr4) # Empty fills with 1 => empty == ones

print(np.linspace(1,16))
print(arr1.reshape((5,2)))
print(arr1)

print(arr3.flatten())