text:
# -*- coding: utf-8 -*-
"""ATMS597_Project3_GroupG_anomalies.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1J8N2kKtayGkRRvv-0YCrpNJQ-y1luNzC

# Install and import packages
"""

# installs python packages needed
!pip install pydap
!pip install netcdf4
!apt-get -qq install libproj-dev proj-data proj-bin libgeos-dev
!pip install --upgrade --force-reinstall shapely --no-binary shapely
!pip install cython
!pip install cartopy

# Commented out IPython magic to ensure Python compatibility.
# import necessary packages
# %pylab inline
import pandas as pd
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.util as cutil

# Commented out IPython magic to ensure Python compatibility.
# the folder content the extreme precipitation day data
# %cd /content/drive/My Drive/ATMS597 Weather Climate Data Analysis/Module 3/Project 3/

"""# Calculate and plot the long term mean"""

# loading the long term mean data
uwnd = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/uwnd.mon.1981-2010.ltm.nc', engine='netcdf4')
vwnd = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/vwnd.mon.1981-2010.ltm.nc', engine='netcdf4')
hgt = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/hgt.mon.1981-2010.ltm.nc', engine='netcdf4')
air = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/air.mon.1981-2010.ltm.nc', engine='netcdf4')
shum = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/shum.mon.1981-2010.ltm.nc', engine='netcdf4')
skt = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/skt.sfc.mon.1981-2010.ltm.nc', engine='netcdf4')
uwnd_sig995 = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface/uwnd.sig995.mon.1981-2010.ltm.nc', engine='netcdf4')
vwnd_sig995 = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface/vwnd.sig995.mon.1981-2010.ltm.nc', engine='netcdf4')
pr_wtr = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface/pr_wtr.eatm.mon.1981-2010.ltm.nc', engine='netcdf4')
prate = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/prate.sfc.mon.ltm.nc', engine='netcdf4')
omega = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/omega.mon.1981-2010.ltm.nc', engine='netcdf4')

xr.Dataset.info(uwnd)

skt.isel(time=[9, 10, 11])

# Figure 1. 250 hPa wind vectors and wind speed
# (1981-2010 OND long term mean)

x = uwnd.lon
y = uwnd.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Caculate wind speed
u = uwnd['uwnd'].sel(level=250.).isel(time=[9, 10, 11]).mean(dim='time')
v = vwnd['vwnd'].sel(level=250.).isel(time=[9, 10, 11]).mean(dim='time')
wspd_uv = sqrt(u**2+v**2)

# Plot wind speed as filled contours
data = wspd_uv
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('250 hPa wind speed (m/s)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
quiver_strides = 4  # change the density of quiver plot
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              uwnd['uwnd'].sel(level=250.).isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides],
              vwnd['vwnd'].sel(level=250.).isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 20, r'$20\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('250 hPa wind vectors and wind speed (1981-2010 OND long term mean)', fontsize=14)

plt.show()
fig.savefig("figure1.png", dpi=300)

# Figure 2. 500 hPa wind vectors and geopotential height
# (1981-2010 OND long term mean)

x = hgt.lon
y = hgt.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot geopotential height as filled contours
data = hgt['hgt'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time')
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('500 hPa geopotential height (m)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              uwnd['uwnd'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides],
              vwnd['vwnd'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 10, r'$10\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('500 hPa wind vectors and geopotential height (1981-2010 OND long term mean)', fontsize=14)

plt.show()
fig.savefig("figure2.png", dpi=300)

# Figure 3. 850 hPa temperature, specific humidity, and winds
# (1981-2010 OND long term mean)

x = air.lon
y = air.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot temperature as filled contours
data = air['air'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time')
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.005, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('850 hPa temperature (degC)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot specific humidity as contours
data = shum['shum'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time')
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
# Option 1: plot colorbar for contour
cs = ax.contour(x_cyc, y, data_cyc, transform=ccrs.PlateCarree(), linewidths=1.5, cmap='Reds')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('850 hPa specific humidity (g/kg)', size=12, rotation=-90, verticalalignment='bottom')
# Option 2: label the number on the contour
# cs = ax.contour(x, y, data, transform=ccrs.PlateCarree(), linewidths=1.5, cmap='Reds');
# ax.clabel(cs, inline=1, fmt='%2.1f', fontsize=14)

# Plot wind barbs, but not all of them
quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              uwnd['uwnd'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides],
              vwnd['vwnd'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.63, 0.83, 5, r'$5\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('850 hPa temperature, specific humidity, and winds (1981-2010 OND long term mean)', fontsize=14)

plt.show()
fig.savefig("figure3.png", dpi=300)

# Figure 4. skin temperature, and surface wind vectors (sig995 level)
# (1981-2010 OND long term mean)

# For skin temperature, lat: 94, lon: 192
x = skt.lon
y = skt.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot skin temperature as filled contours
data = skt['skt'].isel(time=[9, 10, 11]).mean(dim='time')
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('skin temperature (degC)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
# For uwnd_sig995, lat: 73, lon: 144
x = uwnd_sig995.lon
y = uwnd_sig995.lat
quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              uwnd_sig995['uwnd'].isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides],
              vwnd_sig995['vwnd'].isel(time=[9, 10, 11]).mean(dim='time').values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 5, r'$5\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('skin temperature, and surface wind vectors (sig995 level) (1981-2010 OND long term mean)', fontsize=14)
ax.set_ylim(-88.542, 88.542)

plt.show()
fig.savefig("figure4.png", dpi=300)

# Figure 5. total atmospheric column water vapor
# (1981-2010 OND long term mean)

x = pr_wtr.lon
y = pr_wtr.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
data = pr_wtr['pr_wtr'].isel(time=[9, 10, 11]).mean(dim='time')
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('total atmospheric column water vapor (kg/$m^2$)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('total atmospheric column water vapor (1981-2010 OND long term mean)', fontsize=14)

plt.show()
fig.savefig("figure5.png", dpi=300)

# Figure 6. precipitation rate
# (1981-2010 OND long term mean)

x = prate.lon
y = prate.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
data = prate['prate'].isel(time=[9, 10, 11]).mean(dim='time') * 86400  # kg/m^2/s to mm/day
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(0, 17), transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('precipitation rate (mm/day)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('precipitation rate (1981-2010 OND long term mean)', fontsize=14)
ax.set_ylim(-88.542, 88.542)

plt.show()
fig.savefig("figure6.png", dpi=300)

# Figure 7. 500 hPa omega
# (1981-2010 OND long term mean)

x = omega.lon
y = omega.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
data = omega['omega'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time')
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('omega (Pa/s)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('500 hPa omega (1981-2010 OND long term mean)', fontsize=14)
ax.set_ylim(-88.542, 88.542)

plt.show()
fig.savefig("figure7.png", dpi=300)

"""# Calculate extreme precipitation day composite"""

dates_95 = xr.open_dataset('datesexceeded.nc')  # extreme precipitation date

dates = pd.to_datetime(dates_95.time.values)
dates

iyr = 1996
dates_yr = dates[where(dates.year == iyr)]
dates_yr

# loading the daily data for extreme precip dates
years = pd.date_range(start='1996-01-01', end='2019-12-01', freq='AS')
varname = ['uwnd', 'vwnd', 'hgt', 'air', 'shum', 'skt', 'uwnd_sig995', 'vwnd_sig995', 'pr_wtr', 'prate', 'omega']

data_list = {}
for ivar in varname:
    data_list[ivar] = []

for iyr in years.year:
    print('working on '+str(iyr))
    # extreme precip dates for this year
    dates_yr = dates[where(dates.year == iyr)]
    print('how many extreme dates?', dates_yr.size)
    # collect data except for wspd, which is not available in daily dataset
    data_sel = {}
    data_sel['uwnd'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/pressure/uwnd.'+str(iyr)+'.nc', engine='netcdf4').sel(time=dates_yr)
    data_sel['vwnd'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/pressure/vwnd.'+str(iyr)+'.nc', engine='netcdf4').sel(time=dates_yr)
    data_sel['hgt'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/pressure/hgt.'+str(iyr)+'.nc', engine='netcdf4').sel(level=500, time=dates_yr)
    data_sel['air'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/pressure/air.'+str(iyr)+'.nc', engine='netcdf4').sel(level=850, time=dates_yr)
    data_sel['shum'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/pressure/shum.'+str(iyr)+'.nc', engine='netcdf4').sel(level=850, time=dates_yr)
    data_sel['skt'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/surface_gauss/skt.sfc.gauss.'+str(iyr)+'.nc', engine='netcdf4').sel(time=dates_yr)
    data_sel['uwnd_sig995'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/surface/uwnd.sig995.'+str(iyr)+'.nc', engine='netcdf4').sel(time=dates_yr)
    data_sel['vwnd_sig995'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/surface/vwnd.sig995.'+str(iyr)+'.nc', engine='netcdf4').sel(time=dates_yr)
    data_sel['pr_wtr'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/surface/pr_wtr.eatm.'+str(iyr)+'.nc', engine='netcdf4').sel(time=dates_yr)
    data_sel['prate'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/surface_gauss/prate.sfc.gauss.'+str(iyr)+'.nc', engine='netcdf4').sel(time=dates_yr)
    data_sel['omega'] = xr.open_dataset('https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.dailyavgs/pressure/omega.'+str(iyr)+'.nc', engine='netcdf4').sel(level=500, time=dates_yr)
    for ivar in varname:
        data_list[ivar].append(data_sel[ivar])

# concatenate the data and calculate the time average
print('concatenate the data and calculate the time average')
extremes_tavg = {}
for ivar in varname:
    extremes_tavg[ivar] = xr.concat(data_list[ivar], dim='time').mean(dim='time')

# save the extreme precipitation day composite
for ivar in varname:
    filename = 'extremes_tavg_' + ivar
    extremes_tavg[ivar].to_netcdf('./' + filename)

# # we can use the saved data instead of generating the data again
# extremes_tavg = {}
# for ivar in varname:
#     extremes_tavg[ivar] = xr.open_dataset('./extremes_tavg_'+ivar)

extremes_tavg['air']

"""# Plot extreme precipitation day composite"""

# Figure 8. 250 hPa wind vectors and wind speed
# (extreme precipitation day composite)

x = extremes_tavg['uwnd'].lon
y = extremes_tavg['uwnd'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Caculate wind speed
u = extremes_tavg['uwnd'].sel(level=250.)
v = extremes_tavg['vwnd'].sel(level=250.)
wspd_uv = sqrt(u['uwnd']**2 + v['vwnd']**2)

# Plot wind speed as filled contours
data = wspd_uv
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('250 hPa wind speed (m/s)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
quiver_strides = 4  # change the density of quiver plot
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              u['uwnd'][::quiver_strides, ::quiver_strides],
              v['vwnd'][::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 20, r'$20\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('250 hPa wind vectors and wind speed (extreme precipitation day composite)', fontsize=14)

plt.show()
fig.savefig("figure8.png", dpi=300)

# Figure 9. 500 hPa wind vectors and geopotential height
# (extreme precipitation day composite)

x = extremes_tavg['hgt'].lon
y = extremes_tavg['hgt'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot geopotential height as filled contours
data = extremes_tavg['hgt']['hgt']
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('500 hPa geopotential height (m)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              extremes_tavg['uwnd'].sel(level=500.)['uwnd'][::quiver_strides, ::quiver_strides],
              extremes_tavg['vwnd'].sel(level=500.)['vwnd'][::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 10, r'$10\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('500 hPa wind vectors and geopotential height (extreme precipitation day composite)', fontsize=14)

plt.show()
fig.savefig("figure9.png", dpi=300)

# Figure 10. 850 hPa temperature, specific humidity, and winds
# (extreme precipitation day composite)

x = extremes_tavg['air'].lon
y = extremes_tavg['air'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot temperature as filled contours
data = extremes_tavg['air']['air'] - 273.15  # K to degC
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.005, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('850 hPa temperature (degC)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot specific humidity as contours
data = extremes_tavg['shum']['shum'] * 1000  # kg/kg to g/kg
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
# Option 1: plot colorbar for contour
cs = ax.contour(x_cyc, y, data_cyc, transform=ccrs.PlateCarree(), linewidths=1.5, cmap='Reds')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('850 hPa specific humidity (g/kg)', size=12, rotation=-90, verticalalignment='bottom')
# Option 2: label the number on the contour
# cs = ax.contour(x, y, data, transform=ccrs.PlateCarree(), linewidths=1.5, cmap='Reds')
# ax.clabel(cs, inline=1, fmt='%2.1f', fontsize=14)

# Plot wind barbs, but not all of them
quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              extremes_tavg['uwnd'].sel(level=850.)['uwnd'][::quiver_strides, ::quiver_strides],
              extremes_tavg['vwnd'].sel(level=850.)['vwnd'][::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.63, 0.83, 5, r'$5\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('850 hPa temperature, specific humidity, and winds (extreme precipitation day composite)', fontsize=14)

plt.show()
fig.savefig("figure10.png", dpi=300)

# Figure 11. skin temperature, and surface wind vectors (sig995 level)
# (extreme precipitation day composite)

# For skin temperature, lat: 94, lon: 192
x = extremes_tavg['skt'].lon
y = extremes_tavg['skt'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot skin temperature as filled contours
data = extremes_tavg['skt']['skt'] - 273.15  # K to degC
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('skin temperature (degC)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
# For uwnd_sig995, lat: 73, lon: 144
x = extremes_tavg['uwnd_sig995'].lon
y = extremes_tavg['uwnd_sig995'].lat
quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              extremes_tavg['uwnd_sig995']['uwnd'][::quiver_strides, ::quiver_strides],
              extremes_tavg['vwnd_sig995']['vwnd'][::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 5, r'$5\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('skin temperature, and surface wind vectors (sig995 level) (extreme precipitation day composite)', fontsize=14)
ax.set_ylim(-88.542, 88.542)

plt.show()
fig.savefig("figure11.png", dpi=300)

# Figure 12. total atmospheric column water vapor
# (extreme precipitation day composite)

x = extremes_tavg['pr_wtr'].lon
y = extremes_tavg['pr_wtr'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
data = extremes_tavg['pr_wtr']['pr_wtr']
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('total atmospheric column water vapor (kg/$m^2$)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('total atmospheric column water vapor (extreme precipitation day composite)', fontsize=14)

plt.show()
fig.savefig("figure12.png", dpi=300)

# Figure 13. precipitation rate
# (extreme precipitation day composite)

x = extremes_tavg['prate'].lon
y = extremes_tavg['prate'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
data = extremes_tavg['prate']['prate'] * 86400  # kg/m^2/s to mm/day
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(0, 17), transform=ccrs.PlateCarree())
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('precipitation rate (mm/day)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('precipitation rate (extreme precipitation day composite)', fontsize=14)
ax.set_ylim(-88.542, 88.542)

plt.show()
fig.savefig("figure13.png", dpi=300)

# Figure 14. 500 hPa omega
# (extreme precipitation day composite)

x = extremes_tavg['omega'].lon
y = extremes_tavg['omega'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
data = extremes_tavg['omega']['omega']
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('omega (Pa/s)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('500 hPa omega (extreme precipitation day composite)', fontsize=14)

plt.show()
fig.savefig("figure14.png", dpi=300)

"""# Plot extreme precipitation day anomaly"""

# Figure 15. 250 hPa wind vectors and wind speed anomalies
# (compared with 1981-2010 OND long term mean)

x = uwnd.lon
y = uwnd.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot wind speed as filled contours
u = uwnd['uwnd'].sel(level=250.).isel(time=[9, 10, 11]).mean(dim='time')
v = vwnd['vwnd'].sel(level=250.).isel(time=[9, 10, 11]).mean(dim='time')
wspd_uv = sqrt(u**2+v**2)
u_ext = extremes_tavg['uwnd'].sel(level=250.)
v_ext = extremes_tavg['vwnd'].sel(level=250.)
wspd_uv_ext = sqrt(u_ext['uwnd']**2 + v_ext['vwnd']**2)
u_diff = u_ext['uwnd'] - u
v_diff = v_ext['vwnd'] - v
wspd_uv_diff = wspd_uv_ext - wspd_uv

data = wspd_uv_diff
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(-4, 5, 1), transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('250 hPa wind speed (m/s)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              u_diff.values[::quiver_strides, ::quiver_strides],
              v_diff.values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 2, r'$2\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('250 hPa wind vectors and wind speed anomalies', fontsize=14)

plt.show()
fig.savefig("figure15.png", dpi=300)

# Figure 16. 500 hPa wind vectors and geopotential height anomalies
# (compared with 1981-2010 OND long term mean)

x = hgt.lon
y = hgt.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot geopotential height as filled contours
hgt_ltm = hgt['hgt'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time')
hgt_ext = extremes_tavg['hgt']
hgt_diff = hgt_ext['hgt'] - hgt_ltm

data = hgt_diff
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(-50, 51, 10), transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('500 hPa geopotential height (m)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
u = uwnd['uwnd'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time')
v = vwnd['vwnd'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time')
u_ext = extremes_tavg['uwnd'].sel(level=500.)
v_ext = extremes_tavg['vwnd'].sel(level=500.)
u_diff = u_ext['uwnd'] - u
v_diff = v_ext['vwnd'] - v

quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              u_diff.values[::quiver_strides, ::quiver_strides],
              v_diff.values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 1, r'$1\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('500 hPa wind vectors and geopotential height anomalies', fontsize=14)

plt.show()
fig.savefig("figure16.png", dpi=300)

# Figure 17. 850 hPa temperature, specific humidity, and winds anomalies
# (compared with 1981-2010 OND long term mean)

x = air.lon
y = air.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot temperature as filled contours
air_ltm = air['air'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time')
air_ext = extremes_tavg['air'] - 273.15
air_diff = air_ext['air'] - air_ltm

data = air_diff
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(-4, 5, 1), transform=ccrs.PlateCarree(), cmap='coolwarm', alpha=0.7)
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.005, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('850 hPa temperature (degC)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot specific humidity as contours
shum_ltm = shum['shum'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time')
shum_ext = extremes_tavg['shum'] * 1000.
shum_diff = shum_ext['shum'] - shum_ltm

data = shum_diff
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)
# Option 1: plot colorbar for contour
cs = ax.contour(x_cyc, y, data_cyc, np.arange(-1.6, 1.7, 0.4), transform=ccrs.PlateCarree(), linewidths=1.5, cmap='seismic', alpha=1)
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('850 hPa specific humidity (g/kg)', size=12, rotation=-90, verticalalignment='bottom')
# Option 2: label the number on the contour
# cs = ax.contour(x, y, data, transform=ccrs.PlateCarree(), linewidths=1.5, cmap='Reds')
# ax.clabel(cs, inline=1, fmt='%2.1f', fontsize=14)

# Plot wind barbs, but not all of them
u850_ltm = uwnd['uwnd'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time')
v850_ltm = vwnd['vwnd'].sel(level=850.).isel(time=[9, 10, 11]).mean(dim='time')
u850_ext = extremes_tavg['uwnd'].sel(level=850.)
v850_ext = extremes_tavg['vwnd'].sel(level=850.)
u_diff = u850_ext['uwnd'] - u850_ltm
v_diff = v850_ext['vwnd'] - v850_ltm

quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              u_diff.values[::quiver_strides, ::quiver_strides],
              v_diff.values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.75, 0.83, 1, r'$1\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('850 hPa temperature, specific humidity, and winds anomalies', fontsize=14)

plt.show()
fig.savefig("figure17.png", dpi=300)

# Figure 18. skin temperature, and surface wind vectors (sig995 level) anomalies
# (compared with 1981-2010 OND long term mean)

# For skin temperature, lat: 94, lon: 192
x = skt.lon
y = skt.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot skin temperature as filled contours
skt_ltm = skt['skt'].isel(time=[9, 10, 11]).mean(dim='time')
skt_ext = extremes_tavg['skt'] - 273.15
skt_diff = skt_ext['skt'] - skt_ltm

data = skt_diff
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(-8, 9, 2), transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('skin temperature (degC)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()

# Plot wind barbs, but not all of them
# For uwnd_sig995, lat: 73, lon: 144
x = uwnd_sig995.lon
y = uwnd_sig995.lat
usfc_ltm = uwnd_sig995['uwnd'].isel(time=[9, 10, 11]).mean(dim='time')
vsfc_ltm = vwnd_sig995['vwnd'].isel(time=[9, 10, 11]).mean(dim='time')
usfc_ext = extremes_tavg['uwnd_sig995']
vsfc_ext = extremes_tavg['vwnd_sig995']
usfc_diff = usfc_ext['uwnd'] - usfc_ltm
vsfc_diff = vsfc_ext['vwnd'] - vsfc_ltm

quiver_strides = 4
Q = ax.quiver(x[::quiver_strides], y[::quiver_strides],
              usfc_diff.values[::quiver_strides, ::quiver_strides],
              vsfc_diff.values[::quiver_strides, ::quiver_strides])
qk = ax.quiverkey(Q, 0.8, 0.88, 1, r'$1\ \frac{m}{s}$', labelpos='E', coordinates='figure')
ax.set_title('skin temperature, and surface wind vectors (sig995 level) anomalies', fontsize=14)
ax.set_ylim(-88.542, 88.542)

plt.show()
fig.savefig("figure18.png", dpi=300)

# Figure 19. total atmospheric column water vapor anomalies
# (compared with 1981-2010 OND long term mean)

x = pr_wtr.lon
y = pr_wtr.lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
pr_wtr_ltm = pr_wtr['pr_wtr'].isel(time=[9, 10, 11]).mean(dim='time')
pr_wtr_ext = extremes_tavg['pr_wtr']
pr_wtr_diff = pr_wtr_ext['pr_wtr'] - pr_wtr_ltm

data = pr_wtr_diff
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(-6, 7, 1.5), transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('total atmospheric column water vapor (kg/$m^2$)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('total atmospheric column water vapor anomalies', fontsize=14)

plt.show()
fig.savefig("figure19.png", dpi=300)

# Figure 20. precipitation rate anomalies
# (compared with 1981-2010 OND long term mean)

x = extremes_tavg['prate'].lon
y = extremes_tavg['prate'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
prate_ltm = prate['prate'].isel(time=[9, 10, 11]).mean(dim='time')
prate_ext = extremes_tavg['prate']
prate_diff = prate_ext['prate'] - prate_ltm

data = prate_diff * 86400  # kg/m^2/s to mm/day
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(-4, 5, 1), transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('precipitation rate (mm/day)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('precipitation rate anomalies', fontsize=14)
ax.set_ylim(-88.542, 88.542)

plt.show()
fig.savefig("figure20.png", dpi=300)

# Figure 21. 500 hPa omega anomalies
# (compared with 1981-2010 OND long term mean)

x = extremes_tavg['omega'].lon
y = extremes_tavg['omega'].lat

# Create the matplotlib figure and axis
fig = plt.figure(figsize=(20, 8.5))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot total atmospheric column water vapor as filled contours
omega_ltm = omega['omega'].sel(level=500.).isel(time=[9, 10, 11]).mean(dim='time')
omega_ext = extremes_tavg['omega']
omega_diff = omega_ext['omega'] - omega_ltm

data = omega_diff
data_cyc, x_cyc = cutil.add_cyclic_point(data, coord=x)  # Add cyclic point
cs = ax.contourf(x_cyc, y, data_cyc, np.arange(-0.04, 0.05, 0.01), transform=ccrs.PlateCarree(), cmap='coolwarm')
cax, kw = matplotlib.colorbar.make_axes(ax, location='right', pad=0.05, shrink=0.8)
out = fig.colorbar(cs, cax=cax, extend='both')
label = out.set_label('omega (Pa/s)', size=12, rotation=-90, verticalalignment='bottom')
ax.set_global()
ax.coastlines()
ax.set_title('500 hPa omega anomalies', fontsize=14)

plt.show()
fig.savefig("figure21.png", dpi=300)


text:
"""
Plot Yield Plot Histograms
Written by: Dean Keithly
Written on: 11/28/2018
"""
from EXOSIMS.util.read_ipcluster_ensemble import gen_summary
import os
if not 'DISPLAY' in os.environ.keys(): #Check environment for keys
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
else:
    import matplotlib.pyplot as plt
import numpy as np
from cycler import cycler
import math
import datetime
import re

class yieldPlotHistograms(object):
    """Designed to plot yield plot histograms
    """
    _modtype = 'util'

    def __init__(self, args=None):
        """
        Args:
            args is not used
        """
        pass

    def singleRunPostProcessing(self, PPoutpath, folder):
        """Generates a single yield histogram for the run_type
        Args:
            PPoutpath (string) - output path to place data in
            folder (string) - full filepath to folder containing runs
        """
        res = gen_summary(folder)
        self.res = res
        self.dist_plot([res['detected']],PPoutpath=PPoutpath,folder=folder,legtext=[''])

    def multiRunPostProcessing(self, PPoutpath, folders):
        """Generates a yield histograms for the provided run_types
        Args:
            PPoutpath (string) - output path to place data in
            folders (string) - full filepaths to folders containing runs of each run_type
        """
        resList = list()
        for folder in folders:
            resList.append(gen_summary(folder)['detected'])
        self.resList = resList
        self.dist_plot(resList,PPoutpath=PPoutpath,folder=folder)

    def dist_plot(self,res,uniq = True,fig=None,lstyle='--',plotmeans=True,legtext=None,PPoutpath=os.getcwd(),folder=None):
        """
        Args:
            res (list) - list of gen_summary outputs [res1,res2,res3]
            uniq (boolean) - indicates whether to count unique detections (True) or all detections (False)
            fig (integer)
            lstyle (string)
            plotmeans (boolean)
            legtext (list) - list of strings containing legend labels for each element in res
            PPoutpath (string) -
            folder (string) - name of folder containing runs used to generate this plot (used for naming only)
        """
        #Set linewidth and color cycle
        plt.rc('axes',linewidth=2)
        plt.rc('lines',linewidth=2)
        plt.rc('axes',prop_cycle=(cycler('color',['red','purple','blue','black','darkorange','forestgreen'])))

        rcounts = []
        for el in res:
            if uniq:
                rcounts.append(np.array([np.unique(r).size for r in el]))
            else:
                rcounts.append(np.array([len(r) for r in el]))

        bins = range(np.min(np.hstack(rcounts).astype(int)),np.max(np.hstack(rcounts).astype(int))+2)
        bcents = np.diff(bins)/2. + bins[:-1]

        pdfs = []
        for j in range(len(res)):
            pdfs.append(np.histogram(rcounts[j],bins=bins,density=True)[0].astype(float))

        mx = 1.1*np.max(pdfs) #math.ceil(np.max(pdfs)*10)/10#np.round(np.max(pdfs),decimals=1)
        print(mx)

        syms = 'osp^v<>h'
        if fig is None:
            plt.figure()
        else:
            plt.figure(fig)
            plt.gca().set_prop_cycle(None)

        if legtext is None:
            legtext = [None]*len(res)

        for j in range(len(res)):
            leg = legtext[j]
            #c = plt.gca()._get_lines.prop_cycler.next()['color']#before 3.6
            c = plt.gca()._get_lines.prop_cycler.__next__()['color']#after 3.6
            if plotmeans:
                mn = np.mean(rcounts[j])
                plt.plot([mn]*2,[0,mx],'--',color=c)
                if leg is not None:
                    leg += ' ($\\mu = %2.2f$)'%mn
            plt.plot(bcents, pdfs[j], syms[np.mod(j,len(syms))]+lstyle,color=c,label=leg)

        plt.ylim([0,mx])
        if legtext[0] is not None:
            plt.legend()
        plt.xlabel('Unique Detections',weight='bold')
        plt.ylabel('Normalized Yield Frequency',weight='bold')

        date = datetime.datetime.now().strftime("%d_%b_%Y_%H_%M_%S_%f")
        #date = unicode(datetime.datetime.now())
        date = ''.join(c + '_' for c in re.split('-|:| ',date)[0:-1])#Removes seconds from date
        fname = 'yieldPlotHist_' + folder.split('/')[-1] + '_' + date
        plt.savefig(os.path.join(PPoutpath,fname+'.png'))
        plt.savefig(os.path.join(PPoutpath,fname+'.svg'))
        plt.savefig(os.path.join(PPoutpath,fname+'.eps'))
        plt.savefig(os.path.join(PPoutpath,fname+'.pdf'))

ext: py

max_stars_repo_name: jaysonfig/EXOSIMS

max_issues_repo_name: jaysonfig/EXOSIMS

max_forks_repo_name: jaysonfig/EXOSIMS

reason: import numpy


text:
# depricated see jupyter notebook : stats how does error carry forward

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

pi = np.pi
exp = np.exp
sqrt = np.sqrt

def get_cum_dist_from_unif(n_runs=1000,n_data=10000,delta=1.0):
    """a random variable x samples uniformly from (-delta/2,delta/2)
    What is the distribution of the random variable X = sum(x) over n such instances? I couldn't figure it out so this function is to do some experiments."""
    big_x = np.random.uniform(low=-delta/2,high=-delta/2, size=(n_runs,n_data))
    big_x.sum(axis=0)/sqrt(n_runs)
    return big_x

def plot_hist_from_dist(dist,n_runs,b=50):
    n,bins,_ = plt.hist(dist,b)
    plt.show(block=True)
    bin_means = bins[1:] / 2 + bins[:-1] / 2
    gauss = lambda x,sigma,amp:a*exp(-((x)/sigma)**2/2) / sqrt(2*pi*sigma**2) # assumes mean is zero, but not necessarily normalized
    popt,_ = curve_fit(gauss, bin_means, n/n_runs, p0=[1.,1.])
    sigma,amp = popt
    print("sigma = {} \namp {}".format(sigma,amp))
    # plt.show(block=True)
    x = np.linspace(bin_means[0],bin_means[-1],300)
    plt.plot(x, n_runs * gauss(x))
    plt.show(block=True)
    return

n_runs,n_data = 10000,10000
dist = get_cum_dist_from_unif(n_runs,n_data,delta=1.0)
print("got the dist")
plot_hist_from_dist(dist,n_runs,b=25)

max_stars_repo_name: dcxSt/pfb-mod

max_issues_repo_name: dcxSt/pfb-mod

max_forks_repo_name: dcxSt/pfb-mod

reason: import numpy,from scipy


text:
import numpy as np
from tdw.controller import Controller
from tdw.tdw_utils import TDWUtils
from tdw.output_data import OutputData, Meshes, Rigidbodies, Bounds

class PerlinNoiseTerrain(Controller):
    """
    Generate Perlin noise terrain and roll a ball down the terrain.
    """

    def run(self) -> None:
        self.start()

        # Download visual materials into the simulation.
        # Generate Perlin noise terrain.
        terrain_material = "stone_mountain_grey"
        ball_material = "plastic_vinyl_glossy_yellow"
        resp = self.communicate([self.get_add_material(material_name=terrain_material),
                                 self.get_add_material(material_name=ball_material),
                                 {"$type": "perlin_noise_terrain",
                                  "size": {"x": 24, "y": 24},
                                  "subdivisions": 1,
                                  "turbulence": 1.75,
                                  "origin": {"x": 0.5, "y": 0.5},
                                  "visual_material": terrain_material,
                                  "texture_scale": {"x": 4, "y": 2},
                                  "dynamic_friction": 0.25,
                                  "static_friction": 0.4,
                                  "bounciness": 0.2,
                                  "max_y": 10}])

        # Read the mesh data of the terrain.
        # Get a good spot to place the object: The highest position near the center of the terrain mesh.
        object_position: np.array = np.array([0, 0, 0])
        origin: np.array = np.array([0, 0])
        distance = 6
        for i in range(len(resp) - 1):
            r_id = OutputData.get_data_type_id(resp[i])
            if r_id == "mesh":
                meshes = Meshes(resp[i])
                vertices = meshes.get_vertices(0)
                for vertex in vertices:
                    if np.linalg.norm(np.array([vertex[0], vertex[2]]) - origin) > distance:
                        continue
                    if vertex[1] > object_position[1]:
                        object_position = vertex
        # Start the object above the vertex.
        object_position += 0.9

        # Add the ball. Set its visual material and physic material. Request output data.
        object_id = 0
        commands = [self.get_add_object(model_name="sphere",
                                        library="models_flex.json",
                                        position=TDWUtils.array_to_vector3(object_position),
                                        object_id=object_id),
                    {"$type": "set_visual_material",
                     "material_index": 0,
                     "material_name": ball_material,
                     "object_name": "sphere",
                     "id": object_id},
                    {"$type": "set_physic_material",
                     "dynamic_friction": 0.5,
                     "static_friction": 0.2,
                     "bounciness": 0.6,
                     "id": object_id},
                    {"$type": "send_bounds",
                     "frequency": "always"},
                    {"$type": "send_rigidbodies",
                     "frequency": "always"}]
        # Add a third-person camera.
        commands.extend(TDWUtils.create_avatar(position={"x": 12, "y": 6, "z": 12},
                                               look_at=object_id))
        resp = self.communicate(commands)

        # Run the simulation until either the ball stops moving or it falls off the edge of the terrain.
        done = False
        while not done:
            object_position: np.array = np.array([0, 0, 0])
            sleeping: bool = False
            for i in range(len(resp) - 1):
                r_id = OutputData.get_data_type_id(resp[i])
                if r_id == "boun":
                    boun = Bounds(resp[i])
                    for j in range(boun.get_num()):
                        if boun.get_id(j) == object_id:
                            object_position = np.array(boun.get_center(j))
                elif r_id == "rigi":
                    rigi = Rigidbodies(resp[i])
                    for j in range(rigi.get_num()):
                        if rigi.get_id(j) == object_id:
                            sleeping = rigi.get_sleeping(j)
            done = object_position[1] < -0.1 or sleeping
            # Look at the ball.
            resp = self.communicate({"$type": "look_at",
                                     "object_id": object_id})
        # End the simulation.
        self.communicate({"$type": "terminate"})

if __name__ == "__main__":
    c = PerlinNoiseTerrain(launch_build=False)
    c.run()

max_stars_repo_name: gragundier/tdw

max_issues_repo_name: NoelShin/tdw

max_forks_repo_name: NoelShin/tdw


text:
# -*- coding: utf-8 -*-
#
# Copyright (c) 2018, the cclib development team
#
# This file is part of cclib (http://cclib.github.io) and is distributed under
# the terms of the BSD 3-Clause License.

import unittest

import numpy

from cclib.bridge import cclib2pyquante

class PyquanteTest(unittest.TestCase):
    """Tests for the cclib2pyquante bridge in cclib."""

    def test_makepyquante(self):
        from PyQuante.hartree_fock import hf
        atomnos = numpy.array([1, 8, 1],"i")
        a = numpy.array([[-1, 1, 0], [0, 0, 0], [1, 1, 0]], "f")
        pyqmol = cclib2pyquante.makepyquante(a, atomnos)
        en, orbe, orbs = hf(pyqmol)
        ref = -75.824754
        assert abs(en - ref) < 1.0e-6

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

max_stars_repo_name: pstjohn/cclib

max_issues_repo_name: pstjohn/cclib

max_forks_repo_name: pstjohn/cclib


text:
import sys

sys.path.append('core.zip')
import json
import numpy as np
from core.utils.utils import NpEncoder
import requests
import argparse
from core.utils.udecorator import elapsed_time
import time

headers = {
    'User-Agent': 'Mozilla/5.0',
    'content-type': 'application/json',
    'charset': 'UTF-8',
}

# @elapsed_time
def ReindexApiTest(host, port, rc_id, num):
    url = 'http://{host}:{port}/{url}'.format(host=host, port=port, url='rc/reindex')
    d = 100
    nb = num
    np.random.seed(1234)
    xb = np.random.random((nb, d)).astype('float32')
    xb[:, 0] += np.arange(nb) / 1000.

    ids = ['u' + str(i) for i in range(nb)]

    data = json.dumps({
        'rcId': rc_id,
        'ids': ids,
        'vectors': json.loads(json.dumps(xb, cls=NpEncoder)),
    })

    s = requests.session()
    r = s.post(url, data=data, headers=headers)
    return r

if __name__ == '__main__':
    print('test begin...')

    parser = argparse.ArgumentParser()
    parser.add_argument('--host', type=str, default='localhost', required=False, help='')
    parser.add_argument('--port', type=str, default='8088', required=False, help='')
    parser.add_argument('--rc_id', type=str, default='101001101', required=False, help='')
    parser.add_argument('--num', type=int, default=1000, required=False, help='')

    params = parser.parse_args()

    req = ReindexApiTest(params.host, params.port, params.rc_id, params.num)

    print("""
            reindex
            status_code...%s
            content...%s

            """ % (req.status_code, req.content))


text:
import pytest
import numpy as np
from pyvibdmc import DistIt

udu = np.array([[-5.1637122053302426E-002, -0.64283289924786391, 1.5013908299665704],
                [-6.6986645862392238E-002, -1.0712440718151974, 0.63206523297920003],
                [0.68029777946852377, -1.0421568612924319, 1.9715519867047862],
                [-6.2517576960482130E-002, -0.96998368494931764, -1.2761479660404149],
                [-0.10689894504689866, -3.4286135083244317E-003, -1.2309713446779025],
                [-0.74897513808594363, -1.2354360967793472, -1.8873593576677006],
                [2.5686843929439932E-002, 1.6025162419831362, -0.16918517749731218],
                [2.0116888927609659E-002, 1.0626137857410543, 0.63588396614107578],
                [-0.59312555997327554, 2.3143071935596669, -9.0920564389553662E-003]])
uuu = np.array([[-2.5959991041286115E-002, 1.4991768513805555, -0.60084982145893295],
                [0.55551975361128614, 2.2546995583418465, -0.67743590894778982],
                [3.1332763535393070E-002, 1.2106512728103587, 0.32154786324596829],
                [-2.5969936853924985E-002, -1.2699370467222246, -0.99790377379813366],
                [3.1321499633255201E-002, -0.32685620892671469, -1.2092399775892559],
                [0.55548007897032570, -1.7140334262004899, -1.6139337867371595],
                [-2.5940453410975263E-002, -0.22924002701551754, 1.5987489246721986],
                [0.55556956122359069, -0.54066260106652508, 2.2913245852660338],
                [3.1338180636427793E-002, -0.88380042635289757, 0.88768175980027386]])
cds = np.array([udu, uuu])

def test_dist():
    """Unsorted distance matrix descriptor"""
    transformer = DistIt(zs=[8, 1, 1] * 3,
                         method='distance',
                         full_mat=False)
    dist_vec = transformer.run(cds)
    transformer = DistIt(zs=[8, 1, 1] * 3,
                         method='distance',
                         full_mat=True)
    dist_mat = transformer.run(cds)
    assert True

def test_spf_unsorted():
    """Sorted/Unsorted distance matrix descriptor"""
    from pyvibdmc import DistIt
    transformer = DistIt(zs=[8, 1, 1] * 3,
                         method='spf',
                         eq_xyz=cds[0],
                         full_mat=False,
                         )
    dist_mat_sorted = transformer.run(cds)
    assert True

def test_spf():
    """Sorted/Unsorted distance matrix descriptor"""
    from pyvibdmc import DistIt
    transformer = DistIt(zs=[8, 1, 1] * 3,
                         method='spf',
                         eq_xyz=cds[0],
                         full_mat=False,
                         sorted_atoms=[[0], [1, 2], [3], [4, 5], [6], [7, 8]]
                         )
    dist_mat_sorted = transformer.run(cds)
    assert True

def test_coulomb():
    """Sorted/Unsorted distance matrix descriptor"""
    from pyvibdmc import DistIt
    transformer = DistIt(zs=[8, 1, 1] * 3,
                         method='coulomb',
                         full_mat=True,
                         sorted_atoms=[[0], [1, 2], [3], [4, 5], [6], [7, 8]],
                         sorted_groups=[[0, 1, 2], [6, 7, 8]]
                         )
    dist_mat_sorted = transformer.run(np.tile(udu, (10, 1, 1)))
    assert True

def test_sorting():
    # Swap water 1 and water 2
    w1 = [0, 1, 2]
    w2 = [3, 4, 5]
    tot_atms = np.arange(len(udu))
    tot_atms_c = np.copy(tot_atms)
    tot_atms_c[w1] = w2
    tot_atms_c[w2] = w1
    w2w1 = np.copy(udu[tot_atms_c])

    # Swap H1 and H2 on water 1
    cop = np.copy(udu)
    cop[[1, 2]] = cop[[2, 1]]
    w1mod = cop

    # Swap H1 and H2 in water 1 in w2w1
    cop = np.copy(w2w1)
    cop[[1, 2]] = cop[[2, 1]]
    w2w1_w1mod = cop

    # All should have the same sorted spf
    srt_cds = np.array([udu, w2w1, w1mod, w2w1_w1mod])

    from pyvibdmc import DistIt
    transformer = DistIt(zs=[8, 1, 1] * 3,
                         method='spf',
                         eq_xyz=cds[0],
                         full_mat=False,
                         sorted_groups=[[0, 1, 2], [3, 4, 5], [6, 7, 8], ],
                         sorted_atoms=[[0], [1, 2], [3], [4, 5], [6], [7, 8]]
                         )
    dist_mat_sorted = transformer.run(srt_cds)
    assert np.count_nonzero(np.diff(dist_mat_sorted,axis=0)) == 0

def test_sorting_2():
    # Swap water 1 and water 2
    w1 = [0, 1, 2]
    w2 = [3, 4, 5]
    tot_atms = np.arange(len(uuu))
    tot_atms_c = np.copy(tot_atms)
    tot_atms_c[w1] = w2
    tot_atms_c[w2] = w1
    w2w1 = np.copy(uuu[tot_atms_c])

    # Swap H1 and H2 on water 1
    cop = np.copy(uuu)
    cop[[1, 2]] = cop[[2, 1]]
    w1mod = cop

    # Swap H1 and H2 in water 1 in w2w1
    cop = np.copy(w2w1)
    cop[[1, 2]] = cop[[2, 1]]
    w2w1_w1mod = cop

    # All should have the same sorted spf
    srt_cds = np.array([uuu, w2w1, w1mod, w2w1_w1mod])

    from pyvibdmc import DistIt
    transformer = DistIt(zs=[8, 1, 1] * 3,
                         method='spf',
                         eq_xyz=cds[0],
                         full_mat=False,
                         sorted_groups=[[0, 1, 2], [3, 4, 5], [6, 7, 8], ],
                         sorted_atoms=[[0], [1, 2], [3], [4, 5], [6], [7, 8]]
                         )
    dist_mat_sorted = transformer.run(srt_cds)
    assert np.count_nonzero(np.diff(dist_mat_sorted,axis=0)) == 0

max_stars_repo_name: McCoyGroup/pyvibdmc

max_issues_repo_name: McCoyGroup/pyvibdmc

max_forks_repo_name: McCoyGroup/pyvibdmc


text:
import argparse
import os
import time
import numpy as np
import tensorflow as tf
from trainer.model import get_data_generator, get_model_inputs, get_model_output
from trainer.pellet_list import PELLET_LIST
from package_ensemble import EntropyThresholdLayer

MODEL_TFLITE = "models/ab_recognition_i2a_and_amman_data.tflite"
MODEL_KERAS = "models/ab_recognition_i2a_and_amman_data.h5"

parser = argparse.ArgumentParser(
    description='Infer the label using the saved model')
parser.add_argument(
    'image', type=str,
    help='path to the image file showing one pellet with label')
parser.add_argument(
    '--img-size', type=int, default=64,
    help='square size to resize input images to in pixel, default=64')
parser.add_argument(
    '--model-type', choices=['keras', 'tflite'], default='tflite',
    help='keras or tflite')
parser.add_argument(
    '--model', default='models/ab_recognition_i2a_and_amman_data.tflite',
    help='path to model file')

def timeit(method):
    def timed(*args, **kw):
        ts = time.time()
        result = method(*args, **kw)
        te = time.time()
        print('%r  %2.2f ms' % (method.__name__, (te - ts) * 1000))
        return result
    return timed

@timeit
def infer_keras(image_path, img_size, path):
    inputs = np.array([get_model_inputs(image_path, img_size)])
    inputs_gen = get_data_generator().flow(inputs, shuffle=False)

    classifier = tf.keras.models.load_model(
        path, {'EntropyThresholdLayer': EntropyThresholdLayer})
    predictions = classifier.predict(inputs_gen)
    get_model_output(predictions, PELLET_LIST)

@timeit
def infer_tflite(image_path, img_size, path):
    inputs = np.array([get_model_inputs(image_path, img_size)], dtype=np.float32)
    inputs_gen = get_data_generator().standardize(inputs)

    interpreter = tf.lite.Interpreter(model_path=path)
    interpreter.allocate_tensors()

    input_details = interpreter.get_input_details()
    output_details = interpreter.get_output_details()

    interpreter.set_tensor(input_details[0]['index'], inputs_gen)
    interpreter.invoke()

    predictions = interpreter.get_tensor(output_details[0]['index'])
    get_model_output(predictions, PELLET_LIST)

if __name__ == '__main__':
    args = parser.parse_args()
    print("Using model", args.model)

    if args.model_type == 'keras':
        infer_keras(args.image, args.img_size, args.model)
    elif args.model_type == 'tflite':
        infer_tflite(args.image, args.img_size, args.model)

max_stars_repo_name: mpascucci/AST-image-processing

max_issues_repo_name: mpascucci/AST-image-processing

max_forks_repo_name: mpascucci/AST-image-processing


text:
import config
config.load()
"""
The above line MUST be imported first
Global parameters should not mutate at any time
Export path if necessay
export PYTHONPATH=~/dissertation/code; export LD_LIBRARY_PATH=~/venv/lib
"""
import nltk
import config
from nltk.tree import Tree
from tqdm import tqdm
from numba.core import types
from numba.typed import Dict, List
import numpy as np
from numba import njit, prange
from parsing.baseline import prune
from time import time
import os
import multiprocessing as mp
from parsing.contrained import constrained, get_parse_chart

def un_chomsky_normal_form(tree):
    """
    Modification of nltk.treetransforms.un_chomsky_normal_form
    """
    childChar = '|'
    unaryChar = '+'
    nodeList = [(tree, [])]
    while nodeList != []:
        node, parent = nodeList.pop()
        if parent and node not in parent:
            continue
        if isinstance(node, Tree):
            childIndex = node.label().find(childChar)
            if childIndex != -1 and node.label()[:childIndex] != parent.label():
                node.set_label(node.label()[:childIndex])
            elif childIndex != -1:
                nodeIndex = parent.index(node)
                parent.remove(parent[nodeIndex])
                parent.insert(nodeIndex, node[0])
                parent.insert(nodeIndex + 1, node[1])
                node = parent
            else:
                unaryIndex = node.label().find(unaryChar)
                if unaryIndex != -1:
                    newNode = Tree(
                        node.label()[unaryIndex + 1:], [i for i in node]
                    )
                    node.set_label(node.label()[:unaryIndex])
                    node[0:] = [newNode]
            for child in node:
                nodeList.append((child, node))

def transform_int2str(tree, sent, i=0):
    tree.set_label(config.nonterminal_map[int(tree.label())])
    if len(tree) == 1 and not isinstance(tree[0], Tree):
        tree[0] = sent[i]
        return i + 1
    else:
        for subtree in tree:
            i = transform_int2str(subtree, sent, i)
    return i

@njit
def recursive_build(parse_chart, score_chart, i, j, a=-1):
    if a == -1:
        assert(i == 0 and j == len(parse_chart) - 1)
        # print(parse_chart, i, j, a)
        # print(score_chart)
        best_score = -1
        for candidate, score in score_chart[i][j].items():
            if score > best_score:
                best_score = score
                a = candidate
    b, c, k = parse_chart[i][j][a]
    assert(i == j or (i <= k and k < j))
    root = '('+str(a) + ' '
    if i != j:
        root += recursive_build(parse_chart, score_chart, i, k, b)
        root += recursive_build(parse_chart, score_chart, k+1, j, c)
    else:
        root += 'w'
    return root + ')'

# @njit(nogil=True)
# @njit
def get_charts(terminals, r3_p, r1_p, pi_p, r3_lookupC, r1_lookup, prune_cutoff, r3_f, r1_f, pi_f):
    # Passing in the global parameters is necessary
    # t = time()
    #  = args
    # if terminals[0] == -1:
    #     return '()'
    # parse_chart, score_chart = prune(terminals, r3, r1, pi, r3_lookupC, r1_lookup, prune_cutoff)
    constrains = prune(terminals, r3_p, r1_p, pi_p, r3_lookupC, r1_lookup, prune_cutoff)
    if len(constrains[0][len(constrains) - 1]) == 0:
        return '()'

    # parse_chart, score_chart = get_parse_chart(constrains, len(constrains), r3_lookupC)
    # return recursive_build(parse_chart, score_chart, 0, len(parse_chart) - 1)

    marginal = constrained(terminals, r3_f, r1_f, pi_f, r3_lookupC, r1_lookup, constrains)
    if len(marginal[0][len(marginal) - 1]) == 0:
        # parse_chart, score_chart = get_parse_chart(constrains, len(marginal), r3_lookupC)
        return '()'
    else:
        parse_chart, score_chart = get_parse_chart(marginal, len(marginal), r3_lookupC)
    if len(parse_chart[0][len(parse_chart) - 1]) == 0:
        return '()'
    return recursive_build(parse_chart, score_chart, 0, len(parse_chart) - 1)

def process_wrapper(terminals):
    if terminals is None:
        return '()'
    if not config.numba_ready:
        from parsing import prepare_global_param
    return get_charts(List(terminals), config.rule3s_prune,
                        config.rule1s_prune,
                        config.pi_prune,
                        config.rule3s_lookupC,
                        config.rule1s_lookup,
                        config.prune_cutoff,
                        config.rule3s_full,
                        config.rule1s_full,
                        config.pi_full)
from nltk import pos_tag

def prepare_args(sent):
    # uncased = [w.lower() for w in sent]
    uncased = sent
    # print(uncased)
    rules = []
    for wordC, POS in pos_tag(uncased):
        word = wordC.lower()
        flag = False
        local = []
        if word in config.terminal_map.term2int:
            for rule in config.rule1s_lookup[config.terminal_map[word]]:
                rpos = config.nonterminal_map[rule.a]
                idx = rpos.rfind('+')
                rpos = rpos[idx+1:]
                if POS == rpos:
                    flag = True
                    local.append(rule)
            if not flag and POS in config.terminal_map.term2int:
                local.extend(config.rule1s_lookup[config.terminal_map[POS]])
            elif not flag:
                local.extend(config.rule1s_lookup[config.terminal_map[word]])
        elif POS in config.terminal_map.term2int:
            local.extend(config.rule1s_lookup[config.terminal_map[POS]])
        else:
            print(wordC, POS)
            return None
        rules.append(list(map(hash, local)))
    else:
        return rules

def parse_devset(dev_file):
    sents = []
    with open(dev_file, 'r') as f:
        for line in f:
            tree = Tree.fromstring(line)
            sents.append(tree.leaves())
    args = list(map(prepare_args, sents))
    # parses = []
    # t = tqdm(total=len(sents))
    cpu = os.cpu_count()
    with open(config.output_dir + 'baseline_parse.txt', 'w') as f:
        with mp.Pool(cpu-2) as pool:
            # result_futures = list(map(lambda arg: executor.submit(get_charts, *arg), args))
            # for future in concurrent.futures.as_completed(result_futures):
            for i, tree_str in enumerate(tqdm(pool.imap(process_wrapper, args, chunksize=len(sents)//cpu), total=len(sents))):
                if tree_str == '()':
                    f.write('()\n')
                else:
                    tree = Tree.fromstring(tree_str)
                    transform_int2str(tree, sents[i])
                    # un_chomsky_normal_form(tree)
                    tree.un_chomsky_normal_form()
                    parse = tree.pformat(margin=float('inf'))
                    f.write(parse + '\n')
                # parses.append(parse)
                # t.update()
    # return parses

def save(parses):
    new_lined = map(lambda s: s + '\n', parses)
    with open(config.output_dir + 'parse.txt', 'w') as f:
        f.writelines(new_lined)

if __name__ == '__main__':
    mp.set_start_method('fork')
    parse_devset(config.dev_file)

max_stars_repo_name: GavinPHR/Parser

max_issues_repo_name: GavinPHR/Parser

max_forks_repo_name: GavinPHR/Parser

reason: import numpy,from numba


text:
# Copyright 2018 Amazon.com, Inc. or its affiliates. All Rights Reserved.
#
#   Licensed under the Apache License, Version 2.0 (the "License").
#   You may not use this file except in compliance with the License.
#   A copy of the License is located at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   or in the "license" file accompanying this file. This file is distributed
#   on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
#   express or implied. See the License for the specific language governing
#   permissions and limitations under the License.
# ==============================================================================

import pytest

import numpy as np
import mxnet as mx
from scipy.stats import dirichlet as scipy_dirichlet

from mxfusion.components.variables.runtime_variable import add_sample_dimension, array_has_samples, get_num_samples
from mxfusion.components.variables import Variable
from mxfusion.components.distributions import Dirichlet
from mxfusion.util.testutils import numpy_array_reshape
from mxfusion.util.testutils import MockMXNetRandomGenerator

@pytest.mark.usefixtures("set_seed")
class TestDirichletDistribution(object):
    # scipy implementation of dirichlet throws an error if x_i does not sum to one and float32 is
    # not precise enough to have sum close enough to 1 after normalisation so using float64
    @pytest.mark.parametrize("dtype, a, a_is_samples, rv, rv_is_samples, num_samples", [
        (np.float64, np.random.rand(3,2), False, np.random.rand(5, 3, 2), True, 5),
        (np.float64, np.random.rand(10,3,2), True, np.random.rand(10, 3, 2), True, 10),
        (np.float64, np.random.rand(3,2), False, np.random.rand(3, 2), False, 5)
        ])
    def test_log_pdf_with_broadcast(self, dtype, a, a_is_samples, rv, rv_is_samples, num_samples):
        # Add sample dimension if varaible is not samples
        a_mx = mx.nd.array(a, dtype=dtype)
        if not a_is_samples:
            a_mx = add_sample_dimension(mx.nd, a_mx)
        a = a_mx.asnumpy()

        rv_mx = mx.nd.array(rv, dtype=dtype)
        if not rv_is_samples:
            rv_mx = add_sample_dimension(mx.nd, rv_mx)
        rv = rv_mx.asnumpy()

        is_samples_any = a_is_samples or rv_is_samples
        rv_shape = rv.shape[1:]

        n_dim = 1 + len(rv.shape) if is_samples_any and not rv_is_samples else len(rv.shape)
        a_np = np.broadcast_to(a, (num_samples, 3, 2))
        rv_np = numpy_array_reshape(rv, is_samples_any, n_dim)

        # Initialize rand_gen
        rand = np.random.rand(num_samples, *rv_shape)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))

        # Calculate correct Dirichlet logpdf
        r = []
        for s in range(len(rv_np)):
            a = []
            for i in range(len(rv_np[s])):
                a.append(scipy_dirichlet.logpdf(rv_np[s][i]/sum(rv_np[s][i]), a_np[s][i]))
            r.append(a)
        log_pdf_np = np.array(r)

        dirichlet = Dirichlet.define_variable(alpha=Variable(), shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor
        variables = {dirichlet.alpha.uuid: a_mx, dirichlet.random_variable.uuid: rv_mx}
        log_pdf_rt = dirichlet.log_pdf(F=mx.nd, variables=variables)

        assert np.issubdtype(log_pdf_rt.dtype, dtype)
        assert array_has_samples(mx.nd, log_pdf_rt) == is_samples_any
        if is_samples_any:
            assert get_num_samples(mx.nd, log_pdf_rt) == num_samples, (get_num_samples(mx.nd, log_pdf_rt), num_samples)
        assert np.allclose(log_pdf_np, log_pdf_rt.asnumpy())

    @pytest.mark.parametrize("dtype, a, a_is_samples, rv, rv_is_samples, num_samples", [
        (np.float64, np.random.rand(5, 3, 2), True, np.random.rand(5, 3, 2), True, 5),
        (np.float64, np.random.rand(10, 3, 2), True, np.random.rand(10, 3, 2), True, 10),
        ])
    def test_log_pdf_no_broadcast(self, dtype, a, a_is_samples,
                                  rv, rv_is_samples, num_samples):

        a_mx = mx.nd.array(a, dtype=dtype)
        if not a_is_samples:
            a_mx = add_sample_dimension(mx.nd, a_mx)
        a = a_mx.asnumpy()

        rv_mx = mx.nd.array(rv, dtype=dtype)
        if not rv_is_samples:
            rv_mx = add_sample_dimension(mx.nd, rv_mx)
        rv = rv_mx.asnumpy()

        is_samples_any = any([a_is_samples, rv_is_samples])
        rv_shape = rv.shape[1:]

        n_dim = 1 + len(rv.shape) if is_samples_any and not rv_is_samples else len(rv.shape)
        a_np = numpy_array_reshape(a, is_samples_any, n_dim)
        rv_np = numpy_array_reshape(rv, is_samples_any, n_dim)

        rand = np.random.rand(num_samples, *rv_shape)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))

        r = []
        for s in range(len(rv_np)):
            a = []
            for i in range(len(rv_np[s])):
                a.append(scipy_dirichlet.logpdf(rv_np[s][i]/sum(rv_np[s][i]), a_np[s][i]))
            r.append(a)
        log_pdf_np = np.array(r)

        dirichlet = Dirichlet.define_variable(alpha=Variable(), shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor
        variables = {dirichlet.alpha.uuid: a_mx, dirichlet.random_variable.uuid: rv_mx}
        log_pdf_rt = dirichlet.log_pdf(F=mx.nd, variables=variables)

        assert np.issubdtype(log_pdf_rt.dtype, dtype)
        assert array_has_samples(mx.nd, log_pdf_rt) == is_samples_any
        if is_samples_any:
            assert get_num_samples(mx.nd, log_pdf_rt) == num_samples, (get_num_samples(mx.nd, log_pdf_rt), num_samples)
        assert np.allclose(log_pdf_np, log_pdf_rt.asnumpy())

    @pytest.mark.parametrize("dtype, a, a_is_samples, rv_shape, num_samples", [
        (np.float64, np.random.rand(3,2), False, (3, 2), 5)
        ])
    def test_draw_samples_with_broadcast(self, dtype, a, a_is_samples, rv_shape, num_samples):
        a_mx = mx.nd.array(a, dtype=dtype)
        if not a_is_samples:
            a_mx = add_sample_dimension(mx.nd, a_mx)

        rand = np.random.gamma(shape=a, scale=np.ones(a.shape), size=(num_samples,)+rv_shape)
        draw_samples_np = rand / np.sum(rand)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))

        dirichlet = Dirichlet.define_variable(alpha=Variable(), shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor
        variables = {dirichlet.alpha.uuid: a_mx}
        draw_samples_rt = dirichlet.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert np.issubdtype(draw_samples_rt.dtype, dtype)
        assert draw_samples_rt.shape == (5,) + rv_shape
        assert np.allclose(draw_samples_np, draw_samples_rt.asnumpy())

    @pytest.mark.parametrize("dtype, a, a_is_samples, rv_shape, num_samples", [
        (np.float64, np.random.rand(5,3, 2), True, (3, 2), 5)
        ])
    def test_draw_samples_with_broadcast_no_numpy_verification(self, dtype, a, a_is_samples, rv_shape, num_samples):
        a_mx = mx.nd.array(a, dtype=dtype)
        if not a_is_samples:
            a_mx = add_sample_dimension(mx.nd, a_mx)

        dirichlet = Dirichlet.define_variable(alpha=Variable(), shape=rv_shape, dtype=dtype).factor
        variables = {dirichlet.alpha.uuid: a_mx}
        draw_samples_rt = dirichlet.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert np.issubdtype(draw_samples_rt.dtype, dtype)
        assert draw_samples_rt.shape == (5,) + rv_shape

    @pytest.mark.parametrize("dtype, a, a_is_samples, rv_shape, num_samples", [
        (np.float64, np.random.rand(5, 3, 2), True, (3, 2), 5)
        ])
    def test_draw_samples_no_broadcast(self, dtype, a, a_is_samples, rv_shape, num_samples):
        a_mx = mx.nd.array(a, dtype=dtype)
        if not a_is_samples:
            a_mx = add_sample_dimension(mx.nd, a_mx)

        rand = np.random.gamma(shape=a, scale=np.ones(a.shape), size=(num_samples,)+rv_shape)
        draw_samples_np = rand / np.sum(rand)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))

        dirichlet = Dirichlet.define_variable(alpha=Variable(), shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor
        variables = {dirichlet.alpha.uuid: a_mx}
        draw_samples_rt = dirichlet.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert np.issubdtype(draw_samples_rt.dtype, dtype)
        assert np.allclose(draw_samples_np, draw_samples_rt.asnumpy())

max_stars_repo_name: JeremiasKnoblauch/MXFusion

max_issues_repo_name: JeremiasKnoblauch/MXFusion

max_forks_repo_name: JeremiasKnoblauch/MXFusion


text:
import os
import random
import time
import torch
import torch.nn as nn
import numpy as np

from transformers import AdamW
from colbert.utils.runs import Run
from colbert.utils.amp import MixedPrecisionManager

from colbert.training.lazy_batcher import LazyBatcher
from colbert.training.eager_batcher import EagerBatcher
from colbert.parameters import DEVICE

from colbert.modeling.colbert import ColBERT
from colbert.utils.utils import print_message
from colbert.training.utils import print_progress, manage_checkpoints

def train(args):
    random.seed(12345)
    np.random.seed(12345)
    torch.manual_seed(12345)
    if args.distributed:
        torch.cuda.manual_seed_all(12345)

    if args.distributed:
        assert args.bsize % args.nranks == 0, (args.bsize, args.nranks)
        assert args.accumsteps == 1
        args.bsize = args.bsize // args.nranks

        print("Using args.bsize =", args.bsize, "(per process) and args.accumsteps =", args.accumsteps)

    if args.lazy:
        reader = LazyBatcher(args, (0 if args.rank == -1 else args.rank), args.nranks)
    else:
        reader = EagerBatcher(args, (0 if args.rank == -1 else args.rank), args.nranks)

    if args.rank not in [-1, 0]:
        torch.distributed.barrier()

    colbert = ColBERT.from_pretrained('microsoft/codebert-base',
                                      query_maxlen=args.query_maxlen,
                                      doc_maxlen=args.doc_maxlen,
                                      dim=args.dim,
                                      similarity_metric=args.similarity,
                                      mask_punctuation=args.mask_punctuation)

    if args.checkpoint is not None:
        assert args.resume_optimizer is False, "TODO: This would mean reload optimizer too."
        print_message(f"#> Starting from checkpoint {args.checkpoint} -- but NOT the optimizer!")

        checkpoint = torch.load(args.checkpoint, map_location='cpu')

        try:
            colbert.load_state_dict(checkpoint['model_state_dict'])
        except:
            print_message("[WARNING] Loading checkpoint with strict=False")
            colbert.load_state_dict(checkpoint['model_state_dict'], strict=False)

    if args.rank == 0:
        torch.distributed.barrier()

    colbert = colbert.to(DEVICE)
    colbert.train()

    if args.distributed:
        colbert = torch.nn.parallel.DistributedDataParallel(colbert, device_ids=[args.rank],
                                                            output_device=args.rank,
                                                            find_unused_parameters=True)

    optimizer = AdamW(filter(lambda p: p.requires_grad, colbert.parameters()), lr=args.lr, eps=1e-8)
    optimizer.zero_grad()

    amp = MixedPrecisionManager(args.amp)
    criterion = nn.CrossEntropyLoss()
    labels = torch.zeros(args.bsize, dtype=torch.long, device=DEVICE)

    start_time = time.time()
    train_loss = 0.0

    start_batch_idx = 0

    if args.resume:
        assert args.checkpoint is not None
        start_batch_idx = checkpoint['batch']

        reader.skip_to_batch(start_batch_idx, checkpoint['arguments']['bsize'])

    for batch_idx, BatchSteps in zip(range(start_batch_idx, args.maxsteps), reader):
        this_batch_loss = 0.0

        for queries, passages in BatchSteps:
            with amp.context():
                scores = colbert(queries, passages).view(2, -1).permute(1, 0)
                loss = criterion(scores, labels[:scores.size(0)])
                loss = loss / args.accumsteps

            if args.rank < 1:
                print_progress(scores)

            amp.backward(loss)

            train_loss += loss.item()
            this_batch_loss += loss.item()

        amp.step(colbert, optimizer)

        if args.rank < 1:
            avg_loss = train_loss / (batch_idx+1)

            num_examples_seen = (batch_idx - start_batch_idx) * args.bsize * args.nranks
            elapsed = float(time.time() - start_time)

            log_to_mlflow = (batch_idx % 20 == 0)
            Run.log_metric('train/avg_loss', avg_loss, step=batch_idx, log_to_mlflow=log_to_mlflow)
            Run.log_metric('train/batch_loss', this_batch_loss, step=batch_idx, log_to_mlflow=log_to_mlflow)
            Run.log_metric('train/examples', num_examples_seen, step=batch_idx, log_to_mlflow=log_to_mlflow)
            Run.log_metric('train/throughput', num_examples_seen / elapsed, step=batch_idx, log_to_mlflow=log_to_mlflow)

            print_message(batch_idx, avg_loss)
            manage_checkpoints(args, colbert, optimizer, batch_idx+1)

max_stars_repo_name: hueck/ColBERT

max_issues_repo_name: hueck/ColBERT

max_forks_repo_name: hueck/ColBERT


text:
# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

NbrOfNodes = 35
key200 = '  TIME: GANDRA  STEP:   180.000  FRAME:  1.000'
#--------------------------------------------------------------------------
#       File for gain parameter 01
#--------------------------------------------------------------------------
file_g01 = open('validation_005.rsn', 'r')
gain01 = file_g01.readlines()
g01 = pd.Series(gain01)
g01 = g01.replace(r'\n','', regex=True)
g01 = g01.replace(r'\r\n','', regex=True)
g01 = g01.replace(r'\r','', regex=True)
index_Time_g01 = g01[g01.str.contains('TIME', case=False, regex=False)]
index_TimeValues_g01 = index_Time_g01.index.values
#--------------------------------------------------------------------------
G01 = {}
for idx in index_Time_g01.index.values:
    index_start = idx + 1
    index_end = index_start + NbrOfNodes
    tmp_df = g01[index_start:index_end].str.strip()
    tmp_df = tmp_df.str.split('  ',expand=True)
    np.array(tmp_df.values, dtype=float)
    G01[g01[idx]]=np.array(tmp_df.values, dtype=float)

Data_g01 =  np.array([], dtype=np.int64)
Data_g01.shape = (-1, 7)
for key in sorted(G01.keys()):
#    print key
    Data_g01 = np.append(Data_g01,[G01[key][0,:]], axis=0)

Data_g01.sort(axis=0)
#--------------------------------------------------------------------------
#       File for gain parameter 02
#--------------------------------------------------------------------------
file_g02 = open('validation_009.rsn', 'r')
gain02 = file_g02.readlines()
g02 = pd.Series(gain02)
g02 = g02.replace(r'\n','', regex=True)
g02 = g02.replace(r'\r\n','', regex=True)
g02 = g02.replace(r'\r','', regex=True)
index_Time_g02 = g02[g02.str.contains('TIME', case=False, regex=False)]
index_TimeValues_g02 = index_Time_g02.index.values
#--------------------------------------------------------------------------
G02 = {}
for idx in index_Time_g02.index.values:
    index_start = idx + 1
    index_end = index_start + NbrOfNodes
    tmp_df = g02[index_start:index_end].str.strip()
    tmp_df = tmp_df.str.split('  ',expand=True)
    np.array(tmp_df.values, dtype=float)
    G02[g02[idx]]=np.array(tmp_df.values, dtype=float)

Data_g02 =  np.array([], dtype=np.int64)
Data_g02.shape = (-1, 7)
for key in sorted(G02.keys()):
#    print key
    Data_g02 = np.append(Data_g02,[G02[key][0,:]], axis=0)

Data_g02.sort(axis=0)
#--------------------------------------------------------------------------
#       File for gain parameter 03
#--------------------------------------------------------------------------
file_g03 = open('validation_015.rsn', 'r')
gain03 = file_g03.readlines()
g03 = pd.Series(gain03)
g03 = g03.replace(r'\n','', regex=True)
g03 = g03.replace(r'\r\n','', regex=True)
g03 = g03.replace(r'\r','', regex=True)
index_Time_g03 = g03[g03.str.contains('TIME', case=False, regex=False)]
index_TimeValues_g03 = index_Time_g03.index.values
#--------------------------------------------------------------------------
G03 = {}
for idx in index_Time_g03.index.values:
    index_start = idx + 1
    index_end = index_start + NbrOfNodes
    tmp_df = g03[index_start:index_end].str.strip()
    tmp_df = tmp_df.str.split('  ',expand=True)
    np.array(tmp_df.values, dtype=float)
    G03[g03[idx]]=np.array(tmp_df.values, dtype=float)

Data_g03 =  np.array([], dtype=np.int64)
Data_g03.shape = (-1, 7)
for key in sorted(G03.keys()):
#    print key
    Data_g03 = np.append(Data_g03,[G03[key][0,:]], axis=0)

Data_g03.sort(axis=0)
#--------------------------------------------------------------------------
#       File for gain parameter 04
#--------------------------------------------------------------------------
file_g04 = open('cyron2016b_009.png.dat', 'r')
gain04 = file_g04.readlines()
g04 = pd.Series(gain04)
g04 = g04.replace(r'\n','', regex=True)
g04 = g04.replace(r'\r\n','', regex=True)
g04 = g04.replace(r'\r','', regex=True)
#--------------------------------------------------------------------------
for line in g04:
    tmp_df = g04.str.strip()
    tmp_df = tmp_df.str.split('  ',expand=True)
    np.array(tmp_df.values, dtype=float)
    Data_g04 = np.array(tmp_df.values, dtype=float)

#--------------------------------------------------------------------------
#       File for gain parameter 05
#--------------------------------------------------------------------------
file_g05 = open('cyron2016b_005.png.dat', 'r')
gain05 = file_g05.readlines()
g05 = pd.Series(gain05)
g05 = g05.replace(r'\n','', regex=True)
g05 = g05.replace(r'\r\n','', regex=True)
g05 = g05.replace(r'\r','', regex=True)
#--------------------------------------------------------------------------
for line in g05:
    tmp_df = g05.str.strip()
    tmp_df = tmp_df.str.split('  ',expand=True)
    np.array(tmp_df.values, dtype=float)
    Data_g05 = np.array(tmp_df.values, dtype=float)

#--------------------------------------------------------------------------
#       File for gain parameter 06
#--------------------------------------------------------------------------
file_g06 = open('cyron2016b_015.png.dat', 'r')
gain06 = file_g06.readlines()
g06 = pd.Series(gain06)
g06 = g06.replace(r'\n','', regex=True)
g06 = g06.replace(r'\r\n','', regex=True)
g06 = g06.replace(r'\r','', regex=True)
#--------------------------------------------------------------------------
for line in g06:
    tmp_df = g06.str.strip()
    tmp_df = tmp_df.str.split('  ',expand=True)
    np.array(tmp_df.values, dtype=float)
    Data_g06 = np.array(tmp_df.values, dtype=float)

#--------------------------------------------------------------------------
fig = plt.figure()
plt.rcParams.update({'font.size': 12})

#plt.subplot(2,3,3)
plt.plot(Data_g01[:,0],Data_g01[:,3]*1000.0,'b',label='$k_{\sigma}$=0.05/T',linewidth=2.0,markersize=10)
plt.plot(Data_g02[:,0],Data_g02[:,3]*1000.0,'r',label='$k_{\sigma}$=0.09/T',linewidth=2.0,markersize=10)
plt.plot(Data_g03[:,0],Data_g03[:,3]*1000.0,'g',label='$k_{\sigma}$=0.15/T',linewidth=2.0,markersize=10)
plt.plot(Data_g04[:,0]*12.0,Data_g04[:,1],'r--',linewidth=2.0,markersize=10)
plt.plot(Data_g05[:,0]*12.0,Data_g05[:,1],'b--',linewidth=2.0,markersize=10)
plt.plot(Data_g06[:,0]*12.0,Data_g06[:,1],'g--',linewidth=2.0,markersize=10)
plt.xlabel('Time [months]')
plt.ylabel('Radius [mm]')
plt.legend(loc='upper left')
plt.xlim(right=180)
plt.ylim(bottom=10)
#plt.title('Radius')

fig.tight_layout()
plt.show

FIGURENAME = 'validation.eps'
plt.savefig(FIGURENAME)
plt.savefig(fname=FIGURENAME,
            dpi=None,
            facecolor='w',
            edgecolor='w',
            orientation='portrait',
            format=None,
            transparent=False,
            bbox_inches=None,
            pad_inches=0.1,
            frameon=None,
            metadata=None)

plt.close('all')

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

max_stars_repo_name: jdlaubrie/shell-elem

max_issues_repo_name: jdlaubrie/shell-elem

max_forks_repo_name: jdlaubrie/shell-elem


text:
import os
import json
import numpy as np
import re
import xml.etree.ElementTree
from speaksee.data import field, Example, PairedDataset, COCO
from speaksee.utils import nostdout
from itertools import groupby
import pickle as pkl

class COCOEntities(PairedDataset):
    def __init__(self, image_field, det_field, text_field, img_root, ann_root, entities_file, id_root=None,
                 use_restval=True, filtering=False):
        roots = dict()
        roots['train'] = {
            'img': os.path.join(img_root, 'train2014'),
            'cap': os.path.join(ann_root, 'captions_train2014.json')
        }
        roots['val'] = {
            'img': os.path.join(img_root, 'val2014'),
            'cap': os.path.join(ann_root, 'captions_val2014.json')
        }
        roots['test'] = {
            'img': os.path.join(img_root, 'val2014'),
            'cap': os.path.join(ann_root, 'captions_val2014.json')
        }
        roots['trainrestval'] = {
            'img': (roots['train']['img'], roots['val']['img']),
            'cap': (roots['train']['cap'], roots['val']['cap'])
        }

        if id_root is not None:
            ids = {}
            ids['train'] = np.load(os.path.join(id_root, 'coco_train_ids.npy'))
            ids['val'] = np.load(os.path.join(id_root, 'coco_dev_ids.npy'))
            ids['test'] = np.load(os.path.join(id_root, 'coco_test_ids.npy'))
            ids['trainrestval'] = (
                ids['train'],
                np.load(os.path.join(id_root, 'coco_restval_ids.npy')))

            if use_restval:
                roots['train'] = roots['trainrestval']
                ids['train'] = ids['trainrestval']
        else:
            ids = None

        if not filtering:
            dataset_path = 'saved_data/coco/coco_entities_precomp.pkl'
        else:
            dataset_path = 'saved_data/coco/coco_entities_filtered_precomp.pkl'
        if not os.path.isfile(dataset_path):
            with nostdout():
                train_examples, val_examples, test_examples = COCO.get_samples(roots, ids)

            self.train_examples, self.val_examples, self.test_examples = self.get_samples([train_examples, val_examples, test_examples], entities_file, filtering)
            pkl.dump((self.train_examples, self.val_examples, self.test_examples), open(dataset_path, 'wb'), -1)
        else:
            self.train_examples, self.val_examples, self.test_examples = pkl.load(open(dataset_path, 'rb'))

        examples = self.train_examples + self.val_examples + self.test_examples
        super(COCOEntities, self).__init__(examples, {'image': image_field, 'detection': det_field, 'text': text_field})

    @property
    def splits(self):
        train_split = PairedDataset(self.train_examples, self.fields)
        val_split = PairedDataset(self.val_examples, self.fields)
        test_split = PairedDataset(self.test_examples, self.fields)
        return train_split, val_split, test_split

    @classmethod
    def get_samples(cls, samples, entities_file, filtering=False):
        train_examples = []
        val_examples = []
        test_examples = []

        with open(entities_file, 'r') as fp:
            visual_chunks = json.load(fp)

        for id_split, samples_split in enumerate(samples):
            for s in samples_split:
                id_image = str(int(s.image.split('/')[-1].split('_')[-1].split('.')[0]))
                caption = s.text.lower().replace('\t', ' ').replace('\n', '')

                words = caption.strip().split(' ')
                caption_fixed = []
                for w in words:
                    if w not in field.TextField.punctuations and w != '':
                        caption_fixed.append(w)

                det_classes = [None for _ in caption_fixed]
                caption_fixed = ' '.join(caption_fixed)

                for p in field.TextField.punctuations:
                    caption_fixed = caption_fixed.replace(p, '')

                if id_image in visual_chunks:
                    if caption in visual_chunks[id_image]:
                        chunks = visual_chunks[id_image][caption]
                        for chunk in chunks:
                            words = chunk[0].split(' ')
                            chunk_fixed = []
                            for w in words:
                                if w not in field.TextField.punctuations and w != '':
                                    chunk_fixed.append(w)
                            chunk_fixed = ' '.join(chunk_fixed)
                            for p in field.TextField.punctuations:
                                chunk_fixed = chunk_fixed.replace(p, '')

                            sub_str = ' '.join(['_' for _ in chunk_fixed.split(' ')])
                            sub_cap = caption_fixed.replace(chunk_fixed, sub_str).split(' ')
                            for i, w in enumerate(sub_cap):
                                if w == '_':
                                    det_classes[i] = chunk[1]

                        example = Example.fromdict({'image': s.image,
                                                    'detection': (s.image, tuple(det_classes)),
                                                    'text': caption_fixed})

                        det_classes_set = [x[0] for x in groupby(det_classes) if x[0] is not None]
                        chunks_filtered = list(set([c[1] for c in chunks]))
                        if len(det_classes_set) < len(chunks_filtered):
                            pass
                        else:
                            if id_split == 0:
                                train_examples.append(example)
                            elif id_split == 1:
                                if filtering:
                                    if '_' not in example.detection[1]:
                                        val_examples.append(example)
                                else:
                                    val_examples.append(example)
                            elif id_split == 2:
                                if filtering:
                                    if '_' not in example.detection[1]:
                                        test_examples.append(example)
                                else:
                                    test_examples.append(example)

        return train_examples, val_examples, test_examples

class FlickrEntities(PairedDataset):
    def __init__(self, image_field, text_field, det_field, img_root, ann_file, entities_root,
                 precomp_file='saved_data/flickr/flickr_entities_precomp.pkl'):
        if os.path.isfile(precomp_file):
            with open(precomp_file, 'rb') as pkl_file:
                self.train_examples, self.val_examples, self.test_examples = pkl.load(pkl_file)
        else:
            self.train_examples, self.val_examples, self.test_examples = self.get_samples(ann_file, img_root, entities_root)
            with open(precomp_file, 'wb') as pkl_file:
                pkl.dump((self.train_samples, self.val_samples, self.test_samples), pkl_file, protocol=pkl.HIGHEST_PROTOCOL)

        examples = self.train_examples + self.val_examples + self.test_examples
        super(FlickrEntities, self).__init__(examples, {'image': image_field, 'detection': det_field, 'text': text_field})

    @property
    def splits(self):
        train_split = PairedDataset(self.train_examples, self.fields)
        val_split = PairedDataset(self.val_examples, self.fields)
        test_split = PairedDataset(self.test_examples, self.fields)
        return train_split, val_split, test_split

    def get_samples(self, ann_file, img_root, entities_root):
        def _get_sample(d):
            filename = d['filename']
            split = d['split']
            xml_root = xml.etree.ElementTree.parse(os.path.join(entities_root, 'Annotations',
                                                                filename.replace('.jpg', '.xml'))).getroot()
            det_dict = dict()
            id_counter = 1
            for obj in xml_root.findall('object'):
                obj_names = [o.text for o in obj.findall('name')]
                if obj.find('bndbox'):
                    bbox = tuple(int(o.text) for o in obj.find('bndbox'))
                    for obj_name in obj_names:
                        if obj_name not in det_dict:
                            det_dict[obj_name] = {'id': id_counter, 'bdnbox': [bbox]}
                            id_counter += 1
                        else:
                            det_dict[obj_name]['bdnbox'].append(bbox)

            bdnboxes = [[] for _ in range(id_counter - 1)]
            for it in det_dict.values():
                bdnboxes[it['id'] - 1] = tuple(it['bdnbox'])
            bdnboxes = tuple(bdnboxes)

            captions = [l.strip() for l in open(os.path.join(entities_root, 'Sentences',
                                                             filename.replace('.jpg', '.txt')), encoding="utf-8").readlines()]
            outputs = []
            for c in captions:
                matches = prog.findall(c)
                caption = []
                det_ids = []

                for match in matches:
                    for i, grp in enumerate(match):
                        if i in (0, 2):
                            if grp != '':
                                words = grp.strip().split(' ')
                                for w in words:
                                    if w not in field.TextField.punctuations and w != '':
                                        caption.append(w)
                                        det_ids.append(0)
                        elif i == 1:
                            words = grp[1:-1].strip().split(' ')
                            obj_name = words[0].split('#')[-1].split('/')[0]
                            words = words[1:]
                            for w in words:
                                if w not in field.TextField.punctuations and w != '':
                                    caption.append(w)
                                    if obj_name in det_dict:
                                        det_ids.append(det_dict[obj_name]['id'])
                                    else:
                                        det_ids.append(0)

                caption = ' '.join(caption)
                if caption != '' and np.sum(np.asarray(det_ids)) > 0:
                    example = Example.fromdict({'image': os.path.join(img_root, filename),
                                                'detection': (os.path.join(img_root, filename), bdnboxes, det_ids),
                                                'text': caption})
                    outputs.append([example, split])

            return outputs

        train_samples = []
        val_samples = []
        test_samples = []

        prog = re.compile(r'([^\[\]]*)(\[[^\[\]]+\])([^\[\]]*)')
        dataset = json.load(open(ann_file, 'r'))['images']

        samples = []
        for d in dataset:
            samples.extend(_get_sample(d))

        for example, split in samples:
            if split == 'train':
                train_samples.append(example)
            elif split == 'val':
                val_samples.append(example)
            elif split == 'test':
                test_samples.append(example)

        return train_samples, val_samples, test_samples


text:
'''
want to demonstrate a bit of info about ROC curve, AUROC, and PR curves as
    metrics. will use scikit-learn for items

'''

import numpy as np
import sklearn.metrics as skm
import matplotlib.pyplot as plt
np.random.seed(0)
# generate some fake data
ndat = 30
k=0.30 # factor for how much prediction deviates from ytrue. k:[0.0-0.5]
ytrue = np.random.randint(2,size=ndat)
ypred = ((ytrue*k+(1-k)*np.random.rand(ndat))).round(0).astype(int)
yconf = np.random.rand(ndat)/2.1+ypred*0.5

print('data')
print('ytrue:',ytrue)
print('ypred:',ypred)
print('yconf:',yconf.round(2))
CM = skm.confusion_matrix(ytrue,ypred)
print('CM:\n',CM,sep='')

ROC = skm.roc_curve(ytrue,yconf)
AUROC = skm.roc_auc_score(ytrue,yconf)
print('AUROC:',AUROC)

PR = skm.precision_recall_curve(ytrue,yconf)
AP = skm.average_precision_score(ytrue,yconf)
print('AP:',AP)

f,p=plt.subplots(1,2)
p[0].plot(ROC[0],ROC[1],'r.-')
p[0].set_title('ROC')
p[0].set_xlim([0,1])
p[0].set_ylim([0,1])
p[0].set_xlabel('FPR')
p[0].set_ylabel('TPR')
p[0].grid()
p[0].set_aspect('equal')

p[1].plot(PR[1],PR[0],'r.-')
p[1].set_title('PR')
p[1].set_xlim([0,1])
p[1].set_ylim([0,1])
p[1].set_xlabel('recall')
p[1].set_ylabel('precision')
p[1].grid()
p[1].set_aspect('equal')

plt.show()

# import ipdb; ipdb.set_trace()
# eof

max_stars_repo_name: kjgonzalez/codefiles

max_issues_repo_name: kjgonzalez/codefiles

max_forks_repo_name: kjgonzalez/codefiles


text:
import os,datetime
from dateutil import rrule
import scipy as sp
from g5lib import dset

__all__=['tctl','sctl','rhoctl']

class Ctl(dset.NCDset):
    def __init__(self,var='t'):
        name='WOA13 '+var.upper()

        oceanval=os.environ.get('OCEANVAL',
                                '/discover/nobackup/projects/gmao/oceanval/verification')
        data_dir=oceanval+'/WOA13'
        if var=='I':
            ver='01'
        else:
            ver='01v2'
        flist=sp.array([data_dir+'/woa13_decav_'+var+str(mm).zfill(2)+'_'+ver+'.nc' for mm in range(1,13)],dtype='|O')

        rr=rrule.rrule(rrule.MONTHLY,dtstart=datetime.date(2013,1,1),count=12)
        time=sp.array(rr[:],dtype='|O')

        super(Ctl,self).__init__(flist,time=time,name=name,levname='depth')

tctl=Ctl('t')
sctl=Ctl('s')
rhoctl=Ctl('I')

reason: import scipy