text:
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import numpy as np
import array
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
import glob
import struct
import time
import sys
from mpl_toolkits.basemap import Basemap, shiftgrid, addcyclic
from scipy import interpolate
import getopt
import string
import datetime
import scipy.interpolate as interp
from scipy.spatial import cKDTree
import scipy.optimize as optm
import scipy.stats as scstats
import os.path
import math
from scipy.io import netcdf
import calendar
import datetime
sys.path.append('/home/bzhao/python_utils')
import read_utils
import plot_utils
import math_utils
import os.path

#### usage:
## run for 1987 jan whole month:

##     python interpo_aice_hice.py 1987 1 1 30

#####
## this script is modified from /discover/nobackup/zli7/geos5/Sub_seasonal/Obs_to_Tri_interp/20160301_test/interpo_aice_hice.py
## major changes include: write out hice and seaice in different files
#####
#####
## current version of this scripts are only good until year 2014
## for year 2016 has different input files for hice (heff)
## to run for 2016 onwards, gmask and gm need to be used (uncomment line 360 to 366)
#####
def write_state_aice(fname, AICE):

    nx=AICE.shape[1]
    ny=AICE.shape[0]
    nt=1

    ncfile = Dataset(fname,'w')
    ncfile.createDimension('xaxis_1',nx)
    ncfile.createDimension('yaxis_1',ny)
    ncfile.createDimension('Time',nt)

    aice = ncfile.createVariable('AICE',np.dtype('float32').char,('Time','yaxis_1','xaxis_1'))

    aice[:] = AICE
    ncfile.close()
def write_state_hice(fname, HICE):

    nx=HICE.shape[1]
    ny=HICE.shape[0]
    nt=1

    ncfile = Dataset(fname,'w')
    ncfile.createDimension('xaxis_1',nx)
    ncfile.createDimension('yaxis_1',ny)
    ncfile.createDimension('Time',nt)

    hice = ncfile.createVariable('HICE',np.dtype('float32').char,('Time','yaxis_1','xaxis_1'))

    hice[:] = HICE
    ncfile.close()

def nearest_interp(lon, lat, z, LON, LAT, undef=np.nan):

    lon[lon>80.0]=lon[lon>80.0]-360.0
    points=np.array( (lon.flatten(), lat.flatten()) ).swapaxes(0, 1)

    zout=interpolate.NearestNDInterpolator(points, z.flatten())(LON, LAT)

    return zout

def nearest_interp_new(lon, lat, z, LON, LAT):
    lon[lon>80.0]=lon[lon>80.0]-360.0
    xs, ys, zs = lon_lat_to_cartesian(lon.flatten(), lat.flatten())
    xt, yt, zt = lon_lat_to_cartesian(LON.flatten(), LAT.flatten())
    tree = cKDTree(zip(xs, ys, zs))
    #find indices of the nearest neighbors in the flattened array
    #d, inds = tree.query(zip(xt, yt, zt), k = 1)
    #get interpolated 2d field
    #zout = z.flatten()[inds].reshape(LON.shape)

    d, inds = tree.query(zip(xt, yt, zt), k = 10)
    w = 1.0 / d**2
    zout = np.sum(w * z.flatten()[inds], axis=1) / np.sum(w, axis=1)
    zout.shape = LON.shape
    return zout

def lon_lat_to_cartesian(lon, lat, R = 1):
    """
    calculates lon, lat coordinates of a point on a sphere with
    radius R
    """
    lon_r = np.radians(lon)
    lat_r = np.radians(lat)

    x =  R * np.cos(lat_r) * np.cos(lon_r)
    y = R * np.cos(lat_r) * np.sin(lon_r)
    z = R * np.sin(lat_r)
    return x,y,z

def get_grid(): #reads lat lon for tripolar ocean grid
    ##ncfile=Dataset('/gpfsm/dnb42/projects/p17/gvernier/SAND_BOXES/PLOT_ODAS/DATA/grid_spec_720x410x40.nc', "r")
    ncfile=Dataset('/discover/nobackup/yvikhlia/coupled/Forcings/Ganymed/a90x540_o720x410/INPUT/grid_spec.nc',"r")
    #ncfile=Dataset('/gpfsm/dnb42/projects/p17/bzhao/cp100/scratch/INPUT/grid_spec.nc',"r")
    LON     = ncfile.variables['x_T'][:]
    LAT     = ncfile.variables['y_T'][:]
    numlev     = ncfile.variables['num_levels'][:]
    ncfile.close()

    return LON, LAT, numlev

# Read in NSIDC

###
##choose year, month, day for aice
yearstart=int(sys.argv[1])
yearend=int(sys.argv[1])
month=int(sys.argv[2])
daystart=int(sys.argv[3])
day_interval=int(sys.argv[4])

year=yearstart
year1=year
while (year <= yearend):
       pre_year=year-1
       next_year=year+1
       day=daystart

       #if (day == 1):
       #	day1=31
       #	month1=month-1
       #else:
       #day1=day-1
       day1=day
       month1=month
       #if (day == 1 & month == 1):
       #	day1=31
       #	month1=12
       #	year1=year-1

       startdays=datetime.datetime(year,month,day,0,0,0).timetuple().tm_yday
       enddays=startdays+day_interval
       leapday1=366
       leapday2=365
       preyear=[-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1]

       if (day <= 15):
	       kk1=month-2
	       kk2=month-1
       else:
	       kk1=month-1
	       kk2=month
	       print startdays

##linear interpo monthly to daily
       XP=[]
       mm=0

       if (calendar.isleap(year)):
		calenderdays=np.arange(leapday1)
		hice_day=np.zeros((leapday1,276,360))
		nextyear=[367,368,369,370,371,372,373,374,375,376,377,378,379,380,381]
		preXP=[-16]
		nextXP=[381]
		caldayend=381
       else:
		calenderdays=np.arange(leapday2)
		hice_day=np.zeros((leapday2,276,360))
		nextyear=[366,367,368,369,370,371,372,373,374,375,376,377,378,379,380]
		preXP=[-16]
		nextXP=[380]
		caldayend=380
       calenderdays_extend=np.concatenate((preyear,calenderdays,nextyear))
       while (mm <= 11):
               nn=mm+1
               dum=datetime.datetime(year,nn,15,0,0,0).timetuple().tm_yday
               XP.append(dum)
               mm=mm+1
       XP_extend=np.concatenate((preXP,XP,nextXP))
       	# read in hice
       gmask=np.genfromtxt('/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/io.dat_360_276.output',dtype='int32',delimiter=2)
       ncfile=netcdf.netcdf_file('/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/giomass_grid.nc','r')
       dum=ncfile.variables['LON'][:]
       lon2=np.zeros(np.shape(dum))
       lon2[:]=dum
       lon2 = np.transpose(lon2)
       dum=ncfile.variables['LAT'][:]
       lat2=np.zeros(np.shape(dum))
       lat2[:]=dum
       lat2 = np.transpose(lat2)
       ncfile.close()
       print 'lon2 min: ', lon2.min()
       print 'lon2 max: ', lon2.max()
       print 'lon2 shape ', lon2.shape
       print 'lat2 shape ', lat2.shape
       if (year == 2016):
      		 filein2='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d' %year
      		 #filein2pre='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d.nc' %pre_year
      		 filein2pre='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d' %year
       		 filein2next='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d' %year
       		 ss=np.fromfile(filein2,dtype='float32')
       		 hice_p_thisyear=np.reshape(ss,(12, lon2.shape[0], lon2.shape[1]))
       		 ss_next=np.fromfile(filein2next,dtype='float32')
       		 hice_p_next=np.reshape(ss_next,(12, lon2.shape[0], lon2.shape[1]))
                 ss_pre=np.fromfile(filein2pre,dtype='float32')
                 hice_p_pre=np.reshape(ss_pre,(12, lon2.shape[0], lon2.shape[1]))

                 #ncfilepre=netcdf.netcdf_file(filein2pre,'r')
                 #dum=ncfilepre.variables['lon_scaler'][:]
                 #lon2pre=np.zeros(np.shape(dum))
                 #lon2pre[:]=dum

                 #dum=ncfilepre.variables['lat_scaler'][:]
                 #lat2pre=np.zeros(np.shape(dum))
                 #lat2pre[:]=dum

                 #dum=ncfilepre.variables['heff'][:]
                 #hice_p_pre=np.zeros(np.shape(dum))
                 #ice_p_pre[:]=dum
                 #ncfilepre.close()
       if (year == 2015):
		 filein2='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d.nc' %year

                 filein2pre='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d.nc' %pre_year
       		 filein2next='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d' %next_year
		 ncfilethis=netcdf.netcdf_file(filein2,'r')
		 dum=ncfilethis.variables['heff'][:]
	         hice_p_thisyear=np.zeros(np.shape(dum))
       		 hice_p_thisyear[:]=dum
		 ncfilepre=netcdf.netcdf_file(filein2pre,'r')
		 dum=ncfilepre.variables['heff'][:]
	         hice_p_pre=np.zeros(np.shape(dum))
       		 hice_p_pre[:]=dum
		 ss_next=np.fromfile(filein2next,dtype='float32')
       		 hice_p_next=np.reshape(ss_next,(12, lon2.shape[0], lon2.shape[1]))
		 ncfilethis.close()
		 ncfilepre.close()
       if (year < 2015):
                 filein2='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d.nc' %year

                 filein2pre='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d.nc' %pre_year
                 filein2next='/gpfsm/dnb04/projects/p71/aogcm/g5odas/obs/raw/GIOMAS/RAW/heff/heff.H%d.nc' %next_year
                 ncfilethis=netcdf.netcdf_file(filein2,'r')
                 dum=ncfilethis.variables['heff'][:]
                 hice_p_thisyear=np.zeros(np.shape(dum))
                 hice_p_thisyear[:]=dum
                 ncfilepre=netcdf.netcdf_file(filein2pre,'r')
                 dum=ncfilepre.variables['heff'][:]
                 hice_p_pre=np.zeros(np.shape(dum))
                 hice_p_pre[:]=dum
                 ncfilenext=netcdf.netcdf_file(filein2next,'r')
                 dum=ncfilenext.variables['heff'][:]
                 hice_p_next=np.zeros(np.shape(dum))
                 hice_p_next[:]=dum
                 ncfilethis.close()
                 ncfilepre.close()
                 ncfilenext.close()

       hice_p_3year=np.concatenate((hice_p_pre,hice_p_thisyear,hice_p_next))
       hice_p_3year=np.ma.masked_where(hice_p_3year>9999, hice_p_3year)
       hice_p_extend=hice_p_3year[11:25,:,:]
       print 'hice_p_3year shape ', hice_p_3year.shape

       #for n in range(1,13):
#		print n
#      		hice_p[n-1]=np.ma.masked_where(gmask==0, hice_p[n-1])
 #      lon2=np.ma.masked_where(gmask==0, lon2)
 #      lat2=np.ma.masked_where(gmask==0, lat2)
       #print 'calenderdays: ', calenderdays[:]
       #print 'XP: ', XP[:]

       ii=0
       while (ii <= 275):
               jj=0
               while (jj <= 359):
                       dum=np.interp(calenderdays_extend[15:caldayend],XP_extend[:],hice_p_extend[:,ii,jj])
                       hice_day[:,ii,jj]=dum
                       jj=jj+1
               ii=ii+1

       ######
       startdaysforhice=startdays-1
       enddaysforhice=enddays-1
       days=startdaysforhice
       print 'days: ', days
       while (days <= enddaysforhice):
         yyyymmdd=year*10000+month*100+day
         print yyyymmdd

         #### check if aice input data exist
         aice_file = '/discover/nobackup/zli7/geos5/Sub_seasonal/Obs_to_Tri_interp/INPUT/%d/gmao-obs-%d.nc' %(year,yyyymmdd)
	 fileout_hice='OUTPUT/%d/hice/seaice_hice_%d.nc' % (year,yyyymmdd)
         if os.path.isfile(aice_file):
	 	yyyymmdd1=yyyymmdd
	 else:
		yyyymmdd=year1*10000+month1*100+day1
          	print 'aice missing'
		print 'aice_file=', aice_file
	        ###
	 filein1='/discover/nobackup/zli7/geos5/Sub_seasonal/Obs_to_Tri_interp/INPUT/%d/gmao-obs-%d.nc' %(year,yyyymmdd)
	 fileout_aice='OUTPUT/%d/aice/seaice_aice_%d.nc' % (year,yyyymmdd)
#	 ncfile=Dataset(filein1,'r')
#	 dum=ncfile.variables['lon'][:]
#	 lon1=np.zeros(np.shape(dum))
#	 lon1[:]=dum
#	 lon1[lon1>80.0]=lon1[lon1>80.0]-360.0
#	 #lon1[lon1<=0.0]=lon1[lon1<=0.0]+360.0
#	 print 'lon1 ',lon1.shape
#
#	 dum=ncfile.variables['lat'][:]
#	 lat1=np.zeros(np.shape(dum))
#	 lat1[:]=dum
#
#	 dum=ncfile.variables['value'][:]
#	 aice_p=np.zeros(np.shape(dum))
#	 aice_p[:]=dum
#	 aice_p=np.squeeze(aice_p)
#	 ncfile.close()
#	 print aice_p.shape
#
#		 ######get tripolar LON LAT
#	 LON, LAT, numlevels = get_grid()
#
#	 ####find mininum lat in obs to maskout artifactual seaice after interpo
#	 print 'AA ', lon1.min(), lon1.max()
#	 print 'AA ', lat1.min(), lat1.max()
#	 ## divide by lon zone and mask zone by zone
#	 lon_zone=np.linspace(-280,80,37)
#	 min_lat_ind=np.where(aice_p > 0)
#	 #print min_lat_ind
#	 min_lat_N=np.linspace(0,36,37)
#	 min_lat_S=np.linspace(0,36,37)
#	 for ii in min_lat_ind:
#	 	yy=lat1[ii]
#	 	xx=lon1[ii]
#	 	yy_N=yy[yy >= 0] #yy_N: non zero points' lat value in NH
#	 	yy_S=yy[yy < 0]
#	 	yy_S_abs=abs(yy_S)
#	 for mm in range(1,37):
#		Izone_N=np.where((xx>lon_zone[mm-1])& (xx<=lon_zone[mm]) & (yy >= 0))
#	 	min_lat_N[mm-1]=np.min(yy_N[Izone_N]) #min_lat_N:min lat value (mask threadshould) in each lon zone in NH
#	 	#min_lat_N[mm-1]=np.min(yy[Izone_N]) #min_lat_N:min lat value (mask threadshould) in each lon zone in NH
#	 	#min_lat_S[mm-1]=np.min(yy_S_abs[Izone])
#	 	min_lat_S[mm-1]=np.min(yy_S_abs)   #SH is masked out using one value(not by zones)
#	 print 'maskout min lat in NH='+str(min_lat_N[0:36])
#	 #print 'maskout min lat in SH='+str(-min_lat_S[0])
#	 #min_lat_N[:]=44.75
#	 #min_lat_S=yy_S_abs.min()
#
#	 ###Interpo and maskout by zones
#	 aice_out = nearest_interp(lon1, lat1, aice_p, LON, LAT, undef=0.)
#	 #aice_out = nearest_interp_new(lon1, lat1, aice_p, LON, LAT)
#	 aice_out=np.squeeze(aice_out)
#	 for mm in range(1,37):
#		I=np.where((LON>lon_zone[mm-1])& (LON<=lon_zone[mm]) & (LAT>-min_lat_S[mm-1]) & (LAT<min_lat_N[mm-1]))
#	        #print mm, lon_zone[mm-1], lon_zone[mm], -min_lat_S[mm-1], min_lat_N[mm-1]
#	 	#print LON.min(), LON.max()
#	 	#print LAT.min(), LAT.max()
#	 	#print I
#	 	aice_out[I]=0.
#	 	#print aice_out.max()
#	 	#print aice_out.min()
# 	 write_state_aice(fileout_aice, aice_out) #aice_out)

	 ######get tripolar LON LAT
 	 LON, LAT, numlevels = get_grid()
         hice_temp = hice_day[days,:,:].flatten()
         hice_in = hice_temp[hice_temp<1000]
         lon_in = lon2.flatten()[hice_temp<1000]
         lat_in = lat2.flatten()[hice_temp<1000]
         print 'lon_in  shape ', lon_in.shape
         print 'lat_in  shape ', lat_in.shape
         print 'hice_in  shape ', hice_in.shape
         #gm=gmask.flatten()[hice_temp<1000]
         #hice_in = hice_temp[gm>0]
         #lon_in = lon_in[gm>0]
         #lat_in = lat_in[gm>0]
         #print 'lon_in  shape 1 ', lon_in.shape
         #print 'lat_in  shape 1', lat_in.shape
         #print 'hice_in  shape 1', hice_in.shape
         hice_out = nearest_interp(lon_in, lat_in, hice_in, LON, LAT, undef=np.nan)
         #hice_out = nearest_interp_new(lon_in, lat_in, hice_in, LON, LAT)

         write_state_hice(fileout_hice, hice_out) #hice_out)
         days = days + 1
         day = day + 1

         ###write out nc file
       year = year +1

reason: import numpy,import scipy,from scipy


text:
"""
optimized_expression.py

Finds the optimal expression pattern that maximizes coverage of synaptic partners
while minimizing coverage of neighbors.

@author Christopher Brittin
@date 2019 April 14

"""
import sys
sys.path.append('./volumetric_analysis')
sys.path.append('.')
import argparse
import numpy as np
from random import sample
from tqdm import tqdm
from random import random
import matplotlib.pyplot as plt

from mat_loader import MatLoader
from cam.expression import Matrix
import aux

class Optimize(Matrix):
    def __init__(self,fgenes,fexp):
        Matrix.__init__(self,fgenes,fexp)

    def reduced_expression(self,cells):
        idx = [self.cells[c] for c in cells]
        E = self.E[idx,:]
        jdx = np.where(E.any(axis=0))[0]
        E = E[:,jdx]
        genes = [self.gene_idx[j] for j in jdx]
        return genes,E.T

def build_synapse_matices(cells,synapses,synpartners=set([])):
    n = len(cells)
    m = len(synapses)
    cells_idx = {cells[i]:i for i in range(n)}

    P = np.zeros((n,m))
    A = np.zeros((n,m))
    jdx = 0
    for cont in synapses:
        partners = set(synapses[cont]['partners'])
        neighbors = set(synapses[cont]['neighbors'])
        nonsyn = neighbors - partners - synpartners
        for c in partners: P[cells_idx[c],jdx] = 1
        for c in nonsyn: A[cells_idx[c],jdx] = 1
        jdx += 1

    return P,A

def anneal(sol):
    old_cost = cost(sol)
    T = 1.0
    T_min = 0.00001
    alpha = 0.9
    while T > T_min:
        i = 1
        while i <= 100:
            new_sol = neighbor(sol)
            new_cost = cost(new_sol)
            ap = acceptance_probability(old_cost, new_cost, T)
            if ap > random():
                sol = new_sol
                old_cost = new_cost
            i += 1
        T = T*alpha
    return sol, cost

def syncost(E,S,A,gamma = [0.5,0.4,0.1]):
    syn = np.sum(np.dot(E,S))
    adj = np.sum(np.dot(E,A))
    sparse = np.sum(E != 0)
    J = gamma[0] * syn - gamma[1] * adj - gamma[2] * sparse
    return J

def adjcost(E,S,A,gamma = [0.5,0.4,0.1]):
    syn = np.sum(np.dot(E,S))
    adj = np.sum(np.dot(E,A))
    sparse = np.sum(E != 0)
    J = -gamma[0] * syn + gamma[1] * adj - gamma[2] * sparse
    return J

def acceptance_probability(new_cost,old_cost,T):
    p =  np.exp((new_cost - old_cost) / T)
    return min(p,1)

def perturb(coord,E):
    rdx = np.random.randint(0,len(coord))
    (a,b) = coord[rdx]
    E[a,b] = int(not E[a,b])
    return E

def run_optimization(E0,S,A,funcCost,iters=100):
    (n,m) = E0.shape
    E = np.zeros(E0.shape)
    coord = np.nonzero(E0)
    coord = list(zip(coord[0],coord[1]))
    k = int(len(coord) / 2)
    for (a,b) in sample(coord,k): E[a,b] = 1
    if not iters: iters = 10*len(coord)

    old_cost = funcCost(E,S,A)
    T = 1.0
    T_min = 0.00001
    alpha = 0.9
    cost_rec = [old_cost]
    while T > T_min:
        i = 1
        while i <= iters:
            #Perturb E
            rdx = np.random.randint(0,len(coord))
            (a,b) = coord[rdx]
            E[a,b] = int(not E[a,b])
            #rdx = np.random.randint(0,n)
            #old_row = E[rdx,:]
            #nonzero = np.nonzero(old_row)
            #E[rdx,nonzero] = 0

            #Compute new cost
            new_cost = funcCost(E,S,A)
            #Decide to accept perturbation
            ap = acceptance_probability(new_cost,old_cost,T)
            if ap > random():
                old_cost = new_cost
                cost_rec.append(new_cost)
            else:
                E[a,b] = int(not(E[a,b]))
                #E[rdx,nonzero] = 1
            i += 1
        T *= alpha
    return E,cost_rec

cam = 'mat/cam_isoforms.txt'

FOUT = 'cam_analysis/results/cam_optimize_%d_v2.csv'

if __name__=='__main__':

    parser = argparse.ArgumentParser(description=__doc__,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument('matrix',
                        action = 'store',
                        help = 'Path to matrix file')

    parser.add_argument('deg',
                        action = 'store',
                        type= int,
                        help = 'Conserved degree')

    params = parser.parse_args()

    _iter = 100
    M = MatLoader()
    M.load_left()
    C = M.load_consensus_graphs(params.deg)
    S = M.load_consensus_chemical_synapse(params.deg)

    e = Optimize(cam,params.matrix)
    e.load_genes()
    e.load_cells(sorted(C.A.nodes()))
    e.assign_expression()
    e.binarize()
    gene_count = {g:[0]*_iter for g in e.genes}

    for j in range(_iter):
        print('Iter: %d/%d'%(j+1,_iter))
        for cell in tqdm(M.left,desc='Optimize'):
            if not C.A.has_node(cell): continue
            if not C.C.has_node(cell): continue
            neigh = sorted(C.A.neighbors(cell))
            genes,E = e.reduced_expression(neigh)
            P,A = build_synapse_matices(neigh,S[cell],set(C.C.neighbors(cell)))
            #P,A = build_synapse_matices(neigh,S[cell])
            Eopt,cost_rec = run_optimization(E,P,A,syncost)
            idx = set(np.nonzero(Eopt)[0])
            _genes = [genes[i] for i in idx]
            for g in _genes: gene_count[g][j] += 1

    aux.write.from_dict(FOUT%params.deg,gene_count)
    #sorted_genes = sorted(gene_count.items(), key=lambda kv: kv[1])
    #for g in sorted_genes: print(g)

    """
    print(np.sum(E),np.sum(Eopt))
    print(np.sum(np.dot(E,P)),np.sum(np.dot(E,A)))
    print(np.sum(np.dot(Eopt,P)),np.sum(np.dot(Eopt,A)))
    idx = set(np.nonzero(Eopt)[0])
    _genes = [genes[i] for i in idx]
    print(len(genes),len(_genes))
    print(sorted(_genes))
    plt.figure()
    plt.plot(cost_rec,'b-')
    plt.show()
    """

ext: py

reason: import numpy


text:
__all__ = ["processGrid"]

import copy
import cupy as cp
import numpy as np
import logging
from LevelSetPy.Utilities.matlab_utils import *
from LevelSetPy.BoundaryCondition import addGhostPeriodic

logger = logging.getLogger(__name__)

def processGrid(gridIn, data=None, sparse_flag=False):
    """
     processGrid: Construct a grid data structure, and check for consistency.

       gridOut = processGrid(gridIn, data)

     Processes all the various types of grid argument allowed.

     Inp.t Parameters:

       gridIn: A scalar, a vector, or a structure.

         Scalar: It is assumed to be the dimension.  See below for default
         settings for other grid fields.

         Vector: It contains the number of grid nodes in each dimension.  See
         below for default settings for other fields.

         Structure: It must contain some subset of the following fields
         (where each vector has length equal to the number of dimensions):

    	      gridIn.dim: Positive integer scalar, dimension of the grid.

    	      gridIn.min: Double vector specifying the lower left corner of the grid.

    	      gridIn.max: Double vector specifying the upper right corner of the grid.

    	      gridIn.N: Positive integer vector specifying the number of grid
    	      nodes in each dimension.

    	      gridIn.dx: Positive double vector specifying the grid spacing in
    	      each dimension.

    	      gridIn.vs: Cell vector, each element is a vector of node locations
    	      for that dimension.

    	      gridIn.xs: Cell vector, each element is an array of node locations
    	      (result of calling ndgrid on vs).

    	      gridIn.bdry: Cell vector of function handles pointing to boundary
    		    condition generating functions for each dimension.

           gridIn.bdryData: Cell vector of data structures for the boundary
    		    condition generating functions.

           gridIn.axis: Vector specifying computational domain bounds in a
           format suitable to pass to the axis() command (only defined for 2D
           and 3D grids, otherwise grid.axis == []).

           gridIn.shape: Vector specifying grid node count in a format suitable
           to pass to the reshape() command (usually grid.N', except for 1D
           grids).

         If any of the following fields are scalars, they are replicated
         gridIn.dim times: min, max, N, dx, bdry, bdryData.

         In general, it is not necessary to supply the fields: vs, xs, axis, shape.

         If one of N or dx is supplied, the other is inferred.
         If both are supplied, consistency is checked.

         Dimensional consistency is checked on all fields.

         Default settings (only used if value is not given or inferred)
    	      min   = zeros(dim, 1)
    	      max   = ones(dim, 1)
    	      N     = 101
    	      bdry  = periodic

       data: Double array.  Optional.  If present, the data array is checked
       for consistency with the grid.

       sparse_flag: Whether to make a sparse or dense grid. Default = False

     Output Parameters:

       gridOut: the full structure described for gridIn above.

     Copyright 2004 Ian M. Mitchell (mitchell@cs.ubc.ca).
     This software is used, copied and distributed under the licensing
       agreement contained in the file LICENSE in the top directory of
       the distribution.

     Ian Mitchell, 1/22/03
      new version  5/13/03 added fields dim, dx, vs, xs, bdry.
      new version  1/13/04 added field bdryData.
      new version  2/09/04 added field shape.
      new version  8/23/12 fixed some floating point problems with N and dx.

     Lekan Molux, Python, 7/26/2021
     ----------------------------------------------------------------------------
    """
    defaultMin = 0
    defaultMax = 1
    defaultN = 101
    defaultBdry = addGhostPeriodic
    defaultBdryData = None

    # This is just to avoid attempts to allocate 100 dimensional arrays.
    maxDimension = 5

    if not isinstance(gridIn, Bundle):
        if len(gridIn) == 1:
            gridOut = Bundle(dict(dim = gridIn))
        elif(ndims(gridIn) == 2):
            # Should be a vector of node counts.
            if(gridIn.shape[1] != 1):
                error('gridIn vector must be a column vector')
            else:
               gridOut = Bundle(dict(dim =  len(gridIn), N = gridIn))
        else:
            error('Unknown format for gridIn parameter')
    else:
        gridOut = copy.copy(gridIn)

    # Now we should have a partially complete structure in gridOut.

    if(isfield(gridOut, 'dim')):
        if(gridOut.dim > maxDimension):
            warn('Grid dimension > {}, may be dangerously large'.format(maxDimension))
        if(gridOut.dim < 0):
            error('Grid dimension must be positive')
    else:
        error('Grid structure must contain dimension')

    # Process grid boundaries.
    if(isfield(gridOut, 'min')):
        if(not isColumnLength(gridOut.min, gridOut.dim)):
            if(isscalar(gridOut.min)):
                gridOut.min = gridOut.min * np.ones((gridOut.dim, 1))
            else:
                error('min field is not column vector of length dim or a scalar')
    else:
        gridOut.min = defaultMin * np.ones((gridOut.dim, 1))

    if(isfield(gridOut, 'max')):
        if(not isColumnLength(gridOut.max, gridOut.dim)):
            if(isscalar(gridOut.max)):
                gridOut.max = gridOut.max * np.ones((gridOut.dim, 1))
            else:
                error('max field is not column vector of length dim or a scalar')
    else:
        gridOut.max = defaultMin * np.ones((gridOut.dim, 1))

    if(np.any(gridOut.max <= gridOut.min)):
        error('max bound must be strictly greater than min bound in all dimensions')

    # Check N field if necessary.  If N is missing but dx is present, we will
    # determine N later.
    if(isfield(gridOut, 'N')):
        if(np.any(gridOut.N <= 0)):
            error('number of grid cells must be strictly positive')
        if(not isColumnLength(gridOut.N, gridOut.dim)):
            if(isscalar(gridOut.N)):
                gridOut.N *= np.ones((gridOut.dim, 1)).astype(np.int64)
            else:
                error('N field is not column vector of length dim or a scalar')

    # Check dx field if necessary.  If dx is missing but N is present, infer
    # dx.  If both are present, we will check for consistency later.  If
    # neither are present, use the defaults.
    if isfield(gridOut, 'dx'):
        if(np.any(gridOut.dx <= 0)):
            error('grid cell size dx must be strictly positive')
        if(not isColumnLength(gridOut.dx, gridOut.dim)):
            if(isscalar(gridOut.dx)):
                gridOut.dx *= np.ones((gridOut.dim, 1), np.float64)
        else:
            error('dx field is not column vector of length dim or a scalar')
    elif isfield(gridOut, 'N'):
        gridOut.dx = np.divide(gridOut.max - gridOut.min,  gridOut.N-1)
    else:
        logger.warn('Neither fields dx nor dN is present, so use default N and infer dx')
        gridOut.N = defaultN * ones(gridOut.dim, 1).astype(np.int64)
        gridOut.dx = np.divide(gridOut.max - gridOut.min, gridOut.N-1)

    if isfield(gridOut, 'vs'):
        if(iscell(gridOut.vs)):
            if(not isColumnLength(gridOut.vs, gridOut.dim)):
                error(f'vs field is not column cell vector of length dim: {gridOut.dim}')
            else:
                for i in range(gridOut.dim):
                    if(not isColumnLength(gridOut.vs[i], gridOut.N[i])):
                        error('vs cell entry is not correctly sized vector')
        else:
            error('vs field is not a cell vector')
    else:
        gridOut.vs = cell(gridOut.dim)
        for i in range(gridOut.dim):
            gridOut.vs[i] = np.expand_dims(np.linspace(gridOut.min[i,0].item(), gridOut.max[i,0].item(), num=gridOut.N[i,0].item()), 1)
    # Now we can check for consistency between dx and N, based on the size of
    # the vectors in vs.  Note that if N is present, it will be a vector.  If
    # N is not yet a field, set it to be consistent with the size of vs.

    if isfield(gridOut, 'N'):
        for i in range(gridOut.dim):
            if(gridOut.N[i,0] != len(gridOut.vs[i])):
                error(f'Inconsistent grid size in dimension {i}')
    else:
        gridOut.N = np.zeros((gridOut.dim, 1), dtype=np.float64)

    for i in range(gridOut.dim):
        gridOut.N[i,0] = len(gridOut.vs[i])

    if(isfield(gridOut, 'xs')):
        if(iscell(gridOut.xs)):
            if(not isColumnLength(gridOut.xs, gridOut.dim)):
                error(f'xs field is not column cell vector of length/dim: {gridOut.dim}')
            else:
                if(gridOut.dim > 1):
                    for i in range(gridOut.dim):
                        if(np.any(gridOut.xs[i]) != gridOut.N.T):
                            error('xs cell entry is not correctly sized array')
                else:
                    if(len(gridOut.xs[0]) != gridOut.N):
                        error('xs cell entry is not correctly sized array')
        else:
            error('xs field is not a cell vector')
    else:
        gridOut.xs = np.meshgrid(*gridOut.vs, indexing='ij', sparse=sparse_flag)

    if isfield(gridOut, 'bdry'):
        if(iscell(gridOut.bdry) or isinstance(gridOut.bdry, np.ndarray)):
            if(not isColumnLength(gridOut.bdry, gridOut.dim)):
                error(f'bdry field is not column cell vector of length dim: {gridOut.dim}')
            else:
                pass
        else:
            if(isscalar(gridOut.bdry)):
                bdry = gridOut.bdry
                gridOut.bdry = [bdry for _ in range(gridOut.dim)]
            else:
                error('bdry field is not a cell vector or a scalar')
    else:
        gridOut.bdry = np.zeros((gridOut.dim, 1)).fill(defaultBdry)

    if(isfield(gridOut,'bdryData')):
        if(iscell(gridOut.bdryData)):
            if(not isColumnLength(gridOut.bdryData, gridOut.dim)):
                error(f'bdryData field is not column cell vector of length dim: {gridOut.dim}')
            else:
                logger.warn('Maybe not worth checking that entries are structures')
        else:
            if(isscalar(gridOut.bdryData)):
                bdryData = gridOut.bdryData
                gridOut.bdryData = [bdryData for i in range(gridOut.dim)]
            else:
                error('bdryData field is not a cell vector or a scalar')
    else:
        gridOut.bdryData = [defaultBdryData for i in range(gridOut.dim)]

    if((gridOut.dim == 2) or (gridOut.dim == 3)):
        if(isfield(gridOut, 'axis')):
            for i in range(gridOut.dim):
                if(gridOut.axis.take(2 * i) != gridOut.min[i]):
                    error('axis and min fields do not agree')
                if(gridOut.axis.take(2 * i) is not gridOut.max[i]):
                    error('axis and max fields do not agree')
        else:
            gridOut.axis = np.zeros((1, 2 * gridOut.dim), dtype=np.float64)
            for i in range(gridOut.dim):
                gridOut.axis[0, 2*i : 2*i+2] = [ gridOut.min[i,0].item(), gridOut.max[i,0].item() ]
    else:
        gridOut.axis = []

    Nshape = gridOut.N.squeeze()
    Nshape = tuple([x.item() for x in Nshape])
    if(isfield(gridOut, 'shape')):
        if(gridOut.dim == 1):
            if(np.any(gridOut.shape != (Nshape + (1,)) )):
                error('shape and N fields do not agree')
        else:
            if(np.any(gridOut.shape != gridOut.N.T)):
                error('shape and N fields do not agree')
    else:
        if(gridOut.dim == 1):
            gridOut.shape = (Nshape + (1,))
        else:
            gridOut.shape = Nshape

    # check data parameter for consistency
    if data:
        if(ndims(data) != len(gridOut.shape)):
            error('data parameter does not agree in dimension with grid')
        if(np.any(size(data) != gridOut.shape)):
            error('data parameter does not agree in array size with grid')

    return  gridOut

max_stars_repo_name: robotsorcerer/LevelSetPy

max_issues_repo_name: robotsorcerer/LevelSetPy

max_forks_repo_name: robotsorcerer/LevelSetPy

reason: import numpy,import cupy


text:
import trimesh
import numpy as np

def scale_to_unit_sphere(mesh):
    if isinstance(mesh, trimesh.Scene):
        mesh = mesh.dump().sum()

    vertices = mesh.vertices - mesh.bounding_box.centroid
    distances = np.linalg.norm(vertices, axis=1)
    vertices /= np.max(distances)

    return trimesh.Trimesh(vertices=vertices, faces=mesh.faces)

def scale_to_unit_cube(mesh):
    if isinstance(mesh, trimesh.Scene):
        mesh = mesh.dump().sum()

    vertices = mesh.vertices - mesh.bounding_box.centroid
    vertices *= 2 / np.max(mesh.bounding_box.extents)

    return trimesh.Trimesh(vertices=vertices, faces=mesh.faces)

voxel_points = dict()

def get_raster_points(voxel_resolution):
    if voxel_resolution in voxel_points:
        return voxel_points[voxel_resolution]

    points = np.meshgrid(
        np.linspace(-1, 1, voxel_resolution),
        np.linspace(-1, 1, voxel_resolution),
        np.linspace(-1, 1, voxel_resolution)
    )
    points = np.stack(points)
    points = np.swapaxes(points, 1, 2)
    points = points.reshape(3, -1).transpose().astype(np.float32)

    voxel_points[voxel_resolution] = points
    return points

def check_voxels(voxels):
    block = voxels[:-1, :-1, :-1]
    d1 = (block - voxels[1:, :-1, :-1]).reshape(-1)
    d2 = (block - voxels[:-1, 1:, :-1]).reshape(-1)
    d3 = (block - voxels[:-1, :-1, 1:]).reshape(-1)

    max_distance = max(np.max(d1), np.max(d2), np.max(d3))
    return max_distance < 2.0 / voxels.shape[0] * 3**0.5 * 1.1

def sample_uniform_points_in_unit_sphere(amount, mesh_bounds):
    ellipsiod_dim = mesh_bounds*1.5 # bigger than the shape
    x,y,z = ellipsiod_dim

    sample = np.empty((0,3))
    while len(sample) < amount:
        unit_sphere_points = np.random.uniform(-ellipsiod_dim,
                                            ellipsiod_dim,
                                            size=(amount * 2 + 20, 3))
                                # use ellipsoid bigger than the dimensions
        # check if sample points are inside the ellipsoid
        ellipsiod_check = unit_sphere_points[:,0]**2/x**2 + \
                        unit_sphere_points[:,1]**2/y**2 + \
                        unit_sphere_points[:,2]**2/z**2
        unit_sphere_points = unit_sphere_points[ellipsiod_check < 1]

        ## unit sphere, original
        # unit_sphere_points = np.random.uniform(-1, 1, size=(amount * 2 + 20, 3))
        # unit_sphere_points = unit_sphere_points[np.linalg.norm(unit_sphere_points, axis=1) < 1]

        sample = np.vstack((sample, unit_sphere_points))
    return sample[:amount]


text:
import random
import numpy as np

def binarycode(s):
    size = s
    c = []
    for x in range(size):
        if random.random()<0.5:
            c.append("0")
        else:
            c.append("1")
    return c

def randomChange(c, s, x, p, rand):
    size = s
    code = c
    spot = random.randint(0, size-1)
    #if random.random()<0.5:
    #    code[spot]="0"
    #else:
    #    code[spot]="1"
    #print(code)
    if (code[spot] == "1"):
        code[spot] = "0"
    else :
        code[spot] = "1"
    #print("Individuo ", x, " fue mutado en: ", spot, " y quedó como ", code[spot], " random ",rand, " probabilidad ", p )
    #print(code)
    return code

def toDecimals(c, ll, ul, size):
    question = c
    answer = 0

    answer = ll+((question/(pow(2,size)-1))*(ul-ll))

    return answer

def evaluation (q):

    return pow(q, 2)

def torneo (ps, p, dp):
    torneos = pobSize // 2
    poblacion = p
    pobDecimal = dp
    padres = []
    rango = (len(poblacion))
    contador = 1

    for x in range (torneos):
        tmpa = random.randrange(rango)
        tmpb = random.randrange(rango)
        winner  = 0

        while tmpa == tmpb:
            tmpb = random.randrange(rango)

        if (abs(pobDecimal[tmpa]) < abs(pobDecimal[tmpb])):
            padres.append(poblacion[tmpa])
            winner = tmpa
        elif (abs(pobDecimal[tmpa]) > abs(pobDecimal[tmpb])):
            padres.append(poblacion[tmpb])
            winner = tmpb
        else :
            padres.append(poblacion[tmpa])
            winner = tmpa

        #print("Torneo ",contador,": Individuo ", tmpa, " vs ", tmpb, " Ganador : ", winner)
        contador+=1
        winner2 = winner

        while (winner2 == winner ):
            tmpa = random.randrange(rango)
            tmpb = random.randrange(rango)

            while tmpa == tmpb:
                tmpb = random.randrange(rango)

            if (abs(pobDecimal[tmpa]) < abs(pobDecimal[tmpb])):
                winner2 = tmpa
            elif (abs(pobDecimal[tmpa]) > abs(pobDecimal[tmpb])):
                winner2 = tmpb
            else :
                winner2 = tmpa
        padres.append(poblacion[winner2])
        #print("Torneo ",contador,": Individuo ", tmpa, " vs ", tmpb, " Ganador : ", winner2)
        contador+=1

    return padres

def cruza (p, ps, s, pc):
    padres = p
    size = ps // 2
    rango = (len(padres)-1)+1
    longitud = s
    hijos = []
    arrTemp1 = []
    arrTemp2 = []
    contador = 0
    for x in range (size):
        r = random.random()
        if (r<pc):
            arrTemp1=[]
            arrTemp2=[]
            puntoCruza = random.randrange(rango) + 1
            switch = True
            for y in range (longitud):
                if puntoCruza == y:
                    switch = False
                if (switch):
                    arrTemp1.append(padres[contador][y])
                    arrTemp2.append(padres[contador+1][y])
                else:
                    arrTemp1.append(padres[contador+1][y])
                    arrTemp2.append(padres[contador][y])
            #print("Pareja - ", (x+1)," punto de cruza - ", puntoCruza, ", numero aleatorio - ", r, ", probabilidad - ", pc)
            hijos.append(arrTemp1)
            hijos.append(arrTemp2)
        else:
            hijos.append(padres[contador])
            hijos.append(padres[contador+1])
            #print("Pareja - ", (x+1)," sin cruza, numero aleatorio - ", r, ", probabilidad - ", pc)
        contador+=2
    return hijos

def mutacion(h, s, p):
    hijos = h
    longitudCadena = s
    probMutacion = p
    size = len(hijos)
    for x in range (size):
        rand = random.random()
        if (rand < probMutacion):
            hijos[x] = randomChange(hijos[x], longitudCadena, x, p, rand)

    return hijos

size = 12
pobSize = 10
upplimit = 2
lowlimit = -2
probabilidadCruza = 0.6
probabilidadMutacion = 0.1
minimo = 100
condicionParo = 20

mejor = []
pob = []
parsedpob = []
decimalpob = []
padres = []
hijos = []

test = np.array([])

code = binarycode(size)
#print (code)
for x in range (pobSize):
    temp = binarycode(size)
    pob.append(temp)
    parsedpob.append(int("".join(temp), 2))
    decimalpob.append(toDecimals(parsedpob[x], lowlimit, upplimit, size))
    #a =  np.array([])
    #test = np.append(a, x)
#print (pobSize)
#print (pob)
#print (parsedpob)
#print (decimalpob)

#for x in range (pobSize):
 #   print("Individuo ", x, " - ","".join(pob[x]),  " ", parsedpob[x], " ", round(decimalpob[x], 3), " " , round(evaluation(decimalpob[x]), 3))

for x in range (pobSize):
    if (evaluation(decimalpob[x]) < minimo):
        mejor = decimalpob[x]
        minimo = evaluation(decimalpob[x])

print("Generación inicial : Mejor valor ",mejor, " - F(x): ",minimo)

for v in range (condicionParo):

    padres = torneo(size, pob, decimalpob)

    #for x in range (pobSize):
    #   print("Padre ",x," - ","".join(padres[x]))

    #print("--Cruza--")
    hijos = cruza(padres, pobSize, size, probabilidadCruza)

    #for x in range (pobSize):
    #   print("Hijo ",x," - ","".join(hijos[x]))

    #print ("--Mutación--")

    hijos = mutacion(hijos, size, probabilidadMutacion)

    pob = hijos

    for x in range (pobSize):
        parsedpob[x]=(int("".join(pob[x]), 2))
        decimalpob[x]=(toDecimals(parsedpob[x], lowlimit, upplimit, size))
        if (evaluation(decimalpob[x]) < minimo):
            mejor = decimalpob[x]
            minimo = evaluation(decimalpob[x])

    print("Generación ", (v+1),": Mejor valor ",mejor, " - F(x): ",minimo)

print ("El mejor valor fue ", mejor, " - F(x): ",minimo, "  -  " , round(evaluation(decimalpob[x]), 5))
#for x in range (pobSize):

#for x in range (pobSize):
 #   print("Hijo ",x," - ","".join(hijos[x]))

#print (test)
#code = randomChange(code, size)
#print (code)
#parsedValue = int("".join(pob[5]), 2)
#print (parsedValue)
#print(round(toDecimals(273,lowlimit, upplimit, size),3))

max_stars_repo_name: LeoCruzG/Genetic-AI-exercises

max_issues_repo_name: LeoCruzG/Genetic-AI-exercises

max_forks_repo_name: LeoCruzG/Genetic-AI-exercises


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

#----------------------------------------------------------------------------------------------------------
"""
Author : Mohamed

Title : Implementation des equations de saint-venant

Description : Simulation des vagues dans un domaine [0,1] avec les equations de saint-venant

"""
#-----------------------------------------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt

#Params
g = 9.8

def F(U):
	"""
	U[0] : hauteur
	U[1] : debit
	"""

	return np.array([U[1],(U[1]**2)/U[0] + g*(U[0]**2)/2])

def f(U_g,U_d):
	"""
	U_g/d : vecteur a gauche/droite
	"""

	#vitesse
	v_g = U_g[1]/U_g[0]
	v_d = U_d[1]/U_d[0]

	c = np.amax([np.abs(v_g)+np.sqrt(g*U_g[0]), np.abs(v_d)+np.sqrt(g*U_d[0])])

	return (F(U_g)+F(U_d))/2 - (c/2)*(U_d - U_g)

#----------------------------------------------------------------------------------------------------------
# Mayage
#-------------------------

N = 100
dx = 1.0/N

x = np.linspace(0, 1, num=N+1)  # le nombre d'element est N+1 car on a N intervalle

#[End]----------------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------------------------
# Conditions initiales et de bords
#-------------------------

v_0 = np.zeros(N+1) #vitesse nulle a t = 0

#Condition discontinue sur les hauteurs
h_0 = np.ones(N+1)
h_0[N//2:] = np.ones(len(h_0[N//2:])) * 10

#Condition continue
# h_0 = np.linspace(0,1,N+1)
# h_0 = np.exp(h_0)

U_0 = np.array([ h_0, h_0*v_0 ])

#[End]----------------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------------------------
# Definition des tableaux
#-------------------------

U_tab = np.array([U_0]) #tableau de vecteur qui va contenir les vecteurs U (qui contient les donnees de tout le maillage) a tout instant.

time = np.zeros(1) #Tableau des instants

#[End]----------------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------------------------
# Implementation
#-------------------------

def lambda1(v,h):
	return v - np.sqrt(g*h)

def lambda2(v,h):
	return v + np.sqrt(g*h)

#debut de l'evolution
n = 0 #--> n designe le temps

while n < 100:

	h_n = U_tab[n][0] #tableau hauteurs h_n
	q_n = U_tab[n][1] #tableau debits q_n

	v_n = q_n/h_n #tableau vitesses v_n

	m = 2 * np.max([np.abs(lambda1(v_n,h_n)), np.abs(lambda2(v_n,h_n))])
	dt = dx/m

	time = np.append(time, time[n]+dt) #Ajout dans le tableau de temps

	# U_nplus1 = U_tab[n]  #--> fait gaffe la, car U_nplus1 prend la meme adress que U_tab[n] !!!!!
	U_nplus1 = np.array([np.ones(N+1), np.ones(N+1)]) #Tableau qui contient les valeurs de vecteurs U a l'instant n+1 sur tout le maillage

	for i in range(N+1): # i -> designe l'endroit sur le maillage

		# if i < 3 or i > N - 3:  #Près de bords -> vitesse nulle
		# 	u_in = np.array([h_n[i], 0]) #u_in vecteur a l'endroit i mais a vitesse nulle

		# 	u_ig = np.array([h_n[(i-1)%(N+1)], 0]) #definition de u_gauche a l'instant n vitesse nulles aux bords
		# 	u_id = np.array([h_n[(i+1)%(N+1)], 0])

		# 	u_i_nplus1 = u_in - (dt/dx)*(f(u_in, u_id)-f(u_ig, u_in))

		# 	#Stocker les params dans le nouveau tableau de vecteur a l'instant n+1 a l'endroit i
		# 	U_nplus1[0][i] = u_i_nplus1[0]
		# 	U_nplus1[1][i] = 0 #Vitesse nulle

		# else: #Loin de bords

	#Periodisation :
		u_in = np.array([h_n[i], q_n[i]]) #u_in vecteur a l'endroit i qui contient les donnees h_n[i],q_n[i] a l'instant n et a l'endroit i

		u_ig = np.array([h_n[(i-1)%(N+1)], q_n[(i-1)%(N+1)]]) #definition de u_gauche a l'instant n avec condition periodique
		u_id = np.array([h_n[(i+1)%(N+1)], q_n[(i+1)%(N+1)]])

		u_i_nplus1 = u_in - (dt/dx)*(f(u_in, u_id)-f(u_ig, u_in)) #def de vecteur U a l'instant n+1 a l'endroit i

		#Stocker les params dans le nouveau tableau de vecteur a l'instant n+1 a l'endroit i
		U_nplus1[0][i] = u_i_nplus1[0]
		U_nplus1[1][i] = u_i_nplus1[1]

	U_tab = np.concatenate((U_tab, [U_nplus1])) #Ajout des nouveaux donnees de l'instant n+1

	n = n+1 #incrementation

#[End]----------------------------------------------------------------------------------------------------

#----------------------------------------------------------------------------------------------------------
# Animation de plot
#-------------------------

fig = plt.figure()

plt.ion() # mode interactif
plt.figure(1) # nouvelle figure

plt.plot(x, U_tab[0][0])
plt.show()

plt.pause(0.01) # nécessaire (bug matplotlib)
plt.hold(False) # Ne pas superposer les figures pour animer le mouvement

for j in range(1,n):

    plt.plot(x, U_tab[j][0])  # on redessine
    plt.ylim([0, 11]) #limits sur l'axis y
    plt.pause(0.0001) # on attend un peu

plt.ioff()        # fin du mode interactif
plt.show()        #nécessaire

#[End]----------------------------------------------------------------------------------------------------

max_stars_repo_path: saint_venant_periodic boundaries.py

max_issues_repo_path: saint_venant_periodic boundaries.py

max_forks_repo_path: saint_venant_periodic boundaries.py


text:
from pycuda.compiler import SourceModule
import pycuda.gpuarray as garray
import numpy as np

size_of_curandStateXORWOW = 12

def curand_setup(num_threads, seed):
    """
    Setup curand seed
    """
    func = get_curand_int_func()
    grid = ( (int(num_threads)-1)/128 + 1,1)
    block = (128,1,1)

    # curandStateXORWOW is a struct of size 12 4bytes
    state = garray.empty(num_threads*size_of_curandStateXORWOW, np.int32)

    func.prepared_call( grid, block, state.gpudata, num_threads, seed)
    return state

def get_curand_int_func():
    code = """
#include "curand_kernel.h"
extern "C" {
__global__ void
rand_setup(curandStateXORWOW_t* state, int size, unsigned long long seed)
{
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    int total_threads = blockDim.x * gridDim.x;

    for(int i = tid; i < size; i+=total_threads)
    {
        curand_init(seed, i, 0, &state[i]);
    }
}
}
    """
    mod = SourceModule(code, no_extern_c = True)
    func = mod.get_function("rand_setup")
    func.prepare('PiL')#[np.intp, np.int32, np.uint64])
    return func

max_stars_repo_name: MangoHaha/neurodriver

max_issues_repo_name: MangoHaha/neurodriver

max_forks_repo_name: MangoHaha/neurodriver

reason: import numpy,import pycuda,from pycuda


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

from __future__ import absolute_import, division, print_function, unicode_literals

import tensorflow as tf

from tensorflow.keras.layers import Dense, Flatten, Conv2D, MaxPooling2D, UpSampling2D,Conv2DTranspose, Concatenate
from tensorflow.keras import Model, regularizers
import re
import csv
import numpy as np
import matplotlib.pyplot as plt
from glob import glob
from time import time
import matplotlib.animation as animation
import os
import yaml
from encoder_decoder import autoencoder
import shutil
import argparse

stream = open("config.yaml", 'r')
config = yaml.load(stream)
print ("Using the following settings:")
for key, value in config.items():
    print (key + " : " + str(value))
chip_size_x = config["chip_size_x"]
chip_size_y = config["chip_size_y"]
tile_size =  float(config["tile_size"])
time_steps = config["num_time_steps"]
percent_valid_split = config["percent_valid_split"]
learning_rate_val = config["learning_rate"]
decay_steps_num =  config["decay_steps"]
decay_rate_val = config["decay_rate"]
epochs_num = config["epochs"]

# x_data is input image frame, t_data is time, y_data is golden output image frame
def train_model(x_data, t_data, y_data):
    # Create an instance of the model
    model = autoencoder()
    initial_learning_rate = learning_rate_val
    lr_schedule = tf.keras.optimizers.schedules.ExponentialDecay(
        initial_learning_rate,
        decay_steps=decay_steps_num,
        decay_rate=decay_rate_val,
        staircase=True)
    model.compile(optimizer=tf.keras.optimizers.Adam(learning_rate=lr_schedule),
                  loss='mse',
                  metrics=['mse', 'mae', 'mape'])

    st = time()
    history = model.fit(x=[x_data,t_data],y=y_data, epochs=epochs_num,
                validation_split=percent_valid_split,
                shuffle=True)
    et = time()
    tt = et-st
    print("Elapsed time: %03d:%02d:%05.2f"%(int(tt/3600),int(tt/60)%60,tt%60))
    return model

def predict_temperature(model, normalization_data, test_data_dir):
    data_runs = [f for f in glob(test_data_dir+"**/Run_*_contour*", recursive=False)]
    if (len(data_runs) == 0): exit("Could not find transient test runs with Run_\d_contour filenames")

    power_map = np.zeros((len(data_runs),time_steps,chip_size_x,chip_size_y))
    temp_map = np.zeros((len(data_runs),time_steps,chip_size_x,chip_size_y))

    num_test_data = 0;
    for i, run in enumerate(data_runs):
        print("Reading data: ",run)
        run_num = re.search(r'Run_\d+_contour', run)
        if (run_num == None):
            exit("Directory paths must match structure and naming convention of data set to UMN in Transient_runs directory")
        run_num = run_num.group()
        run_num = re.findall(r'\d+', run_num)
        run_num = int(run_num[0])
        for frame in range(time_steps):
            fname = run+"/ml_raw_data_%d_%d.csv"%(run_num,frame+1)
            power_map[i,frame,...],temp_map[i,frame,...] = read_image(fname)
        num_test_data+=1

    max_x = normalization_data[0]
    min_x = normalization_data[1]
    max_t = normalization_data[2]
    min_t = normalization_data[3]
    max_y = normalization_data[4]
    min_y = normalization_data[5]
    Writer = animation.writers['html']
    writer = Writer(fps=5, metadata=dict(artist='Me'), bitrate=1800)
    im = []
    figs= []
    anis = []
    num_im = num_test_data
    predicted_video = np.zeros((num_im,time_steps,chip_size_x,chip_size_y))
    err_video = np.zeros((num_im,time_steps,chip_size_x,chip_size_y))
    #get_ipython().run_line_magic('matplotlib', 'notebook')
    for im_num in range(num_im):
        def updatefig1(frame,im_num):
            im[(im_num)*3].set_array(temp_map[im_num,frame,...])
            im[(im_num)*3+1].set_array(predicted_video[im_num,frame,...])

        def updatefig2(frame,im_num):
            im[(im_num)*3+2].set_array(err_video[im_num,frame,...])

        for frame in range(time_steps):
            data_pt = np.array([[frame]])/max_t
            in_data = power_map[im_num:im_num+1,0,...,np.newaxis]
            in_data = (in_data-min_x)/(max_x-min_x)
            predicted = model.predict((in_data,data_pt))
            predicted_video[im_num,frame,...] = (np.squeeze(predicted)*(max_y-min_y))+min_y
            err_video[im_num,frame,...] = predicted_video[im_num,frame,...] - temp_map[im_num,frame,...]
        print("Writing result of test: ", data_runs[im_num])
        test_data_num = re.search(r'Run_\d+_contour', data_runs[im_num])
        if (test_data_num == None):
            exit("Directory paths must match structure and naming convention of data set to UMN in Transient_runs directory")
        test_data_num = test_data_num.group();
        test_data_num = re.findall(r'\d+', test_data_num)
        test_data_num = int(test_data_num[0])
        fig, axes = plt.subplots(1,2)
        max_val = np.max((np.max(temp_map[im_num,...]), np.max(predicted_video[im_num,...])))
        max_err = np.max(err_video[im_num,...])
        min_err = np.min(err_video[im_num,...])
        im.append(axes[0].imshow(temp_map[im_num,0], vmin=0, vmax=max_val))
        im.append(axes[1].imshow(predicted_video[im_num,0,...], vmin=0, vmax=max_val))
        fig.colorbar(im[3*(im_num)+1], ax=axes.ravel().tolist())

        anis.append(animation.FuncAnimation(fig, updatefig1, frames=range(time_steps), fargs=(im_num,),
                                      interval=150 ))
        anis[2*(im_num)].save("output_plots/run_contours_%d.html"%test_data_num, writer=writer)
        plt.close()
        fig2 = plt.figure()
        im.append(plt.imshow(err_video[im_num,0,...], vmin=min_err, vmax=max_err))
        plt.colorbar()
        anis.append(animation.FuncAnimation(fig2, updatefig2, frames=range(time_steps), fargs=(im_num,),
                                      interval=150 ))
        anis[2*(im_num)+1].save("output_plots/err_contours_%d.html"%test_data_num, writer=writer)
        plt.close()
        dir = os.path.join("./","output_plots","output_plots")
        if not os.path.exists(dir):
            os.mkdir(dir)
            os.rename("output_plots/run_contours_%d_frames"%test_data_num,"output_plots/output_plots/run_contours_%d_frames"%test_data_num )
            os.rename("output_plots/err_contours_%d_frames"%test_data_num,"output_plots/output_plots/err_contours_%d_frames"%test_data_num )
        else:
            if (os.path.exists("output_plots/output_plots/run_contours_%d_frames"%test_data_num)
            and os.path.exists("output_plots/output_plots/err_contours_%d_frames"%test_data_num)):
                shutil.rmtree("output_plots/output_plots/run_contours_%d_frames"%test_data_num )
                shutil.rmtree("output_plots/output_plots/err_contours_%d_frames"%test_data_num )
            os.rename("output_plots/run_contours_%d_frames"%test_data_num,"output_plots/output_plots/run_contours_%d_frames"%test_data_num )
            os.rename("output_plots/err_contours_%d_frames"%test_data_num,"output_plots/output_plots/err_contours_%d_frames"%test_data_num )

def read_data(data_dir):
    data_runs = [f for f in glob(data_dir+"**/Run_*_contour*", recursive=False)]
    if (len(data_runs) == 0): exit("Could not find transient training runs with Run_\d_contour filenames")

    power_map = np.zeros((len(data_runs),time_steps,chip_size_x,chip_size_y))
    temp_map = np.zeros((len(data_runs),time_steps,chip_size_x,chip_size_y))

    num_images = 0;
    for i, run in enumerate(data_runs):
        print("Reading data: ",run)
        run_num = re.search(r'Run_\d+_contour', run)
        if (run_num == None):
            exit("Directory paths must match structure and naming convention of data set to UMN in Transient_runs directory")
        run_num = run_num.group()
        run_num = re.findall(r'\d+', run_num)
        run_num = int(run_num[0])
        for frame in range(time_steps):
            fname = run+"/ml_raw_data_%d_%d.csv"%(run_num,frame+1)
            power_map[i,frame,...],temp_map[i,frame,...] = read_image(fname)
            num_images+=1

    count = 0
    x_data = np.zeros((num_images,chip_size_x,chip_size_y,1))
    y_data = np.zeros((num_images,chip_size_x,chip_size_y,1))
    t_data = np.zeros((num_images,1))
    for im_num in range(len(data_runs)):
        for frame in range(time_steps):
            x_data[count,:,:,0] = power_map[im_num,time_steps-1,...]
            y_data[count,:,:,0] = temp_map[im_num,frame,...]
            t_data[count,...] = frame
            count +=1
    return [x_data, y_data, t_data]

def process_data(x_data, y_data, t_data):
    indices = np.arange(x_data.shape[0])
    np.random.shuffle(indices)

    x_data = x_data[indices]
    y_data = y_data[indices]
    t_data = t_data[indices]

    min_x = np.min(x_data)
    max_x = np.max(x_data)
    x_data = (x_data-min_x)/(max_x-min_x)

    max_t =time_steps-1
    min_t = 0
    t_data = (t_data)/max_t

    min_y=np.min(y_data)
    max_y = np.max(y_data)
    y_data = (y_data-min_y)/(max_y-min_y)
    normalization_data = [max_x, min_x, max_t, min_t, max_y, min_y]
    return x_data, y_data, t_data, normalization_data

def read_image(fname):
   power_map = np.zeros((chip_size_x,chip_size_y))
   temp_map = np.zeros((chip_size_x,chip_size_y))
   with open(fname) as csvfile:
       readCSV = csv.reader(csvfile, delimiter=',')
       for row in readCSV:
           x = int(np.round(float(row[1])/tile_size))
           y = int(np.round(float(row[2])/tile_size))
           dyn_pow = float(row[3])
           leak_pow = float(row[4])
           alpha = float(row[5])
           power_map[x,y] = alpha*dyn_pow + leak_pow
           temp_map[x,y] = float(row[7])
   return power_map,temp_map

def main():

    parser = argparse.ArgumentParser(description="Training and inference flow for transient thermal analysis")
    parser.add_argument("-train_data_path",
                        help="Provide path to the directory containing Transient_runs and Steady_run folders with training data",
                        required=True)
    parser.add_argument("-test_data_path",
                       help="Provide path to the directory containing Transient_runs and Steady_run folders with test data",
                       required=True)
    parser.add_argument("-output_plot_dir",
                        help="Provide path to generate the output plots",
                       required=True)
    args = parser.parse_args()
    train_data_dir =  args.train_data_path
    test_data_dir =  args.test_data_path
    output_plot_dir = args.output_plot_dir

    data_train = read_data(train_data_dir)
    #data_test = read_data(test_data_dir)
    x_data, y_data, t_data, normalization_data = process_data(data_train[0],
            data_train[1], data_train[2])
    model = train_model(x_data, t_data, y_data)
    predict_temperature(model, normalization_data, test_data_dir)
    #plot_test_results(predicted_video, err_video, temp_map, output_plot_dir)

if __name__ == '__main__':
    main()

max_stars_repo_name: VidyaChhabria/ThermEDGe-and-IREDGe

max_issues_repo_name: VidyaChhabria/ThermEDGe-and-IREDGe

max_forks_repo_name: VidyaChhabria/ThermEDGe-and-IREDGe


text:
# %%[markdown]
# Author: Nelson Liu
#
# Email: [nliu@uncharted.software](mailto:nliu@uncharted.software)

# %%[markdown]
# Content:
# * Load data of the full EMMAA Covid-19 graph
# * Filter statements by a given document DOI/PMID/PMCID
# * Get the list of nodes and edges associated with these statements
# * Get documents referenced by these edges
# * Filter statements by these documents
# * Get the expanded list of nodes and edges
# * Generate the associated output files (`nodes.jsonl, edges.jsonl, G.pkl, ...`)

# %%
import json
import pickle
import time
import numpy as np
import scipy as sp
import networkx as nx
import matplotlib as mpl
import matplotlib.pyplot as plt

import emmaa_lib as emlib
import importlib
# importlib.reload(emlib)

# %%
np.random.seed(0)

# %%[markdown]
# # Define the document of interest
# doc = {'id_type': 'DOI', 'id': '10.1016/j.immuni.2020.04.003'.upper()}
doc = {'id_type': 'PMID', 'id': '32325025'}

# %%[markdown]
# # Load Statement Data
statements_all = {}
with open('./data/covid19-snapshot_dec8-2020/source/latest_statements_covid19.json', 'r') as x:
    statements_all = json.load(x)

x = None
del x

# %%[markdown]
# # Filter Statements by Referenced Document
s_doc = set([i for i, s in enumerate(statements_all) for ev in s['evidence'] if ('text_refs' in ev.keys()) if (doc['id_type'] in ev['text_refs'].keys()) if ev['text_refs'][doc['id_type']] == doc['id']])

# statements_all = None
# del statements_all

# %%
print(f"{len(s_doc)} statements found to reference the document with {doc['id_type']} = {doc['id']}")
# 7 statements found to reference the document with DOI = 10.1016/J.IMMUNI.2020.04.003

# %%[markdown]
# # Load the Nodes and Edges Data of the Model Graph

nodes_model = emlib.load_jsonl('./dist/v3/nodes_model.jsonl', remove_preamble = True)
edges_model = emlib.load_jsonl('./dist/v3/edges_model.jsonl', remove_preamble = True)

# %%[markdown]
# # Filter the Node and Edge List

x = [statements_all[s]['matches_hash'] for s in s_doc]
e_doc = [e for e, edge in enumerate(edges_model) if edge['statement_id'] in x]
n_doc = set([edges_model[e]['source_id'] for e in e_doc] + [edges_model[e]['target_id'] for e in e_doc])

x = None
del x

# %%
print(f"{len(n_doc)} nodes and {len(e_doc)} edges found to be associated with the document with {doc['id_type']} = {doc['id']}")
# 10 nodes and 8 edges found to be associated with the document with DOI = 10.1016/J.IMMUNI.2020.04.003

# %%[markdown]
# # Get Documents Referenced by the Edge List

# %%

e_exp = set([e for n in n_doc for e in nodes_model[n]['edge_ids_source']] + [e for n in n_doc for e in nodes_model[n]['edge_ids_target']])
n_exp = set([edges_model[e]['source_id'] for e in e_exp] + [edges_model[e]['target_id'] for e in e_exp])

x = set([edges_model[e]['statement_id'] for e in e_exp])
s_exp = [s for s, statement in enumerate(statements_all) if statement['matches_hash'] in x]
docs_exp = list(set([ev['text_refs']['PMID'] for s in s_exp for ev in statements_all[s]['evidence'] if ('text_refs' in ev.keys()) if 'PMID' in ev['text_refs'].keys()]))

x = None
del x

# %%
print(f"The previous {len(n_doc)} nodes are linked to {len(e_exp)} edges, which are linked to {len(n_exp)} nodes and reference {len(docs_exp)} documents in total.")
# The previous 10 nodes are linked to 12715 edges, which are linked to 4924 nodes and reference 9252 documents in total.

# %%
# # Plot Results

# %%

nodes_doc = [nodes_model[n] for n in n_doc]
edges_doc = [edges_model[e] for e in e_doc]

G_doc = emlib.generate_nx_object(nodes_doc, edges_doc)

__ = emlib.generate_nx_layout(G = G_doc, layout = 'kamada_kawai', plot = True,
    plot_atts = {'verticalalignment': 'top','font_size': 10, 'with_labels': True, 'arrows': True,
    'node_size': [10 * node[1] for node in G_doc.degree()], 'width': 1,
    'labels': {node: G_doc.nodes[node]['name'] for node in G_doc.nodes()},
    'cmap': 'cool'})
# fig.savefig('./figures/v3/boutique_subgraph_layout_degree.png', dpi = 150)

# %%

nodes_exp = [nodes_model[n] for n in n_exp]
edges_exp = [edges_model[e] for e in e_exp]

G_exp = emlib.generate_nx_object(nodes_exp, edges_exp)

__, __, fig, __ = emlib.generate_nx_layout(G = G_exp, layout = 'spring', layout_atts = {'k': 0.5}, plot = True,
    plot_atts = {'verticalalignment': 'top','font_size': 10, 'with_labels': False, 'arrows': False,
    'node_size': [node[1] for node in G_exp.degree()], 'width': 1,
    'labels': {node: G_exp.nodes[node]['name'] for node in G_exp.nodes()},
    'cmap': 'cool'})
# fig.savefig('./figures/v3/boutique_subgraph_layout_degree.png', dpi = 150)

# %%
# Save `nodes`
preamble = {
    'model_id': '<int> unique model ID that is present in all related distribution files',
    'id': '<int> unique node ID that is referenced by other files',
    'name': '<str> unique human-interpretable name of this node (from the `name` attribute in `latest_statements_covid19.jsonl`)',
    'db_refs': '<dict> database references of this node (from the `db_refs` attribute in `latest_statements_covid19.jsonl`)',
    'grounded': '<bool> whether this node is grounded to any database',
    'edge_ids_source': '<list of int> ID of edges that have this node as a source',
    'edge_ids_target': '<list of int> ID of edges that have this node as a target',
    'out_degree': '<int> out-degree of this node',
    'in_degree': '<int> in-degree of this node',
}
emlib.save_jsonl(nodes_doc, './dist/v3/doc/nodes_doc.jsonl', preamble = preamble)

# %%
# Save `edges`
preamble = {
    'model_id': '<int> unique model ID that is present in all related distribution files',
    'id': '<int> unique edge ID that is referenced by other files',
    'type': '<str> type of this edge (from `type` attribute in `latest_statements_covid19.jsonl`)',
    'belief': '<float> belief score of this edge (from `belief` attribute in `latest_statements_covid19.jsonl`)',
    'statement_id': '<str> unique statement id (from `matches_hash` in `latest_statements_covid19.jsonl`)',
    'source_id': '<int> ID of the source node (as defined in `nodes.jsonl`)' ,
    'target_id': '<int> ID of the target node (as defined in `nodes.jsonl`)',
    'tested': '<bool> whether this edge is tested'
}
emlib.save_jsonl(edges_doc, './dist/v3/doc/edges_doc.jsonl', preamble = preamble)

preamble = None
del preamble

# %%
# Generate a `G.pkl`

G_doc = emlib.generate_nx_object(nodes_doc, edges_doc)
with open('./dist/v3/doc/G_doc.pkl', 'wb') as x:
    pickle.dump(G_doc, x)

map_ids_edges = {edge['id']: i for i, edge in enumerate(edges_doc)}
# x = [max([edges_doc[map_ids_edges[edge_id]]['belief'] for edge_id in node[1]]) if len(node[1]) > 0 else 0.0 for node in G_doc.nodes.data('edge_ids_target')]
# __, __, fig, __ = emlib.generate_nx_layout(G = G_doc, layout = 'spring', layout_atts = {'k': 0.1}, plot = True, plot_atts = {'node_color': x, 'vmin': 0.0, 'vmax': 1.0, 'cmap': 'cool'})
__, __, fig, __ = emlib.generate_nx_layout(G = G_doc, layout = 'spring', layout_atts = {'k': 0.1}, plot = True, plot_atts = {'width': 1, 'node_size': 100, 'labels': {node: G_doc.nodes[node]['name'] for node in G_doc.nodes()}, 'cmap': 'cool'})
fig.savefig('./figures/v3/doc_subgraph_layout_degree.png', dpi = 150)

G_doc = x = fig = None
del x, fig, G_doc

# %%[markdown]
# # Impose Ontology and Generate Hyperedges

# %%
%%time

with open('./data/bio_ontology_v1.5.json', 'r') as x:
    G_onto_JSON = json.load(x)

# Remove 'xref' links
G_onto_JSON['links'] = [link for link in G_onto_JSON['links'] if link['type'] != 'xref']

# Generate a namespace list common to the model graph and the ontology
namespaces_priority = ['FPLX', 'UPPRO', 'HGNC', 'UP', 'CHEBI', 'GO', 'MESH', 'MIRBASE', 'DOID', 'HP', 'EFO']
namespaces, namespaces_count = emlib.generate_ordered_namespace_list(namespaces_priority, G_onto_JSON, nodes_doc)

# Reduce 'db_refs' of each model node to a single entry by namespace priority
nodes_doc, __ = emlib.reduce_nodes_db_refs(nodes_doc, namespaces)

# Calculate in-ontology paths
emlib.calculate_onto_root_path(nodes_doc, G_onto_JSON)

# Extract Ontological Categories
ontocats_doc = emlib.extract_ontocats(nodes_doc, G_onto_JSON)

# Generate Hyperedges
hyperedges_doc = emlib.generate_hyperedges(nodes_doc, edges_doc, ontocats_doc)

# Compute Layout using Ontological Categories and Hyperedges
__ = emlib.generate_onto_layout(nodes_doc, ontocats_doc, hyperedges_doc, plot = True)

# %%[markdown]
# # Save Outputs

# %%
# Save layout of `nodes` as `nodeLayout`
preamble = {
    'model_id': '<int> unique model ID that is present in all related distribution files',
    'id': '<int> unique node ID that is defined in `nodes.jsonl`',
    'x': '<float> position of the node in the graph layout',
    'y': '<float> position of the node in the graph layout',
    'z': '<float> position of the node in the graph layout',
}
emlib.save_jsonl(nodes_doc, './dist/v3/doc/nodeLayout_doc.jsonl', preamble = preamble)

# Save new attributes of `nodes` as `nodeAtts`
preamble = {
    'model_id': '<int> unique model ID that is present in all related distribution files',
    'id': '<int> unique node ID that is defined in `nodes.jsonl`',
    'db_ref_priority': '<str> database reference from `db_refs` of `nodes.jsonl`, that is used by the INDRA ontology v1.5',
    'grounded_onto': '<bool> whether this model node is grounded to something that exists within the ontology',
    'ontocat_level': '<int> level of the ontology node/category to which this model node was mapped (`-1` if not mappable, `0` if root)',
    'ontocat_ids': '<array of int> ordered list of ontological category IDs (see `ontocats.jsonl`) to which this node is mapped (order = root-to-leaf)',
    'cluster_ids': '<array of int> (placeholder for cluster)'
}
emlib.save_jsonl(nodes_doc, './dist/v3/doc/nodeAtts_doc.jsonl', preamble = preamble)

# Save `ontocats`
preamble = {
    'model_id': '<int> unique model ID that is present in all related distribution files',
    'id': '<int> unique ID for the ontological category that is referenced by other files',
    'ref': '<str> unique reference ID of this category (as given by the INDRA Ontology v1.5)',
    'name': '<str> name of this category (as given by the INDRA Ontology v1.5)',
    'size': '<int> number of model nodes that were mapped to this category and its children',
    'level': '<int> number of hops to reach the local root (`0` if root)',
    'parent_id': '<int> ID of the parent of this category in the ontology',
    'children_ids': '<array of int> unordered list of the child category IDs',
    'node_ids': '<array of int> unordered list of IDs from model nodes in the membership of this category',
    'node_ids_direct': '<array of int> node_ids but only model nodes which were directly mapped to this category and not any of the child categories',
    'hyperedge_ids': '<array of int> unordered list of hyperedge IDs (see `hyperedges.jsonl`) that are within this category',
}
emlib.save_jsonl(ontocats_doc, './dist/v3/doc/ontocats_doc.jsonl', preamble = preamble)

# Save layout of ontocats as `ontocatLayout`
preamble = {
    'model_id': '<int> unique model ID that is present in all related distribution files',
    'id': '<int> unique ID for the ontological category that is referenced by other files',
    'x': '<float> position of the node in the graph layout',
    'y': '<float> position of the node in the graph layout',
    'z': '<float> position of the node in the graph layout'
}
emlib.save_jsonl(ontocats_doc, './dist/v3/doc/ontocatLayout_doc.jsonl', preamble = preamble)

# Save `hyperedges`
preamble = {
    'model_id': '<int> unique model ID that is present in all related distribution files',
    'id': '<int> unique hyperedge ID that is referenced by other files',
    'level': '<float> ontological level of this hyperedge (taken from the source ontocat)',
    'size': '<int> number of model edges that is aggregated here',
    'edge_ids': '<array of int> unordered list of edge ID of the underlying model edges (see `edges.jsonl`)',
    'source_type': '<str> object type of the source (only `ontocat`)',
    'source_id': '<int> ID of the source (defined in `ontocats.jsonl`)' ,
    'target_type': '<str> object type of the target (either `ontocat` or `node`)',
    'target_id': '<int> ID of the target (defined in either `ontocats.jsonl` or `nodes.jsonl`)',
}
emlib.save_jsonl(hyperedges_doc, './dist/v3/doc/hyperedges_doc.jsonl', preamble = preamble)

# %%

max_stars_repo_name: uncharted-aske/research

max_issues_repo_name: uncharted-aske/research

max_forks_repo_name: uncharted-aske/research

reason: import numpy,import scipy,import networkx


text:
import flood_simulator_stats as fss
from itertools import cycle
import numpy as np

from constants import *

cf=fss.StatsExplorer("../sim_results/naiveflood_10000", 0)
mf=fss.StatsExplorer("../sim_results/mpr_10000", 0)

def get_markers ():
        markers = cycle(('s','+', '.', 'o', '*'))
        colors  = cycle(('blue', 'red', 'green', 'black'))
        linestyles = cycle(('solid', 'dashed', 'dotted'))
        return (markers, colors, linestyles)

def adjustFigAspect(fig,aspect=1):
    '''
    Adjust the subplot parameters so that the figure has the correct
    aspect ratio.
    '''
    xsize,ysize = fig.get_size_inches()
    minsize = min(xsize,ysize)
    xlim = .4*minsize/xsize
    ylim = .4*minsize/ysize
    if aspect < 1:
        xlim *= aspect
    else:
        ylim /= aspect
    fig.subplots_adjust(left=.5-xlim, right=.5+xlim, bottom=.5-ylim, top=.5+ylim)

cf_s, mf_s = "solid", "dashed"

pc = {
        4 : 0.15,
        5:  .18,
        6: .21,
        7: .24,
        8: .27,
        9: .295,
        10: .32,
        11: .345,
        12: .37,
        13: .395,
        14: .42,
        15: .445,
        16: .47,
        17: .49,
        18: .51,
        19: .535,
        20: .56,
        22: .59,
        24: .63,
        26: .66,
        28: .68,
        30: .71,
        32: .73,
        34: .76,
        36: .78,
        38: .81,
        40: .83,
        42: .84,
        44: .86,
        46: .87,
        48: .88,
        50: .89
        }

PAYLOAD_SIZE = 32
NODE_ID_SIZE = 2
SLOC = 8
NUMBER_NAME_COMP = 3
TTL = 4

def nb_tries_per_message(pc, mtr=5):
    return (1.-(pc**mtr)*(mtr+1.)+(pc**(mtr+1))*mtr)/(1.-pc)

def total_number_transmitted_messages(density, ttl):
    d=density
    return np.sum(cf.exps[d].get_means()[0:ttl])*nb_tries_per_message(pc[d])

# -ICN_CNT_SEC_LEN: We don't need the ICN signature since the interest is never forwarded
# 2+1+NODE_ID_SIZE: NameComponent TLVs for /ndb/node_id
def beacon_size(sloc):
    return HDR_802154_LEN + ICN_CNT_HDR_LEN - ICN_CNT_SEC_LEN + 2 + 1 + NODE_ID_SIZE + sloc

def name_size(number_of_nodes):
    return NUMBER_NAME_COMP + np.ceil(np.log2(number_of_nodes))

def geo_energy_capacity(number_of_nodes, density, change_frequency, payload_size=PAYLOAD_SIZE, sloc=SLOC, ttl=TTL):
    nsize = name_size(number_of_nodes)
    p = pc[density]
    ntr = nb_tries_per_message(p)
    res = float(AA_J)
    control = ((beacon_size(sloc) * ntr * COST_TX_UJ + CRYPTO_UJ) + density * (beacon_size(sloc) * ntr * COST_RX_UJ + CRYPTO_UJ)) / J_IN_UJ
    res = res - control
    denominator = control/change_frequency + ttl * (4*CRYPTO_UJ + ntr*(icn_geo_int_frame_size(nsize, sloc)+icn_cnt_frame_size(nsize, payload_size))*(COST_RX_UJ+COST_TX_UJ)) / J_IN_UJ
    return res/denominator

def mpr_energy_capacity(number_of_nodes, density, change_frequency, payload_size=PAYLOAD_SIZE, ttl=TTL):
    nsize = name_size(number_of_nodes)
    p = pc[density]
    ntr = nb_tries_per_message(p)
    total = min(total_number_transmitted_messages(density, ttl), number_of_nodes*ntr/2)
    res = float(AA_J)
    control = (float(total*2*CRYPTO_UJ)/ntr + total*icn_int_frame_size(nsize)*(COST_RX_UJ+COST_TX_UJ)) / J_IN_UJ
    res = res - control
    denominator = control/change_frequency + ttl * (4*CRYPTO_UJ + ntr*(icn_int_frame_size(nsize)+icn_cnt_frame_size(nsize, payload_size))*(COST_RX_UJ+COST_TX_UJ)) / J_IN_UJ
    return res/denominator

DENSITIES = np.arange(4, 15, 1)
NODES     = np.arange(10, 750)
X,Y = np.meshgrid(DENSITIES, NODES)
colorbar = (0.9, 1.8)

max_stars_repo_name: marceleng/geographic-icthings

max_issues_repo_name: marceleng/geographic-icthings

max_forks_repo_name: marceleng/geographic-icthings


text:
"""

Tests module catalog

# Author: Vladan Lucic
# $Id$
"""
from __future__ import unicode_literals
from __future__ import absolute_import

__version__ = "$Revision$"

from copy import copy, deepcopy
import unittest
import os.path

import numpy
import numpy.testing as np_test
import scipy

from pyto.analysis.catalog import Catalog
from pyto.analysis.test import common

class TestCatalog(np_test.TestCase):
    """
    """

    def setUp(self):
        """
        """

        #self.relative_dir = 'all_catalogs'

    def makeInstance(self):
        """
        Returns an istance of Catalog()
        """

        curr_dir, base = os.path.split(__file__)
        rel_dir =  os.path.join(curr_dir, common.rel_catalogs_dir)
        catal = Catalog()
        catal.read(dir=rel_dir, catalog=r'catalog_[0-9]*\.',
                   type='distributed')

        return catal

    def testRead(self):
        """
        Tests read
        """

        ###################################################
        #
        # Absolute path to catalogs directory
        #

        curr_dir, base = os.path.split(os.path.abspath(__file__))
        abs_dir =  os.path.join(curr_dir, common.rel_catalogs_dir)
        catal = Catalog()
        catal.read(dir=abs_dir, catalog=r'catalog_[0-9]*\.',
                   type='distributed')
        #print 'catal._db: ', catal._db

        # test feature
        np_test.assert_equal(catal._db['feature']['exp_1'], 'good')
        np_test.assert_equal(catal._db['feature']['exp_2'], 'bad')
        np_test.assert_equal(catal._db['feature']['exp_3'], 'good')
        np_test.assert_equal(catal._db['feature']['exp_4'], 'bad')

        # test results
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_b',
                                               'results_1.dat'))
        np_test.assert_equal(
            os.path.realpath(catal._db['results_file']['exp_1']), desired)
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_b',
                                               'results_2.dat'))
        np_test.assert_equal(
            os.path.realpath(catal._db['results_file']['exp_2']), desired)
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_a',
                                               'results_3.dat'))
        np_test.assert_equal(
            os.path.realpath(catal._db['results_file']['exp_3']), desired)
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_a',
                                      'catalogs_c', 'res', 'results_4.dat'))
        np_test.assert_equal(
            os.path.realpath(catal._db['results_file']['exp_4']), desired)

        ###################################################
        #
        # Absolute path to catalogs directory, test identifiers
        #

        catal = Catalog()
        catal.read(dir=abs_dir, catalog=r'catalog_[0-9]*\.',
                   type='distributed', identifiers=['exp_2', 'exp_3'])

        # test feature
        np_test.assert_equal(catal._db['feature'].get('exp_1', 'foo'), 'foo')
        np_test.assert_equal(catal._db['feature']['exp_2'], 'bad')
        np_test.assert_equal(catal._db['feature']['exp_3'], 'good')
        np_test.assert_equal(catal._db['feature'].get('exp_4', 'foo'), 'foo')

        ###################################################
        #
        # Absolute path to catalogs directory
        #

        curr_dir, base = os.path.split(__file__)
        rel_dir =  os.path.join(curr_dir, common.rel_catalogs_dir)
        catal = Catalog()
        catal.read(dir=rel_dir, catalog=r'catalog_[0-9]*\.',
                   type='distributed')
        #print 'catal._db: ', catal._db

        # test results
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_b',
                                               'results_1.dat'))
        np_test.assert_equal(catal._db['results_file']['exp_1'], desired)
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_b',
                                               'results_2.dat'))
        np_test.assert_equal(catal._db['results_file']['exp_2'], desired)
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_a',
                                               'results_3.dat'))
        np_test.assert_equal(catal._db['results_file']['exp_3'], desired)
        desired = os.path.realpath(os.path.join(curr_dir, 'catalogs_a',
                                      'catalogs_c', 'res', 'results_4.dat'))
        np_test.assert_equal(catal._db['results_file']['exp_4'], desired)

    def testMakeGroups(self):
        """
        Tests makeGroups
        """

        # read db
        #curr_dir, base = os.path.split(__file__)
        #rel_dir =  os.path.join(curr_dir, common.rel_catalogs_dir)
        #catal = Catalog()
        #catal.read(dir=rel_dir, catalog=r'catalog_[0-9]*\.',
        #           type='distributed')
        #print 'catal._db: ', catal._db

        # use all identifiers
        catal = self.makeInstance()
        catal.makeGroups()
        #for name in catal.getProperties():
        #    print "property: ", name, "  value: ", getattr(catal, name)
        desired = {'first' : {'exp_2':'bad', 'exp_3':'good'},
                   'second' : {'exp_1':'good', 'exp_4':'bad'}}
        np_test.assert_equal(catal.feature, desired)

        # specify all identifiers
        catal = self.makeInstance()
        catal.makeGroups(include=['exp_2', 'exp_3', 'exp_1', 'exp_4'])
        desired = {'first' : {'exp_2':'bad', 'exp_3':'good'},
                   'second' : {'exp_1':'good', 'exp_4':'bad'}}
        np_test.assert_equal(catal.feature, desired)

        # specify some identifiers
        catal = self.makeInstance()
        catal.makeGroups(include=['exp_2', 'exp_3', 'exp_4'])
        desired = {'first' : {'exp_2':'bad', 'exp_3':'good'},
                   'second' : {'exp_4':'bad'}}
        np_test.assert_equal(catal.feature, desired)

        # specify some identifiers
        catal = self.makeInstance()
        catal.makeGroups(include=['exp_2', 'exp_3'])
        desired = {'first' : {'exp_2':'bad', 'exp_3':'good'}}
        np_test.assert_equal(catal.feature, desired)

        # some properties not in all experiments
        catal = self.makeInstance()
        catal._db['partial'] = {'exp_2' : 2, 'exp_3' : 3}
        catal._db['categ'] = {'exp_2' : 2, 'exp_3' : 3}
        # print catal._db.keys()
        # failed to find name that comes first in _db.keys()
        catal.makeGroups()
        desired = {'first' : {'exp_2':2, 'exp_3':3}}
        np_test.assert_equal(catal.partial, desired)
        desired = {'first' : {'exp_2':'bad', 'exp_3':'good'},
                   'second' : {'exp_1':'good', 'exp_4':'bad'}}
        np_test.assert_equal(catal.feature, desired)

        # use all identifiers
        catal = self.makeInstance()
        catal.makeGroups()
        #for name in catal.getProperties():
        #    print "property: ", name, "  value: ", getattr(catal, name)
        desired = {'first' : {'exp_2':'bad', 'exp_3':'good'},
                   'second' : {'exp_1':'good', 'exp_4':'bad'}}
        np_test.assert_equal(catal.feature, desired)

        # put everything in one (dummy) group
        catal = self.makeInstance()
        catal.makeGroups(
            feature=None, singleGroupName='ggg', singleFeature='fff')
        desired = set(
            ['results_file', 'results_obj', 'cleft_columns_obj', 'feature',
             'pixel_size', 'category', 'fff'])
        np_test.assert_equal(set(catal.getProperties()), desired)
        desired = {'ggg' : {'exp_2':'ggg', 'exp_3':'ggg',
                    'exp_1':'ggg', 'exp_4':'ggg'}}
        np_test.assert_equal(catal.fff, desired)
        desired = {'ggg' : {
            'exp_2':'bad', 'exp_3':'good', 'exp_1':'good', 'exp_4':'bad'}}
        np_test.assert_equal(catal.feature, desired)

    def testGetProperties(self):
        """
        Tests getProperties()
        """

        # make instance
        catal = self.makeInstance()
        catal.makeGroups()

        # test
        props = catal.getProperties()
        desired = set(
            ['results_file', 'results_obj', 'cleft_columns_obj', 'feature',
             'pixel_size', 'category'])
        np_test.assert_equal(set(props), desired)

    def testAdd(self):
        """
        Tests add()
        """

        # read catalog
        curr_dir, base = os.path.split(os.path.abspath(__file__))
        abs_dir =  os.path.join(curr_dir, common.rel_catalogs_dir)
        catal = Catalog()
        catal.read(dir=abs_dir, catalog=r'catalog_[0-9]*\.',
                   type='distributed')

        # add new property
        new_values = {'exp_1' : 'new_1', 'exp_3' : 'new_3'}
        catal.add(name='new', values=new_values)
        np_test.assert_equal(catal._db['new']['exp_1'], 'new_1')
        np_test.assert_equal(catal._db['new'].get('exp_2', "none"), 'none')
        np_test.assert_equal(catal._db['new']['exp_3'], 'new_3')
        np_test.assert_equal(catal._db['new'].get('exp_4', "none"), 'none')
        #print catal._db

        # add an existing property, overwrite
        new_values = {'exp_2' : 'new_2', 'exp_3' : 'new_33'}
        catal.add(name='new', values=new_values, overwrite=True)
        np_test.assert_equal(catal._db['new']['exp_1'], 'new_1')
        np_test.assert_equal(catal._db['new']['exp_2'], 'new_2')
        np_test.assert_equal(catal._db['new']['exp_3'], 'new_33')
        np_test.assert_equal(catal._db['new'].get('exp_4', "none"), 'none')

        # add an existing property, not overwrite
        new_values = {'exp_2' : 'new_2', 'exp_5' : 'new_55'}
        np_test.assert_raises(
            ValueError, catal.add,
            **{'name':'new', 'values':new_values, 'overwrite':False})

        # test other features
        np_test.assert_equal(catal._db['feature']['exp_1'], 'good')
        np_test.assert_equal(catal._db['feature']['exp_2'], 'bad')
        np_test.assert_equal(catal._db['feature']['exp_3'], 'good')
        np_test.assert_equal(catal._db['feature']['exp_4'], 'bad')

    def testPool(self):
        """
        Test pool()
        """

        # make instance
        catal = self.makeInstance()
        catal.makeGroups(feature='category')

        # test
        catal.pool(categories=['first', 'second'], name='together')
        desired = {'exp_2' : 3.2, 'exp_3' : 3.3}
        np_test.assert_equal(catal.pixel_size['first'], desired)
        desired = {'exp_1' : 3.1, 'exp_4' : 3.4}
        np_test.assert_equal(catal.pixel_size['second'], desired)
        desired = {'exp_1' : 3.1, 'exp_2' : 3.2, 'exp_3' : 3.3, 'exp_4' : 3.4}
        np_test.assert_equal(catal.pixel_size['together'], desired)

if __name__ == '__main__':
    suite = unittest.TestLoader().loadTestsFromTestCase(TestCatalog)
    unittest.TextTestRunner(verbosity=2).run(suite)

reason: import numpy,import scipy


text:
# Copyright (c) Meta Platforms, Inc. and affiliates.
#
# This source code is licensed under the MIT license found in the
# LICENSE file in the root directory of this source tree.

import collections
import time
import logging

from typing import Callable, Optional, Sequence, Tuple, Union

import numpy as np
import torch

import rlmeta.core.remote as remote
import rlmeta.utils.data_utils as data_utils
import rlmeta_extension.nested_utils as nested_utils

from rlmeta.core.launchable import Launchable
from rlmeta.core.segment_tree import SumSegmentTree, MinSegmentTree
from rlmeta.core.server import Server
from rlmeta.core.types import Tensor, NestedTensor
from rlmeta_extension import CircularBuffer

class ReplayBuffer(remote.Remotable, Launchable):

    def __init__(self,
                 capacity: int,
                 collate_fn: Optional[Callable[[Sequence[NestedTensor]],
                                               NestedTensor]] = None):
        self._buffer = CircularBuffer(capacity)
        if collate_fn is not None:
            self._collate_fn = collate_fn
        else:
            self._collate_fn = data_utils.stack_tensors

    def __len__(self) -> int:
        return self.size

    def __getitem__(self, index: Union[int, Tensor]) -> NestedTensor:
        data = self._buffer[index]
        if not isinstance(index, int):
            data = nested_utils.collate_nested(self._collate_fn, data)
        return data

    @property
    def size(self) -> int:
        return len(self._buffer)

    @property
    def capacity(self) -> int:
        return self._buffer.capacity()

    @property
    def cursor(self) -> int:
        return self._buffer.cursor()

    def init_launching(self) -> None:
        pass

    def init_execution(self) -> None:
        pass

    @remote.remote_method(batch_size=None)
    def get_size(self) -> int:
        return self.size

    @remote.remote_method(batch_size=None)
    def get_capacity(self) -> int:
        return self.capacity

    @remote.remote_method(batch_size=None)
    def append(self, data: NestedTensor) -> None:
        self._append(data)

    @remote.remote_method(batch_size=None)
    def extend(self, data: Sequence[NestedTensor]) -> None:
        self._extend(data)

    @remote.remote_method(batch_size=None)
    def sample(self, batch_size: int) -> NestedTensor:
        index = np.random.randint(0, self.size, size=batch_size)
        return self.__getitem__(index)

    def _append(self, data: NestedTensor) -> int:
        return self._buffer.append(data)

    def _extend(self, data: Sequence[NestedTensor]) -> np.ndarray:
        return self._buffer.extend(data)

class PrioritizedReplayBuffer(ReplayBuffer):

    def __init__(self,
                 capacity: int,
                 alpha: float,
                 beta: float,
                 eps: float = 1e-8,
                 collate_fn: Optional[Callable[[Sequence[NestedTensor]],
                                               NestedTensor]] = None,
                 priority_type: np.dtype = np.float32) -> None:
        super().__init__(capacity, collate_fn)

        assert alpha > 0
        assert beta >= 0
        self._alpha = alpha
        self._beta = beta
        self._eps = eps

        self._priority_type = priority_type
        self._sum_tree = SumSegmentTree(capacity, dtype=priority_type)
        self._min_tree = MinSegmentTree(capacity, dtype=priority_type)

        self._max_priority = 1.0

    def __getitem__(
        self, index: Union[int, Tensor]
    ) -> Tuple[NestedTensor, Union[float, torch.Tensor]]:
        data = super().__getitem__(index)
        weight = self._compute_weight(index)
        return data, weight

    @property
    def alpha(self) -> float:
        return self._alpha

    @property
    def beta(self) -> float:
        return self._beta

    @property
    def eps(self) -> float:
        return self._eps

    @property
    def priority_type(self) -> np.dtype:
        return self._priority_type

    @property
    def max_priority(self) -> float:
        return self._max_priority

    @remote.remote_method(batch_size=None)
    def append(self,
               data: NestedTensor,
               priority: Optional[Union[float, Tensor]] = None) -> None:
        self._append(data, priority)

    @remote.remote_method(batch_size=None)
    def extend(self,
               data: Sequence[NestedTensor],
               priority: Optional[Union[float, Tensor]] = None) -> None:
        self._extend(data, priority)

    @remote.remote_method(batch_size=None)
    def sample(
            self,
            batch_size: int) -> Tuple[NestedTensor, torch.Tensor, torch.Tensor]:
        data, weight, index = self._sample(batch_size)
        return data, weight, torch.from_numpy(index)

    @remote.remote_method(batch_size=None)
    def update_priority(self, index: Union[int, Tensor],
                        priority: Union[float, Tensor]) -> None:
        self._update_priority(index, priority)

    def warm_up(self):
        capacity = self.get_capacity()
        width = len(str(capacity)) + 1
        cur_size = self.get_size()
        while cur_size < capacity:
            time.sleep(1)
            cur_size = self.get_size()
            logging.info("Warming up replay buffer: " +
                         f"[{cur_size: {width}d} / {capacity} ]")

    def _init_priority(self, index) -> None:
        priority = self._max_priority**self.alpha
        self._sum_tree[index] = priority
        self._min_tree[index] = priority

    def _update_priority(self, index: Union[int, Tensor],
                         priority: Union[float, Tensor]) -> None:
        priority += self.eps
        if isinstance(priority, float):
            self._max_priority = max(self._max_priority, priority)
        else:
            self._max_priority = max(self._max_priority, priority.max().item())
        priority = priority**self.alpha

        if isinstance(priority, np.ndarray):
            priority = priority.astype(self.priority_type)
        elif isinstance(priority, torch.Tensor):
            priority = priority.to(
                data_utils.numpy_dtype_to_torch(self.priority_type))

        self._sum_tree[index] = priority
        self._min_tree[index] = priority

    def _compute_weight(
            self, index: Union[int, Tensor]) -> Union[float, torch.Tensor]:
        p = self._sum_tree[index]
        if isinstance(p, np.ndarray):
            p = torch.from_numpy(p)
        p_min = self._min_tree.query(0, self.capacity)

        # Importance sampling weight formula:
        #   w_i = (p_i / sum(p) * N) ^ (-beta)
        #   weight_i = w_i / max(w)
        #   weight_i = (p_i / sum(p) * N) ^ (-beta) /
        #       ((min(p) / sum(p) * N) ^ (-beta))
        #   weight_i = ((p_i / sum(p) * N) / (min(p) / sum(p) * N)) ^ (-beta)
        #   weight_i = (p_i / min(p)) ^ (-beta)
        return (p / p_min)**(-self.beta)

    def _append(self,
                data: NestedTensor,
                priority: Optional[Union[float, Tensor]] = None) -> int:
        index = super()._append(data)
        if priority is None:
            self._init_priority(index)
        else:
            self._update_priority(index, priority)
        return index

    def _extend(self,
                data: Sequence[NestedTensor],
                priority: Optional[Union[float, Tensor]] = None) -> np.ndarray:
        index = super()._extend(data)
        if priority is None:
            self._init_priority(index)
        else:
            self._update_priority(index, priority)
        return index

    def _sample(
            self,
            batch_size: int) -> Tuple[NestedTensor, torch.Tensor, np.ndarray]:
        p_sum = self._sum_tree.query(0, self.capacity)
        mass = np.random.uniform(0.0, p_sum,
                                 size=batch_size).astype(self.priority_type)
        index = self._sum_tree.scan_lower_bound(mass)
        data, weight = self.__getitem__(index)
        return data, weight, index

class RemoteReplayBuffer(remote.Remote):

    def __init__(self,
                 target: ReplayBuffer,
                 server_name: str,
                 server_addr: str,
                 name: Optional[str] = None,
                 prefetch: int = 0,
                 timeout: float = 60) -> None:
        super().__init__(target, server_name, server_addr, name, timeout)
        self._prefetch = prefetch
        self._futures = collections.deque()

    @property
    def prefetch(self) -> Optional[int]:
        return self._prefetch

    def sample(
        self, batch_size: int
    ) -> Union[NestedTensor, Tuple[NestedTensor, torch.Tensor]]:
        if len(self._futures) > 0:
            ret = self._futures.popleft().result()
        else:
            ret = self.client.sync(self.server_name, "sample", batch_size)

        while len(self._futures) < self.prefetch:
            fut = self.client.async_(self.server_name, "sample", batch_size)
            self._futures.append(fut)

        return ret

    async def async_sample(
        self, batch_size: int
    ) -> Union[NestedTensor, Tuple[NestedTensor, torch.Tensor]]:
        if len(self._futures) > 0:
            ret = await self._futures.popleft()
        else:
            ret = await self.client.async_(self.server_name, "sample",
                                           batch_size)

        while len(self._futures) < self.prefetch:
            fut = self.client.async_(self.server_name, "sample", batch_size)
            self._futures.append(fut)

        return ret

    def warm_up(self):
        capacity = self.get_capacity()
        width = len(str(capacity)) + 1
        cur_size = self.get_size()
        while cur_size < capacity:
            time.sleep(1)
            cur_size = self.get_size()
            logging.info("Warming up replay buffer: " +
                         f"[{cur_size: {width}d} / {capacity} ]")

ReplayBufferLike = Union[ReplayBuffer, RemoteReplayBuffer]

def make_remote_replay_buffer(target: ReplayBuffer,
                              server: Server,
                              name: Optional[str] = None,
                              prefetch: int = 0,
                              timeout: float = 60):
    return RemoteReplayBuffer(target, server.name, server.addr, name, prefetch,
                              timeout)

max_stars_repo_name: stjordanis/rlmeta

max_issues_repo_name: stjordanis/rlmeta

max_forks_repo_name: stjordanis/rlmeta


text:
import timeit
import numpy as np
import fasterraster as fr
from pathlib import Path

NTESTS = 10

# Load a .bil file containing a DEM
fname = Path('./test_data/dem.bil')
dem = fr.read(fname)

# regular python implementation of hillshade function
# from https://www.neonscience.org/resources/learning-hub/tutorials/create-hillshade-py
def py_hillshade(dem, cell_size, azimuth=330, altitude=30):
    azimuth = 360.0 - azimuth

    dem = dem / cell_size
    x, y = np.gradient(dem)
    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
    aspect = np.arctan2(-x, y)
    azimuthrad = azimuth*np.pi/180.
    altituderad = altitude*np.pi/180.

    shaded = np.sin(altituderad)*np.sin(slope)
    + np.cos(altituderad)*np.cos(slope) * np.cos(
    (azimuthrad - np.pi/2.) - aspect)

    return 255*(shaded + 1)/2

# Time hillshade computation using regular python
time = timeit.timeit(lambda: py_hillshade(dem.raster, dem.XDIM), number=NTESTS)
print(f'python hillshade averaged {time/NTESTS:.3f} seconds')

# Time hillshade computation using open-MP for numt-threads
num_threads = [1,2,4,8]
for numt in num_threads:
    time = timeit.timeit(lambda: fr.hillshade_mp_faster(dem.raster, 330, 30, dem.XDIM, numt), number=10)
    print(f'hillshade averaged {time/NTESTS:.3f} seconds for {numt} threads')

max_stars_repo_name: asenogles/fasterraster

max_issues_repo_name: asenogles/fasterraster

max_forks_repo_name: asenogles/fasterraster