text:
# -*- coding: utf-8 -*-

import time
import tensorflow as tf
import numpy as np

class TimeSeriesModel:
    def _verbose_func(avg_cost):
        return np.exp(avg_cost)

    def _optimizer_getter(lr):
        return tf.train.GradientDescentOptimizer(lr)

    def __init__(self, is_training, config, result_map=None, loss_map=None, input_proj=None, use_peepholes=True, state_is_tuple=False, optimizer_getter=_optimizer_getter, verbose_func=_verbose_func):
        self.verbose_func=verbose_func
        self.batch_size=config.batch_size
        self.num_steps=config.input_shape[0]
        self.num_units=config.num_units
        self.last_units=config.last_units

        self.input_data=tf.placeholder(tf.float32, [self.batch_size]+config.input_shape)
        self.targets=tf.placeholder(tf.float32, [self.batch_size]+config.target_shape)
        self.initial_state=tf.placeholder(tf.float32, [self.batch_size, 2*self.num_units*config.num_layers])

        def lstm_cell():
            return tf.nn.rnn_cell.LSTMCell(self.num_units, use_peepholes=use_peepholes, cell_clip=0.95, state_is_tuple=state_is_tuple)
        attn_cell=lstm_cell
        if is_training and config.keep_prob<1:
            def attn_cell():
                return tf.nn.rnn_cell.DropoutWrapper(lstm_cell(), output_keep_prob=config.keep_prob)
        self.cell=tf.nn.rnn_cell.MultiRNNCell([attn_cell() for _ in range(config.num_layers)], state_is_tuple=state_is_tuple)

        if input_proj is None:
            inputs=self.input_data
        else:
            with tf.variable_scope("input_proj", reuse=True):
                inputs=input_proj(self.input_data)

        outputs, state=tf.nn.static_rnn(self.cell, tf.unstack(inputs, num=self.num_steps, axis=1), initial_state=self.initial_state)
        if is_training and loss_map is not None:
            the_inputs=tf.reshape(inputs, [-1]+inputs.get_shape()[2:].as_list())

        self.output=tf.reshape(tf.concat(values=outputs, axis=1), [-1, self.num_units])#[batch_size*self.num_steps, self.num_units]
        self.last_w=tf.get_variable("last_w", [self.num_units, self.last_units])
        self.last_b=tf.get_variable("last_b", [self.last_units])
        self.logits=tf.matmul(self.output, self.last_w)+self.last_b
        self.final_state=state

        if result_map is None:
            self.result=tf.nn.softmax(self.logits)
        else:
            with tf.variable_scope("result_map", reuse=True):
                self.result=result_map(self.logits)
        self.result=tf.reshape(self.result, [self.batch_size, -1, self.last_units])#[batch_size, self.num_steps, self.last_units]

        if is_training:
            if loss_map is None:
                self.loss=tf.contrib.seq2seq.sequence_loss([self.logits], [tf.reshape(self.targets, [-1])], [tf.ones([self.batch_size*self.num_steps])])
            else:
                with tf.variable_scope("loss_map", reuse=True):
                    self.loss=loss_map(self.logits, self.targets, the_inputs)

            self.cost=tf.reduce_sum(self.loss)/self.batch_size

            self.lr=tf.Variable(config.learning_rate, trainable=False)
            tvars=tf.trainable_variables()
            grads, _=tf.clip_by_global_norm(tf.gradients(self.cost, tvars), config.max_grad_norm)
            self.optimizer=optimizer_getter(self.lr)
            self.train_op=self.optimizer.apply_gradients(zip(grads, tvars))

    def assign_lr(self, session, lr_value):
        session.run(tf.assign(self.lr, lr_value))

    def run_epoch(self, session, data_iterator, eval_op, verbose=-1):
        start_time=time.time()
        costs=0.0
        iters=0
        state=self.cell.zero_state(self.batch_size, tf.float32).eval()
        for step, (x, y) in enumerate(data_iterator):
            cost, state, _=session.run([self.cost, self.final_state, eval_op], feed_dict={self.input_data:x, self.targets:y, self.initial_state:state})
            costs+=cost
            iters+=self.num_steps

            if verbose>=0 and step%verbose==0:
                print("step {0} avg perplexity: {1:.3f} current perplexity: {2:.3f} speed: {3:.0f} wps".format(step, self.verbose_func(costs/iters), self.verbose_func(cost), iters*self.batch_size/(time.time()-start_time)))

        return self.verbose_func(costs/(iters if iters!=0 else 1))

    def run_training(self, session, data_iterator, verbose=-1):
        start_time=time.time()
        costs=0.0
        iters=0
        for step, (x, y, state) in enumerate(data_iterator):
            cost, state, _=session.run([self.cost, self.final_state, self.train_op], feed_dict={self.input_data:x, self.targets:y, self.initial_state:state})
            costs+=cost
            iters+=self.num_steps

            if verbose>=0 and step%verbose==0:
                print("step {0} avg perplexity: {1:.3f} current perplexity: {2:.3f} speed: {3:.0f} sps".format(step, self.verbose_func(costs/iters), self.verbose_func(cost/self.num_steps), iters*self.batch_size/(time.time()-start_time)))

        return self.verbose_func(costs/(iters if iters!=0 else 1))

    def run_predicting(self, session, data_iterator):
        start_time=time.time()
        for step, (input, state) in enumerate(data_iterator):
            if state is None:
                state=self.initial_state.eval()
            result, state=session.run([self.result, self.final_state], feed_dict={self.input_data:input, self.initial_state:state})
            yield result, state
            print("step {0} finished. time: {1:.3f} sec".format(step, time.time()-start_time))

max_stars_repo_name: xiaoshenxian/mlxsx

max_issues_repo_name: xiaoshenxian/mlxsx

max_forks_repo_name: xiaoshenxian/mlxsx


text:
#!/usr/bin/python -W ignore
from PIL import Image
import numpy as np
import os
from sklearn import linear_model

trainingCancerDir = "../images/extract-automatic/cancer/"
trainingBloodDir = "../images/extract-automatic/blood/"
trainingOtherDir = "../images/extract-automatic/other/"
testCancerDir = "../images/extract-automatic/cancer-test/"
testBloodDir = "../images/extract-automatic/blood-test/"
testOtherDir = "../images/extract-automatic/other-test/"

def getY(numCancerSamples, numBloodSamples, numOtherSamples):
    '''Returns a numpy array of Y values'''
    cancerArray = np.ones(numCancerSamples)
    bloodArray = np.zeros(numBloodSamples)
    otherArray = np.full(numOtherSamples, 2)
    return np.concatenate((cancerArray,bloodArray,otherArray))

def readImage(path):
    # look at the jpegs
    if path.endswith(".jpg"):

        # open each image (max size is known to be 30x30)
        image = Image.open(path)

        # grab the pixel values into a 2d numpy array
        pixels = np.array(image)

        # reshape the 2d array into a 1d vector of RGB values
        r = pixels[:,:,0].reshape(-1)
        g = pixels[:,:,1].reshape(-1)
        b = pixels[:,:,2].reshape(-1)

        if len(r) < 900:
            rgbfill = np.full(900-len(r), -1)
            r = np.append(r, rgbfill)
            g = np.append(g, rgbfill)
            b = np.append(b, rgbfill)

        return np.concatenate((r,g,b))

    return np.array([0,0])

def getTrainingData():
    trainingImageList = []

    numCancerSamples = 0
    numBloodSamples = 0
    numOtherSamples = 0

    # read the training data
    for dirname in (trainingCancerDir, trainingBloodDir, trainingOtherDir):
        for filename in os.listdir(dirname):
            p = readImage(dirname+filename)
            if p.size > 0:
                # add the vector to our list of training data
                trainingImageList.append(p)

                # store the number of positive samples we've got ready for later
                if dirname == trainingCancerDir:
                    numCancerSamples += 1

                # store the number of positive samples we've got ready for later
                if dirname == trainingBloodDir:
                    numBloodSamples += 1

                # store the number of positive samples we've got ready for later
                if dirname == trainingOtherDir:
                    numOtherSamples += 1

    # convert training data into numpy array
    X = np.array(trainingImageList)

    # generate labels
    y = getY(numCancerSamples,numBloodSamples,numOtherSamples)

    return (X, y)

def runTestData():
    testCancerExamples = []
    testBloodExamples = []
    testOtherExamples = []

    for filename in os.listdir(testCancerDir):
        p = readImage(testCancerDir+filename)
        if p.size > 0:
            testCancerExamples.append({"image":p,"filename":filename});

    for filename in os.listdir(testBloodDir):
        p = readImage(testBloodDir+filename)
        if p.size > 0:
            testBloodExamples.append({"image":p,"filename":filename})

    for filename in os.listdir(testOtherDir):
        p = readImage(testOtherDir+filename)
        if p.size > 0:
            testOtherExamples.append({"image":p,"filename":filename})

    print "Test cancer cells (should be 1.0):"
    correct = 0
    wrong = 0
    for testCancerExample in testCancerExamples:
        prediction = logreg.predict(testCancerExample["image"])[0]
        if prediction == 1.0:
            correct += 1
        else:
            wrong += 1
        print "  " + testCancerExample["filename"] + " " + str(prediction)
    print correct, "correct,", wrong, "wrong,", "Accuracy =", str(correct*100.0/(correct+wrong))

    print "Test blood cells (should be 0.0):"
    correct = 0
    wrong = 0
    for testBloodExample in testBloodExamples:
        prediction = logreg.predict(testBloodExample["image"])[0]
        if prediction == 0.0:
            correct += 1
        else:
            wrong += 1
        print "  " + testBloodExample["filename"] + " " + str(prediction)
    print correct, "correct,", wrong, "wrong,", "Accuracy =", str(correct*100.0/(correct+wrong))

    print "Test other cells (should be 2.0):"
    correct = 0
    wrong = 0
    for testOtherExample in testOtherExamples:
        prediction = logreg.predict(testOtherExample["image"])[0]
        if prediction == 2.0:
            correct += 1
        else:
            wrong += 1
        print "  " + testOtherExample["filename"] + " " + str(prediction)
    print correct, "correct,", wrong, "wrong,", "Accuracy =", str(correct*100.0/(correct+wrong))

if __name__ == "__main__":

    X, y = getTrainingData()

    print "X = " + str(X.shape[0]) + " x " + str(X.shape[1]) + " with " + str(len(y)) + " labels"

    # train the model
    logreg = linear_model.LogisticRegression()
    logreg.fit(X, y)

    print "Training complete"

    runTestData()

ext: py

max_stars_repo_name: mangoeatscicada/Tackling-Cancer

max_issues_repo_name: mangoeatscicada/Tackling-Cancer

max_forks_repo_name: mangoeatscicada/Tackling-Cancer

reason: import numpy


text:
import os
import numpy as np
import math
import abc

from gym import spaces
import pybullet as p
from gym_pybullet_drones.envs.single_agent_rl.obstacles.LandingZone import LandingZone
from gym_pybullet_drones.envs.single_agent_rl import BaseSingleAgentAviary

POSITIVE_REWARD = 1
NEGATIVE_REWARD = -1

class DenseReward:
    """Dense reward for a drone.

    The reward is the distance to the closest obstacle.

    """

    def __init__(self, aviary, scale):
        self.aviary = aviary
        self.scale = scale

    def _calculateReward():
        return 0

    def _getDroneStateVector(self, drone_id=0):
        """Get the drone state vector.

        Parameters
        ----------
        drone_id : int
            Drone id.

        Returns
        -------
        np.array
            Drone state vector.

        """
        # state = p.getLinkState(drone_id, 0)
        state = self.aviary._getDroneStateVector(drone_id)
        return np.array(state)

    def calculateReward():
        return self._calculateReward() * self.scale

class BoundsReward(DenseReward):
    """Checks rewards for out of bounds"""

    def __init__(self, aviary, scale, bounds, useTimeScaling=False):
        super().__init__(aviary, scale)
        self.landing_zone = landing_zone
        # Defined as [[x_high, y_high, z_high], [x_low, y_low, z_low]]
        self.bounds = bounds
        self.XYZ_IDX = [0, 1, 2]
        self.useTimeScaling = useTimeScaling

    def _calculateReward(self):
        state = self._getDroneStateVector(0)
        position = state[0:3]

        for dim_idx in self.XYZ_IDX:
            if (
                state[dim_idx] > self.bounds[0, dim_idx]
                or state[dim_idx] < self.bounds[1, dim_idx]
            ):
                self.aviary.completeEpisode = True
                if self.useTimeScaling:
                    return NEGATIVE_REWARD * (
                        1
                        - (
                            (self.aviary.step_counter / self.aviary.SIM_FREQ)
                            / self.aviary.EPISODE_LEN_SEC
                        )
                    )
                return NEGATIVE_REWARD

class LandingReward(DenseReward):
    """Calculate the landing reward."""

    def __init__(self, aviary, scale, landing_zone: LandingZone):
        super().__init__(aviary, scale)
        self.landing_zone = landing_zone

    def _calculateReward(self):
        state = self._getDroneStateVector(0)
        position = state[0:3]
        velocity = state[10:13]
        target_position = self.landing_zone.xyz
        target_velocity = np.asarray([0, 0, 0])
        pos_dist = np.linalg.norm(position - target_position)
        vel_dist = np.linalg.norm(velocity - target_velocity)

        if pos_dist < 0.1 and vel_dist < 0.1:
            self.aviary.landing_frames += 1
            self.aviary.completeEpisode = True
            if self.aviary.landing_frames >= 10:
                return 2240
            else:
                return 80
        elif pos_dist < 1:
            if vel_dist < 1:
                inv_vel_dist = 2 - vel_dist
            else:
                inv_vel_dist = 1 / vel_dist

            return 2 - inv_vel_dist + 0.85 * inv_vel_dist
        # Penalize if out of bounds
        return 1 / pos_dist

    ################################################################################

max_stars_repo_name: trangml/gym-pybullet-drones

max_issues_repo_name: trangml/gym-pybullet-drones

max_forks_repo_name: trangml/gym-pybullet-drones


text:
#!/usr/bin/python
"""
PlotFundamentalCouplings.py

Plot an illustration of the fundamental (Lagrangian-level)
couplings which correspond to the regions where good
discrimination is possible.

BJK 15/09/2017
"""

import numpy as np

#----- Matplotlib parameters ------

import matplotlib.pyplot as pl
import matplotlib as mpl
import matplotlib.colors as colors

font = {'family' : 'sans-serif',
        'size'   : 18}
mpl.rc('font', **font)

mpl.rcParams['xtick.major.size'] = 8
mpl.rcParams['xtick.major.width'] = 1
mpl.rcParams['xtick.minor.size'] = 3
mpl.rcParams['xtick.minor.width'] = 1
mpl.rcParams['ytick.major.size'] = 8
mpl.rcParams['ytick.major.width'] = 1
mpl.rcParams['ytick.minor.size'] = 3
mpl.rcParams['ytick.minor.width'] = 1
mpl.rcParams['axes.linewidth'] = 1.5

mpl.rc('font', **font)
mpl.rc('text.latex', preamble=[r'\usepackage{color}', r'\usepackage{amssymb}'])

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

#---Functions---

#Calculate fundamental Lagrangian couplings (see Eq. 1)
#for a given value of f and rnp = lambda_n/lambda_p
#Note: lne denotes lambda_{n,e} etc. (coupling of
#even operator to neutrons).
#Returns couplings in format (lpe, lpo, lne, lno)

#Positive Solution
def CalculateFundamental_plus(rnp, f, lne, lno = 1.0):

    theta = np.arccos(f)

    lpe = (np.cos(theta)*lne - np.sin(theta)*lno)*1.0/rnp
    lpo = (np.sin(theta)*lne + np.cos(theta)*lno)*1.0/rnp

    return (lpe, lpo, lne, lno)

#Negative Solution
def CalculateFundamental_minus(rnp, f, lne, lno = 1.0):

    theta = -np.arccos(f)

    lpe = (np.cos(theta)*lne - np.sin(theta)*lno)*1.0/rnp
    lpo = (np.sin(theta)*lne + np.cos(theta)*lno)*1.0/rnp

    return (lpe, lpo, lne, lno)

#----Calculations----

# Generate 100 samples each with a different value of lne
# (lno is fixed to 1)
# These 100 samples will correspond to 200 fundamental couplings
# because there is a positive and negative solution for each.
Np = 100
points = np.zeros((2*Np, 4))
lne_vals = np.linspace(-1.0, 1.0, Np)

# Specify central values (and spread) for r_np and f
rnp0 = 0.75
f0 = -0.99

delta_rnp = 0.05
delta_f = 0.005

#Calculate fundamental couplings (for values of r_np and f near the central values)
for i in range(Np):
    points[2*i,:] = CalculateFundamental_plus(rnp0+delta_rnp*(2.0*np.random.rand()-1.0), f0+delta_f*(2.0*np.random.rand()-1.0), lne_vals[i])
    points[2*i+1,:] = CalculateFundamental_minus(rnp0+delta_rnp*(2.0*np.random.rand()-1.0), f0+delta_f*(2.0*np.random.rand()-1.0), lne_vals[i])

#----Plotting----

fig, ax = pl.subplots(1,1, figsize=(8,6))

scat = ax.scatter(points[:,0], points[:,2], c=points[:,1],vmin=-1.5, vmax=-1.0)

# Add diagonal line for central value of r_np
ax.plot([-2.0, 2.0], [2.0*0.75, -2.0*0.75], 'k--', linewidth=2.0)

# Some grid lines
ax.axvline(0, linestyle='-', color='grey', alpha=0.5)
ax.axhline(0, linestyle='-', color='grey', alpha=0.5)
ax.set_xlim(-1.75,1.75)
ax.set_ylim(-1.75,1.75)

# Add colorbar
cbar = fig.colorbar(scat,ax=ax)
cbar.ax.set_ylabel("$\lambda_{p,o}$", fontsize=20)

# Add labels
ax.set_xlabel("$\lambda_{p,e}$", fontsize=20)
ax.set_ylabel("$\lambda_{n,e}$", fontsize=20)
ax.text(-1.55,-1.0, r"$f \sim - 0.99$", ha="left")
ax.text(-1.55,-1.25, r"$\lambda_n/\lambda_p \sim 0.75$", ha="left")
ax.text(-1.55,-1.5, r"$\lambda_{n,o} = 1$", ha="left")

# Output to file
pl.tight_layout()
pl.savefig("../plots/FundamentalCouplings.pdf", bbox_inches="tight")
#pl.show()

max_stars_repo_name: LBJ-Wade/AntiparticleDM

max_issues_repo_name: LBJ-Wade/AntiparticleDM

max_forks_repo_name: LBJ-Wade/AntiparticleDM


text:
r"""Solve Stokes equations using a coupled formulation

The Stokes equations are in strong form

.. math::

    \nabla^2 u - \nabla p &= f \\
    \nabla \cdot u &= h \\
    u(x, y, z=\pm 1) &= 0

where :math:`f` and :math:`g` are given functions of space.
In addition we require :math:`\int p d\ = 0`, which is achieved by
fixing the coefficient :math:`\hat{p}_{0, 0, 0} = 0`.

We use a tensorproductspace with Fourier expansions in the x- and
y-directions and a composite Chebyshev or Legendre basis in the z-direction
for ``u`` and a regular (no boundary conditions) Chebyshev or Legendre basis
for ``p``.

When both Fourier wavenumbers are zero the assembled coefficient matrix has
two nullspaces. One of these are removed by enforcing the global constraint
on the pressure. The second is removed by fixing :math:`\hat{p}_{0, 0, N-1} = 0`.

"""
import os
import sys
import numpy as np
from sympy import symbols, sin, cos
from shenfun import *

x, y, z = symbols("x,y,z", real=True)

# Some right hand side (manufactured solution)
uex = sin(2*y)*(1-z**2)
uey = sin(2*x)*(1-z**2)
uez = sin(2*z)*(1-z**2)
pe = -0.1*sin(2*x)*cos(4*y)
fx = uex.diff(x, 2) + uex.diff(y, 2) + uex.diff(z, 2) - pe.diff(x, 1)
fy = uey.diff(x, 2) + uey.diff(y, 2) + uey.diff(z, 2) - pe.diff(y, 1)
fz = uez.diff(x, 2) + uez.diff(y, 2) + uez.diff(z, 2) - pe.diff(z, 1)
h = uex.diff(x, 1) + uey.diff(y, 1) + uez.diff(z, 1)

N = (20, 20, 20)
family = sys.argv[-1] if len(sys.argv) == 2 else 'Legendre'
K0 = FunctionSpace(N[0], 'Fourier', dtype='D', domain=(0, 2*np.pi))
K1 = FunctionSpace(N[1], 'Fourier', dtype='d', domain=(0, 2*np.pi))
SD = FunctionSpace(N[2], family, bc=(0, 0))
ST = FunctionSpace(N[2], family)

TD = TensorProductSpace(comm, (K0, K1, SD), axes=(2, 0, 1))
Q = TensorProductSpace(comm, (K0, K1, ST), axes=(2, 0, 1))
V = VectorSpace(TD)
VQ = CompositeSpace([V, Q])

up = TrialFunction(VQ)
vq = TestFunction(VQ)

u, p = up
v, q = vq

# Assemble blocks of complete matrix
if family.lower() == 'chebyshev':
    A = inner(v, div(grad(u)))
    G = inner(v, -grad(p))
else:
    A = inner(grad(v), -grad(u))
    G = inner(div(v), p)
D = inner(q, div(u))

# Create block matrix
M = BlockMatrix(A+G+D)

# Get f and h on quad points
fh = Array(VQ, buffer=(fx, fy, fz, h))
f_, h_ = fh

fh_hat = Function(VQ)
f_hat, h_hat = fh_hat
f_hat = inner(v, f_, output_array=f_hat)
h_hat = inner(q, h_, output_array=h_hat)

# Solve problem using integral constraint on pressure
up_hat = M.solve(fh_hat, constraints=((3, 0, 0), (3, N[2]-1, 0)))
up = up_hat.backward()
u_, p_ = up

# Exact solution
ux, uy, uz = Array(V, buffer=(uex, uey, uez))
pe = Array(Q, buffer=pe)

error = [comm.reduce(np.linalg.norm(ux-u_[0])),
         comm.reduce(np.linalg.norm(uy-u_[1])),
         comm.reduce(np.linalg.norm(uz-u_[2])),
         comm.reduce(np.linalg.norm(pe-p_))]

if comm.Get_rank() == 0:
    print('Error    u          v          w        p')
    print('     %2.4e %2.4e %2.4e %2.4e' %(error[0], error[1], error[2], error[3]))
    assert np.all(abs(np.array(error)) < 1e-8), error

if 'pytest' not in os.environ:
    import matplotlib.pyplot as plt
    plt.figure()
    X = Q.local_mesh(True)
    plt.contourf(X[0][:, :, 0], X[1][:, :, 0], u_[2, :, :, 6], 100)
    plt.figure()
    plt.quiver(X[0][:, :, 0], X[1][:, :, 0], u_[0, :, :, 6], u_[1, :, :, 6])
    plt.figure()
    l, m = 2, 2
    plt.spy(M.diags((l, m), 'csr'), markersize=2, color='k') # The matrix for Fourier given wavenumber
    plt.title('Block matrix: l, m = ({}, {})'.format(l, m))
    plt.xticks([])
    plt.yticks([])
    plt.figure()
    plt.contourf(X[0][:, :, 0], X[1][:, :, 0], u_[0, :, :, 6], 100)
    plt.show()

max_stars_repo_name: spectralDNS/shenfun

max_issues_repo_name: spectralDNS/shenfun

max_forks_repo_name: spectralDNS/shenfun

reason: import numpy,from sympy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jun  8 10:07:02 2021

Instead of using the uncertainity for each region, we will use the total uncertainty

@author: ccamargo
"""

import os
import xarray as xr
import numpy as np
import warnings
import sys
sys.path.append("/Users/ccamargo/Documents/py_scripts/")
import utils_SL as sl

import utils_SLE as sle
import pandas as pd

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.gridspec import GridSpec
import cmocean as cm
from cartopy import crs as ccrs, feature as cfeature

from scipy.integrate import simps,trapz#,cumtrapz,romb

#%%# % % get ANT dataset
# IMBIE 2018 Antartic Dataset
#This spreadsheet contains the IMBIE-2018 datasets for Antarctica (worksheet 1),
#Antarctic Peninsula (worksheet 2), East Antarctica (worksheet 3) and West Antarctica (worksheet 4).
#Each worksheet includes data on monthly cumulative ice sheet mass changes and their estimated uncertainty.
#The data are expressed in units of mass (Gigatons – columns B and C)
#and in units of equivalent mean global sea level rise (millimetres – columns D and E).

df=pd.read_excel('/Volumes/LaCie_NIOZ/data/barystatic/original/ANT_IMBIE.xlsx',
                  #sheet_name=[0,1] # load first and second sheet as a dict of df
                  sheet_name=0 # sheet 1
                  )

df.columns
i,j=df.shape
data=np.zeros((i,j))
c=[None]*(j+1)
for j,col in enumerate(df.columns):
    data[:,j]=np.array(df[col])
    c[j]=col

ant=data[:,4] # uncertainity
years=np.array(df.Year)
tstring=sle.from_tdec_to_tstring(years)

#% % Now get data for each region:
sheets=[1,2,3]
icesheets=['AP','EA','WA']
ant_is=np.zeros((len(sheets),len(ant)))
glb=[0,0,0]
for l,n in enumerate(sheets):
    print(n)
    df=pd.read_excel('/Volumes/LaCie_NIOZ/data/barystatic/original/ANT_IMBIE.xlsx',
                      #sheet_name=[0,1] # load first and second sheet as a dict of df
                      sheet_name=n # sheet 1
                      )
    df.columns
    i,j=df.shape
    data=np.zeros((i,j))
    c=[None]*(j+1)
    for j,col in enumerate(df.columns):
        data[:,j]=np.array(df[col])
        c[j]=col

    ant_is[l,:]=data[:,4] # sea level contr uncertianity (mm)

#%%  make regional time series
ice_in=np.array(ant_is)
path_to_mask='/Volumes/LaCie_NIOZ/data/barystatic/masks/barystatic_mask2.nc'
# path_to_mask='/Volumes/LaCie_NIOZ/data/barystatic/masks/rignot_basins/ANT_Basins_IMBIE2_v1.6/ANT_Basins_IMBIE2_v1.6.nc'
# path_to_mask='/Volumes/LaCie_NIOZ/PhD/Barystatic/antarctica/rignot2019/basins/AIS_basins_v2.nc'
path_to_save='/Volumes/LaCie_NIOZ/data/barystatic/intrinsic_unc/',

#% %
#% % open ANT mask:
ds=xr.open_dataset(path_to_mask)
ds
# on IMBIE we don't have island separation, so we will use 'mask_noisland':
mask=np.array(ds.mask10) # using zwally_IceSheets_1x1_etopo_cor.nc
print(ds.code_mask10)
# we want codes 10-12
mask[np.where(mask>12)]='nan'
mask[np.where(mask<10)]='nan'
np.unique(mask)
grid=np.zeros((mask.shape))
#% %
# mask=np.array(ds.mask_noisland)
mask=np.hstack(mask)
lat=np.array(ds.lat)
lon=np.array(ds.lon)
llon,llat=np.meshgrid(lon,lat)
#llat,llon=np.meshgrid(lat,lon) THIS IS WRONG!!!!!!!!!!!!!!!
llat=np.hstack(llat)
llon=np.hstack(llon)

df=pd.DataFrame(llon,columns=['lon'])
df['lat']=llat
df['mask']=mask

#%% get grid area
grid_area=sl.get_grid_area(grid)
df['area']=np.hstack(grid_area)
area_per_grid=np.hstack(grid_area)
#%% test for one time step
k=1
ice_slc=np.zeros((len(llat)))

# AIS
mean_load=np.array(ant[k])
region_area=df['area'].where(np.isfinite(df['mask'])).sum()
load_per_area=np.array(mean_load/region_area)
ice_slc[np.where(np.isfinite(mask))]=load_per_area*area_per_grid[np.where(np.isfinite(mask))]

print(np.sum(ant[k]))
print(np.sum(ice_slc))

#%% loop over time
# use mean value of AIS
# codes=[2,4,3]
data_reg=np.zeros((len(years),len(ds.lat),len(ds.lon)))

# Loop over time:
for k in range(0,len(years)):
    print(k)
    #Divide global number by total area of each region to obtain the mean number:
    ice_slc=np.zeros((len(llat)))

    # AIS
    mean_load=np.array(ant[k])
    region_area=df['area'].where(np.isfinite(df['mask'])).sum()
    load_per_area=np.array(mean_load/region_area)
    ice_slc[np.where(np.isfinite(mask))]=load_per_area*area_per_grid[np.where(np.isfinite(mask))]

    data_reg[k,:,:]=np.array(ice_slc).reshape((len(ds.lat),len(ds.lon)))
df['ice_slc']=ice_slc

#%% make dataset

da=xr.Dataset(data_vars={'mask':(('lat','lon'),mask.reshape(180,360)),
                          'ant_slc_reg':(('time','lat','lon'),-data_reg),
                          'ant_glb':(('time'),-ant)

                        },
                          coords={'lat':ds.lat,
                                  'lon':ds.lon,

                                  'time':years,
                                  'reg':['AIS']
                                  })

#% %
da.attrs['mask_exp']='Mask for Antarctica Basisns: 0 = land + ocean; 3 = Ant Peninsula; 1 = West Ant; 2 = East Ant'
da.attrs['Metadata0']='ant_glb: global values of cumlative sea-level change for entire Antarctica (sheet 0 of IMBIE dataset)'
da.attrs['Metadata1']='ant_slc_reg: Sea-level change UNCERTAINTY from Regional Cumulative sea level contribution (mm), done by redistributing the global values of cumulative slc'
da.attrs['source'] = 'IMBIE dataset 2018 (imbie_dataset-2018_07_23).'

da.attrs['script']='intrinsic_unc-IMBIE-prep_v2'

# da.to_netcdf('/Users/ccamargo/Documents/PhD/Barystatic/imbie/imbie_ant_reg.nc')
da.to_netcdf(path_to_save[0]+'AIS_IMB.nc')

#%% GRE
#%% Greenland

# # IMBIE 2019 Greendland Dataset
# #This spreadsheet contains the IMBIE-2019 datasets for Greenland,
# #which includes data on the annual rate of change and cumulative change in Greenland’s ice sheet mass,
# #its surface mass balance and ice discharge anomalies, and their estimated uncertainty.

# # Sheet 2: equivalent mean global sea level rise (mm/yr) and sea-level rise (mm)
# # and in units of equivalent mean global sea level rise (millimetres per year – sheet 2, columns B, C, F, G, J and K,
# #and millimetres – sheet 2, columns D, E, H, I, L and M).
df=pd.read_excel('/Volumes/LaCie_NIOZ/data/barystatic/original/GRE_IMBIE.xlsx',
                  #sheet_name=[0,1] # load first and second sheet as a dict of df
                  sheet_name=1 # sheet 2
                  )
i,j=df.shape
data=np.zeros((i,j))
c=[None]*(j+1)
for j,col in enumerate(df.columns):
    data[:,j]=np.array(df[col])
    c[j]=col

years=np.array(df.Year)
print(years[144])

# From 1992 onwards UNCERTAINTIES

gre_tot=data[144:len(years),4] #Cumulative ice sheet mass change (mm sea level)

years=years[144:len(years)]
#% %
plt.plot(years,gre_tot,'green',label='total')

#%% Now using the Greenland Mask
#open regional grid mask:
# ds=xr.open_dataset('/Volumes/LaCie_NIOZ/PhD/Barystatic/greenland/mouginout2019/GrIS_division/GIS_basins.nc')
# ds=xr.open_dataset('/Volumes/LaCie_NIOZ/data/barystatic/masks/Mouginot_GrIS_division/Greenland_Basins_PS_v1.4.2.shp.nc')
ds=xr.open_dataset(path_to_mask)
ds
# 0:NW, 1:CE, 2:CW, 3:SW, 4:NO, 5:NE, 6:SE
ds.mask10.plot()
maskgre=np.array(ds.mask10)
# we want from 2 to 8
maskgre[np.where(maskgre>8)]='nan'
maskgre[np.where(maskgre<2)]='nan'
# remove 2 so we have it from 0-6
maskgre=maskgre-2
np.unique(maskgre[np.isfinite(maskgre)])

lat=np.array(ds.lat)
lon=np.array(ds.lon)
llon,llat=np.meshgrid(lon,lat)
#llat,llon=np.meshgrid(lat,lon) THIS IS WRONG!!!!!!!!!!!!!!!
llat=np.hstack(llat)
llon=np.hstack(llon)

df=pd.DataFrame(llon,columns=['lon'])
df['lat']=llat
df['mask']=np.hstack(maskgre)

df['area']=np.hstack(grid_area)

area_per_grid=np.hstack(grid_area)
#%% GRE TOT
data_reg = np.zeros((len(gre_tot),len(lat),len(lon)))

for k in range(0,len(years)):
    print(k)

    # GIS
    ice_slc=np.zeros((len(llat)))
    mean_load=np.array(gre_tot[k])
    #Divide global number by total area of each region to obtain the mean number:
    region_area=df['area'].where(np.isfinite(df['mask'])).sum()
    load_per_area=np.array(mean_load/region_area)
    ice_slc[np.isfinite(df['mask'])]=load_per_area*area_per_grid[np.isfinite(df['mask'])]

    #% %
    data_reg[k,:,:]=np.array(ice_slc).reshape((len(ds.lat),len(ds.lon)))

#%%

#% %
da=xr.Dataset(data_vars={'gre_slc_tot':(('time','lat','lon'),-data_reg),
                         'gre_glb':(('time'),-gre_tot),
                         'mask':(('lat','lon'),maskgre),
#                         'maskdyn':(('lat','lon'),maskdyn)  ,

                         },
                         coords={'lat':lat,
                                 'lon':lon,
                                 'time':years,
                                 'reg':['TOT']
                                 })

#%%
da.attrs['maskgre_exp']='Greenland mask: 0:NW, 1:CE, 2:CW, 3:SW, 4:NO, 5:NE, 6:SE Greenland Drainage division from Mouginot and Rignot, 2019'
da.attrs['Metadata']='gre_glb: global values of cumlative sea-level change UNCERTAINTY for entire ice sheet (column E of excel)'
da.attrs['units']='mm'
da.attrs['regional']='Sea-level change UNC from Regional Cumulative sea level contribution (mm), done by redistributing the global values of cumulative slc'

da.attrs['source'] = 'IMBIE dataset 2019 (imbie_dataset_greenland_dynamics-2020_02_28).'
da.attrs['script']='intrinsic_unc-IMBIE-prep_v2'
# da.to_netcdf('/Users/ccamargo/Documents/PhD/Barystatic/imbie/imbie_gre_reg.nc')
# da.to_netcdf('/Volumes/LaCie_NIOZ/PhD/Barystatic/imbie/imbie_gre_reg.nc')
da.to_netcdf(path_to_save[0]+'GIS_IMB.nc')

#%% Standardize
#%%
import datetime as dt
import utils_SLE_v2 as sle
#%%
#%% Description
# from the barystatic_test_units.py we saw that:
# datasets from IMBIE, Rignot 2019, MOUGINOT 2019 are in mm of SL

flist=sl.get_filelist(path_to_save[0],ext='*.nc')
flist=[f for f in flist if f.split('.')[-2].split('_')[-1]=='IMB']

print(flist)

names=['AIS_IMB',
       'GIS_IMB']
variables = ['ant_slc_reg','gre_slc_tot']
#%%
# time range: months 1993-2020.08
time,time2=sl.make_date('01-01-1993',(27*12)+8)
tdec,tdec0=sl.get_dec_time(time2,ns=1e-6)
t=[dt.datetime.utcfromtimestamp(t.astype(int) * 1e-6) for t in time2]
ty=np.array([dt.datetime.utcfromtimestamp(t.astype(int) * 1e-6).timetuple().tm_year
                                for t in time2])
tm=np.array([dt.datetime.utcfromtimestamp(t.astype(int) * 1e-6).timetuple().tm_mon
                                for t in time2])

# yearly
timey=np.arange(1993,2021)
idx=np.zeros((len(timey)))
for i,year in enumerate(timey):
    idx[i]=np.where(ty==year)[0][0]

SL_mm=np.zeros((len(names),len(time),180,360))
SL_mm.fill('nan')
SL_EWH=np.full_like(SL_mm,np.nan)

SL_mm_y=np.zeros((len(names),len(timey),180,360))
SL_mm_y.fill('nan')
SL_EWH_y=np.full_like(SL_mm_y,np.nan)

for iname,f in enumerate(flist):
    ds=xr.open_dataset(f)
    tdec_local=np.array(ds.time)
    time_local,time2_local=sl.make_date('01-01-'+str(int(tdec_local[0])),len(tdec_local))
    ds['time']=time_local
    ds_y=ds.groupby('time.year').mean(dim='time')

    name=variables[iname]
    lon=np.array(ds.lon)
    lat=np.array(ds.lat)
    timey_local=np.array(ds_y.year)

    # find index for yearly:
    idx_local=np.zeros((len(timey)))
    idx_local.fill('nan')
    # idx_local=np.full_like(timey,np.nan)
    for i,year in enumerate(timey):
        if year<=np.max(timey_local) and \
            year>=np.min(timey_local):
            idx_local[i]=np.where(timey_local==year)[0][0]

    idx_0y=np.array(np.where(np.isfinite(idx_local)))[0,0]
    idx_00y=int(idx_local[np.isfinite(idx_local)][0])
    idx_1y=np.array(np.where(np.isfinite(idx_local)))[0,-1]+1

    # find index for monthly:
    ty_local=np.array([dt.datetime.utcfromtimestamp(t.astype(int) * 1e-6).timetuple().tm_year
                                    for t in time2_local])
    tm_local=np.array([dt.datetime.utcfromtimestamp(t.astype(int) * 1e-6).timetuple().tm_mon
                                    for t in time2_local])

    idx_local_m=np.zeros((len(ty)))
    idx_local_m.fill('nan')
    t2_local=[dt.datetime.utcfromtimestamp(t.astype(int) * 1e-6) for t in time2_local]
    for i,iy in enumerate(t):
        year=iy.year
        if np.any(ty_local==year):
            month=iy.month
            if np.any(tm_local==month):
                # print(year)
                # print(month)
                jdx=np.array(np.where(year==ty_local))[0]
                jdx2=np.array(np.where(month==tm_local))[0]
                for jd in jdx2:
                    if np.any(jdx==jd):
                        j=jdx[jdx==jd]
                        idx_local_m[i]=j

    # Get data
    data=np.array(ds[name]) # mmm of SL height
    datay=np.array(ds_y[name])
    # Put data in main array
    SL_mm_y[iname,idx_0y:idx_1y,:,:]=datay[idx_00y:len(timey_local),:,:]
    for i,j in enumerate(idx_local_m):
        if np.isfinite(j):
            j=int(j)
            SL_mm[iname,i,:,:]=data[j,:,:]
    # Transform data:
    for i in range(len(data)):
        data[i,:,:]=sle.height_to_EWH(data[i,:,:]).reshape(180,360)
    for i in range(len(datay)):
        datay[i,:,:]=sle.height_to_EWH(datay[i,:,:]).reshape(180,360)
    # Put transformed data in main array:
    SL_EWH_y[iname,idx_0y:idx_1y,:,:]=datay[idx_00y:len(timey_local),:,:]
    for i,j in enumerate(idx_local_m):
        if np.isfinite(j):
            j=int(j)
            SL_EWH[iname,i,:,:]=data[j,:,:]

#%%
da=xr.Dataset(data_vars={'SL_mm':(('name','time','lat','lon'),SL_mm),
                              'SL_EWH':(('name','time','lat','lon'),SL_EWH),
                              'SL_mm_y':(('name','year','lat','lon'),SL_mm_y),
                              'SL_EWH_y':(('name','year','lat','lon'),SL_EWH_y),
                              'mask_gre':(('lat','lon'),maskgre),
                              'mask_ant':(('lat','lon'),mask.reshape(180,360)),

                            },
                              coords={'lat':ds.lat,
                                      'lon':ds.lon,
                                      'time':time2,
                                      'tdec':tdec,
                                      'year':timey,
                                      'name':names})
da['SL_mm'].attrs['units']='mm of sea level'
da['SL_mm_y'].attrs['units']='mm of sea level'
da['SL_EWH'].attrs['units']='mm of Equivalent Water Thickness'
da['SL_EWH_y'].attrs['units']='mm of Equivalent Water Thickness'
da['SL_mm'].attrs['long_name']='Monthly ocean mass change in mm of sea level height'
da['SL_EWH'].attrs['long_name']='Monthly ocean mass change in mm of equivalent water thickness'
da['SL_mm_y'].attrs['long_name']='Yearly averages of ocean mass change in mm of sea level height'
da['SL_EWH_y'].attrs['long_name']='Yearly averages of ocean mass change in mm of equivalent water thickness'

da['name'].attrs['long_name']='Name of the source region/contribution of ocean mass change'
#% %
da.attrs['metadata']='Ocean mass sea level changes UNCERTANTIES (in mm and EWH), from IMBIE ANT and GRE'
da.attrs['sources']='IMBIE Team, 2020; 2018, Nature '

da.attrs['script']='intrinsic_unc-IMBIE-prep.py'
da.attrs['Comment']='different sources of barystatic contributions standardized to use as input for the Sl Equation'
da.attrs['Author']='Carolina M.L. Camargo'
da.attrs['reference']='Camargo et al. 2021'
da.attrs['date_created']=str(dt.datetime.now())

path='/Volumes/LaCie_NIOZ/data/barystatic/intrinsic_unc/use/'
da.to_netcdf(path+'IMBIE_v2.nc')

max_stars_repo_name: carocamargo/barystaticSLC

max_issues_repo_name: carocamargo/barystaticSLC

max_forks_repo_name: carocamargo/barystaticSLC

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
R"""
Created on Fri Dec 18 22:59:59 2020

@author: Christian

In this example, the node positions are optimized by .
Three nodes are fixed in the example, the top most node, and the two bottom
nodes.
The connectivity between nodes is also fixed.

"""

"""
Import functions and set the random seed for reproducability.
"""
import naturalize as nat
import numpy as np
import matplotlib.pyplot as plt
from functions import ftest, fitness_basic, fitness_Complex, environment, plotIndividual

dataName = 'Analysis data.P'

xlims = np.array([0,2])
ylims = np.array([0,2])
trussMat = 10       # can be 10 or 11
forces = np.array([10000., 1000., 0.])
env = environment(xlims, ylims, trussMat, forces)

"""
Set up the gene pool. Two genes will be used, one for the x position of each
node, and one for the y position of each node.
"""

llims = [np.array([0,0,0,0]), np.array([0,0,0,0])]
ulims = [np.array([1, 1, 1, 1]),np.array([2, 2, 2, 2])]
pool = nat.BasicGenePool(llims, ulims)

"""
Define population size. A fairly large sample size is used. We use a mutation
rate of about 1/N_gene.
"""
helper = nat.AlgorithmHelper(ftest, pool, fitness_Complex, environment = env)

Ngen = 1000
Npop = 25
Ncouples = 12
Nsurvive = 0
mutateThresold = 0.065

"""
Define the algorithm, analysis, and a recoder.
"""
algorithm = nat.GeneticAlgorithm(helper, Npop, Ncouples, Nsurvive, mutateThresold)
recorder = nat.basicRecorder(5, 5)
analysis = nat.Analysis(algorithm, recorder)
best = analysis.runAnalysis(Ngen)

"""
Save the analyisis to a pickle and plot the output data. It's often useful to
see how the results of the analysis are changing over time.'
"""
data = analysis.getRecorderData()
nat.pickleData(data, dataName)
plt.plot(data.genNumber, data.bestScores)
plotIndividual(data.bestIndividuals[-1], env)

max_stars_repo_path: examples/Structural Engineering Examples/Ex2 - OpenSees Optimization of a Truss/Ex2.2 - Truss Node Positions/main.py

max_stars_repo_name: cslotboom/Naturalize

max_issues_repo_name: cslotboom/Naturalize

max_forks_repo_name: cslotboom/Naturalize


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

"""
递归神经网络测试
"""

import numpy as np
from NeuralNetwork.RNN import tree_rnn
from NeuralNetwork.ANN.activators import IdentityActivator

def gradient_check():
    """
    梯度检查
    """
    # 设计一个误差函数，取所有节点输出项之和
    error_function = lambda o: o.sum()

    rnn = tree_rnn.RecursiveLayer(2, 2, IdentityActivator(), 1e-3)

    # 计算forward值
    x, d = tree_rnn.data_set()
    rnn.forward(x[0], x[1])
    rnn.forward(rnn.root, x[2])

    # 求取sensitivity map
    sensitivity_array = np.ones((rnn.node_width, 1),
                                dtype=np.float64)
    # 计算梯度
    rnn.backward(sensitivity_array)

    # 检查梯度
    epsilon = 10e-4
    for i in range(rnn.W.shape[0]):
        for j in range(rnn.W.shape[1]):
            rnn.W[i, j] += epsilon
            rnn.reset_state()
            rnn.forward(x[0], x[1])
            rnn.forward(rnn.root, x[2])
            err1 = error_function(rnn.root.data)
            rnn.W[i, j] -= 2 * epsilon
            rnn.reset_state()
            rnn.forward(x[0], x[1])
            rnn.forward(rnn.root, x[2])
            err2 = error_function(rnn.root.data)
            expect_grad = (err1 - err2) / (2 * epsilon)
            rnn.W[i, j] += epsilon
            print('weights(%d,%d): expected - actural %.4e - %.4e' % (
                i, j, expect_grad, rnn.W_grad[i, j]))
    return rnn

if __name__ == '__main__':
    gradient_check()

max_stars_repo_name: neohope/MLStudy

max_issues_repo_name: neohope/MLStudy

max_forks_repo_name: neohope/MLStudy


text:
from functools import partial
from unittest import TestCase
import os
import os.path as osp

import numpy as np

from datumaro.components.extractor import (
    AnnotationType, Bbox, DatasetItem, Label, LabelCategories, Points, Polygon,
    PolyLine,
)
from datumaro.components.project import Dataset
from datumaro.plugins.cvat_format.converter import CvatConverter
from datumaro.plugins.cvat_format.extractor import CvatImporter
from datumaro.util.image import Image
from datumaro.util.test_utils import (
    TestDir, compare_datasets, test_save_and_load,
)

from .requirements import Requirements, mark_requirement

DUMMY_IMAGE_DATASET_DIR = osp.join(osp.dirname(__file__),
    'assets', 'cvat_dataset', 'for_images')

DUMMY_VIDEO_DATASET_DIR = osp.join(osp.dirname(__file__),
    'assets', 'cvat_dataset', 'for_video')

class CvatImporterTest(TestCase):
    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_detect_image(self):
        self.assertTrue(CvatImporter.detect(DUMMY_IMAGE_DATASET_DIR))

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_detect_video(self):
        self.assertTrue(CvatImporter.detect(DUMMY_VIDEO_DATASET_DIR))

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_load_image(self):
        expected_dataset = Dataset.from_iterable([
            DatasetItem(id='img0', subset='train',
                image=np.ones((8, 8, 3)),
                annotations=[
                    Bbox(0, 2, 4, 2, label=0, z_order=1,
                        attributes={
                            'occluded': True,
                            'a1': True, 'a2': 'v3'
                        }),
                    PolyLine([1, 2, 3, 4, 5, 6, 7, 8],
                        attributes={'occluded': False}),
                ], attributes={'frame': 0}),
            DatasetItem(id='img1', subset='train',
                image=np.ones((10, 10, 3)),
                annotations=[
                    Polygon([1, 2, 3, 4, 6, 5], z_order=1,
                        attributes={'occluded': False}),
                    Points([1, 2, 3, 4, 5, 6], label=1, z_order=2,
                        attributes={'occluded': False}),
                ], attributes={'frame': 1}),
        ], categories={
            AnnotationType.label: LabelCategories.from_iterable([
                ['label1', '', {'a1', 'a2'}],
                ['label2'],
            ])
        })

        parsed_dataset = Dataset.import_from(DUMMY_IMAGE_DATASET_DIR, 'cvat')

        compare_datasets(self, expected_dataset, parsed_dataset)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_load_video(self):
        expected_dataset = Dataset.from_iterable([
            DatasetItem(id='frame_000010', subset='annotations',
                image=255 * np.ones((20, 25, 3)),
                annotations=[
                    Bbox(3, 4, 7, 1, label=2,
                        id=0,
                        attributes={
                            'occluded': True,
                            'outside': False, 'keyframe': True,
                            'track_id': 0
                        }),
                    Points([21.95, 8.00, 2.55, 15.09, 2.23, 3.16],
                        label=0,
                        id=1,
                        attributes={
                            'occluded': False,
                            'outside': False, 'keyframe': True,
                            'track_id': 1, 'hgl': 'hgkf',
                        }),
                ], attributes={'frame': 10}),
            DatasetItem(id='frame_000013', subset='annotations',
                image=255 * np.ones((20, 25, 3)),
                annotations=[
                    Bbox(7, 6, 7, 2, label=2,
                        id=0,
                        attributes={
                            'occluded': False,
                            'outside': True, 'keyframe': True,
                            'track_id': 0
                        }),
                    Points([21.95, 8.00, 9.55, 15.09, 5.23, 1.16],
                        label=0,
                        id=1,
                        attributes={
                            'occluded': False,
                            'outside': True, 'keyframe': True,
                            'track_id': 1, 'hgl': 'jk',
                        }),
                    PolyLine([7.85, 13.88, 3.50, 6.67, 15.90, 2.00, 13.31, 7.21],
                        label=2,
                        id=2,
                        attributes={
                            'occluded': False,
                            'outside': False, 'keyframe': True,
                            'track_id': 2,
                        }),
                ], attributes={'frame': 13}),
            DatasetItem(id='frame_000016', subset='annotations',
                image=Image(path='frame_0000016.png', size=(20, 25)),
                annotations=[
                    Bbox(8, 7, 6, 10, label=2,
                        id=0,
                        attributes={
                            'occluded': False,
                            'outside': True, 'keyframe': True,
                            'track_id': 0
                        }),
                    PolyLine([7.85, 13.88, 3.50, 6.67, 15.90, 2.00, 13.31, 7.21],
                        label=2,
                        id=2,
                        attributes={
                            'occluded': False,
                            'outside': True, 'keyframe': True,
                            'track_id': 2,
                        }),
                ], attributes={'frame': 16}),
        ], categories={
            AnnotationType.label: LabelCategories.from_iterable([
                ['klhg', '', {'hgl'}],
                ['z U k'],
                ['II']
            ]),
        })

        parsed_dataset = Dataset.import_from(DUMMY_VIDEO_DATASET_DIR, 'cvat')

        compare_datasets(self, expected_dataset, parsed_dataset)

class CvatConverterTest(TestCase):
    def _test_save_and_load(self, source_dataset, converter, test_dir,
            target_dataset=None, importer_args=None, **kwargs):
        return test_save_and_load(self, source_dataset, converter, test_dir,
            importer='cvat',
            target_dataset=target_dataset, importer_args=importer_args, **kwargs)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_save_and_load(self):
        src_label_cat = LabelCategories(attributes={'occluded', 'common'})
        for i in range(10):
            src_label_cat.add(str(i))
        src_label_cat.items[2].attributes.update(['a1', 'a2', 'empty'])

        source_dataset = Dataset.from_iterable([
            DatasetItem(id=0, subset='s1', image=np.zeros((5, 10, 3)),
                annotations=[
                    Polygon([0, 0, 4, 0, 4, 4],
                        label=1, group=4,
                        attributes={ 'occluded': True, 'common': 't' }),
                    Points([1, 1, 3, 2, 2, 3],
                        label=2,
                        attributes={ 'a1': 'x', 'a2': 42, 'empty': '',
                            'unknown': 'bar' }),
                    Label(1),
                    Label(2, attributes={ 'a1': 'y', 'a2': 44 }),
                ]
            ),
            DatasetItem(id=1, subset='s1',
                annotations=[
                    PolyLine([0, 0, 4, 0, 4, 4],
                        label=3, id=4, group=4),
                    Bbox(5, 0, 1, 9,
                        label=3, id=4, group=4),
                ]
            ),

            DatasetItem(id=2, subset='s2', image=np.ones((5, 10, 3)),
                annotations=[
                    Polygon([0, 0, 4, 0, 4, 4], z_order=1,
                        label=3, group=4,
                        attributes={ 'occluded': False }),
                    PolyLine([5, 0, 9, 0, 5, 5]), # will be skipped as no label
                ]
            ),

            DatasetItem(id=3, subset='s3', image=Image(
                path='3.jpg', size=(2, 4))),
        ], categories={ AnnotationType.label: src_label_cat })

        target_label_cat = LabelCategories(
            attributes={'occluded'}) # unable to represent a common attribute
        for i in range(10):
            target_label_cat.add(str(i), attributes={'common'})
        target_label_cat.items[2].attributes.update(['a1', 'a2', 'empty', 'common'])
        target_dataset = Dataset.from_iterable([
            DatasetItem(id=0, subset='s1', image=np.zeros((5, 10, 3)),
                annotations=[
                    Polygon([0, 0, 4, 0, 4, 4],
                        label=1, group=4,
                        attributes={ 'occluded': True, 'common': 't' }),
                    Points([1, 1, 3, 2, 2, 3],
                        label=2,
                        attributes={ 'occluded': False, 'empty': '',
                            'a1': 'x', 'a2': 42 }),
                    Label(1),
                    Label(2, attributes={ 'a1': 'y', 'a2': 44 }),
                ], attributes={'frame': 0}
            ),
            DatasetItem(id=1, subset='s1',
                annotations=[
                    PolyLine([0, 0, 4, 0, 4, 4],
                        label=3, group=4,
                        attributes={ 'occluded': False }),
                    Bbox(5, 0, 1, 9,
                        label=3, group=4,
                        attributes={ 'occluded': False }),
                ], attributes={'frame': 1}
            ),

            DatasetItem(id=2, subset='s2', image=np.ones((5, 10, 3)),
                annotations=[
                    Polygon([0, 0, 4, 0, 4, 4], z_order=1,
                        label=3, group=4,
                        attributes={ 'occluded': False }),
                ], attributes={'frame': 0}
            ),

            DatasetItem(id=3, subset='s3', image=Image(
                    path='3.jpg', size=(2, 4)),
                attributes={'frame': 0}),
        ], categories={ AnnotationType.label: target_label_cat })

        with TestDir() as test_dir:
            self._test_save_and_load(source_dataset,
                partial(CvatConverter.convert, save_images=True), test_dir,
                target_dataset=target_dataset)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_allow_undeclared_attrs(self):
        source_dataset = Dataset.from_iterable([
            DatasetItem(id=0, annotations=[
                Label(0, attributes={ 'x': 4, 'y': 2 }),
                Bbox(1, 2, 3, 4, label=0, attributes={ 'x': 1, 'y': 1 }),
            ]),
        ], categories=[ ('a', '', {'x'}) ])

        target_label_cat = LabelCategories(attributes={'occluded'})
        target_label_cat.add('a', attributes={'x'})
        target_dataset = Dataset.from_iterable([
            DatasetItem(id=0, annotations=[
                Label(0, attributes={ 'x': 4, 'y': 2 }),
                Bbox(1, 2, 3, 4, label=0,
                    attributes={ 'x': 1, 'y': 1, 'occluded': False }),
            ], attributes={'frame': 0}),
        ], categories={ AnnotationType.label: target_label_cat })

        with TestDir() as test_dir:
            self._test_save_and_load(source_dataset,
                partial(CvatConverter.convert, allow_undeclared_attrs=True),
                test_dir, target_dataset=target_dataset)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_relative_paths(self):
        source_dataset = Dataset.from_iterable([
            DatasetItem(id='1', image=np.ones((4, 2, 3))),
            DatasetItem(id='subdir1/1', image=np.ones((2, 6, 3))),
            DatasetItem(id='subdir2/1', image=np.ones((5, 4, 3))),
        ])

        target_dataset = Dataset.from_iterable([
            DatasetItem(id='1', image=np.ones((4, 2, 3)),
                attributes={'frame': 0}),
            DatasetItem(id='subdir1/1', image=np.ones((2, 6, 3)),
                attributes={'frame': 1}),
            DatasetItem(id='subdir2/1', image=np.ones((5, 4, 3)),
                attributes={'frame': 2}),
        ], categories=[])

        with TestDir() as test_dir:
            self._test_save_and_load(source_dataset,
                partial(CvatConverter.convert, save_images=True), test_dir,
                target_dataset=target_dataset, require_images=True)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_save_dataset_with_cyrillic_and_spaces_in_filename(self):
        label_categories = LabelCategories(attributes={'occluded'})
        for i in range(10):
            label_categories.add(str(i))
        label_categories.items[2].attributes.update(['a1', 'a2', 'empty'])

        source_dataset = Dataset.from_iterable([
            DatasetItem(id='кириллица с пробелом',
                subset='s1', image=np.zeros((5, 10, 3)),
                annotations=[
                    Label(1),
                ]
            ),
        ], categories={
            AnnotationType.label: label_categories,
        })

        target_dataset = Dataset.from_iterable([
            DatasetItem(id='кириллица с пробелом',
                subset='s1', image=np.zeros((5, 10, 3)),
                annotations=[
                    Label(1),
                ], attributes={'frame': 0}
            ),
        ], categories={
            AnnotationType.label: label_categories,
        })

        with TestDir() as test_dir:
            self._test_save_and_load(source_dataset,
                partial(CvatConverter.convert, save_images=True), test_dir,
                target_dataset=target_dataset, require_images=True)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_can_save_and_load_image_with_arbitrary_extension(self):
        expected = Dataset.from_iterable([
            DatasetItem('q/1', image=Image(path='q/1.JPEG',
                data=np.zeros((4, 3, 3))), attributes={'frame': 1}),
            DatasetItem('a/b/c/2', image=Image(path='a/b/c/2.bmp',
                data=np.zeros((3, 4, 3))), attributes={'frame': 2}),
        ], categories=[])

        with TestDir() as test_dir:
            self._test_save_and_load(expected,
                partial(CvatConverter.convert, save_images=True),
                test_dir, require_images=True)
            self.assertTrue(osp.isfile(
                osp.join(test_dir, 'images', 'q', '1.JPEG')))
            self.assertTrue(osp.isfile(
                osp.join(test_dir, 'images', 'a', 'b', 'c', '2.bmp')))

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_preserve_frame_ids(self):
        expected_dataset = Dataset.from_iterable([
            DatasetItem(id='some/name1', image=np.ones((4, 2, 3)),
                attributes={'frame': 40}),
        ], categories=[])

        with TestDir() as test_dir:
            self._test_save_and_load(expected_dataset,
                CvatConverter.convert, test_dir)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_reindex(self):
        source_dataset = Dataset.from_iterable([
            DatasetItem(id='some/name1', image=np.ones((4, 2, 3)),
                attributes={'frame': 40}),
        ])

        expected_dataset = Dataset.from_iterable([
            DatasetItem(id='some/name1', image=np.ones((4, 2, 3)),
                attributes={'frame': 0}),
        ], categories=[])

        with TestDir() as test_dir:
            self._test_save_and_load(source_dataset,
                partial(CvatConverter.convert, reindex=True), test_dir,
                target_dataset=expected_dataset)

    @mark_requirement(Requirements.DATUM_GENERAL_REQ)
    def test_inplace_save_writes_only_updated_data(self):
        with TestDir() as path:
            # generate initial dataset
            dataset = Dataset.from_iterable([
                DatasetItem(1, subset='a'),
                DatasetItem(2, subset='b'),
                DatasetItem(3, subset='c', image=np.ones((2, 2, 3))),
            ])
            dataset.export(path, 'cvat', save_images=True)
            os.unlink(osp.join(path, 'a.xml'))
            os.unlink(osp.join(path, 'b.xml'))
            os.unlink(osp.join(path, 'c.xml'))
            self.assertFalse(osp.isfile(osp.join(path, 'images', '2.jpg')))
            self.assertTrue(osp.isfile(osp.join(path, 'images', '3.jpg')))

            dataset.put(DatasetItem(2, subset='a', image=np.ones((3, 2, 3))))
            dataset.remove(3, 'c')
            dataset.save(save_images=True)

            self.assertTrue(osp.isfile(osp.join(path, 'a.xml')))
            self.assertFalse(osp.isfile(osp.join(path, 'b.xml')))
            self.assertTrue(osp.isfile(osp.join(path, 'c.xml')))
            self.assertTrue(osp.isfile(osp.join(path, 'images', '2.jpg')))
            self.assertFalse(osp.isfile(osp.join(path, 'images', '3.jpg')))

max_stars_repo_name: certiware/posemaro

max_issues_repo_name: certiware/posemaro

max_forks_repo_name: certiware/posemaro


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

import sys, os, logging
# logging.disable(sys.maxsize)

## Set dir
testdir = os.path.dirname(os.path.abspath(__file__))
datdir = testdir+'/dat/'
outdir = testdir+'/out/'
if not os.path.exists(outdir):
    os.makedirs(outdir)

import numpy as np
from astropy.io import fits
from astropy.wcs import WCS

## Local
sys.path.insert(0, testdir+'/..') ## astylo path
from astylo.iolib import read_fits, write_fits
from astylo.astrolib import fixwcs
from astylo.ipro import (improve, islice, icrop, iswarp, iconvolve,
                          sextract, wclean, interfill, hextract, hswarp)

print('\n TEST improve ')
print('--------------')
imp = improve(datdir+'M83', verbose=True)
print('raw: ', imp.im[0,0,:])
imp.rand_norm(datdir+'M83_unc')
print('add sym rand: ', imp.im[0,0,:])
imp.rand_splitnorm([datdir+'M83_unc', datdir+'M83_unc'])
print('add split rand', imp.im[0,0,:])
imp.crop(outdir+'M83_imp_crop', sizval=(0.005,0.009), cenpix=(7,8))
print('See out/M83_imp_crop.fits [Done]')

print('\n TEST islice ')
print('-------------')
slc = islice(datdir+'M83', filSL=outdir+'M83_inv_sqrt',
             filUNC=datdir+'M83_unc', slicetype='inv_sq', postfix='_test_')
print(slc.filenames()[0])
if input('Clean slices (y/n): ')=='y':
    slc.clean()

print('\n TEST icrop ')
print('------------')
crp = icrop(datdir+'M83', filOUT=outdir+'M83_icrop',
            sizpix=(3,6), cenval=(204.2529675, -29.8656962),
            filUNC=[datdir+'M83_unc', datdir+'M83_unc'],
            dist='splitnorm', verbose=True)
print('See out/M83_icrop.fits [Done]')

print('\n TEST iswarp ')
print('-------------')
hdr_ref = fixwcs(datdir+'M82_09_L86').header
swp = iswarp(refheader=hdr_ref,
             center=None, pixscale=6.,
             tmpdir=outdir+'swp/')
swp.combine(datdir+'M82_08_L86', uncpdf='norm',
            filOUT=outdir+'M82_08_L86_rep', tmpdir=outdir+'swp/')
print('Reproject M82_08_L86 to M82_09_L86 (pixscale=6") [Done]')

swp_cube = iswarp((datdir+'M82_09_L86', datdir+'M82_04_SL1'),
                  center='9:55:51,69:40:45', pixscale=None,
                  tmpdir=outdir+'swp_cube/')
swp_cube.combine(datdir+'M82_04_SL1', uncpdf='norm',
                 filOUT=outdir+'M82_04_SL1_rep')
print('Reproject M82_04_SL1 to M82_09_L86 (recenter to M82 center) [Done]')

print('\n TEST iconvolve ')
print('----------------')
## See also idlib/conv_prog.pro & idlib/convolve_image.pro
convdir = outdir+'conv/' # see also idlib/convolve_image.pro
if not os.path.exists(convdir):
    os.makedirs(convdir)
path_ker = datdir
path_idl = testdir+'/../idlib/'
csv_ker = outdir+'kernelist' # see also idlib/conv_prog.pro

irs_ker = []
psf = [2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6.]
psf_ref = 'Gauss_06.0'
for p in psf:
    irs_ker.append(path_ker+'Kernel_HiRes_Gauss_0'+str(p)+'_to_'+psf_ref)
conv_cube = iconvolve(datdir+'M82_04_SL1',
                      kfile=irs_ker, klist=csv_ker,
                      convdir=convdir, filOUT=outdir+'M82_04_SL1'+'_conv')
conv_cube.do_conv(path_idl, verbose=False)
print('Convolve M82_04_SL1 [Done]')

irac_ker = path_ker+'Kernel_HiRes_IRAC_8.0_to_Gauss_06.0'
conv = iconvolve(datdir+'M82_IRAC4',
                 kfile=irac_ker, klist=csv_ker,
                 filOUT=outdir+'M82_IRAC4'+'_conv')
conv.do_conv(path_idl)
print('Convolve M82_IRAC4 [Done]')

print('\n TEST sextract ')
print('---------------')

print('\n TEST wclean ')
print('-------------')

print('\n TEST interfill ')
print('----------------')
data = read_fits(datdir+'IC10_SL2').data[0]
hdr = fixwcs(datdir+'IC10_SL2').header
newdata = interfill(data, axis=0)
write_fits(outdir+'IC10_fillgap', hdr, newdata)
print('interfill test [Done]')

print('\n TEST hextract ')
print('---------------')
hextract(datdir+'M82_08_L86', outdir+'M82_hextract', 20, 40, 32, 42)
print('hextract test [Done]')

print('\n TEST hswarp ')
print('-------------')
print('Reproject M82_08_L86 to M82_09_L86')
old = read_fits(datdir+'M82_08_L86')
oldimage = old.data
oldheader = old.header
refheader = read_fits(datdir+'M82_09_L86').header
hswp = hswarp(oldimage, oldheader, refheader, keepedge=True,
              tmpdir=outdir+'hswp/', verbose=True)
# print('hswarp image: ', hswp.image)
# print('hswarp image header: ', hswp.header)
print('See hswp/coadd.fits [Done]')

if input('Clean tmp files (y/n): ')=='y':
    swp.clean()
    swp_cube.clean()
    conv_cube.clean()
    conv.clean()
    conv.clean(outdir+'hswp/')

max_stars_repo_name: kxxdhdn/astylo

max_issues_repo_name: kxxdhdn/astylo

max_forks_repo_name: kxxdhdn/astylo

reason: import numpy,from astropy


text:
# -*- coding: utf-8 -*-
"""
Created on Sun July 20 22:19:53 2021
©2020. Triad National Security, LLC. All rights reserved. This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National Nuclear Security Administration. All rights in the program are reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear Security Administration. The Government is granted for itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare derivative works, distribute copies to the public, perform publicly and display publicly, and to permit others to do so.
@author: Reshad
"""
import os
import random
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cv2
from tqdm import tqdm_notebook, tnrange
from itertools import chain
from skimage.io import imread, imshow, concatenate_images
from skimage.transform import resize
from skimage.morphology import label
from sklearn.model_selection import train_test_split

import tensorflow as tf

from keras.models import Model, load_model
from keras.layers import Input, BatchNormalization, Activation, Dense, Dropout
from keras.layers.core import Lambda, RepeatVector, Reshape
from keras.layers.convolutional import Conv2D, Conv2DTranspose
from keras.layers.pooling import MaxPooling2D, GlobalMaxPool2D
from keras.layers.merge import concatenate, add
from keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
from keras.optimizers import Adam
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img

# Set some parameters
im_width = 128
im_height = 128
border = 5
path_train = r'D:\job\LAL\data\Figures_Only-20210704T024044Z-001\deep learning\training/'
#path_train = r'D:\job\LAL\data\data\400_figures\deep-learning/'

# Get and resize train images and masks
def get_data(path, train=True):
    ids = next(os.walk(path + "images"))[2]
    X = np.zeros((len(ids), im_height, im_width, 1), dtype=np.float32)
    if train:
        y = np.zeros((len(ids), im_height, im_width, 1), dtype=np.float32)
    print('Getting and resizing images ... ')
    for n, id_ in tqdm_notebook(enumerate(ids), total=len(ids)):
        # Load images
        img = load_img(path + '/images/' + id_, grayscale=True)
        x_img = img_to_array(img)
        x_img = resize(x_img, (128, 128, 1), mode='constant', preserve_range=True)

        # Load masks
        if train:
            mask = img_to_array(load_img(path + '/masks/' + id_, grayscale=True))
            mask = resize(mask, (128, 128, 1), mode='constant', preserve_range=True)

        # Save images
        X[n, ..., 0] = x_img.squeeze() / 255
        if train:
            y[n] = mask / 255
    print('Done!')
    if train:
        return X, y
    else:
        return X

X_train, y_train = get_data(path_train, train=True)

# Check if training data looks all right
ix = random.randint(0, len(X_train))
has_mask = y_train[ix].max() > 0

fig, ax = plt.subplots(1, 2, figsize=(20, 10))

ax[0].imshow(X_train[ix, ..., 0], cmap='seismic', interpolation='bilinear')
if has_mask:
    ax[0].contour(y_train[ix].squeeze(), colors='k', levels=[0.5])
ax[0].set_title('Input Images (dummy color)')

ax[1].imshow(y_train[ix].squeeze(), interpolation='bilinear', cmap='gray')
ax[1].set_title('Ground-truth');

################## Unet #########################3

def conv2d_block(input_tensor, n_filters, kernel_size=3, batchnorm=True):
    # first layer
    x = Conv2D(filters=n_filters, kernel_size=(kernel_size, kernel_size), kernel_initializer="he_normal",
               padding="same")(input_tensor)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation("relu")(x)
    # second layer
    x = Conv2D(filters=n_filters, kernel_size=(kernel_size, kernel_size), kernel_initializer="he_normal",
               padding="same")(x)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation("relu")(x)
    return x

def get_unet(input_img, n_filters=16, dropout=0.5, batchnorm=True):
    # contracting path
    c1 = conv2d_block(input_img, n_filters=n_filters * 1, kernel_size=3, batchnorm=batchnorm)
    p1 = MaxPooling2D((2, 2))(c1)
    p1 = Dropout(dropout * 0.5)(p1)

    c2 = conv2d_block(p1, n_filters=n_filters * 2, kernel_size=3, batchnorm=batchnorm)
    p2 = MaxPooling2D((2, 2))(c2)
    p2 = Dropout(dropout)(p2)

    c3 = conv2d_block(p2, n_filters=n_filters * 4, kernel_size=3, batchnorm=batchnorm)
    p3 = MaxPooling2D((2, 2))(c3)
    p3 = Dropout(dropout)(p3)

    c4 = conv2d_block(p3, n_filters=n_filters * 8, kernel_size=3, batchnorm=batchnorm)
    p4 = MaxPooling2D(pool_size=(2, 2))(c4)
    p4 = Dropout(dropout)(p4)

    c5 = conv2d_block(p4, n_filters=n_filters * 16, kernel_size=3, batchnorm=batchnorm)

    # expansive path
    u6 = Conv2DTranspose(n_filters * 8, (3, 3), strides=(2, 2), padding='same')(c5)
    u6 = concatenate([u6, c4])
    u6 = Dropout(dropout)(u6)
    c6 = conv2d_block(u6, n_filters=n_filters * 8, kernel_size=3, batchnorm=batchnorm)

    u7 = Conv2DTranspose(n_filters * 4, (3, 3), strides=(2, 2), padding='same')(c6)
    u7 = concatenate([u7, c3])
    u7 = Dropout(dropout)(u7)
    c7 = conv2d_block(u7, n_filters=n_filters * 4, kernel_size=3, batchnorm=batchnorm)

    u8 = Conv2DTranspose(n_filters * 2, (3, 3), strides=(2, 2), padding='same')(c7)
    u8 = concatenate([u8, c2])
    u8 = Dropout(dropout)(u8)
    c8 = conv2d_block(u8, n_filters=n_filters * 2, kernel_size=3, batchnorm=batchnorm)

    u9 = Conv2DTranspose(n_filters * 1, (3, 3), strides=(2, 2), padding='same')(c8)
    u9 = concatenate([u9, c1], axis=3)
    u9 = Dropout(dropout)(u9)
    c9 = conv2d_block(u9, n_filters=n_filters * 1, kernel_size=3, batchnorm=batchnorm)

    outputs = Conv2D(1, (1, 1), activation='sigmoid')(c9)
    model = Model(inputs=[input_img], outputs=[outputs])
    return model

input_img = Input((im_height, im_width, 1), name='img')
model = get_unet(input_img, n_filters=16, dropout=0., batchnorm=True)

model.compile(optimizer=Adam(), loss="binary_crossentropy", metrics=["accuracy"])

data_gen_args = dict(horizontal_flip=True,
                     vertical_flip=True)
image_datagen = ImageDataGenerator(**data_gen_args)
mask_datagen = ImageDataGenerator(**data_gen_args)

seed = 2018
bs = 32

image_generator = image_datagen.flow(X_train, seed=seed, batch_size=bs, shuffle=True)
mask_generator = mask_datagen.flow(y_train, seed=seed, batch_size=bs, shuffle=True)

# Just zip the two generators to get a generator that provides augmented images and masks at the same time
train_generator = zip(image_generator, mask_generator)

callbacks = [
    EarlyStopping(patience=10, verbose=1),
    ReduceLROnPlateau(factor=0.1, patience=3, min_lr=0.000001, verbose=1),
    ModelCheckpoint('model-unet.h5', verbose=1, save_best_only=True, save_weights_only=True)
]

results = model.fit_generator(train_generator, steps_per_epoch=(len(X_train) // bs), epochs=100, callbacks=callbacks,
                              validation_data=(X_train, y_train))

model.save("model-unet-new1-without-tagg.h5")

path_test = r'D:\job\LAL\data\data-training\test/'
X_valid, y_valid = get_data(path_test, train=True)

# Evaluate on validation set (this must be equals to the best log_loss)
model.evaluate(X_valid, y_valid, verbose=1)

# Predict on train, val and test
preds_train = model.predict(X_train, verbose=1)
preds_val = model.predict(X_valid, verbose=1)

# Threshold predictions
#preds_train_t = (preds_train > 0.5).astype(np.uint8)
preds_val_t = (preds_val > 0.5).astype(np.uint8)

def plot_sample(X, y, preds, binary_preds, ix=None):
    if ix is None:
        ix = random.randint(0, len(X))

    has_mask = y[ix].max() > 0

    fig, ax = plt.subplots(1, 4, figsize=(20, 10))
    ax[0].imshow(X[ix, ..., 0], cmap='seismic')
    if has_mask:
        ax[0].contour(y[ix].squeeze(), colors='k', levels=[0.5])
    ax[0].set_title('Original-Dummy-colors')

    ax[1].imshow(y[ix].squeeze())
    ax[1].set_title('Ground-truth Mask')

    ax[2].imshow(preds[ix].squeeze(), vmin=0, vmax=1)
    if has_mask:
        ax[2].contour(y[ix].squeeze(), colors='k', levels=[0.5])
    ax[2].set_title('Mask Predicted')

    ax[3].imshow(binary_preds[ix].squeeze(), vmin=0, vmax=1)
    if has_mask:
        ax[3].contour(y[ix].squeeze(), colors='k', levels=[0.5])
    ax[3].set_title('Predicted binary');

def plot_sample1(X, y, binary_preds, ix=None):
    if ix is None:
        ix = random.randint(0, len(X))

    has_mask = y[ix].max() > 0

    fig, ax = plt.subplots(1, 2, figsize=(20, 10))
    ax[0].imshow(y[ix].squeeze())
    ax[0].set_title('Ground-truth Mask')
    ax[1].imshow(binary_preds[ix].squeeze(), vmin=0, vmax=1)
    if has_mask:
        ax[1].contour(y[ix].squeeze(), colors='k', levels=[0.5])
    ax[1].set_title('Predicted binary');

plot_sample(X_valid, y_valid, preds_val, preds_val_t, ix=2)
plot_sample1(X_valid, y_valid, preds_val_t, ix=2)

max_stars_repo_name: GoFigure-LANL/figure-segmentation

max_issues_repo_name: GoFigure-LANL/figure-segmentation

max_forks_repo_name: GoFigure-LANL/figure-segmentation


text:
# This file is part of the pyMOR project (https://www.pymor.org).
# Copyright pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)

import hashlib
import os
import sys
from pprint import pformat
from functools import wraps

import hypothesis
import numpy as np
from pickle import dump, load
from pkg_resources import resource_filename, resource_stream

from pymor.algorithms.basic import almost_equal, relative_error

def runmodule(filename):
    import pytest

    sys.exit(pytest.main(sys.argv[1:] + [filename]))

def check_results(test_name, params, results, *args):
    params = str(params)
    tols = (1e-13, 1e-13)
    keys = {}
    for arg in args:
        if isinstance(arg, tuple):
            assert len(arg) == 2
            tols = arg
        else:
            keys[arg] = tols

    assert results is not None
    assert set(keys.keys()) <= set(results.keys()), \
        f'Keys {set(keys.keys()) - set(results.keys())} missing in results dict'
    results = {k: np.asarray(results[k]) for k in keys.keys()}
    assert all(v.dtype != object for v in results.values())

    basepath = resource_filename('pymortests', 'testdata/check_results')
    arg_id = hashlib.sha1(params.encode()).hexdigest()
    filename = resource_filename('pymortests', f'testdata/check_results/{test_name}/{arg_id}')
    testname_dir = os.path.join(basepath, test_name)

    def _dump_results(fn, res):
        with open(fn, 'wb') as f:
            f.write((params + '\n').encode())
            res = {k: v.tolist() for k, v in res.items()}
            dump(res, f, protocol=2)

    try:
        with resource_stream('pymortests', f'testdata/check_results/{test_name}/{arg_id}') as f:
            f.readline()
            old_results = load(f)
    except FileNotFoundError:
        if not os.path.exists(testname_dir):
            os.mkdir(testname_dir)
        _dump_results(filename, results)
        assert False, \
            f'No results found for test {test_name} ({params}), saved current results. Remember to check in {filename}.'

    for k, (atol, rtol) in keys.items():
        if not np.all(np.allclose(old_results[k], results[k], atol=atol, rtol=rtol)):
            abs_errs = np.abs(results[k] - old_results[k])
            rel_errs = abs_errs / np.abs(old_results[k])
            _dump_results(filename + '_changed', results)
            assert False, (f'Results for test {test_name}({params}, key: {k}) have changed.\n'
                           f'(maximum error: {np.max(abs_errs)} abs / {np.max(rel_errs)} rel).\n'
                           f'Saved new results in {filename}_changed')

def assert_all_almost_equal(U, V, product=None, sup_norm=False, rtol=1e-14, atol=1e-14):
    cmp_array = almost_equal(U, V, product=product, sup_norm=sup_norm, rtol=rtol, atol=atol)
    too_large_relative_errors = dict((i, relative_error(u, v, product=product))
                                     for i, (u, v, f) in enumerate(zip(U, V, cmp_array)) if not f)
    assert np.all(cmp_array), f'Relative errors for not-equal elements:{pformat(too_large_relative_errors)}'

def might_exceed_deadline(deadline=-1):
    """For hypothesis magic to work properly this must be the topmost decorator on test function"""
    def _outer_wrapper(func):
        @wraps(func)
        def _inner_wrapper(*args, **kwargs):
            dl = deadline
            if os.environ.get('PYMOR_ALLOW_DEADLINE_EXCESS', False):
                dl = None
            elif dl == -1:
                dl = hypothesis.settings.default.deadline.total_seconds() * 1e3
            assert dl is None or dl > 1
            return hypothesis.settings(deadline=dl)(func)(*args, **kwargs)
        return _inner_wrapper
    return _outer_wrapper

max_stars_repo_name: kinnala/pymor

max_issues_repo_name: kinnala/pymor


text:
'''Tests for util module'''

from unittest.mock import patch, MagicMock
from memn2n import util

import numpy as np

FILE_CONTENTS = '''1 Mary got the milk.
2 John moved to the bedroom.
3 Daniel journeyed to the office.
4 John grabbed the apple there.
5 John got the football.
6 John journeyed to the garden.
7 Mary left the milk.
8 John left the football.
9 Daniel moved to the garden.
10 Daniel grabbed the football.
11 Mary moved to the hallway.
12 Mary went to the kitchen.
13 John put down the apple there.
14 John picked up the apple.
15 Sandra moved to the hallway.
16 Daniel left the football there.
17 Daniel took the football.
18 John travelled to the kitchen.
19 Daniel dropped the football.
20 John dropped the apple.
21 John grabbed the apple.
22 John went to the office.
23 Sandra went back to the bedroom.
24 Sandra took the milk.
25 John journeyed to the bathroom.
26 John travelled to the office.
27 Sandra left the milk.
28 Mary went to the bedroom.
29 Mary moved to the office.
30 John travelled to the hallway.
31 Sandra moved to the garden.
32 Mary moved to the kitchen.
33 Daniel took the football.
34 Mary journeyed to the bedroom.
35 Mary grabbed the milk there.
36 Mary discarded the milk.
37 John went to the garden.
38 John discarded the apple there.
39 Where was the apple before the bathroom? 	office	38 25 22
40 Sandra travelled to the bedroom.
41 Daniel moved to the bathroom.
42 Where was the apple before the hallway? 	office	38 30 26
43 Sandra got the milk.
44 Daniel travelled to the garden.
45 Where was the apple before the hallway? 	office	38 30 26
46 Sandra went back to the bathroom.
47 Daniel took the apple there.
48 Mary went back to the hallway.
49 Daniel went to the hallway.
50 Sandra went to the kitchen.
51 Mary journeyed to the bedroom.
52 Sandra journeyed to the hallway.
53 Daniel put down the apple.
54 Daniel put down the football there.
55 Sandra journeyed to the garden.
56 Where was the football before the garden? 	bathroom	54 44 41
57 Mary travelled to the office.
58 Sandra dropped the milk.
59 Where was the football before the garden? 	bathroom	54 44 41
1 Mary got the milk.
2 John moved to the bedroom.
3 Daniel journeyed to the office.
4 John grabbed the apple there.
5 John got the football.
6 John journeyed to the garden.
7 Mary left the milk.
8 John left the football.
9 Daniel moved to the garden.
10 Daniel grabbed the football.
11 Mary moved to the hallway.
12 Mary went to the kitchen.
13 John put down the apple there.
14 John picked up the apple.
15 Sandra moved to the hallway.
16 Daniel left the football there.
17 Daniel took the football.
18 John travelled to the kitchen.
19 Daniel dropped the football.
20 John dropped the apple.
21 John grabbed the apple.
22 John went to the office.
23 Sandra went back to the bedroom.
24 Sandra took the milk.
25 John journeyed to the bathroom.
26 John travelled to the office.
27 Sandra left the milk.
28 Mary went to the bedroom.
29 Mary moved to the office.
30 John travelled to the hallway.
31 Sandra moved to the garden.
32 Mary moved to the kitchen.
33 Daniel took the football.
34 Mary journeyed to the bedroom.
35 Mary grabbed the milk there.
36 Mary discarded the milk.
37 John went to the garden.
38 John discarded the apple there.
39 Where was the apple before the bathroom? 	office	38 25 22'''

def test_load_dataset():
    '''Test dataset load'''

    m = MagicMock(name='open')
    m().__enter__().__iter__.return_value = iter(FILE_CONTENTS.splitlines())
    with patch('memn2n.util.open', m):
        batch = list(util.load_dataset('filename'))

        assert len(batch) == 6
        sample = batch[0]
        assert len(sample) == 3
        assert len(sample[0]) == 38
        assert sample[1] == 'Where was the apple before the bathroom ?'
        assert sample[2] == 'office'

def test_memory_capacity_calc():
    m = MagicMock(name='open')
    m().__enter__().__iter__.return_value = iter(FILE_CONTENTS.splitlines())
    with patch('memn2n.util.open', m):
        dataset = list(util.load_dataset('filename'))
        capacity = util.calc_memory_capacity_for(dataset)
        assert capacity == 59 - 5

def test_sentence_length_calc():
    m = MagicMock(name='open')
    m().__enter__().__iter__.return_value = iter(FILE_CONTENTS.splitlines())
    with patch('memn2n.util.open', m):
        dataset = list(util.load_dataset('filename'))
        capacity = util.calc_sentence_length_for(dataset)
        assert capacity == 8

def test_vectorize_dataset():
    from nltk.tokenize import word_tokenize as tokenize
    m = MagicMock(name='open')
    m().__enter__().__iter__.return_value = iter(FILE_CONTENTS.splitlines())
    with patch('memn2n.util.open', m):
        dataset = list(util.load_dataset('filename'))

        memory_size, sentence_length = 20, 15
        # fake vocabulary to make sure we get correct memories
        word2idx = {
           'John': 1,
           'discarded': 2,
           'the': 3,
           'apple': 4,
           'there': 5,
           '.': 6
        }
        fcts, q, a = util.vectorize_dataset(dataset, word2idx, memory_size, sentence_length)

        assert fcts.shape == (6, memory_size, sentence_length)
        assert q.shape == (6, sentence_length)
        # last memory from particular sample
        assert np.all(fcts[-1][-1] == [1, 2, 3, 4, 5, 6] + [0, 0, 0, 0, 0, 0, 0, 0, 0])


text:
import numpy as np
np.seterr(all='ignore')
import pyfits
import atpy

from hyperion.util.constants import c, kpc
from hyperion.util.integrate import integrate_loglog

# Set list of filters
filt_name = ['B', 'V', 'J', 'H', 'K', '12', '25']

# Set wavelengths in microns
filt_wav = np.array([0.44, 0.55, 1.235, 1.662, 2.159, 12., 25.])

# Read in spectral information
ts = atpy.TableSet('SKY/models/sky_model_spectra.fits', verbose=False)

# Extract default table
t_default = ts['DEFAULT FLUXES']

# Find number of spectral types
n_spec = len(t_default)

# Compile list of spectral types
spectral_types = [t_default['Type'][ispec].strip().upper() for ispec in range(n_spec)]

sky_loaded = {}

def setup_sources(model, model_name, sky_model, spectral_type=None):

    global sky_loaded

    # Read in SKY model maps
    if sky_model in sky_loaded:
        print "[%s] Reading %s (cached)" % (model_name, sky_model)
        hdulist = sky_loaded[sky_model]
    else:
        print "[%s] Reading %s" % (model_name, sky_model)
        hdulist = pyfits.open(sky_model)
        sky_loaded[sky_model] = hdulist

    # Set cylindrical grid parameters
    r_wall = hdulist[1].data * kpc
    z_wall = hdulist[2].data * kpc
    p_wall = hdulist[3].data
    model.set_cylindrical_polar_grid(r_wall, z_wall, p_wall)

    # Loop through spectral types
    for ispec in range(n_spec):

        if spectral_type is not None and t_default['Type'][ispec].strip().upper() != spectral_type:
            continue
        else:
            spec_type = t_default['Type'][ispec].strip().upper()

        # Check if full spectrum exists
        if spec_type in ts.tables:
            nu = ts[spec_type]['nu']
            fnu = ts[spec_type]['fnu']
        else:
            nu = c / (filt_wav * 1.e-4)
            fnu = np.array([t_default[band][ispec] for band in filt_name])

        # Sort spectrum
        order = np.argsort(nu)
        nu, fnu = nu[order], fnu[order]

        # Find total luminosity of single object
        l_single = integrate_loglog(nu, fnu)

        # Find number of sources in each cell
        number = hdulist[0].data[ispec, :, :, :]

        # Set up source
        source = model.add_map_source()
        source.luminosity = np.sum(number) * l_single
        source.spectrum = (nu.astype(np.float64), fnu.astype(np.float64))
        source.map = number

    return model

max_stars_repo_name: hyperion-rt/paper-galaxy-rt-model

max_issues_repo_name: hyperion-rt/paper-galaxy-rt-model

max_forks_repo_name: hyperion-rt/paper-galaxy-rt-model


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""A module that implements the Gated Recurrent Unit layer.
"""

import theano
import theano.tensor as tensor

from theanolm.network.weightfunctions import get_submatrix
from theanolm.network.basiclayer import BasicLayer

class GRULayer(BasicLayer):
    """Gated Recurrent Unit Layer

    K. Cho et al. (2014)
    Learning Phrase Representations Using RNN Encoder-Decoder for Statistical
    Machine Translation
    Proc. 2014 Conference on Empiricial Methods in Natural Language Processing
    """

    def __init__(self, *args, **kwargs):
        """Initializes the parameters used by this layer.

        The weight matrices are concatenated so that they can be applied in a
        single parallel matrix operation. The same thing for bias vectors.
        """

        super().__init__(*args, **kwargs)

        input_size = sum(x.output_size for x in self._input_layers)
        output_size = self.output_size

        # Add state variables to be passed between time steps.
        self.hidden_state_index = self._network.add_recurrent_state(output_size)

        # Initialize the parameters.
        num_gates = 2
        # layer input weights for each gate and the candidate state
        self._init_weight('layer_input/W', (input_size, output_size),
                          scale=0.01, count=num_gates+1)
        # hidden state input weights for each gate and the candidate state
        self._init_weight('step_input/W', (output_size, output_size),
                          count=num_gates+1, split_to_devices=False)
        # biases for each gate and the candidate state
        self._init_bias('layer_input/b', output_size, [None] * (num_gates + 1))

        self.output = None

    def create_structure(self):
        """Creates the symbolic graph of this layer.

        The input is always 3-dimensional: the first dimension is the time step,
        the second dimension are the sequences, and the third dimension is the
        layer input. When processing mini-batches, all dimensions can have size
        greater than one.

        The function can also be used to create a structure for generating the
        probability distribution of the next word. Then the input is still
        3-dimensional, but the size of the first dimension (time steps) is 1,
        and the state outputs from the previous time step are read from
        ``self._network.recurrent_state_input``.

        Saves the recurrent state in the Network object. There's just one state
        in a GRU layer, h_(t). ``self.output`` will be set to the same hidden
        state output, which is also the actual output of this layer.
        """

        layer_input = tensor.concatenate([x.output for x in self._input_layers],
                                         axis=2)
        num_time_steps = layer_input.shape[0]
        num_sequences = layer_input.shape[1]

        # Compute the gate and candidate state pre-activations, which don't
        # depend on the state input from the previous time step.
        layer_input_preact = self._tensor_preact(layer_input, 'layer_input')
        if self._reverse_time:
            # Shift the input left by two time steps in the backward layer of a
            # bidirectional layer. The target word is the next input word, so
            # otherwise predicting it will be trivial.
            preact_dim = layer_input_preact.shape[2]
            padding = tensor.zeros([2, num_sequences, preact_dim])
            layer_input_preact = tensor.concatenate(
                [layer_input_preact[2:], padding],
                axis=0)

        # Weights of the hidden state input of each time step have to be applied
        # inside the loop.
        hidden_state_weights = self._get_param('step_input/W')

        if self._network.mode.minibatch:
            sequences = [self._network.mask, layer_input_preact]
            non_sequences = [hidden_state_weights]
            initial_hidden_state = tensor.zeros(
                (num_sequences, self.output_size), dtype=theano.config.floatX)

            hidden_state_output, _ = theano.scan(
                self._create_time_step,
                sequences=sequences,
                outputs_info=[initial_hidden_state],
                non_sequences=non_sequences,
                go_backwards=self._reverse_time,
                name='gru_steps',
                n_steps=num_time_steps,
                profile=self._profile,
                strict=True)

            self.output = hidden_state_output
            if self._reverse_time:
                self.output = self.output[::-1]
        elif self._reverse_time:
            raise RuntimeError("Text generation and lattice decoding are not "
                               "possible with bidirectional layers.")
        else:
            hidden_state_input = \
                self._network.recurrent_state_input[self.hidden_state_index]

            hidden_state_output = self._create_time_step(
                self._network.mask[0],
                layer_input_preact[0],
                hidden_state_input[0],
                hidden_state_weights)

            # Create a new axis for time step with size 1.
            hidden_state_output = hidden_state_output[None, :, :]

            self._network.recurrent_state_output[self.hidden_state_index] = \
                hidden_state_output
            self.output = hidden_state_output

    def _create_time_step(self, mask, x_preact, h_in, h_weights):
        """The GRU step function for theano.scan(). Creates the structure of one
        time step.

        The inputs do not contain the time step dimension. ``mask`` is a vector
        containing a boolean mask for each sequence. ``x_preact`` is a matrix
        containing the preactivations for each sequence. ``C_in`` and ``h_in``,
        as well as the outputs, are matrices containing the state vectors for
        each sequence.

        The required affine transformations have already been applied to the
        input prior to creating the loop. The transformed inputs and the mask
        that will be passed to the step function are vectors when processing a
        mini-batch - each value corresponds to the same time step in a different
        sequence.

        :type mask: Variable
        :param mask: a symbolic vector that masks out sequences that are past
                     the last word

        :type x_preact: Variable
        :param x_preact: concatenation of the input x_(t) pre-activations
                         computed using the gate and candidate state weights and
                         biases; shape is (the number of sequences, state size *
                         3)

        :type h_in: Variable
        :param h_in: h_(t-1), hidden state output of the previous time step;
                     shape is (the number of sequences, state size)

        :type h_weights: Variable
        :param h_weights: concatenation of the gate and candidate state weights
                          to be applied to h_(t-1); shape is (state size, state
                          size * 3)

        :rtype: Variable
        :returns: h_(t), the hidden state output
        """

        # pre-activation of the gates
        h_preact = tensor.dot(h_in, h_weights)
        preact_gates = get_submatrix(h_preact, 0, self.output_size, 1)
        preact_gates += get_submatrix(x_preact, 0, self.output_size, 1)

        # reset and update gates
        r = tensor.nnet.sigmoid(
            get_submatrix(preact_gates, 0, self.output_size))
        u = tensor.nnet.sigmoid(
            get_submatrix(preact_gates, 1, self.output_size))

        # pre-activation of the candidate state
        preact_candidate = get_submatrix(h_preact, 2, self.output_size)
        preact_candidate *= r
        preact_candidate += get_submatrix(x_preact, 2, self.output_size)

        # hidden state output
        h_candidate = tensor.tanh(preact_candidate)
        h_out = (1.0 - u) * h_in + u * h_candidate

        # Apply the mask. None creates a new axis with size 1, causing the mask
        # to be broadcast to all the outputs.
        h_out = tensor.switch(mask[:, None], h_out, h_in)

        return h_out

max_stars_repo_name: vasuneralla/theanolm

reason: import theano,from theano


text:
import numpy as np
import random

"""
canyon 是自己写的 optimizer, 可以无缝替换到 controller.py 中 Hive.BeeHive 的位置
思路就是模仿加农炮打靶, 根据 Reward 不断优化自己的 action
"""

class Canyon(object):
    def __init__(self, lower, upper, fun):
        self.lower = lower
        self.upper = upper
        self.horizen = len(self.lower)
        self.evaluate = fun
        self.solution = None

    def run(self):
        epsilon = 0.1  # 只要 action 的变化量大于 epsilon 就说明没有收敛
        action_delta = 1  # action 更新的变化量
        reward_episode = 0
        render = False
        left_action = self.lower
        right_action = self.upper

        cur_action = []
        for tmp in range(self.horizen):
            action_0 = left_action[tmp][0] + random.random() * (right_action[tmp][0] - left_action[tmp][0])
            action_1 = left_action[tmp][1] + random.random() * (right_action[tmp][1] - left_action[tmp][1])
            cur_action.append( [action_0, action_1] )
        cur_action = np.array(cur_action)

        label_tmp = []

        while(action_delta > epsilon):
            action_delta = 0
            for step in range(self.horizen):  # 依次更新每一个 action
                left_reward = self.evaluate(left_action)
                right_reward = self.evaluate(right_action)
                if left_reward > right_reward:
                    action_delta += ((cur_action[step][0] - right_action[step][0])**2+(cur_action[step][1] - right_action[step][1])**2)**0.5
                    left_action[step] = cur_action[step]
                    cur_action[step] = (right_action[step] + cur_action[step]) / 2
                else:
                    action_delta += ((cur_action[step][0] - right_action[step][0])**2+(cur_action[step][1] - right_action[step][1])**2)**0.5
                    right_action[step] = cur_action[step]
                    cur_action[step] = (left_action[step] + cur_action[step]) / 2
        self.solution = cur_action
        return cur_action

max_stars_repo_name: WeixiongLin/Quanser-robots

max_issues_repo_name: WeixiongLin/Quanser-robots

max_forks_repo_name: WeixiongLin/Quanser-robots


text:
# Author: Behrad Soleimani <behrad@umd.edu>
# Author: Proloy Das <proloy@umd.edu>
"Deviance calculation"

import numpy as np
import warnings
from scipy import linalg

from .opt.m_step import calculate_ss

def sample_path_bias(q, a, x_bar, zeroed_index, n_eigenmodes):
    """Computes the bias in the deviance

    Parameters
    ----------
    q:  ndarray of shape (n_sources*mo, n_sources*mo)
    a:  ndarray of shape (n_sources*mo, n_sources*order*mo)
    x_bar:  ndarray of shape (t, n_sources*mo)
    idx_src: source index

    Returns
    -------
    bias

    """
    t, dxm = x_bar.shape
    _, dtot = a.shape
    p = dtot // dxm

    bias = 0
    qd = np.diag(q)
    cx = np.zeros((t - p, dtot))

    for idx_src in range(dxm):
        ai = a[idx_src]

        xi = x_bar[p:, idx_src]

        for k in range(p):
            cx[:, k * dxm:(k + 1) * dxm] = x_bar[p - 1 - k:t - 1 - k]

        # gradient of log - likelihood
        ldot = cx.T.dot(xi - cx.dot(ai)) / qd[idx_src]
        # hessian of log - likelihood
        ldotdot = -cx.T.dot(cx) / qd[idx_src]

        # if zeroed_index is not None:
        #     x_index, y_index = zeroed_index
        #     if idx_src in x_index:
        #         removed_idx = list(np.asanyarray(y_index)[np.asanyarray(x_index) == idx_src])
        #         ldot = np.delete(ldot, removed_idx)
        #         ldotdot = np.delete(ldotdot, removed_idx, axis=0)
        #         ldotdot = np.delete(ldotdot, removed_idx, axis=1)

        # FIX removing cross-talk components (that forced to be zero)
        for l in range(0, dxm, n_eigenmodes):
            for u in range(n_eigenmodes):
                for v in range(n_eigenmodes):
                    if v != u and idx_src == l + v:
                        removed_idx = list(range(l + u, dtot, dxm))
                        if zeroed_index is not None:
                            x_index, y_index = zeroed_index
                            if idx_src in x_index:
                                removed_idx.extend(list(np.asanyarray(y_index)[np.asanyarray(x_index) == idx_src]))
                        ldot = np.delete(ldot, removed_idx)
                        ldotdot = np.delete(ldotdot, removed_idx, axis=0)
                        ldotdot = np.delete(ldotdot, removed_idx, axis=1)

        bias += ldot.dot(np.linalg.solve(ldotdot, ldot))
    return bias

def bias_by_idx(idx_src, q, a, x_bar, s_bar, b, m, p, zeroed_index=None):
    """Computes the bias in the deviance (proloy@umd.edu)

    Parameters
    ----------
    q:  ndarray of shape (n_sources*mo, n_sources*mo)
    a:  ndarray of shape (n_sources*mo, n_sources*order*mo)
    x_bar:  ndarray of shape (t, n_sources*mo)
    idx_src: source index

    Returns
    -------
    bias

    """
    warnings.filterwarnings('always')
    _, dtot = a.shape

    ### These uses the whole distribution.
    s1, s2, s3, n = calculate_ss(x_bar, s_bar, b, m, p)

    ai = a[idx_src]  # in python slicing returns 1d array, so transpose is meaningless.
    qi = q[idx_src, idx_src]

    ldot = np.empty((dtot))
    ldotdot = np.empty((dtot, dtot))

    temp1 = s2.dot(ai)
    temp1 -= s1[idx_src]

    ldot[:] = - temp1
    ldot[:] /= qi

    ldotdot[:, :] = - s2 / qi

    if zeroed_index is not None:
        x_index, y_index = zeroed_index
        if idx_src in x_index:
            removed_idx = list(np.asanyarray(y_index)[np.asanyarray(x_index) == idx_src])
            ldot = np.delete(ldot, removed_idx)
            ldotdot = np.delete(ldotdot, removed_idx, axis=0)
            ldotdot = np.delete(ldotdot, removed_idx, axis=1)

    try:
        c, low = linalg.cho_factor(-ldotdot)
        temp = linalg.cho_solve((c, low), ldot)
        bias = n * ldot.dot(temp)
    except linalg.LinAlgError:
        warnings.warn('source-index {:d} ldotdot is not negative definite: '
                      'setting positive eigenvalues equal to zero, '
                      'result may not be accurate.'.format(idx_src), RuntimeWarning, stacklevel=2)
        e, v = linalg.eigh(-ldotdot)
        temp = v.dot(ldot)
        idx = e > 0
        bias = np.sum(temp[idx] ** 2 / e[idx])
        bias *= n

    return bias

def debias_deviances(dev_raw, bias_f, bias_r):
    d = dev_raw.copy()
    bias_mat = bias_r - bias_f
    d[bias_r != 0] += bias_mat[bias_r != 0]
    np.fill_diagonal(d, 0)
    d[d < 0] = 0
    return d

max_stars_repo_name: BabadiLab/NLGC

max_issues_repo_name: BabadiLab/NLGC

max_forks_repo_name: BabadiLab/NLGC