text:
import os
import dotenv
from pathlib import Path
import random

from scipy import stats

import nibabel as nib
import matplotlib.pyplot as plt
from msm.run import run_msm, prepare_darrays

dotenv.load_dotenv()
FSL_PATH = os.getenv("FSL_PATH")
FSL_CONFIG_PATH = os.getenv("FSL_CONFIG_PATH")

# Load spherical mesh produced with `mris_convert`
spherical_mesh = "./data/lh.sphere.gii"

# load data
subject_input = "sub-07"
subject_reference = "sub-04"
data_path = "./data/"

# parameters
random.seed(42)
lambd = [0.075, 0.1, 0.15]
# if true will run MSM, if false will use an existing run_MSM output
train_mode = True

sessions = [
    ["ses-00"],
    ["ses-00", "ses-01", "ses-02", "ses-04"],
    ["ses-00", "ses-01"],
    ["ses-00", "ses-01", "ses-02"],
]

for session in sessions:
    # data - dictionary of type: input_fname -> reference_fname

    data = {}

    data_input = [
        os.path.join(data_path, file)
        for file in os.listdir(data_path)
        if subject_input in file
        and "lh" in file
        and (any(ses in file for ses in session))
    ]

    data_reference = [
        os.path.join(data_path, file)
        for file in os.listdir(data_path)
        if subject_reference in file
        and "lh" in file
        and (any(ses in file for ses in session))
    ]

    for input_fname in data_input:
        for reference_fname in data_reference:
            reference_contrast = reference_fname.split(subject_reference)[1]
            input_contrast = input_fname.split(subject_input)[1]
            if reference_contrast == input_contrast:
                data[input_fname] = reference_fname

    print(len(data))
    # split data to train and test

    train_data_input = random.sample(list(data.keys()), int(len(data) * 0.8))

    train_data = {
        input_fname: data[input_fname] for input_fname in train_data_input
    }

    test_data = {
        input_fname: data[input_fname]
        for input_fname in list(data.keys())
        if input_fname not in train_data_input
    }

    cross_correlation = {}
    # cross-correlation without deformation
    base_cross_correlation = {}
    for lam in lambd:

        # replacing lambda values in config file
        with open(f"{FSL_CONFIG_PATH}", "r") as f:
            list_of_lines = f.readlines()
            list_of_lines[3] = f"--lambda={lam},{lam},{lam},{lam}\n"

        with open(f"{FSL_CONFIG_PATH}", "w") as f:
            f.writelines(list_of_lines)

        # running MSM with train data
        output_dir = "test_outputs_lambda_" + str(len(data)) + "_" + str(lam)

        if train_mode and not os.path.isdir(output_dir):
            transformed_mesh, transformed_func = run_msm(
                in_data_list=list(train_data.keys()),
                in_mesh=spherical_mesh,
                ref_data_list=list(train_data.values()),
                debug=False,
                verbose=True,
                output_dir=output_dir,
            )

        # testing MSM with test data

        deformed = os.path.join(output_dir, "transformed_in_mesh.surf.gii")
        cross_correlation[lam] = 0
        base_cross_correlation[lam] = 0
        # Load the coordsys from the mesh associated to the data to make
        # sure it is well specified
        mesh = nib.load(spherical_mesh)
        coordsys = mesh.darrays[0].coordsys

        for input_fname, reference_fname in test_data.items():
            data_input = nib.load(input_fname)
            data_input.darrays = prepare_darrays(data_input.darrays, coordsys)
            # copy dataset in order to cope with a bug of MSM
            extra_data = nib.load(input_fname)
            data_input.darrays.extend(
                prepare_darrays(extra_data.darrays, coordsys)
            )
            filename = str(Path(output_dir) / "input_test.func.gii")
            data_input.to_filename(filename)
            input_fname_preprocessed = filename

            input_contrast_name = input_fname.split(subject_input)[1]
            test_transformed = os.path.join(
                output_dir,
                f"test_transformed_and_projected{input_contrast_name}",
            )
            cmd = " ".join(
                [
                    f"{FSL_PATH}/fsl/bin/msmresample",
                    f"{deformed} ",
                    test_transformed,
                    f"-labels {input_fname_preprocessed} ",
                    f"-project {spherical_mesh}",
                ]
            )
            exit_code = os.system(cmd)
            if exit_code != 0:
                raise RuntimeError(f"Failed to run MSM with command:\n{cmd}")

            # compute cross_correlation
            transformed_data_path = f"{test_transformed}.func.gii"
            transformed_data = nib.load(transformed_data_path).darrays[0].data
            reference_data = nib.load(reference_fname).darrays[0].data

            cross_correlation[lam] += stats.pearsonr(
                transformed_data, reference_data
            )[0] / len(test_data)
            base_cross_correlation[lam] += stats.pearsonr(
                data_input, reference_data
            )[0] / len(test_data)

    print(cross_correlation)
    print(base_cross_correlation)
    f1 = plt.figure()
    plt.xlabel("Lambda")
    plt.ylabel("Pearson correlation")
    plt.title(f"Test data, data number: {len(data)}")
    plt.plot(
        lambd,
        list(cross_correlation.values()),
        "-o",
        label="Pearson correlation after transformation",
    )
    plt.plot(
        lambd,
        list(base_cross_correlation.values()),
        "-o",
        label="Pearson correlation before transformation",
    )
    plt.legend()
    plt.savefig(str(len(data)) + "lambda_optimization_test_data.png")

max_stars_repo_name: wingrune/Parietal

max_issues_repo_name: wingrune/Parietal

max_forks_repo_name: wingrune/Parietal

reason: from scipy


text:
#!/usr/bin/env python
# coding: utf-8

# In[30]:

import numpy as np
from astropy.io import fits
from astropy.table import Table, vstack
from astropy.wcs import WCS
import os
import time
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from matplotlib import cm
import pandas as pd
pd.options.display.max_columns = 250
pd.options.display.max_rows = 250
import healpy as hp
from copy import copy, deepcopy
import logging, traceback
import sys
logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)

# In[31]:

# cd to code directory
get_ipython().run_line_magic('cd', '/Users/tparsota/Documents/BAT_SCRIPTS/NITRATES/NITRATES_GIT/NITRATES')
#/storage/work/j/jjd330/local/bat_data/BatML

# In[32]:

from event2dpi_funcs import det2dpis, mask_detxy
from flux_models import Cutoff_Plaw_Flux, Plaw_Flux, get_eflux_from_model
from sqlite_funcs import get_conn
from dbread_funcs import get_info_tab
from do_manage2 import im_dist, get_rate_res_fnames, get_peak_res_fnames, get_out_res_fnames,                    get_merged_csv_df, get_merged_csv_df_wpos
from hp_funcs import ang_sep
from coord_conv_funcs import theta_phi2imxy, imxy2theta_phi, convert_imxy2radec,                            convert_radec2thetaphi, convert_radec2imxy
from do_llh_inFoV4realtime2 import parse_bkg_csv
from LLH import LLH_webins
from minimizers import NLLH_ScipyMinimize_Wjacob

# In[33]:

ebins0 = np.array([15.0, 24.0, 35.0, 48.0, 64.0])
ebins0 = np.append(ebins0, np.logspace(np.log10(84.0), np.log10(500.0), 5+1))[:-1]
ebins0 = np.round(ebins0, decimals=1)[:-1]
ebins1 = np.append(ebins0[1:], [350.0])
nebins = len(ebins0)
print("Number of ebins: ", nebins)

# directory with results
work_dir = '/Users/tparsota/Documents/BAT_SCRIPTS/NITRATES/F646018360/'
#'/storage/work/j/jjd330/local/bat_data/realtime_workdir/F646018360/'

conn = get_conn(os.path.join(work_dir,'results.db'))
info_tab = get_info_tab(conn)
# trigger time stored in results.db (sqlite DB)
# not much else stored there (results used to be stored there but not anymore)
trigger_time = info_tab['trigtimeMET'][0]
print("trigger time: ", trigger_time)

# default file names made by do_data_setup.py
evfname = os.path.join(work_dir,'filter_evdata.fits')
ev_data = fits.open(evfname)[1].data
# GTI extensions added to event fits file by do_data_setup.py
GTI_PNT = Table.read(evfname, hdu='GTI_POINTING') # when the analysis can be run (good data time and pointing)
GTI_SLEW = Table.read(evfname, hdu='GTI_SLEW') # slewing times
# the final set of detectors to mask (disabled dets, hot/cold dets, and dets with glitches)
dmask = fits.open(os.path.join(work_dir,'detmask.fits'))[0].data
attfile = fits.open(os.path.join(work_dir,'attitude.fits'))[1].data

# number of detectors being used
ndets = np.sum(dmask==0)
print ("Ndets: ", np.sum(dmask==0))

t_end = trigger_time + 1e3
t_start = trigger_time - 1e3
mask_vals = mask_detxy(dmask, ev_data)
bl_dmask = (dmask==0.)

# get rid of events:
# far away from trigger time
# from bad dets
# with bad event flags
bl_ev = (ev_data['EVENT_FLAGS']<1)&        (ev_data['ENERGY']<=500.)&(ev_data['ENERGY']>=14.)&        (mask_vals==0.)&(ev_data['TIME']<=t_end)&        (ev_data['TIME']>=t_start)

print("Number of events passing cuts: ", np.sum(bl_ev))
ev_data0 = ev_data[bl_ev]

# In[34]:

print (GTI_PNT)
print()
for row in GTI_PNT:
    print (row['START'] - trigger_time, row['STOP'] - trigger_time)

# In[35]:

bkg_fname = os.path.join(work_dir,'bkg_estimation.csv')
bkg_df = pd.read_csv(bkg_fname)

# In[36]:

bkg_df.head()

# In[37]:

# plotting the rate per detector in each energy bin for the diffuse models
nr = 3
nc = 3
nplt = 1
fig = plt.figure(dpi=80, figsize=(14,6))
for j in range(9):
    ax = fig.add_subplot(nr,nc,nplt)
    nplt+=1
    try:
        rate_name = 'bkg_rate_' + str(j)
        rate = bkg_df[rate_name]
    except:
        rate_name = 'Background_bkg_rate_' + str(j)
        rate = bkg_df[rate_name]
    err0 = rate - bkg_df['err_'+rate_name]
    err1 = rate + bkg_df['err_'+rate_name]
    plt.fill_between(bkg_df['dt'], err0, err1, alpha=.5)
    plt.plot(bkg_df['dt'], rate, 'o--')
    plt.grid(True)

# In[38]:

# plotting time hist of event data
tbins = np.arange(-40*1.024, 40*1.024, 0.064*4*4*1) # tbins with multiples of 64ms (64ms * 4 * 4 = 1.024s)
dt = tbins[1] - tbins[0]
tax = (tbins[1:] + tbins[:-1])/2.
ntbins = len(tax)
# choosing which tbins to estimate a flat bkg from for the plot
# bkg_bl = (np.abs(tax)>8.0)
bkg_bl = ((tax)<-10.0)&(tax>-40.0)

nc = 1
nr = 1
nplt = 1

fig = plt.figure(dpi=100, figsize=(8,3*nr))
ax = fig.add_subplot(nr,nc,nplt)
nplt += 1

h=plt.hist(ev_data0['TIME'] - trigger_time, bins=tbins,           histtype='step', label='data')[0]

bkg_mean = np.mean(h[bkg_bl])
bkg_std = np.std(h[bkg_bl])

print("time with min counts, max counts")
print (tbins[np.argmin(h)], tbins[np.argmax(h)])
print("counts at min time, max time")
print(np.min(h), np.max(h))
print("bkg mean, bkg std")
print(bkg_mean, bkg_std)

plt.axhline(bkg_mean)
plt.grid(True)
plt.xlim(-22, 26)
# plt.xlim(-4, 4.)
plt.ylim(np.around(.95*np.min(h[(h>0)]), decimals=-1),         np.round(1.01*np.max(h), decimals=-1))
plt.xlabel('t - trig_time (s)')
plt.ylabel('Counts 15-350 keV')
plt.legend(loc='lower left')

# In[39]:

# plotting time hist with same bins as last cell but for each Ebin
nc = 1
nr = nebins
nplt = 1
fig = plt.figure(dpi=100, figsize=(8,3.5*nr))

for ei in range(nebins):
    ax = fig.add_subplot(nr,nc,nplt)
    nplt += 1
    ebl = (ev_data0['ENERGY']>=ebins0[ei])&            (ev_data0['ENERGY']<ebins1[ei])
    h=plt.hist(ev_data0[ebl]['TIME'] - trigger_time, bins=tbins,               histtype='step')[0]
    bkg_mean = np.mean(h[bkg_bl])
    bkg_std = np.std(h[bkg_bl])
    ttl = '%.1f - %.1f keV' %(ebins0[ei],ebins1[ei])
    print(ttl)
    print("time with min counts, max counts")
    print (tbins[np.argmin(h)], tbins[np.argmax(h)])
    print("counts at min time, max time")
    print(np.min(h), np.max(h))
    print("bkg mean, bkg std")
    print(bkg_mean, bkg_std)
    print
    plt.axhline(bkg_mean)
    plt.title(ttl)
    plt.grid(True)
#     plt.xlim(-2, 2)
    plt.xlim(-22, 26)
    plt.ylim(np.around(.9*np.min(h[(h>0)]), decimals=-1),             np.round(1.1*np.max(h), decimals=-1))

# In[40]:

# getting attitude information at trigger time
att_ind = np.argmin(np.abs(attfile['TIME'] - trigger_time))
att_quat = attfile['QPARAM'][att_ind]
print("QUATERNION: ", att_quat)
pnt_ra, pnt_dec = attfile['POINTING'][att_ind,:2]
print("Pointing RA, Dec")
print(pnt_ra, pnt_dec)
plt.plot(attfile['TIME']-trigger_time, attfile['POINTING'], 'o')
plt.grid(True)
plt.xlim(-50,150)
plt.legend(['ra','dec','roll'])

# In[41]:

# if you have some ra, dec of interest (like a gbm localization or something)
# here's how to do conversions into detector coordinates (theta, phi) and (imx, imy)
ra_interest, dec_interest = 233.117, -26.213
print("RA of interest, Dec of interest")
print(ra_interest, dec_interest)
theta_interest, phi_interest = convert_radec2thetaphi(ra_interest, dec_interest, att_quat)
print("theta, phi")
print (theta_interest, phi_interest)
# imx, imy only valid at theta < 90 deg
imx_interest, imy_interest = convert_radec2imxy(ra_interest, dec_interest, att_quat)
print("imx, imy")
print (imx_interest, imy_interest)

# In[42]:

#%%time
# getting split rate analysis results

# get the file names
res_rate_fnames = get_rate_res_fnames(work_dir)
print("%d split rate result files"%(len(res_rate_fnames)))
# read files and merge into one Pandas Dataframe
res_rate_tab = get_merged_csv_df(res_rate_fnames, work_dir, ignore_index=True)
print("Merged split rate results table has %d rows"%(len(res_rate_tab)))
res_rate_tab['dt'] = res_rate_tab['time'] - trigger_time

# In[43]:

res_rate_tab.sort_values('TS', ascending=False).head(64)

# In[44]:

#%%time
# getting out of FoV analysis results

res_out_fnames = get_out_res_fnames(work_dir)
print("%d out of FoV result files"%(len(res_out_fnames)))

res_out_tab = get_merged_csv_df_wpos(res_out_fnames, attfile, direc=work_dir, ignore_index=True)
print("Merged out of FoV results table has %d rows"%(len(res_out_tab)))
res_out_tab['dt'] = res_out_tab['time'] - trigger_time

# In[45]:

#%%time
# getting in FoV analysis results

res_peak_fnames = get_peak_res_fnames(work_dir)
print("%d in FoV peaks result files"%(len(res_peak_fnames)))
# read files and merge in single dataframe, also convert detector coords into RA, Dec
res_peak_tab = get_merged_csv_df_wpos(res_peak_fnames, attfile, direc=work_dir, ignore_index=True)
print("Merged in FoV peaks results table has %d rows"%(len(res_peak_tab)))
res_peak_tab['dt'] = res_peak_tab['time'] - trigger_time

# In[46]:

# getting the max TS for each square/time seed combo
idx = res_peak_tab.groupby(['squareID','timeID'])['TS'].transform(max) == res_peak_tab['TS']
res_peak_maxSq_tab = res_peak_tab[idx]
print(len(res_peak_maxSq_tab))

# In[47]:

# initing flux model to calculate flux/fluences
flux_params = {'A':1.0, 'gamma':0.5, 'Epeak':1e2}
flux_mod = Cutoff_Plaw_Flux(E0=100.0)

# In[48]:

#%%time
# calculating fluence for each row based on the best fit spectral parameters
fluncs = np.zeros(len(res_peak_maxSq_tab))
flux_pars = {'A':1.0, 'Epeak':1e2, 'gamma':0.5}
i = 0
for ind, row in res_peak_maxSq_tab.iterrows():
    flux_pars['A'] = row['A']
    flux_pars['gamma'] = row['gamma']
    flux_pars['Epeak'] = row['Epeak']
    fluncs[i] = get_eflux_from_model(flux_mod, flux_pars, 1e1, 1e3)*row['dur'] # fluence from 10keV to 1MeV
    i+=1
#     res_peak_maxSq_tab.loc[ind]['fluence'] = get_eflux_from_model(flux_mod, flux_pars, 1e1, 1e3)*row['dur']
res_peak_maxSq_tab['fluence'] = fluncs

# In[49]:

# max TS peak result
print (np.max(res_peak_maxSq_tab['TS']))
idx = res_peak_maxSq_tab['TS'].idxmax()
row = res_peak_maxSq_tab.loc[idx]
max_TS_timeID = row['timeID']
row

# In[50]:

res_peak_maxSq_tab.sort_values('TS', ascending=False).head(64)

# In[51]:

# max TS out result
print (np.max(res_out_tab['TS']))
idx = res_out_tab['TS'].idxmax()
row = res_out_tab.loc[idx]
max_TSout_timeID = row['timeID']
row

# In[52]:

res_out_tab.sort_values('TS', ascending=False).head(64)

# In[53]:

# get best out of FoV result for each healpix pixel for a certrain time bin
bl = np.isclose(res_out_tab['timeID'],max_TSout_timeID)
idx = res_out_tab[bl].groupby(['hp_ind'])['TS'].transform(max) == res_out_tab[bl]['TS']
res_hpmax_tab = res_out_tab[bl][idx]

# In[54]:

#%%time
# get fluence for each row of res_hpmax_tab
fluncs = np.zeros(len(res_hpmax_tab))
flux_pars = {'A':1.0, 'Epeak':1e2, 'gamma':0.5}
# for i in range(len(res_peak_maxSq_tab)):
i = 0
for ind, row in res_hpmax_tab.iterrows():
    flux_pars['A'] = row['A']
    flux_pars['gamma'] = row['gamma']
    flux_pars['Epeak'] = row['Epeak']
    fluncs[i] = get_eflux_from_model(flux_mod, flux_pars, 1e1, 1e3)*row['dur']
    i+=1
res_hpmax_tab['fluence'] = fluncs

# In[55]:

res_hpmax_tab.sort_values('TS', ascending=False).head(64)

# In[56]:

def get_dlogl_peak_out(res_peak_tab, res_out_tab):
    '''
    returns DeltaLLH_peak and DeltaLLH_out for the time bin with the max TS
    '''

    idx = res_peak_tab['TS'].idxmax()
    row = res_peak_tab.loc[idx]

    timeID = row['timeID']
    imdists = im_dist(row['imx'], row['imy'], res_peak_tab['imx'], res_peak_tab['imy'])
    bld = (imdists>0.012)&(res_peak_tab['timeID']==timeID)

    dlogl_peak = np.nanmin(res_peak_tab[bld]['nllh']) - row['nllh']

    blo = (res_out_tab['timeID']==timeID)

    dlogl_out = np.nanmin(res_out_tab[blo]['nllh']) - row['nllh']

    return dlogl_peak, dlogl_out

def get_dlogls_inout(res_tab, res_out_tab):
    '''
    returns DeltaLLH_peak, DeltaLLH_out, timeID for each time bin
    '''

    dlogls = []
    dlogls_in_out = []
    timeIDs = []
    for timeID, df in res_tab.groupby('timeID'):
        idx = df['TS'].idxmax()
        row = df.loc[idx]
        imdists = im_dist(row['imx'], row['imy'], df['imx'], df['imy'])
        bld = (imdists>0.012)
        try:
            dlogls.append(np.nanmin(df[bld]['nllh']) - row['nllh'])
        except Exception as E:
            print(E)
            dlogls.append(np.nan)
        blo = np.isclose(res_out_tab['timeID'],timeID,rtol=1e-9,atol=1e-3)
        dlogls_in_out.append(np.nanmin(res_out_tab[blo]['nllh']) - row['nllh'])
        timeIDs.append(timeID)
    return dlogls, dlogls_in_out, timeIDs

# In[57]:

dlogl_peak, dlogl_out = get_dlogl_peak_out(res_peak_tab, res_out_tab)
print("\Delta LLH_peak = %.3f"%(dlogl_peak))
print("\Delta LLH_out = %.3f"%(dlogl_out))

# In[58]:

# %%time
dlogls_peak, dlogls_out, list_timeIDs = get_dlogls_inout(res_peak_tab, res_out_tab)
for i in range(len(list_timeIDs)):
    print( "timeID = %d"%(list_timeIDs[i]))
    print( "\Delta LLH_peak = %.3f"%(dlogls_peak[i]))
    print ("\Delta LLH_out = %.3f"%(dlogls_out[i]))
    print()

# In[59]:

ext: py

max_stars_repo_name: parsotat/NITRATES

max_issues_repo_name: parsotat/NITRATES

max_forks_repo_name: parsotat/NITRATES

reason: import numpy,from astropy


text:
from typing import Optional

import h5py
import numpy as np

from ef.output import OutputWriter
from ef.particle_interaction_model import Model
from ef.util.serializable_h5 import tree_to_hdf5

class OutputWriterHistory(OutputWriter):
    def __init__(self, prefix: str, suffix: str):
        self.prefix: str = prefix
        self.suffix: str = suffix
        self.h5file: h5py.File = h5py.File(f"{prefix}history{suffix}", 'w')

    def __del__(self):
        self.h5file.close()

    def write(self, sim: 'Simulation', name: Optional[str] = None) -> None:
        if name is not None:
            return  # don't write fields_without_particles etc.
        if self.h5file.get('history') is None:
            self.init_file(sim, self.h5file)
        t = sim.time_grid.current_node // sim.time_grid.node_to_save
        h = self.h5file['history']
        for p in sim.particle_arrays:
            for i, id_ in enumerate(p.ids):
                h['particles/position'][id_, t] = p.positions[i]
                h['particles/momentum'][id_, t] = p.momentums[i]
                h['particles/mass'][id_] = p.mass
                h['particles/charge'][id_] = p.charge
        if sim.particle_interaction_model == Model.PIC:
            h['field/potential'][t] = sim.potential

    def init_file(self, sim: 'Simulation', h5file: h5py.File) -> None:
        h = h5file.create_group('history')
        n_particles = sum(s.initial_number_of_particles +
                          s.particles_to_generate_each_step * sim.time_grid.total_nodes
                          for s in sim.particle_sources)
        n_time = (sim.time_grid.total_nodes - 1) // sim.time_grid.node_to_save + 1
        h['time'] = np.linspace(0, sim.time_grid.total_time, n_time)
        h['particles/ids'] = np.arange(n_particles)
        h['particles/coordinates'] = [np.string_('x'), np.string_('y'), np.string_('z')]
        h.create_dataset('particles/position', (n_particles, n_time, 3))
        h['particles/position'].dims[0].label = 'id'
        h['particles/position'].dims.create_scale(h['particles/ids'], 'ids')
        h['particles/position'].dims[0].attach_scale(h['particles/ids'])
        h['particles/position'].dims[1].label = 'time'
        h['particles/position'].dims.create_scale(h['time'], 'time')
        h['particles/position'].dims[1].attach_scale(h['time'])
        h['particles/position'].dims[2].label = 'coordinates'
        h['particles/position'].dims.create_scale(h['particles/coordinates'], 'coordinates')
        h['particles/position'].dims[2].attach_scale(h['particles/coordinates'])
        h.create_dataset('particles/momentum', (n_particles, n_time, 3))
        h['particles/momentum'].dims[0].label = 'id'
        h['particles/momentum'].dims.create_scale(h['particles/ids'], 'ids')
        h['particles/momentum'].dims[0].attach_scale(h['particles/ids'])
        h['particles/momentum'].dims[1].label = 'time'
        h['particles/momentum'].dims.create_scale(h['time'], 'time')
        h['particles/momentum'].dims[1].attach_scale(h['time'])
        h['particles/momentum'].dims[2].label = 'coordinate'
        h['particles/momentum'].dims.create_scale(h['particles/coordinates'], 'coordinates')
        h['particles/momentum'].dims[2].attach_scale(h['particles/coordinates'])
        h.create_dataset('particles/mass', (n_particles,))
        h['particles/mass'].dims[0].label = 'id'
        h['particles/mass'].dims.create_scale(h['particles/ids'], 'ids')
        h['particles/mass'].dims[0].attach_scale(h['particles/ids'])
        h.create_dataset('particles/charge', (n_particles,))
        h['particles/charge'].dims[0].label = 'id'
        h['particles/charge'].dims.create_scale(h['particles/ids'], 'ids')
        h['particles/charge'].dims[0].attach_scale(h['particles/ids'])
        if sim.particle_interaction_model == Model.PIC:
            h.create_dataset('field/potential', (n_time, *sim.potential.n_nodes))
            h['field/potential'].dims[0].label = 'time'
            h['field/potential'].dims.create_scale(h['time'], 'time')
            h['field/potential'].dims[0].attach_scale(h['time'])
            for i, c in enumerate('xyz'):
                h[f'field/{c}'] = np.linspace(0, sim.electric_field.size[i], sim.electric_field.n_nodes[i])
                h['field/potential'].dims[i + 1].label = c
                h['field/potential'].dims.create_scale(h[f'field/{c}'], c)
                h['field/potential'].dims[i + 1].attach_scale(h[f'field/{c}'])
        tree_to_hdf5(sim.tree, self.h5file.create_group('simulation'))

reason: import numpy


text:
#
#  2021/06/11 Yamashita
#
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
from matplotlib.colors import ListedColormap
from .cutil import ColUtils

class val2col():
    """カラーマップの設定（グラデーション）

    Parameters:
    ----------
    tmin: float
        カラーマップの下限
    tmax: float
        カラーマップの上限
    tstep: float
        カラーマップのラベルを描く間隔
    cmap: str
        色テーブルの名前
    ----------
    """
    def __init__(self, tmin=0., tmax=20., tstep=2, cmap='jet'):
        self.tmin = tmin
        self.tmax = tmax
        self.tstep = tstep
        self.cmap = cmap
        try:
            self.cm = plt.get_cmap(self.cmap)
        except:
            cutils = ColUtils(cmap)  # 色テーブルの選択
            self.cm = cutils.get_ctable()  # 色テーブルの取得

    def conv(self, val):
        """データをカラーに変換

        Parameters:
        ----------
        val: float
            データ
        ----------
        Returns:
        ----------
        cmap
            カラーマップ
        ----------
        """
        n = (val - self.tmin) / (self.tmax - self.tmin) * self.cm.N
        n = max(min(n, self.cm.N), 0)
        return self.cm(int(n))

    def colorbar(self,
                 fig=None,
                 anchor=(0.35, 0.24),
                 size=(0.3, 0.02),
                 fmt="{f:.0f}",
                 label=True):
        """カラーバーを描く

        Parameters:
        ----------
        fig: matplotlib Figure
            プロット領域を作成した際の戻り値
        anchor: tuple(float, float)
            カラーバーの位置
        size: tuple(float, float)
            カラーバーの大きさ
        fmt: str
            カラーバーの目盛り線ラベルの書式
        label: bool
            カラーバーの目盛り線ラベルを描くかどうか
        ----------
        """
        if fig is None:
            raise Exception('fig is needed')
        ax = fig.add_axes(anchor + size)
        gradient = np.linspace(0, 1, self.cm.N)
        gradient_array = np.vstack((gradient, gradient))
        ticks = list()
        labels = list()
        ll = np.arange(self.tmin, self.tmax, self.tstep)
        for t in ll:
            ticks.append((t - self.tmin) / (self.tmax - self.tmin) * self.cm.N)
            if label:
                labels.append(fmt.format(f=t))
        # カラーバーを描く
        ax.imshow(gradient_array, aspect='auto', cmap=self.cm)
        ax.yaxis.set_major_locator(mticker.NullLocator())
        ax.yaxis.set_minor_locator(mticker.NullLocator())
        ax.set_xticks(ticks)
        if label:
            ax.set_xticklabels(labels)
        else:
            ax.xaxis.set_major_formatter(mticker.NullFormatter())
            ax.xaxis.set_minor_formatter(mticker.NullFormatter())
        #ax.set_axis_off()

    def clabel(self,
               fig=None,
               anchor=(0.34, 0.24),
               size=(0.1, 0.02),
               text=None,
               ha='right',
               va='bottom',
               fontsize=14):
        """カラーバーにラベルを付ける

        Parameters:
        ----------
        fig: matplotlib Figure
            プロット領域を作成した際の戻り値
        anchor: tuple(float, float)
            カラーバーラベルの位置
        size: tuple(float, float)
            カラーバーラベルの大きさ
        text: str
            カラーバーラベルに描く文字列
        ha: str
            カラーバーラベルの水平位置の揃え方
        va: str
            カラーバーラベルの鉛直位置の揃え方
        fontsize: int
            カラーバーラベルの文字サイズ
        ----------
        """
        if fig is None:
            raise Exception('fig is needed')
        ax = fig.add_axes(anchor + size)
        ax.text(0.0, 0.0, text, ha=ha, va=va, fontsize=fontsize)
        ax.set_axis_off()

class collevs():
    """カラーマップの設定（指定したレベル）

    Parameters:
    ----------
    clevs: list(float, float,,,)
        カラーマップの境の値
    tmax: list(float, float, float,,,)
        カラーマップの色
    ----------
    """

    def __init__(self, clevs=[], ccols=[]):
        self.clevs_min = clevs[0]
        self.clevs = clevs
        self.ccols = ccols
        if len(clevs) != len(ccols) + 1:
            print(len(clevs), len(ccols))
            raise ValueError('len(clevs) must be len(ccols)+1')

    def conv(self, val):
        """データをカラーに変換

        Parameters:
        ----------
        val: float
            データ
        ----------
        Returns:
        ----------
        cmap
            カラーマップ
        ----------
        """
        if val < float(self.clevs[0]):
            col = 'gray'
            #col = self.ccols[0]
        for n in np.arange(len(self.ccols)):
            if val >= float(self.clevs[n]) and val < float(self.clevs[n + 1]):
                col = self.ccols[n]
        return col

    def colorbar(self,
                 fig=None,
                 anchor=(0.35, 0.24),
                 size=(0.3, 0.02),
                 fontsize=11):
        """カラーバーを描く

        Parameters:
        ----------
        fig: matplotlib Figure
            プロット領域を作成した際の戻り値
        anchor: tuple(float, float)
            カラーバーの位置
        size: tuple(float, float)
            カラーバーの大きさ
        ----------
        """
        if fig is None:
            raise Exception('fig is needed')
        ax = fig.add_axes(anchor + size)
        gradient = np.linspace(0, 1, len(self.clevs) - 1)
        gradient_array = np.vstack((gradient, gradient))
        ticks = list()
        labels = list()
        ll = self.clevs[:-1]
        for n, t in enumerate(ll):
            ticks.append(n - 0.5)
            if t >= 1.:
                labels.append("{f:.0f}".format(f=t))
            else:
                labels.append("{f:.1f}".format(f=t))
        # カラーバーを描く
        cmap = ListedColormap(self.ccols)
        ax.imshow(gradient_array, aspect='auto', cmap=cmap)
        ax.yaxis.set_major_locator(mticker.NullLocator())
        ax.yaxis.set_minor_locator(mticker.NullLocator())
        ax.set_xticks(ticks)
        ax.set_xticklabels(labels, fontsize=fontsize)
        #ax.set_axis_off()


text:
from posixpath import dirname
import matplotlib.pyplot as plt
from PIL import Image
import numpy as np
import datetime
import os
import sys
import inspect

from yolov5 import model_adapter as yol
from dtlpy.ml import train_utils
import dtlpy as dl

def get_globals():
    model = dl.models.get(model_name='yolo-v5')
    snapshot = model.snapshots.get('pretrained-yolo-v5-s')
    model.snapshots.list().to_df()
    return model, snapshot

def create_sample_dataset(project):
    try:
        dataset = project.datasets.get(dataset_name='FruitImage')
    except dl.exceptions.NotFound:
        dataset = project.datasets.create('FruitImage')
        _ = dataset.items.upload(local_path='../../assets/sample_datasets/FruitImage/items/train/*',
                                 local_annotations_path='../../assets/sample_datasets/FruitImage/json/train',
                                 remote_path='train')
        _ = dataset.items.upload(local_path='../../assets/sample_datasets/FruitImage/items/test/*',
                                 local_annotations_path='../../assets/sample_datasets/FruitImage/json/test',
                                 remote_path='test')
        dataset.labels = ['apple', 'banana', 'orange', 'mixed']

    dataset.to_df()

    return dataset

def create_model_and_snapshot(project):
    try:
        model = yol.model_creation(project_name=project.name, env=args.env)
    except dl.exceptions.BadRequest:
        model = project.models.get('yolo-v5')
    snapshot = yol.snapshot_creation(model, yolo_size='small', env=args.env)
    return model, snapshot

def pretrained_inference(item_id, model, snapshot):
    adapter = model.build()
    adapter.load_from_snapshot(snapshot=snapshot)  # local_path=snapshot.bucket.local_path
    item = dl.items.get(item_id=item_id)
    image = np.asarray(Image.open(item.download()))

    # use the adapter to inference without uploading annotations
    all_annotations = adapter.predict_items([item], with_upload=False, min_score=0.1)

    item_annotations = all_annotations[0]  # predict items returns a list for each of the input items
    print(item_annotations.to_df())

    # add the items to the image
    annotated_image = item.annotations.show(image=image, thickness=5)
    plt.imshow(annotated_image)

def train_on_new_dataset(model, project, snapshot, dataset):

    partitions = {dl.SnapshotPartitionType.TRAIN: dl.Filters(field='dir', values='/train'),
                  dl.SnapshotPartitionType.VALIDATION: dl.Filters(field='dir', values='/test')}

    snapshot_name = f'trained-{dataset.name}'
    try:
        new_snapshot = model.snapshots.get(snapshot_name=snapshot_name)
    except dl.exceptions.NotFound:
        cloned_dataset = train_utils.prepare_dataset(dataset,
                                                     filters=None,
                                                     partitions=partitions)
        bucket = project.buckets.create(bucket_type=dl.BucketType.ITEM, model_name=model.name, snapshot_name=snapshot_name)
        new_snapshot = snapshot.clone(snapshot_name=snapshot_name,
                                      dataset_id=cloned_dataset.id,
                                      bucket=bucket)

    new_snapshot.configuration.update({'batch_size': 16,
                                       'start_epoch': 0,
                                       'max_epoch': 5,
                                       'data_yaml_fname': 'data.yaml',
                                       'hyp_yaml_fname': 'hyp.finetune.yaml',
                                       'id_to_label_map': {ind: label for ind, label in
                                                                         enumerate(dataset.instance_map)}
                                         })

    adapter = model.build()
    adapter.load_from_snapshot(snapshot=new_snapshot,
                               configuration=new_snapshot.configuration)

    root_path, data_path, output_path = adapter.prepare_training(root_path=os.path.join('tmp', new_snapshot.id))

    print("Training {!r} with snapshot {!r} on data {!r}".format(model.name, new_snapshot.id, data_path))
    adapter.train(data_path=data_path,
                  output_path=output_path)

    adapter.snapshot.bucket.list_content().to_df()

    adapter.save_to_snapshot(local_path=output_path,
                             replace=True)

    adapter.snapshot.bucket.list_content().to_df()

def main(args, **kwargs):
    project = dl.projects.get(args.project)

    model, snapshot = get_globals()
    # model, snapshot = create_model_and_snapshot(project)

    if args.mode == 'train':
        if not args.dataset:
            dataset = create_sample_dataset(project)
        else:
            dataset = project.datasets.get(None, args.dataset)
        dataset.to_df()
        train_on_new_dataset(model=model,
                             project=project,
                             snapshot=snapshot,
                             dataset=dataset)
    elif args.mode == 'inference':
        pretrained_inference(item_id=args.item,
                             model=model,
                             snapshot=snapshot)

if __name__ == '__main__':
    import argparse
    dl.setenv('prod')
    parser = argparse.ArgumentParser(description='local runner for model management testing')

    parser.add_argument('--env', '-e', default='prod', help='dtlpy env')
    parser.add_argument('--project', '-p', default='distillator', help='dtlpy project name', required=True)
    parser.add_argument('--dataset', '-d', default='', help='dtlpy dataset id')
    parser.add_argument('--item', '-i', default='6205097d8ea6ad0abe7b90ba', help='dtlpy single item id')
    parser.add_argument('--mode', '-m', default='inference', help='inference or training')

    args = parser.parse_args()

    main(args)
    print('Done!')

max_stars_repo_name: dataloop-ai/dtlpy-documentation

max_issues_repo_name: dataloop-ai/dtlpy-documentation

max_forks_repo_name: dataloop-ai/dtlpy-documentation


text:
#!/usr/bin/env python3
# Name: Bryan Thornlow
# Date: 2/1/2018
# compareDatabases.py

import sys
import os
import datetime
import numpy
from numpy import random
import gzip
import math
import re

##########################
##### MAIN FUNCTIONS #####
##########################

def checkClusters():
    trioTo3P = {}
    with open('filtering/data/allRelevantNodesMNKPval.txt') as f:
        for line in f:
            splitLine = (line.strip()).split('\t')
            trioTo3P[str(splitLine[0])+'_'+str(splitLine[1])+'_'+str(splitLine[2])] = splitLine[8:]

    goodTrios = {}
    with open('filtering/data/combinedCatOnlyBestWithPValsFinalReportWithInfSites.txt') as f:
        for line in f:
            splitLine = (line.strip()).split('\t')
            goodTrios[str(splitLine[0])+'_'+str(splitLine[3])+'_'+str(splitLine[6])] = False

    myOutString = ''
    bp1 = {}
    bp2 = {}
    myOKs = {28881:True,28882:True,28883:True,28280:True,28281:True,28282:True}
    counter = 0
    with open('filtering/data/combinedCatOnlyBestWithPValsFinalReportWithInfSites.txt') as f:
        for line in f:
            splitLine = (line.strip()).split('\t')
            myTrio = str(splitLine[0])+'_'+str(splitLine[3])+'_'+str(splitLine[6])
            myStart = int(splitLine[1].split(',')[0][1:])
            myEnd = int(splitLine[2].split(',')[1][:-1])
            mySeq = splitLine[16]
            mySites = toInt(splitLine[15].split(','))
            myA = []
            myB = []
            for i in range(0,len(mySites)):
                if mySeq[i] == 'A':
                    myA.append(mySites[i])
                elif mySeq[i] == 'B':
                    myB.append(mySites[i])
            if max(myA)-min(myA) > 20 and max(myB)-min(myB) > 20:
                myOutString += line.strip()+'\t'+joiner(trioTo3P[myTrio])+'\n'
                if myStart == 0 or myEnd == 29903:
                    bp1[splitLine[0]] = True
                else:
                    bp2[splitLine[0]] = True

            if str(splitLine[0])+'_'+str(splitLine[3])+'_'+str(splitLine[6]) in goodTrios:
                counter += 1
                goodTrios[str(splitLine[0])+'_'+str(splitLine[3])+'_'+str(splitLine[6])] = True

    for k in goodTrios:
        if goodTrios[k] == False:
            print(k.split('_')[0])

    print(len(bp1),len(bp2))
    print(counter)
    open('filtering/data/combinedCatOnlyBestWithPValsFinalReportWithInfSitesNoClusters.txt','w').write(myOutString)

##########################
#### HELPER FUNCTIONS ####
##########################

def toInt(myList):
    myReturn = []
    for k in myList:
        myReturn.append(int(k))
    return(myReturn)

def joiner(entry):
    newList = []
    for k in entry:
        newList.append(str(k))
    return('\t'.join(newList))

def joinerC(entry):
    newList = []
    for k in entry:
        newList.append(str(k))
    return(','.join(newList))

#########################
##### FUNCTION CALL #####
#########################

def main():
    checkClusters()

if __name__ == "__main__":
    """
    Calls main when program is run by user.
    """
    main();
    raise SystemExit

max_stars_repo_name: mrkylesmith/usher

max_issues_repo_name: mrkylesmith/usher

max_forks_repo_name: mrkylesmith/usher

reason: import numpy,from numpy


text:
import random
import scipy.io as sio

def Site(choice):
    return {
        1 : 'booth',
        2 : 'office',
        3 : 'meeting',
        4 : 'lecture'}.get(choice, 'error')    # error is default if choice not found

def get_impulse():

  Room = Site(random.randint(1,4))
  channel = str(random.randint(0,1))
  rir_no = str(random.randint(1,3))

  file_name = ['AIR_1_4/air_binaural_' + Room + '_' + channel +
              '_1_' + rir_no + '.mat'][0]

  mat =  sio.loadmat(file_name)
  Fs =  float(mat['air_info'][0]['fs'][0][0][0])
  sound = mat['h_air'][0]

  return sound, Fs

max_stars_repo_name: andimarafioti/AIAMI

reason: import scipy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov  4 09:55:30 2019

Calculate cumulative growth of wt_depth areas and marine/tidal w/ slr

@author: kbefus
"""

import os,sys
import glob
import numpy as np
import geopandas as gpd
from scipy.spatial import cKDTree as KDTree
from shapely.ops import unary_union
from shapely import speedups
speedups.enable()
res_dir = r'/mnt/data2/CloudStation'
code_dir = os.path.join(res_dir,r'ca_slr/scripts')
sys.path.insert(1,code_dir)

from cgw_model.cgw_utils import cgw_general_utils as cgu
from cgw_model.cgw_utils import cgw_feature_utils as cfu
from shapely.ops import shared_paths,snap,linemerge

#%%

research_dir_main = os.path.join(res_dir,'ca_slr')
research_dir = r'/mnt/762D83B545968C9F'
data_dir = os.path.join(research_dir,'data')

active_date = '29Oct19'
output_dir = 'outputs_fill_gdal_{}'.format(active_date)
shp_dir = os.path.join(data_dir,output_dir)

model_types = ['model_lmsl_noghb','model_mhhw_noghb']

sealevel_elevs = np.hstack([np.arange(0,2.25,.25),2.5,3.,5.])# m
Kh_vals = [0.1,1.,10.]
dirname_fmt = 'slr_{0:3.2f}_m_Kh{1:3.1f}mday'

utm10n = 3717
utm11n = 3718
cell_spacing = 10

ncrs =  {'init': 'epsg:{}'.format(utm10n)}
scrs = {'init': 'epsg:{}'.format(utm11n)}
crs_dict = {'norca':ncrs,'paca':ncrs,'sfbay':ncrs,
            'cenca':ncrs,'soca':scrs}

marine_type = 'marine/tidal'
swi_type = 'swi_footprint'
fresh_type = 'fresh_footprint'
feature_types = [fresh_type,swi_type,marine_type]

all_footprints = []
min_area = 1 # smaller than 1e2 takes all

overwrite_bool = False
save_dir = os.path.join(research_dir_main,'results','no_ghb','wt_analysis')
sv_fname = os.path.join(save_dir,'wtdepth_areas_bycounty_6Dec19.csv')

#linear_resp_bool = True
out_list = []

wt_col = 'wtdepth'
for linear_resp_bool in [False, True]:
    for model_type in model_types:
        datum_type = model_type.split('_')[1].upper()
        scenario_type = '_'.join(model_type.split('_')[1:])
        print('---------- {} --------------'.format(scenario_type))
        county_dirs = glob.glob(os.path.join(shp_dir,model_type,'shp','*'))
        county_dirs = [idir for idir in county_dirs if os.path.isdir(idir) and 'Error' not in idir]
        for Kh in Kh_vals:
            print('------------ Kh = {} ---------------'.format(Kh))
            kh_dir = 'Kh{0:3.2f}mday'.format(Kh)
            kh_dir=kh_dir.replace('.','p')

            for cdir in county_dirs:
                county_name = os.path.basename(cdir)
                print('----------- {} -------------'.format(county_name))
                county_merged = {}
                for sl in sealevel_elevs:
                    print('--- SL = {} ----'.format(sl))
                    if linear_resp_bool and sl==0:
                        # Load original, not lin, modeled output for sl=0
                        temp_fname = '{0}_{1}_slr{2:3.2f}m_Kh{3:3.2f}mday_emergent'.format(county_name,scenario_type,sl,Kh)
                        temp_fname = temp_fname.replace('.','p')
                        shp_name = os.path.join(cdir,kh_dir,'{}.shp'.format(temp_fname))
                        shp_df = gpd.read_file(shp_name)
                        out_type = '_'.join([scenario_type,'linresp'])
                    else:
                        temp_fname = '{0}_{1}_slr{2:3.2f}m_Kh{3:3.2f}mday_emergent'.format(county_name,scenario_type,sl,Kh)
                        if linear_resp_bool:
                           kh_dir2 = '_'.join(['linresponse',kh_dir])
                           temp_fname = '{}_lin'.format(temp_fname)
                           out_type = '_'.join([scenario_type,'linresp'])
                        else:
                            kh_dir2 = kh_dir
                            out_type = scenario_type

                        temp_fname = temp_fname.replace('.','p')
                        shp_name = os.path.join(cdir,kh_dir2,'{}.shp'.format(temp_fname))
                        shp_df = gpd.read_file(shp_name)

#                    unique_types = shp_df[wt_col].unique()

                    shp_df['area_km2'] = shp_df.area/1e6
                    gdf = shp_df.groupby(['wtdepth'])['area_km2'].agg('sum')
                    unique_types = gdf.index.values.unique()

                    for temp_type in unique_types:
                        # new way
                        out_list.append([out_type,Kh,county_name,sl,temp_type,gdf[temp_type]])

                        # old way
#                        temp_shp = shp_df[shp_df[wt_col]==temp_type].unary_union
#                        if temp_type not in list(county_merged.keys()):
#                            # Add key and merge features
#                            county_merged.update({temp_type:temp_shp})
#
#                        else:
#                            county_merged[temp_type] = unary_union([county_merged[temp_type],temp_shp])
#
#                        out_list.append([out_type,Kh,county_name,sl,temp_type,county_merged[temp_type].area/1e6])

            sv_cols = ['Scenario','Kh_mday','County','Sea_level_m','Category','Area_km2']
            sv_df = gpd.pd.DataFrame(out_list,columns=sv_cols)
            if os.path.isfile(sv_fname) and not overwrite_bool:
                sv_df_orig = gpd.pd.read_csv(sv_fname)
                sv_df = gpd.pd.concat([sv_df_orig,sv_df],ignore_index=True)

            sv_df.to_csv(sv_fname,index=False)

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Oct  9 11:22:22 2018

@author: michaelek
"""
from time import time
import numpy as np
import pandas as pd
import xarray as xr
from hydrointerp.io.raster import save_geotiff
from hydrointerp.util import grp_ts_agg
from nzmetservice import select_bounds, to_df
from scipy import ndimage

pd.options.display.max_columns = 10

#####################################
### Parameters

min_lat = -47
max_lat = -40
min_lon = 166
max_lon = 175
nc1 = r'N:\met_service\forecasts\wrf_hourly_precip_nz4kmN-NCEP_2018090918.nc'
nc2 = r'N:\nasa\precip\gpm_3IMERGHHL\gpm_3IMERGHHL_v07_20190101-20190108.nc4'

from_crs = 4326
to_crs = 2193
grid_res = 1000
interp_fun = 'cubic'

time_col = 'time'
x_col = 'longitude'
y_col = 'latitude'
data_col = 'precip_rate'
interp_fun = 'cubic'
digits = 2

ts_resample_code = '4H'

point_path = r'N:\met_service\point_test1.shp'
point_site_col = 'site_id'
site_test = 5

####################################
### Import

ms1 = select_bounds(nc1, min_lat, max_lat, min_lon, max_lon)

ms_df = to_df(ms1, True).dropna().reset_index()
ms_df1 = ms_df[(ms_df.time <= '2018-09-24 12:00') & (ms_df.time >= '2018-09-24')]

### Resample in time
ms_df2 = grp_ts_agg(ms_df1, ['longitude', 'latitude'], 'time', ts_resample_code, 'left', 'right')[data_col].sum().reset_index()

### Interp

new_df = interp_to_grid(ms_df2, 'time', 'longitude', 'latitude', 'precip_rate', grid_res, from_crs, to_crs, interp_fun)

new_ds = interp_to_grid(ms_df2, 'time', 'longitude', 'latitude', 'precip_rate', grid_res, from_crs, to_crs, interp_fun, output='xarray')

new_points = interp_to_points(ms_df2, 'time', 'longitude', 'latitude', 'precip_rate', point_shp, point_site_col, from_crs)

###################################
### Testing

nc2 = r'N:\nasa\precip\gpm_3IMERGHHL_v07_20190101-20190106.nc4'
time_name = 'time'
x_name = 'lon'
y_name = 'lat'
data_name = 'precipitationCal'
from_crs = 4326
to_crs = 2193
grid_res = 1000
bbox=None
order=1
extrapolation='constant'
cval=np.nan
digits = 2
min_lat = -48
max_lat = -41
min_lon = 170
max_lon = 178
min_val=0
method='linear'
point_data = r'N:\met_service\point_test1.shp'

grid1 = xr.open_dataset(nc2)

grid = grid1.where((grid1[y_name] >= min_lat) & (grid1[y_name] <= max_lat) & (grid1[x_name] >= min_lon) & (grid1[x_name] <= max_lon), drop=True).copy()
grid1.close()

grid2d = grid.isel(time=280)[data_name]

grid2d.plot.pcolormesh(x='lon', y='lat')

start1 = time()
interp1 = grid_to_grid(grid, time_name, x_name, y_name, data_name, grid_res, from_crs, to_crs, bbox, order, extrapolation, cval, digits, min_val)
end1 = time()
end1 - start1

output2d = interp1.isel(time=280)[data_name]

output2d.plot.pcolormesh(x='x', y='y')

df = grid.to_dataframe().reset_index()

x = da1[x_name].values
y = da1[y_name].values

xy_orig_pts = np.dstack(np.meshgrid(x, y)).reshape(-1, 2)

da2 = da1.transpose(x_name, y_name, time_name)
ar1 = da2.values
ar1[:2] = np.nan

np.nan_to_num(ar1, False)

points_x, points_y = np.broadcast_arrays(xinterp.reshape(-1,1), yinterp)
coord = np.vstack((points_x.flatten()*(len(xgrid)-1), points_y.flatten()*(len(ygrid)-1)))

grouped = df1.groupby([x_name, y_name, time_name])[data_name].mean()

grouped = df1.set_index([x_name, y_name, time_name])[data_name]

# create an empty array of NaN of the right dimensions
shape = tuple(len(i) for i in grouped.index.levels)
arr = np.full(shape, np.nan)

# fill it using Numpy's advanced indexing
arr[tuple(grouped.index.labels)] = grouped.values.flat

df = pd.DataFrame({'x': [1, 2, 1, 3, 1, 2, 3, 1, 2],
                   'y': [1, 1, 2, 2, 1, 1, 1, 2, 2],
                   'z': [1, 1, 1, 1, 2, 2, 2, 2, 2],
                   'value': [1, 2, 3, 4, 5, 6, 7, 8, 9]})

grouped = df.groupby(['z', 'y', 'x'])['value'].mean()

a = np.arange(12.).reshape((4, 3))

dx, dy = 0.4, 0.4
xmax, ymax = 2, 4
x = np.arange(-xmax, xmax, dx)
y = np.arange(-ymax, ymax, dy)
X, Y = np.meshgrid(x, y)
Z = np.exp(-(2*X)**2 - (Y/2)**2)

interp2b = grid_to_points(grid, time_name, x_name, y_name, data_name, point_data, from_crs, to_crs, order, extrapolation, fill_val, digits, min_val)

interp2

d2 = df[df.time == pd.Timestamp('2019-01-06T20:29:59.999000')].copy()

interp2c = points_to_grid(d2, time_name, x_name, y_name, data_name, grid_res, from_crs, to_crs, bbox, method, fill_val, digits, min_val)

output2d = interp2c.isel(time=0)[data_name]

output2d.plot.pcolormesh(x='x', y='y')

grid1 = xr.open_dataset(nc1)
grid = grid1.copy()
grid1.close()

nc2 = r'N:\nasa\precip\gpm_3IMERGHH\gpm_3IMERGHH_v07_20170101-20171231.nc4'
nc2 = r'N:\nasa\precip\gpm_3IMERGHH\gpm_3IMERGHH_v07_20180101-20180701.nc4'

grid1 = xr.open_dataset(nc2)
da = grid1[data_name].copy()
grid1.close()

da1 = da * 0.5

da2 = da1.resample(time='D', closed='right', label='left').sum('time')

max_day = da2.groupby('time').sum().to_dataframe().sort_values('precipitationCal')

da3 = da2.loc['2018-01-15':'2018-01-22']

#grid2d.plot.pcolormesh(x='lon', y='lat')

plt1 = da3.plot(x='lon', y='lat', col='time', col_wrap=4)

save1 = r'E:\ecan\git\hydrointerp\hydrointerp\datasets\nasa_gpm_2017-07-20.nc'
save2 = r'E:\ecan\git\hydrointerp\hydrointerp\datasets\nasa_gpm_2017-07-20.tif'

da4 = da1.loc['2017-07-20':'2017-07-21 00:00']
da5 = da4.resample(time='H', closed='right', label='right').sum('time')
da5.attrs = grid1[data_name].attrs

da5.to_dataset().to_netcdf(save1)

#da5 = xr.open_dataset(save1)

df5 = da5.resample(time='D', closed='right', label='left').sum('time').to_dataframe().reset_index()

save_geotiff(df5, 4326, 'precipitationCal', 'lon', 'lat', export_path=save2)

a = np.arange(50, step=2).reshape((5,5))
a[2,2] = 2

ndimage.gaussian_filter(a, sigma=1, order=0)
ndimage.gaussian_gradient_magnitude(a, sigma=1)

ndimage.median_filter(a, size=3)
ndimage.percentile_filter(a, percentile=100, size=2)
ndimage.uniform_filter(a, size=3)

date0 = '2018-01-17'

date1 = pd.Timestamp(date0)
date2 = date1 + pd.DateOffset(days=1)

da4 = da1.loc[date1:date2]
da5 = da4.resample(time='D', closed='right', label='left').sum('time')
da6 = da5.copy()
da6['time'] = da6['time'].to_series() + pd.DateOffset(hours=1)

#da6.data = ndimage.percentile_filter(da5, percentile=100, size=2)
da6.data = ndimage.median_filter(da5, size=2)

da7 = da5.combine_first(da6)

plt1 = da7.plot(x='lon', y='lat', col='time', col_wrap=2)

save2 = r'E:\ecan\git\hydrointerp\hydrointerp\datasets\nasa_gpm_2018-01-18.tif'
save_geotiff(da6.to_dataframe().reset_index(), 4326, 'precipitationCal', 'lon', 'lat', export_path=save2)

from nasadap import Nasa

product = '3IMERGDF'
product = '3IMERGHH'
#product2c = '3IMERGHHE'
product = '3B42'
username = 'Dryden'
password = 'NasaData4me'
from_date = '2018-04-15'
to_date = '2018-04-17'
mission = 'gpm'
mission = 'trmm'
dataset_type = ['precipitationCal', 'randomError', 'precipitationQualityIndex']
dataset_type = ['precipitation']
min_lat=-49
max_lat=-33
min_lon=165
max_lon=180
cache_dir = r'\\fs02\GroundWaterMetData$\nasa\cache\nz'

ge = Nasa(username, password, mission, cache_dir)
ds2 = ge.get_data(product, dataset_type, from_date, to_date, min_lat, max_lat, min_lon, max_lon, 60)
ds3 = ds2.copy().load()

ds3['time'] = ds2.time.to_index() + pd.DateOffset(hours=12)
del ds2

da5 = ds3['precipitationCal'].resample(time='D', closed='right', label='left').sum('time').isel(time=slice(1, 4))
ea1 = ds3['randomError'].resample(time='D', closed='right', label='left').mean('time').isel(time=slice(1, 4))
qa1 = ds3['precipitationQualityIndex'].resample(time='D', closed='right', label='left').mean('time').isel(time=slice(1, 4))

ea2 = ea1/da5

plt1 = da5.plot(x='lon', y='lat', col='time', col_wrap=3)

da6 = da5.where(qa1 > 0.4).fillna(400)
da6.name = 'precip'

plt2 = da6.plot(x='lon', y='lat', col='time', col_wrap=3)

bad_spot1 = [168.75, -44.45]

day1 = '2018-01-17'
qa1.sel(time='2018-01-17', lat=-44.45, lon=168.75)
qa1.loc['2018-01-17'].sel(lon=168.75)

ea1.loc['2018-01-17'].sel(lon=168.75)

(qa1 < 0.26).sum(['lat', 'lon'])

qa2 = da5.where(qa1 > 0.26)
qa2.name = 'precip'

df1 = qa2.to_dataframe().reset_index()

### A 'precipitationQualityIndex' of > 0.4 is the ideal setting for GPM data

da5 = ds3['precipitation'].resample(time='D', closed='right', label='left').sum('time').isel(time=slice(1, 4))

plt1 = da5.plot(x='lon', y='lat', col='time', col_wrap=3)

da5 = ds3['precipitationCal'].resample(time='D', closed='right', label='left').sum('time').isel(time=slice(1, 4))
qa1 = ds3['precipitationQualityIndex'].resample(time='D', closed='right', label='left').mean('time').isel(time=slice(1, 4))

da6 = da5.where(qa1 > 0.4)
da6.name = 'precip'

ds7 = grid_interp_na(da6.to_dataset(), time_name, x_name, y_name, 'precip', min_val=0)

#da7 = ds7['precipitationCal'].resample(time='D', closed='right', label='left').sum('time').isel(time=slice(1, 4))

plt1 = da5.plot(x='lon', y='lat', col='time', col_wrap=3)
plt2 = ds7['precip'].plot(x='lon', y='lat', col='time', col_wrap=3)

max_stars_repo_name: mullenkamp/hydrointerp

max_issues_repo_name: mullenkamp/hydrointerp

max_forks_repo_name: mullenkamp/hydrointerp


text:
from numpy import ones
from data import generateFakeSamples, generateLatentPoints, getRealSamples
from generator import createGenerator
from discriminator import createDiscriminator
from keras.models import Sequential
from matplotlib import pyplot

def createGan(discriminator: Sequential, generator: Sequential):
    discriminator.trainable = False
    model = Sequential()
    model.add(generator)
    model.add(discriminator)
    model.compile(loss='binary_crossentropy', optimizer='adam')
    return model

def evaluatePerformance(epoch, discriminator: Sequential, generator: Sequential, latentDim, n=100):
    xReal, yReal = getRealSamples(n)
    _, accReal = discriminator.evaluate(xReal, yReal, verbose=0)

    xFake, yFake = generateFakeSamples(generator, latentDim, n)
    _, accFake = discriminator.evaluate(xFake, yFake, verbose=0)

    print(epoch, accReal, accFake)

    pyplot.scatter(xReal[:, 0], xReal[:, 1], color='red')
    pyplot.scatter(xFake[:, 0], xFake[:, 1], color='blue')
    filename = 'eval/generated_plot_e%03d.png' % (epoch+1)
    pyplot.savefig(filename)
    pyplot.close()

def train(discriminator:Sequential, generator:Sequential, gan:Sequential, latentDim, numEpochs=10000, numBatch=128, evalFreq=2000):
    halfBatch = int(numBatch / 2)
    for i in range(numEpochs):
        xReal, yReal = getRealSamples(halfBatch)
        discriminator.train_on_batch(xReal, yReal)

        xFake, yFake = generateFakeSamples(generator, latentDim, halfBatch)
        discriminator.train_on_batch(xFake, yFake)

        xGan = generateLatentPoints(latentDim, halfBatch)
        yGan = ones((halfBatch, 1))
        gan.train_on_batch(xGan, yGan)

        if (i+1) % evalFreq == 0:
            evaluatePerformance(i, discriminator, generator, latentDim)

if __name__ == '__main__':
    discriminator = createDiscriminator()
    generator = createGenerator()
    gan = createGan(discriminator, generator)
    gan.summary()

reason: from numpy


text:
import io
import pathlib

import h5py
import PIL.Image
import numpy as np
from torch.utils.data import Dataset
from tqdm import tqdm

# For high performance PIL operations use
# pillow-simd instead of pillow
# https://python-pillow.org/pillow-perf/
# https://github.com/uploadcare/pillow-simd
import pkg_resources
try:
    pkg_resources.get_distribution("Pillow-SIMD")
except pkg_resources.DistributionNotFound:
    import warnings
    warnings.warn("\033[93mPillow-SIMD not installed. PIL operations could be faster using SIMD/AVX instructions https://github.com/uploadcare/pillow-simd\033[0m")

def make_dataset(hf):

    samples = []
    for label, class_ in enumerate(tqdm(sorted(hf), leave=False)):

        for image in sorted(hf[class_]):
            samples.append( (hf[class_][image].name, label) )

    return samples

def class_to_idx(hf):

    class_to_idx = {}
    for label, class_ in enumerate(sorted(hf)):
        class_to_idx[class_] = label

    return class_to_idx

def array_loader(array):
    return PIL.Image.fromarray(array)

def binary_loader(array):
    img_bytes = io.BytesIO(array)
    return PIL.Image.open(img_bytes).convert('RGB')

mode_to_loader = {
    "binary": binary_loader,
    "array": array_loader,
}

class ImageHDF5Dataset(Dataset):

    def __init__(self, h5file,
                 transform=None,  target_transform=None):
        super(ImageHDF5Dataset, self).__init__()

        self.h5file = pathlib.Path(h5file)
        indexfile = self.h5file.with_suffix('.idx.npy')

        with h5py.File(h5file, 'r') as hf:
            # TODO remove next two lines after adding mode to existing files
            if 'mode' not in hf.attrs:
                self.loader = array_loader
            else:
                self.loader = mode_to_loader[hf.attrs['mode']]
        self.transform = transform
        self.target_transform = target_transform
        # Compile image & label list
        if indexfile.exists():
            self.samples = [(x, int(y)) for x, y in np.load(indexfile)]
        else:
            print("Compiling dataset indices")
            with h5py.File(self.h5file, 'r') as hf:
                self.samples = make_dataset(hf)
                np.save(indexfile, self.samples)
            # self.class_to_idx = class_to_idx(hf)

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

    def __getitem__(self, index):
        # File descriptor must be opened here because h5 does
        # support processing inheriting the file descriptor
        # https://gist.github.com/bkj/f448025fdef08c0609029489fa26ea2a
        path, target = self.samples[index]
        with h5py.File(self.h5file, 'r') as hf:
            sample = self.loader(hf[path][()])
        if self.transform is not None:
            sample = self.transform(sample)
        if self.target_transform is not None:
            target = self.target_transform(target)

        return sample, target

# class ImageHDF5Dataset(Dataset):

#     def __init__(self, h5file, loader=PIL_loader,
#                  transform=None,  target_transform=None):
#         super(ImageHDF5Dataset, self).__init__()

#         h5file = pathlib.Path(h5file)
#         indexfile = h5file.with_suffix('.idx.npy')

#         self.hf = h5py.File(h5file, 'r')
#         self.loader = loader
#         self.transform = transform
#         self.target_transform = target_transform
#         # Compile image & label list
#         if indexfile.exists():
#             self.samples = [(x, int(y)) for x, y in np.load(indexfile)]
#         else:
#             print("Compiling dataset indices")
#             self.samples = make_dataset(self.hf)
#             np.save(indexfile, self.samples)
#         self.class_to_idx = class_to_idx(self.hf)

#     def __len__(self):
#         return len(self.samples)

#     def __getitem__(self, index):
#         path, target = self.samples[index]
#         sample = self.loader(self.hf[path].value)
#         if self.transform is not None:
#             sample = self.transform(sample)
#         if self.target_transform is not None:
#             target = self.target_transform(target)

#         return sample, target

#     def __del__(self):
#         if hasattr(self, 'hf'):
#             self.hf.close()


text:
#!/usr/bin/env python3
# Author:   Michael E. Rose <michael.ernst.rose@gmail.com>
"""Combines faculty information for faculty members
with Scopus IDs and PHD.
"""

from collections import Counter
from glob import glob
from math import log10
from os.path import basename, splitext

import pandas as pd
from numpy import nan

SOURCE_FOLDER = "./source_files/"
SCOPUS_FILE = "./mapping_files/persons.csv"
MAPPING_FILE = "./mapping_files/institutions.csv"
TARGET_FILE = "./hasselback.csv"

def complete_year(y):
    """Complete a year by adding missing century digits."""
    # Keep the following try-except until all numbers are correct
    try:
        y = int(y)
    except ValueError:
        return None
    try:
        digits = int(log10(y)) + 1
    except ValueError:  # y is 0
        return str(2000)
    if digits == 2:
        return str(1900 + y)
    elif digits == 1:
        return str(2000 + y)
    elif digits == 3:
        return str(1000 + y)
    else:
        return str(y)

def read_hasselback_file(fname, scopus, aff_map, degrees=('PHD',)):
    """Read and compile file, subsetted to relevant degress only."""
    year = splitext(basename(fname))[0].split("-")[0]
    df = pd.read_csv(fname).dropna(subset=["dep"])
    df = df[df['degree'].isin(degrees)]
    # Aggregate entries
    df['grad_year'] = df['grad_year'].apply(complete_year)
    df['school'] = df['school'].apply(lambda x: aff_map.get(x, x))
    id_cols = ['name', 'school', 'grad_year']
    df['ID'] = df[id_cols].fillna("").apply(lambda l: ";".join(l), axis=1)
    df['ID'] = df['ID'].str.replace(";;", ";").str.replace("nan", "").str.strip(";")
    # Aggregate departments
    df['dep'] = df['dep'].apply(lambda x: aff_map.get(x))
    before = df.shape[0]
    df = df.dropna(subset=["dep"]).set_index('ID').drop('grad_year', axis=1)
    after = df.shape[0]
    print(f">>> {year}: Dropping {before-after} names because of missing "
          "department data")
    # Merge with Scopus
    df = df.merge(scopus, "left", left_index=True, right_index=True)
    identified = df[~df['scopus_id'].isnull()].copy()
    identified = identified.set_index("scopus_id").add_prefix(year + "_")
    unidentified = df[df['scopus_id'].isnull()][["dep"]].copy()
    return identified, unidentified

def main():
    # Read in
    scopus = pd.read_csv(SCOPUS_FILE, index_col=0)
    mapping = pd.read_csv(MAPPING_FILE, index_col=0).dropna(subset=['our_name'])
    mapping = mapping["our_name"].to_dict()
    files = glob(SOURCE_FOLDER + "*.csv")
    df, unidentified = read_hasselback_file(files[0], scopus[["scopus_id"]], mapping)
    for file in files[1:]:
        new, new_un = read_hasselback_file(file, scopus[["scopus_id"]], mapping)
        df = df.merge(new, "outer", left_index=True, right_index=True)
        unidentified = unidentified.append(new_un)

    # Add time-invariant information
    scopus = scopus.drop_duplicates(subset="scopus_id")
    scopus = scopus.set_index("scopus_id")
    df = df.merge(scopus[["scopus_name", 'grad_year']], "left",
                  left_index=True, right_index=True)

    # Maintenance file
    unidentified = unidentified[~unidentified.index.duplicated()]
    unidentified.to_csv("./mapping_files/unmapped.csv", index_label="faculty")
    print(f">>> {unidentified.shape[0]:,} individuals from "
          f"{unidentified['dep'].nunique():,} different universities without ID")

    # Write out
    df.index = df.index.astype(int)
    df.to_csv(TARGET_FILE, index_label="scopus_id")

    # Statistics
    dep = set()
    dep_cols = [c for c in df.columns if c.endswith("dep")]
    for dep_col in dep_cols:
        dep.update(set(df[dep_col]))
    print(f">>> {df.shape[0]:,} individuals from {len(dep):,} different "
          "universities with ID")
    all_deps = dep.union(set(unidentified["dep"]))
    stats = {"N_of_Hasselback_fac_scopus": f"{df.shape[0]:,}",
             "N_of_Hasselback_dep_scopus": len(dep),
             "N_of_Hasselback_fac": f"{df.shape[0] + unidentified.shape[0]:,}",
             "N_of_Hasselback_dep": len(all_deps)}
    print(stats)

if __name__ == '__main__':
    main()

max_stars_repo_name: Michael-E-Rose/HasselbackFacultyRoster

max_issues_repo_name: Michael-E-Rose/HasselbackFacultyRoster

max_forks_repo_name: Michael-E-Rose/HasselbackFacultyRoster


text:
# -*- coding: utf-8 -*-
"""Numpy dataset"""

def _numpy_generator(data_dict, subset):
    for i in subset:
        datum = {k: v[i] for k, v in data_dict.items()}
        yield datum

def load_numpy(dataset, splits=None, shuffle=True, seed=0):
    """Loads dataset from numpy array

    It is assumed that the numpy arrays contain structures with the
    same number of atoms, and stacked (with the same axis being the
    index of the structure).

    Args:
        dataset: a dictionary of numpy arrays
        splits (dict): key-val pairs specifying the ratio of subsets
        shuffle (bool): shuffle the dataset (only used when splitting)
        seed (int): random seed for shuffling
    """
    import numpy as np
    import tensorflow as tf
    from pinn.io.base import split_list

    dtypes = {
        k: ('int32' if v.dtype.kind in ('u', 'i') else tf.keras.backend.floatx())
        for k, v in dataset.items()}
    shapes = {k: v.shape[1:] for k, v in dataset.items()}

    def generator_fn(subset): return tf.data.Dataset.from_generator(
        lambda: _numpy_generator(dataset, subset), dtypes, shapes)

    indices = list(range(dataset['elems'].shape[0]))
    if splits is None:
        return generator_fn(indices)
    else:
        subsets = split_list(indices, splits=splits, shuffle=shuffle, seed=seed)
        splitted = {k: generator_fn(v) for k,v in subsets.items()}
        return splitted

max_stars_repo_name: yqshao/PiNN

max_issues_repo_name: yqshao/pinn

max_forks_repo_name: yqshao/pinn


text:
#!/usr/bin/env python
# -*- coding:utf-8 -*-
#
# @Author: José Sánchez-Gallego (gallegoj@uw.edu)
# @Date: 2018-07-13
# @Filename: dapall.py
# @License: BSD 3-clause (http://www.opensource.org/licenses/BSD-3-Clause)
#
# @Last modified by: José Sánchez-Gallego (gallegoj@uw.edu)
# @Last modified time: 2018-11-14 12:03:58

import distutils
import os

import astropy.io.fits

import marvin
from marvin.core.exceptions import MarvinError
from marvin.utils.general import get_dapall_file, map_dapall

__all__ = ['DAPallMixIn']

class DAPallMixIn(object):
    """A mixin that provides access to DAPall paremeters.

    Must be used in combination with `.MarvinToolsClass` and initialised
    before `~.DAPallMixIn.dapall` can be called.

    `DAPallMixIn` uses the `.MarvinToolsClass.data_origin` of the object to
    determine how to obtain the DAPall information. However, if the object
    contains a ``dapall`` attribute with the path to a DAPall file, that file
    will be used.

    """

    __min_dapall_version__ = distutils.version.StrictVersion('2.1.0')

    @property
    def dapall(self):
        """Returns the contents of the DAPall data for this target."""

        if (not self._dapver or
                distutils.version.StrictVersion(self._dapver) < self.__min_dapall_version__):
            raise MarvinError('DAPall is not available for versions before MPL-6.')

        if hasattr(self, '_dapall') and self._dapall is not None:
            return self._dapall

        if self.data_origin == 'file':
            try:
                dapall_data = self._get_dapall_from_file()
            except IOError:
                marvin.log.debug('cannot find DAPall file. Trying remote request.')
                dapall_data = self._get_from_api()
        elif self.data_origin == 'db':
            dapall_data = self._get_dapall_from_db()
        else:
            dapall_data = self._get_dapall_from_api()

        self._dapall = dapall_data

        return self._dapall

    def _get_dapall_from_file(self):
        """Uses DAPAll file to retrieve information."""

        daptype = self.bintype.name + '-' + self.template.name

        dapall_path = get_dapall_file(self._drpver, self._dapver)

        assert dapall_path is not None, 'cannot build DAPall file.'

        if not os.path.exists(dapall_path):
            raise MarvinError('cannot find DAPall file in the system.')

        dapall_hdu = astropy.io.fits.open(dapall_path)

        header = dapall_hdu[0].header
        dapall_table = dapall_hdu[-1].data

        dapall_row = dapall_table[(dapall_table['PLATEIFU'] == self.plateifu) &
                                  (dapall_table['DAPTYPE'] == daptype)]

        assert len(dapall_row) == 1, 'cannot find matching row in DAPall.'

        return map_dapall(header, dapall_row[0])

    def _get_dapall_from_db(self):
        """Uses the DB to retrieve the DAPAll data."""

        dapall_data = {}

        daptype = self.bintype.name + '-' + self.template.name

        mdb = marvin.marvindb

        if not mdb.isdbconnected:
            raise MarvinError('No DB connected')

        datadb = mdb.datadb
        dapdb = mdb.dapdb

        dapall_row = mdb.session.query(dapdb.DapAll).join(
            dapdb.File, datadb.PipelineInfo, datadb.PipelineVersion).filter(
                mdb.datadb.PipelineVersion.version == self._dapver,
                dapdb.DapAll.plateifu == self.plateifu,
                dapdb.DapAll.daptype == daptype).first()

        if dapall_row is None:
            raise MarvinError('cannot find a DAPall match for this target in the DB.')

        for col in dapall_row.__table__.columns.keys():
            if col != 'pk' and '_pk' not in col:
                dapall_data[col] = getattr(dapall_row, col)

        return dapall_data

    def _get_dapall_from_api(self):
        """Uses the API to retrieve the DAPall data."""

        url = marvin.config.urlmap['api']['dapall']['url']

        url_full = url.format(name=self.plateifu,
                              bintype=self.bintype.name,
                              template=self.template.name)

        try:
            response = self._toolInteraction(url_full)
        except Exception as ee:
            raise MarvinError('found a problem while getting DAPall: {0}'.format(str(ee)))

        if response.results['error'] is not None:
            raise MarvinError('found a problem while getting DAPall: {}'
                              .format(str(response.results['error'])))

        data = response.getData()

        return data['dapall_data']

max_stars_repo_name: margudo/marvin

max_issues_repo_name: margudo/marvin

max_forks_repo_name: margudo/marvin

reason: import astropy


text:
r"""
Incompressible Stokes flow with Navier (slip) boundary conditions, flow driven
by a moving wall and a small diffusion for stabilization.

This example demonstrates the use of `no-penetration` boundary conditions as
well as `edge direction` boundary conditions together with Navier or slip
boundary conditions.

Find :math:`\ul{u}`, :math:`p` such that:

.. math::
    \int_{\Omega} \nu\ \nabla \ul{v} : \nabla \ul{u}
    - \int_{\Omega} p\ \nabla \cdot \ul{v}
    + \int_{\Gamma_1} \beta \ul{v} \cdot (\ul{u} - \ul{u}_d)
    + \int_{\Gamma_2} \beta \ul{v} \cdot \ul{u}
    = 0
    \;, \quad \forall \ul{v} \;,

    \int_{\Omega} \mu \nabla q \cdot \nabla p
    + \int_{\Omega} q\ \nabla \cdot \ul{u}
    = 0
    \;, \quad \forall q \;,

where :math:`\nu` is the fluid viscosity, :math:`\beta` is the slip
coefficient, :math:`\mu` is the (small) numerical diffusion coefficient,
:math:`\Gamma_1` is the top wall that moves with the given driving velocity
:math:`\ul{u}_d` and :math:`\Gamma_2` are the remaining walls. The Navier
conditions are in effect on both :math:`\Gamma_1`, :math:`\Gamma_2` and are
expressed by the corresponding integrals in the equations above.

The `no-penetration` boundary conditions are applied on :math:`\Gamma_1`,
:math:`\Gamma_2`, except the vertices of the block edges, where the `edge
direction` boundary conditions are applied. Optionally, Dirichlet boundary
conditions can be applied on the inlet, see the code below.

The mesh is created by ``gen_block_mesh()`` function - try different mesh
dimensions and resolutions below. For large meshes use the ``'ls_i'`` linear
solver - PETSc + petsc4py is needed in that case.
"""
import numpy as nm

from sfepy.discrete.fem.meshio import UserMeshIO
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.homogenization.utils import define_box_regions

# Mesh dimensions.
dims = nm.array([3, 1, 0.5])

# Mesh resolution: increase to improve accuracy.
shape = [11, 15, 15]

def mesh_hook(mesh, mode):
    """
    Generate the block mesh.
    """
    if mode == 'read':
        mesh = gen_block_mesh(dims, shape, [0, 0, 0], name='user_block',
                              verbose=False)
        return mesh

    elif mode == 'write':
        pass

filename_mesh = UserMeshIO(mesh_hook)

regions = define_box_regions(3, 0.5 * dims)
regions.update({
    'Omega' : 'all',
    'Edges_v' : ("""(r.Near *v r.Bottom) +v
                    (r.Bottom *v r.Far) +v
                    (r.Far *v r.Top) +v
                    (r.Top *v r.Near)""", 'edge'),
    'Gamma1_f' : ('copy r.Top', 'face'),
    'Gamma2_f' : ('r.Near +v r.Bottom +v r.Far', 'face'),
    'Gamma_f' : ('r.Gamma1_f +v r.Gamma2_f', 'face'),
    'Gamma_v' : ('r.Gamma_f -v r.Edges_v', 'face'),
    'Inlet_f' : ('r.Left -v r.Gamma_f', 'face'),
})

fields = {
    'velocity' : ('real', 3, 'Omega', 1),
    'pressure' : ('real', 1, 'Omega', 1),
}

def get_u_d(ts, coors, region=None):
    """
    Given stator velocity.
    """
    out = nm.zeros_like(coors)
    out[:] = [1.0, 1.0, 0.0]

    return out

functions = {
    'get_u_d' : (get_u_d,),
}

variables = {
    'u' : ('unknown field', 'velocity', 0),
    'v' : ('test field',    'velocity', 'u'),
    'u_d' : ('parameter field', 'velocity',
             {'setter' : 'get_u_d'}),
    'p' : ('unknown field', 'pressure', 1),
    'q' : ('test field',    'pressure', 'p'),
}

# Try setting the inlet velocity by un-commenting the ebcs.
ebcs = {
    ## 'inlet' : ('Inlet_f', {'u.0' : 1.0, 'u.[1, 2]' : 0.0}),
}

lcbcs = {
    'walls' : ('Gamma_v', {'u.all' : 'no_penetration'},
               'normals_Gamma.vtk'),
    'edges' : ('Edges_v', {'u.all' : 'edge_direction'},
               'edges_Edges.vtk'),
}

materials = {
    'm' : ({
        'nu' : 1e-3,
        'beta' : 1e-2,
        'mu' : 1e-10,
    },),
}

equations = {
    'balance' :
    """dw_div_grad.5.Omega(m.nu, v, u)
     - dw_stokes.5.Omega(v, p)
     + dw_surface_dot.5.Gamma1_f(m.beta, v, u)
     + dw_surface_dot.5.Gamma2_f(m.beta, v, u)
     =
     + dw_surface_dot.5.Gamma1_f(m.beta, v, u_d)""",
    'incompressibility' :
    """dw_laplace.5.Omega(m.mu, q, p)
     + dw_stokes.5.Omega(u, q) = 0""",
}

solvers = {
    'ls_d' : ('ls.scipy_direct', {}),
    ## 'ls_i' : ('ls.petsc', {
    ##     'method' : 'bcgsl', # ksp_type
    ##     'precond' : 'ilu', # pc_type
    ##     'eps_a' : 0.0, # abstol
    ##     'eps_r' : 1e-12, # rtol
    ##     'eps_d' : 1e10, # Divergence tolerance.
    ##     'i_max' : 2500, # maxits
    ## }),
    'newton' : ('nls.newton', {
        'i_max' : 1,
        'eps_a'      : 1e-10,
    }),
}

options = {
    'nls' : 'newton',
    'ls' : 'ls_d',
}

max_stars_repo_name: vondrejc/sfepy

max_issues_repo_name: vondrejc/sfepy

max_forks_repo_name: vondrejc/sfepy


text:
'''
output : Gaussian Blurring (Image Smoothing)

Syntax :  cv2.GaussianBlur(src, ksize, sigmaX[, dst[, sigmaY[, borderType]]]) ->dst

Parameters :
src	    input image; the image can have any number of channels, which are processed independently, but the depth should be CV_8U, CV_16U, CV_16S, CV_32F or CV_64F.
dst	    output image of the same size and type as src.
ksize	Gaussian kernel size. ksize.width and ksize.height can differ but they both must be positive and odd. Or, they can be zero's and then they are computed from sigma.
sigmaX	Gaussian kernel standard deviation in X direction.
sigmaY	Gaussian kernel standard deviation in Y direction; if sigmaY is zero, it is set to be equal to sigmaX, if both sigmas are zeros, they are computed from ksize.width and ksize.height, respectively (see getGaussianKernel for details); to fully control the result regardless of possible future modifications of all this semantics, it is recommended to specify all of ksize, sigmaX, and sigmaY.
borderType	pixel extrapolation method, see BorderTypes. BORDER_WRAP is not supported.

Reference : https://docs.opencv.org/master/d4/d86/group__imgproc__filter.html#gaabe8c836e97159a9193fb0b11ac52cf1
'''

import cv2
import numpy as np
import matplotlib.pyplot as plt

"""
We’ll see in the code below, when the size of our kernel increases so will the amount of blurring that is applied to our output image.
However, the blurring will appear to be more “natural” and will preserve edges in our image better than simple average smoothing:
"""

img = cv2.imread('../Images and videos/Gaussian-noise.jpg')

gaussian_3x3 = cv2.GaussianBlur(img, (3,3), sigmaX=0, sigmaY=0)
gaussian_5x5 = cv2.GaussianBlur(img, (5,5), sigmaX=0, sigmaY=0)
gaussian_7x7 = cv2.GaussianBlur(img, (7,7), sigmaX=0, sigmaY=0)
gaussian_9x9 = cv2.GaussianBlur(img, (9,9), sigmaX=0, sigmaY=0)
gaussian_11x11 = cv2.GaussianBlur(img, (11,11), sigmaX=0, sigmaY=0)
images = [img, gaussian_3x3, gaussian_5x5, gaussian_7x7, gaussian_9x9, gaussian_11x11]
titles = ['Original', 'Gaussian Smoothening - 3x3','Gaussian Smoothening - 5x5','Gaussian Smoothening - 7x7','Gaussian Smoothening - 9x9','Gaussian Smoothening - 11x11']

for i in range(0,6):
    plt.subplot(2,3,(i+1))
    plt.imshow(images[i],'gray')
    plt.title(titles[i])
    plt.axis('off')
plt.show()

# Callback Function which helps to get the position of trackbar
def nothing(x):
    pass

# 1.trackbar name  2.window name 3.minimum value 4.maximum value 5.callback function
cv2.namedWindow('trackbars')
cv2.createTrackbar('sigmaX','trackbars',0,50,nothing)
cv2.createTrackbar('sigmaY','trackbars',0,50,nothing)
cv2.createTrackbar('ksize','trackbars',0,30,nothing)

while True:

    # Read the image
    image = cv2.imread('../Images and videos/flower_images/F10.jpg')

    # Get the trackbar position
    sigmaX = cv2.getTrackbarPos('sigmaX','trackbars')
    sigmaY = cv2.getTrackbarPos('sigmaY','trackbars')
    ksize = cv2.getTrackbarPos('ksize','trackbars')

    odd_ksize = [i for i in range(60) if i%2!=0]
    filtered = cv2.GaussianBlur(image,ksize = (odd_ksize[ksize],odd_ksize[ksize]), sigmaX=sigmaX,sigmaY=sigmaY,borderType=cv2.BORDER_CONSTANT)

    cv2.imshow('frame',image)
    cv2.imshow('Gaussian',filtered)

    # Press the 'ESC' to exit
    k = cv2.waitKey(1) & 0xFF
    if k ==27:
        break

cv2.destroyAllWindows()

max_stars_repo_path: src/07 - Blurring And Smoothing/03-Gaussian_filter.py

max_stars_repo_name: Robofied/Awesome-Computer-Vision-Guide

max_issues_repo_name: Robofied/Computer-Vision-Tutorials

max_forks_repo_name: Robofied/Computer-Vision-Tutorials


text:
"""
========
Nan Test
========

Example: simple line plots with NaNs inserted.
"""
import numpy as np
import matplotlib.pyplot as plt

t = np.arange(0.0, 1.0 + 0.01, 0.01)
s = np.cos(2 * 2*np.pi * t)
t[41:60] = np.nan

plt.subplot(2, 1, 1)
plt.plot(t, s, '-', lw=2)

plt.xlabel('time (s)')
plt.ylabel('voltage (mV)')
plt.title('A sine wave with a gap of NaNs between 0.4 and 0.6')
plt.grid(True)

plt.subplot(2, 1, 2)
t[0] = np.nan
t[-1] = np.nan
plt.plot(t, s, '-', lw=2)
plt.title('Also with NaN in first and last point')

plt.xlabel('time (s)')
plt.ylabel('more nans')
plt.grid(True)

plt.tight_layout()
plt.show()


text:
import os
import shutil
import sys

import numpy as np
from scipy import sparse

import matplotlib.pyplot as plt

import seaborn as sn
sn.set()

import pandas as pd

import tensorflow as tf
from tensorflow.contrib.layers import apply_regularization, l2_regularizer

import bottleneck as bn
from model import MultiDAE, MultiVAE
from metrics import *
from utils import *

def test_model_k_rounds(p_dims, noise, test_data, num_rounds, num_users, num_items, topkItens, dcg_gt, dataset):
    N_test = test_data.shape[0]
    idxlist_test = range(N_test)

    batch_size_test = 2000
    tf.reset_default_graph()
    vae = MultiVAE(p_dims, lam=0.0)
    saver, logits_var, _, _, _ = vae.build_graph(noise)
    arch_str = "I-%s-I" % ('-'.join([str(d) for d in vae.dims[1:-1]]))
    # the total number of gradient updates for annealing
    total_anneal_steps = 200000
    # largest annealing parameter
    anneal_cap = 0.2
    chkpt_dir = './chkpt/'+dataset+'/VAE_anneal{}K_cap{:1.1E}/{}'.format(
        total_anneal_steps/1000, anneal_cap, arch_str)
    print("chkpt directory: %s" % chkpt_dir)

    preds_k_rounds = []

    for k in range(num_rounds):
        preds = []
        with tf.Session(config=tf.ConfigProto(device_count={'GPU': 0})) as sess:
        #with tf.Session() as sess:
            saver.restore(sess, '{}/model'.format(chkpt_dir))
            for bnum, st_idx in enumerate(range(0, N_test, batch_size_test)):
                end_idx = min(st_idx + batch_size_test, N_test)
                X = test_data[idxlist_test[st_idx:end_idx]]
            #    X = test_data

                if sparse.isspmatrix(X):
                    X = X.toarray()
                X = X.astype('float32')

                pred_val = sess.run(logits_var, feed_dict={vae.input_ph: X})
                # exclude examples from training and validation (if any)
                pred_val[X.nonzero()] = -np.inf
                preds.extend(pred_val)
        preds_k_rounds.append(preds)
    preds_k_rounds = np.array(preds_k_rounds)
    print(preds_k_rounds.shape)

    preds_k_rounds_filter = np.zeros((num_rounds,num_users,num_items))
    for user in range(num_users):
        preds_k_rounds_filter[:,user,:] = preds_k_rounds[:,user,sorted(topkItens[user,:])]

    #preds_k_rounds = preds_k_rounds[:,:,sorted(topkItens)]
    preds_k_rounds_filter = np.array(preds_k_rounds_filter)
    print(preds_k_rounds_filter.shape)

    ufairs, ndcgs = Fairness_at_k_rounds(preds_k_rounds_filter, dcg_gt)
    return ufairs,ndcgs

def main(dataset, test_file):
    #dataset = "netflix"
    #dataset = "ml-20m"

    DATA_DIR = '../data/'+dataset+'/'
    pro_dir = os.path.join(DATA_DIR, 'pro_sg')

    unique_sid = list()
    with open(os.path.join(pro_dir, 'unique_sid.txt'), 'r') as f:
        for line in f:
            unique_sid.append(line.strip())

    n_items = len(unique_sid)
    print("NUMERO ITENS: "+str(n_items))

    p_dims = [200, 600, n_items]
    # the total number of gradient updates for annealing
    total_anneal_steps = 200000
    # largest annealing parameter
    anneal_cap = 0.2

    # Test data
    #test_data_tr, test_data_te = load_tr_te_data(
    #    os.path.join(pro_dir, 'test_tr.csv'),
    #    os.path.join(pro_dir, 'test_te.csv'),
    #    n_items)

    test_data = load_test_data(
        os.path.join(pro_dir, test_file),
        n_items)

    N_test = test_data.shape[0]
    idxlist_test = range(N_test)

    batch_size_test = 2000

    tf.reset_default_graph()
    vae = MultiVAE(p_dims, lam=0.0)
    saver, logits_var, _, _, _ = vae.build_graph(0)

    arch_str = "I-%s-I" % ('-'.join([str(d) for d in vae.dims[1:-1]]))

    chkpt_dir = './chkpt/'+dataset+'/VAE_anneal{}K_cap{:1.1E}/{}'.format(
        total_anneal_steps/1000, anneal_cap, arch_str)
    print("chkpt directory: %s" % chkpt_dir)

    n100_list, r20_list, r50_list = [], [], []
    preds = []

    with tf.Session(config=tf.ConfigProto(device_count={'GPU': 0})) as sess:
    #with tf.Session() as sess:
        saver.restore(sess, '{}/model'.format(chkpt_dir))
        for bnum, st_idx in enumerate(range(0, N_test, batch_size_test)):
            end_idx = min(st_idx + batch_size_test, N_test)
            X = test_data[idxlist_test[st_idx:end_idx]]
            #X = test_data

            if sparse.isspmatrix(X):
                X = X.toarray()
            X = X.astype('float32')

            pred_val = sess.run(logits_var, feed_dict={vae.input_ph: X})
            # exclude examples from training and validation (if any)
            pred_val[X.nonzero()] = -np.inf
            preds.extend(pred_val)

    num_users = 1000
    num_items = 100
    num_rounds = 100

    print("USUARIOS: "+str(test_data.shape[0]))
    preds = np.array(preds)
    np.random.seed(98765)
    rnd_users = np.random.choice(range(test_data.shape[0]), num_users, replace=False)
    preds = preds[rnd_users,:]
    #preds = preds[1000:1100,:]
    #plot_sorted_preds(preds)
    #print(sorted(preds))
    topkItens = np.argsort(preds)[:,-num_items:]
    print(topkItens.shape)

    predsn = []
    for user in range(topkItens.shape[0]):
        predsn.append(preds[user,sorted(topkItens[user,:])])
    predsn = np.array(predsn)

    print(predsn.shape)
    dcg_gt = dcg_k_users(predsn)
    #subSet = test_data_tr.T[topkItens,0]
    test_data = test_data[rnd_users,:]

    print("NO NOISE")
    ufairs, ndcgs = test_model_k_rounds(p_dims, 0, test_data, num_rounds, num_users, num_items, topkItens, dcg_gt, dataset)
    print(np.mean(ufairs),np.std(ufairs),np.mean(ndcgs),np.std(ndcgs))
    print("NORMAL NOISE STD 0.5")
    ufairs_n05, ndcgs_n05 = test_model_k_rounds(p_dims, 1, test_data, num_rounds, num_users, num_items, topkItens, dcg_gt, dataset)
    print(np.mean(ufairs_n05),np.std(ufairs_n05),np.mean(ndcgs_n05),np.std(ndcgs_n05))
    print("NORMAL NOISE STD 1.0")
    ufairs_n10, ndcgs_n10 = test_model_k_rounds(p_dims, 2, test_data, num_rounds, num_users, num_items, topkItens, dcg_gt, dataset)
    print(np.mean(ufairs_n10),np.std(ufairs_n10),np.mean(ndcgs_n10),np.std(ndcgs_n10))
    print("NORMAL NOISE STD 2.0")
    ufairs_n20, ndcgs_n20 = test_model_k_rounds(p_dims, 3, test_data, num_rounds, num_users, num_items, topkItens, dcg_gt, dataset)
    print(np.mean(ufairs_n20),np.std(ufairs_n20),np.mean(ndcgs_n20),np.std(ndcgs_n20))
    print("UNIFORM NOISE")
    ufairs_unif, ndcgs_unif = test_model_k_rounds(p_dims, 4, test_data, num_rounds, num_users, num_items, topkItens, dcg_gt, dataset)
    print(np.mean(ufairs_unif),np.std(ufairs_unif),np.mean(ndcgs_unif),np.std(ndcgs_unif))

    plot_comparison([ufairs,ufairs_n05, ufairs_n10,ufairs_n20, ufairs_unif],[1-ndcgs,1-ndcgs_n05, 1-ndcgs_n10,1-ndcgs_n20,1-ndcgs_unif],['original','N(std=0.5)','N(std=1.0)','N(std=2.0)','uniform'], dataset, test_file)

if __name__ == "__main__":
    main(sys.argv[1], sys.argv[2])

max_stars_repo_name: rcaborges/variational-fairness

max_issues_repo_name: rcaborges/variational-fairness

max_forks_repo_name: rcaborges/variational-fairness


text:
from __future__ import absolute_import
from __future__ import print_function
from __future__ import division

import numpy as np

import nngen.operator as operator

from . import util

def Flatten(visitor, node, no_transpose=False):

    srcs = []

    for src in node.input:
        src_obj = visitor.visit(src)
        srcs.append(src_obj)

    input = srcs[0]

    for attribute in node.attribute:
        if attribute.name == 'axis':
            axis = attribute.i

    if (not no_transpose and
        input.get_layout() is not None and input.get_onnx_layout() is not None and
            input.get_layout() != input.get_onnx_layout()):

        perm = [input.get_layout().index(l) for l in input.get_onnx_layout()]
        onnx_perm = [i for i, l in enumerate(input.get_onnx_layout())]

        input = operator.transpose(input, perm)
        input.transpose_onnx_perm = onnx_perm

        input.layout = input.get_onnx_layout()

    name = util.get_name(node)

    if axis == 0:
        shape = [1, -1]
    else:
        shape = [int(np.prod(input.shape[:axis])), -1]

    c = operator.reshape(input, shape, name=name)

    return c

max_stars_repo_name: hirayamy/nngen

max_issues_repo_name: hirayamy/nngen

max_forks_repo_name: hirayamy/nngen


text:
# Copyright 2020 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.
# ==============================================================================
# This code is a small data loader function for our cifar experiments.

import os
import numpy as np
import torch
import torchvision
import enum

class DatasetTypes(enum.Enum):
    stl10 = 'STL10'

def generate_data(args, model_cfg):
    print("Loading dataset {} from {}".format(args.dataset, args.data_path))
    dataset = getattr(torchvision.datasets, args.dataset)
    path = os.path.join(args.data_path, args.dataset.lower())

    if args.dataset == DatasetTypes.stl10:
        train_set = dataset(
            root=path, split="train", download=False, transform=model_cfg.transform_train
        )
        num_classes = 10
        # this is a manual mapping of STL10 classes to CIFAR10 classes
        # CIFAR10 classes: {automobile, bird, cat, deer, dog, frog, horse, ship}
        # STL10 classes: {airplane, bird, car, cat, deer, dog, horse, monkey, ship, truck}
        cls_mapping = np.array([0, 2, 1, 3, 4, 5, 7, 6, 8, 9])
        train_set.labels = cls_mapping[train_set.labels]

        test_set = dataset(
            root=path, split="test", download=False, transform=model_cfg.transform_test
        )
        test_set.labels = cls_mapping[test_set.labels]
    else:
        train_set = dataset(
            root=path, train=True, download=False, transform=model_cfg.transform_train
        )
        # zero indexing so the max target is one less than the number of classes
        num_classes = max(train_set.targets) + 1

        test_set = dataset(
            root=path, train=False, download=False, transform=model_cfg.transform_test
        )

    loaders = {
        "train": torch.utils.data.DataLoader(
            train_set,
            batch_size=args.batch_size,
            shuffle=True,
            num_workers=args.num_workers,
            pin_memory=True,
        ),
        "test": torch.utils.data.DataLoader(
            test_set,
            batch_size=args.batch_size,
            shuffle=False,
            num_workers=args.num_workers,
            pin_memory=True,
        ),
    }
    num_data = len(loaders["train"].dataset)
    print("Number of data points: ", num_data)
    return loaders, num_classes, num_data