text:
from __future__ import division,print_function

from ..calcs import Calc

from numba import njit,guvectorize,vectorize,prange
import numba.cuda as cuda
from PIL import Image
from pyproj import Geod
from six import iteritems
import scipy.interpolate as interp
import scipy.signal as signal
import numpy as np
from tqdm import tqdm

__all__=["Scene","Renderer_35mm","Renderer_Composite","land_model"]

@njit
def _get_water(rs,ind):
    n_v = rs.shape[0]
    m = rs.shape[1]
    ind[:] = m
    for i in range(n_v):
        for j in range(m):
            if rs[i,j] <= 0:
                ind[i] = j
                break

@vectorize(["b1(f8,f8,f8)"])
def _get_bounds(a,a_min,a_max):
    return a >= a_min and a < a_max

@njit
def _get_bounds_block(a,a_mins,a_maxs,mask):
    n = a.shape[0]
    m = a_mins.shape[0]
    for i in range(m):
        mask[i] = False
        for j in range(n):
            if a[j] >= a_mins[i] and a[j] < a_maxs[i]:
                mask[i] = True
                break

@njit
def _get_bounds_sum(a,a_min,a_max):
    n = a.shape[0]
    count = 0
    for i in range(n):
        if a[i] >= a_min and a[i] < a_max:
            count += 1

    return count

@njit
def _get_vertical_mask(rh,h_min,h_max,inds,ds,d,sky,mask):
    n = rh.shape[0]

    for i in range(n):
        mask[i] = False
        if rh[i] >= h_min and rh[i] < h_max:
            if sky[i]: # if ray has hit sky (not hit anything) it has hit image
                mask[i] = True
            else:
                if ds[inds[i]] < d: # if ray hits land or water before image no change required
                    mask[i] = False
                else:
                    mask[i] = True

@cuda.jit("void(f8,f8[:,::1],f8[:],i4[:],b1[:],b1[:],b1[:])")
def _ray_crossing_gpu(h_min,rs,heights,inds,water,land,sky):
    i = cuda.grid(1)
    n_d = rs.shape[1]

    if i < rs.shape[0]:
        hit = False
        water[i] = False
        land[i] = False
        sky[i] = False
        inds[i] = -1
        for j in range(n_d):
            if rs[i,j] <= heights[j]:
                hit = True
                if heights[j] > h_min:
                    land[i] = True
                    inds[i] = j
                else:
                    water[i] = True
                    inds[i] = j

                break

        sky[i] = not hit

@njit(parallel=True)
def _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky):
    n_v = rs.shape[0]
    n_d = rs.shape[1]

    water[:] = False
    land[:] = False
    sky[:] = False
    inds[:] = -1

    for i in prange(n_v):
        hit = False

        for j in range(n_d):
            if rs[i,j] <= heights[j]:
                hit = True
                if heights[j] > h_min:
                    land[i] = True
                    inds[i] = j
                else:
                    water[i] = True
                    inds[i] = j

                break

        sky[i] = not hit

@njit
def _update_png_data(i,j,png_data,k,l,img_png_data):
    n = i.shape[0]
    m = j.shape[0]

    for a in range(n):
        ii = i[a]
        kk = k[a]
        for b in range(m):
            alpha = img_png_data[kk,l[b],3] / (255.0)
            png_data[ii,j[b],:] = (1.0 - alpha)*png_data[ii,j[b],:] + alpha*img_png_data[kk,l[b],:3]

@njit
def _update_png_data_slice(j,png_data,l,img_png_data):
    m = j.shape[0]

    for b in range(m):
       if img_png_data[l[b],3] > 0:
            alpha = img_png_data[l[b],3] / (255.0)
            png_data[j[b],:] =  (1.0 - alpha)*png_data[j[b],:] + alpha*img_png_data[l[b],:3]

def _check_gps(lat,lon):
    if np.any(lat < -90) or np.any(lat > 90):
        raise ValueError("latitude must be between -90 and 90")
    if np.any(lon < -180) or np.any(lon > 180):
        raise ValueError("lonitude must be between -180 and 180")

def _defualt_cfunc(d,heights,background_color):
    alpha = np.asarray(np.exp(-d/200000.0)) # exponential attenuation of light due to atmospheric scattering
    out = np.zeros(alpha.shape+(3,),dtype=np.uint8)

    out[:,0] = 50*alpha + background_color[0]*(1-alpha)
    out[:,1] = 150*alpha + background_color[1]*(1-alpha)
    out[:,2] = background_color[2]*(1-alpha)

    return out

def _render_cpu(png_data,h_min,rs,ds,h_angles,surface_color,background_color,terrain_args,image_args,disp=False):

    lat_obs,lon_obs,terrain,cfunc,cfunc_args = terrain_args
    img_datas,ray_heights = image_args

    n_v = rs.shape[0]
    n_z = rs.shape[1]

    water = np.zeros(n_v,dtype=np.bool)
    land = np.zeros(n_v,dtype=np.bool)
    sky = np.zeros(n_v,dtype=np.bool)
    inds = np.zeros(n_v,dtype=np.int32)

    v_mask = np.zeros(n_v,dtype=np.bool)
    h_mask = np.zeros(h_angles.shape[0],dtype=np.bool)

    h_mins = np.array([h_px.min() for _,h_px,_,_ in img_datas])
    h_maxs = np.array([h_px.max() for _,h_px,_,_ in img_datas])
    img_mask = np.zeros_like(h_mins,dtype=np.bool)

    if terrain.has_data: # render land model

        h_angle_max = h_angles.max()
        if disp:
            h_iter = tqdm(h_angles)
        else:
            h_iter = iter(h_angles)

        for i,h_angle in enumerate(h_iter):

            heights = terrain.get_terrain(lat_obs,lon_obs,h_angle,ds)

            _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky)

            png_data[i,water,:] = surface_color
            png_data[i,sky,:] = background_color

            if np.any(land):
                iland, = np.where(land)
                land_inds = inds[land]
                land_inds_m = land_inds-1
                land_inds_m[land_inds_m<0] = 0

                _heights = heights[land_inds]
                dx = ds[land_inds] - ds[land_inds_m]
                dy = (rs[iland,land_inds] - rs[iland,land_inds_m])
                dh = (_heights - heights[land_inds_m])

                t_ray = np.stack((dx,dy),axis=-1)
                t_land = np.stack((dx,dh),axis=-1)
                t_ray = (t_ray.T/np.linalg.norm(t_ray,axis=1))
                t_land = (t_land.T/np.linalg.norm(t_land,axis=1))
                n_ray = t_ray - 2*np.sum(t_ray*t_land,axis=0)*t_land

                try:
                    png_data[i,land,:] = cfunc(ds[land_inds],_heights,n_ray,*cfunc_args)
                except TypeError:
                    png_data[i,land,:] = cfunc(ds[land_inds],_heights,*cfunc_args)

            _get_bounds(h_angle,h_mins,h_maxs,out=img_mask)
            if np.any(img_mask):
                img_indx = np.argwhere(img_mask).ravel()
                for I in img_indx:
                    img_png_data,h_px,v_px,d = img_datas[I]

                    rh = ray_heights[d]
                    _get_vertical_mask(rh,v_px[0],v_px[-1],inds,ds,d,sky,v_mask)

                    j = np.argwhere(v_mask).ravel()
                    k = np.searchsorted(h_px,h_angle)
                    l = np.searchsorted(v_px,rh[j])
                    _update_png_data_slice(j,png_data[i,...],l,img_png_data[k,...])

    else: # if no terrain to render, just render sky and sphere surface

        _get_water(rs,inds)
        np.less(inds,n_z,out=water)
        png_data[:,water,:] = surface_color
        np.logical_not(water,out=sky)
        png_data[:,sky,:] = background_color

        _get_bounds_block(h_angles,h_mins,h_maxs,img_mask)
        if np.any(img_mask):
            img_indx = np.argwhere(img_mask).ravel()
            for I in img_indx:
                img_png_data,h_px,v_px,d = img_datas[I]

                rh = ray_heights[d]
                _get_bounds(h_angles,h_px[0],h_px[-1],out=h_mask)
                _get_vertical_mask(rh,v_px[0],v_px[-1],inds,ds,d,sky,v_mask)
                if np.any(v_mask):
                    i = np.argwhere(h_mask).ravel()
                    j = np.argwhere(v_mask).ravel()
                    k = np.searchsorted(h_px,h_angles[i])
                    l = np.searchsorted(v_px,rh[j])

                    _update_png_data(i,j,png_data,k,l,img_png_data)

def _render_gpu(png_data,h_min,rs,ds,h_angles,surface_color,background_color,terrain_args,image_args,disp=False):

    lat_obs,lon_obs,terrain,cfunc,cfunc_args = terrain_args
    img_datas,ray_heights = image_args

    if terrain.has_data: # render land model
        n_v = rs.shape[0]
        n_h = h_angles.shape[0]
        n_z = rs.shape[1]

        water = np.zeros(n_v,dtype=np.bool)
        land = np.zeros(n_v,dtype=np.bool)
        sky = np.zeros(n_v,dtype=np.bool)
        inds = np.zeros(n_v,dtype=np.int32)

        water_dev = cuda.to_device(water)
        land_dev = cuda.to_device(land)
        sky_dev = cuda.to_device(sky)
        inds_dev = cuda.to_device(inds)
        heights_dev = cuda.device_array_like(ds)
        ds_dev = cuda.to_device(ds)
        rs_dev = cuda.to_device(rs)

        stream = cuda.stream()

        v_mask = np.zeros(n_v,dtype=np.bool)
        h_mask = np.zeros(h_angles.shape[0],dtype=np.bool)

        h_mins = np.array([h_px.min() for _,h_px,_,_ in img_datas])
        h_maxs = np.array([h_px.max() for _,h_px,_,_ in img_datas])
        img_mask = np.zeros_like(h_mins,dtype=np.bool)

        nth = 1024
        nbk = max(n_v//nth,1)+1
        print(n_v,nth,nbk)

        h_angle_max = h_angles.max()

        if disp:
            h_iter = tqdm(h_angles)
        else:
            h_iter = iter(h_angles)

        for i,h_angle in enumerate(h_iter):
            heights = terrain.get_terrain(lat_obs,lon_obs,h_angle,ds)

            cuda.to_device(heights,to=heights_dev)

            _ray_crossing_gpu[nbk,nth](h_min,rs_dev,heights_dev,inds_dev,water_dev,land_dev,sky_dev)

            water_dev.copy_to_host(ary=water, stream=stream)
            inds_dev.copy_to_host(ary=inds, stream=stream)
            land_dev.copy_to_host(ary=land, stream=stream)
            sky_dev.copy_to_host(ary=sky, stream=stream)

            png_data[i,water,:] = surface_color
            png_data[i,sky,:] = background_color

            if np.any(land):
                iland, = np.where(land)
                land_inds = inds[land]
                land_inds_m = land_inds-1
                land_inds_m[land_inds_m<0] = 0

                _heights = heights[land_inds]
                dx = ds[land_inds] - ds[land_inds_m]
                dy = (rs[iland,land_inds] - rs[iland,land_inds_m])
                dh = (_heights - heights[land_inds_m])

                t_ray = np.stack((dx,dy),axis=-1)
                t_land = np.stack((dx,dh),axis=-1)
                t_ray = (t_ray.T/np.linalg.norm(t_ray,axis=1))
                t_land = (t_land.T/np.linalg.norm(t_land,axis=1))
                n_ray = t_ray - 2*np.sum(t_ray*t_land,axis=0)*t_land

                try:
                    png_data[i,land,:] = cfunc(ds[land_inds],_heights,n_ray,*cfunc_args)
                except TypeError:
                    png_data[i,land,:] = cfunc(ds[land_inds],_heights,*cfunc_args)

            _get_bounds(h_angle,h_mins,h_maxs,out=img_mask)
            if np.any(img_mask):
                img_indx = np.argwhere(img_mask).ravel()
                for I in img_indx:
                    img_png_data,h_px,v_px,d = img_datas[I]

                    rh = ray_heights[d]
                    _get_vertical_mask(rh,v_px[0],v_px[-1],inds,ds,d,sky,v_mask)

                    j = np.argwhere(v_mask).ravel()
                    k = np.searchsorted(h_px,h_angle)
                    l = np.searchsorted(v_px,rh[j])
                    _update_png_data_slice(j,png_data[i,...],l,img_png_data[k,...])

    else: # if no terrain to render, just render sky and sphere surface

        n_v = rs.shape[0]
        n_z = rs.shape[1]

        water = np.zeros(n_v,dtype=np.bool)
        sky = np.zeros(n_v,dtype=np.bool)
        inds = np.zeros(n_v,dtype=np.int32)

        v_mask = np.zeros(n_v,dtype=np.bool)
        h_mask = np.zeros(h_angles.shape[0],dtype=np.bool)

        h_mins = np.array([h_px.min() for _,h_px,_,_ in img_datas])
        h_maxs = np.array([h_px.max() for _,h_px,_,_ in img_datas])
        img_mask = np.zeros_like(h_mins,dtype=np.bool)

        nth = 1024
        nbk = max(n_v//nth,1)

        _get_water(rs,inds)
        np.less(inds,n_z,out=water)
        png_data[:,water,:] = surface_color
        np.logical_not(water,out=sky)
        png_data[:,sky,:] = background_color

        _get_bounds_block(h_angles,h_mins,h_maxs,img_mask)
        if np.any(img_mask):
            img_indx = np.argwhere(img_mask).ravel()
            for I in img_indx:
                img_png_data,h_px,v_px,d = img_datas[I]

                rh = ray_heights[d]
                _get_bounds(h_angles,h_px[0],h_px[-1],out=h_mask)
                _get_vertical_mask(rh,v_px[0],v_px[-1],inds,ds,d,sky,v_mask)
                if np.any(v_mask):
                    i = np.argwhere(h_mask).ravel()
                    j = np.argwhere(v_mask).ravel()
                    k = np.searchsorted(h_px,h_angles[i])
                    l = np.searchsorted(v_px,rh[j])

                    _update_png_data(i,j,png_data,k,l,img_png_data)

def _prep_scene(scene,h_angles,lat_obs,lon_obs,geod,sol):
    img_datas = []
    ray_heights = {}
    dh_angle = np.diff(h_angles).min()

    for image,(im,image_data) in iteritems(scene._image_dict):
        px_width,px_height = im.size

        for lat,lon,img_height,width,height,heading in image_data:
            horz_pixel_pos = np.linspace(-width/2,width/2,px_width)
            v_px = img_height + np.linspace(0,height,px_height)

            f_az,b_az,dist = geod.inv(lon_obs,lat_obs,lon,lat)

            if heading is None:
                alpha = 0.0
            else:
                b_az = b_az%360
                heading = heading%360
                alpha = np.deg2rad(heading - b_az)

            h_px=np.arctan(horz_pixel_pos*np.abs(np.cos(alpha))/dist)

            np.rad2deg(h_px,out=h_px)
            h_px += (f_az % 360)
            np.mod(h_px,360,out=h_px)

            if dist not in ray_heights:
                rh = sol(dist)
                ray_heights[dist] = rh
            else:
                rh = ray_heights[dist]

            n_h = max(int((h_px[-1]-h_px[0])/dh_angle),2)
            n_v = max(_get_bounds_sum(rh,v_px[0],v_px[-1]),2)

            try:
                new_im = im.resize((n_h,n_v),Image.LANCZOS)
            except ValueError:
                continue
            im_data = np.array(new_im)
            im_data = im_data[::-1,:,:].transpose((1,0,2)).copy()
            h_px = np.linspace(h_px[0],h_px[-1],n_h)
            v_px = np.linspace(v_px[0],v_px[-1],n_v)
            img_datas.append((im_data,h_px,v_px,dist))

    img_datas.sort(key=lambda x:-x[-1])

    return img_datas,ray_heights

@njit
def is_sorted(a):
    for i in range(a.size-1):
         if a[i+1] < a[i] :
               return False
    return True

class Renderer_35mm(object):
    """
    Object used to set a camera angle and position to render a scene.
    """
    def __init__(self,calc,h_obs,lat_obs,lon_obs,direction,max_distance,
                 distance_res=10,vert_obs_angle=0.0,vert_res=1000,
                 focal_length=2000,atol=1.1e-7,rtol=1.1e-7,aspect_ratio=None,ellps="sphere"):

        """
        Parameters
        ----------

        calc: calcs object
            object which can be used to calculate the trajectory of light rays

        h_obs: float
            height of renderer in meters

        lat_obs: float
            latitude of renderer in degrees between -90 and 90

        lon_obs: float
            longitude of renderer in degrees between -180 adn 180

        direction: tuple or float
            tuple: tuple containing a latitude and longitude point to point renderer towards.
            float: heading to point the renderer in.

        max_distance: float
            maximum distance away from the renderer to calculate the light ray path

        distance_res: float
            distance between points where the light ray path is checked for intersection point.

        vert_obs_angle: float
            vertical tilt of the renderer.

        vert_res: int
            resolution of pixels on the verticle

        focal_length: float
            focal length in milimeters to calculate the field of view of the camera. This, along with the vertial resolution, sets the resolution on the horizontal.

        atol: float
            absolute tolerance of ode solver

        rtol: float
            relative tolerance of the ode solver

        aspect_ratio: tuple
            aspect ratio for image

        ellps: str,
            type of geoid to use for calculating gps, default is "sphere"

        Note
        ----
        The focal length sets the verticle field of view of the frame using the standard 24 mm for 35 mm film.

        """
        if not isinstance(calc,Calc):
            raise ValueError("expecting calculator to be instance of Calc base.")

        self._calc = calc
        self._geod = Geod(ellps=ellps)

        self._atol=atol
        self._rtol=rtol

        self._lat_obs = float(lat_obs)
        self._lon_obs = float(lon_obs)
        _check_gps(self._lat_obs,self._lon_obs)
        if type(direction) is tuple:
            if len(direction) != 2:
                raise ValueError("direction must be either heading or tuple containing latitude and lonitude respectively.")

            lat_dir,lon_dir = direction
            _check_gps(lat_dir,lon_dir)
            f_az,b_az,dist = self._geod.inv(lon_obs,lat_obs,lon_dir,lat_dir)

        else:
            f_az = float(direction)

        self._h_obs = float(h_obs)

        f_az = f_az%360

        vert_res = int(vert_res)

        if aspect_ratio is None:
            dx=18
            horz_res = int(vert_res*1.5)
        else:
            dx = float(aspect_ratio[0])*12.0/float(aspect_ratio[1])
            r = float(aspect_ratio[0])/aspect_ratio[1]
            horz_res = int(vert_res*r)

        self._vert_res = vert_res
        self._horz_res = horz_res
        self._focal_length = focal_length

        y_grid = np.linspace(-12,12,vert_res)
        x_grid = np.linspace(-dx,dx,horz_res)
        self._v_angles = np.rad2deg(np.arctan(y_grid/focal_length))+vert_obs_angle
        self._h_angles = np.rad2deg(np.arctan(x_grid/focal_length))+f_az
        self._vfov = self._v_angles.max()-self._v_angles.min()

        self._ds = np.arange(0.0,max_distance,distance_res,dtype=np.float64)
        sol = self._calc.solve_ivp(max_distance,h_obs,alpha=self._v_angles,atol=atol,rtol=rtol,dense_output=True)
        self._rs = np.ascontiguousarray(sol.sol(self._ds)[:vert_res])

        self._sol = interp.interp1d(self._ds,self._rs,axis=-1,copy=False,bounds_error=True)

    @property
    def vfov(self):
        """Vertial field of view for this particular renderer"""
        return self._vfov

    @property
    def v_angles(self):
        """Vertical angular scale of the image frame."""
        v_angles = self._v_angles[...]
        v_angles.setflags(write=False)
        return v_angles

    @property
    def h_angles(self):
        """Horizontal angular scale of the image frame."""
        h_angles = self._h_angles[...]
        h_angles.setflags(write=False)
        return h_angles

    @property
    def calc(self):
        """Calc object used to calculate verticle rays."""
        return self._calc

    @property
    def geod(self):
        """geod object used to distances between latitude and longitude coordinates."""
        return self._geod

    @property
    def ds(self):
        """distance points."""
        ds = self._ds[...]
        ds.setflags(write=False)
        return ds

    def set_location(self,lat_obs,lon_obs,direction):
        """This function can be used to change the location and heading of the renderer.

        Parameters
        ----------
        lat_obs: float
            new latitude of observer
        lon_obs: float
            new longitude of observer
        direction: tuple or float
            tuple: tuple containing a latitude and longitude point to point renderer towards.
            float: heading to point the renderer in.

        """
        self._lat_obs = float(lat_obs)
        self._lon_obs = float(lon_obs)
        _check_gps(self._lat_obs,self._lon_obs)
        if type(direction) is tuple:
            if len(direction) != 2:
                raise ValueError("direction must be either heading or tuple containing latitude and lonitude respectively.")

            lat_dir,lon_dir = direction
            _check_gps(lat_dir,lon_dir)
            f_az,b_az,dist = self._geod.inv(lon_obs,lat_obs,lon_dir,lat_dir)

        else:
            f_az = float(direction)

        f_az = f_az%360

        x_grid = np.linspace(-18,18,self._horz_res)
        self._h_angles = np.rad2deg(np.arctan(x_grid/self._focal_length))+f_az

    def ray_diagram(self,ax,heading=None,heights=None,style="sphere_top",
                    eye_level=True,linewidth_rays=0.2,linewidth_earth=0.2,R0=6371008,h_min=0.01):
        """Side profile view showing the rays trajectories.

        Parameters
        ----------
        ax: matplotlib Axes object
            this function plots lines on a graph, you can specify which graph to plot to by passing in the axes object for the graph you want to plot to.

        heading: scalar, optional
            the heading to take slice of tarrain for ray paths. default is the center of the frame for renderer

        heights: array_like, (N,), optional
            list of values which represent the elevation profile along ray trajectory (in addition to elevation)

        style: str, optional
            style to plot the graph. "flat": plot rays on flat plane, "sphere_top": plot rays as if the earth is filling away from the observer, "sphere_side": plot rays with the 'bulge' in the middle.

        eye_level: bool, optional
            color the ray which is closest to eye level (0 degrees initial angle) orange.

        linewidth_rays: float, optional
            linewidth used for plotting rays in diagram

        linewidth_earth: float, optional
            linewidth used for plotting surface of the earth, both water and land.

        h_min: float, optional
            minimum values which should count as water when calculating hit locations.

        """
        d = self._ds

        if len(d) == 0 or not is_sorted(d):
            raise ValueError("array 'd' must contain distance values in ascending order.")

        d_max = self._ds.max()

        if heading is None:
            h_angle = self.h_angles.mean()
        else:
            h_angle = heading

        rs = self._rs
        if heights is None:
            heights = np.zeros_like(d)

        n_v = self.v_angles.size

        water = np.zeros(n_v,dtype=np.bool)
        land = np.zeros(n_v,dtype=np.bool)
        sky = np.zeros(n_v,dtype=np.bool)
        inds = np.zeros(n_v,dtype=np.int32)

        i_horz = np.abs(self.v_angles).argmin()

        if style == "sphere_top":
            c = np.cos(np.pi/2-d/R0)
            s = np.sin(np.pi/2-d/R0)

            _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky)
            rs = rs + R0
            heights = heights + R0
            for i in range(n_v):
                if i == i_horz and eye_level:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0,color="orange",linewidth=linewidth_rays)
                    continue

                if water[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0,color="blue",linewidth=linewidth_rays)
                    continue

                if land[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0,color="green",linewidth=linewidth_rays)
                    continue

                if sky[i]:
                    ax.plot(rs[i,:]*c,rs[i,:]*s-R0,color="cyan",linewidth=linewidth_rays)
                    continue

            ax.plot(c*heights,s*heights-R0,color="green",linewidth=linewidth_earth)
            ax.plot(c*R0,s*R0-R0,color="blue",linewidth=linewidth_earth)
        elif style == "sphere_side":
            c = np.cos(np.pi/2-d/R0+d_max/(2*R0))
            c -= c.min()
            s = np.sin(np.pi/2-d/R0+d_max/(2*R0))
            _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky)
            rs = rs + R0
            heights = heights + R0
            for i in range(n_v):
                if i == i_horz and eye_level:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0*s.min(),color="orange",linewidth=linewidth_rays)
                    continue

                if water[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0*s.min(),color="blue",linewidth=linewidth_rays)
                    continue

                if land[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0*s.min(),color="green",linewidth=linewidth_rays)
                    continue

                if sky[i]:
                    ax.plot(rs[i,:]*c,rs[i,:]*s-R0*s.min(),color="cyan",linewidth=linewidth_rays)
                    continue

            ax.plot(c*heights,s*heights-R0*s.min(),color="green",linewidth=linewidth_earth)
            ax.plot(c*R0,s*R0-R0*s.min(),color="blue",linewidth=linewidth_earth)
        elif style == "flat":
            _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky)

            for i in range(n_v):
                if i == i_horz and eye_level:
                    i_max = inds[i]
                    ax.plot(d[:i_max],rs[i,:i_max],color="orange",linewidth=linewidth_rays)
                    continue

                if water[i]:
                    i_max = inds[i]
                    ax.plot(d[:i_max],rs[i,:i_max],color="blue",linewidth=linewidth_rays)
                    continue

                if land[i]:
                    i_max = inds[i]
                    ax.plot(d[:i_max],rs[i,:i_max],color="green",linewidth=linewidth_rays)
                    continue

                if sky[i]:
                    ax.plot(d,rs[i,:],color="cyan",linewidth=linewidth_rays)
                    continue

            ax.plot(d,heights,color="green",linewidth=linewidth_earth)
            ax.plot(d,np.zeros_like(d),color="blue",linewidth=linewidth_earth)
        else:
            raise ValueError

    def render_scene(self,scene,image_name,surface_color=None,background_color=None,cfunc=None,cfunc_args=None,
        disp=False,eye_level=False,postprocess=None,h_min=0.01,gpu=False):
        """Render a scene object for the renderer's given field of view and direction.

        Parameters
        ----------
        scene: Scene object
            object which contains data which the renderer can extract and render

        image_name: str
            name for the image

        background_color: array_like, optional
            3 element array containing an color via RGB color code (numbered 0 to 255)
            default value: `[135,206,250]`

        surface_color: array_like, optional
            3 element array containing an color via RGB color code (numbered 0 to 255)
            default value: `[0,80,120]`

        cfunc: callable, optional
            Callable function which is used to color the elevation data. The function takes
            in two arrays, the first is distances from the observer    the second is the elevation
            above the surface of the earth.

        cfunc_args: array_like, optional
            extra arguments to pass into `cfunc`

        disp: bool, optional
            when rendering topographical data this will print out the heading slice which has been rendered.

        eye_level: bool, optional
            when rendering the image, an orange line is placed at eye level in the image.

        post_process: callable, optional
            function which processes the final image before saving it.

        h_min: float, optional
            minimum value for ray to count as crossing water.

        gpu: bool, optional
            utilize GPU to calculate 3D land model rendering.

        """
        if surface_color is None:
            surface_color = np.array([0,80,120],dtype=np.uint8)
        else:
            surface_color = np.fromiter(surface_color,dtype=np.uint8)

        if background_color is None:
            background_color = np.array([135,206,250],dtype=np.uint8)
        else:
            background_color = np.fromiter(background_color,dtype=np.uint8)

        if cfunc is None and cfunc_args is None:
            cfunc = _defualt_cfunc
            cfunc_args = (background_color,)
        elif cfunc is None and cfunc_args is not None:
            raise Exception("cfunc_args given without cfunc.")

        if cfunc_args is None:
            cfunc_args = ()

        img_datas,ray_heights = _prep_scene(scene,self._h_angles,self._lat_obs,self._lon_obs,self._geod,self._sol)

        land_model = scene._land_model

        png_data = np.empty((len(self._h_angles),self._rs.shape[0],3),dtype=np.uint8)
        png_data[...] = 0

        terrain_args = (self._lat_obs,self._lon_obs,land_model,cfunc,cfunc_args)
        image_args = (img_datas,ray_heights)

        if gpu:
            _render_gpu(png_data,h_min,self._rs,self._ds,self._h_angles,surface_color,background_color,terrain_args,image_args,disp)
        else:
            _render_cpu(png_data,h_min,self._rs,self._ds,self._h_angles,surface_color,background_color,terrain_args,image_args,disp)

        if eye_level:
            i_horz = np.argwhere(np.abs(self._v_angles)<(self._vfov/800.0)).ravel()
            png_data[:,i_horz,0] = 255
            png_data[:,i_horz,1] = 100
            png_data[:,i_horz,2] = 0

        png_data = png_data.transpose((1,0,2))
        png_data = png_data[::-1,:,:]
        im = Image.fromarray(png_data,mode="RGB")

        if postprocess is not None:
            im = postprocess(im)

        im.save(image_name)

class Renderer_Composite(object):
    """
    Object used to set a camera angle and position to render a scene for a user specified horizontal field of view.
    """
    def __init__(self,calc,h_obs,lat_obs,lon_obs,max_distance,
                 distance_res=10,vert_obs_angle=0.0,vert_res=1000,
                 focal_length=2000,atol=1.1e-7,rtol=1.1e-7,ellps="sphere"):
        """
        Parameters
        ----------

        calc: calcs object
            object which can be used to calculate the trajectory of light rays

        h_obs: float
            height of renderer in meters

        lat_obs: float
            latitude of renderer in degrees between -90 and 90

        lon_obs: float
            longitude of renderer in degrees between -180 adn 180

        max_distance: float
            maximum distance away from the renderer to calculate the light ray path

        distance_res: float
            distance between points where the light ray path is checked for intersection point.

        vert_obs_angle: float
            vertical tilt of the renderer.

        vert_res: int
            resolution of pixels on the verticle

        focal_length: float
            focal length in milimeters to calculate the field of view of the camera. This, along with the vertial resolution, sets the resolution on the horizontal.

        atol: float
            absolute tolerance of ode solver

        rtol: float
            relative tolerance of the ode solver

        """

        if not isinstance(calc,Calc):
            raise ValueError("expecting calculator to be instance of Calc base.")

        self._calc = calc
        self._atol=atol
        self._rtol=rtol

        self._lat_obs = float(lat_obs)
        self._lon_obs = float(lon_obs)
        self._h_obs = float(h_obs)
        _check_gps(self._lat_obs,self._lon_obs)

        self._calc = calc
        self._geod = Geod(ellps=ellps)

        vert_res = int(vert_res)
        vert_min = np.rad2deg(np.arctan(-12.0/focal_length))+vert_obs_angle
        vert_max = np.rad2deg(np.arctan( 12.0/focal_length))+vert_obs_angle

        self._v_angles,self._dangles = np.linspace(vert_min,vert_max,vert_res,retstep=True)
        self._vfov = self._v_angles.max() - self._v_angles.min()

        self._ds = np.arange(0.0,max_distance,distance_res,dtype=np.float64)
        sol = self._calc.solve_ivp(max_distance,h_obs,alpha=self._v_angles,atol=atol,rtol=rtol,dense_output=True)
        self._rs = sol.sol(self._ds)[:vert_res].copy()

        self._sol = interp.interp1d(self._ds,self._rs,axis=-1,copy=False,bounds_error=True)

    @property
    def vfov(self):
        """Vertial field of view for this particular renderer"""
        return self._vfov

    @property
    def v_angles(self):
        """Vertical angular scale of the image frame."""
        v_angles = self._v_angles[...]
        v_angles.setflags(write=False)
        return v_angles

    @property
    def calc(self):
        """Calc object used to calculate verticle rays."""
        return self._calc

    @property
    def geod(self):
        """geod object used to distances between latitude and longitude coordinates."""
        return self._geod

    @property
    def ds(self):
        """distance points."""
        ds = self._ds[...]
        ds.setflags(write=False)
        return ds

    def set_location(self,lat_obs,lon_obs):
        """This function can be used to change the position of the renderer.

        Parameters
        ----------
        lat_obs: float
            new latitude of observer
        lon_obs: float
            new longitude of observer

        """
        self._lat_obs = float(lat_obs)
        self._lon_obs = float(lon_obs)
        _check_gps(self._lat_obs,self._lon_obs)

    def ray_diagram(self,ax,heading,heights=None,style="sphere_top",
                    eye_level=True,linewidth_rays=0.2,linewidth_earth=0.2,R0=6371008,h_min=0.01):
        """Side profile view showing the rays trajectories.

        Parameters
        ----------
        ax: matplotlib Axes object
            this function plots lines on a graph, you can specify which graph to plot to by passing in the axes object for the graph you want to plot to.

        heading: scalar
            the heading to take slice of tarrain for ray paths.

        heights: array_like, (N,), optional
            list of values which represent the elevation profile along ray trajectory (in addition to elevation)

        style: str, optional
            style to plot the graph. "flat": plot rays on flat plane, "sphere_top": plot rays as if the earth is filling away from the observer, "sphere_side": plot rays with the 'bulge' in the middle.

        eye_level: bool, optional
            color the ray which is closest to eye level (0 degrees initial angle) orange.

        linewidth_rays: float, optional
            linewidth used for plotting rays in diagram

        linewidth_earth: float, optional
            linewidth used for plotting surface of the earth, both water and land.

        h_min: float, optional
            minimum values which should count as water when calculating hit locations.

        """
        d = self._ds

        if len(d) == 0 or not is_sorted(d):
            raise ValueError("array 'd' must contain distance values in ascending order.")

        d_max = self._ds.max()

        h_angle = heading

        rs = self._rs
        if heights is None:
            heights = np.zeros_like(d)

        n_v = self.v_angles.size

        water = np.zeros(n_v,dtype=np.bool)
        land = np.zeros(n_v,dtype=np.bool)
        sky = np.zeros(n_v,dtype=np.bool)
        inds = np.zeros(n_v,dtype=np.int32)

        i_horz = np.abs(self.v_angles).argmin()

        if style == "sphere_top":
            c = np.cos(np.pi/2-d/R0)
            s = np.sin(np.pi/2-d/R0)

            _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky)
            rs = rs + R0
            heights = heights + R0
            for i in range(n_v):
                if i == i_horz and eye_level:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0,color="orange",linewidth=linewidth_rays)
                    continue

                if water[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0,color="blue",linewidth=linewidth_rays)
                    continue

                if land[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0,color="green",linewidth=linewidth_rays)
                    continue

                if sky[i]:
                    ax.plot(rs[i,:]*c,rs[i,:]*s-R0,color="cyan",linewidth=linewidth_rays)
                    continue

            ax.plot(c*heights,s*heights-R0,color="green",linewidth=linewidth_earth)
            ax.plot(c*R0,s*R0-R0,color="blue",linewidth=linewidth_earth)
        elif style == "sphere_side":
            c = np.cos(np.pi/2-d/R0+d_max/(2*R0))
            c -= c.min()
            s = np.sin(np.pi/2-d/R0+d_max/(2*R0))
            _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky)
            rs = rs + R0
            heights = heights + R0
            for i in range(n_v):
                if i == i_horz and eye_level:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0*s.min(),color="orange",linewidth=linewidth_rays)
                    continue

                if water[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0*s.min(),color="blue",linewidth=linewidth_rays)
                    continue

                if land[i]:
                    i_max = inds[i]
                    ax.plot(rs[i,:i_max]*c[:i_max],
                            rs[i,:i_max]*s[:i_max]-R0*s.min(),color="green",linewidth=linewidth_rays)
                    continue

                if sky[i]:
                    ax.plot(rs[i,:]*c,rs[i,:]*s-R0*s.min(),color="cyan",linewidth=linewidth_rays)
                    continue

            ax.plot(c*heights,s*heights-R0*s.min(),color="green",linewidth=linewidth_earth)
            ax.plot(c*R0,s*R0-R0*s.min(),color="blue",linewidth=linewidth_earth)
        elif style == "flat":
            _ray_crossing_cpu(h_min,rs,heights,inds,water,land,sky)

            for i in range(n_v):
                if i == i_horz and eye_level:
                    i_max = inds[i]
                    ax.plot(d[:i_max],rs[i,:i_max],color="orange",linewidth=linewidth_rays)
                    continue

                if water[i]:
                    i_max = inds[i]
                    ax.plot(d[:i_max],rs[i,:i_max],color="blue",linewidth=linewidth_rays)
                    continue

                if land[i]:
                    i_max = inds[i]
                    ax.plot(d[:i_max],rs[i,:i_max],color="green",linewidth=linewidth_rays)
                    continue

                if sky[i]:
                    ax.plot(d,rs[i,:],color="cyan",linewidth=linewidth_rays)
                    continue

            ax.plot(d,heights,color="green",linewidth=linewidth_earth)
            ax.plot(d,np.zeros_like(d),color="blue",linewidth=linewidth_earth)
        else:
            raise ValueError

    def render_scene(self,scene,image_names,heading_mins,heading_maxs,surface_color=None,background_color=None,
        cfunc=None,cfunc_args=None,disp=False,eye_level=False,postprocess=None,h_min=0.01,gpu=False):
        """Renders a composites over a very wide horizontal field.

        Parameters
        ----------
        scene: Scene object
            object which contains data which the renderer can extract and render

        image_names: array_like (n,) or 'str'
            name(s) for the image(s) being rendered.

        heading_mins: array_like (n,) or float
            minimum starting point(s) for composite.

        heading_maxs: array_like (n,) or float
            maximum starting point(s) for composite, pair with 'heading_mins'.

        background_color: array_like, optional
            3 element array containing an color via RGB color code (numbered 0 to 255)
            default value: `[135,206,250]`

        surface_color: array_like, optional
            3 element array containing an color via RGB color code (numbered 0 to 255)
            default value: `[0,80,120]`

        cfunc: callable, optional
            Callable function which is used to color the elevation data. The function takes
            in two arrays, the first is distances from the observer    the second is the elevation
            above the surface of the earth.

        cfunc_args: array_like, optional
            extra arguments to pass into `cfunc`

        disp: bool, optional
            when rendering topographical data this will print out the heading slice which has been rendered.

        eye_level: bool, optional
            when rendering the image, an orange line is placed at eye level in the image.

        post_process: callable, optional
            function which processes the final image before saving it.

        h_min: float, optional
            minimum value for ray to count as crossing water.

        gpu: bool, optional
            utilize GPU to calculate 3D land model rendering.

        """
        if surface_color is None:
            surface_color = np.array([0,80,120])
        else:
            surface_color = np.fromiter(surface_color,dtype=np.uint8)

        if background_color is None:
            background_color = np.array([135,206,250])
        else:
            background_color = np.fromiter(background_color,dtype=np.uint8)

        if cfunc is None and cfunc_args is None:
            cfunc = _defualt_cfunc
            cfunc_args = (background_color,)
        elif cfunc is None and cfunc_args is not None:
            raise Exception("cfunc_args given without cfunc.")

        if cfunc_args is None:
            cfunc_args = ()

        n_v = self._rs.shape[0]

        land_model = scene._land_model

        heading_mins = np.atleast_1d(heading_mins).ravel()
        heading_maxs = np.atleast_1d(heading_maxs).ravel()
        image_names =  np.atleast_1d(image_names).ravel()
        tup = np.broadcast_arrays(heading_mins,heading_maxs,image_names)

        for heading_min,heading_max,image_name in zip(*tup):
            print(heading_min,heading_max,image_name)
            h_angles = np.arange(heading_min,heading_max,self._dangles)

            img_datas,ray_heights = _prep_scene(scene,h_angles,self._lat_obs,self._lon_obs,self._geod,self._sol)

            png_data = np.empty((len(h_angles),n_v,3),dtype=np.uint8)
            png_data[...] = 0

            terrain_args = (self._lat_obs,self._lon_obs,land_model,cfunc,cfunc_args)
            image_args = (img_datas,ray_heights)

            if gpu:
                _render_gpu(png_data,h_min,self._rs,self._ds,h_angles,surface_color,background_color,terrain_args,image_args,disp)
            else:
                _render_cpu(png_data,h_min,self._rs,self._ds,h_angles,surface_color,background_color,terrain_args,image_args,disp)

            if eye_level:
                i_horz = np.argwhere(np.abs(self._v_angles)<(self._vfov/800.0)).ravel()
                png_data[:,i_horz,0] = 255
                png_data[:,i_horz,1] = 100
                png_data[:,i_horz,2] = 0

            png_data = png_data.transpose((1,0,2))
            png_data = png_data[::-1,:,:]
            im = Image.fromarray(png_data,mode="RGB")
            if postprocess is not None:
                im = postprocess(im)
            im.save(image_name)

class Scene(object):
    """
    Simple wrapper which keeps track of data used to render an image.
    """
    def __init__(self):
        """
        This function initialized `Scene` object.
        """
        self._land_model = land_model()
        self._image_dict = {}

    @property
    def land_model(self):
        """Object used to generate elevation along great circle."""
        return self._land_model

    def add_elevation_model(self,*args):
        """Add terrain data to the interpolated model.

        Parameters
        ----------
        args: tuple
            tuple which contains elevation data:
            if len(args) == 3: args = (lats,lons,elevation) which contains the arguments for scipy.interpolate.RegularGridInterpolator.
            if len(args) == 2: args = (points,elevation) which contains the arguments for scipy.interpolate.LinearNDInterpolator.
            if len(args) == 1: args = object that is instance of `land_model`.

        """
        if len(args)==1 and isinstance(args[0],land_model):
            self._land_model += args[0]
        else:
            self._land_model.add_elevation_data(*args)

    def add_image(self,image,image_pos,dimensions,direction=None):
        """Add image to scene.

        Parameters
        ----------
        image: str
            string which contains path to image file.

        image_pos: array_like
            either has (h_obj,lat,lon) or (lat,lon) h_obj is height above the earth's surface

        dimensions: tuple
            contains dimensions of the image is in meters. If either one has value `-1`
            then that dimension is determined by the resolution of the picture.
        """
        if type(image) is str:
            if image in self._image_dict:
                im = self._image_dict[image][0]
            else:
                im = Image.open(image)
                self._image_dict[image] = (im,[])
        else:
            raise ValueError("image must be a string containing name of file.")

        px_width,px_height = im.size

        aspect = float(px_height)/px_width

        if len(image_pos) == 2:
            lat,lon = image_pos
            h = self._land_model(lat,lon)
        elif len(image_pos) == 3:
            h,lat,lon = image_pos
        else:
            raise ValueError("expecting image_pos to contain gps coordinates and optionally height of image.")

        if direction is None:
            heading = None
        elif type(direction) is tuple:
            if len(direction) != 2:
                raise ValueError("direction must be either heading or tuple containing latitude and lonitude respectively.")

            lat_dir,lon_dir = direction
            _check_gps(lat_dir,lon_dir)
            heading,b_az,dist = self._geod.inv(lon_obs,lat_obs,lon_dir,lat_dir)

        else:
            heading = float(direction)

        try:
            width,height = dimensions
        except ValueError:
            raise ValueError("dimensions must contain only length and width of object.")

        if width == -1 and height == -1:
            raise ValueError("at least one height must be specifided to deduce size.")

        if width == -1:
            width = height/aspect

        if height == -1:
            height = width*aspect

        self._image_dict[image][1].append((lat,lon,h,width,height,heading))

class land_model(object):
    """
    This class is used to grab slices of the tarrain alon great circle
    slices of the earth.
    """
    def __init__(self,ellps='sphere'):
        """
        Simple wrapper which does the bilinear interpolation of terrain data.

        Parameters
        ----------
        ellps: str, optional
            String representing which pyproj datum to use for the data being interpolated.

        """
        self._geod = Geod(ellps=ellps)
        self._terrain_list = []
        self._d_min = np.inf
        self._d_max = 0
        self._h_max = 0

    def get_terrain(self,lat,lon,heading,dist):
        """
        Sample terrain along geodesic defined by an intial latitude and longitude point along a heading at different distances.

        Parameters
        ----------
        lat: array_like
            latitude(s) which starting point(s) begin at
        lon: array_like
            longitude(s) which starting point(s) begin at
        heading: array_like
            heading(s) which determine which direction to go on the spheriod
        dist: array_like
            distance(s) in meters of how far to go along the heading

        Notes
        -----
        This function follows NumPy's broadcasting rules for handing arguments of different shapes.

        """
        lon,lat,az = self._geod.fwd(*np.broadcast_arrays(lon,lat,heading,dist))

        return self.__call__(lat,lon)

    def __call__(self,lat,lon):
        lat = np.atleast_1d(lat)
        lon = np.atleast_1d(lon)
        heights = np.zeros_like(lon)

        coors = np.stack(np.broadcast_arrays(lat,lon),axis=-1)

        for terrain in self._terrain_list:
            try:
                heights += terrain(coors)
            except TypeError:
                heights += terrain.ev(*np.broadcast_arrays(lat,lon))

        return np.squeeze(heights)

    @property
    def has_data(self):
        """flag if True the terrain object has land else it has no land."""
        return len(self._terrain_list) > 0

    def add_elevation_data(self,*args,smooth=False):
        """Add terrain data to the interpolated model.

        Parameters
        ----------
        args: tuple
            tuple which contains elevation data:
            if len(args) == 3: args = (lats,lons,elevation) which contains the regular grid values and elevation values on that grid.
            if len(args) == 2: args = (points,elevation) which contains unstructured collection of points for the elevation data.
        smooth: bool, optional
            If True, then use 3rd order bivariate spline to interpolate data, otherwise use faster linear interpolation. default is False.

        """
        if len(args) == 3:
            if smooth:
                self._terrain_list.append(interp.RectBivariateSpline(*args))
            else:
                lats,lons,elevation = args
                self._terrain_list.append(interp.RegularGridInterpolator((lats,lons),elevation,bounds_error=False,fill_value=0.0,method="linear"))
        elif len(args) == 2:
            points,elevation = args
            if smooth:
                self._terrain_list.append(interp.SmoothBivariateSpline(points[:,0],points[:,1],elevation))
            else:
                self._terrain_list.append(interp.LinearNDInterpolator(points,elevation,fill_value=0.0))
        else:
            raise ValueError("can't interpret arguments.")

    def __iadd__(self,other):
        if isinstance(other,land_model):
            self._terrain_list += other._terrain_list
            return self
        else:
            return NotImplemented

    def __add__(self,other):
        if isinstance(other,land_model):
            new = land_model()
            new += other
            return new
        else:
            return NotImplemented

reason: import numpy,import scipy,import numba,from numba


text:
"""Definitions of monomial orderings. """

from __future__ import print_function, division

from typing import Optional

__all__ = ["lex", "grlex", "grevlex", "ilex", "igrlex", "igrevlex"]

from sympy.core import Symbol
from sympy.core.compatibility import iterable

class MonomialOrder(object):
    """Base class for monomial orderings. """

    alias = None  # type: Optional[str]
    is_global = None  # type: Optional[bool]
    is_default = False

    def __repr__(self):
        return self.__class__.__name__ + "()"

    def __str__(self):
        return self.alias

    def __call__(self, monomial):
        raise NotImplementedError

    def __eq__(self, other):
        return self.__class__ == other.__class__

    def __hash__(self):
        return hash(self.__class__)

    def __ne__(self, other):
        return not (self == other)

class LexOrder(MonomialOrder):
    """Lexicographic order of monomials. """

    alias = "lex"
    is_global = True
    is_default = True

    def __call__(self, monomial):
        return monomial

class GradedLexOrder(MonomialOrder):
    """Graded lexicographic order of monomials. """

    alias = "grlex"
    is_global = True

    def __call__(self, monomial):
        return (sum(monomial), monomial)

class ReversedGradedLexOrder(MonomialOrder):
    """Reversed graded lexicographic order of monomials. """

    alias = "grevlex"
    is_global = True

    def __call__(self, monomial):
        return (sum(monomial), tuple(reversed([-m for m in monomial])))

class ProductOrder(MonomialOrder):
    """
    A product order built from other monomial orders.

    Given (not necessarily total) orders O1, O2, ..., On, their product order
    P is defined as M1 > M2 iff there exists i such that O1(M1) = O2(M2),
    ..., Oi(M1) = Oi(M2), O{i+1}(M1) > O{i+1}(M2).

    Product orders are typically built from monomial orders on different sets
    of variables.

    ProductOrder is constructed by passing a list of pairs
    [(O1, L1), (O2, L2), ...] where Oi are MonomialOrders and Li are callables.
    Upon comparison, the Li are passed the total monomial, and should filter
    out the part of the monomial to pass to Oi.

    Examples
    ========

    We can use a lexicographic order on x_1, x_2 and also on
    y_1, y_2, y_3, and their product on {x_i, y_i} as follows:

    >>> from sympy.polys.orderings import lex, grlex, ProductOrder
    >>> P = ProductOrder(
    ...     (lex, lambda m: m[:2]), # lex order on x_1 and x_2 of monomial
    ...     (grlex, lambda m: m[2:]) # grlex on y_1, y_2, y_3
    ... )
    >>> P((2, 1, 1, 0, 0)) > P((1, 10, 0, 2, 0))
    True

    Here the exponent `2` of `x_1` in the first monomial
    (`x_1^2 x_2 y_1`) is bigger than the exponent `1` of `x_1` in the
    second monomial (`x_1 x_2^10 y_2^2`), so the first monomial is greater
    in the product ordering.

    >>> P((2, 1, 1, 0, 0)) < P((2, 1, 0, 2, 0))
    True

    Here the exponents of `x_1` and `x_2` agree, so the grlex order on
    `y_1, y_2, y_3` is used to decide the ordering. In this case the monomial
    `y_2^2` is ordered larger than `y_1`, since for the grlex order the degree
    of the monomial is most important.
    """

    def __init__(self, *args):
        self.args = args

    def __call__(self, monomial):
        return tuple(O(lamda(monomial)) for (O, lamda) in self.args)

    def __repr__(self):
        from sympy.core import Tuple

        return self.__class__.__name__ + repr(Tuple(*[x[0] for x in self.args]))

    def __str__(self):
        from sympy.core import Tuple

        return self.__class__.__name__ + str(Tuple(*[x[0] for x in self.args]))

    def __eq__(self, other):
        if not isinstance(other, ProductOrder):
            return False
        return self.args == other.args

    def __hash__(self):
        return hash((self.__class__, self.args))

    @property
    def is_global(self):
        if all(o.is_global is True for o, _ in self.args):
            return True
        if all(o.is_global is False for o, _ in self.args):
            return False
        return None

class InverseOrder(MonomialOrder):
    """
    The "inverse" of another monomial order.

    If O is any monomial order, we can construct another monomial order iO
    such that `A >_{iO} B` if and only if `B >_O A`. This is useful for
    constructing local orders.

    Note that many algorithms only work with *global* orders.

    For example, in the inverse lexicographic order on a single variable `x`,
    high powers of `x` count as small:

    >>> from sympy.polys.orderings import lex, InverseOrder
    >>> ilex = InverseOrder(lex)
    >>> ilex((5,)) < ilex((0,))
    True
    """

    def __init__(self, O):
        self.O = O

    def __str__(self):
        return "i" + str(self.O)

    def __call__(self, monomial):
        def inv(l):
            if iterable(l):
                return tuple(inv(x) for x in l)
            return -l

        return inv(self.O(monomial))

    @property
    def is_global(self):
        if self.O.is_global is True:
            return False
        if self.O.is_global is False:
            return True
        return None

    def __eq__(self, other):
        return isinstance(other, InverseOrder) and other.O == self.O

    def __hash__(self):
        return hash((self.__class__, self.O))

lex = LexOrder()
grlex = GradedLexOrder()
grevlex = ReversedGradedLexOrder()
ilex = InverseOrder(lex)
igrlex = InverseOrder(grlex)
igrevlex = InverseOrder(grevlex)

_monomial_key = {
    "lex": lex,
    "grlex": grlex,
    "grevlex": grevlex,
    "ilex": ilex,
    "igrlex": igrlex,
    "igrevlex": igrevlex,
}

def monomial_key(order=None, gens=None):
    """
    Return a function defining admissible order on monomials.

    The result of a call to :func:`monomial_key` is a function which should
    be used as a key to :func:`sorted` built-in function, to provide order
    in a set of monomials of the same length.

    Currently supported monomial orderings are:

    1. lex       - lexicographic order (default)
    2. grlex     - graded lexicographic order
    3. grevlex   - reversed graded lexicographic order
    4. ilex, igrlex, igrevlex - the corresponding inverse orders

    If the ``order`` input argument is not a string but has ``__call__``
    attribute, then it will pass through with an assumption that the
    callable object defines an admissible order on monomials.

    If the ``gens`` input argument contains a list of generators, the
    resulting key function can be used to sort SymPy ``Expr`` objects.

    """
    if order is None:
        order = lex

    if isinstance(order, Symbol):
        order = str(order)

    if isinstance(order, str):
        try:
            order = _monomial_key[order]
        except KeyError:
            raise ValueError(
                "supported monomial orderings are 'lex', 'grlex' and 'grevlex', got %r"
                % order
            )
    if hasattr(order, "__call__"):
        if gens is not None:

            def _order(expr):
                return order(expr.as_poly(*gens).degree_list())

            return _order
        return order
    else:
        raise ValueError(
            "monomial ordering specification must be a string or a callable, got %s"
            % order
        )

class _ItemGetter(object):
    """Helper class to return a subsequence of values."""

    def __init__(self, seq):
        self.seq = tuple(seq)

    def __call__(self, m):
        return tuple(m[idx] for idx in self.seq)

    def __eq__(self, other):
        if not isinstance(other, _ItemGetter):
            return False
        return self.seq == other.seq

def build_product_order(arg, gens):
    """
    Build a monomial order on ``gens``.

    ``arg`` should be a tuple of iterables. The first element of each iterable
    should be a string or monomial order (will be passed to monomial_key),
    the others should be subsets of the generators. This function will build
    the corresponding product order.

    For example, build a product of two grlex orders:

    >>> from sympy.polys.orderings import grlex, build_product_order
    >>> from sympy.abc import x, y, z, t

    >>> O = build_product_order((("grlex", x, y), ("grlex", z, t)), [x, y, z, t])
    >>> O((1, 2, 3, 4))
    ((3, (1, 2)), (7, (3, 4)))

    """
    gens2idx = {}
    for i, g in enumerate(gens):
        gens2idx[g] = i
    order = []
    for expr in arg:
        name = expr[0]
        var = expr[1:]

        def makelambda(var):
            return _ItemGetter(gens2idx[g] for g in var)

        order.append((monomial_key(name), makelambda(var)))
    return ProductOrder(*order)

ext: py

max_stars_repo_name: msgoff/sympy

max_issues_repo_name: msgoff/sympy

max_forks_repo_name: msgoff/sympy

reason: from sympy


text:
from typing import Dict, Tuple, Any, List
import logging
import copy

from overrides import overrides
import torch
from torch.nn.modules import Dropout
import numpy

from allennlp.common.checks import check_dimensions_match, ConfigurationError
from allennlp.data import TextFieldTensors, Vocabulary
from allennlp.modules import Seq2SeqEncoder, TextFieldEmbedder, Embedding, InputVariationalDropout
from allennlp.modules.matrix_attention.bilinear_matrix_attention import BilinearMatrixAttention
from allennlp.modules import FeedForward
from allennlp.models.model import Model
from allennlp.nn import InitializerApplicator, Activation
from allennlp.nn.util import get_text_field_mask
from allennlp.nn.util import get_lengths_from_binary_sequence_mask
from allennlp.training.metrics import F1Measure

logger = logging.getLogger(__name__)

@Model.register("graph_parser")
class GraphParser(Model):
    """
    A Parser for arbitrary graph structures.

    # Parameters

    vocab : `Vocabulary`, required
        A Vocabulary, required in order to compute sizes for input/output projections.
    text_field_embedder : `TextFieldEmbedder`, required
        Used to embed the `tokens` `TextField` we get as input to the model.
    encoder : `Seq2SeqEncoder`
        The encoder (with its own internal stacking) that we will use to generate representations
        of tokens.
    tag_representation_dim : `int`, required.
        The dimension of the MLPs used for arc tag prediction.
    arc_representation_dim : `int`, required.
        The dimension of the MLPs used for arc prediction.
    tag_feedforward : `FeedForward`, optional, (default = None).
        The feedforward network used to produce tag representations.
        By default, a 1 layer feedforward network with an elu activation is used.
    arc_feedforward : `FeedForward`, optional, (default = None).
        The feedforward network used to produce arc representations.
        By default, a 1 layer feedforward network with an elu activation is used.
    pos_tag_embedding : `Embedding`, optional.
        Used to embed the `pos_tags` `SequenceLabelField` we get as input to the model.
    dropout : `float`, optional, (default = 0.0)
        The variational dropout applied to the output of the encoder and MLP layers.
    input_dropout : `float`, optional, (default = 0.0)
        The dropout applied to the embedded text input.
    edge_prediction_threshold : `int`, optional (default = 0.5)
        The probability at which to consider a scored edge to be 'present'
        in the decoded graph. Must be between 0 and 1.
    initializer : `InitializerApplicator`, optional (default=`InitializerApplicator()`)
        Used to initialize the model parameters.
    """

    def __init__(
        self,
        vocab: Vocabulary,
        text_field_embedder: TextFieldEmbedder,
        encoder: Seq2SeqEncoder,
        tag_representation_dim: int,
        arc_representation_dim: int,
        tag_feedforward: FeedForward = None,
        arc_feedforward: FeedForward = None,
        pos_tag_embedding: Embedding = None,
        dropout: float = 0.0,
        input_dropout: float = 0.0,
        edge_prediction_threshold: float = 0.5,
        initializer: InitializerApplicator = InitializerApplicator(),
        **kwargs,
    ) -> None:
        super().__init__(vocab, **kwargs)

        self.text_field_embedder = text_field_embedder
        self.encoder = encoder
        self.edge_prediction_threshold = edge_prediction_threshold
        if not 0 < edge_prediction_threshold < 1:
            raise ConfigurationError(
                f"edge_prediction_threshold must be between "
                f"0 and 1 (exclusive) but found {edge_prediction_threshold}."
            )

        encoder_dim = encoder.get_output_dim()

        self.head_arc_feedforward = arc_feedforward or FeedForward(
            encoder_dim, 1, arc_representation_dim, Activation.by_name("elu")()
        )
        self.child_arc_feedforward = copy.deepcopy(self.head_arc_feedforward)

        self.arc_attention = BilinearMatrixAttention(
            arc_representation_dim, arc_representation_dim, use_input_biases=True
        )

        num_labels = self.vocab.get_vocab_size("labels")
        self.head_tag_feedforward = tag_feedforward or FeedForward(
            encoder_dim, 1, tag_representation_dim, Activation.by_name("elu")()
        )
        self.child_tag_feedforward = copy.deepcopy(self.head_tag_feedforward)

        self.tag_bilinear = BilinearMatrixAttention(
            tag_representation_dim, tag_representation_dim, label_dim=num_labels
        )

        self._pos_tag_embedding = pos_tag_embedding or None
        self._dropout = InputVariationalDropout(dropout)
        self._input_dropout = Dropout(input_dropout)

        representation_dim = text_field_embedder.get_output_dim()
        if pos_tag_embedding is not None:
            representation_dim += pos_tag_embedding.get_output_dim()

        check_dimensions_match(
            representation_dim,
            encoder.get_input_dim(),
            "text field embedding dim",
            "encoder input dim",
        )
        check_dimensions_match(
            tag_representation_dim,
            self.head_tag_feedforward.get_output_dim(),
            "tag representation dim",
            "tag feedforward output dim",
        )
        check_dimensions_match(
            arc_representation_dim,
            self.head_arc_feedforward.get_output_dim(),
            "arc representation dim",
            "arc feedforward output dim",
        )

        self._unlabelled_f1 = F1Measure(positive_label=1)
        self._arc_loss = torch.nn.BCEWithLogitsLoss(reduction="none")
        self._tag_loss = torch.nn.CrossEntropyLoss(reduction="none")
        initializer(self)

    @overrides
    def forward(
        self,  # type: ignore
        tokens: TextFieldTensors,
        pos_tags: torch.LongTensor = None,
        metadata: List[Dict[str, Any]] = None,
        arc_tags: torch.LongTensor = None,
    ) -> Dict[str, torch.Tensor]:

        """
        # Parameters

        tokens : TextFieldTensors, required
            The output of `TextField.as_array()`.
        pos_tags : torch.LongTensor, optional (default = None)
            The output of a `SequenceLabelField` containing POS tags.
        metadata : List[Dict[str, Any]], optional (default = None)
            A dictionary of metadata for each batch element which has keys:
                tokens : `List[str]`, required.
                    The original string tokens in the sentence.
        arc_tags : torch.LongTensor, optional (default = None)
            A torch tensor representing the sequence of integer indices denoting the parent of every
            word in the dependency parse. Has shape `(batch_size, sequence_length, sequence_length)`.

        # Returns

        An output dictionary.
        """
        embedded_text_input = self.text_field_embedder(tokens)
        if pos_tags is not None and self._pos_tag_embedding is not None:
            embedded_pos_tags = self._pos_tag_embedding(pos_tags)
            embedded_text_input = torch.cat([embedded_text_input, embedded_pos_tags], -1)
        elif self._pos_tag_embedding is not None:
            raise ConfigurationError("Model uses a POS embedding, but no POS tags were passed.")

        mask = get_text_field_mask(tokens)
        embedded_text_input = self._input_dropout(embedded_text_input)
        encoded_text = self.encoder(embedded_text_input, mask)

        float_mask = mask.float()
        encoded_text = self._dropout(encoded_text)

        # shape (batch_size, sequence_length, arc_representation_dim)
        head_arc_representation = self._dropout(self.head_arc_feedforward(encoded_text))
        child_arc_representation = self._dropout(self.child_arc_feedforward(encoded_text))

        # shape (batch_size, sequence_length, tag_representation_dim)
        head_tag_representation = self._dropout(self.head_tag_feedforward(encoded_text))
        child_tag_representation = self._dropout(self.child_tag_feedforward(encoded_text))
        # shape (batch_size, sequence_length, sequence_length)
        arc_scores = self.arc_attention(head_arc_representation, child_arc_representation)
        # shape (batch_size, num_tags, sequence_length, sequence_length)
        arc_tag_logits = self.tag_bilinear(head_tag_representation, child_tag_representation)
        # Switch to (batch_size, sequence_length, sequence_length, num_tags)
        arc_tag_logits = arc_tag_logits.permute(0, 2, 3, 1).contiguous()

        minus_inf = -1e8
        minus_mask = (1 - float_mask) * minus_inf
        arc_scores = arc_scores + minus_mask.unsqueeze(2) + minus_mask.unsqueeze(1)

        arc_probs, arc_tag_probs = self._greedy_decode(arc_scores, arc_tag_logits, mask)

        output_dict = {"arc_probs": arc_probs, "arc_tag_probs": arc_tag_probs, "mask": mask}

        if metadata:
            output_dict["tokens"] = [meta["tokens"] for meta in metadata]

        if arc_tags is not None:
            arc_nll, tag_nll = self._construct_loss(
                arc_scores=arc_scores, arc_tag_logits=arc_tag_logits, arc_tags=arc_tags, mask=mask
            )
            output_dict["loss"] = arc_nll + tag_nll
            output_dict["arc_loss"] = arc_nll
            output_dict["tag_loss"] = tag_nll

            # Make the arc tags not have negative values anywhere
            # (by default, no edge is indicated with -1).
            arc_indices = (arc_tags != -1).float()
            tag_mask = float_mask.unsqueeze(1) * float_mask.unsqueeze(2)
            one_minus_arc_probs = 1 - arc_probs
            # We stack scores here because the f1 measure expects a
            # distribution, rather than a single value.
            self._unlabelled_f1(
                torch.stack([one_minus_arc_probs, arc_probs], -1), arc_indices, tag_mask
            )

        return output_dict

    @overrides
    def make_output_human_readable(
        self, output_dict: Dict[str, torch.Tensor]
    ) -> Dict[str, torch.Tensor]:
        arc_tag_probs = output_dict["arc_tag_probs"].cpu().detach().numpy()
        arc_probs = output_dict["arc_probs"].cpu().detach().numpy()
        mask = output_dict["mask"]
        lengths = get_lengths_from_binary_sequence_mask(mask)
        arcs = []
        arc_tags = []
        for instance_arc_probs, instance_arc_tag_probs, length in zip(
            arc_probs, arc_tag_probs, lengths
        ):

            arc_matrix = instance_arc_probs > self.edge_prediction_threshold
            edges = []
            edge_tags = []
            for i in range(length):
                for j in range(length):
                    if arc_matrix[i, j] == 1:
                        edges.append((i, j))
                        tag = instance_arc_tag_probs[i, j].argmax(-1)
                        edge_tags.append(self.vocab.get_token_from_index(tag, "labels"))
            arcs.append(edges)
            arc_tags.append(edge_tags)

        output_dict["arcs"] = arcs
        output_dict["arc_tags"] = arc_tags
        return output_dict

    def _construct_loss(
        self,
        arc_scores: torch.Tensor,
        arc_tag_logits: torch.Tensor,
        arc_tags: torch.Tensor,
        mask: torch.Tensor,
    ) -> Tuple[torch.Tensor, torch.Tensor]:
        """
        Computes the arc and tag loss for an adjacency matrix.

        # Parameters

        arc_scores : `torch.Tensor`, required.
            A tensor of shape (batch_size, sequence_length, sequence_length) used to generate a
            binary classification decision for whether an edge is present between two words.
        arc_tag_logits : `torch.Tensor`, required.
            A tensor of shape (batch_size, sequence_length, sequence_length, num_tags) used to generate
            a distribution over edge tags for a given edge.
        arc_tags : `torch.Tensor`, required.
            A tensor of shape (batch_size, sequence_length, sequence_length).
            The labels for every arc.
        mask : `torch.Tensor`, required.
            A mask of shape (batch_size, sequence_length), denoting unpadded
            elements in the sequence.

        # Returns

        arc_nll : `torch.Tensor`, required.
            The negative log likelihood from the arc loss.
        tag_nll : `torch.Tensor`, required.
            The negative log likelihood from the arc tag loss.
        """
        float_mask = mask.float()
        arc_indices = (arc_tags != -1).float()
        # Make the arc tags not have negative values anywhere
        # (by default, no edge is indicated with -1).
        arc_tags = arc_tags * arc_indices
        arc_nll = (
            self._arc_loss(arc_scores, arc_indices)
            * float_mask.unsqueeze(1)
            * float_mask.unsqueeze(2)
        )
        # We want the mask for the tags to only include the unmasked words
        # and we only care about the loss with respect to the gold arcs.
        tag_mask = float_mask.unsqueeze(1) * float_mask.unsqueeze(2) * arc_indices

        batch_size, sequence_length, _, num_tags = arc_tag_logits.size()
        original_shape = [batch_size, sequence_length, sequence_length]
        reshaped_logits = arc_tag_logits.view(-1, num_tags)
        reshaped_tags = arc_tags.view(-1)
        tag_nll = (
            self._tag_loss(reshaped_logits, reshaped_tags.long()).view(original_shape) * tag_mask
        )

        valid_positions = tag_mask.sum()

        arc_nll = arc_nll.sum() / valid_positions.float()
        tag_nll = tag_nll.sum() / valid_positions.float()
        return arc_nll, tag_nll

    @staticmethod
    def _greedy_decode(
        arc_scores: torch.Tensor, arc_tag_logits: torch.Tensor, mask: torch.Tensor
    ) -> Tuple[torch.Tensor, torch.Tensor]:
        """
        Decodes the head and head tag predictions by decoding the unlabeled arcs
        independently for each word and then again, predicting the head tags of
        these greedily chosen arcs independently.

        # Parameters

        arc_scores : `torch.Tensor`, required.
            A tensor of shape (batch_size, sequence_length, sequence_length) used to generate
            a distribution over attachments of a given word to all other words.
        arc_tag_logits : `torch.Tensor`, required.
            A tensor of shape (batch_size, sequence_length, sequence_length, num_tags) used to
            generate a distribution over tags for each arc.
        mask : `torch.Tensor`, required.
            A mask of shape (batch_size, sequence_length).

        # Returns

        arc_probs : `torch.Tensor`
            A tensor of shape (batch_size, sequence_length, sequence_length) representing the
            probability of an arc being present for this edge.
        arc_tag_probs : `torch.Tensor`
            A tensor of shape (batch_size, sequence_length, sequence_length, sequence_length)
            representing the distribution over edge tags for a given edge.
        """
        # Mask the diagonal, because we don't self edges.
        inf_diagonal_mask = torch.diag(arc_scores.new(mask.size(1)).fill_(-numpy.inf))
        arc_scores = arc_scores + inf_diagonal_mask
        # shape (batch_size, sequence_length, sequence_length, num_tags)
        arc_tag_logits = arc_tag_logits + inf_diagonal_mask.unsqueeze(0).unsqueeze(-1)
        # Mask padded tokens, because we only want to consider actual word -> word edges.
        minus_mask = (1 - mask).to(dtype=torch.bool).unsqueeze(2)
        arc_scores.masked_fill_(minus_mask, -numpy.inf)
        arc_tag_logits.masked_fill_(minus_mask.unsqueeze(-1), -numpy.inf)
        # shape (batch_size, sequence_length, sequence_length)
        arc_probs = arc_scores.sigmoid()
        # shape (batch_size, sequence_length, sequence_length, num_tags)
        arc_tag_probs = torch.nn.functional.softmax(arc_tag_logits, dim=-1)
        return arc_probs, arc_tag_probs

    @overrides
    def get_metrics(self, reset: bool = False) -> Dict[str, float]:
        metrics = {}
        precision, recall, f1_measure = self._unlabelled_f1.get_metric(reset)
        metrics["precision"] = precision
        metrics["recall"] = recall
        metrics["f1"] = f1_measure
        return metrics

reason: import numpy


text:
import tensorflow as tf
import numpy as np
from tensorflow.contrib.rnn import LSTMCell, LSTMStateTuple
from tensorflow.python.ops.rnn import bidirectional_dynamic_rnn as bi_rnn
import networkx as nx

class STNE(object):
    ############
    def encoder(self,X,layers):
        for i in range(layers - 1):
            name_W = 'encoder_W_' + str(i)
            name_b = 'encoder_b_' + str(i)
            X = tf.nn.tanh(tf.matmul(X, self.W[name_W]) + self.b[name_b])
        return X

    def decoder(self,X,layers):
        for i in range(layers - 1):
            name_W = 'decoder_W_' + str(i)
            name_b = 'decoder_b_' + str(i)
            X = tf.nn.tanh(tf.matmul(X, self.W[name_W]) + self.b[name_b])
        return X

    def make_autoencoder_loss(self,X_new,X_re):
        def get_autoencoder_loss(X, newX):
            return tf.reduce_sum(tf.pow((newX - X), 2))

        def get_reg_loss(weights, biases):
            reg = tf.add_n([tf.nn.l2_loss(w) for w in weights.values()])
            reg += tf.add_n([tf.nn.l2_loss(b) for b in biases.values()])
            return reg
        loss_autoencoder = get_autoencoder_loss(X_new, X_re)
        loss_reg = get_reg_loss(self.W, self.b)
        return self.config.alpha * loss_autoencoder + self.config.reg * loss_reg
    ################
    def construct_traget_neighbors(self, nx_G, X, mode='EMN'):
        # construct target neighbor feature matrix
        X_target = np.zeros(X.shape)
        nodes = nx_G.nodes()
        if mode == 'OWN':
            # autoencoder for reconstructing itself
            return X
        elif mode == 'EMN':
            # autoencoder for reconstructing Elementwise Median Neighbor
            for node in nodes:
                neighbors = list(nx_G.neighbors(node))
                if len(neighbors) == 0:
                    X_target[node] = X[node]
                else:
                    temp = np.array(X[node])
                    for n in neighbors:
                        # # if FLAGS.weighted:
                        # #     # weighted sum
                        # #     # temp = np.vstack((temp, X[n] * edgeWeight))
                        # #     pass
                        # else:
                        temp = np.vstack((temp, X[n]))
                    temp = np.median(temp, axis=0)
                    X_target[node] = temp
            return X_target
    #################
    def __init__(self,config, hidden_dim, nx_G, X_1,node_num, fea_dim, seq_len, attention_size,
                 depth=1, node_fea=None, node_fea_trainable=False):
        self.node_num, self.fea_dim, self.seq_len = node_num, fea_dim, seq_len
        self.attention_size = attention_size
        self.nx_G = nx_G
        self.X_1 = X_1
        self.config = config
        self.keep_prob = tf.placeholder(tf.float32, name='keep_prob')
        self.input_seqs = tf.placeholder(tf.int32, shape=(None, self.seq_len), name='input_seq')
        print(self.input_seqs,'AZAZAZAAAAAAAAAAAAAAA')
        self.dropout = tf.placeholder(tf.float32, name='dropout')
        X_target = self.construct_traget_neighbors(self.nx_G, self.X_1, mode='EMN')
        print(X_target, '222222222222222222222222')
        X_target = tf.constant(X_target, dtype=tf.float32)
        self.layers = len(config.struct)
        struct = config.struct
        self.W = {}
        self.b = {}
        # encode module
        for i in range(self.layers - 1):
            name_W = 'encoder_W_' + str(i)
            name_b = 'encoder_b_' + str(i)
            print(struct[i], '55555555555555555555')
            print(struct[i + 1], '6666666666666666666')
            self.W[name_W] = tf.get_variable(name_W, [struct[i], struct[i + 1]],
                                             initializer=tf.contrib.layers.xavier_initializer())
            print(self.W[name_W], 'AAAAAAAAAAAAAAAAAAAAAAA')
            self.b[name_b] = tf.get_variable(name_b, [struct[i + 1]], initializer=tf.zeros_initializer())
            print(self.b[name_b], 'SSSSSSSSSSSSSSSSSSSSS')
        # decode module
        struct.reverse()
        for i in range(self.layers - 1):
            name_W = 'decoder_W_' + str(i)
            name_b = 'decoder_b_' + str(i)
            self.W[name_W] = tf.get_variable(name_W, [struct[i], struct[i + 1]],
                                             initializer=tf.contrib.layers.xavier_initializer())
            self.b[name_b] = tf.get_variable(name_b, [struct[i + 1]], initializer=tf.zeros_initializer())
        config.struct.reverse()
        ############## define input ###################
        self.Y1 = self.encoder(self.X_1,self.layers)
        self.X1_reconstruct = self.decoder(self.Y1, self.layers)
        self.loss_autoencoder_1 = self.make_autoencoder_loss(X_target, self.X1_reconstruct)
        input_seq_embed = tf.nn.embedding_lookup(self.Y1, self.input_seqs, name='input_embed_lookup')
        print(input_seq_embed, '44444444444444444NNNNNNNNNNNNNNNNNN4')

        # encoder
        encoder_cell_fw_0 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim), output_keep_prob=1 - self.dropout)
        encoder_cell_bw_0 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim), output_keep_prob=1 - self.dropout)
        if depth == 1:
            encoder_cell_fw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_fw_0])
            encoder_cell_bw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_bw_0])
        else:
            encoder_cell_fw_1 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim), output_keep_prob=1 - self.dropout)
            encoder_cell_bw_1 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim), output_keep_prob=1 - self.dropout)

            encoder_cell_fw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_fw_0] + [encoder_cell_fw_1] * (depth - 1))
            encoder_cell_bw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_bw_0] + [encoder_cell_bw_1] * (depth - 1))

        encoder_outputs, encoder_final = bi_rnn(encoder_cell_fw_all, encoder_cell_bw_all, inputs=input_seq_embed,
                                                dtype=tf.float32)
        ######################Attention
        self.encoder_output = tf.concat(encoder_outputs, 2)
        input_shape = self.encoder_output.shape  # (batch_size, sequence_length, hidden_size)
        sequence_size = input_shape[1].value
        hidden_size = input_shape[2].value
        attention_w = tf.Variable(tf.truncated_normal([hidden_size, self.attention_size], stddev=0.1),
                                  name='attention_w')
        attention_b = tf.Variable(tf.constant(0.1, shape=[self.attention_size]), name='attention_b')
        attention_u = tf.Variable(tf.truncated_normal([self.attention_size], stddev=0.1), name='attention_u')
        z_list = []
        for t in range(sequence_size):
            u_t = tf.tanh(tf.matmul(self.encoder_output[:, t, :], attention_w) + tf.reshape(attention_b, [1, -1]))
            z_t = tf.matmul(u_t, tf.reshape(attention_u, [-1, 1]))
            z_list.append(z_t)
        # Transform to batch_size * sequence_size

        attention_z = tf.concat(z_list, axis=1)
        alpha = tf.nn.softmax(attention_z)

        # print((self.encoder_output * tf.reshape(self.alpha, [-1, sequence_size, 1])), '++++++++++++++++++++++')
        attention_output = self.encoder_output * tf.reshape(alpha, [-1, sequence_size, 1])

        final_output = tf.nn.dropout(attention_output, self.keep_prob)
        ######################
        c_fw_list, h_fw_list, c_bw_list, h_bw_list = [], [], [], []
        for d in range(depth):
            (c_fw, h_fw) = encoder_final[0][d]
            (c_bw, h_bw) = encoder_final[1][d]
            c_fw_list.append(c_fw)
            h_fw_list.append(h_fw)
            c_bw_list.append(c_bw)
            h_bw_list.append(h_bw)

        decoder_init_state = tf.concat(c_fw_list + c_bw_list, axis=-1), tf.concat(h_fw_list + h_bw_list, axis=-1)
        decoder_cell = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim * 2), output_keep_prob=1 - self.dropout)
        decoder_init_state = LSTMStateTuple(
            tf.layers.dense(decoder_init_state[0], units=hidden_dim * 2, activation=None),
            tf.layers.dense(decoder_init_state[1], units=hidden_dim * 2, activation=None))
        encoder_output_T = tf.transpose(final_output, [1, 0, 2])
        new_state = decoder_init_state
        outputs_list = []
        for i in range(seq_len):
            new_output, new_state = decoder_cell(tf.zeros(shape=tf.shape(encoder_output_T)[1:]), new_state)  # None
            outputs_list.append(new_output)

        decoder_outputs = tf.stack(outputs_list, axis=0)  # seq_len * batch_size * hidden_dim
        decoder_outputs = tf.transpose(decoder_outputs, [1, 0, 2])  # batch_size * seq_len * hidden_dim
        self.decoder_outputs = decoder_outputs
        output_preds = tf.layers.dense(decoder_outputs, units=self.node_num, activation=None)
        loss_ce = tf.nn.sparse_softmax_cross_entropy_with_logits(labels=self.input_seqs, logits=output_preds)
        self.loss_ce = tf.reduce_mean(loss_ce, name='loss_ce')
        self.toll_loss=config.alpha*self.loss_ce+config.reg*self.loss_autoencoder_1
        self.train_op = tf.train.RMSPropOptimizer(config.sg_learning_rate).minimize(self.toll_loss)

class STNEConv(object):
    def conv_pool(self, in_tensor, filter_size, num_filters, s_length, embedding_size=256):
        with tf.name_scope("conv-maxpool-%s" % filter_size):
            # Convolution Layer
            filter_shape = [filter_size, embedding_size, 1, num_filters]
            W = tf.Variable(tf.truncated_normal(filter_shape, stddev=0.1), name="W")
            b = tf.Variable(tf.constant(0.1, shape=[num_filters]), name="b")
            # None, seq_len, word_dim, 1
            conv = tf.nn.conv2d(in_tensor, W, strides=[1, 1, 1, 1], padding="VALID", name="conv")
            h = tf.nn.relu(tf.nn.bias_add(conv, b), name="relu")
            # Maxpooling over the outputs
            pooled = tf.nn.max_pool(
                h,
                ksize=[1, s_length - filter_size + 1, 1, 1],
                strides=[1, 1, 1, 1],
                padding='VALID',
                name="pool")
            # print(pooled)
        return pooled  # None, 1, 1, num_filters

    def __init__(self, hidden_dim, node_num, fea_dim, seq_len, contnt_len, num_filters, word_dim,
                 vocab_size, attention_size, depth=1, filter_sizes=[2, 4, 8]):
        self.node_num, self.fea_dim = node_num, fea_dim
        self.attention_size = attention_size
        self.seq_len = seq_len
        self.keep_prob = tf.placeholder(tf.float32, name='keep_prob')
        self.input_seqs = tf.placeholder(tf.int32, shape=(None, self.seq_len), name='input_seq')
        self.input_seq_content = tf.placeholder(tf.int32, shape=(None, self.seq_len, contnt_len),
                                                name='input_seq_content')
        self.dropout_rnn = tf.placeholder(tf.float32, name='dropout_rnn')
        self.dropout_word = tf.placeholder(tf.float32, name='dropout_word')
        self.word_embeds_W = tf.Variable(initial_value=tf.random_uniform(shape=(vocab_size, word_dim)),
                                         name='content_embed', trainable=True)

        contnt_embeds = tf.nn.embedding_lookup(self.word_embeds_W, self.input_seq_content, name='input_content_embed')
        contnt_embeds = tf.reshape(contnt_embeds, [-1, contnt_len, word_dim, 1])
        print(contnt_embeds, '1111111111111111111111')
        pooled = []
        for fsize in filter_sizes:
            # batch*seq_len, 1, num_filters
            tmp = self.conv_pool(contnt_embeds, fsize, num_filters, contnt_len, word_dim)
            pooled.append(tf.reshape(tmp, [-1, self.seq_len, num_filters]))
        input_seq_embed = tf.concat(pooled, axis=-1)  # batch, seq_len, num_filters*len(filter_sizes)
        input_seq_embed = tf.nn.dropout(input_seq_embed, keep_prob=1 - self.dropout_word)
        print(input_seq_embed, '22222222222222222222222')
        encoder_cell_fw_0 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim), output_keep_prob=1 - self.dropout_rnn)
        encoder_cell_bw_0 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim), output_keep_prob=1 - self.dropout_rnn)
        if depth == 1:
            encoder_cell_fw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_fw_0])
            encoder_cell_bw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_bw_0])
        else:
            encoder_cell_fw_1 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim),
                                                              output_keep_prob=1 - self.dropout_rnn)
            encoder_cell_bw_1 = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim),
                                                              output_keep_prob=1 - self.dropout_rnn)

            encoder_cell_fw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_fw_0] + [encoder_cell_fw_1] * (depth - 1))
            encoder_cell_bw_all = tf.contrib.rnn.MultiRNNCell([encoder_cell_bw_0] + [encoder_cell_bw_1] * (depth - 1))

        encoder_outputs, encoder_final = bi_rnn(encoder_cell_fw_all, encoder_cell_bw_all, inputs=input_seq_embed,
                                                dtype=tf.float32)
        ############################
        self.encoder_output = tf.concat(encoder_outputs, 2)
        input_shape = self.encoder_output.shape  # (batch_size, sequence_length, hidden_size)
        sequence_size = input_shape[1].value
        hidden_size = input_shape[2].value
        attention_w = tf.Variable(tf.truncated_normal([hidden_size, self.attention_size], stddev=0.1),
                                  name='attention_w')
        attention_b = tf.Variable(tf.constant(0.1, shape=[self.attention_size]), name='attention_b')
        attention_u = tf.Variable(tf.truncated_normal([self.attention_size], stddev=0.1), name='attention_u')
        z_list = []
        for t in range(sequence_size):
            u_t = tf.tanh(tf.matmul(self.encoder_output[:, t, :], attention_w) + tf.reshape(attention_b, [1, -1]))
            z_t = tf.matmul(u_t, tf.reshape(attention_u, [-1, 1]))
            z_list.append(z_t)
        # Transform to batch_size * sequence_size

        attention_z = tf.concat(z_list, axis=1)
        self.alpha = tf.nn.softmax(attention_z)

        # print((self.encoder_output * tf.reshape(self.alpha, [-1, sequence_size, 1])), '++++++++++++++++++++++')
        attention_output = self.encoder_output * tf.reshape(self.alpha, [-1, sequence_size, 1])

        self.final_output = tf.nn.dropout(attention_output, self.keep_prob)
        ####################################
        c_fw_list, h_fw_list, c_bw_list, h_bw_list = [], [], [], []

        for d in range(depth):
            (c_fw, h_fw) = encoder_final[0][d]
            (c_bw, h_bw) = encoder_final[1][d]
            c_fw_list.append(c_fw)
            h_fw_list.append(h_fw)
            c_bw_list.append(c_bw)
            h_bw_list.append(h_bw)

        decoder_init_state = tf.concat(c_fw_list + c_bw_list, axis=-1), tf.concat(h_fw_list + h_bw_list, axis=-1)
        decoder_cell = tf.contrib.rnn.DropoutWrapper(LSTMCell(hidden_dim * 2), output_keep_prob=1 - self.dropout_rnn)
        decoder_init_state = LSTMStateTuple(
            tf.layers.dense(decoder_init_state[0], units=hidden_dim * 2, activation=None),
            tf.layers.dense(decoder_init_state[1], units=hidden_dim * 2, activation=None))

        # self.encoder_output = tf.concat(encoder_outputs, axis=-1)
        # encoder_output_T = tf.transpose(self.encoder_output, [1, 0, 2])  # h
        encoder_output_T = tf.transpose(self.final_output, [1, 0, 2])
        new_state = decoder_init_state
        outputs_list = []
        for i in range(seq_len):
            new_output, new_state = decoder_cell(tf.zeros(shape=tf.shape(encoder_output_T)[1:]), new_state)  # None
            outputs_list.append(new_output)

        decoder_outputs = tf.stack(outputs_list, axis=0)  # seq_len * batch_size * hidden_dim
        decoder_outputs = tf.transpose(decoder_outputs, [1, 0, 2])  # batch_size * seq_len * hidden_dim
        ###########################
        input_shape = decoder_outputs.shape  # (batch_size, sequence_length, hidden_size)
        sequence_size = input_shape[1].value
        hidden_size = input_shape[2].value
        attention_w = tf.Variable(tf.truncated_normal([hidden_size, self.attention_size], stddev=0.1),
                                  name='attention_w')
        attention_b = tf.Variable(tf.constant(0.1, shape=[self.attention_size]), name='attention_b')
        attention_u = tf.Variable(tf.truncated_normal([self.attention_size], stddev=0.1), name='attention_u')
        z_list = []
        for t in range(sequence_size):
            u_t = tf.tanh(tf.matmul(decoder_outputs[:, t, :], attention_w) + tf.reshape(attention_b, [1, -1]))
            z_t = tf.matmul(u_t, tf.reshape(attention_u, [-1, 1]))
            z_list.append(z_t)
        # Transform to batch_size * sequence_size

        attention_z = tf.concat(z_list, axis=1)
        self.alpha = tf.nn.softmax(attention_z)

        # print((self.encoder_output * tf.reshape(self.alpha, [-1, sequence_size, 1])), '++++++++++++++++++++++')
        attention_output = decoder_outputs * tf.reshape(self.alpha, [-1, sequence_size, 1])

        final_output_1 = tf.nn.dropout(attention_output, self.keep_prob)
        #####################
        self.decoder_output = final_output_1
        # decoder_outputs, _ = dynamic_rnn(decoder_cell, inputs=self.encoder_output, initial_state=decoder_init_state)
        output_preds = tf.layers.dense(final_output_1, units=self.node_num, activation=None)
        print(self.input_seqs, '@@@@@@@@@@@@@@@@@@@@')
        print(output_preds, '$$$$$$$$$$$$$$$$$$$$')
        loss_ce = tf.nn.sparse_softmax_cross_entropy_with_logits(labels=self.input_seqs, logits=output_preds)
        self.loss_ce = tf.reduce_mean(loss_ce, name='loss_ce')

        self.global_step = tf.Variable(1, name="global_step", trainable=False)

reason: import numpy,import networkx