text:
import gym

from gym import spaces, logger
from gym.utils import seeding
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.patches import Rectangle, Circle
from brl_gym.envs.util import fig2np

# Crosswalk is 4 x 4 space, [0, 0] x [4, 4]
# Driver is approaching the crosswalk from y=0
# Driver actions are (velocity, steering angle).

# 3 pedestrians, randomly picked to be on the right or left.
# Pedestrians start anywhere vertically, y in [1, 6]. x is either 0 or 4.
# And go straight to their goals, with fixed velocity,
# uniformly sampled from [0.2 m/s, 0.8  m/s].
# Done when the car reaches y = 4.5.
# Angles are ccw, with 0 rad correponding to y-axis, the initial direction of the car.

#colors = ['#C7980A', '#F4651F', '#82D8A7', '#CC3A05', '#575E76', '#156943']
#colors = ['#88CCEE','#DDCC77','#117733','#332288','#44AA99','#661100']
#colors = ['#377eb8','#4daf4a','#984ea3','#ff7f00','#ffff33','#a65628', '#e41a1c']
from matplotlib import cm
tab10 = cm.get_cmap('tab10')
colors = np.vstack([tab10.colors[:3], tab10.colors[4:6], tab10.colors[7], tab10.colors[3]])

def check_intersect(seg1, seg2):
    p0, p1 = seg1
    q0, q1 = seg2

    # Vectors
    u = p1 - p0
    v = q1 - q0

    # Vector perpendicular to v
    v_perp = np.array([-v[1], v[0]])

    # Vector q0 -> p0
    w = p0 - q0

    # Check parallel
    denom = np.dot(v_perp, u)
    if np.abs(denom) <= 1e-3:
        return False

    # Find s such that dot(v_perp, w + s*u) = 0
    s = -np.dot(v_perp,w)/denom

    # Intersection = w + s*u
    i = p0 + s * u
    if s < 0 or s > 1:
        return False

    # Find r for which i = q0 + r * v
    r = (i - q0)[0] / v[0]

    if r < 0 or r > 1:
        return False

    return True

def min_dist_small(seg1, seg2, thresh=0.65):
    p0, p1 = seg1
    q0, q1 = seg2
    x = p1 - p0
    x_normalized = x / np.dot(x,x)

    seg1 = np.vstack(seg1)
    seg2 = np.vstack(seg2)

    t = np.linspace(0, 1.0, 11)
    q = q0 + (q1-q0).reshape(-1,2)*t.reshape(-1,1)
    d = p1 - q
    xp = x_normalized.reshape(-1,2)*np.dot(d, x).reshape(-1,1)
    dist = np.linalg.norm(d - xp, axis=1)

    return np.any(dist < thresh)
    # for t in np.linspace(0, 1.0, 5):
    #     q = q0 + (q1 - q0) * t
    #     d = p1 - q
    #     xp = np.dot(d, x) / np.dot(x,x) * x
    #     dist = np.linalg.norm(d - xp)

    # if t in [0.0, 0.5, 1.0]:
    #     print("dist", dist)
    #     plt.figure()
    #     plt.plot(seg1[:,0], seg1[:,1], 'r')
    #     plt.plot(seg2[:,0], seg2[:,1], 'b')
    #     # debug
    #     plt.plot([q[0],(p1-xp)[0]], [q[1], (p1-xp)[1]], 'r--')
    #     plt.show()
    #     plt.close()

    #     if dist < thresh:
    #         return True
    # return False

def draw_all(seg1, seg2, seg3, seg4):

    seg1 = np.vstack(seg1)
    seg2 = np.vstack(seg2)
    seg3 = np.vstack(seg3)
    seg4 = np.vstack(seg4)

    print(seg1)
    print(seg2)
    print(seg3)
    print(seg4)

    plt.figure()
    plt.plot(seg1[:,0], seg1[:,1],'r')
    plt.plot(seg2[:,0], seg2[:,1],'g')
    plt.plot(seg3[:,0], seg3[:,1],'b')
    plt.plot(seg4[:,0], seg4[:,1],'k')

    plt.scatter(seg1[0,0], seg1[0,1], c='r', s=30)
    plt.scatter(seg2[0,0], seg2[0,1], c='g', s=30)
    plt.scatter(seg3[0,0], seg3[0,1], c='b', s=30)
    plt.scatter(seg4[0,0], seg4[0,1], c='k', s=30)
    plt.xlim(0.0,3.5)
    plt.ylim(0.0,4.5)

    plt.show()

# Crosswalk but with (velocity, steering angle) control
class CrossWalkVelEnv(gym.Env):
    def __init__(self, timestep=0.05):
        self.action_space = spaces.Box(np.array([0.3, -0.3]), np.array([0.6, 0.3]))
        self.car_length = 0.4 # Length of the MuSHR car
        self.num_pedestrians = 3
        self.timestep = timestep
        self.observation_space = spaces.Box(low=np.ones(8+self.num_pedestrians*4)*-10.0,
                                            high=np.ones(8+self.num_pedestrians*4)*10.0,
                                            dtype=np.float32)
        self.x_limit = np.array([0, 3.5], dtype=np.float32)
        self.y_limit = np.array([0, 4], dtype=np.float32)
        self.car_y_goal = 4.5
        self.x_left = 0.0
        self.x_right = 3.5
        self.y_starts = np.array([1.5, 3.8])
        self.ped_speed_limits = np.array([0.5, 0.6], dtype=np.float32)
        self.time_penalty = 1./(10. / 0.05)

    def _get_init_goal(self):
        while True:
            goal_xs = np.zeros(self.num_pedestrians, dtype=np.float32)

            # Choose left or right for each pedestrian. # 0 for left, 1 for right
            #sides = np.random.choice(2, size=self.num_pedestrians)
            #if np.all(sides == 0) or np.all(sides == 1):
            #    continue
            sides = np.array([0, 0, 1])
            if np.random.random() > 0.5:
                sides[1] = 1

            peds = np.zeros((self.num_pedestrians, 2), dtype=np.float32)
            peds[sides == 0, 0] = self.x_left
            peds[sides == 1, 0] = self.x_right

            # Choose init-y
            peds[:, 1] = np.random.uniform(self.y_starts[0], self.y_starts[1], size=self.num_pedestrians)

            # Goals
            goals = np.zeros((self.num_pedestrians, 2), dtype=np.float32)
            goals[sides == 0, 0] = self.x_right
            goals[sides == 1, 0] = self.x_left
            goals[:, 1] = np.random.uniform(self.y_starts[0], self.y_starts[1], size=self.num_pedestrians)

            # If too close, pass
            close = False
            pairs = np.array([(0,1), (0,2), (1,2)])
            g = goals[pairs][:,:,1]
            p = peds[pairs][:,:,1]
            if np.any(np.abs(g[:,0] - g[:,1]) < 0.35) or np.any(np.abs(p[:,0] - p[:,1]) < 0.35):
                close = True

            # for pair in [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]:
            #     if (np.abs(goals[pair[0],1] - goals[pair[1],1]) < 0.35
            #        or np.abs(peds[pair[0],1] - peds[pair[1],1]) < 0.35):
            #         close = True
            if close:
                continue

            # Check if paths cross
            intersect = False
            for pair in pairs:
                ped1 = (peds[pair[0]], goals[pair[0]])
                ped2 = (peds[pair[1]], goals[pair[1]])
                if min_dist_small(ped1, ped2):
                    intersect = True
                    break

            if intersect:
                continue
            else:
                break
        # draw_all((peds[0], goals[0]),(peds[1], goals[1]), (peds[2], goals[2]), (peds[3], goals[3]))
        return goals, peds

    def reset(self):
        self.goals, self.pedestrians = self._get_init_goal()

        self.t = 0

        self.random_delays = np.random.choice(20, size=self.num_pedestrians)

        # Pedestrians have fixed speed, but can change directions
        self.pedestrian_speeds = np.random.uniform(
            self.ped_speed_limits[0], self.ped_speed_limits[1], size=self.num_pedestrians)
        self.pedestrian_angles = self._get_pedestrian_angles()

        # Agent's initial position, speed, angle
        self.pose = np.array([np.random.uniform(1.2, 2.3), 0.0, np.random.uniform(-0.1, 0.1)])
        self.speed = 0.0
        self.steering_angle = 0.0
        self.car_front = self.pose[:2] + \
                         self.car_length * np.array([-np.sin(self.pose[2]), np.cos(self.pose[2])])

        ped_speeds = self.pedestrian_speeds.copy()
        for i in range(self.num_pedestrians):
            if self.random_delays[i] < self.t:
                ped_speeds[i] = 0.0

        self.pedestrian_directions = self._get_pedestrian_directions(self.pedestrian_angles, ped_speeds)

        self.fig = None
        self.car = None
        return self.get_obs()

    def _get_pedestrian_angles(self, add_noise=True):
        # Angles are directed straight to the goals
        diff = self.goals - self.pedestrians
        angles = -1.0*np.arctan2(diff[:,0], diff[:,1])
        if add_noise:
            angles += np.random.normal(size=self.pedestrians.shape[0], scale=0.5)
        return angles

    def _get_actual_action(self, action):
        action = np.clip(action, self.action_space.low, self.action_space.high)
        action *= 0.8
        if action[0] < 0.3:
            action[0] = 0.0
        return action

    def step(self, action):
        self.t += 1
        assert isinstance(action, np.ndarray)
        action[0] += np.clip(np.random.normal(size=1)*0.1, -0.1, 0.1)
        action[1] += np.clip(np.random.normal(size=1)*0.05, -0.05, 0.05)
        action = np.clip(action, self.action_space.low, self.action_space.high)
        action *= 0.8
        if action[0] < 0.3:
            action[0] = 0.0

        self.speed = action[0]
        theta = self.pose[2] + np.pi / 2.0

        # move car
        next_theta = theta + self.speed / self.car_length * np.tan(1e-5+action[1]) * self.timestep
        self.pose[0] += self.car_length / (1e-5+np.tan(action[1])) * (np.sin(next_theta) - np.sin(theta))
        self.pose[1] += self.car_length / (1e-5+np.tan(action[1])) * (-np.cos(next_theta) + np.cos(theta))
        self.pose[2] = next_theta - np.pi / 2.0
        self.car_front = self.pose[:2] + \
                         self.car_length * np.array([-np.sin(self.pose[2]), np.cos(self.pose[2])])

        # move pedestrians
        ped_speeds = self.pedestrian_speeds.copy()

        for i in range(self.num_pedestrians):
            if self.pedestrians[i,0] >= self.x_limit[1]:
                ped_speeds[i] = 0.0
            if self.pedestrians[i,0] <= self.x_limit[0]:
                ped_speeds[i] = 0.0
            if self.random_delays[i] > self.t:
                ped_speeds[i] = 0.0

        # for i in [2]:
        #     if self.random_delays[i] > self.t:
        #         ped_speeds[i] = 0.0
        #     if self.pedestrians[i,0] <= self.x_limit[0]:
        #         ped_speeds[i] = 0.0
        self.pedestrian_angles = self._get_pedestrian_angles()
        self.pedestrian_directions = self._get_pedestrian_directions(self.pedestrian_angles, ped_speeds)
        self.pedestrians += self.pedestrian_directions

        obs = self.get_obs()
        reward, done = self._get_reward(obs)
        return obs, reward, done, dict(
                pedestrians=self.pedestrians,
                pedestrian_speeds = ped_speeds,
                pedestrian_angles = self.pedestrian_angles)

    def _get_reward(self, obs):
        done = False
        pose = obs[:2]
        angle = obs[7]
        car_front = obs[2:4]
        speed = obs[6]
        pedestrians = obs[8:8+self.num_pedestrians*2].reshape(-1, 2)
        next_pedestrians = obs[8+self.num_pedestrians*2:8+self.num_pedestrians*4].reshape(-1,2)

        # Penalize time & distance to goal. y = 4.0
        reward = -self.time_penalty - np.abs(self.car_y_goal - pose[1]) * 0.05

        # angle penalty
        reward -= np.abs(angle) * 0.05

        collision_dist = 0.65
        dist = np.linalg.norm(pose - pedestrians, axis=1)
        front_dist = np.linalg.norm(car_front - pedestrians, axis=1)
        next_dist = np.linalg.norm(pose - next_pedestrians, axis=1)

        reward *= 0.5

        # Collision
        collision = False
        # print(np.around(dist,2), np.around(next_dist, 2))
        for i, (d, fd, nd) in enumerate(zip(dist, front_dist, next_dist)):
            if (d < collision_dist or fd < collision_dist) and nd < d:
                collision = True
                break

        if collision:
            done = True
            reward = -10
        elif car_front[0] <= self.x_limit[0] + 0.5 \
                or car_front[0] >= self.x_limit[1] -0.5 \
                or car_front[1] <= -0.5:
            reward = -10.0
            done = True
        elif pose[1] >= self.car_y_goal:
            reward = 10.0
            done = True

        return reward, done

    def _get_pedestrian_directions(self, angles, ped_speeds):
        speeds = ped_speeds.reshape(-1,1).copy()
        speeds[speeds < 0.1] = 0.5
        speeds += np.clip(np.random.normal(size=speeds.shape)*0.1, -.1, 0.1)
        speeds = np.clip(speeds, 0, 1.0)

        return self.timestep*speeds \
                * np.array([-np.sin(angles),
                            np.cos(angles)]).transpose()

    def get_obs(self):
        car_direction = self.car_front - self.pose[:2]

        obs = np.concatenate([
            self.pose[:2].ravel(), # 2
            self.car_front.ravel(), # 4
            (self.car_front + car_direction).ravel(), # 6
            [self.speed, self.pose[2]], # 8
            self.pedestrians.ravel(), # +4*2 = 16
            (self.pedestrians + 1.0/self.timestep*self.pedestrian_directions).ravel()]) # + 4*2 = 24
        return obs

    def _visualize(self, show=False, filename=None, nparray=False, head_only=False, transparent=True):
        # if filename is None:
        #     plt.ion()
        #     self.fig = plt.figure()
        fig = plt.figure()

        # Draw goal
        plt.plot([0, 4], [self.car_y_goal, self.car_y_goal], linewidth=3, color='r')

        # Draw boundaries
        plt.plot([0, 0], [-1, 5], linewidth=2, color='k', linestyle='--')
        plt.plot([3.5, 3.5], [-1, 5], linewidth=2, color='k', linestyle='--')

        plt.axis('square')
        plt.axis('off')
        plt.xlim((-0.1, 3.6))
        plt.ylim((-0.5, self.car_y_goal+0.1))

        # Car
        # Pose of the end of the car
        car = self.pose.copy()
        angle = car[2]
        if self.car is not None:
            self.car.remove()
        car_front = self.car_front.copy()

        # self.car = Rectangle((car[0]-0.33/2.0, car[1]-self.car_length/2.0),
        #             0.33, self.car_length, angle=np.rad2deg(angle), color=colors[-1])
        # else:
        self.car = Circle([car[0], car[1]], radius=self.car_length, color=colors[-1], zorder=25)
        plt.plot([car[0], car_front[0]], [car[1], car_front[1]], color="k", lw=3, zorder=30)

        plt.gca().add_patch(self.car)

        # pedestrians
        for i in range(len(self.pedestrians)):
            ped = Circle(self.pedestrians[i], radius=0.33, color=colors[i], zorder=20)
            plt.gca().add_patch(ped)

            xy = self.pedestrians[i]
            dxdy = self.pedestrian_directions[i]

            #rect = Rectangle((xy[0], xy[1]),
            #        0.1, 1.0, angle=np.rad2deg(self.pedestrian_angles[i]), color=colors[i])
            #plt.gca().add_patch(rect)
        """
        for i in range(2):
            goal = Circle(self.goals[i], radius=0.25, color='g', zorder=10)
            plt.gca().add_patch(goal)
        """

        if nparray:
            raise NotImplementedError
            return fig2np(fig)

        if filename is not None:
            plt.savefig(filename, bbox_inches='tight', transparent=transparent)

        if show:
            fig.canvas.draw()
        plt.close()

    def _visualize_history(self, states, actions):
        pass

    def render(self, mode='human'):
        self._visualize(show=True)

if __name__ == "__main__":
    env = CrossWalkVelEnv()

    obs = env.reset()


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jul 24 09:57:01 2018

@author: yichen
"""

import pandas as pd
import numpy as np
import os
import re
import  difflib
from collections import Counter
import logging
import os
from gensim import corpora, utils
from gensim.models.wrappers.dtmmodel import DtmModel
import numpy as np
import scipy
import string

data = pd.read_csv('tc_record_log.csv')
data = data.loc[data['Volume']>=96]
data = data[['ArticleID','Date_Published_in_Print']]
data = data.drop_duplicates().dropna()

publish_year = []
for i in data['Date_Published_in_Print']:
    publish_year.append(i[:4])
data['Date_Published_in_Print'] = publish_year

files = os.listdir('./Data')
files_39_htm = []
for i in files:
    if (('39_' in i) or ('38_' in i)) and ('.htm' in i) and ('_g' not in i):
        files_39_htm.append(i)

articleid = list(data['ArticleID'])

final_files = []
final_articleid = []
problem = []
for i in files_39_htm:
    try:
        this = int(i[3:-4])
        if this in articleid:
            final_files.append(i)
            final_articleid.append(this)
    except:
        problem.append(i)

no_htm = list(set(articleid) - set(final_articleid))

documents = []
indexs = []
for index,file in enumerate(final_files):
    try:
        f=open("./Data/"+file,'r')
        document= f.read()
        documents.append(document)
        indexs.append(final_articleid[index])
    except:
        pass

documents_final = [0] * data.shape[0]
for i,j in enumerate(indexs):
    where = list(data['ArticleID']).index(j)
    documents_final[where] = documents[i]

data['document'] = documents_final
data = data.loc[data['document']!=0]
data = data.drop_duplicates().dropna()

data = data.reindex()

references = []
for i in data['document']:
    reference = re.findall(r'<P style=.*>.*\(.*\).* <I>Teachers College Record',i)
    references.append(reference)

for i,j in enumerate(references):
    if len(j) != 0:
        this = []
        for k in j:
            index_begin = k.index(')')
            index_end = k.index('<I>')
            potential = k[index_begin+1:index_end].replace('.','').strip()
            translator=str.maketrans('','',string.punctuation)
            potential = potential.translate(translator).lower()
            if len(potential) != 0 :
                this.append(potential)
            else:
                this.append(k)
        references[i] = this

data_1 = pd.read_csv('tc_record_log.csv')
title = list(data_1['Title'])
title = [str(i) for i in title]

Articleid = list(data_1['ArticleID'])

for i,a in enumerate(references):
    if len(a) != 0:
        for j,b in enumerate(a):
            references[i][j] = Articleid[title.index(difflib.get_close_matches(b, title,1,0)[0])]

final_references = []
for i in references:
    final_references += i

dic = Counter(final_references)

reference = [0] * data.shape[0]
for key,value in dic.items():
    try:
        where =  list(data['ArticleID']).index(key)
        reference[where] = value
    except:
        pass
data['reference'] = reference

data.to_csv('data_reference_doc.csv',index=False)

dic_time = Counter(data['Date_Published_in_Print'])

time = list(data['Date_Published_in_Print'])
time_step = []

dic_time = Counter(time_step)

for i in time:
    if int(i) >= 1994 and int(i) <= 1995:
        time_step.append(1)
    elif int(i) >= 1999 and int(i) <= 2003:
        time_step.append(2)
    elif int(i) >= 2004 and int(i) <= 2008:
        time_step.append(3)
    elif int(i) >= 2009 and int(i) <= 2013:
        time_step.append(4)
    else:
        time_step.append(5)

data['time_step'] = time_step
data = data.sort_values(by=['Date_Published_in_Print'])

data.to_csv('data_reference_doc_DIM.csv',index=False)

time_seq = [276, 535, 864, 482, 545]

dtm_path = '/Users/yichen/dtm-master/dtm/'

document = pd.read_csv('document.csv')
del document['Unnamed: 0']

document = document.iloc[0]
documents = []
wrong_index = []
for i in range(2702):
    try:
        this_doc = document[i].split(" ")
        this_doc = [x for x in this_doc if x != '']
        documents.append(this_doc)
        wrong_index.append(0)
    except:
        wrong_index.append(1)

data['wrong'] = wrong_index
data = data[data['wrong']==0]

data['document'] = document

time_seq = [276, 534, 864, 482, 544]

class DTMcorpus(corpora.textcorpus.TextCorpus):

    def get_texts(self):
        return self.input

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

corpus = DTMcorpus(documents)
model = DtmModel(dtm_path, corpus, time_seq, num_topics=2,id2word=corpus.dictionary, initialize_lda=True, model='fixed')

## document topic propotion

np.save('gamma',model.gamma_)
np.save('influence',model.influences_time)
np.save('model_result',model)

ext: py

reason: import numpy,import scipy


text:
""" find peaks in a periodic signal such as ECG
    Assume:
      1. cycle is stable and have one strong period signal
      2. target peaks have enough features
    Proprosal:
      1. peak value features
      2. cycle priors
      3. cycle transition priors
      4. autodetect to detech peaks or valleys
    Implemention:
      1. merge findpeaks results
      2. analysis peak results
      3. adjust findpeaks function parameters
"""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import itertools
import numpy as np
import findpeaks2
from filter import BandpassFilter
from filter import NotchFilter
from collections import deque
from collections import namedtuple
from absl import app

import pdb

class FindPeaksMaster(object):

    def __init__(self, fs=500):
        self._fs = fs
        self._bandpass_filter = BandpassFilter(fs)
        self._notch_filter = NotchFilter(fs)

        self._workers = []
        threses = [0.2]
        min_dists = [int(0.6*fs)]
        directions = ["up", "down"]
        buffer_sizes = [4*fs]
        combinations = itertools.product(threses, min_dists, directions, buffer_sizes)
        for thres, min_dist, direction, buffer_size in combinations:
            findpeaks = FindPeaks(fs=fs,
                                  thres=thres,
                                  min_dist=min_dist,
                                  direction=direction,
                                  buffer_size=buffer_size)
            worker = FindPeaksWorker(findpeaks, nsamples=6)
            self._workers.append(worker)
        self._peak_index_buffer = deque()
        self._peak_interval_buffer = deque()
        self._peak_value_buffer = deque()
        self._counter = 0
        self._counter_maximum = 1*fs

    def findpeaks(self, value):
        value = self._notch_filter.filter([value])[0]
        value = self._bandpass_filter.filter([value])[0]
        self._counter += 1
        for worker in self._workers:
            worker.findpeaks(value)
        if self._counter < self._counter_maximum:
            return
        self._counter = 0
        losses = []
        for worker in self._workers:
            losses.append(worker.loss)
        best_worker_idx = np.argmin(losses)
        best_worker = self._workers[best_worker_idx]
        self._peak_index_buffer.extend(best_worker.peak_indices)
        self._peak_interval_buffer.extend(best_worker.peak_intervals)
        self._peak_value_buffer.extend(best_worker.peak_values)
        for worker in self._workers:
            worker.peak_indices
            worker.peak_values
            worker.peak_intervals

    @property
    def peak_indices(self):
        peak_indices = list(self._peak_index_buffer)
        self._peak_index_buffer.clear()
        return peak_indices

    @property
    def peak_intervals(self):
        intervals = list(self._peak_interval_buffer)
        self._peak_interval_buffer.clear()
        return intervals

    @property
    def peak_values(self):
        peak_values = list(self._peak_value_buffer)
        self._peak_value_buffer.clear()
        return peak_values

class FindPeaksWorker(object):

    def __init__(self, instance, nsamples=60):
        """ master fetch indices must be more than 1 miniutes one time"""
        self._instance = instance
        self._peak_indices_for_loss = deque(maxlen=nsamples)
        self._peak_index_buffer = deque()
        self._peak_interval_buffer = deque()
        self._peak_value_buffer = deque()
        self._last_peak_index = None

    def findpeaks(self, value):
        self._instance.findpeaks(value)

    @property
    def mean_and_std_intervals(self):
        intervals = np.diff(self._peak_indices_for_loss)
        fs = self._instance.fs
        hist,bin_edges = np.histogram(intervals, bins=5, range=(0.5*fs, 1.5*fs))
        ind = np.argmax(hist)
        bin_st, bin_ed = bin_edges[ind], bin_edges[ind+1]
        cluster = []
        for interval in intervals:
            if interval >= bin_st and interval <= bin_ed:
                cluster.append(interval)
        mean = np.median(cluster)
        #std = np.std(np.abs(np.diff(cluster)))
        std = np.std(cluster)
        return mean, std

    @property
    def loss(self):
        '''
        set loss to be error rate
        '''
        loss = 1000000000
        peak_indices = self._instance.peak_indices
        self._peak_indices_for_loss.extend(peak_indices)

        if self._last_peak_index is None:
            if len(peak_indices) > 1:
                diff = peak_indices[1] - peak_indices[0]
                self._last_peak_index = peak_indices[0] - diff
            elif len(peak_indices) == 1:
                self._last_peak_index = peak_indices[0] - self._instance.fs
            else:
                return loss

        if self._peak_index_buffer:
            last_peak_index = self._peak_index_buffer[-1]
        else:
            last_peak_index = self._last_peak_index

        self._peak_interval_buffer.extend(np.diff([last_peak_index] + peak_indices))
        self._peak_index_buffer.extend(peak_indices)

        peak_values = self._instance.peak_values
        self._peak_value_buffer.extend(peak_values)

        if len(self._peak_indices_for_loss) >= 2:
            #mean, std = self.mean_and_std_intervals
            intervals = np.diff(self._peak_indices_for_loss)
            #diff_abs = np.abs(np.diff(intervals))
            #mean = np.median(diff_abs)
            #diff = diff_abs - mean
            #loss = np.sum(diff * diff) / len(diff)
            #inerval_diffs_abs = np.abs(np.diff(intervals))
            #count = 0
            #for interval in intervals:
            #    if np.abs(interval - mean) >= 5*std:
            #        count += 1
            #loss = std + 100*count
            loss = np.std(intervals)
        return loss

    @property
    def peak_indices(self):
        peak_indices = list(self._peak_index_buffer)
        self._peak_index_buffer.clear()
        if len(peak_indices) > 0:
            self._last_peak_index = peak_indices[-1]
        return peak_indices

    @property
    def peak_intervals(self):
        intervals = list(self._peak_interval_buffer)
        self._peak_interval_buffer.clear()
        return intervals

    @property
    def peak_values(self):
        peak_values = list(self._peak_value_buffer)
        self._peak_value_buffer.clear()
        return peak_values

class FindPeaks(object):
    """ online findpeaks for periodic signals """
    def __init__(self, fs=500, thres=0.65, min_dist=300, direction="up", buffer_size=None):
        if buffer_size is None:
            buffer_size = int(fs * 2)
        self.fs = fs
        self._signal_buffer = deque(maxlen=buffer_size)
        self._peak_value_buffer = deque()
        self._peak_index_buffer = deque()
        self._peaks_buffer = deque()
        self._raw_extremas_buffer = deque()
        self._raw_peaks_buffer = deque()
        self._thres_abs_buffer = deque()
        self._signal_index = -1 # starts from 0
        self._update_counter = 0
        self.thres = thres
        self.min_dist = min_dist
        self.direction = direction

    @property
    def raw_extremas(self):
        extremas = list(self._raw_extremas_buffer)
        self._raw_extremas_buffer.clear()
        return extremas

    @property
    def raw_peaks(self):
        peaks = list(self._raw_peaks_buffer)
        self._raw_peaks_buffer.clear()
        return peaks

    @property
    def peaks(self):
        peaks = list(self._peaks_buffer)
        self._peaks_buffer.clear()
        return peaks

    @property
    def thres_abs(self):
        thres = list(self._thres_abs_buffer)
        self._thres_abs_buffer.clear()
        return thres

    @property
    def peak_values(self):
        """ get the peak values """
        values = list(self._peak_value_buffer)
        self._peak_value_buffer.clear()
        return values
        #while self._peak_value_buffer:
        #    index = self._peak_value_buffer.popleft()
        #    yield index

    @property
    def peak_indices(self):
        """ get the peak indexes """
        indexes = list(self._peak_index_buffer)
        self._peak_index_buffer.clear()
        return indexes
        #while self._peak_index_buffer:
        #    index = self._peak_index_buffer.popleft()
        #    yield index

    def clear(self):
        """ clear signal buffer and reset signal index"""
        self._signal_buffer.clear()
        self._signal_index = -1

    def findpeaks(self, value):
        """ find peaks in signal buffer,
            add peak indexes and peak values to buffer
        """
        self._signal_index += 1
        self._update_counter += 1
        signal = self._signal_buffer
        signal.append(value)
        if len(signal) != signal.maxlen:
            return
        update_thres = signal.maxlen / 2
        if self._update_counter < update_thres:
            return
        self._update_counter = 0

        peaks, raw_extremas, raw_peaks, thres_abs = \
            findpeaks2.findpeaks2(np.asarray(signal),
                                  thres=self.thres,
                                  min_dist=self.min_dist,
                                  direction=self.direction)
        index_offset = self._signal_index - signal.maxlen + 1
        chosen_index_start = int(signal.maxlen / 4)
        chosen_index_end = signal.maxlen - int(signal.maxlen / 4)
        chosen_dist = chosen_index_end - chosen_index_start

        def filter_peaks(peaks):
            ret = []
            for peak in peaks:
                peak_index = peak.peak_index
                if index_offset == 0 and peak_index < chosen_index_end:
                    ret.append(peak)
                    continue
                if peak_index < chosen_index_start or peak_index >= chosen_index_end:
                    continue
                peak_index = peak_index + index_offset
                args = list(peak)
                args[0] = peak_index
                new_peak = findpeaks2.Peak(*args)
                ret.append(new_peak)
            return ret

        def filter_extremas(peaks):
            ret = []
            for peak in peaks:
                peak_index = peak.index
                if index_offset == 0 and peak_index < chosen_index_end:
                    ret.append(peak)
                    continue
                if peak_index < chosen_index_start or peak_index >= chosen_index_end:
                    continue
                peak_index = peak_index + index_offset
                args = list(peak)
                args[0] = peak_index
                new_peak = findpeaks2.Extrema(*args)
                ret.append(new_peak)
            return ret

        # add raw peaks to buffer
        chosen_raw_extremas = filter_extremas(raw_extremas)
        self._raw_extremas_buffer.extend(chosen_raw_extremas)

        # add raw peaks to buffer
        chosen_raw_peaks = filter_peaks(raw_peaks)
        self._raw_peaks_buffer.extend(chosen_raw_peaks)

        # add peaks to buffer
        chosen_peaks = filter_peaks(peaks)
        self._peaks_buffer.extend(chosen_peaks)

        # add thres_abs to buffer
        self._thres_abs_buffer.extend([thres_abs for _ in range(chosen_dist)])

        # add peak.peak_index, peak.peak_value to buffer
        self._peak_index_buffer.extend(list(map(lambda peak:peak.peak_index, chosen_peaks)))
        self._peak_value_buffer.extend(list(map(lambda peak:peak.peak_value, chosen_peaks)))

        """
        # add peak_index from peaks to peak_index_buffer
        ipks = list(map(lambda peak:peak.peak_index, peaks))
        for ipk in ipks:
            if ipk_offset == 0 and ipk < chosen_ipks_end:
                self._peak_index_buffer.append(ipk)
                self._peak_value_buffer.append(signal[ipk])
                continue
            if ipk < chosen_ipks_start or ipk >= chosen_ipks_end:
                continue
            self._peak_index_buffer.append(ipk + ipk_offset)
            self._peak_value_buffer.append(signal[ipk])
        """

def main(argv):
    del argv
    from scipy.io import loadmat
    filepath = "/home/guo/data/sleepstage/database/ISRUC_Sleep/" + \
        "ExtractedChannels/subgroupI-Extractedchannels/subject1.mat"
    data = loadmat(filepath)
    ecg = data["X2"].flatten()[:5000]
    autopeaks = AutoPeaks(thres=0.65, min_dist=80, fs=200)
    list(map(autopeaks.findpeaks, ecg))
    ipks = autopeaks.peak_indexes
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(ecg)
    ax.plot(ipks, ecg[ipks], 'o')
    plt.show()

if __name__ == "__main__":
    app.run(main)

max_stars_repo_name: qixiuai/BCGHeart

max_issues_repo_name: qixiuai/BCGHeart

max_forks_repo_name: qixiuai/BCGHeart

reason: import numpy,from scipy


text:
import cv2
import numpy as np
from __utils__.general import show_image

def scaling(image, scale=1, dst=None):
    """
    :param image: Gray image or RGB image
    :param scale: The scale of image.
    :param dst: The image after scaling
    :return:
    """
    # Converting RGB to Gray image
    gray = image if len(image.shape) == 2 else cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)
    height, width = gray.shape

    # Get the size of image after scaling
    dst_height = int(height * scale)
    dst_width = int(width * scale)

    # Init the 2D matrix
    out = np.zeros(shape=(dst_height, dst_width), dtype=np.uint8)

    for x in range(0, dst_height - 1):
        for y in range(0, dst_width - 1):
            x_dst = x + 0.5
            y_dst = y + 0.5
            float_x = x_dst / scale
            float_y = y_dst / scale
            x_source = int(float_x + 0.5)
            y_source = int(float_y + 0.5)
            out[x][y] = gray[x_source][y_source]
    return out

def opencv_scaling(image, scale=1):
    height, width = image.shape[:2]
    new_height = scale * height
    new_width = scale * width
    res = cv2.resize(image, (new_width, new_height), interpolation=cv2.INTER_CUBIC)
    return res

if __name__ == '__main__':
    image_path = '../../asserts/images/bird.png'
    image = cv2.imread(image_path)

    # scaled_image = scaling(image, 4)
    scaled_image = opencv_scaling(image, 3)

    show_image(scaled_image)

max_stars_repo_name: jerry-le/computer-vision

max_issues_repo_name: jerry-le/computer-vision

max_forks_repo_name: jerry-le/computer-vision

reason: import numpy


text:
import subprocess as sp
from pathlib import Path
import json
import numpy as np

def _read_info(path):
    stdout_data = sp.check_output([
        'ffprobe', "-loglevel", "panic",
        str(path), '-print_format', 'json', '-show_format', '-show_streams'
    ])
    return json.loads(stdout_data.decode('utf-8'))

def _read_audio(path,
                filename,
                seek_time=None,
                duration=None,
                sample_rate=None,
                channels=None):
    if duration is None:
        target_size = None
        query_duration = None
    else:
        target_size = int(sample_rate * duration)
        query_duration = float((target_size + 1) / sample_rate)

    command = ['ffmpeg', '-y']
    command += ['-loglevel', 'panic']
    if seek_time:
        command += ['-ss', str(seek_time)]
    command += ['-i', path]

    if query_duration is not None:
        command += ['-t', str(query_duration)]
    command += ['-f', 'f32le']
    if sample_rate is not None:
        command += ['-ar', str(sample_rate)]
    command += ['-ac', '2']
    command += [filename]

    sp.run(command, check=True)
    wav = np.fromfile(filename, dtype=np.float32)

class AudioFile:
    def __init__(self, perf, details):
        self.perf = perf
        self.vocal = details['vocal']
        self.back = details['background']
        self.mix = details['original_mix']

        self.original_mean = details['original_mix_mean']
        self.original_std = details['original_mix_std']

        self.mean = details['mix_mean']
        self.std = details['mix_std']

        self.duration = details['duration']
        self.vocal_start = details['vocal_start']
        self.background_start = details['background_start']

        self.scaler = details['scaler']

max_stars_repo_name: groadabike/DAMP-VSEP-Singles

max_issues_repo_name: groadabike/DAMP-VSEP-Singles

max_forks_repo_name: groadabike/DAMP-VSEP-Singles


text:
"""Optimizing input init points to object surface."""

import os
import pdb
import tqdm
import argparse
import numpy as np
import torch
import torch.nn.functional as F

from src import config

from defense import SORDefense
from defense import repulsion_loss

def str2bool(v):
    return v.lower() in ("yes", "true", "t", "1")

parser = argparse.ArgumentParser(
    description='Extract meshes from occupancy process.'
)
parser.add_argument('--config', type=str,
                    default='configs/convonet_3plane_mn40.yaml',
                    help='Path to config file.')
parser.add_argument('--sample_npoint', type=int, default=1024,
                    help='Re-sample points number per mesh.')
parser.add_argument('--padding_scale', type=float, default=0.9,
                    help='Used in pre-processing point clouds,'
                         'padding in unit cube')
parser.add_argument('--data_root', type=str, default='',
                    help='Path to point cloud npz file.')
parser.add_argument('--train', type=str2bool, default=False,
                    help='whether defend training data')

parser.add_argument('--init_sigma', type=float, default=0.01,
                    help='sigma for normal dist used in ori_init')
parser.add_argument('--iterations', type=int, default=200,
                    help='Optimization iterations.')
parser.add_argument('--batch_size', type=int, default=192,
                    help='Batch process points.')
parser.add_argument('--lr', type=float, default=0.001,
                    help='lr in optimization')
parser.add_argument('--rep_weight', type=float, default=500.,
                    help='loss weight for repulsion term')

parser.add_argument('--sor', type=str2bool, default=True,
                    help='whether use SOR before reconstruction')
parser.add_argument('--sor_k', type=int, default=2,
                    help='KNN in SOR')
parser.add_argument('--sor_alpha', type=float, default=1.1,
                    help='Threshold = mean + alpha * std')

args = parser.parse_args()
cfg = config.load_config(args.config, 'configs/default.yaml')

device = torch.device("cuda")

args.threshold = cfg['test']['threshold']
args.input_npoint = cfg['data']['pointcloud_n']

# Model
model = config.get_model(cfg, device=device, dataset=None)
model.load_state_dict(torch.load(cfg['test']['model_file']))

# Generator
generator = config.get_generator(model, cfg, device=device)

# model and generator not updated
model.eval()
for p in model.parameters():
    p.requires_grad = False

def normalize_batch_pc(points):
    """points: [batch, K, 3]"""
    centroid = torch.mean(points, dim=1)  # [batch, 3]
    points -= centroid[:, None, :]  # center, [batch, K, 3]
    dist = torch.sum(points ** 2, dim=2) ** 0.5  # [batch, K]
    max_dist = torch.max(dist, dim=1)[0]  # [batch]
    points /= max_dist[:, None, None]
    return points

def sor_process(pc):
    """Use SOR to pre-process pc.
    Inputs:
        pc: [N, K, 3]

    Returns list of [K_i, 3]
    """
    N = len(pc)
    batch_size = 32
    sor_pc = []
    sor_defense = SORDefense(k=args.sor_k, alpha=args.sor_alpha)
    for i in range(0, N, batch_size):
        input_pc = pc[i:i + batch_size]  # [B, K, 3]
        input_pc = torch.from_numpy(input_pc).float().cuda()
        output_pc = sor_defense(input_pc)
        # to np array list
        output_pc = [
            one_pc.detach().cpu().numpy().
            astype(np.float32) for one_pc in output_pc
        ]
        sor_pc.append(output_pc)
    pc = []
    for i in range(len(sor_pc)):
        pc += sor_pc[i]  # list of [k, 3]
    assert len(pc[0].shape) == 2 and pc[0].shape[1] == 3
    return pc

def preprocess_pc(pc, num_points=None, padding_scale=1.):
    """Center and scale to be within unit cube.
    Inputs:
        pc: np.array of [K, 3]
        num_points: pick a subset of points as OccNet input.
        padding_scale: padding ratio in unit cube.
    """
    # normalize into unit cube
    center = np.mean(pc, axis=0)  # [3]
    centered_pc = pc - center
    max_dim = np.max(centered_pc, axis=0)  # [3]
    min_dim = np.min(centered_pc, axis=0)  # [3]
    scale = (max_dim - min_dim).max()
    scaled_centered_pc = centered_pc / scale * padding_scale

    # select a subset as ONet input
    if num_points is not None and \
            scaled_centered_pc.shape[0] > num_points:
        idx = np.random.choice(
            scaled_centered_pc.shape[0], num_points,
            replace=False)
        pc = scaled_centered_pc[idx]
    else:
        pc = scaled_centered_pc

    # to torch tensor
    # torch_pc is ONet input, have fixed number of points
    # torch_all_pc is for initializing the defense point cloud
    torch_pc = torch.from_numpy(pc).\
        float().cuda().unsqueeze(0)
    torch_all_pc = torch.from_numpy(scaled_centered_pc).\
        float().cuda().unsqueeze(0)
    return torch_all_pc, torch_pc

def init_points(pc):
    """Initialize points to be optimized.

    Args:
        pc (tensor): input (adv) pc, [B, N, 3]
    """
    with torch.no_grad():
        B = len(pc)
        # init points from ori_points with noise
        # if not enough points in input pc, randomly duplicate
        # else select a subset
        if isinstance(pc, list):  # after SOR
            idx = [
                torch.randint(
                    0, len(one_pc),
                    (args.sample_npoint,)).long().cuda() for one_pc in pc
            ]
        else:
            idx = torch.randint(
                0, pc.shape[1], (B, args.sample_npoint)).long().cuda()
        points = torch.stack([
            pc[i][idx[i]] for i in range(B)
        ], dim=0).float().cuda()

        # add noise
        noise = torch.randn_like(points) * args.init_sigma
        points = torch.clamp(
            points + noise,
            min=-0.5 * args.padding_scale,
            max=0.5 * args.padding_scale)
    return points

def optimize_points(opt_points, z, c,
                    rep_weight=1.,
                    iterations=1000,
                    printing=False):
    """Optimization process on point coordinates.

    Args:
        opt_points (tensor): input init points to be optimized
        z (tensor): latent code
        c (tensor): feature vector
        iterations (int, optional): opt iter. Defaults to 1000.
        printing (bool, optional): print info. Defaults to False.
    """
    # 2 losses in total
    # Geo-aware loss enforces occ_value = occ_threshold by BCE
    # Dist-aware loss pushes points uniform by repulsion loss
    opt_points = opt_points.float().cuda()
    opt_points.requires_grad_()
    B, K = opt_points.shape[:2]

    # GT occ for surface
    with torch.no_grad():
        occ_threshold = torch.ones(
            (B, K)).float().cuda() * args.threshold

    opt = torch.optim.Adam([opt_points], lr=args.lr)

    # start optimization
    for i in range(iterations + 1):
        # 1. occ = threshold
        occ_value = generator.model.decode(opt_points, c).logits
        occ_loss = F.binary_cross_entropy_with_logits(
            occ_value, occ_threshold, reduction='none')  # [B, K]
        occ_loss = torch.mean(occ_loss)
        occ_loss = occ_loss * K

        # 2. repulsion loss
        rep_loss = torch.tensor(0.).float().cuda()
        if rep_weight > 0.:
            rep_loss = repulsion_loss(opt_points)  # [B]
            rep_loss = torch.mean(rep_loss)
            rep_loss = rep_loss * rep_weight

        loss = occ_loss + rep_loss
        opt.zero_grad()
        loss.backward()
        opt.step()
        if printing and i % 100 == 0:
            print('iter {}, loss {:.4f}'.format(i, loss.item()))
            print('occ loss: {:.4f}, '
                  'rep loss: {:.4f}\n'
                  'occ value mean: {:.4f}'.
                  format(occ_loss.item(),
                         rep_loss.item(),
                         torch.sigmoid(occ_value).mean().item()))
    opt_points.detach_()
    opt_points = normalize_batch_pc(opt_points)
    return opt_points.detach().cpu().numpy()

def get_save_name(path, train=args.train):
    """Saving name containing model settings."""
    sub_path = path.split('/')
    save_name = sub_path[-1]
    save_name = 'convonet_opt-' + save_name
    save_folder = path[:path.rindex(sub_path[-1])]
    save_folder = os.path.join(save_folder, 'ConvONet-Opt')
    if not os.path.exists(save_folder):
        os.makedirs(save_folder)
    save_path = os.path.join(save_folder, save_name)
    return save_path

def defend_point_cloud(pc):
    """Apply defense to input point clouds.

    Args:
        pc (tensor): [num_data, K, 3]
    """
    # possible SOR preprocessor
    if args.sor:
        with torch.no_grad():
            pc = sor_process(pc)
    torch.cuda.empty_cache()

    opt_pc = np.zeros(
        (len(pc), args.sample_npoint, 3),
        dtype=np.float32)

    # batch process
    for idx in tqdm.trange(0, len(pc), args.batch_size):
        # prepare for input
        with torch.no_grad():
            batch_pc = pc[idx:idx + args.batch_size]  # [B, K, 3]
            # preprocess
            batch_proc_pc = [
                preprocess_pc(
                    one_pc, num_points=args.input_npoint,
                    padding_scale=args.padding_scale) for one_pc in batch_pc
            ]
            # the selected input_n points from batch_pc after preprocess
            # sel_pc are for ONet input and have fixed number of points
            batch_proc_sel_pc = torch.cat([
                one_pc[1] for one_pc in batch_proc_pc
            ], dim=0).float().cuda()
            # proc_pc may have different num_points because of SOR
            # they're used for initializing the defense point clouds
            try:
                batch_proc_pc = torch.cat([
                    one_pc[0] for one_pc in batch_proc_pc
                ], dim=0).float().cuda()
            except RuntimeError:
                batch_proc_pc = [
                    one_pc[0][0] for one_pc in batch_proc_pc
                ]  # list of [num, 3]

            # get latent feature vector c
            # c is [B, c_dim (typically 512)]
            c = generator.model.encode_inputs(batch_proc_sel_pc)

            # z is of no use
            z = None

        # init points and optimize
        points = init_points(batch_proc_pc)
        points.requires_grad_()
        points = optimize_points(points, z, c,
                                 rep_weight=args.rep_weight,
                                 iterations=args.iterations,
                                 printing=True)
        opt_pc[idx:idx + args.batch_size] = points

    return opt_pc

def defend_npz_test_data(path=args.data_root):
    '''Apply defense to all pc in a npz file.
    Then save it for evaluation.
    '''
    # load data
    npz = np.load(path)
    test_pc = npz['test_pc'][..., :3]  # [..., K, 3]
    test_label = npz['test_label']
    try:
        target_label = npz['target_label']
    except KeyError:
        target_label = None

    # defend
    def_test_pc = defend_point_cloud(test_pc)

    # save new npz file
    save_path = get_save_name(path)
    if target_label is None:
        np.savez(save_path,
                 test_pc=def_test_pc.astype(np.float32),
                 test_label=test_label.astype(np.uint8))
    else:
        np.savez(save_path,
                 test_pc=def_test_pc.astype(np.float32),
                 test_label=test_label.astype(np.uint8),
                 target_label=target_label.astype(np.uint8))
    print('defense result saved to {}'.format(save_path))

def defend_npz_train_test_data(path=args.data_root):
    '''Apply defense to all pc in a npz file.
    Then save it for evaluation.
    '''
    # load data
    npz = np.load(path)
    train_pc = npz['train_pc'][..., :3]  # [num_data, K, 3]
    test_pc = npz['test_pc'][..., :3]
    train_label = npz['train_label']
    test_label = npz['test_label']

    # defend
    def_train_pc = defend_point_cloud(train_pc)
    def_test_pc = defend_point_cloud(test_pc)

    # save new npz file
    save_path = get_save_name(path)
    np.savez(save_path,
             train_pc=def_train_pc.astype(np.float32),
             train_label=train_label.astype(np.uint8),
             test_pc=def_test_pc.astype(np.float32),
             test_label=test_label.astype(np.uint8))
    print('defense result saved to {}'.format(save_path))

if __name__ == '__main__':
    if args.train:
        # defend both train and test data
        # for hybrid training usage
        defend_npz_train_test_data(args.data_root)
    else:
        # if given a dir, then defend all the files in it
        data_root = args.data_root
        if os.path.isdir(data_root):
            all_files = os.listdir(data_root)
            for file in all_files:
                one_file = os.path.join(data_root, file)
                if os.path.isfile(one_file):
                    defend_npz_test_data(one_file)
        else:
            defend_npz_test_data(args.data_root)

max_stars_repo_name: code-roamer/IF-Defense

max_issues_repo_name: code-roamer/IF-Defense


text:
# Code based on an original given in
# https://github.com/agaidus/census_data_extraction/blob/master/census_mapper.py

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point
from numpy.random import RandomState, uniform

def generate_random_points_in_polygon(poly, num_points, seed=None):
	"""
	Create a list of randomly generated points within a polygon.

	Parameters
	----------
	poly : Polygon
	num_points : int
		The number of random points to create within the polygon
	seed : int, optional
		A random seed

	Returns
	-------
	List
	"""
	min_x, min_y, max_x, max_y = poly.bounds
	points = []
	i = 0
	while len(points) < num_points:
		s = RandomState(seed + i) if seed else RandomState(seed)
		random_point = Point([s.uniform(min_x, max_x), s.uniform(min_y, max_y)])
		if random_point.within(poly):
			points.append(random_point)
		i += 1
	return points

def generate_points_in_areas(gdf, values, points_per_unit=1, seed=None):
	"""
	Create a GeoSeries of random points in polygons.

	Parameters
	----------
	gdf : GeoDataFrame
		The areas in which to create points
	values : str or Series
		The [possibly scaled] number of points to create in each area
	points_per_unit : numeric, optional
		The rate to scale the values in point generation.
	seed : int, optional
		A random seed

	Returns
	-------
	GeoSeries
	"""
	geometry = gdf.geometry
	if isinstance(values, str) and values in gdf.columns:
		values = gdf[values]
	new_values = (values / points_per_unit).astype(int)
	g = gpd.GeoDataFrame(data={'vals': new_values}, geometry=geometry)
	a = g.apply(lambda row: tuple(generate_random_points_in_polygon(row['geometry'], row['vals'], seed)), 1)
	b = gpd.GeoSeries(a.apply(pd.Series).stack(), crs=geometry.crs)
	b.name = 'geometry'
	return b

reason: from numpy


text:
from __future__ import absolute_import
import re
from copy import copy

import numpy as nm

from sfepy.base.base import (as_float_or_complex, get_default, assert_,
                             Container, Struct, basestr, goptions)
from sfepy.base.compat import in1d

# Used for imports in term files.
from sfepy.terms.extmods import terms
import six
from six.moves import range
from functools import reduce

_match_args = re.compile('^([^\(\}]*)\((.*)\)$').match
_match_virtual = re.compile('^virtual$').match
_match_state = re.compile('^state(_[_a-zA-Z0-9]+)?$').match
_match_parameter = re.compile('^parameter(_[_a-zA-Z0-9]+)?$').match
_match_material = re.compile('^material(_[_a-zA-Z0-9]+)?$').match
_match_material_opt = re.compile('^opt_material(_[_a-zA-Z0-9]+)?$').match
_match_material_root = re.compile('(.+)\.(.*)').match
_match_ts = re.compile('^ts$').match

def get_arg_kinds(arg_types):
    """
    Translate `arg_types` of a Term to a canonical form.

    Parameters
    ----------
    arg_types : tuple of strings
        The term argument types, as given in the `arg_types` attribute.

    Returns
    -------
    arg_kinds : list of strings
        The argument kinds - one of 'virtual_variable', 'state_variable',
        'parameter_variable', 'opt_material', 'ts', 'user'.
    """
    arg_kinds = []
    for ii, arg_type in enumerate(arg_types):
        if _match_virtual(arg_type):
            arg_kinds.append('virtual_variable')

        elif _match_state(arg_type):
            arg_kinds.append('state_variable')

        elif _match_parameter(arg_type):
            arg_kinds.append('parameter_variable')

        elif _match_material(arg_type):
            arg_kinds.append('material')

        elif _match_material_opt(arg_type):
            arg_kinds.append('opt_material')
            if ii > 0:
                msg = 'opt_material at position %d, must be at 0!' % ii
                raise ValueError(msg)

        elif _match_ts(arg_type):
            arg_kinds.append('ts')

        else:
            arg_kinds.append('user')

    return arg_kinds

def get_shape_kind(integration):
    """
    Get data shape kind for given integration type.
    """
    if integration == 'surface':
        shape_kind = 'surface'

    elif integration in ('volume', 'plate', 'surface_extra'):
        shape_kind = 'volume'

    elif integration == 'point':
        shape_kind = 'point'

    else:
        raise NotImplementedError('unsupported term integration! (%s)'
                                  % integration)

    return shape_kind

def split_complex_args(args):
    """
    Split complex arguments to real and imaginary parts.

    Returns
    -------
    newargs : dictionary
        Dictionary with lists corresponding to `args` such that each
        argument of numpy.complex128 data type is split to its real and
        imaginary part. The output depends on the number of complex
        arguments in 'args':

          - 0: list (key 'r') identical to input one

          - 1: two lists with keys 'r', 'i' corresponding to real
            and imaginary parts

          - 2: output dictionary contains four lists:

            - 'r' - real(arg1), real(arg2)
            - 'i' - imag(arg1), imag(arg2)
            - 'ri' - real(arg1), imag(arg2)
            - 'ir' - imag(arg1), real(arg2)
    """
    newargs = {}
    cai = []

    for ii, arg in enumerate(args):
        if isinstance(arg, nm.ndarray) and (arg.dtype == nm.complex128):
            cai.append(ii)

    if len(cai) > 0:
        newargs['r'] = list(args[:])
        newargs['i'] = list(args[:])

        arg1 = cai[0]
        newargs['r'][arg1] = args[arg1].real.copy()
        newargs['i'][arg1] = args[arg1].imag.copy()

        if len(cai) == 2:
            arg2 = cai[1]
            newargs['r'][arg2] = args[arg2].real.copy()
            newargs['i'][arg2] = args[arg2].imag.copy()

            newargs['ri'] = list(args[:])
            newargs['ir'] = list(args[:])
            newargs['ri'][arg1] = newargs['r'][arg1]
            newargs['ri'][arg2] = newargs['i'][arg2]
            newargs['ir'][arg1] = newargs['i'][arg1]
            newargs['ir'][arg2] = newargs['r'][arg2]

        elif len(cai) > 2:
            raise NotImplementedError('more than 2 complex arguments! (%d)'
                                      % len(cai))

    else:
        newargs['r'] = args[:]

    return newargs

def create_arg_parser():
    from pyparsing import Literal, Word, delimitedList, Group, \
         StringStart, StringEnd, Optional, nums, alphas, alphanums

    inumber = Word("+-" + nums, nums)

    history = Optional(Literal('[').suppress() + inumber
                       + Literal(']').suppress(), default=0)("history")
    history.setParseAction(lambda str, loc, toks: int(toks[0]))

    variable = Group(Word(alphas, alphanums + '._') + history)

    derivative = Group(Literal('d') + variable\
                       + Literal('/').suppress() + Literal('dt'))

    trace = Group(Literal('tr') + Literal('(').suppress() + variable \
                  + Literal(')').suppress())

    generalized_var = derivative | trace | variable

    args = StringStart() + delimitedList(generalized_var) + StringEnd()

    return args

class ConnInfo(Struct):

    def get_region(self, can_trace=True):
        if self.is_trace and can_trace:
            return self.region.get_mirror_region()
        else:
            return self.region

    def get_region_name(self, can_trace=True):
        if self.is_trace and can_trace:
            reg = self.region.get_mirror_region()
        else:
            reg = self.region

        if reg is not None:
            return reg.name
        else:
            return None

class Terms(Container):

    @staticmethod
    def from_desc(term_descs, regions, integrals=None):
        """
        Create terms, assign each term its region.
        """
        from sfepy.terms import term_table

        terms = Terms()
        for td in term_descs:
            try:
                constructor = term_table[td.name]
            except:
                msg = "term '%s' is not in %s" % (td.name,
                                                  sorted(term_table.keys()))
                raise ValueError(msg)

            try:
                region = regions[td.region]
            except IndexError:
                raise KeyError('region "%s" does not exist!' % td.region)

            term = Term.from_desc(constructor, td, region, integrals=integrals)
            terms.append(term)

        return terms

    def __init__(self, objs=None):
        Container.__init__(self, objs=objs)

        self.update_expression()

    def insert(self, ii, obj):
        Container.insert(self, ii, obj)
        self.update_expression()

    def append(self, obj):
        Container.append(self, obj)
        self.update_expression()

    def update_expression(self):
        self.expression = []
        for term in self:
            aux = [term.sign, term.name, term.arg_str,
                   term.integral_name, term.region.name]
            self.expression.append(aux)

    def __mul__(self, other):
        out = Terms()
        for name, term in self.iteritems():
            out.append(term * other)

        return out

    def __rmul__(self, other):
        return self * other

    def __add__(self, other):
        if isinstance(other, Term):
            out = self.copy()
            out.append(other)

        elif isinstance(other, Terms):
            out = Terms(self._objs + other._objs)

        else:
            raise ValueError('cannot add Terms with %s!' % other)

        return out

    def __radd__(self, other):
        return self + other

    def __sub__(self, other):
        if isinstance(other, Term):
            out = self + (-other)

        elif isinstance(other, Terms):
            out = self + (-other)

        else:
            raise ValueError('cannot subtract Terms with %s!' % other)

        return out

    def __rsub__(self, other):
        return -self + other

    def __pos__(self):
        return self

    def __neg__(self):
        return -1.0 * self

    def setup(self):
        for term in self:
            term.setup()

    def assign_args(self, variables, materials, user=None):
        """
        Assign all term arguments.
        """
        for term in self:
            term.assign_args(variables, materials, user)

    def get_variable_names(self):
        out = []
        for term in self:
            out.extend(term.get_variable_names())
        return list(set(out))

    def get_material_names(self):
        out = []
        for term in self:
            out.extend(term.get_material_names())
        return list(set(out))

    def get_user_names(self):
        out = []
        for term in self:
            out.extend(term.get_user_names())
        return list(set(out))

class Term(Struct):
    name = ''
    arg_types = ()
    arg_shapes = {}
    integration = 'volume'
    geometries = ['1_2', '2_3', '2_4', '3_4', '3_8']

    @staticmethod
    def new(name, integral, region, **kwargs):
        from sfepy.terms import term_table

        arg_str = _match_args(name)
        if arg_str is not None:
            name, arg_str = arg_str.groups()

        else:
            raise ValueError('bad term syntax! (%s)' % name)

        if name in term_table:
            constructor = term_table[name]

        else:
            msg = "term '%s' is not in %s" % (name, sorted(term_table.keys()))
            raise ValueError(msg)

        obj = constructor(name, arg_str, integral, region, **kwargs)
        return obj

    @staticmethod
    def from_desc(constructor, desc, region, integrals=None):
        from sfepy.discrete import Integrals

        if integrals is None:
            integrals = Integrals()

        integral = integrals.get(desc.integral)
        obj = constructor(desc.name, desc.args, integral, region)
        obj.sign = desc.sign

        return obj

    def __init__(self, name, arg_str, integral, region, **kwargs):
        self.name = name
        self.arg_str = arg_str
        self.region = region
        self._kwargs = kwargs
        self._integration = self.integration
        self.sign = 1.0

        self.set_integral(integral)

    def __mul__(self, other):
        try:
            mul = as_float_or_complex(other)

        except ValueError:
            raise ValueError('cannot multiply Term with %s!' % other)

        out = self.copy(name=self.name)

        out.sign = mul * self.sign

        return out

    def __rmul__(self, other):
        return self * other

    def __add__(self, other):
        if isinstance(other, Term):
            out = Terms([self, other])

        else:
            out = NotImplemented

        return out

    def __sub__(self, other):
        if isinstance(other, Term):
            out = Terms([self, -1.0 * other])

        else:
            out = NotImplemented

        return out

    def __pos__(self):
        return self

    def __neg__(self):
        out = -1.0 * self
        return out

    def get_str(self):
        return ('%+.2e * %s.%d.%s(%s)'
                % (self.sign, self.name, self.integral.order,
                   self.region.name, self.arg_str))

    def set_integral(self, integral):
        """
        Set the term integral.
        """
        self.integral = integral
        if self.integral is not None:
            self.integral_name = self.integral.name

    def setup(self):
        self.function = Struct.get(self, 'function', None)

        self.step = 0
        self.dt = 1.0
        self.is_quasistatic = False
        self.has_region = True

        self.setup_formal_args()

        if self._kwargs:
            self.setup_args(**self._kwargs)

        else:
            self.args = []

    def setup_formal_args(self):
        self.arg_names = []
        self.arg_steps = {}
        self.arg_derivatives = {}
        self.arg_traces = {}

        parser = create_arg_parser()
        self.arg_desc = parser.parseString(self.arg_str)

        for arg in self.arg_desc:
            trace = False
            derivative = None

            if isinstance(arg[1], int):
                name, step = arg

            else:
                kind = arg[0]
                name, step = arg[1]
                if kind == 'd':
                    derivative = arg[2]
                elif kind == 'tr':
                    trace = True

            match = _match_material_root(name)
            if match:
                name = (match.group(1), match.group(2))

            self.arg_names.append(name)
            self.arg_steps[name] = step
            self.arg_derivatives[name] = derivative
            self.arg_traces[name] = trace

    def setup_args(self, **kwargs):
        self._kwargs = kwargs

        self.args = []
        for arg_name in self.arg_names:
            if isinstance(arg_name, basestr):
                self.args.append(self._kwargs[arg_name])

            else:
                self.args.append((self._kwargs[arg_name[0]], arg_name[1]))

        self.classify_args()
        self.check_args()

    def assign_args(self, variables, materials, user=None):
        """
        Check term argument existence in variables, materials, user data
        and assign the arguments to terms. Also check compatibility of
        field and term regions.
        """
        if user is None:
            user = {}

        user.setdefault('ts', Struct())

        kwargs = {}
        for arg_name in self.arg_names:
            if isinstance(arg_name, basestr):
                if arg_name in variables.names:
                    kwargs[arg_name] = variables[arg_name]

                elif arg_name in user:
                    kwargs[arg_name] = user[arg_name]

                else:
                    raise ValueError('argument %s not found!' % arg_name)

            else:
                arg_name = arg_name[0]
                if arg_name in materials.names:
                    kwargs[arg_name] = materials[arg_name]

                else:
                    raise ValueError('material argument %s not found!'
                                     % arg_name)

        self.setup_args(**kwargs)

    def classify_args(self):
        """
        Classify types of the term arguments and find matching call
        signature.

        A state variable can be in place of a parameter variable and
        vice versa.
        """
        self.names = Struct(name='arg_names',
                            material=[], variable=[], user=[],
                            state=[], virtual=[], parameter=[])

        # Prepare for 'opt_material' - just prepend a None argument if needed.
        if isinstance(self.arg_types[0], tuple):
            arg_types = self.arg_types[0]

        else:
            arg_types = self.arg_types

        if len(arg_types) == (len(self.args) + 1):
            self.args.insert(0, (None, None))
            self.arg_names.insert(0, (None, None))

        if isinstance(self.arg_types[0], tuple):
            assert_(len(self.modes) == len(self.arg_types))
            # Find matching call signature using variable arguments - material
            # and user arguments are ignored!
            matched = []
            for it, arg_types in enumerate(self.arg_types):
                arg_kinds = get_arg_kinds(arg_types)
                if self._check_variables(arg_kinds):
                    matched.append((it, arg_kinds))

            if len(matched) == 1:
                i_match, arg_kinds = matched[0]
                arg_types = self.arg_types[i_match]
                self.mode = self.modes[i_match]

            elif len(matched) == 0:
                msg = 'cannot match arguments! (%s)' % self.arg_names
                raise ValueError(msg)

            else:
                msg = 'ambiguous arguments! (%s)' % self.arg_names
                raise ValueError(msg)

        else:
            arg_types = self.arg_types
            arg_kinds = get_arg_kinds(self.arg_types)
            self.mode = Struct.get(self, 'mode', None)

            if not self._check_variables(arg_kinds):
                raise ValueError('cannot match variables! (%s)'
                                 % self.arg_names)

        # Set actual argument types.
        self.ats = list(arg_types)

        for ii, arg_kind in enumerate(arg_kinds):
            name = self.arg_names[ii]
            if arg_kind.endswith('variable'):
                names = self.names.variable

                if arg_kind == 'virtual_variable':
                    self.names.virtual.append(name)

                elif arg_kind == 'state_variable':
                    self.names.state.append(name)

                elif arg_kind == 'parameter_variable':
                    self.names.parameter.append(name)

            elif arg_kind.endswith('material'):
                names = self.names.material

            else:
                names = self.names.user

            names.append(name)

        self.n_virtual = len(self.names.virtual)
        if self.n_virtual > 1:
            raise ValueError('at most one virtual variable is allowed! (%d)'
                             % self.n_virtual)

        self.set_arg_types()

        self.setup_integration()

    def _check_variables(self, arg_kinds):
        for ii, arg_kind in enumerate(arg_kinds):
            if arg_kind.endswith('variable'):
                var = self.args[ii]
                check = {'virtual_variable' : var.is_virtual,
                         'state_variable' : var.is_state_or_parameter,
                         'parameter_variable' : var.is_state_or_parameter}
                if not check[arg_kind]():
                    return False

        else:
            return True

    def set_arg_types(self):
        pass

    def check_args(self):
        """
        Common checking to all terms.

        Check compatibility of field and term regions.
        """
        vns = self.get_variable_names()
        for name in vns:
            field = self._kwargs[name].get_field()
            if field is None:
                continue

            if not nm.all(in1d(self.region.vertices,
                               field.region.vertices)):
                msg = ('%s: incompatible regions: (self, field %s)'
                       + '(%s in %s)') %\
                       (self.name, field.name,
                        self.region.vertices, field.region.vertices)
                raise ValueError(msg)

    def get_variable_names(self):
        return self.names.variable

    def get_material_names(self):
        out = []
        for aux in self.names.material:
            if aux[0] is not None:
                out.append(aux[0])
        return out

    def get_user_names(self):
        return self.names.user

    def get_virtual_name(self):
        if not self.names.virtual:
            return None

        var = self.get_virtual_variable()
        return var.name

    def get_state_names(self):
        """
        If variables are given, return only true unknowns whose data are of
        the current time step (0).
        """
        variables = self.get_state_variables()
        return [var.name for var in variables]

    def get_parameter_names(self):
        return copy(self.names.parameter)

    def get_conn_key(self):
        """The key to be used in DOF connectivity information."""
        key = (self.name,) + tuple(self.arg_names)
        key += (self.integral_name, self.region.name)

        return key

    def get_conn_info(self):
        vvar = self.get_virtual_variable()
        svars = self.get_state_variables()
        pvars = self.get_parameter_variables()

        all_vars = self.get_variables()

        dc_type = self.get_dof_conn_type()
        tgs = self.get_geometry_types()

        v_tg = None
        if vvar is not None:
            field = vvar.get_field()
            if field is not None:
                if vvar.name in tgs:
                    v_tg = tgs[vvar.name]

                else:
                    v_tg = None

        else:
            # No virtual variable -> all unknowns are in fact known parameters.
            pvars += svars
            svars = []

        region = self.get_region()
        if region is not None:
            is_any_trace = reduce(lambda x, y: x or y,
                                  list(self.arg_traces.values()))
            if is_any_trace:
                region.setup_mirror_region()

        vals = []
        aux_pvars = []
        for svar in svars:
            # Allow only true state variables.
            if not svar.is_state():
                aux_pvars.append(svar)
                continue

            field = svar.get_field()
            is_trace = self.arg_traces[svar.name]

            if svar.name in tgs:
                ps_tg = tgs[svar.name]
            else:
                ps_tg = v_tg

            val = ConnInfo(virtual=vvar,
                           state=svar,
                           primary=svar,
                           has_virtual=True,
                           has_state=True,
                           is_trace=is_trace,
                           dc_type=dc_type,
                           v_tg=v_tg,
                           ps_tg=ps_tg,
                           region=region,
                           all_vars=all_vars)
            vals.append(val)

        pvars += aux_pvars
        for pvar in pvars:
            field = pvar.get_field()
            is_trace = self.arg_traces[pvar.name]

            if pvar.name in tgs:
                ps_tg = tgs[pvar.name]
            else:
                ps_tg = v_tg

            val = ConnInfo(virtual=vvar,
                           state=None,
                           primary=pvar.get_primary(),
                           has_virtual=vvar is not None,
                           has_state=False,
                           is_trace=is_trace,
                           dc_type=dc_type,
                           v_tg=v_tg,
                           ps_tg=ps_tg,
                           region=region,
                           all_vars=all_vars)
            vals.append(val)

        if vvar and (len(vals) == 0):
            # No state, parameter variables, just the virtual one.
            val = ConnInfo(virtual=vvar,
                           state=vvar.get_primary(),
                           primary=vvar.get_primary(),
                           has_virtual=True,
                           has_state=False,
                           is_trace=False,
                           dc_type=dc_type,
                           v_tg=v_tg,
                           ps_tg=v_tg,
                           region=region,
                           all_vars=all_vars)
            vals.append(val)

        return vals

    def get_args_by_name(self, arg_names):
        """
        Return arguments by name.
        """
        out = []
        for name in arg_names:
            try:
                ii = self.arg_names.index(name)

            except ValueError:
                raise ValueError('non-existing argument! (%s)' % name)

            out.append(self.args[ii])

        return out

    def get_args(self, arg_types=None, **kwargs):
        """
        Return arguments by type as specified in arg_types (or
        self.ats). Arguments in **kwargs can override the ones assigned
        at the term construction - this is useful for passing user data.
        """
        ats = self.ats
        if arg_types is None:
            arg_types = ats
        args = []

        region_name, iorder = self.region.name, self.integral.order
        for at in arg_types:
            ii = ats.index(at)
            arg_name = self.arg_names[ii]
            if isinstance(arg_name, basestr):
                if arg_name in kwargs:
                    args.append(kwargs[arg_name])

                else:
                    args.append(self.args[ii])

            else:
                mat, par_name = self.args[ii]
                if mat is not None:
                    mat_data = mat.get_data((region_name, iorder), par_name)

                else:
                    mat_data = None

                args.append(mat_data)

        return args

    def get_kwargs(self, keys, **kwargs):
        """Extract arguments from **kwargs listed in keys (default is
        None)."""
        return [kwargs.get(name) for name in keys]

    def get_arg_name(self, arg_type, full=False, join=None):
        """
        Get the name of the argument specified by `arg_type.`

        Parameters
        ----------
        arg_type : str
            The argument type string.
        full : bool
            If True, return the full name. For example, if the name of a
            variable argument is 'u' and its time derivative is
            requested, the full name is 'du/dt'.
        join : str, optional
            Optionally, the material argument name tuple can be joined
            to a single string using the `join` string.

        Returns
        -------
        name : str
            The argument name.
        """
        try:
            ii = self.ats.index(arg_type)
        except ValueError:
            return None

        name = self.arg_names[ii]
        if full:
            # Include derivatives.
            if self.arg_derivatives[name]:
                name = 'd%s/%s' % (name, self.arg_derivatives[name])

        if (join is not None) and isinstance(name, tuple):
            name = join.join(name)

        return name

    def setup_integration(self):
        self.has_geometry = True

        self.geometry_types = {}
        if isinstance(self.integration, basestr):
            for var in self.get_variables():
                self.geometry_types[var.name] = self.integration

        else:
            if self.mode is not None:
                self.integration = self._integration[self.mode]

            if self.integration is not None:
                for arg_type, gtype in six.iteritems(self.integration):
                    var = self.get_args(arg_types=[arg_type])[0]
                    self.geometry_types[var.name] = gtype

        gtypes = list(set(self.geometry_types.values()))

        if 'surface_extra' in gtypes:
            self.dof_conn_type = 'volume'

        elif len(gtypes):
            self.dof_conn_type = gtypes[0]

    def get_region(self):
        return self.region

    def get_geometry_types(self):
        """
        Returns
        -------
        out : dict
            The required geometry types for each variable argument.
        """
        return self.geometry_types

    def get_dof_conn_type(self):
        return Struct(name='dof_conn_info', type=self.dof_conn_type,
                      region_name=self.region.name)

    def get_assembling_cells(self, shape=None):
        """
        Return the assembling cell indices into a DOF connectivity.
        """
        cells = nm.arange(shape[0], dtype=nm.int32)

        return cells

    def time_update(self, ts):
        if ts is not None:
            self.step = ts.step
            self.dt = ts.dt
            self.is_quasistatic = ts.is_quasistatic

        if 'ts' in self._kwargs:
            self._kwargs['ts'].update(ts)

    def advance(self, ts):
        """
        Advance to the next time step. Implemented in subclasses.
        """

    def get_vector(self, variable):
        """Get the vector stored in `variable` according to self.arg_steps
        and self.arg_derivatives. Supports only the backward difference w.r.t.
        time."""

        name = variable.name
        return variable(step=self.arg_steps[name],
                        derivative=self.arg_derivatives[name])

    def get_variables(self, as_list=True):

        if as_list:
            variables = self.get_args_by_name(self.names.variable)

        else:
            variables = {}
            for var in self.get_args_by_name(self.names.variable):
                variables[var.name] = var

        return variables

    def get_virtual_variable(self):
        aux = self.get_args_by_name(self.names.virtual)
        if len(aux) == 1:
            var = aux[0]

        else:
            var = None

        return var

    def get_state_variables(self, unknown_only=False):
        variables = self.get_args_by_name(self.names.state)

        if unknown_only:
            variables = [var for var in variables
                         if (var.kind == 'unknown') and
                         (self.arg_steps[var.name] == 0)]

        return variables

    def get_parameter_variables(self):
        return self.get_args_by_name(self.names.parameter)

    def get_materials(self, join=False):
        materials = self.get_args_by_name(self.names.material)

        for mat in materials:
            if mat[0] is None:
                materials.remove(mat)

        if join:
            materials = list(set(mat[0] for mat in materials))

        return materials

    def get_qp_key(self):
        """
        Return a key identifying uniquely the term quadrature points.
        """
        return (self.region.name, self.integral.order)

    def get_physical_qps(self):
        """
        Get physical quadrature points corresponding to the term region
        and integral.
        """
        from sfepy.discrete.common.mappings import get_physical_qps, PhysicalQPs

        if self.integration == 'point':
            phys_qps = PhysicalQPs()

        else:
            phys_qps = get_physical_qps(self.region, self.integral)

        return phys_qps

    def get_mapping(self, variable, get_saved=False, return_key=False):
        """
        Get the reference mapping from a variable.

        Notes
        -----
        This is a convenience wrapper of Field.get_mapping() that
        initializes the arguments using the term data.
        """
        integration = self.geometry_types[variable.name]
        is_trace = self.arg_traces[variable.name]

        if is_trace:
            region = self.region.get_mirror_region()

        else:
            region = self.region

        out = variable.field.get_mapping(region,
                                         self.integral, integration,
                                         get_saved=get_saved,
                                         return_key=return_key)

        return out

    def get_data_shape(self, variable):
        """
        Get data shape information from variable.

        Notes
        -----
        This is a convenience wrapper of FieldVariable.get_data_shape() that
        initializes the arguments using the term data.
        """
        integration = self.geometry_types[variable.name]
        is_trace = self.arg_traces[variable.name]

        if is_trace:
            region = self.region.get_mirror_region()

        else:
            region = self.region

        out = variable.get_data_shape(self.integral, integration, region.name)
        return out

    def get(self, variable, quantity_name, bf=None, integration=None,
            step=None, time_derivative=None):
        """
        Get the named quantity related to the variable.

        Notes
        -----
        This is a convenience wrapper of Variable.evaluate() that
        initializes the arguments using the term data.
        """
        name = variable.name

        step = get_default(step, self.arg_steps[name])
        time_derivative = get_default(time_derivative,
                                      self.arg_derivatives[name])
        integration = get_default(integration, self.geometry_types[name])

        data = variable.evaluate(mode=quantity_name,
                                 region=self.region, integral=self.integral,
                                 integration=integration,
                                 step=step, time_derivative=time_derivative,
                                 is_trace=self.arg_traces[name], bf=bf)
        return data

    def check_shapes(self, *args, **kwargs):
        """
        Check term argument shapes at run-time.
        """
        from sfepy.base.base import output
        from sfepy.mechanics.tensors import dim2sym

        dim = self.region.dim
        sym = dim2sym(dim)

        def _parse_scalar_shape(sh):
            if isinstance(sh, basestr):
                if sh == 'D':
                    return dim

                elif sh == 'D2':
                    return dim**2

                elif sh == 'S':
                    return sym

                elif sh == 'N': # General number.
                    return nm.inf

                elif sh == 'str':
                    return 'str'

                else:
                    return int(sh)

            else:
                return sh

        def _parse_tuple_shape(sh):
            if isinstance(sh, basestr):
                return tuple((_parse_scalar_shape(ii.strip())
                              for ii in sh.split(',')))

            else:
                return (int(sh),)

        arg_kinds = get_arg_kinds(self.ats)

        arg_shapes_list = self.arg_shapes
        if not isinstance(arg_shapes_list, list):
            arg_shapes_list = [arg_shapes_list]

        # Loop allowed shapes until a match is found, else error.
        allowed_shapes = []
        prev_shapes = {}
        actual_shapes = {}
        for _arg_shapes in arg_shapes_list:
            # Unset shapes are taken from the previous iteration.
            arg_shapes = copy(prev_shapes)
            arg_shapes.update(_arg_shapes)
            prev_shapes = arg_shapes

            allowed_shapes.append(arg_shapes)

            n_ok = 0
            for ii, arg_kind in enumerate(arg_kinds):
                if arg_kind in ('user', 'ts'):
                    n_ok += 1
                    continue

                arg = args[ii]
                key = '%s:%s' % (self.ats[ii], self.arg_names[ii])

                if self.mode is not None:
                    extended_ats = self.ats[ii] + ('/%s' % self.mode)

                else:
                    extended_ats = self.ats[ii]

                try:
                    sh = arg_shapes[self.ats[ii]]

                except KeyError:
                    sh = arg_shapes[extended_ats]

                if arg_kind.endswith('variable'):
                    n_el, n_qp, _dim, n_en, n_c = self.get_data_shape(arg)
                    actual_shapes[key] = (n_c,)
                    shape = _parse_scalar_shape(sh[0] if isinstance(sh, tuple)
                                                else sh)
                    if nm.isinf(shape):
                        n_ok += 1

                    else:
                        n_ok += shape == n_c

                elif arg_kind.endswith('material'):
                    if arg is None: # Switched-off opt_material.
                        n_ok += sh is None
                        continue

                    if sh is None:
                        continue

                    prefix = ''
                    if isinstance(sh, basestr):
                        aux = sh.split(':')
                        if len(aux) == 2:
                            prefix, sh = aux

                    if sh == 'str':
                        n_ok += isinstance(arg, basestr)
                        continue

                    shape = _parse_tuple_shape(sh)
                    ls = len(shape)

                    aarg = nm.array(arg, ndmin=1)
                    actual_shapes[key] = aarg.shape

                    # Substiture general dimension 'N' with actual value.
                    iinfs = nm.where(nm.isinf(shape))[0]
                    if len(iinfs):
                        shape = list(shape)
                        for iinf in iinfs:
                            shape[iinf] = aarg.shape[-ls+iinf]
                        shape = tuple(shape)

                    if (ls > 1) or (shape[0] > 1):
                        # Array.
                        n_ok += shape == aarg.shape[-ls:]
                        actual_shapes[key] = aarg.shape[-ls:]

                    elif (ls == 1) and (shape[0] == 1):
                        # Scalar constant.
                        from numbers import Number
                        n_ok += isinstance(arg, Number)

                else:
                    n_ok += 1

            if n_ok == len(arg_kinds):
                break

        else:
            term_str = self.get_str()
            output('allowed argument shapes for term "%s":' % term_str)
            output(allowed_shapes)
            output('actual argument shapes:')
            output(actual_shapes)
            raise ValueError('wrong arguments shapes for "%s" term! (see above)'
                             % term_str)

    def standalone_setup(self):
        from sfepy.discrete import create_adof_conns, Variables

        conn_info = {'aux' : self.get_conn_info()}
        adcs = create_adof_conns(conn_info, None)

        variables = Variables(self.get_variables())
        variables.set_adof_conns(adcs)

        materials = self.get_materials(join=True)

        for mat in materials:
            mat.time_update(None, [Struct(terms=[self])])

    def call_get_fargs(self, args, kwargs):
        try:
            fargs = self.get_fargs(*args, **kwargs)

        except (RuntimeError, ValueError):
            terms.errclear()
            raise

        return fargs

    def call_function(self, out, fargs):
        try:
            status = self.function(out, *fargs)

        except (RuntimeError, ValueError):
            terms.errclear()
            raise

        if status:
            terms.errclear()
            raise ValueError('term evaluation failed! (%s)' % self.name)

        return status

    def eval_real(self, shape, fargs, mode='eval', term_mode=None,
                  diff_var=None, **kwargs):
        out = nm.empty(shape, dtype=nm.float64)

        if mode == 'eval':
            status = self.call_function(out, fargs)
            # Sum over elements but not over components.
            out1 = nm.sum(out, 0).squeeze()
            return out1, status

        else:
            status = self.call_function(out, fargs)

            return out, status

    def eval_complex(self, shape, fargs, mode='eval', term_mode=None,
                     diff_var=None, **kwargs):
        rout = nm.empty(shape, dtype=nm.float64)

        fargsd = split_complex_args(fargs)

        # Assuming linear forms. Then the matrix is the
        # same both for real and imaginary part.
        rstatus = self.call_function(rout, fargsd['r'])
        if (diff_var is None) and len(fargsd) >= 2:
            iout = nm.empty(shape, dtype=nm.float64)
            istatus = self.call_function(iout, fargsd['i'])

            if mode == 'eval' and len(fargsd) >= 4:
                irout = nm.empty(shape, dtype=nm.float64)
                irstatus = self.call_function(irout, fargsd['ir'])
                riout = nm.empty(shape, dtype=nm.float64)
                ristatus = self.call_function(riout, fargsd['ri'])

                out = (rout - iout) + (riout + irout) * 1j
                status = rstatus or istatus or ristatus or irstatus

            else:
                out = rout + 1j * iout
                status = rstatus or istatus

        else:
            out, status = rout + 0j, rstatus

        if mode == 'eval':
            out1 = nm.sum(out, 0).squeeze()
            return out1, status

        else:
            return out, status

    def evaluate(self, mode='eval', diff_var=None,
                 standalone=True, ret_status=False, **kwargs):
        """
        Evaluate the term.

        Parameters
        ----------
        mode : 'eval' (default), or 'weak'
            The term evaluation mode.

        Returns
        -------
        val : float or array
            In 'eval' mode, the term returns a single value (the
            integral, it does not need to be a scalar), while in 'weak'
            mode it returns an array for each element.
        status : int, optional
            The flag indicating evaluation success (0) or failure
            (nonzero). Only provided if `ret_status` is True.
        iels : array of ints, optional
            The local elements indices in 'weak' mode. Only provided in
            non-'eval' modes.
        """
        if standalone:
            self.standalone_setup()

        kwargs = kwargs.copy()
        term_mode = kwargs.pop('term_mode', None)

        if mode in ('eval', 'el_eval', 'el_avg', 'qp'):
            args = self.get_args(**kwargs)
            self.check_shapes(*args)

            emode = 'eval' if mode == 'el_eval' else mode
            _args = tuple(args) + (emode, term_mode, diff_var)
            fargs = self.call_get_fargs(_args, kwargs)

            shape, dtype = self.get_eval_shape(*_args, **kwargs)

            if dtype == nm.float64:
                val, status = self.eval_real(shape, fargs, mode, term_mode,
                                             **kwargs)

            elif dtype == nm.complex128:
                val, status = self.eval_complex(shape, fargs, mode, term_mode,
                                                **kwargs)

            else:
                raise ValueError('unsupported term dtype! (%s)' % dtype)

            val *= self.sign
            out = (val,)

        elif mode == 'weak':
            varr = self.get_virtual_variable()
            if varr is None:
                raise ValueError('no virtual variable in weak mode! (in "%s")'
                                 % self.get_str())

            if diff_var is not None:
                varc = self.get_variables(as_list=False)[diff_var]

            args = self.get_args(**kwargs)
            self.check_shapes(*args)

            _args = tuple(args) + (mode, term_mode, diff_var)
            fargs = self.call_get_fargs(_args, kwargs)

            n_elr, n_qpr, dim, n_enr, n_cr = self.get_data_shape(varr)
            n_row = n_cr * n_enr

            if diff_var is None:
                shape = (n_elr, 1, n_row, 1)

            else:
                n_elc, n_qpc, dim, n_enc, n_cc = self.get_data_shape(varc)
                n_col = n_cc * n_enc

                shape = (n_elr, 1, n_row, n_col)

            if varr.dtype == nm.float64:
                vals, status = self.eval_real(shape, fargs, mode, term_mode,
                                              diff_var, **kwargs)

            elif varr.dtype == nm.complex128:
                vals, status = self.eval_complex(shape, fargs, mode, term_mode,
                                                 diff_var, **kwargs)

            else:
                raise ValueError('unsupported term dtype! (%s)'
                                 % varr.dtype)

            if not isinstance(vals, tuple):
                vals *= self.sign
                iels = self.get_assembling_cells(vals.shape)

            else:
                vals = (self.sign * vals[0],) + vals[1:]
                iels = None

            out = (vals, iels)

        if goptions['check_term_finiteness']:
            assert_(nm.isfinite(out[0]).all(),
                    msg='"%s" term values not finite!' % self.get_str())

        if ret_status:
            out = out + (status,)

        if len(out) == 1:
            out = out[0]

        return out

    def assemble_to(self, asm_obj, val, iels, mode='vector', diff_var=None):
        """
        Assemble the results of term evaluation.

        For standard terms, assemble the values in `val` corresponding to
        elements/cells `iels` into a vector or a CSR sparse matrix `asm_obj`,
        depending on `mode`.

        For terms with a dynamic connectivity (e.g. contact terms), in
        `'matrix'` mode, return the extra COO sparse matrix instead. The extra
        matrix has to be added to the global matrix by the caller. By default,
        this is done in :func:`Equations.evaluate()
        <sfepy.discrete.equations.Equations.evaluate()>`.
        """
        import sfepy.discrete.common.extmods.assemble as asm

        vvar = self.get_virtual_variable()
        dc_type = self.get_dof_conn_type()

        extra = None

        if mode == 'vector':
            if asm_obj.dtype == nm.float64:
                assemble = asm.assemble_vector

            else:
                assert_(asm_obj.dtype == nm.complex128)
                assemble = asm.assemble_vector_complex
                for ii in range(len(val)):
                    if not(val[ii].dtype == nm.complex128):
                        val[ii] = nm.complex128(val[ii])

            if not isinstance(val, tuple):
                dc = vvar.get_dof_conn(dc_type)
                assert_(val.shape[2] == dc.shape[1])

                assemble(asm_obj, val, iels, 1.0, dc)

            else:
                vals, rows, var = val
                if var.eq_map is not None:
                    eq = var.eq_map.eq

                    rows = eq[rows]
                    active = (rows >= 0)
                    vals, rows = vals[active], rows[active]

                # Assumes no repeated indices in rows!
                asm_obj[rows] += vals

        elif mode == 'matrix':
            if asm_obj.dtype == nm.float64:
                assemble = asm.assemble_matrix

            else:
                assert_(asm_obj.dtype == nm.complex128)
                assemble = asm.assemble_matrix_complex

            svar = diff_var
            tmd = (asm_obj.data, asm_obj.indptr, asm_obj.indices)

            if ((asm_obj.dtype == nm.complex128)
                and (val.dtype == nm.float64)):
                val = val.astype(nm.complex128)

            sign = 1.0
            if self.arg_derivatives[svar.name]:
                if not self.is_quasistatic or (self.step > 0):
                    sign *= 1.0 / self.dt

                else:
                    sign = 0.0

            if not isinstance(val, tuple):
                rdc = vvar.get_dof_conn(dc_type)

                is_trace = self.arg_traces[svar.name]
                cdc = svar.get_dof_conn(dc_type, is_trace=is_trace)
                assert_(val.shape[2:] == (rdc.shape[1], cdc.shape[1]))

                assemble(tmd[0], tmd[1], tmd[2], val, iels, sign, rdc, cdc)

            else:
                from scipy.sparse import coo_matrix

                vals, rows, cols, rvar, cvar = val
                if rvar.eq_map is not None:
                    req, ceq = rvar.eq_map.eq, cvar.eq_map.eq

                    rows, cols = req[rows], ceq[cols]
                    active = (rows >= 0) & (cols >= 0)
                    vals, rows, cols = vals[active], rows[active], cols[active]

                extra = coo_matrix((sign * vals, (rows, cols)),
                                   shape=asm_obj.shape)

        else:
            raise ValueError('unknown assembling mode! (%s)' % mode)

        return extra

max_stars_repo_name: lokik/sfepy

max_issues_repo_name: lokik/sfepy

max_forks_repo_name: lokik/sfepy


text:
#! /usr/bin/python
# -*- coding: utf-8 -*-
import sys
import os
path_to_script = os.path.dirname(os.path.abspath(__file__))
sys.path.append(os.path.join(path_to_script, "../extern/pyseg_base/src/"))
sys.path.append(os.path.join(path_to_script, "../extern/sed3/"))
sys.path.append(os.path.join(path_to_script, "../src/"))
import numpy as np

from lisa import organ_segmentation
import imcut.dcmreaddata as dcmr
import os.path

def main():

    dcmdir = os.path.join(path_to_script,'./../sample_data/matlab/examples/sample_data/DICOM/digest_article/')
    #data3d, metadata = dcmr.dcm_read_from_dir(dcmdir)
    oseg = organ_segmentation.OrganSegmentation(dcmdir, working_voxelsize_mm = 4, manualroi=False)

    oseg.interactivity()

main()

max_stars_repo_name: mjirik/lisa

max_issues_repo_name: mjirik/lisa

max_forks_repo_name: mjirik/lisa