text:
# coding: utf-8

# In[1]:

import math
import random

import gym
import numpy as np

import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
from torch.distributions import Normal
import visdom

vis = visdom.Visdom()
vis.close()
# In[2]:

from IPython.display import clear_output
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

#
Use CUDA

# In[3]:

use_cuda = torch.cuda.is_available()
device   = torch.device("cuda" if use_cuda else "cpu")
print(device)

#
Replay Buffer

# In[5]:

class ReplayBuffer:
    def __init__(self, capacity):
        self.capacity = capacity
        self.buffer = []
        self.position = 0

#    def push(self,mem):
#        if len(self.buffer) < self.capacity:
#            self.buffer.append(None)
#        self.buffer[self.position] = mem
#        self.position = (self.position + 1) % self.capacity
#

    def push(self, state, action, reward, next_state, done):
        if len(self.buffer) < self.capacity:
            self.buffer.append(None)
        self.buffer[self.position] = (state, action, reward, next_state, done)
        self.position = (self.position + 1) % self.capacity

    def sample(self, batch_size,seq_size=1):
        batch = random.sample(self.buffer, batch_size)
        state, action, reward, next_state, done = map(np.stack, zip(*batch))
        return state, action, reward, next_state, done
#        batch_data = random.sample(self.buffer,batch_size)
#        state, action, reward, next_state, done = [],[],[],[],[]
#        for data in batch_data:
#            start = random.randint(0,len(data)-seq_size)
#            aa = data[start:start+seq_size]
#            t0,t1,t2,t3,t4 = list(zip(*aa))
#
#            state.append(torch.stack(t0,dim=1))
#            action.append(torch.stack(t1,dim=1))
#            reward.append(torch.stack(t2,dim=0))
#            next_state.append(torch.stack(t3,dim=1))
#            done.append(torch.tensor(t4))
#
#        state  = torch.cat(state).reshape((batch_size,seq_size,-1)).to(device)
#        action = torch.cat(action).reshape((batch_size,seq_size,-1)).to(device)
#        reward = torch.cat(reward).reshape((batch_size,seq_size,-1)).to(device)
#        next_state = torch.cat(next_state).reshape((batch_size,seq_size,-1)).to(device)
#        done = torch.cat(done).reshape((batch_size,seq_size,-1)).to(device)
#        return state, action, reward, next_state, done

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

#
Ornstein-Uhlenbeck process

# Adding time-correlated noise to the actions taken by the deterministic policy

# wiki

# In[12]:

class OUNoise(object):
    def __init__(self, action_space, mu=0.0, theta=0.15, max_sigma=0.3, min_sigma=0.3, decay_period=100000):
        self.mu           = mu
        self.theta        = theta
        self.sigma        = max_sigma
        self.max_sigma    = max_sigma
        self.min_sigma    = min_sigma
        self.decay_period = decay_period
        self.action_dim   = action_space
        self.low          = -1
        self.high         = 1
        self.reset()

    def reset(self):
        self.state = np.ones(self.action_dim) * self.mu

    def evolve_state(self):
        x  = self.state
        dx = self.theta * (self.mu - x) + self.sigma * np.random.randn(self.action_dim)
        self.state = x + dx
        return self.state

    def get_action(self, action, t=0):
        ou_state = torch.from_numpy(self.evolve_state()).float()
        self.sigma = self.max_sigma - (self.max_sigma - self.min_sigma) * min(1.0, t / self.decay_period)
        return torch.clamp(action + ou_state, self.low, self.high)

#https://github.com/vitchyr/rlkit/blob/master/rlkit/exploration_strategies/ou_strategy.py

#
 Continuous control with deep reinforcement learning

#
Arxiv

# In[18]:

class ValueNetwork(nn.Module):
    def __init__(self, num_inputs, num_actions, hidden_size, init_w=3e-3):
        super(ValueNetwork, self).__init__()

        self.linear1 = nn.Linear(num_inputs + num_actions, hidden_size)
        self.linear2 = nn.Linear(hidden_size, hidden_size)
        self.linear3 = nn.Linear(hidden_size, 1)

        self.linear3.weight.data.uniform_(-init_w, init_w)
        self.linear3.bias.data.uniform_(-init_w, init_w)

    def forward(self, state, action):
        x = torch.cat([state, action], 1)
        x = F.relu(self.linear1(x))
        x = F.relu(self.linear2(x))
        x = self.linear3(x)
        return x

class PolicyNetwork(nn.Module):
    def __init__(self, num_inputs, num_actions, hidden_size, init_w=3e-3):
        super(PolicyNetwork, self).__init__()

        self.linear1 = nn.Linear(num_inputs, hidden_size)
        self.linear2 = nn.Linear(hidden_size, hidden_size)
        self.linear3 = nn.Linear(hidden_size, num_actions)

        self.linear3.weight.data.uniform_(-init_w, init_w)
        self.linear3.bias.data.uniform_(-init_w, init_w)

    def forward(self, state):
        x = F.relu(self.linear1(state))
        x = F.relu(self.linear2(x))
        x = torch.tanh(self.linear3(x))
        return x

    def get_action(self, state):
#        state  = torch.FloatTensor(state).unsqueeze(0).to(device)
        action = self.forward(state.to(device))
        return action.detach().cpu()

#
DDPG Update

# In[19]:

def ddpg_update(batch_size,
           gamma = 0.99,
           min_value=-np.inf,
           max_value=np.inf,
           soft_tau=1e-2):

    state, action, reward, next_state, done = replay_buffer.sample(batch_size)

    state      = torch.FloatTensor(state).reshape((batch_size,-1)).to(device)
    next_state = torch.FloatTensor(next_state).reshape((batch_size,-1)).to(device)
    action     = torch.FloatTensor(action).reshape((batch_size,-1)).to(device)
    reward     = torch.FloatTensor(reward).unsqueeze(1).to(device)
    done       = torch.FloatTensor(np.float32(done)).unsqueeze(1).to(device)

    policy_loss = value_net(state, policy_net(state))
    policy_loss = -policy_loss.mean()

    next_action    = target_policy_net(next_state)
    target_value   = target_value_net(next_state, next_action.detach())
    expected_value = reward + (1.0 - done) * gamma * target_value
    expected_value = torch.clamp(expected_value, min_value, max_value)

    value = value_net(state, action)
    value_loss = value_criterion(value, expected_value.detach())

    policy_optimizer.zero_grad()
    policy_loss.backward()
    policy_optimizer.step()

    value_optimizer.zero_grad()
    value_loss.backward()
    value_optimizer.step()

    for target_param, param in zip(target_value_net.parameters(), value_net.parameters()):
            target_param.data.copy_(
                target_param.data * (1.0 - soft_tau) + param.data * soft_tau
            )

    for target_param, param in zip(target_policy_net.parameters(), policy_net.parameters()):
            target_param.data.copy_(
                target_param.data * (1.0 - soft_tau) + param.data * soft_tau
            )

    return policy_loss.detach().cpu(),value_loss.detach().cpu()

# In[]:

import os
import data_manager
import pandas as pd
import numpy as np
import random
import datetime

class env_stock():
    def __init__(self,count_max,view_seq,init_money):
        self.count_max= count_max
        self.view_seq = view_seq
        self.win1= vis.line(Y=torch.Tensor([0]))
        self.win2= vis.line(Y=torch.Tensor([0]))
        self.init_money = init_money

    def reset(self):
        stock_code = np.random.choice(['005930' , '000270','000660','005380','005490','009240','009540'])

        self.prev_action = 0
        self.count = 0
        self.balance = self.init_money
        self.num_stocks = 0
        self.sum_action = 0

        chart_data = data_manager.load_chart_data(
            os.path.join('./',
                         'data/chart_data/{}.csv'.format(stock_code)))
        prep_data = data_manager.preprocess(chart_data)
        training_data = data_manager.build_training_data(prep_data)

        # 기간 필터링
        start = random.randint(self.view_seq ,(len(training_data)-self.count_max-200))

        training_data = training_data[start-self.view_seq :start+self.count_max+200]

#        training_data = training_data[(training_data['date'] >= self.start_date) &
#                                      (training_data['date'] <= self.end_date)]
        training_data = training_data.dropna()

        # 차트 데이터 분리
        features_chart_data = ['date', 'open', 'high', 'low', 'close', 'volume']
        self.chart_data = training_data[features_chart_data]
        self.chart_data['date']= pd.to_datetime(self.chart_data.date).astype(np.int64)/1e12

        # 학습 데이터 분리
        features_training_data = [
            'open_lastclose_ratio', 'high_close_ratio', 'low_close_ratio',
            'close_lastclose_ratio', 'volume_lastvolume_ratio',
            'close_ma5_ratio', 'volume_ma5_ratio',
            'close_ma10_ratio', 'volume_ma10_ratio',
            'close_ma20_ratio', 'volume_ma20_ratio',
            'close_ma60_ratio', 'volume_ma60_ratio',
            'close_ma120_ratio', 'volume_ma120_ratio'
        ]
        training_data = training_data[features_training_data]

        self.data = torch.from_numpy(training_data.values).float()

        state = self.data[self.count :self.count+self.view_seq].view(1,-1)
        state = torch.cat([state, torch.Tensor([self.sum_action]).view(1,-1)],dim=1)
        return state
        pass

    def step(self,action):
        #continuouse action space
        r_t = d_t= 0
        action = action.item()
        self.sum_action = np.clip(self.sum_action + action, 0,100)
        aa = round(self.sum_action)

        delta_a = (self.prev_action - aa)

        self.prev_action = aa

        cost = self.data[self.count,1]*delta_a*0.99
        self.pocket += cost

#        r_t = cost
#        r_t = cost
        if self.count+1 == self.count_max:
#            self.pocket += self.data[self.count,1]*(self.prev_action*quantize)
            d_t = 1
            r_t = self.pocket

        else:
            self.count +=1

        state = self.data[self.count :self.count+self.view_seq].view(1,-1)
        state = torch.cat([state, torch.Tensor([self.sum_action]).view(1,-1)],dim=1)

        return state ,r_t ,d_t,0
        pass

    def vis(self):
        self.win1= vis.line(Y=self.data[:,1],win=self.win1,opts=dict(title='price'))
        self.win2= vis.line(Y=self.data[:,2],win=self.win2,opts=dict(title='vol'))

# In[25]:

STEP_MAX =400
WINDOW_SIZE = 2

env = env_stock(STEP_MAX,WINDOW_SIZE,10000000)
env.reset()

#env = NormalizedActions(gym.make("Pendulum-v0"))
ou_noise = OUNoise(1,theta=0.1)

state_dim  = 15*WINDOW_SIZE+1
action_dim = 1
hidden_dim = 256

value_net  = ValueNetwork(state_dim, action_dim, hidden_dim).to(device)
policy_net = PolicyNetwork(state_dim, action_dim, hidden_dim).to(device)

target_value_net  = ValueNetwork(state_dim, action_dim, hidden_dim).to(device)
target_policy_net = PolicyNetwork(state_dim, action_dim, hidden_dim).to(device)

for target_param, param in zip(target_value_net.parameters(), value_net.parameters()):
    target_param.data.copy_(param.data)

for target_param, param in zip(target_policy_net.parameters(), policy_net.parameters()):
    target_param.data.copy_(param.data)

value_lr  = 1e-4
policy_lr = 1e-5

value_optimizer  = optim.Adam(value_net.parameters(),  lr=value_lr)
policy_optimizer = optim.Adam(policy_net.parameters(), lr=policy_lr)

clip = 1
torch.nn.utils.clip_grad_norm_(value_net.parameters(),clip)
torch.nn.utils.clip_grad_norm_(policy_net.parameters(),clip)

value_criterion = nn.MSELoss()

replay_buffer_size = 1000000
replay_buffer = ReplayBuffer(replay_buffer_size)

# In[28]:

max_frames  = 12000000
max_steps   = env.count_max
frame_idx   = 0
rewards     = []
batch_size  = 128
#batch_size  = 5

# In[29]:

win_p = vis.line(Y=torch.Tensor([0]),opts=dict(title='policy'))
win_v = vis.line(Y=torch.Tensor([0]), opts=dict(title='value'))
win_r = vis.line(Y=torch.Tensor([0]), opts=dict(title='reward'))
win_a = vis.line(Y=torch.Tensor([0]), opts=dict(title='action'))
win_rew = vis.line(Y=torch.Tensor([0]), opts=dict(title='reward_each_epi'))
win_pocket = vis.line(Y=torch.Tensor([0]))

print(max_steps)

while frame_idx < max_frames:
    state = env.reset()

    mem = []
    ou_noise.reset()
    episode_reward = 0
    traj = []
    traj_reward = []
    traj_pocket = []

    for step in range(max_steps):
        action = policy_net.get_action(state)
        action = ou_noise.get_action(action, step)

        traj.append(env.sum_action)

        next_state, reward, done, _ = env.step(action*1)
        traj_pocket.append(env.pocket.item())
        replay_buffer.push(state, action, reward, next_state, done)

        if len(replay_buffer) > batch_size and frame_idx%10 == 0:
            ploss,vloss=ddpg_update(batch_size)
            if frame_idx%1000 ==0:
                win_p = vis.line(X=torch.Tensor([frame_idx]), Y=ploss.view(1) , win = win_p, update='append')
                win_v = vis.line(X=torch.Tensor([frame_idx]), Y=vloss.view(1) , win = win_v,  update='append')

        state = next_state
        traj_reward.append(reward)
        episode_reward += reward
        frame_idx += 1

#        if frame_idx % max(1000, max_steps + 1) == 0:
#            plot(frame_idx, rewards)

        if done:
#            replay_buffer.push(mem)
            win_r = vis.line(X=torch.Tensor([frame_idx]), Y=torch.Tensor([episode_reward]).clamp(-1e5,1e5).view(-1) , win = win_r , update='append')
            win_pocket = vis.line(Y=torch.Tensor(traj_pocket).view(-1) , win = win_pocket ,opts=dict(title='pocket_state each step'))
            win_rew = vis.line(Y=torch.Tensor(traj_reward).view(-1) , win = win_rew,opts=dict(title='reward each step'))
            win_a = vis.line(Y=torch.Tensor(traj) , win = win_a,opts=dict(title='sum action each step'))
            env.vis()

            break

    rewards.append(episode_reward)


text:
#!/usr/bin/env python
"""
Generic python script.
"""
__author__ = "Alex Drlica-Wagner"

import functools
import matplotlib.ticker as mticker
import numpy as np

import cartopy.crs as ccrs
from cartopy.mpl.geoaxes import GeoAxes
from cartopy.geodesic import Geodesic
from shapely.geometry.polygon import Polygon

from cartosky.utils import setdefaults
from cartosky import healpix

def default_transform(self, func, *args, transform=None, **kwargs):
    """
    Makes ``transform=cartopy.crs.PlateCarree()`` the default
    for cartopy plots. This means you no longer have to
    pass ``transform=cartopy.crs.PlateCarree()`` if your data
    coordinates are longitude and latitude.
    Note
    ----
    This function wraps %(methods)s for `~proplot.axes.CartopyAxes`.
    Comes from ProPlot
    """
    # Apply default transform
    # TODO: Do some cartopy methods reset backgroundpatch or outlinepatch?
    # Deleted comment reported this issue
    if transform is None:
        transform = ccrs.PlateCarree()
    result = func(self, *args, transform=transform, **kwargs)
    return result

def default_crs(self, func, *args, crs=None, **kwargs):
    """
    Fixes the `~cartopy.mpl.geoaxes.GeoAxes.set_extent` bug associated with
    tight bounding boxes and makes ``crs=cartopy.crs.PlateCarree()`` the
    default for cartopy plots.
    Note
    ----
    This function wraps %(methods)s for `~proplot.axes.CartopyAxes`.
    Comes from ProPlot
    """
    # Apply default crs
    if crs is None:
        crs = ccrs.PlateCarree()
    try:
        result = func(self, *args, crs=crs, **kwargs)
    except TypeError as err:  # duplicate keyword args, i.e. crs is positional
        if not args:
            raise err
        args, crs = args[:-1], args[-1]
        result = func(self, *args, crs=crs, **kwargs)

    return result

def default_gridlines(self, func, *args, **kwargs):
    """
    Default gridlines options.
    """
    fmt = mticker.FuncFormatter(lambda v, pos: '{:g}'.format(v))
    defaults = dict(linestyle=':', xformatter=fmt, yformatter=fmt, draw_labels=True,
                    crs=ccrs.PlateCarree())
    setdefaults(kwargs, defaults)
    result = func(self, *args, **kwargs)
    return result

def ignore_warnings(self, func, *args, **kwargs):
    """
    Turn off warnings for this function.
    """
    import warnings

    with warnings.catch_warnings():
        warnings.simplefilter("ignore", UserWarning)
        result = func(self, *args, **kwargs)

    return result

def _wrapper_decorator(driver):
    """
    Generate generic wrapper decorator and dynamically modify the docstring
    to list methods wrapped by this function. Also set `__doc__` to ``None`` so
    that ProPlot fork of automodapi doesn't add these methods to the website
    documentation. Users can still call help(ax.method) because python looks
    for superclass method docstrings if a docstring is empty.
    """
    driver._docstring_orig = driver.__doc__ or ''
    driver._methods_wrapped = []

    def decorator(func):
        # Define wrapper and suppress documentation
        # We only document wrapper functions, not the methods they wrap
        @functools.wraps(func)
        def _wrapper(self, *args, **kwargs):
            return driver(self, func, *args, **kwargs)
        _wrapper.__doc__ = None
        return _wrapper
    return decorator

_default_transform = _wrapper_decorator(default_transform)
_default_crs = _wrapper_decorator(default_crs)
_default_gridlines = _wrapper_decorator(default_gridlines)
_ignore_warnings = _wrapper_decorator(ignore_warnings)

class SkyAxes(GeoAxes):
    """ Subclass cartopy.GeoAxes """

    def tissot(self, ra, dec, radius=1.0, n_samples=100, **kwargs):
        """
        Add Tissot's indicatrix to the axes.
        https://en.wikipedia.org/wiki/Tissot%27s_indicatrix

        Parameters
        ----------
        ra    : right asencion of circle center (deg)
        dec   : declination of circle center (deg)
        radius: angular radius (deg)
        n_samples : number of samples
        ``**kwargs`` are passed to `SkyAxes.add_geometries`.

        Returns
        -------
        feat  : feature artist
        """
        ra = np.asarray(ra).flatten()
        dec = np.asarray(dec).flatten()

        if ra.shape != dec.shape:
            msg = "ra and dec must have the same shape."
            raise ValueError(msg)

        globe = self.projection.globe

        params = globe.to_proj4_params()
        a, b = params['a'], params['b']
        geod = Geodesic(radius=a, flattening=1 - b/a)
        radius_m = a * np.radians(radius)

        geoms = []
        for x, y in zip(ra, dec):
            circle = geod.circle(x, y, radius_m, n_samples=n_samples)
            geoms.append(Polygon(circle))

        return self.add_geometries(geoms, ccrs.Geodetic(globe), **kwargs)

    def tissot_indicatrices(self, radius=5.0, num=6, **kwargs):
        """
        Add Tissot's indicatrices to the axes.
        https://en.wikipedia.org/wiki/Tissot%27s_indicatrix

        Parameters
        ----------
        radius : angular radius of circles (deg)
        ``**kwargs`` are passed to `Skymap.tissot`.

        Returns
        -------
        feat : feature artist
        """
        ra = np.linspace(-180, 180, num, endpoint=False)
        dec = np.linspace(-80, 80, num)
        return self.tissot(*np.meshgrid(ra, dec), radius=radius, **kwargs)

    # Change gridlines default to ls=':'
    gridlines = _default_gridlines(
        GeoAxes.gridlines
    )

    plot = _default_transform(
        GeoAxes.plot
    )
    scatter = _default_transform(
        GeoAxes.scatter
    )
    fill_between = _default_transform(
        GeoAxes.fill_between
    )
    fill_betweenx = _default_transform(
        GeoAxes.fill_betweenx
    )
    contour = _default_transform(
        GeoAxes.contour
    )
    contourf = _default_transform(
        GeoAxes.contourf
    )
    pcolor = _default_transform(
        GeoAxes.pcolor
    )
    pcolormesh = _default_transform(
        GeoAxes.pcolormesh
    )
    quiver = _default_transform(
        GeoAxes.quiver
    )
    streamplot = _default_transform(
        GeoAxes.streamplot
    )
    barbs = _default_transform(
        GeoAxes.barbs
    )
    tripcolor = _default_transform(
        GeoAxes.tripcolor
    )
    tricontour = _default_transform(
        GeoAxes.tricontour
    )
    tricontourf = _default_transform(
        GeoAxes.tricontourf
    )
    get_extent = _ignore_warnings(
        GeoAxes.get_extent
    )
    set_extent = _default_crs(
        GeoAxes.set_extent
    )
    set_xticks = _default_crs(
        GeoAxes.set_xticks
    )
    set_yticks = _default_crs(
        GeoAxes.set_yticks
    )

    def hpxmap(self, hpxmap, pixel=None, nside=None, xsize=800,
               lonra=None, latra=None, badval=healpix.UNSEEN, smooth=None,
               **kwargs):
        """ Draw a healpix map with pcolormesh.

        Parameters
        ----------
        hpxmap: input healpix or HealSparse map
        pixel:  explicit pixel indices in RING scheme (required for partial healpix maps)
        nside:  explicit nside of the map (required for partial healpix maps) if
                passed while visualizing a HealSparse map it will doegrade the map to this nside.
        xsize:  resolution of the output image
        lonra:  longitude range [-180,180] (deg)
        latra:  latitude range [-90,90] (deg)
        badval: set of values considered "bad"
        smooth: gaussian smoothing kernel (deg)
        kwargs: passed to pcolormesh

        Returns
        -------
        im,lon,lat,values : mpl image with pixel longitude, latitude (deg), and values
        """
        healpix.check_hpxmap(hpxmap, pixel, nside)
        hpxmap = healpix.masked_array(hpxmap, badval)

        if smooth:
            # To smooth we need the full map...
            # It'd be good to check we aren't going to blow anything up
            hpxmap = healpix.create_map(hpxmap, pixel, nside, badval)
            pixel, nside = None, None
            hpxmap = healpix.masked_array(hpxmap, badval)
            hpxmap = healpix.smooth(hpxmap, sigma=smooth)

        vmin, vmax = np.percentile(hpxmap.compressed(), [2.5, 97.5])

        defaults = dict(rasterized=True, vmin=vmin, vmax=vmax)
        setdefaults(kwargs, defaults)

        lon, lat, values = healpix.hpx2xy(hpxmap, pixel=pixel, nside=nside,
                                          xsize=xsize,
                                          lonra=lonra, latra=latra)

        im = self.pcolormesh(lon, lat, values, **kwargs)
        self._sci(im)
        return im, lon, lat, values

    def hpxbin(self, lon, lat, nside=256, **kwargs):
        """
        Create a healpix histogram of the counts.

        Like `hexbin` from matplotlib

        Parameters:
        -----------
        lon : input longitude (deg)
        lat : input latitude (deg)
        nside : heaplix nside resolution
        kwargs : passed to SkyAxes.hpxmap

        Returns:
        --------
        hpxmap, im : healpix map and image
        """
        hpxmap = healpix.hpxbin(lon, lat, nside)
        return hpxmap, self.hpxmap(hpxmap, **kwargs)

    def hspmap(self, hspmap, pixel=None, nside=None, xsize=800,
               lonra=None, latra=None, badval=healpix.UNSEEN, smooth=None,
               **kwargs):
        """ Draw a healpix map with pcolormesh.

        Parameters
        ----------
        hpxmap: input healpix or HealSparse map
        pixel:  explicit pixel indices in RING scheme (required for partial healpix maps)
        nside:  explicit nside of the map (required for partial healpix maps) if
                passed while visualizing a HealSparse map it will doegrade the map to this nside.
        xsize:  resolution of the output image
        lonra:  longitude range [-180,180] (deg)
        latra:  latitude range [-90,90] (deg)
        badval: set of values considered "bad"
        smooth: gaussian smoothing kernel (deg)
        kwargs: passed to pcolormesh

        Returns
        -------
        im,lon,lat,values : mpl image with pixel longitude, latitude (deg), and values
        """
        # ADW: probably still not the best way to do this...
        import healsparse as hsp
        import healpy as hp

        if not isinstance(hspmap, hsp.HealSparseMap):
            # Assume that it is a healpix map in RING ordering
            hpxmap = hp.reorder(hspmap, r2n=True)
            hspmap = hsp.HealSparseMap(healpix_map=hpxmap, nside_coverage=nside)

        lon, lat, values = healpix.hsp2xy(hspmap, xsize=xsize, lonra=lonra, latra=latra)

        if smooth:
            msg = "Healsparse smoothing not implemented."
            raise Exception(msg)

        vmin, vmax = np.percentile(values.compressed(), [2.5, 97.5])

        defaults = dict(rasterized=True, vmin=vmin, vmax=vmax)
        setdefaults(kwargs, defaults)

        im = self.pcolormesh(lon, lat, values, **kwargs)
        self._sci(im)
        return im, lon, lat, values

ext: py

max_stars_repo_name: erykoff/cartosky

max_issues_repo_name: erykoff/cartosky

max_forks_repo_name: erykoff/cartosky

reason: import numpy


text:
# Copyright 2021 United States Government as represented by the Administrator of the National Aeronautics and Space
# Administration.  No copyright is claimed in the United States under Title 17, U.S. Code. All Other Rights Reserved.

"""
This module implements a basic point object for GIANT, used to represent objects that cannot be traced but that exist in
a scene.

Description
-----------

:class:`.Point` objects in GIANT are essentially just a size 3 numpy array containing the location of the point in the
current frame.  It contains no geometry and will not produce anything visible when rendering a scene.  It is primarily
used to represent the light source (the sun) in the scene, though you could also use it if you are doing unresolved
navigation using :mod:`.unresolved`.
"""

from typing import Union

import numpy as np

from giant.rotations import Rotation

from giant._typing import ARRAY_LIKE

class Point:
    """
    Represents a single, unrenderable point.

    The point is just a shell around a size 3 numpy array to keep track of the location (but not orientation) of
    something in a scene that you don't actually want to render when doing ray tracing.  Therefore, this is most
    frequently used to represent the location of the sun in the scene.  Technically you can also use this if you are
    only doing :mod:`.unresolved` relative navigation; however, this is not recommended since no knowledge of the size
    of the object is available.

    To use this class simply provide the initial location as a length 3 array like object.  It can then be
    rotated/translated like other objects in a scene (though again note that it won't be rendered, and its orientation
    is not tracked, only the position).  If you are looking for something that will track both position and orientation,
    and can be rendered then we recommend checking out the :class:`.Ellipsoid` class instead.
    """

    def __init__(self, location: ARRAY_LIKE):
        """
        :param location: The starting location in the current frame for the point as a size 3 array like object.
        """

        self.location: np.ndarray = np.asarray(location).reshape(3).astype(np.float64)
        """
        The location of the point as a length 3 double array.
        """

    def rotate(self, rotation: Union[Rotation, ARRAY_LIKE]):
        """
        This rotates the point into a new frame in place.

        Only the location is rotated, since the point itself doesn't have any orientation.

        :param rotation: The rotation from the previous frame to the new frame.  For possible inputs see the
                         :class:`.Rotation` documentation.
        """

        if isinstance(rotation, Rotation):
            self.location = np.matmul(rotation.matrix, self.location)

        else:
            self.location = np.matmul(Rotation(rotation).matrix, self.location)

    def translate(self, translation: ARRAY_LIKE):
        """
        This translates the location of the point.

        :param translation: an array like object of size 3
        :raises ValueError: if the provided translation is not size 3
        """

        if np.size(translation) == 3:
            trans_array = np.asarray(translation).astype(np.float64)

            self.location += trans_array.reshape(3)

        else:
            raise ValueError("You have entered an improperly sized translation.\n"
                             "Only length 3 translations are allowed.\n"
                             "You entered {0}".format(np.size(translation)))

max_stars_repo_name: nasa/giant

max_issues_repo_name: nasa/giant

max_forks_repo_name: nasa/giant


text:
###############################
### IMPORT REQUIRED MODULES ###
###############################
import pandas as pd
import numpy as np
import tensorflow as tf

from tensorflow import keras
from tensorflow.keras import layers
from tensorflow.keras.layers import LSTM, GRU, TimeDistributed, RepeatVector, BatchNormalization, Conv1D, MaxPooling1D, GlobalAveragePooling2D
from tensorflow.keras.layers import Dense, Dropout, Reshape, Flatten, Input, GaussianNoise, Attention, Bidirectional, Embedding, Dot, GlobalMaxPooling2D
from tensorflow.keras.layers import Concatenate, Multiply, Add, Activation, ReLU, ConvLSTM2D, Conv2D, MaxPooling2D, Conv3D, MaxPooling3D
from tensorflow.keras.models import Sequential
from tensorflow.keras import Model
from Transformer import TransformerBlock

###########################
### DEFINE DF-TAR MODEL ###
###########################
def DF_TAR(n_steps, length, n_districts, n_features, conv_hidden=256, gru_hidden=256, fc_hidden=1024, batchnorm=True):
    # Input Layers
    risk = Input(shape=[length, n_districts, n_features['risk']])
    E_static = Input(shape=[length, n_districts, n_features['E_static']])
    E_dynamic = Input(shape=[length, n_districts, n_features['E_dynamic']])
    c = Input(shape=[length, n_districts, n_features['c']])

    # Convolutional Block for Static Environmental Features
    h_static = Conv2D(conv_hidden, 9, padding='same', activation='relu')(E_static)
    h_static = MaxPooling2D(3, padding='same')(h_static)
    h_static = Conv2D(conv_hidden, 4, padding='same', activation='relu')(h_static)
    h_static = MaxPooling2D(3, padding='same')(h_static)
    h_static = Conv2D(conv_hidden, 3, padding='same', activation='relu')(h_static)
    h_static = MaxPooling2D(3, padding='same')(h_static)
    h_static = Conv2D(conv_hidden, 3, padding='same', activation='relu')(h_static)
    h_static = GlobalMaxPooling2D()(h_static)
    h_static = Flatten()(h_static)
    h_static = Dense(conv_hidden)(h_static)
    h_static = RepeatVector(length)(h_static)

    # Recurrent Block for Historical Risk Scores
    h_risk = Reshape([length, n_districts * n_features['risk']])(risk)
    h_risk = TransformerBlock(n_districts * n_features['risk'], n_features['risk'], 256)(h_risk)
    h_risk = Concatenate()([h_risk, h_static])
    h_risk = GRU(gru_hidden)(h_risk)

    # Recurrent Block for Dynamic Environmental Features
    h_dynamic = Reshape([length, n_districts * n_features['E_dynamic']])(E_dynamic)
    h_dynamic = TransformerBlock(n_districts * n_features['E_dynamic'], n_features['E_dynamic'], 256)(h_dynamic)
    h_dynamic = Concatenate()([h_dynamic, h_static])
    h_dynamic = GRU(gru_hidden)(h_dynamic)

    # Recurrent Block for Dangerous Driving Statistics
    h_case = Reshape([length, n_districts * n_features['c']])(c)
    h_case = TransformerBlock(n_districts * n_features['c'], n_features['c'], 256)(h_case)
    h_case = Concatenate()([h_case, h_static])
    h_case = GRU(gru_hidden)(h_case)

    # Fusion Block
    h_concat = Concatenate()([h_dynamic, h_risk, h_case])
    tanh = Activation('tanh')(h_concat)
    sigmoid = Activation('sigmoid')(h_concat)
    h_fusion = Multiply()([tanh, sigmoid])

    # Fully-Connected Block for Final Prediction
    f = Dense(fc_hidden)(h_fusion)
    if batchnorm:
        f = BatchNormalization()(f)
    f = Activation('relu')(f)

    f = Dense(fc_hidden)(f)
    if batchnorm:
        f = BatchNormalization()(f)
    f = Activation('relu')(f)

    f = Dense(fc_hidden)(f)
    if batchnorm:
        f = BatchNormalization()(f)
    f = Activation('relu')(f)

    # Predicted Risk Scores optimized by MAE
    y_mae = Dense(n_steps * n_districts)(f)
    y_mae = Reshape([n_steps, n_districts], name='Y_MAE')(y_mae)

    # Predicted Risk Scores optimized by MSE
    y_mse = Dense(n_steps * n_districts)(f)
    y_mse = Reshape([n_steps, n_districts], name='Y_MSE')(y_mse)

    return Model(inputs=[risk, E_static, E_dynamic, c], outputs=[y_mae, y_mse])


text:
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name:        Raytheon Company
# Contractor Address:     6825 Pine Street, Suite 340
#                         Mail Stop B8
#                         Omaha, NE 68106
#                         402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
#----------------------------------------------------------------------------
# This software is in the public domain, furnished "as is", without technical
# support, and with no warranty, express or implied, as to its usefulness for
# any purpose.
#
# SnowAmtQPFPoPWxCheck
#
# Author: Jay Smith, WFO Fairbanks, jay.smith@noaa.gov, 907-458-3721
# Version: 1.0.0, 09/14/2006 - Initial version
#          1.0.1, 10/12/2006 - Added PoP/QPF check at request of DSPAC
#          1.0.2, 10/18/2006 - Changed PoP/QPF check to treat the PoP as
#                 floating. Instead of checking each individual PoP grid
#                 against its corresponding QPF grid, the max of all the
#                 PoP grids overlapping a QPF grid will be checked.
#          1.1.0, 01/25/2007 - Added options to choose which checks to run.
#                 Reorganized code so that each check is its own method.
#                 Added a check for QPF and Wx. Added highlighting for the
#                 created temporary grids.
#          1.1.1, 02/01/2007 - Changed the SnowAmt/Wx check to return
#                 consistent results for SnowAmt > 0 and Wx grid containing
#                 S, SW, or IP regardless of whether the frozen precip is
#                 mixed with freezing and/or liquid precip.
#          1.2.0, 02/13/2007 - Added a configuration option to provide a CWA
#                 edit area to run the procedure over. A bad edit area or no
#                 edit area will result in running over the whole domain.
#                 Modified the SnowAmt/Wx and QPF/Wx checks to handle two
#                 cases. Case 1: The SnowAmt/QPF grid is 6-hr long and starts
#                 at 00, 06, 12, or 18 UTC. Then only one of the corresponding
#                 Wx grids has to meet the consistency rule. Case 2: The
#                 SnowAmt/QPF grid does not meet the case 1 definition. Then
#                 all of the corresponding Wx grids must meet the consistency
#                 rule.
# The procedure performs the following checks:
# 1. If SnowAmt present and >= 0.5 inches, then corresponding QPF grids
#    must add up to 0.01 inches.
# 2. If SnowAmt >= 0.1 inches, then there are two cases:
#    a. If the SnowAmt grid is exactly 6 hours long and starts at 00, 06, 12,
#       or 18 UTC, then at least one of the corresponding Wx grids must have
#       S, SW, or IP.
#    b. If the SnowAmt grid does not adhere to the time constraints listed in
#       in the previous paragraph, then all of the corresponding Wx grids
#       must have S, SW, or IP. This more stringent test is required because
#       with grids offset from the NDFD time constraints, it's possible for
#       the GFE to evaluate the grids as consistent using an "any"
#       criteria but have the NDFD flag those same grids as inconsistent.
# 3. If QPF > 0, then at least one of the corresponding PoP grids must be > 0
# 4. If QPF > 0, then there are two cases:
#    a. If the QPF grid is exactly 6 hours long and starts at 00, 06, 12, or 18
#       UTC, then at least one of the corresponding Wx grids must have R, RW,
#       S, SW, RS, IP, L, ZR, ZL.
#    b. If the QPF grid does not adhere to the time constraints listed in the
#       previous paragraph, then all corresponding Wx grids must contain a
#       precipitating weather type. This more stringent test is required
#       because with grids offset from the NDFD time constraints, it's
#       possible for the GFE to evaluate grids as consistent using an "any"
#       criteria but have the NDFD flag those same grids as inconsistent.
# For all of the checks above, if the initial threshold is not exceeded, then
# the two grids are consistent by definition. In other words:
# 1. If SnowAmt < 0.5, then SnowAmt and QPF are always consistent.
# 2. If SnowAmt < 0.1, then SnowAmt and Wx are always consistent.
# 3. If QPF = 0, then QPF and PoP are always consistent.
# 4. If QPF = 0, then QPF and Wx are always consistent.
# For the Wx checks above, only the Wx type is considered.
#
# ****** NOTE NOTE NOTE NOTE ******
# At this time, the check for two 6-hour QPF grids vs. one 12-hr PoP grid
# is not implemented because neither of those grid definitions is implemented
# in the GFE baseline. I don't know how to do a check on grids that don't
# exist.
# ****** NOTE NOTE NOTE NOTE ******
#
# If discrepancies are found, then the "bad" grids will be highlighted.
# Temporary grids showing where the discrepancies occur will be created and
# also highlighted.
#
# Dealing with QPF and SnowAmt is always a pain, because they are "cumulative"
# elements. This procedure will account for the possibility that the SnowAmt and
# QPF grids are not the same duration. It will also account for the possibilty
# that the SnowAmt and QPF grids are not aligned on either or both ends.
# The only sane way to handle either situation is to believe that the QPF
# accumulation happens uniformally across the grid's duration and to use
# the proportional amount of the QPF that corresponds the SnowAmt grid's
# duration. Some examples:
# 1. The QPF grid is 3 hours long and there are 3, 1-hour, SnowAmt grids.
# Each SnowAmt grid will be compared to 1/3 the value of the QPF grid.
# 2. The last two hours of a 3-hour QPF grid overlaps a 2-hour SnowAmt grid.
# The SnowAmt grid will be compared to 2/3 the value of the QPF grid.
# 3. Two 3-hour QPF grids align with one 6-hour SnowAmt grid. The first QPF
# grid will be compared to the SnowAmt grid. If the consistency check passes
# on that comparison, the program will continue. If the consistency check
# fails, then the sum of the two QPF grids will be compared to the SnowAmt
# grid.
# 4. The last four hours of a 6-hour QPF grid and the first two hours of a
# 3-hour QPF grid overlap a 6-hour SnowAmt grid. The SnowAmt grid will be
# compared to 2/3 of the first QPF grid. If the consistency check passes,
# the program will continue. If the consistency check fails, then 2/3 of the
# first QPF grid will be added to 2/3 of the second QPF grid and that QPF
# sum will be compared against the SnowAmt grid.
#
# Confused yet? Of course, all of these gyrations can be avoided if the
# QPF and SnowAmt grids are aligned and of the same duration.
#
# Unfortunately, the GFE does not provide a way to deal with proportional
# amounts of the accumulative grids, so I have done this.
#
# I've written this code such that it's optimized to minimize memory usage
# (at least I think I've done that). As a result, it's not particularly
# optimized for ifpServer database access. In fact, I retrieve the various
# grids from the ifpServer database many times during the procedure's run.
# This will have an impact on how fast the procedure runs (it'll run slower
# than if I had optimized for ifpServer database access). The choice to favor
# memory optimization comes from my belief that there are still "memory leak"
# problems in the GFE and that the consequences of those problems will be most
# manifest when this procedure is most likely to be run (near the end of the
# shift). Funky memory problems are a prime cause of funky application
# behavior like application crashes or spontaneous logouts. So, this procedure
# basically reads a grid into memory, keeps it as long as it's needed, and
# then discards it.
#
# Finally, this procedure is also intended to provide an example to other
# developers of how to write and document code. I have reservations as to how
# well I've succeeded at that task. The code is heavily documented, probably
# excessively so. Also, it's not as well as organized as it could be. As you
# look through the various methods, it should become quickly apparent that
# there is a lot of repeated code. I've consciously left the code this way in
# the hopes that it will be easier to understand by more novice programmers
# and because the code hasn't quite grown to the point where updating the
# repeating code is onerous or overly error-prone. It would be better to
# capture the repeating code in separate methods, but keeping track of the
# where you are in the code becomes harder the more you have to jump around
# from method to method. As with all things, there are trade-offs involved.
# ----------------------------------------------------------------------------

##
# This is an absolute override file, indicating that a higher priority version
# of the file will completely replace a lower priority version of the file.
##

MenuItems = ["Consistency"]

VariableList = []
VariableList.append(('Check_Cleanup', 'Check', 'radio', ['Check', 'Cleanup']))
VariableList.append(('Run SnowAmt/QPF Check?', ['Yes'], 'check', ['Yes']))
VariableList.append(('Run SnowAmt/Wx Check?', ['Yes'], 'check', ['Yes']))
VariableList.append(('Run QPF/PoP Check?', ['Yes'], 'check', ['Yes']))
VariableList.append(('Run QPF/Wx Check?', ['Yes'], 'check', ['Yes']))
VariableList.append(('If "Cleanup" is selected, then only cleanup actions will run.\nNo checks will be made, regardless of the above settings.', '', 'label'))

#### Config section
# Both the QPF and SnowAmt grids have values which are floating point
# numbers. This means comparisons must use a tolerance value. In other
# words, 0.5 may be represented in machine numbers as 0.49999999999 or
# 0.500000000001. By specifying a tolerance value, we account for the
# vagaries of machine representation of floating point numbers while
# keeping the precision of the comparisons to acceptable levels. Depending
# on the comparison being done, the tolerance value will be added to or
# subtracted from the comparison value to allow for machine error in the
# floating point number representation.
# By default in the GFE, QPF precision is to the nearest one-hundredth while
# SnowAmt precision is to the nearest tenth.
qpfTol = 0.00001  # 1/100,000 tolerance vs 1/100 precision
snowAmtTol = 0.0001  # 1/10,000 tolerance vs 1/10 precision
# Inconsistent grid highlight color. One size fits all. To turn off
# highlighting, set the variable to the empty string, ''.
inconGridColor = 'red'
# Temporary grid highlight color. One size fits all. To turn off highlighting,
# set the variable to the empty string, ''.
tempGridColor = 'orange'
# Name of CWA edit area to use instead of running the procedure over the
# whole domain. Set to the empty string, '', if you want the procedure to
# always run over the whole domain. If the procedure has a problem with the
# edit area you provide, it will run over the whole domain. You should probably
# choose an edit area that is slightly larger than your entire CWA. It's
# possible that when mapping your GFE grids to NDFD grids that the NDFD thinks
# some GFE grid cells are in your CWA that the GFE does not think are in your
# CWA. Using an edit area slightly larger than the CWA, like the ISC_Send_Area
# which is the mask used when sending grids to the NDFD, should eliminate the
# possibibilty of the NDFD intermittently flagging CWA border "points" as
# inconsistent. Note: running the procedure over a subset of the entire GFE
# domain does not really provide any performance gains. Given the way the
# underlying array data structure works, calculations are almost always made
# at every single grid point first and then a mask is applied to limit the
# meaningful results to the edit area. For the purposes of this procedure, the
# values outside the edit area are set to the appropriate "consistent" result.
# The real benefit of this option is it limits the inconsistent results to the
# areas the forecaster really cares about, which should lessen the workload of
# using this procedure. Marine Offices: Make sure the edit area provided
# includes your marine zones.
cwaEditArea = 'ISC_Send_Area'
#### Config section end

import SmartScript
from numpy import *

class Procedure (SmartScript.SmartScript):
    def __init__(self, dbss):
        SmartScript.SmartScript.__init__(self, dbss)

    def __cleanup(self, timeRange):
        # Remove any temporary grids created previously.
        for element in (
            'SnowAmtQPFInconsistent', 'SnowAmtWxInconsistent',
            'QPFPoPInconsistent', 'QPFWxInconsistent'):
            try:
                # From SmartScript
                self.unloadWE('Fcst', element, 'SFC')
            except:
                # A failure is almost certainly no grids to unload.
                pass
        # Turn off any highlights. From SmartScript
        self.highlightGrids('Fcst', 'SnowAmt', 'SFC', timeRange, inconGridColor, on=0)
        self.highlightGrids('Fcst', 'QPF', 'SFC', timeRange, inconGridColor, on=0)
        self.highlightGrids('Fcst', 'Wx', 'SFC', timeRange, inconGridColor, on=0)
        self.highlightGrids('Fcst', 'PoP', 'SFC', timeRange, inconGridColor, on=0)
        return

    def __checkConfigValueTypes(self):
        import types
        message = ''
        badValues = False
        if not type(inconGridColor) is types.StringType:
            message = '%sThe "inconGridColor" variable is not defined as a string value. Please contact your IFPS focal point to fix this problem.\n' % message
            badValues = True
        if not type(tempGridColor) is types.StringType:
            message = '%sThe "tempGridColor" variable is not defined as a string value. Please contact your IFPS focal point to fix this problem.\n' % message
            badValues = True
        if not type(cwaEditArea) is types.StringType:
            message = '%sThe "cwaEditArea" variable is not defined as a string value. Please contact your IFPS focal point to fix this problem.\n' % message
            badValues = True
        if badValues:
            message = '%sYou will not be able to run the procedure until the problem is corrected.' % message
            # The next two commands are from SmartScript
            self.statusBarMsg(message, 'U')
            self.cancel()
        return

    def _runSnowAmtQPFCheck(self, timeRange):
        # This method implements the check that if SnowAmt >= 0.5, then
        # QPF must be >= 0.01.

        # There can be a significant difference between the values stored
        # in memory and the values returned from the database. This is because
        # when values are saved, the element's precision (as defined in
        # serverConfig.py/localConfig.py) is enforced. Values in memory do not
        # have the element's precision enforced; in fact, they have the
        # machine precision of the underlying data type.
        # If there are locks, post an urgent message and return from the method.
        message = ''
        # lockedByMe is from SmartScript
        if self.lockedByMe('QPF', 'SFC'):
            message = '%sYou have the QPF grid locked. Please save the QPF grid.\n' % message
        if self.lockedByMe('SnowAmt', 'SFC'):
            message = '%sYou have the SnowAmt grid locked. Please save the SnowAmt grid.\n' % message
        # lockedByOther is from SmartScript
        if self.lockedByOther('QPF', 'SFC'):
            message = '%sThe QPF grid is locked by someone else. Please have that person save the QPF grid.\n' % message
        if self.lockedByOther('SnowAmt', 'SFC'):
            message = '%sThe SnowAmt grid is locked by someone else. Please have that person save the SnowAmt grid.\n' % message
        if message:
            message = '%sThe SnowAmt/QPF Check was not run.' % message
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have locked grid problems.
            return

        # Make sure there are actually SnowAmt grids in the time range.
        # The self.getGrids command will return None if there are no grids
        # in the time range for mode='First' and noDataError=0. The None
        # variable cannot be iterated over. Rather than trap in a try/except,
        # I'll just check for the condititon. This may not be the most
        # Pythonic way of doing things, but it allows me to avoid having
        # a bunch of code indented beneath a try statement. If no SnowAmt
        # grids are found, post an urgent message and return from the method.
        # getGrids is from SmartScript
        snowAmtInfoList = self.getGridInfo('Fcst', 'SnowAmt', 'SFC', timeRange)
        if [] == snowAmtInfoList:
            message = 'There are no SnowAmt grids in the time range you selected.\nThe SnowAmt/QPF Check did not run.'
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have missing grid problems.
            return

        # getGridInfo is from SmartScript
        # One might ask why I don't just return the result of self.getGrids
        # to a variable and iterate over that. I'm trying to minimize the
        # memory footprint of the procedure. Reading all the grids into a
        # variable could be a fairly large memory hit. The construct below
        # only reads one SnowAmt grid at a time into memory, the one that's
        # being checked. By using the cache=0 switch on all the self.getGrids
        # command, I prevent the GFE from saving the grids into memory for me.
        # The Python builtin command enumerate loops over an iterable object
        # and returns a 2-tuple containing the current index of the
        # iteration and the object at that index. In cases where I need
        # both the index and the object, I think this construct is more
        # elegant than:
        # for i in xrange(len(iterableObject)):
        #     object = iterableObject[i]
        snowAmtGrids = self.getGrids('Fcst', 'SnowAmt', 'SFC',
            timeRange, mode='List', noDataError=0,cache=0)
        for snowAmtIndex, snowAmtGrid in enumerate(snowAmtGrids):
            # greater_equal is from Numeric. For the given array and
            # threshold, a new array of the same dimensions as the input
            # array is returned. The new array has the value 1 where the
            # input array was greater than or equal to the threshold and
            # has the value 0 elsewhere.
            halfInchMask = greater_equal(snowAmtGrid, 0.5 - snowAmtTol)
            gridTR = snowAmtInfoList[snowAmtIndex].gridTime()
            # zeros is from Numeric. It creates an array of all zeros for
            # the given dimensions and numeric type.
            qpfSum = self.empty()
            qpfGrids = self.getGrids(
                'Fcst', 'QPF', 'SFC', gridTR, mode='List', noDataError=0,
                cache=0)
            if qpfGrids is None:
                message = '''There are no QPF grids in time range %s.
The SnowAmt/QPF Check skipped the time range.''' % gridTR
                self.statusBarMsg(message, 'U')
                continue
            qpfInfoList = self.getGridInfo('Fcst', 'QPF', 'SFC', gridTR)
            for qpfIndex, qpfGrid in enumerate(qpfGrids):
                snowAmtGridStartTime = gridTR.startTime().unixTime()
                qpfGridTR = qpfInfoList[qpfIndex].gridTime()
                qpfGridStartTime = qpfGridTR.startTime().unixTime()
                fraction = 1.0
                if qpfGridStartTime < snowAmtGridStartTime:
                    diff = snowAmtGridStartTime - qpfGridStartTime
                    fraction -= (float(diff) / qpfGridTR.duration())
                snowAmtGridEndTime = gridTR.endTime().unixTime()
                qpfGridEndTime = qpfGridTR.endTime().unixTime()
                if qpfGridEndTime > snowAmtGridEndTime:
                    diff = qpfGridEndTime - snowAmtGridEndTime
                    fraction -= (float(diff) / qpfGridTR.duration())
                # For some reason, the construct:
                # qpfSum = qpfSum + (qpfGrid * fraction)
                # doesn't assign the expression evaluation back to qpfSum.
                # Thus, I use a temporary variable.
                qpfTemp = qpfSum + (qpfGrid * fraction)
                qpfSum = qpfTemp
                del qpfTemp
                # less is from Numeric. It behaves analogously to greater_equal,
                # described above.
                qpfMask = less(qpfSum, 0.01 + qpfTol)
                # The following is the "truth" table for the logical
                # comparison.
                # SnowAmt >= 0.5, 1; SnowAmt < 0.5, 0
                # QPF < 0.01, 1; QPF >= 0.01, 0
                # SnowAmt >= 0.5 (1) and QPF < 0.01 (1) = 1 (Bad result)
                # SnowAmt >= 0.5 (1) and QPF >= 0.01 (0) = 0 (Good result)
                # SnowAmt < 0.5 (0) and QPF < 0.01 (1) = 0 (Good result)
                # SnowAmt < 0.5 (0) and QPF >= 0.01 (0) = 0 (Good result)
                # logical_and is from Numeric
                consistMask = logical_and(halfInchMask, qpfMask)
                # Now, apply the CWA mask. There's an assumption here that
                # all offices will use a mask and provide a valid one, which
                # means this step does something meaningful. If that assumption
                # does not hold, then the next statement doesn't actually
                # change anything, even though each and every grid point has a
                # comparison check made.
                # where is from Numeric. The first argument is a mask.
                # The second argument is/are the value/values to use at the
                # array points where the mask is one. The third argument
                # is/are the value/values to use at the array points
                # where the mask is zero. For this comparison, I want
                # the values of consistMask where self.cwaMask is one and
                # I want the "good result", which is zero, where
                # self.cwaMask is zero.
                consistMask[logical_not(self.cwaMask)] = 0
                # ravel and sometrue are from Numeric.
                if not sometrue(ravel(consistMask)):
                    # This is the good result, even though it may not be
                    # intuitive. The ravel function reduces the rank of the
                    # array by one. Since we had a 2-d array, the ravel
                    # function creates a 1-d array (a vector) such that
                    # reading the 2-d array from left-to-right, top-to-
                    # bottom returns the same values as reading the 1-d
                    # array from left-to-right. The sometrue function
                    # performs a logical or on subsequent element pairs
                    # in the 1-d array and returns the final result. If
                    # there's no inconsistency, the result will be 0.
                    # Thus, negating the sometrue result gives us the
                    # positive outcome. Phew.
                    # Since QPF is an accumulative element, we don't need
                    # to continue the loop once the QPF sum meets the
                    # threshold.
                    break
            else:
                # This block will only execute if the for loop runs to
                # completion, i.e., the break statement is not executed.
                # So, if we get here, we have an inconsistency and need to
                # highlight the appropriate grids.
                if inconGridColor:
                    self.highlightGrids(
                        'Fcst', 'SnowAmt', 'SFC', gridTR, inconGridColor)
                    self.highlightGrids(
                        'Fcst', 'QPF', 'SFC', gridTR, inconGridColor)
                # createGrid is from SmartScript
                # Since this block of code only executes if the for loop
                # runs to completion, then the value of consistMask from
                # the for loop will contain all of the inconsistencies.
                self.createGrid(
                    'Fcst', 'SnowAmtQPFInconsistent', 'SCALAR', consistMask,
                    gridTR, descriptiveName='SnowAmtQPFInconsistent',
                    minAllowedValue=0, maxAllowedValue=1, units='Good/Bad')
                if tempGridColor:
                    self.highlightGrids(
                        'Fcst', 'SnowAmtQPFInconsistent', 'SFC', gridTR,
                        tempGridColor)
                self.inconsistent = True
        # While not required, I like to terminate my methods with a return
        # statement to make it clear this is where the method ends.
        return

    def _runSnowAmtWxCheck(self, timeRange):
        # This implements the check that if SnowAmt >= 0.1, then the Wx grid
        # must contain S, SW, or IP, regardless of whether or not there is
        # any freezing or liquid types. Finally, the check does not look at
        # anything other than the Wx type. In other words, the check will be
        # okay if SnowAmt != 0 and Wx has Chc:S:- or Def:SW:-- or Lkly:S:+.

        # There can be a significant difference between the values stored
        # in memory and the values returned from the database. This is because
        # when values are saved, the element's precision (as defined in
        # serverConfig.py/localConfig.py) is enforced. Values in memory do not
        # have the element's precision enforced; in fact, they have the
        # machine precision of the underlying data type.
        # If there are locks, post an urgent message and return from the method.
        message = ''
        # lockedByMe is from SmartScript
        if self.lockedByMe('Wx', 'SFC'):
            message = '%sYou have the Wx grid locked. Please save the Wx grid.\n' % message
        if self.lockedByMe('SnowAmt', 'SFC'):
            message = '%sYou have the SnowAmt grid locked. Please save the SnowAmt grid.\n' % message
        # lockedByOther is from SmartScript
        if self.lockedByOther('Wx', 'SFC'):
            message = '%sThe Wx grid is locked by someone else. Please have that person save the Wx grid.\n' % message
        if self.lockedByOther('SnowAmt', 'SFC'):
            message = '%sThe SnowAmt grid is locked by someone else. Please have that person save the SnowAmt grid.\n' % message
        if message:
            message = '%sThe SnowAmt/Wx Check was not run.' % message
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have locked grid problems.
            return

        # Make sure there are actually SnowAmt grids in the time range.
        # The self.getGrids command will return None if there are no grids
        # in the time range for noDataError=0. The None
        # variable cannot be iterated over. Rather than trap in a try/except,
        # I'll just check for the condititon. This may not be the most
        # Pythonic way of doing things, but it allows me to avoid having
        # a bunch of code indented beneath a try statement. If no SnowAmt
        # grids are found, post an urgent message and return from the method.
        # getGrids is from SmartScript
        snowAmtInfoList = self.getGridInfo('Fcst', 'SnowAmt', 'SFC', timeRange)
        if [] == snowAmtInfoList:
            message = 'There are no SnowAmt grids in the time range you selected.\nThe SnowAmt/Wx Check did not run.'
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have missing grid problems.
            return

        snowAmtGrids = self.getGrids(
            'Fcst', 'SnowAmt', 'SFC', timeRange, mode='List', noDataError=0,
            cache=0)
        for snowAmtIndex, snowAmtGrid in enumerate(snowAmtGrids):
            nonZeroMask = greater_equal(snowAmtGrid, 0.1 - snowAmtTol)
            gridTR = snowAmtInfoList[snowAmtIndex].gridTime()

            wxInfoList = self.getGridInfo('Fcst', 'Wx', 'SFC', gridTR)
            if [] == wxInfoList:
                message = '''There are no Wx grids in time range %s.
The SnowAmt/Wx Check skipped the time range.''' % gridTR
                self.statusBarMsg(message, 'U')
                continue
            # There are two cases, which I'll capture in individual methods
            # If the SnowAmt grid is exactly 6 hours long and starts at
            # 00, 06, 12, or 18 UTC, then only one overlapping Wx grid needs
            # to match. Otherwise, all overlapping Wx grids need to match.
            if gridTR.duration() / 3600 == 6 and \
               gridTR.startTime().hour in (0, 6, 12, 18):
                self._snowAmtWxCheckLocked(nonZeroMask, gridTR, wxInfoList)
            else:
                self._snowAmtWxCheckUnlocked(nonZeroMask, gridTR, wxInfoList)
        return

    def _snowAmtWxCheckLocked(self, nonZeroMask, gridTR, wxInfoList):
        # The "Locked" comes from the idea that if the SnowAmt grid meets
        # the duration and start time constraints, then it's been "locked".
        # I need to capture the consistency masks for each individual Wx grid
        # just in case I end up with inconsistencies.
        consistMaskList = []
        for wxIndex, wxGrid in enumerate(self.getGrids(
            'Fcst', 'Wx', 'SFC', gridTR, mode='List', noDataError=0,
            cache=0)):
            # wxMask is from SmartScript
            sMask = self.wxMask(wxGrid, ':S:')
            swMask = self.wxMask(wxGrid, ':SW:')
            ipMask = self.wxMask(wxGrid, ':IP:')
            snowMask = logical_or(logical_or(sMask, swMask), ipMask)
            del (sMask, swMask, ipMask)
            wxMask = logical_not(snowMask)
            # "Truth" table for the logical comparison follows
            # SnowAmt >= 0.1, 1; SnowAmt < 0.1, 0
            # Wx has S, SW, or IP, 0; Wx doesn't have S, SW, or IP, 1
            # SnowAmt >= 0.1 (1) and Wx has (0) = 0 (Good result)
            # SnowAmt >= 0.1 (1) and Wx doesn't have (1) = 1 (Bad result)
            # SnowAmt < 0.1 (0) and Wx has (0) = 0 (Good result)
            # SnowAmt < 0.1 (0) and Wx doesn't have (1) = 0 (Good result)
            #
            consistMask = logical_and(nonZeroMask, wxMask)
            consistMask[logical_not(self.cwaMask)] = 0
            consistMaskList.append(consistMask)
            if not sometrue(ravel(consistMask)):
                # There were no inconsistencies with this Wx grid. Since only
                # one needs to be consistent, we don't need to do any more
                # checks.
                break
        else:
            # This block will only execute if the for loop runs to
            # completion, i.e., the break statement is not executed.
            # So, if we get here, we have an inconsistency and need to
            # highlight the appropriate grids.
            if inconGridColor:
                self.highlightGrids(
                    'Fcst', 'SnowAmt', 'SFC', gridTR, inconGridColor)
                self.highlightGrids(
                    'Fcst', 'Wx', 'SFC', gridTR, inconGridColor)
            # createGrid is from SmartScript
            for index in xrange(len(wxInfoList)):
                # Create temporary grids for each Wx grid. Limit the start and
                # end times of the temporary grids so that they don't extend
                # beyond the start and end times of the corresponding SnowAmt
                # grid.
                wxGridTR = wxInfoList[index].gridTime()
                tempGridStartTime = wxGridTR.startTime().unixTime()
                if tempGridStartTime < gridTR.startTime().unixTime():
                    tempGridStartTime = gridTR.startTime().unixTime()
                tempGridEndTime = wxGridTR.endTime().unixTime()
                if tempGridEndTime > gridTR.endTime().unixTime():
                    tempGridEndTime = gridTR.endTime().unixTime()
                tempGridDur = (tempGridEndTime - tempGridStartTime) / 3600
                offset = (tempGridStartTime - \
                          self.timeRange0_1.startTime().unixTime()) / 3600
                # Because the time range may be different for the temporary
                # grid, I need to create and use that time range when
                # creating the temporary grid.
                tempGridTR = self.createTimeRange(
                    offset, offset+tempGridDur, 'Zulu')
                self.createGrid(
                    'Fcst', 'SnowAmtWxInconsistent', 'SCALAR',
                    consistMaskList[index], tempGridTR,
                    descriptiveName='SnowAmtWxInconsistent',
                    minAllowedValue=0, maxAllowedValue=1, units='Good/Bad')
                if tempGridColor:
                    self.highlightGrids(
                        'Fcst', 'SnowAmtWxInconsistent', 'SFC', gridTR,
                        tempGridColor)
            self.inconsistent = True
        return

    def _snowAmtWxCheckUnlocked(self, nonZeroMask, gridTR, wxInfoList):
        # The "Unlocked" comes from the idea that if the SnowAmt grid does
        # not meet the duration and start time constraints, then it's been
        # left "unlocked".
        for wxIndex, wxGrid in enumerate(self.getGrids(
            'Fcst', 'Wx', 'SFC', gridTR, mode='List', noDataError=0,
            cache=0)):
            # wxMask is from SmartScript
            sMask = self.wxMask(wxGrid, ':S:')
            swMask = self.wxMask(wxGrid, ':SW:')
            ipMask = self.wxMask(wxGrid, ':IP:')
            snowMask = logical_or(logical_or(sMask, swMask), ipMask)
            del (sMask, swMask, ipMask)
            wxMask = logical_not(snowMask)
            # "Truth" table for the logical comparison follows
            # SnowAmt >= 0.1, 1; SnowAmt < 0.1, 0
            # Wx has S, SW, or IP, 0; Wx doesn't have S, SW, or IP, 1
            # SnowAmt >= 0.1 (1) and Wx has (0) = 0 (Good result)
            # SnowAmt >= 0.1 (1) and Wx doesn't have (1) = 1 (Bad result)
            # SnowAmt < 0.1 (0) and Wx has (0) = 0 (Good result)
            # SnowAmt < 0.1 (0) and Wx doesn't have (1) = 0 (Good result)
            #
            # All Wx grids overlapping the SnowAmt grid must be consistent.
            consistMask = logical_and(nonZeroMask, wxMask)
            consistMask[logical_not(self.cwaMask)] = 0
            if sometrue(ravel(consistMask)):
                # I'll highlight the SnowAmt grids and Wx grids in
                # gridTR as I did with QPF. However, I'll make
                # temporary grids here using the Wx grid's time
                # range but, the temporary grid cannot start before
                # the start of the corresponding SnowAmt grid nor can
                # it end after the end of the corresponding SnowAmt grid.
                wxGridTR = wxInfoList[wxIndex].gridTime()
                tempGridStartTime = wxGridTR.startTime().unixTime()
                if tempGridStartTime < gridTR.startTime().unixTime():
                    # Clip to start of SnowAmt grid
                    tempGridStartTime = gridTR.startTime().unixTime()
                tempGridEndTime = wxGridTR.endTime().unixTime()
                if tempGridEndTime > gridTR.endTime().unixTime():
                    # Clip to end of SnowAmtGrid
                    tempGridEndTime = gridTR.endTime().unixTime()
                tempGridDur = (tempGridEndTime - tempGridStartTime) / 3600
                offset = (tempGridStartTime - \
                          self.timeRange0_1.startTime().unixTime()) / 3600
                # Since either the front or end of the Wx grid's
                # time range may have been clipped, create a time
                # range using those values.
                tempGridTR = self.createTimeRange(
                    offset, offset+tempGridDur, 'Zulu')
                self.createGrid(
                    'Fcst', 'SnowAmtWxInconsistent', 'SCALAR', consistMask,
                    tempGridTR, descriptiveName='SnowAmtWxInconsistent',
                    minAllowedValue=0, maxAllowedValue=1, units='Good/Bad')
                if tempGridColor:
                    self.highlightGrids(
                        'Fcst', 'SnowAmtWxInconsistent', 'SFC', gridTR,
                        tempGridColor)
                if inconGridColor:
                    self.highlightGrids(
                        'Fcst', 'SnowAmt', 'SFC', gridTR, inconGridColor)
                    self.highlightGrids(
                        'Fcst', 'Wx', 'SFC', wxGridTR, inconGridColor)
                self.inconsistent = True
        return

    def _runQPFPoPCheck(self, timeRange):
        # This method implements the check that if any QPF grid is non zero
        # then one of the corresponding floating PoP grids must also be non
        # zero.

        # There can be a significant difference between the values stored
        # in memory and the values returned from the database. This is because
        # when values are saved, the element's precision (as defined in
        # serverConfig.py/localConfig.py) is enforced. Values in memory do not
        # have the element's precision enforced; in fact, they have the
        # machine precision of the underlying data type.
        # If there are locks, post an urgent message and return from the method.
        message = ''
        # lockedByMe is from SmartScript
        if self.lockedByMe('QPF', 'SFC'):
            message = '%sYou have the QPF grid locked. Please save the QPF grid.\n' % message
        if self.lockedByMe('PoP', 'SFC'):
            message = '%sYou have the PoP grid locked. Please save the PoP grid.\n' % message
        # lockedByOther is from SmartScript
        if self.lockedByOther('QPF', 'SFC'):
            message = '%sThe QPF grid is locked by someone else. Please have that person save the QPF grid.\n' % message
        if self.lockedByOther('PoP', 'SFC'):
            message = '%sThe PoP grid is locked by someone else. Please have that person save the PoP grid.\n' % message
        if message:
            message = '%sThe QPF/PoP Check was not run.' % message
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have locked grid problems.
            return

        # Make sure there are actually QPF grids in the time range.
        # The self.getGrids command will return None if there are no grids
        # in the time range for mode='First' and noDataError=0. The None
        # variable cannot be iterated over. Rather than trap in a try/except,
        # I'll just check for the condititon. This may not be the most
        # Pythonic way of doing things, but it allows me to avoid having
        # a bunch of code indented beneath a try statement. If no SnowAmt
        # grids are found, post an urgent message and return from the method.
        # getGrids is from SmartScript
        qpfInfoList = self.getGridInfo('Fcst', 'QPF', 'SFC', timeRange)
        if [] == qpfInfoList:
            message = 'There are no QPF grids in the time range you selected.\nThe QPF/PoP Check did not run.'
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have missing grid problems.
            return
        qpfGrids = self.getGrids(
            'Fcst', 'QPF', 'SFC', timeRange, mode='List', noDataError=0,
            cache=0)
        for qpfIndex, qpfGrid in enumerate(qpfGrids):
            gridTR = qpfInfoList[qpfIndex].gridTime()

            popGrid = self.getGrids(
                'Fcst', 'PoP', 'SFC', gridTR, mode='Max', noDataError=0,
                cache=0)
            if popGrid is None:
                message = '''There are no PoP grids in time range %s.
The QPF/PoP Check skipped the time range.''' % gridTR
                self.statusBarMsg(message, 'U')
                continue
            qpfNonZeroMask = greater(qpfGrid, qpfTol)
            popZeroMask = equal(popGrid, 0)
            # popZeroMask = 1 if PoP = 0; popZeroMask = 0 if PoP != 0
            # qpfNonZeroMask = 1 if QPF > 0; qpfNonZeroMask = 0 if QPF = 0
            # PoP = 0 (1) and QPF = 0 (0) => 0 (Good result)
            # PoP != 0 (0) and QPF = 0 (0) => 0 (Good result)
            # PoP != 0 (0) and QPF > 0 (1) => 0 (Good result)
            # PoP = 0 (1) and QPF > 0 (1) => 1 (Bad result)
            consistMask = logical_and(qpfNonZeroMask, popZeroMask)
            consistMask[logical_not(self.cwaMask)] = 0
            if sometrue(ravel(consistMask)):
                # The good result is if the logical_and returns zeros
                # for every grid point, that is "none true". So, if
                # the sometrue method evaluates True, there are
                # inconsistencies.
                self.createGrid(
                    'Fcst', 'QPFPoPInconsistent', 'SCALAR', consistMask, gridTR,
                    descriptiveName='QPFPoPInconsistent',
                    minAllowedValue=0, maxAllowedValue=1, units='Good/Bad')
                if tempGridColor:
                    self.highlightGrids(
                        'Fcst', 'QPFPoPInconsistent', 'SFC', gridTR,
                        tempGridColor)
                if inconGridColor:
                    self.highlightGrids(
                        'Fcst', 'QPF', 'SFC', gridTR, inconGridColor)
                    self.highlightGrids(
                        'Fcst', 'PoP', 'SFC', gridTR, inconGridColor)
                self.inconsistent = True

            ##### Edited by Rob Radzanowski (WFO-CTP) 03-16-2009 to add missing NDFD check for QPF=0 & PoP > 50
            ##### which is causing unexplained yellow banners due to lack of checking for this error.
            qpfZeroMask = equal(qpfGrid, 0)
            popGrid = self.getGrids(
                'Fcst', 'PoP', 'SFC', gridTR, mode='Max', noDataError=0, cache=0)
            popGreater50Mask = greater(popGrid, 50)
            # popGreater50Mask = 1 if PoP > 50; popGreater50Mask = 0 if PoP <= 50
            # qpfZeroMask = 0 if QPF > 0; qpfZeroMask = 1 if QPF = 0
            # PoP > 50 (1) and QPF > 0 (0) => 0 (Good result)
            # PoP > 50 (1) and QPF = 0 (1) => 1 (Bad result)
            # PoP <= 50 (0) and QPF > 0 (0) => 0 (Good/Irrelevant result)
            # PoP <= 50 (0) and QPF = 0 (1) => 0 (Good result)

            consistMask2 = logical_and(qpfZeroMask, popGreater50Mask)
            consistMask2[logical_not(self.cwaMask)] = 0
            if sometrue(ravel(consistMask2)):
                # The good result is if the logical_and returns zeros
                # for every grid point, that is "none true". So, if
                # the sometrue method evaluates True, there are
                # inconsistencies.
                self.createGrid(
                    'Fcst', 'QPFPoPInconsistent', 'SCALAR', consistMask2, gridTR,
                    descriptiveName='QPFPoPInconsistent',
                    minAllowedValue=0, maxAllowedValue=1, units='Good/Bad')

                if tempGridColor:
                    self.highlightGrids('Fcst', 'QPFPoPInconsistent', 'SFC', gridTR, tempGridColor)
                if inconGridColor:
                    self.highlightGrids('Fcst', 'QPF', 'SFC', gridTR, inconGridColor)
                    self.highlightGrids('Fcst', 'PoP', 'SFC', gridTR, inconGridColor)
                self.inconsistent = True
        return

    def _runQPFWxCheck(self, timeRange):
        # This method implements the check that if QPF non zero, then the
        # corresponding Wx grids must contain a precipitable Wx type. Note:
        # the method only checks the Wx type, no cov/prob, no inten, etc.

        # There can be a significant difference between the values stored
        # in memory and the values returned from the database. This is because
        # when values are saved, the element's precision (as defined in
        # serverConfig.py/localConfig.py) is enforced. Values in memory do not
        # have the element's precision enforced; in fact, they have the
        # machine precision of the underlying data type.
        # If there are locks, post an urgent message and return from the method.
        message = ''
        # lockedByMe is from SmartScript
        if self.lockedByMe('QPF', 'SFC'):
            message = '%sYou have the QPF grid locked. Please save the QPF grid.\n' % message
        if self.lockedByMe('Wx', 'SFC'):
            message = '%sYou have the Wx grid locked. Please save the Wx grid.\n' % message
        # lockedByOther is from SmartScript
        if self.lockedByOther('QPF', 'SFC'):
            message = '%sThe QPF grid is locked by someone else. Please have that person save the QPF grid.\n' % message
        if self.lockedByOther('Wx', 'SFC'):
            message = '%sThe Wx grid is locked by someone else. Please have that person save the Wx grid.\n' % message
        if message:
            message = '%sThe QPF/Wx Check was not run.' % message
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have locked grid problems.
            return

        # Make sure there are actually QPF grids in the time range.
        # I'll just check for the condititon. If no SnowAmt
        # grids are found, post an urgent message and return from the method.
        qpfInfoList = self.getGridInfo('Fcst', 'QPF', 'SFC', timeRange)
        if [] == qpfInfoList:
            message = 'There are no QPF grids in the time range you selected.\nThe QPF/PoP Check did not run.'
            self.statusBarMsg(message, 'U')
            # I return instead of aborting because the user may have asked for
            # other tests that do not have missing grid problems.
            return
        for qpfIndex, qpfGrid in enumerate(self.getGrids(
            'Fcst', 'QPF', 'SFC', timeRange, mode='List', noDataError=0,
            cache=0)):
            qpfNonZeroMask = greater(qpfGrid, qpfTol)
            gridTR = qpfInfoList[qpfIndex].gridTime()
            wxInfoList = self.getGridInfo('Fcst', 'Wx', 'SFC', gridTR)
            if [] == wxInfoList:
                message = '''There are no Wx grids in time range %s.
The QPF/Wx Check skipped the time range.''' % gridTR
                self.statusBarMsg(message, 'U')
                continue
            # There are two cases. If the QPF grid is exactly 6 hours long and
            # starts at 00, 06, 12, or 18 UTC, then only one of the
            # corresponding Wx grids needs to be consistent. Otherwise, all the
            # corresponding Wx grids need to be consistent.
            if gridTR.duration() / 3600 == 6 and gridTR.startTime().hour in (0, 6, 12, 18):
                self._qpfWxCheckLocked(qpfNonZeroMask, gridTR, wxInfoList)
            else:
                self._qpfWxCheckUnlocked(qpfNonZeroMask, gridTR, wxInfoList)
        return

    def _qpfWxCheckLocked(self, qpfNonZeroMask, gridTR, wxInfoList):
        # The "Locked" comes from the idea that if the QPF grid is
        # exactly 6 hours long and starts at 00, 06, 12, or 18 UTC, then it
        # is "locked".
        consistMaskList = []
        for wxIndex, wxGrid in enumerate(self.getGrids(
            'Fcst', 'Wx', 'SFC', gridTR, mode='List', noDataError=0,
            cache=0)):
            # wxMask is from SmartScript
            sMask = self.wxMask(wxGrid, ':S:')
            swMask = self.wxMask(wxGrid, ':SW:')
            ipMask = self.wxMask(wxGrid, ':IP:')
            snowMask = logical_or(logical_or(sMask, swMask), ipMask)
            del (sMask, swMask, ipMask)
            rMask = self.wxMask(wxGrid, ':R:')
            rwMask = self.wxMask(wxGrid, ':RW:')
            lMask = self.wxMask(wxGrid, ':L:')
            zlMask = self.wxMask(wxGrid, ':ZL:')
            zrMask = self.wxMask(wxGrid, ':ZR:')
            # logical_or is from Numeric
            rainMask = logical_or(
                rMask, logical_or(
                    rwMask, logical_or(
                        lMask, logical_or(zlMask, zrMask))))
            del (rMask, rwMask, lMask, zlMask, zrMask)
            precipMask = logical_or(snowMask, rainMask)
            del (snowMask, rainMask)
            wxMask = logical_not(precipMask)
            # QPF >= 0.01, 1; QPF < 0.01, 0
            # Wx has precip, 0; Wx doesn't have precip, 1
            # QPF >= 0.01 (1) and Wx has (0) = 0 (Good result)
            # QPF >= 0.01 (1) and Wx doesn't have (1) = 1 (Bad result)
            # QPF < 0.01 (0) and Wx has (0) = 0 (Good result)
            # QPF < 0.01 (0) and Wx doesn't have (1) = 0 (Good result)
            consistMask = logical_and(qpfNonZeroMask, wxMask)
            consistMask[logical_not(self.cwaMask)] = 0
            consistMaskList.append(consistMask)
            if not sometrue(ravel(consistMask)):
                # There were no inconsistencies with this Wx grid. Since only
                # one needs to be consistent, we don't need to do any more
                # checks.
                break
        else:
            # This block will only execute if the for loop runs to
            # completion, i.e., the break statement is not executed.
            # So, if we get here, we have an inconsistency and need to
            # highlight the appropriate grids.
            if inconGridColor:
                self.highlightGrids(
                    'Fcst', 'QPF', 'SFC', gridTR, inconGridColor)
                self.highlightGrids(
                    'Fcst', 'Wx', 'SFC', gridTR, inconGridColor)
            # createGrid is from SmartScript
            for index in xrange(len(wxInfoList)):
                # Create temporary grids for each Wx grid. Limit the time
                # range of the temporary grid so that it doesn't start any
                # earlier or any later than the corresponding QPF grid.
                wxGridTR = wxInfoList[index].gridTime()
                tempGridStartTime = wxGridTR.startTime().unixTime()
                if tempGridStartTime < gridTR.startTime().unixTime():
                    tempGridStartTime = gridTR.startTime().unixTime()
                tempGridEndTime = wxGridTR.endTime().unixTime()
                if tempGridEndTime > gridTR.endTime().unixTime():
                    tempGridEndTime = gridTR.endTime().unixTime()
                tempGridDur = (tempGridEndTime - tempGridStartTime) / 3600
                offset = (tempGridStartTime - \
                          self.timeRange0_1.startTime().unixTime()) / 3600
                # Since the temporary grid could have a different time range
                # than the Wx grid, I need to create and use that time range
                # when creating the temporary grid.
                tempGridTR = self.createTimeRange(
                    offset, offset+tempGridDur, 'Zulu')
                self.createGrid(
                    'Fcst', 'QPFWxInconsistent', 'SCALAR',
                    consistMaskList[index], tempGridTR,
                    descriptiveName='QPFWxInconsistent',
                    minAllowedValue=0, maxAllowedValue=1, units='Good/Bad')
                if tempGridColor:
                    self.highlightGrids(
                        'Fcst', 'QPFWxInconsistent', 'SFC', gridTR,
                        tempGridColor)
            self.inconsistent = True
        return

    def _qpfWxCheckUnlocked(self, qpfNonZeroMask, gridTR, wxInfoList):
        # The "Unlocked" comes from the idea that if the QPF grid is not
        # exactly 6 hours long and starting at 00, 06, 12, or 18 UTC, then it
        # is "unlocked".
        for wxIndex, wxGrid in enumerate(self.getGrids(
            'Fcst', 'Wx', 'SFC', gridTR, mode='List', noDataError=0,
            cache=0)):
            # wxMask is from SmartScript
            sMask = self.wxMask(wxGrid, ':S:')
            swMask = self.wxMask(wxGrid, ':SW:')
            ipMask = self.wxMask(wxGrid, ':IP:')
            snowMask = logical_or(logical_or(sMask, swMask), ipMask)
            del (sMask, swMask, ipMask)
            rMask = self.wxMask(wxGrid, ':R:')
            rwMask = self.wxMask(wxGrid, ':RW:')
            lMask = self.wxMask(wxGrid, ':L:')
            zlMask = self.wxMask(wxGrid, ':ZL:')
            zrMask = self.wxMask(wxGrid, ':ZR:')
            # logical_or is from Numeric
            rainMask = logical_or(
                rMask, logical_or(
                    rwMask, logical_or(
                        lMask, logical_or(zlMask, zrMask))))
            del (rMask, rwMask, lMask, zlMask, zrMask)
            precipMask = logical_or(snowMask, rainMask)
            del (snowMask, rainMask)
            wxMask = logical_not(precipMask)
            # QPF >= 0.01, 1; QPF < 0.01, 0
            # Wx has precip, 0; Wx doesn't have precip, 1
            # QPF >= 0.01 (1) and Wx has (0) = 0 (Good result)
            # QPF >= 0.01 (1) and Wx doesn't have (1) = 1 (Bad result)
            # QPF < 0.01 (0) and Wx has (0) = 0 (Good result)
            # QPF < 0.01 (0) and Wx doesn't have (1) = 0 (Good result)
            #
            # All Wx grids overlapping the SnowAmt grid must be consistent.
            consistMask = logical_and(qpfNonZeroMask, wxMask)
            consistMask[logical_not(self.cwaMask)] = 0
            if sometrue(ravel(consistMask)):
                wxGridTR = wxInfoList[wxIndex].gridTime()
                tempGridStartTime = wxGridTR.startTime().unixTime()
                if tempGridStartTime < gridTR.startTime().unixTime():
                    # Clip to start of QPF grid
                    tempGridStartTime = gridTR.startTime().unixTime()
                tempGridEndTime = wxGridTR.endTime().unixTime()
                if tempGridEndTime > gridTR.endTime().unixTime():
                    # Clip to end of QPF Grid
                    tempGridEndTime = gridTR.endTime().unixTime()
                tempGridDur = (tempGridEndTime - tempGridStartTime) / 3600
                offset = (tempGridStartTime - \
                          self.timeRange0_1.startTime().unixTime()) / 3600
                # Since either the front or end of the Wx grid's
                # time range may have been clipped, create a time
                # range using those values.
                tempGridTR = self.createTimeRange(
                    offset, offset+tempGridDur, 'Zulu')
                self.createGrid(
                    'Fcst', 'QPFWxInconsistent', 'SCALAR', consistMask,
                    tempGridTR, descriptiveName='QPFWxInconsistent',
                    minAllowedValue=0, maxAllowedValue=1, units='Good/Bad')
                if tempGridColor:
                    self.highlightGrids(
                        'Fcst', 'QPFWxInconsistent', 'SFC', gridTR,
                        tempGridColor)
                if inconGridColor:
                    self.highlightGrids(
                        'Fcst', 'QPF', 'SFC', gridTR, inconGridColor)
                    self.highlightGrids(
                        'Fcst', 'Wx', 'SFC', wxGridTR, inconGridColor)
                self.inconsistent = True
        return

    def _calcTolerance(self, gridInfo):
        precision = gridInfo.gridParmInfo.getPrecision()
        return pow(10, -precision)

    def execute(self, timeRange, varDict):
        # Make sure the configuration values are the correct types.
        self.__checkConfigValueTypes()
        # createTimeRange is from SmartScript
        timeRange0_240 = self.createTimeRange(0, 241, 'Zulu')
        checkCleanup = varDict.get('Check_Cleanup', 'Check')
        self.__cleanup(timeRange0_240)
        if checkCleanup == 'Cleanup':
            message = 'SnowQPFPoPWxCheck complete.'
            self.statusBarMsg(message, 'R')
            self.cancel()
        if timeRange.endTime().unixTime() - timeRange.startTime().unixTime() < \
           3600: # No time range selected, use create a 0 to 240 hour range
            timeRange = timeRange0_240

        # If the user has a time range swept out, send an informational
        # message.
        if (timeRange.startTime().unixTime() != timeRange0_240.startTime().unixTime()) or \
           (timeRange.endTime().unixTime() != timeRange0_240.endTime().unixTime()) or \
           (timeRange.duration() != timeRange0_240.duration()):
            message = 'The SnowAmtQPFPoPWxCheck procedure did not run over the 0 to 240 hour time period,\nit ran over %s. This may be what you desired.' % str(timeRange)
            self.statusBarMsg(message, 'S')

        # I'll need to know the unix time of 00Z so I can determine the
        # start time of temporary grids later. I'll need this in more than
        # one of the methods called later, so this will become an instance
        # variable, i.e., prefixed with "self." I also need an instance
        # variable that flags whether or not there were inconsistent grids.
        self.timeRange0_1 = self.createTimeRange(0, 1, 'Zulu')
        self.inconsistent = False

        # A CWA edit area can be provided in the configuration section.
        # Attempt to encode that edit area as a Numeric Python mask so that
        # the later checks are limited to the edit area. The GFE is not very
        # friendly if the encoding fails. The GFE will send a nasty message
        # to the user, but continue executing the procedure. No trappable
        # error is thrown. As of this writing, the GFE appears to create an
        # array of shape (0, 0) if the encoding cannot be done, so I will
        # check for that and, if I find it, then set the edit area to the
        # domain.
        # encodeEditArea comes from SmartScript. For the points that are in
        # the edit area, a value of one is assigned. Otherwise, a value of
        # zero is assigned.
        if cwaEditArea:
            self.cwaMask = self.encodeEditArea(cwaEditArea)
            if self.cwaMask.shape == (0, 0):
                # Use the getGridInfo command to get information about the
                # SnowAmt grid. From this, the grid size can be extracted. I
                # could use getGridInfo on any valid GFE grid.
                # getGridInfo is from SmartScript
                snowAmtInfoList = self.getGridInfo(
                    'Fcst', 'SnowAmt', 'SFC', timeRange)
                # I painfully discovered that the array shape is (y, x)
                gridSize = (snowAmtInfoList[0].gridLocation().gridSize().y,
                            snowAmtInfoList[0].gridLocation().gridSize().x)
                # ones is from Numeric. It creates an array of the given size
                # and data type where all values are one.
                self.cwaMask = ones(gridSize, Int)
                message = \
'''The procedure was not able to use the CWA edit area, %s, provided
in the configuration. You should inform the person responsible for procedures
of this problem. The procedure ran over the whole domain.''' % cwaEditArea
                self.statusBarMsg(message, 'S')
        else:
            snowAmtInfoList = self.getGridInfo(
                'Fcst', 'SnowAmt', 'SFC', timeRange)
            gridSize = (snowAmtInfoList[0].gridLocation().gridSize().y,
                        snowAmtInfoList[0].gridLocation().gridSize().x)
            self.cwaMask = ones(gridSize, Int)

        # Based on the user's input, run the appropriate checks.
        # By making each of these options a checkbox with only one option in
        # the VariableList above, if an option is unchecked then an empty
        # list, [], will be what's in varDict. If an option is checked then a
        # list with the value "Yes", ["Yes"], will be what's in varDict. In
        # Python, a conditional expression can be whether or not a data
        # structure is empty. In these cases, an empty data structure,
        # e.g., an empty list, an empty tuple, an empty dictionary,
        # conditionally test to False while non empty data structures
        # conditionally test to True. In the if statements below, every varDict
        # lookup returns a list: either [] or ["Yes"]. I think the constructs
        # below or more elegant and easier to understand.
        if varDict['Run SnowAmt/QPF Check?']:
            # Call the SnowAmt/QPF check method
            self._runSnowAmtQPFCheck(timeRange)
        if varDict['Run SnowAmt/Wx Check?']:
            # Call the SnowAmt/Wx check method
            self._runSnowAmtWxCheck(timeRange)
        if varDict['Run QPF/PoP Check?']:
            # Call the QPF/PoP check method
            self._runQPFPoPCheck(timeRange)
        if varDict['Run QPF/Wx Check?']:
            # Call the QPF/Wx check method
            self._runQPFWxCheck(timeRange)
        message = 'SnowAmtQPFPoPWxCheck complete.'
        if self.inconsistent:
            message = '%s Inconsistencies found! Grids highlighted %s and %s.' % (
                message, inconGridColor, tempGridColor)
            self.statusBarMsg(message, 'S')
        else:
            self.statusBarMsg(message, 'R')

reason: from numpy


text:
import numpy as np
import os
from maskgen.image_wrap import openImageFile,ImageWrapper
import cv2

"""
Convert a gray image into a single color channel.
"""

def transform(img, source, target, **kwargs):
    channel_map = {
        "red":0,
        "green":1,
        "blue":2
    }
    donor = kwargs['mask'] if 'mask' in kwargs else source
    channel_name = kwargs['channel'] if 'channel' in kwargs else "green"
    img = openImageFile(donor)
    color_im = np.zeros((img.size[1],img.size[0],3),dtype=np.uint8)
    color_im [:,:,channel_map[channel_name]]  = img.image_array
    ImageWrapper(color_im).save(target)
    return None,None

def operation():
    return {'name':'Blend',
          'category':'Layer',
          'description':'Convert a gray image into a single color channel',
          'software':'OpenCV',
          'version':cv2.__version__,
          'arguments': {
              "channel" : {
                  "type":"list",
                  "defaultvalue":"green",
                  "values": ["red","green","blue"],
                  "description": "which channel is set to the gray image"
              },
              "mode": { "type": "text",  "defaultvalue":"Color" }
          },
          'transitions': [
              'image.image'
          ]
          }

def suffix():
    return '.png'

max_stars_repo_name: j-h-m/Media-Journaling-Tool

max_issues_repo_name: j-h-m/Media-Journaling-Tool

max_forks_repo_name: j-h-m/Media-Journaling-Tool


text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn.decomposition

def scale(df, center=True, scale=True):
    ''' compatible with R scale() '''
    if center:
        df = df - df.mean()
    if scale:
        df = df / df.std()
    return df

class PCA:
    def __init__(self, df, npc=2, method='svd'):
        self.df = df
        self.npc = npc
        self.method = method

    def svd(self):
        ''' SVD decomposite X to U*diag(S)*Wt '''
        self.U, self.diagS, self.Wt = np.linalg.svd(self.df, full_matrices=False)

        # singular values are already sorted
        assert np.all(self.diagS[:-1] >= self.diagS[1:])

        ''' T = X*W = U*S '''
        pcXW = np.dot(self.df, self.Wt.T[:,:self.npc])
        pcUdS = self.U[:,:self.npc] * self.diagS[:self.npc]
        # these two ways return same results
        assert np.allclose(pcXW, pcUdS)
        self.score = pcUdS

    def eig(self):
        # calculate eigenvalues(l) and eigenvectors(w) of the covariance matrix
        C = np.cov(self.df.T)
        self.l,self.w = np.linalg.eig(C)

        # sort eigenvectors by eigenvalue in descending order
        self.w = self.w[:,np.argsort(self.l)[::-1]]

        ''' T = X*W '''
        self.score = np.dot(self.df, self.w[:,:self.npc])

    def pc(self):
        if self.method == 'svd':
            self.svd()
        else:
            self.eig()
        self.score[:,1] = -self.score[:,1]
        return self.score

def plot(df):
    for name in df['Species'].unique():
        cond = df['Species'] == name
        plt.plot(df.loc[cond,'PC1'], df.loc[cond,'PC2'], 'o', label=name)
    plt.grid(True)
    plt.legend(framealpha=0.5)

def main():
    np.set_printoptions(precision=4, threshold=30)

    # famous Fisher iris data
    df = pd.read_csv('data/iris.csv')
    print(df)

    # use only numerical columns
    values = scale(df.iloc[:,:4], center=True, scale=True)
    print('DataFrame is a (%d x %d) matrix of iris.csv:\n' % (values.shape[0], values.shape[1]))

    plt.figure(figsize=(12,4))

    plt.subplot(1,3,1)
    plt.title('Singular value decomposition')
    pca_svd = PCA(values, method='svd')
    score = pd.DataFrame(pca_svd.pc(), columns=['PC1', 'PC2'])
    score = pd.concat([df, score], axis=1)
    plot(score)

    plt.subplot(1,3,2)
    plt.title('Eigenvalue decomposition')
    pca_eig = PCA(values, method='eig')
    score = pd.DataFrame(pca_eig.pc(), columns=['PC1', 'PC2'])
    score = pd.concat([df, score], axis=1)
    plot(score)

    plt.subplot(1,3,3)
    plt.title('Principal Component Analysis of scikit-learn')
    pca_sklearn = sklearn.decomposition.PCA(n_components=2)
    score = pca_sklearn.fit_transform(values)
    score = pd.DataFrame(score, columns=['PC1', 'PC2'])
    score = pd.concat([df, score], axis=1)
    plot(score)

    assert np.allclose(score[['PC1', 'PC2']], pca_svd.pc())

    plt.tight_layout()
    plt.savefig('doc/pca.png')

    # assert two princople components based on SVD and EigenDecomposition
    print('PCA(svd) == PCA(eig) is %s' % np.allclose(pca_svd.pc(), pca_eig.pc()))

if __name__ == '__main__':
    main()

max_stars_repo_name: taneishi/pca

max_issues_repo_name: taneishi/pca

max_forks_repo_name: taneishi/pca