text:
#!/usr/bin/python

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
import gym
import threading,queue

EP_MAX = 500
EP_LEN = 200
GAMMA = 0.9
A_LR = 0.0001
C_LR = 0.0002
MIN_BATCH_SIZE = 32
UPDATE_STEPS = 10
EPSION = 0.2
GAME = 'Pendulum-v0'
S_DIM,A_DIM = 3,1
N_WORKER = 2

class PPO(object):

    def __init__(self,coord,data_queue,rolling_event,update_event):
        self.sess = tf.Session()
        self.tfs = tf.placeholder(tf.float32,[None,S_DIM],'state')
        self.coord = coord
        self.data_queue = data_queue
        self.rolling_event = rolling_event
        self.update_event = update_event

        #critic
        with tf.variable_scope('critic'):
            l1 = tf.layers.dense(self.tfs,100,tf.nn.relu)
            self.v = tf.layers.dense(l1,1)
            self.tfdc_r = tf.placeholder(tf.float32,[None,1],'discounted_r')
            self.advantage = self.tfdc_r - self.v
            self.closs = tf.reduce_mean(tf.square(self.advantage))
            self.ctrain_op = tf.train.AdamOptimizer(C_LR).minimize(self.closs)

        #Actor
        pi,pi_params = self._build_anet('pi',trainable=True)
        oldpi,oldpi_params = self._build_anet('oldpi',trainable=False)

        with tf.variable_scope('sample_action'):
            self.sample_op = tf.squeeze(pi.sample(1),axis=0)
        with tf.variable_scope('update_oldpi'):
            self.update_oldpi_op = [oldp.assign(p) for p,oldp in zip(pi_params,oldpi_params)]

        self.tfa = tf.placeholder(tf.float32,[None,A_DIM],'action')
        self.tfadv = tf.placeholder(tf.float32,[None,1],'advantage')
        with tf.variable_scope('loss'):
            with tf.variable_scope('surrogate'):
                ratio = pi.prob(self.tfa) / oldpi.prob(self.tfa)
                surr = ratio * self.tfadv

                #clipped surrogate objection
                self.aloss = -tf.reduce_mean(tf.minimum(surr,tf.clip_by_value(ratio,1. - EPSION,1 + EPSION)*self.tfadv))

        with tf.variable_scope('atrain'):
            self.atrain_op = tf.train.AdamOptimizer(A_LR).minimize(self.aloss)
        self.sess.run(tf.global_variables_initializer())

    def update(self):
        global GLOBAL_UPDATE_COUNTER

        while not self.coord.should_stop():
            if GLOBAL_EP < EP_MAX:
                self.update_event.wait()
                self.sess.run(self.update_oldpi_op)

                data = [self.data_queue.get() for _ in range(self.data_queue.qsize())]

                s = np.vstack([item[0] for item in data])
                a = np.vstack([item[1] for item in data])
                r = np.vstack([item[2] for item in data])
                #data = np.vstack(data)
                #s,a,r = data[:,:S_DIM],data[:, S_DIM:S_DIM + A_DIM], data[:, -1:]
                adv = self.sess.run(self.advantage,{self.tfs:s,self.tfdc_r:r})
                #update actor
                [self.sess.run(self.atrain_op,{self.tfs:s ,self.tfa: a,self.tfadv: adv}) for _ in range(UPDATE_STEPS)]
                #update critic
                [self.sess.run(self.ctrain_op,{self.tfs: s, self.tfdc_r: r}) for _ in range(UPDATE_STEPS)]

                self.update_event.clear()
                GLOBAL_UPDATE_COUNTER = 0
                self.rolling_event.set()

    def _build_anet(self, name, trainable):
        with tf.variable_scope(name):
            l1 = tf.layers.dense(self.tfs, 100, tf.nn.relu,trainable=trainable)
            mu = 2 * tf.layers.dense(l1, A_DIM, tf.nn.tanh, trainable=trainable)
            sigma = tf.layers.dense(l1,A_DIM,tf.nn.softplus,trainable=trainable)

            norm_dist = tf.distributions.Normal(loc=mu,scale=sigma)
        params = tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES,scope=name)
        return norm_dist, params

    def choose_action(self,s):
        s = s[np.newaxis,:]
        a = self.sess.run(self.sample_op,{self.tfs:s})[0]
        return np.clip(a,-2,2)

    def get_v(self,s):
        if s.ndim < 2:
            s = s[np.newaxis,:]
        return self.sess.run(self.v,{self.tfs:s})[0,0]

class Worker(object):
    def __init__(self,wid,global_ppo,coord,rolling_event,update_event,data_queue,max_ep = EP_MAX,ep_len = EP_LEN,mini_batch_size = MIN_BATCH_SIZE,gamma=GAMMA):
        self.wid = wid
        self.env = gym.make(GAME).unwrapped
        self.ppo = global_ppo
        self.coord = coord
        self.rolling_event = rolling_event
        self.update_event = update_event
        self.data_queue = data_queue
        self.max_ep = max_ep
        self.ep_len = ep_len
        self.mini_batch_size = mini_batch_size
        self.gamma = gamma

    def work(self):
        global GLOBAL_EP, GLOBAL_UPDATE_COUNTER
        while not self.coord.should_stop():
            s = self.env.reset()
            ep_r = 0
            buffer_s,buffer_a,buffer_r = [],[],[]

            for t in range(self.ep_len):
                if not self.rolling_event.is_set():
                    self.rolling_event.wait()
                    buffer_s,buffer_a,buffer_r = [],[],[]
                a = self.ppo.choose_action(s)
                s_, r, done, _ = self.env.step(a)

                buffer_s.append(s)
                buffer_a.append(a)
                buffer_r.append((r+8)/8)
                s = s_
                ep_r += r

                GLOBAL_UPDATE_COUNTER += 1
                if t == self.ep_len -1 or GLOBAL_UPDATE_COUNTER >= self.mini_batch_size:
                    v_s_ = self.ppo.get_v(s_)
                    discounted_r = []
                    for r in buffer_r[::-1]:
                        v_s_ = r + self.gamma * v_s_
                        discounted_r.append(v_s_)
                    discounted_r.reverse()

                    #bs, ba, br = np.vstack(buffer_s),np.vstack(buffer_a),np.array(discounted_r)[:,np.newaxis]
                    print("Worker[%i] collects %d samples!"%(self.wid,len(buffer_s)))
                    for index in range(len(buffer_s)):
                        self.data_queue.put((buffer_s[index],buffer_a[index],discounted_r[index]))
                    buffer_s,buffer_a,buffer_r = [],[],[]

                    if GLOBAL_UPDATE_COUNTER >= self.mini_batch_size:
                        self.update_event.set()
                        self.rolling_event.clear()

                    if GLOBAL_EP >= self.max_ep:
                        self.coord.request_stop()
                        break
            GLOBAL_EP += 1
            print("%d-%d  Worker[%i] reward: %.2f"%(GLOBAL_EP,self.max_ep,self.wid,ep_r))

GLOBAL_UPDATE_COUNTER,GLOBAL_EP = 0 , 0
def main():
    env = gym.make('Pendulum-v0').unwrapped
    COORD = tf.train.Coordinator()
    UPDATE_EVENT, ROLLING_EVENT = threading.Event(),threading.Event()
    QUEUE = queue.Queue()
    UPDATE_EVENT.clear()
    ROLLING_EVENT.set()
    GLOBAL_PPO = PPO(COORD,data_queue=QUEUE,rolling_event = ROLLING_EVENT,update_event = UPDATE_EVENT)

    workers = [Worker(wid=i,global_ppo = GLOBAL_PPO,coord=COORD,rolling_event = ROLLING_EVENT,update_event = UPDATE_EVENT, data_queue = QUEUE) for i in range(N_WORKER)]

    threads = []
    for worker in workers:
        t = threading.Thread(target=worker.work,args=())
        t.start()
        threads.append(t)

    threads.append(threading.Thread(target=GLOBAL_PPO.update,args=()))
    threads[-1].start()
    COORD.join(threads)

    env = gym.make('Pendulum-v0').unwrapped
    while True:
        s = env.reset()
        for t in range(300):
            env.render()
            s = env.step(GLOBAL_PPO.choose_action(s))[0]

if __name__ == "__main__":
    main()

max_stars_repo_name: xiaoyangzai/DeepReinforcementLearning

max_issues_repo_name: xiaoyangzai/DeepReinforcementLearning

max_forks_repo_name: xiaoyangzai/DeepReinforcementLearning


text:
### Vectorized equivalents of freqs, frequtils, etc.

from music21 import *
import numpy as np

tones = [n + a for n in 'CDEFGABcdefgab' for a in ['-', '#', '']]
keys = map(key.Key, tones)
chords = [chord.Chord(map(k.pitches.__getitem__, [0, 2, 4])) for k in keys] # TODO does not include diminished
notes = [n for n in tones if n[0].isupper()]
K, C, N = len(keys), len(chords), len(notes)

C = F = 7*7*4
N = T = 7*7

# vectorize these
# may also use music21.roman later

def Func(chord, key=key.Key('C')):
    try:
        num = int(chord.root().diatonicNoteNum - key.tonic.diatonicNoteNum) % 7
        acc = int(chord.root().ps - key.pitches[num].ps + 3) % 12 # assuming only 5 different types of accidentals occur
        qual = ['major','minor','diminished','augmented', 'other'].index(chord.quality) % 4
        return np.ravel_multi_index((num, acc, qual), (7, 7, 4))
    except AttributeError:
        return -1

def Chord(func, key=key.Key('C')):
    num, acc, qual = np.unravel_index(func, (7,7,4))
    root = Pitch(np.ravel_multi_index((num, acc), (7, 7)), key)
    # TODO

def Tone(note, key=key.Key('C')): # if no key, use absolute value
    try:
        num = int(note.diatonicNoteNum - key.tonic.diatonicNoteNum) % 7
        acc = int(note.ps - key.pitches[num].ps + 3) % 12
        return np.ravel_multi_index((num, acc), (7, 7)) # store accidental
    except:
        return -1

def Pitch(tone, key=key.Key('C')):
    try:
        num, acc = np.unravel_index(tone, (7,7))
        p = key.pitches[num]
        return pitch.Pitch(p.name[0], accidental=p.alter+acc-3)
    except:
        return None

abs2func_tone = np.array([[Tone(Pitch(t, k)) for t in xrange(T)] for k in keys])
func2abs_tone = np.array([[Tone(Pitch(t), k) for t in xrange(T)] for k in keys])
def transpose_tone(hist, key, tofunc=True):
    'Transpose histogram of absolute tones to tones in a given key.'
    if tofunc:
        arr = abs2func_tone[keys.index(key)]
    else:
        arr = func2abs_tone[keys.index(key)]
    return hist[arr]*(arr!=-1)

# matrices converting absolute notes/chords to functions in a key
C = F
# chord2func = np.array([[np.arange(F)==Func(Chord(f, k)) for f in xrange(F)] for k in keys])
# func2chord = np.array([[np.arange(F)==Func(Chord(f), k) for f in xrange(F)] for k in keys])
# pitch2tone = np.array([[np.arange(T)==Tone(Pitch(t, k)) for t in xrange(T)] for k in keys])
# tone2pitch = np.array([[np.arange(T)==Tone(Pitch(t), k) for t in xrange(T)] for k in keys])

def Key(key):
    return keys.index(key)

class Sample:
    '''
    A processed sample chorale containing a list of chords, with optionally specified harmonic velocity.
    '''
    def __init__(self, filename, fromCorpus=False):
        self.filename = filename
        if fromCorpus:
            self.s = corpus.parse(filename)
        else:
            self.s = converter.parse(filename)
        # self.cs = self.s.flat.notesAndRests.stream().chordify()

        self.key = self.s.analyze('krumhansl')
        self.vel = None # TODO change this, perhaps use qualities other than vel, such as measure ends
        self.matrix = None

    # function converting sample into TxN matrix of notes + Tx1 array of absolute durations
    def get_matrix(self):
        if self.matrix is None:
            s = self.s.flat.notesAndRests.stream()
            endtimes = s._uniqueOffsetsAndEndTimes(endTimesOnly=True)
            l = len(endtimes)
            self.matrix = np.zeros((l,T))
            self.ts = np.zeros(l)
            t0 = 0
            for i, t in enumerate(endtimes):
                a = s.getElementsByOffset(t0, t, includeEndBoundary=False)
                for n in a:
                    if isinstance(n, chord.Chord):
                        for p in n.pitches:
                            self.matrix[i, Tone(p)] += t - t0
                    elif isinstance(n, note.Note):
                        self.matrix[i, Tone(n.pitch)] += t - t0
                self.ts[i] = (t0 + t)/2
                t0 = t
        return self.matrix, self.ts

    def measures(self):
        'Iterates through each measure of sample.'
        for m in self.cs:
            if isinstance(m, stream.Measure):
                yield m

    def chords(self):
        'Iterates through each chord of sample.'
        for m in self.measures():
            for c in m.notes: # cannot call notes directly apparently
                yield c

    def notes(self, voice):
        'Iterates through each note in sample at a given voice.'
        for p in self.s.parts[voice]:
            for m in p.getElementsByClass('Measure'):
                for n in m.notes:
                    yield n

    # iterable overrides, may be made more efficient
    def __iter__(self):
        return self.chords()
    def __len__(self):
        return len(self.cs)
    def __getitem__(self, key):
        return list(self.chords())[key] # can also take slice objects

    def __hash__(self):
        return hash(self.filename)
    def __eq__(self, other):
        return isinstance(other, Sample) and self.filename == other.filename
    def __ne__(self, other):
        return not self.__eq__(other)

    pass

class Stats:
    '''
    Stores counts of chords, notes, transitions, etc.
    '''
    def __init__(self, sample=None):
        self.cs = np.zeros((2,C)) # np.zeros((C,))
        self.ts = np.zeros((2,F,F)) # np.zeros((C,C)) # no of transitions between each two functions, for each key mode
        self.ns = np.zeros((T,))
        self.es = np.zeros((2,T,F)) # instances of each tone in a chord, for each key mode
        self.vs = np.zeros((T,T))
        self.ks = np.zeros((K,))
        # self.ts = np.zeros((2,F,F)) # no of transitions between each two functions, for each key mode
        # self.ns = np.zeros((2,F,T)) # instances of each tone in a chord, for each key mode
        # self.cs may store initial chords for each key

        if sample is not None:
            # construct stats from sample matrix
            self.ks[Key(sample.key)] += 1
            # construct note histogram
            m = sample.get_matrix()[0]
            self.ns += transpose_tone(m.sum(axis=0), sample.key) # remove later
            # TODO handle the rest of the arrays
            c1 = None
            for t in xrange(m.shape[0]):
                # infer which chord is active
                notes = [Pitch(n) for n in m[t].nonzero()[0]]
                if len(notes) == 0:
                    continue
                c = chord.Chord(notes)
                # if t % 50 == 0:
                #     print c, m[t].sum(), Func(c, sample.key), Chord(Func(c, sample.key), sample.key)
                if c1 is None:
                    self.cs[int(sample.key.mode == 'major'), Func(c,sample.key)] += m[t].sum()
                else: # may distribute 'other' chords to every other mode
                    self.ts[int(sample.key.mode == 'major'), Func(c1,sample.key), Func(c,sample.key)] += m[t].sum()
                for n in notes:
                    self.es[int(sample.key.mode == 'major'), Tone(n,sample.key), Func(c,sample.key)] += m[t,Tone(n)]
                c1 = c

    def __add__(self, other):
        res = Stats()
        res.cs = self.cs + other.cs
        res.ts = self.ts + other.ts
        res.ns = self.ns + other.ns
        res.vs = self.vs + other.vs
        res.es = self.es + other.es
        res.ks = self.ks + other.ks
        return res

    def __radd__(self, other):
        return other.__add__(self)

    # add methods returning NxK etc. probability matrices looking up entries for Tone(n,k) etc.

class Env:
    '''
    Program environment containing probabilities for each chord, note and transition, acquired from sample chord sequences.
    '''
    def __init__(self):
        self.stats = Stats()
        self.samples = dict()

    def train(self, filenames, verbose=False):
        'Train probabilities from given iterator of filenames, for example corpus.getBachChorales().'
        for f in filenames:
            s = Sample(f)
            if s not in self.samples:
                if verbose:
                    print 'Processing %s...' % f
                self.process(s)

    def process(self, sample):
        'Update probabilities based on a sample sequence of chords.'
        if sample in self.samples:
            return

        # construct Stats from Sample
        # may also look at small windows in sample
        self.samples[sample] = Stats(sample)
        self.stats += self.samples[sample]

    # TODO fill these below, or move some to Stats

    def kprob(self, k):
        'kprob(k) -> Return probability of key k.'
        return (self.stats.ks[k]+1) / (self.stats.ks+1).sum()

    def cprob(self, c, k):
        'cprob(c, k) -> Return probability of chord c occurring in key k.'
        return self.cfreq[True, Func(c, k)]

    def nprob(self, n, c, k, v = None):
        'nprob(n, c, k[, v]) -> Return probability of note n occurring in voice v of chord c in key k.'
        return self.nfreq[True, Tone(n, k), k.mode, Func(c, k), v] # should normalize over v if not none

    def tprob(self, c1, c, k, vel = None):
        'tprob(c1, c, k[, vel]) -> Return probability of chord c changing to chord c1 in key k under harmonic velocity vel.'
        if c == c1:
            return self.cprob(c, k) # perhaps dependent on vel
        f, f1 = Func(c, k), Func(c1, k)
        return self.tfreq[True, f1, f, s]

    def vprob(self, n1, n, c1, c, k, v = None, vel = None):
        '''vprob(n1, n, c1, c, k[, v, vel]) ->
        Return probability of note n in chord c changing to note n1 in chord c1
        in key k and voice v under harmonic velocity vel.'''
        if n == n1 and c == c1:
            return self.nprob(n, c, k, v) # perhaps dependent on vel
        f, f1 = Func(c, k), Func(c1, k)
        t, t1 = Tone(n, k), Tone(n1, k)
        return self.vfreq[True, t1, t, f1, f, v, s]

table = Env()

ext: py

max_stars_repo_name: meagtan/harmonizer

max_issues_repo_name: meagtan/harmonizer

max_forks_repo_name: meagtan/harmonizer

reason: import numpy


text:
import matplotlib.pyplot as plt
import numpy as np
import csv

def is_number(s):
    try:
        float(s)
        return True
    except ValueError:
        pass

    try:
        import unicodedata
        unicodedata.numeric(s)
        return True
    except (TypeError, ValueError):
        pass

    return False

def process_columns(age, weight):
    pairs = list(zip(age, weight))
    filtered = [(int(t[0]), float(t[1])) for t in pairs
                if is_number(t[0]) and is_number(t[1])]
    flat = list(zip(*filtered))
    return flat

with open('baby_data.csv', newline='') as csvfile:
    reader = csv.reader(csvfile, delimiter=',')
    columns = list(zip(*reader))
    perrin_data = process_columns(columns[1], columns[4])
    tobin_data = process_columns(columns[7], columns[10])

# with plt.xkcd():
if True:
    # Based on "Stove Ownership" from XKCD by Randall Munroe
    # https://xkcd.com/418/

    fig = plt.figure()
    ax = fig.add_axes((0.1, 0.2, 0.8, 0.7))
    # ax.spines['right'].set_color('none')
    # ax.spines['top'].set_color('none')
    ax.set_xticks(np.arange(0, 100, step=7))
    # ax.set_yticks([])
    ax.set_xlim([0, 30])
    ax.set_ylim([4, 14])

    # ax.annotate(
    #     'THE DAY I REALIZED\nI COULD COOK BACON\nWHENEVER I WANTED',
    #     xy=(70, 1), arrowprops=dict(arrowstyle='->'), xytext=(15, -10))

    ax.scatter(perrin_data[0], perrin_data[1])
    ax.scatter(tobin_data[0], tobin_data[1])
    # ax.plot(tobin_weight)
    ax.legend(['Perrin', 'Tobin'])
    ax.grid(True)

    ax.set_xlabel('Age (days)')
    ax.set_ylabel('Weight (lbs)')

    # fig.text(
    #     0.5, 0.05,
    #     '"Stove Ownership" from xkcd by Randall Munroe',
    #     ha='center')

    plt.show()


text:
from numpy.testing import assert_allclose

from blocks.extensions import Timing, FinishAfter
from blocks.utils.testing import MockMainLoop

def test_timing():
    epochs = 2
    main_loop = MockMainLoop(delay_time=0.1,
                             extensions=[Timing(prefix='each'),
                                         Timing(prefix='each_second',
                                                every_n_epochs=2),
                                         FinishAfter(after_n_epochs=epochs)])
    main_loop.run()
    iterations = int(main_loop.log.status['iterations_done'] / epochs)
    assert_allclose(
        (main_loop.log[iterations]['each_time_train_this_epoch'] +
         main_loop.log[iterations]['each_time_train_this_epoch']) / 2,
        main_loop.log.current_row['each_second_time_train_this_epoch'],
        atol=1e-2)
    assert 'each_time_read_data_this_epoch' in main_loop.log[iterations]
    assert 'each_second_time_read_data_this_epoch' in main_loop.log[iterations]

max_stars_repo_name: KIKOcaoyue/blocks

max_issues_repo_name: KIKOcaoyue/blocks

max_forks_repo_name: KIKOcaoyue/blocks

reason: from numpy


text:
"""
test_data.py

Author: Matthew Yu, Array Lead (2020).
Contact: matthewjkyu@gmail.com
Created: 8/16/20
Last Modified: 8/16/20
Description: Parses, cleans, and displays cell data collected from the 9-14-19 test.
    Rough and dirty script to grab some data to use for simulator testing. Use at your own risk.
"""
import sys
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

def round_num(x, prec=2, base=.05):
  return round(base * round(float(x)/base),prec)

f = open("Array Testing 9-14-19/Array2_Test2_9_14_19_fixed.dat", "r")
# remove first line with column headers
f.readline()

cycle = []
temp1 = []
temp2 = []

while True:
    # only read every five lines
    line = f.readline()

    if not line:
        break
    try:
        cycle.append(float(line.strip('\n').split(' ')[0]))
        temp1.append(float(line.strip('\n').split(' ')[1]))
        temp2.append(float(line.strip('\n').split(' ')[2]))
    except Exception as e:
        print(e)

    # number of lines to skip
    for i in range(30):
        line = f.readline()

fig, ax = plt.subplots(2, 1)
fig.set_size_inches(9, 10)

ax[0].set_xlabel('Cycle')
ax[0].set_ylabel('Temp 1')
ax[0].set_xlim([0, 5000])
ax[0].set_ylim([0, 150])
ax[0].scatter(cycle, temp1, color='r', marker='o', s=[1])
ax[0].title.set_text("Array temp 1 as a function of cycle")

ax[0].grid()

ax[1].set_xlabel('Cycle')
ax[1].set_ylabel('Temp 2')
ax[1].set_xlim([0, 5000])
ax[1].set_ylim([0, 150])
ax[1].scatter(cycle, temp2, color='r', marker='o', s=[1])
ax[1].title.set_text("Array temp 2 as a function of cycle")
ax[1].grid()

plt.tight_layout()

fig.savefig("test.png")
plt.show()

f = open("reformatted_data.txt", "w")
for idx in range(0, len(cycle)):
    line = "[" + str(idx) + ", " + str(1000) + ", " + str(round_num(temp1[idx], 1, .5)) + "],\n"
    f.write(line)

max_stars_repo_name: lhr-solar/Array-Sim-Training

max_issues_repo_name: lhr-solar/Array-Sim-Training

max_forks_repo_name: lhr-solar/Array-Sim-Training


text:
# Area under the curve splitting criterion

import numpy as np
from sklearn.linear_model import LinearRegression  #LogisticRegression
from sklearn.metrics import accuracy_score, roc_auc_score
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
import sys
sys.path.insert(0, '../dt')

from node_types import Predicate

class Split_auc:
  def __init__(self, use_svm=False):
    self.b_score = -1.
    self.b_pos = -1
    self.b_val = -1
    self.b_eq = None
    self.feature_mask = []
    self.use_svm = use_svm
    self.EPSILON = 0.00000001
    self.pos_uses = []

  def best(self, data, node):
    self.b_score = 0.
    self.b_pos = -1
    self.b_val = -1
    self.b_eq = None

    # Create a mask to only include features that
    # do not contain the same value in all samples
    self.feature_mask = []
    self.pos_uses = []
    for i, dom in enumerate(data.Xdomains):
      assert(len(dom) >= 1)
      self.feature_mask.append(len(dom) > 1)
      self.pos_uses.append(0)

    cur = node.parent
    while (cur is not None):
      assert(cur.is_predicate())
      assert(cur.predicate.fpos >= 0 and cur.predicate.fpos < len(self.pos_uses))
      self.pos_uses[cur.predicate.fpos] += 1
      cur = cur.parent

    # Compute for each predicate
    for i, dom in enumerate(data.Xdomains):
      assert(len(dom) > 0)
      if (len(dom) == 2):
        # {0,1} --> =1
        # {0,8} --> =8
        self.split_score(data, i, max(dom), True, dom)
      elif (len(dom) > 2):
        # {0,1,2} --> =0 =1 =2  (<1 IS =0; <2 IS =!2)
        # {0,3,4,6} --> =0 =3 =4 <4 (<3 IS =0; <6 IS =!6)
        for idx, val in enumerate(sorted(dom)):
          assert(idx == 0 or val != min(dom))
          also_ineq = (i not in data.Xineqforbidden and
                       idx >= 2 and val != max(dom))
          if (also_ineq):
            self.split_score(data, i, val, False, dom)
          self.split_score(data, i, val, True, dom)

    # Done; return the best predicate
    assert(self.b_eq is not None)
    assert(self.b_pos >= 0)
    assert(self.b_pos < data.Xnames.size)

    numname = None
    if (data.Xnames[self.b_pos] == 'module'):
      assert(self.b_eq)
      assert(self.b_val in data.ModuleIDtoName)
      numname = data.ModuleIDtoName[self.b_val]
    elif (data.Xnames[self.b_pos] == 'action'):
      assert(self.b_eq)
      assert(self.b_val in data.ActionIDtoName)
      numname = data.ActionIDtoName[self.b_val]

    return Predicate(fname=data.Xnames[self.b_pos], fpos=self.b_pos,
                     equality=self.b_eq, number=self.b_val, numberName=numname)

  def split_score(self, data, pos, value, equality, domain):
    pred = Predicate(fname='', fpos=pos, equality=equality, number=value)
    mask = (pred.evaluate_matrix(data.X))
    sat_Y = data.Y[mask]
    if (sat_Y.size == 0 or sat_Y.size == data.Y.size):
      return
    sat_X = data.X[mask]
    sat_X = sat_X[:,self.feature_mask]  # Apply feature mask
    unsat_X = data.X[~mask]
    unsat_X = unsat_X[:,self.feature_mask]  # Apply feature mask
    unsat_Y = data.Y[~mask]

    reg = LinearRegression()
    #reg = LogisticRegression(dual=False, tol=0.0001, C=100.0, solver='liblinear',
                             #fit_intercept=True, random_state=31337)
    clf = None if not self.use_svm else \
          LinearSVC(penalty='l1', tol=0.000001, C=10000.0,
                    dual=False, fit_intercept=True, random_state=42)
    area = {'sat': 0., 'unsat': 0.}

    sides_done = set()
    sides_done_clean = 0
    for (X, Y, name) in [(sat_X, sat_Y, 'sat'), (unsat_X, unsat_Y, 'unsat')]:
      assert(X.shape[0] > 0)
      Ys_present = set()
      for y in Y:
        Ys_present.add(y)
      assert(len(Ys_present) > 0)
      if (len(Ys_present) == 1):  # only one class present
        area[name] = 1.
        sides_done.add(name)
        sides_done_clean += 1
      elif (len(Ys_present) == 2):
        X_tr = StandardScaler().fit_transform(X)
        if (self.use_svm):
          # SVM classifier
          clf.fit(X_tr, Y)
          ac = accuracy_score(normalize=False, y_true=Y,
                              y_pred=clf.predict(X_tr))
          if (ac == Y.size):
            area[name] = 1.
            sides_done.add(name)
          else:
            reg.fit(X_tr, Y)
            area[name] = roc_auc_score(y_true=Y, y_score=reg.predict(X_tr))
        else:
          # Linear regressor
          reg.fit(X_tr, Y)
          area[name] = roc_auc_score(y_true=Y, y_score=reg.predict(X_tr))
      assert(area[name] < 1. + self.EPSILON)

    assert(len(sides_done) <= 2)
    if (len(sides_done) == 2):
      # Making sure solving split wins
      assert(sides_done_clean <= 2)
      if (sides_done_clean == 2):
        # Two clean partitions - best
        area['sat'], area['unsat'] = 4., 4.
      elif (sides_done_clean == 1):
        # One clean one LC partition - second best
        area['sat'], area['unsat'] = 3., 3.
      else:
        # Two LC partitions - third best
        area['sat'], area['unsat'] = 2., 2.
    else:
      # Punishment of obviously unfavourable split
      if (pos not in data.Xineqforbidden):
        assert(pos >= 0 and pos < len(self.pos_uses))
        area['sat'] /= float(1 + self.pos_uses[pos])
        area['unsat'] /= float(1 + self.pos_uses[pos])
      #
      if (equality and (pos not in data.Xineqforbidden) and
          value != min(domain) and value != max(domain)):
        area['sat'] /= 2
        area['unsat'] /= 2
      #
      if (len(domain) > 2):
        sat_dom = 0
        for domv in domain:
          if (equality):
            if (domv == value):
              sat_dom += 1
          else:
            if (domv < value):
              sat_dom += 1
        unsat_dom = len(domain) - sat_dom
        assert(sat_dom >= 1 and unsat_dom >= 1)
        assert((not equality) or sat_dom == 1)
        for name, sz, dm in [('sat', sat_Y.size, sat_dom),
                             ('unsat', unsat_Y.size, unsat_dom)]:
          if (dm > 1 and (name not in sides_done) and
              sz / float(data.Y.size) > 0.99):
            area['sat'] /= float(dm)
            area['unsat'] /= float(dm)

    if (area['sat'] + area['unsat'] > self.b_score + self.EPSILON):
      self.b_score = area['sat'] + area['unsat']
      self.b_pos = pos
      self.b_val = value
      self.b_eq = equality

max_stars_repo_name: ViToSVK/dt-linear

max_issues_repo_name: ViToSVK/dt-linear

max_forks_repo_name: ViToSVK/dt-linear


text:
from ensemble.superknn import Superknn
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.ensemble import RandomForestClassifier
from xgboost import XGBClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.datasets import load_digits
import time
import warnings
warnings.filterwarnings(module='sklearn*', action='ignore', category=DeprecationWarning)

# Prepare data
digits = load_digits()
X, y = digits.data, digits.target
X_train, X_test, y_train, y_test = train_test_split(X, y, shuffle=True, test_size=0.2, random_state=0)
'''实验1：不对model调参'''

# print ('******************** 1.XGB *********************' )
# model = XGBClassifier(n_jobs=1, random_state=0)
# start_time = time.time()
# model.fit(X_train, y_train)
# print ('training(XGB) took %fs!' % (time.time() - start_time))
# start_time = time.time()
# #直接pred
# y_pred = model.predict(X_test)
# print ('predict(XGB-直接) took %fs!' % (time.time() - start_time))
# # #一个个pred
# # start_time = time.time()
# # y_pred2 = []
# # for vector in X_test:
# #     pred_temp = model.predict(vector.reshape(1, -1))
# #     y_pred2.append(pred_temp)
# # print ('predict(XGB-一个个) took %fs!' % (time.time() - start_time))
# print('Final prediction score with optimal(use XGB): [%.8f]' % accuracy_score(y_test, y_pred))
# #
# print ('******************** 2.Extra *********************' )
# model = ExtraTreesClassifier(n_jobs=1, random_state=0)
# start_time = time.time()
# model.fit(X_train, y_train)
# print ('training(EXtra) took %fs!' % (time.time() - start_time))
# start_time = time.time()
# #直接pred
# y_pred = model.predict(X_test)
# print ('predict(Extra-直接) took %fs!' % (time.time() - start_time))
# # #一个个pred
# # start_time = time.time()
# # y_pred2 = []
# # for vector in X_test:
# #     pred_temp = model.predict(vector.reshape(1, -1))
# #     y_pred2.append(pred_temp)
# # print ('predict(Extra-一个个) took %fs!' % (time.time() - start_time))
# print('Final prediction score with optimal(use Extra): [%.8f]' % accuracy_score(y_test, y_pred))
# #
# #
# print ('******************** 3.GBDT *********************' )
# model = GradientBoostingClassifier(random_state=0)
# start_time = time.time()
# model.fit(X_train, y_train)
# print ('training(GBDT) took %fs!' % (time.time() - start_time))
# start_time = time.time()
# #直接pred
# y_pred = model.predict(X_test)
# print ('predict(GBDT-直接) took %fs!' % (time.time() - start_time))
# #一个个pred
# start_time = time.time()
# y_pred2 = []
# for vector in X_test:
#     pred_temp = model.predict(vector.reshape(1, -1))
#     y_pred2.append(pred_temp)
# print ('predict(GBDT-一个个) took %fs!' % (time.time() - start_time))
# print('Final prediction score with optimal(use GBDT): [%.8f]' % accuracy_score(y_test, y_pred))
#
# print ('******************** 4.RF *********************' )
# model = RandomForestClassifier(n_jobs=1, random_state=0)
# start_time = time.time()
# model.fit(X_train, y_train)
# print ('training(RF) took %fs!' % (time.time() - start_time))
# start_time = time.time()
# #直接pred
# y_pred = model.predict(X_test)
# print ('predict(RF-直接) took %fs!' % (time.time() - start_time))
# # #一个个pred
# # start_time = time.time()
# # y_pred2 = []
# # for vector in X_test:
# #     pred_temp = model.predict(vector.reshape(1, -1))
# #     y_pred2.append(pred_temp)
# # print ('predict(RF-一个个) took %fs!' % (time.time() - start_time))
# print('Final prediction score with optimal(use RF): [%.8f]' % accuracy_score(y_test, y_pred))
# #
print ('******************** 5.KNN *********************' )
model = KNeighborsClassifier(n_jobs=1)
start_time = time.time()
model.fit(X_train, y_train)
print ('training(KNN) took %fs!' % (time.time() - start_time))
start_time = time.time()
#直接pred
y_pred = model.predict(X_test)
print ('predict(KNN-直接) took %fs!' % (time.time() - start_time))
#一个个pred
start_time = time.time()
y_pred2 = []
for vector in X_test:
    pred_temp = model.predict(vector.reshape(1, -1))
    y_pred2.append(pred_temp)
print ('predict(KNN-一个个) took %fs!' % (time.time() - start_time))
print('Final prediction score with optimal(use KNN): [%.8f]' % accuracy_score(y_test, y_pred2))

print ('******************** 6.Superknn *********************' )
# Make train/test split
# As usual in machine learning task we have X_train, y_train, and X_test
# X_train, X_test, y_train, y_test = train_test_split(X, y, shuffle=True, test_size=0.2, random_state=0)

models = [ExtraTreesClassifier(n_jobs=7, random_state=0),
    RandomForestClassifier(n_jobs=7, random_state=0),

    XGBClassifier(n_jobs=7, random_state=0)
]

'''方法1：不适用need_proba,这时只有参数alpha,所以不需要调用opmimal_parameters函数'''
# model = Superknn(models=models, metric=accuracy_score, regression=False, random_state=0,
#                  folds=5, needs_proba=False, shuffle=True)
# # #一般来说，此时参数比较少可以一个个试alpha，不需要调用optimal_parametersh函数
# model.opmimal_parameters(X=X_train, y=y_train, eps_size=0.2)
# alphas = np.arange(2, 5, 1)
# for alpha in alphas:
#     print('alpha:', alpha)
#     model.alpha = alpha
#     model.fit(X_train, y_train)
#     y_pred = model.predict(X_test)
#     # print(y_test)
#     # print(y_pred)
#     print('Final prediction score with optimal(use superknn): [%.8f]' % accuracy_score(y_test, y_pred))

'''方法2：使用need_proba,这时只有参数alpha和epsilon，其中alpha一共有models*classes = 4*10=40个,
grid_points = 10个
所以先调用opmimal_parameters函数先找到一个大概的范围'''
model = Superknn(models=models, metric=accuracy_score, regression=False, random_state=0,
                 folds=5, needs_proba=True, n_jobs=20, shuffle=True)
# # model.opmimal_parameters(X=X_train, y=y_train, eps_size=0.1, grid_points=10)
#
#
epsilons = np.arange(0.1, 0.5, 0.01)
model.alpha = 20
for epsilon in epsilons:
    model.epsilon = epsilon

    model.fit(X_train, y_train)

    y_pred = model.predict(X_test)
    print('Final prediction score with optimal(use superknn): [%.8f]\n' % accuracy_score(y_test, y_pred))

max_stars_repo_name: DAIM-ML/autotf

max_issues_repo_name: DAIM-ML/autotf

max_forks_repo_name: DAIM-ML/autotf


text:
"""
Modul, welches die Wordcloud-Image-Funktionen beinhaltet.
"""
import numbers
import random

import numpy as np
from wordcloud import WordCloud

from visuanalytics.analytics.control.procedures.step_data import StepData
from visuanalytics.util import resources

WORDCLOUD_DEFAULT_PARAMETER = {
    "background_color": "white",
    "width": 400,
    "height": 200,
    "collocations": True,
    "max_font_size": None,
    "max_words": 200,
    "contour_width": 0,
    "contour_color": "white",
    "font_path": "fonts/Dosis-Regular.ttf",
    "prefer_horizontal": 0.90,
    "scale": 1,
    "min_font_size": 4,
    "font_step": 1,
    "mode": "RGB",
    "relative_scaling": 0.5,
    "color_func": False,
    "regexp": None,
    "colormap": "viridis",
    "normalize_plurals": True,
    "stopwords": None,
    "repeat": False,
    "mask": None
}
"""Ein Dictionary bestehend aus allen Default-Parametern für die Generierung einer Wordcloud.
   Siehe auch Quelle: https://www.datacamp.com/community/tutorials/wordcloud-python"""

DEFAULT_COLOR_FUNC_VALUES = [0, 0, 0, 40]
"""Ein Array bestehend aus den Default-Werten für die Funktion color_func. Dieses Array stellt die Farbwerte (HSL-Farbschema)
    zur Generierung eines Farbverlaufs mit Schwarz- & Grautönen dar. """

def get_color_func(h, s, l_start, l_end):
    """
    Erstellt das Farbspektrum (mit HSL-Farbschema), in welchem die Wörter der Wordcloud dargestellt werden sollen.

    Die Werte werden jeweils als Integer angegeben.

    Quelle: https://github.com/amueller/word_cloud/blob/master/wordcloud/wordcloud.py (siehe random_color_func)
    Funktion wurde für unsere Zwecke angepasst.

    :param h: (hue) Farbton in Grad auf einem Farbenrad. Mögliche Werte: 0 bis 360
        0 ist Rot, 120 ist Grün und 240 ist Blau.
    :param s: (saturation) Sättigung in Prozent. 0% entspricht einem Grauton, 100% ist die volle Farbe.
    :param l_start: (lightness range start) Start des Helligkeitsbereichs in Prozent: 0% ist Schwarz, 100% ist Weiß.
    :param l_end: (lightness range end) Ende des Helligkeitsbereichs in Prozent: 0% ist Schwarz, 100% ist Weiß.
    :return: Funktion, die den Farbverlauf innerhalb der Wordcloud erstellt
    :rtype: callable
    """

    def color_func(word, font_size, position, orientation, random_state=None, **kwargs):
        return "hsl(%d, %d%%, %d%%)" % (h, s, random.randint(l_start, l_end))

    return color_func

def wordcloud(values: dict, step_data: StepData, prev_paths):
    """Erstellt ein Wordcloud-Bild.

    Der Standard-Farbverlauf bei color_func true ist Grau/Schwarz.
    Die Standard-Farbe ist generell die colormap viridis.

    Quelle bzgl. Verwendung: https://github.com/amueller/word_cloud

    :param values: Image Bauplan des zu erstellenden Bildes
    :param step_data: Daten aus der API
    :param prev_paths: Alle Image Baupläne und somit auch alle Pfade zu den bisher erstellen Bildern
    :return: Den Pfad zum erstellten Bild
    :rtype: str
    """
    wordcloud_parameter = dict(WORDCLOUD_DEFAULT_PARAMETER)
    parameter = values.get("parameter", {})

    for param in parameter:
        if param in wordcloud_parameter:
            if isinstance(wordcloud_parameter[param], bool):
                value = step_data.get_data(parameter[param], {}, bool)
            elif isinstance(wordcloud_parameter[param], numbers.Number):
                value = step_data.get_data(parameter[param], {}, numbers.Number)
            else:
                value = step_data.format(parameter[param])

            wordcloud_parameter[param] = value

    path = resources.get_resource_path(wordcloud_parameter["font_path"])
    wordcloud_parameter["font_path"] = path

    if bool(wordcloud_parameter.get("color_func", False)):
        cfw = list(DEFAULT_COLOR_FUNC_VALUES)

        if "color_func_words" in values:
            cfw_list = step_data.format(values["color_func_words"]).split(" ")

            for idx, c in enumerate(cfw_list):
                cfw[idx] = int(c)

        wordcloud_parameter["color_func"] = get_color_func(*cfw)
    else:
        wordcloud_parameter["color_func"] = None

    if parameter.get("colormap", ""):
        wordcloud_parameter["colormap"] = step_data.format(parameter["colormap"])

    if parameter.get("figure", None) is not None:
        figure = step_data.format(parameter["figure"], {})
        if figure == "circle":
            x0 = wordcloud_parameter["width"]
            y0 = wordcloud_parameter["height"]
            x, y = np.ogrid[:x0, :y0]

            mask = (x - (x0 / 2)) ** 2 + (y - (y0 / 2)) ** 2 > 500 ** 2

            mask = 255 * mask.astype(int)
            wordcloud_parameter["mask"] = mask

    if values.get("use_stopwords", None) is not None:
        try:
            file = resources.get_resource_path("stopwords/stopwords.txt")
            with open(file, "r", encoding='utf-8') as f:
                list_stopwords = f.read().splitlines()
        except IOError:
            list_stopwords = []

        dont_use = step_data.get_data(values.get("stopwords", []), {}, list)
        for each in list_stopwords:
            if each not in dont_use:
                dont_use.append(each)
        wordcloud_parameter["stopwords"] = set(dont_use)
    else:
        dont_use = step_data.get_data(values.get("stopwords", []), {}, list)
        wordcloud_parameter["stopwords"] = set(dont_use)
    if values.get("text", None) is not None:
        wordcloud_image = WordCloud(**wordcloud_parameter).generate(step_data.format(values["text"], {}))
    elif values.get("dict", None) is not None:
        wordcloud_image = WordCloud(**wordcloud_parameter).generate_from_frequencies(
            step_data.get_data(values["dict"], {}, dict))

    image = wordcloud_image.to_image()
    file = resources.new_temp_resource_path(step_data.data["_pipe_id"], "png")
    image.save(file)
    return file

max_stars_repo_name: mxsph/Data-Analytics

max_issues_repo_name: mxsph/Data-Analytics

max_forks_repo_name: visuanalytics/visuanalytics


text:
# mv ~/Downloads/input.txt input-15.txt

import sys
from collections import defaultdict
from itertools import *
from more_itertools import *
# import numpy as np
# import scipy.signal
from heapq import *

lines = [line.strip() for line in open("input-15.txt")]

risks = [list(map(int, list(line))) for line in lines]
width = len(risks[0])
height = len(risks)

new_risks = []
for y in range(height*5):
    row = []
    for x in range(width*5):
        value = risks[y % height][x % width] + (y // height) + (x // width)
        if value > 9:
            value -= 9

        row.append(value)

    new_risks.append(row)
risks = new_risks
width *= 5
height *= 5

visited = [[False]*width for i in range(height)]
dist = [[9999999]*width for i in range(height)]

dist[0][0] = 0
h = []

def check(x, y, fr):
    if x < 0 or y < 0 or x >= width or y >= height or visited[y][x]:
        return

    new = fr + risks[y][x]
    if new < dist[y][x]:
        dist[y][x] = new
        heappush(h, (new, x, y))

for x in range(width):
    for y in range(height):
        heappush(h, (99999, x, y) )
heappush(h, (0, 0, 0) )

x = 0
y = 0
while True:
    while True:
        value, x, y = heappop(h)
        if not visited[y][x]:
            break

    visited[y][x] = True
    check(x + 1, y, dist[y][x])
    check(x - 1, y, dist[y][x])
    check(x, y - 1, dist[y][x])
    check(x, y + 1, dist[y][x])

    if y == height - 1 and x == width - 1:
        print(dist[y][x])
        sys.exit(0)

max_stars_repo_name: lkesteloot/advent-of-code

max_issues_repo_name: lkesteloot/advent-of-code

max_forks_repo_name: lkesteloot/advent-of-code

reason: import numpy,import scipy


text:
"""Example showing how to use the ProblemSpec approach.

Showcases method chaining, and parallel model runs using 2 processors.

The following caveats apply:

1. Functions passed into `sample`, `analyze` and `evaluate` must
   accept a numpy array of `X` values as the first parameter, and return a
   numpy array of results.
2. Parallel evaluation/analysis is only beneficial for long-running models
   or large datasets
3. Currently, results must fit in memory - no on-disk caching is provided.
"""

from SALib.analyze import sobol
from SALib.sample import saltelli
from SALib.test_functions import Ishigami
from SALib import ProblemSpec
import numpy as np

import time

if __name__ == '__main__':

    # Create the SALib Problem specification
    sp = ProblemSpec({
        'names': ['x1', 'x2', 'x3'],
        'groups': None,
        'bounds': [[-np.pi, np.pi]]*3,
        'outputs': ['Y']
    })

    # Single core example
    start = time.perf_counter()
    (sp.sample_saltelli(2**15)
        .evaluate(Ishigami.evaluate)
        .analyze_sobol(calc_second_order=True, conf_level=0.95))
    print("Time taken with 1 core:", time.perf_counter() - start, '\n')

    # Same above example, but passing in specific functions
    # (sp.sample(saltelli.sample, 25000, calc_second_order=True)
    #     .evaluate(Ishigami.evaluate)
    #     .analyze(sobol.analyze, calc_second_order=True, conf_level=0.95))

    # Samples, model results and analyses can be extracted:
    # print(sp.samples)
    # print(sp.results)
    # print(sp.analysis)
    # print(sp.to_df())

    # Can set pre-existing samples/results as needed
    # sp.samples = some_numpy_array
    # sp.set_samples(some_numpy_array)
    #
    # Using method chaining...
    # (sp.set_samples(some_numpy_array)
    #    .set_results(some_result_array)
    #    .analyze_sobol(calc_second_order=True, conf_level=0.95))

    # Parallel example
    start = time.perf_counter()
    (sp.sample(saltelli.sample, 2**15)
         # can specify number of processors to use with `nprocs`
         # this will be capped to the number of detected processors
         # or, in the case of analysis, the number of outputs.
        .evaluate(Ishigami.evaluate, nprocs=2)
        .analyze_sobol(calc_second_order=True, conf_level=0.95, nprocs=2))
    print("Time taken with 2 cores:", time.perf_counter() - start, '\n')

    print(sp)

    # Distributed example
    # Specify itself as servers as an example
    servers = ('localhost:55774',
               'localhost:55775',
               'localhost:55776')

    start = time.perf_counter()
    (sp.sample(saltelli.sample, 2**15)
        .evaluate_distributed(Ishigami.evaluate, nprocs=2, servers=servers, verbose=True)
        .analyze(sobol.analyze, calc_second_order=True, conf_level=0.95))
    print("Time taken with distributed cores:", time.perf_counter() - start, '\n')

    print(sp)

    # To display plots:
    # import matplotlib.pyplot as plt
    # sp.plot()
    # plt.show()

# First-order indices expected with Saltelli sampling:
# x1: 0.3139
# x2: 0.4424
# x3: 0.0

max_stars_repo_name: zjzh/SALib

max_issues_repo_name: zjzh/SALib

max_forks_repo_name: zjzh/SALib


text:
"""General matrix math functions."""
# Author: Ben Gravell

import numpy as np
from numpy import linalg as la
from scipy.linalg import solve_discrete_lyapunov, solve_discrete_are
from functools import reduce

from extramath import quadratic_formula

def vec(A):
    """Return the vectorized matrix A by stacking its columns."""
    return A.reshape(-1, order="F")

def svec(A):
    """Return the symmetric vectorization i.e. the vectorization of the upper triangular part of matrix A."""
    return A[np.triu_indices(A.shape[0])]

def svec2(A):
    """
    Return the symmetric vectorization i.e. the vectorization of the upper triangular part of matrix A
    with off-diagonal entries multiplied by sqrt(2) so that la.norm(A, ord='fro')**2 == np.dot(svec2(A), svec2(A))
    """
    B = A + np.triu(A, 1)*(2**0.5 - 1)
    return B[np.triu_indices(A.shape[0])]

def smat(v):
    """Return the symmetric matricization i.e. the inverse operation of svec of vector v."""
    m = v.size
    n = int(((1+m*8)**0.5 - 1)/2)
    idx_upper = np.triu_indices(n)
    idx_lower = np.tril_indices(n, -1)
    A = np.zeros([n, n])
    A[idx_upper] = v
    A[idx_lower] = A.T[idx_lower]
    return A

def smat2(v):
    """Return the symmetric matricization i.e. the inverse operation of svec2 of vector v."""
    m = v.size
    n = int(((1+m*8)**0.5 - 1)/2)
    idx_upper = np.triu_indices(n)
    idx_lower = np.tril_indices(n, -1)
    A = np.zeros([n, n])
    A[idx_upper] = v
    A[np.triu_indices(n,1)] /= 2**0.5
    A[idx_lower] = A.T[idx_lower]
    return A

def sympart(A):
    """Return the symmetric part of matrix A."""
    return 0.5*(A+A.T)

def is_pos_def(A):
    """Check if matrix A is positive definite."""
    try:
        la.cholesky(A)
        return True
    except np.linalg.LinAlgError:
        return False

def succ(A,B):
    """Check the positive definite partial ordering of A > B."""
    return is_pos_def(A-B)

def psdpart(X):
    """Return the positive semidefinite part of a symmetric matrix."""
    X = sympart(X)
    Y = np.zeros_like(X)
    eigvals, eigvecs = la.eig(X)
    for i in range(X.shape[0]):
        if eigvals[i] > 0:
            Y += eigvals[i]*np.outer(eigvecs[:,i],eigvecs[:,i])
    Y = sympart(Y)
    return Y

def kron(*args):
    """Overload and extend the numpy kron function to take a single argument."""
    if len(args)==1:
        return np.kron(args[0], args[0])
    else:
        return np.kron(*args)

def mdot(*args):
    """Multiple dot product."""
    return reduce(np.dot, args)

def mip(A,B):
    """Matrix inner product of A and B."""
    return np.trace(mdot(A.T, B))

def specrad(A):
    """Spectral radius of matrix A."""
    try:
        return np.max(np.abs(la.eig(A)[0]))
    except np.linalg.LinAlgError:
        return np.nan

def printeigs(A):
    """Print all eigenvalues of matrix A."""
    print(la.eig(A)[0])
    return

def minsv(A):
    """Minimum singular value."""
    return la.svd(A)[1].min()

def solveb(a, b):
    """
    Solve a = bx.
    Similar to MATLAB / operator for square invertible matrices.
    """
    return la.solve(b.T, a.T).T

def lstsqb(a, b):
    """
    Return least-squares solution to a = bx.
    Similar to MATLAB / operator for rectangular matrices.
    If b is invertible then the solution is la.solve(a, b).T
    """
    return la.lstsq(b.T, a.T, rcond=None)[0].T

def dlyap(A, Q):
    """
    Solve the discrete-time Lyapunov equation.
    Wrapper around scipy.linalg.solve_discrete_lyapunov.
    Pass a copy of input matrices to protect them from modification.
    """
    try:
        return solve_discrete_lyapunov(np.copy(A), np.copy(Q))
    except ValueError:
        return np.full_like(Q, np.inf)

def dare_scalar(A, B, Q, R):
    """
    Solve the discrete-time algebraic Riccati equation for the scalar case of
    a single state and a single input.
    In this case the equation is a scalar quadratic equation.
    """
    A, B, Q, R = [float(var) for var in [A, B, Q, R]]

    A2 = A**2
    B2 = B**2

    aa = -B2
    bb = R*(A2-1) + Q*B2
    cc = Q*R

    roots = np.array(quadratic_formula(aa, bb, cc))

    if not(roots[0] > 0 or roots[1] > 0):
        P = None
    else:
        P = roots[roots > 0][0]*np.eye(1)
    return P

def dare_gain_scalar(A, B, Q, R):
    P = dare_scalar(A, B, Q, R)
    K = -B*P*A / (R+B*B*P)
    return P, K

def dare(A, B, Q, R):
    """
    Solve the discrete-time algebraic Riccati equation.
    Wrapper around scipy.linalg.solve_discrete_are.
    Pass a copy of input matrices to protect them from modification.
    """
    # Handle the scalar case more efficiently by solving the ARE exactly
    if A.shape[1] == 1 and B.shape[1] == 1:
        return dare_scalar(A, B, Q, R)*np.eye(1)
    else:
        return solve_discrete_are(np.copy(A), np.copy(B), np.copy(Q), np.copy(R))

def dare_gain(A, B, Q, R):
    """
    Solve the discrete-time algebraic Riccati equation.
    Return the optimal cost-to-go matrix P and associated gains K
    such that u = Kx is the optimal control.
    """
    # Handle the scalar case more efficiently by solving the ARE exactly
    if A.shape[1] == 1 and B.shape[1] == 1:
        P, K = dare_gain_scalar(A, B, Q, R)
        return P*np.eye(1), K*np.eye(1)
    else:
        P = dare(A, B, Q, R)
        K = -la.solve((R + mdot(B.T, P, B)),  mdot(B.T, P, A))
        return P, K

reason: import numpy,from numpy,from scipy


text:
import glob
import pickle

import cv2
import numpy as np

# The directory you downloaded CIFAR-10
# You can download cifar10 data via https://www.kaggle.com/janzenliu/cifar-10-batches-py
data_dir = './data'

all_files = glob.glob(data_dir + '/data_batch' + '*')
test_files = glob.glob(data_dir + '/test_batch' + '*')

all_files = all_files + test_files

def unpickle(file):
    with open(file, 'rb') as fo:
        dict = pickle.load(fo, encoding='bytes')
    return dict

# labal
# 0:airplane, 1:automobile, 2:bird. 3:cat, 4:deer, 5:dog, 6:frog, 7:horse, 8:ship, 9:truck
all_image = []
all_label = []

print(all_files)

for file in all_files:
    print(file)
    ret = unpickle(file)

    for i, arr in enumerate(ret[b'data']):
        img = np.reshape(arr, (3, 32, 32))
        img = img.transpose(1, 2, 0)
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        all_image.append(img)
        all_label.append(ret[b'labels'][i])

all_images = np.array(all_image)
all_labels = np.array(all_label)
print(all_images.shape)

np.savez(data_dir + '/' + 'cifar10.npz', images=all_images, labels=all_labels)


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

from __future__ import unicode_literals

import itertools
import numpy as np
from math import factorial

def _embed(x, order=3, delay=1):
    """Time-delay embedding.

    Parameters
    ----------
    x : 1d-array, shape (n_times)
        Time series
    order : int
        Embedding dimension (order)
    delay : int
        Delay.

    Returns
    -------
    embedded : ndarray, shape (n_times - (order - 1) * delay, order)
        Embedded time-series.
    """
    N = len(x)
    Y = np.empty((order, N - (order - 1) * delay))
    for i in range(order):
        Y[i] = x[i * delay:i * delay + Y.shape[1]]
    return Y.T

def util_pattern_space(time_series, lag, dim):
    """Create a set of sequences with given lag and dimension

    Args:
       time_series: Vector or string of the sample data
       lag: Lag between beginning of sequences
       dim: Dimension (number of patterns)

    Returns:
        2D array of vectors
    """
    n = len(time_series)

    if lag * dim > n:
        raise Exception('Result matrix exceeded size limit, try to change lag or dim.')
    elif lag < 1:
        raise Exception('Lag should be greater or equal to 1.')

    pattern_space = np.empty((n - lag * (dim - 1), dim))
    for i in range(n - lag * (dim - 1)):
        for j in range(dim):
            pattern_space[i][j] = time_series[i + j * lag]

    return pattern_space

def util_standardize_signal(time_series):
    return (time_series - np.mean(time_series)) / np.std(time_series)

def util_granulate_time_series(time_series, scale):
    """Extract coarse-grained time series

    Args:
        time_series: Time series
        scale: Scale factor

    Returns:
        Vector of coarse-grained time series with given scale factor
    """
    n = len(time_series)
    b = int(np.fix(n / scale))
    temp = np.reshape(time_series[0:b*scale], (b, scale))
    cts = np.mean(temp, axis = 1)
    return cts

def util_rolling_window(a, window):
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

def shannon_entropy(time_series):
    """Return the Shannon Entropy of the sample data.

    Args:
        time_series: Vector or string of the sample data

    Returns:
        The Shannon Entropy as float value
    """

    # Check if string
    if not isinstance(time_series, str):
        time_series = list(time_series)

    # Create a frequency data
    data_set = list(set(time_series))
    freq_list = []
    for entry in data_set:
        counter = 0.
        for i in time_series:
            if i == entry:
                counter += 1
        freq_list.append(float(counter) / len(time_series))

    # Shannon entropy
    ent = 0.0
    for freq in freq_list:
        ent += freq * np.log2(freq)
    ent = -ent
    return ent

def sample_entropy(time_series, sample_length, tolerance = None):
    """Calculates the sample entropy of degree m of a time_series.

    This method uses chebychev norm.
    It is quite fast for random data, but can be slower is there is
    structure in the input time series.

    Args:
        time_series: numpy array of time series
        sample_length: length of longest template vector
        tolerance: tolerance (defaults to 0.1 * std(time_series)))
    Returns:
        Array of sample entropies:
            SE[k] is ratio "#templates of length k+1" / "#templates of length k"
            where #templates of length 0" = n*(n - 1) / 2, by definition
    Note:
        The parameter 'sample_length' is equal to m + 1 in Ref[1].

    References:
        [1] http://en.wikipedia.org/wiki/Sample_Entropy
        [2] http://physionet.incor.usp.br/physiotools/sampen/
        [3] Madalena Costa, Ary Goldberger, CK Peng. Multiscale entropy analysis
            of biological signals
    """
    #The code below follows the sample length convention of Ref [1] so:
    M = sample_length - 1;

    time_series = np.array(time_series)
    if tolerance is None:
        tolerance = 0.1*np.std(time_series)

    n = len(time_series)

    #Ntemp is a vector that holds the number of matches. N[k] holds matches templates of length k
    Ntemp = np.zeros(M + 2)
    #Templates of length 0 matches by definition:
    Ntemp[0] = n*(n - 1) / 2

    for i in range(n - M - 1):
        template = time_series[i:(i+M+1)];#We have 'M+1' elements in the template
        rem_time_series = time_series[i+1:]

        searchlist = np.nonzero(np.abs(rem_time_series - template[0]) < tolerance)[0]

        go = len(searchlist) > 0;

        length = 1;

        Ntemp[length] += len(searchlist)

        while go:
            length += 1
            nextindxlist = searchlist + 1;
            nextindxlist = nextindxlist[nextindxlist < n - 1 - i]#Remove candidates too close to the end
            nextcandidates = rem_time_series[nextindxlist]
            hitlist = np.abs(nextcandidates - template[length-1]) < tolerance
            searchlist = nextindxlist[hitlist]

            Ntemp[length] += np.sum(hitlist)

            go = any(hitlist) and length < M + 1

    sampen =  - np.log(Ntemp[1:] / Ntemp[:-1])
    return sampen

def multiscale_entropy(time_series, sample_length, tolerance = None, maxscale = None):
    """Calculate the Multiscale Entropy of the given time series considering
    different time-scales of the time series.

    Args:
        time_series: Time series for analysis
        sample_length: Bandwidth or group of points
        tolerance: Tolerance (default = 0.1*std(time_series))

    Returns:
        Vector containing Multiscale Entropy

    Reference:
        [1] http://en.pudn.com/downloads149/sourcecode/math/detail646216_en.html
    """

    if tolerance is None:
        #we need to fix the tolerance at this level. If it remains 'None' it will be changed in call to sample_entropy()
        tolerance = 0.1*np.std(time_series)
    if maxscale is None:
        maxscale = len(time_series)

    mse = np.zeros(maxscale)

    for i in range(maxscale):
        temp = util_granulate_time_series(time_series, i+1)
        mse[i] = sample_entropy(temp, sample_length, tolerance)[-1]
    return mse

def permutation_entropy(time_series, order=3, delay=1, normalize=False):
    """Permutation Entropy.

    Parameters
    ----------
    time_series : list or np.array
        Time series
    order : int
        Order of permutation entropy
    delay : int
        Time delay
    normalize : bool
        If True, divide by log2(factorial(m)) to normalize the entropy
        between 0 and 1. Otherwise, return the permutation entropy in bit.

    Returns
    -------
    pe : float
        Permutation Entropy

    References
    ----------
    .. [1] Massimiliano Zanin et al. Permutation Entropy and Its Main
        Biomedical and Econophysics Applications: A Review.
        http://www.mdpi.com/1099-4300/14/8/1553/pdf

    .. [2] Christoph Bandt and Bernd Pompe. Permutation entropy — a natural
        complexity measure for time series.
        http://stubber.math-inf.uni-greifswald.de/pub/full/prep/2001/11.pdf

    Notes
    -----
    Last updated (Oct 2018) by Raphael Vallat (raphaelvallat9@gmail.com):
    - Major speed improvements
    - Use of base 2 instead of base e
    - Added normalization

    Examples
    --------
    1. Permutation entropy with order 2

        >>> x = [4, 7, 9, 10, 6, 11, 3]
        >>> # Return a value between 0 and log2(factorial(order))
        >>> print(permutation_entropy(x, order=2))
            0.918

    2. Normalized permutation entropy with order 3

        >>> x = [4, 7, 9, 10, 6, 11, 3]
        >>> # Return a value comprised between 0 and 1.
        >>> print(permutation_entropy(x, order=3, normalize=True))
            0.589
    """
    x = np.array(time_series)
    hashmult = np.power(order, np.arange(order))
    # Embed x and sort the order of permutations
    sorted_idx = _embed(x, order=order, delay=delay).argsort(kind='quicksort')
    # Associate unique integer to each permutations
    hashval = (np.multiply(sorted_idx, hashmult)).sum(1)
    # Return the counts
    _, c = np.unique(hashval, return_counts=True)
    # Use np.true_divide for Python 2 compatibility
    p = np.true_divide(c, c.sum())
    pe = -np.multiply(p, np.log2(p)).sum()
    if normalize:
        pe /= np.log2(factorial(order))
    return pe

def multiscale_permutation_entropy(time_series, m, delay, scale):
    """Calculate the Multiscale Permutation Entropy

    Args:
        time_series: Time series for analysis
        m: Order of permutation entropy
        delay: Time delay
        scale: Scale factor

    Returns:
        Vector containing Multiscale Permutation Entropy

    Reference:
        [1] Francesco Carlo Morabito et al. Multivariate Multi-Scale Permutation Entropy for
            Complexity Analysis of Alzheimer’s Disease EEG. www.mdpi.com/1099-4300/14/7/1186
        [2] http://www.mathworks.com/matlabcentral/fileexchange/37288-multiscale-permutation-entropy-mpe/content/MPerm.m
    """
    mspe = []
    for i in range(scale):
        coarse_time_series = util_granulate_time_series(time_series, i + 1)
        pe = permutation_entropy(coarse_time_series, order=m, delay=delay)
        mspe.append(pe)
    return mspe

# TODO add tests
def composite_multiscale_entropy(time_series, sample_length, scale, tolerance=None):
    """Calculate the Composite Multiscale Entropy of the given time series.

    Args:
        time_series: Time series for analysis
        sample_length: Number of sequential points of the time series
        scale: Scale factor
        tolerance: Tolerance (default = 0.1...0.2 * std(time_series))

    Returns:
        Vector containing Composite Multiscale Entropy

    Reference:
        [1] Wu, Shuen-De, et al. "Time series analysis using
            composite multiscale entropy." Entropy 15.3 (2013): 1069-1084.
    """
    cmse = np.zeros((1, scale))

    for i in range(scale):
        for j in range(i):
            tmp = util_granulate_time_series(time_series[j:], i + 1)
            cmse[i] += sample_entropy(tmp, sample_length, tolerance) / (i + 1)
    return cmse

def weighted_permutation_entropy(time_series, order=3, normalize=False):
    """Calculate the Weighted Permutation Entropy.

    Weighted permutation entropy is based on the regular permutation entropy,
    but puts additional weight on those windows that show a high variability
    in the initial time series.

    Parameters
    ----------
    time_series : list or np.array
        Time series
    order : int
        Order of permutation entropy
    normalize : bool
        If True, divide by log2(factorial(m)) to normalize the entropy
        between 0 and 1. Otherwise, return the permutation entropy in bit.
    Returns
    -------
    pe : float
        Weighted Permutation Entropy
    References
    ----------
    .. [1] Bilal Fadlallah et al. Weighted-permutation entropy: A complexity
    measure for time series incorporating amplitude information
    https://link.aps.org/accepted/10.1103/PhysRevE.87.022911
    """
    x = np.array(time_series)
    hashmult = np.power(order, np.arange(order))
    # Embed x and sort the order of permutations

    embedded = _embed(x, order=order)
    sorted_idx = embedded.argsort(kind='quicksort')
    weights = np.var(util_rolling_window(x, order), 1)
    hashval = (np.multiply(sorted_idx, hashmult)).sum(1)
    mapping = {}
    for i in np.unique(hashval):
        mapping[i] = np.where(hashval == i)[0]
    weighted_counts = dict.fromkeys(mapping)
    for k, v in mapping.items():
        weighted_count = 0
        for i in v:
            weighted_count += weights[i]
        weighted_counts[k] = weighted_count
    # Associate unique integer to each permutations
    # Return the counts
    # Use np.true_divide for Python 2 compatibility
    weighted_counts_array = np.array(list(weighted_counts.values()))
    p = np.true_divide(weighted_counts_array, weighted_counts_array.sum())
    pe = -np.multiply(p, np.log2(p)).sum()
    if normalize:
        pe /= np.log2(factorial(order))
    return pe

max_stars_repo_name: julesjulian/pyEntropy

max_issues_repo_name: julesjulian/pyEntropy

max_forks_repo_name: julesjulian/pyEntropy


text:
#coding=utf-8
#Word2Vec.py
import collections
import math
import os
import random
import zipfile
import numpy as np
import urllib
import tensorflow as tf
import matplotlib.pyplot as plt
#定义下载文本数据的函数，使用urllib.request.urlretrieve
url = 'http://mattmahoney.net/dc/'

def maybe_download(filename, expected_bytes):
    if not os.path.exists(filename):
        filename, _ = urllib.request.urlretrieve(url + filename, filename)
    statinfo = os.stat(filename)
    if statinfo.st_size == expected_bytes:
        print('Found and verified', filename)
    else:
        print(statinfo.st_size)
        raise Exception('Failed to verify' + filename + '. Can you get to it with a browser?')
    return filename

filename = maybe_download('text8.zip', 31344016)

max_stars_repo_name: KangChou/NLPs

max_issues_repo_name: KangChou/NLPs

max_forks_repo_name: KangChou/NLPs


text:
# -*- coding: utf-8 -*-
from distutils.version import StrictVersion
from urllib.error import URLError

import numpy as np
from qe_tools import CONSTANTS
from xmlschema import XMLSchema

from aiida_quantumespresso.utils.mapping import get_logging_container

from .exceptions import XMLParseError
from .versions import get_default_schema_filepath, get_schema_filepath

def raise_parsing_error(message):
    raise XMLParseError(message)

def parser_assert(condition, message, log_func=raise_parsing_error):
    if not condition:
        log_func(message)

def parser_assert_equal(val1, val2, message, log_func=raise_parsing_error):
    if not (val1 == val2):
        msg = f'Violated assertion: {val1} == {val2}'
        if message:
            msg += ' - '
            msg += message
        log_func(msg)

def cell_volume(a1, a2, a3):
    r"""Returns the volume of the primitive cell: :math:`|\vec a_1\cdot(\vec a_2\cross \vec a_3)|`"""
    a_mid_0 = a2[1] * a3[2] - a2[2] * a3[1]
    a_mid_1 = a2[2] * a3[0] - a2[0] * a3[2]
    a_mid_2 = a2[0] * a3[1] - a2[1] * a3[0]

    return abs(float(a1[0] * a_mid_0 + a1[1] * a_mid_1 + a1[2] * a_mid_2))

def parse_xml_post_6_2(xml):
    """Parse the content of XML output file written by `pw.x` and `cp.x` with the new schema-based XML format.

    :param xml: parsed XML
    :returns: tuple of two dictionaries, with the parsed data and log messages, respectively
    """
    e_bohr2_to_coulomb_m2 = 57.214766  # e/a0^2 to C/m^2 (electric polarization) from Wolfram Alpha

    logs = get_logging_container()

    schema_filepath = get_schema_filepath(xml)

    try:
        xsd = XMLSchema(schema_filepath)
    except URLError:

        # If loading the XSD file specified in the XML file fails, we try the default
        schema_filepath_default = get_default_schema_filepath()

        try:
            xsd = XMLSchema(schema_filepath_default)
        except URLError:
            raise XMLParseError(
                f'Could not open or parse the XSD files {schema_filepath} and {schema_filepath_default}'
            )
        else:
            schema_filepath = schema_filepath_default

    # Validate XML document against the schema
    # Returned dictionary has a structure where, if tag ['key'] is "simple", xml_dictionary['key'] returns its content.
    # Otherwise, the following keys are available:
    #
    #  xml_dictionary['key']['$'] returns its content
    #  xml_dictionary['key']['@attr'] returns its attribute 'attr'
    #  xml_dictionary['key']['nested_key'] goes one level deeper.

    # Fix a bug of QE 6.8: the output XML is not consistent with schema, see
    # https://github.com/aiidateam/aiida-quantumespresso/pull/717
    xml_creator = xml.find('./general_info/creator')
    if xml_creator is not None and 'VERSION' in xml_creator.attrib:
        creator_version = xml_creator.attrib['VERSION']
        if creator_version == '6.8':
            root = xml.getroot()
            timing_info = root.find('./timing_info')
            partial_pwscf = timing_info.find("partial[@label='PWSCF'][@calls='0']")
            try:
                timing_info.remove(partial_pwscf)
            except (TypeError, ValueError):
                pass

    xml_dictionary, errors = xsd.to_dict(xml, validation='lax')
    if errors:
        logs.error.append(f'{len(errors)} XML schema validation error(s) schema: {schema_filepath}:')
        for err in errors:
            logs.error.append(str(err))

    xml_version = StrictVersion(xml_dictionary['general_info']['xml_format']['@VERSION'])
    inputs = xml_dictionary.get('input', {})
    outputs = xml_dictionary['output']

    lattice_vectors = [
        [x * CONSTANTS.bohr_to_ang for x in outputs['atomic_structure']['cell']['a1']],
        [x * CONSTANTS.bohr_to_ang for x in outputs['atomic_structure']['cell']['a2']],
        [x * CONSTANTS.bohr_to_ang for x in outputs['atomic_structure']['cell']['a3']],
    ]

    has_electric_field = inputs.get('electric_field', {}).get('electric_potential', None) == 'sawtooth_potential'
    has_dipole_correction = inputs.get('electric_field', {}).get('dipole_correction', False)

    if 'occupations' in inputs.get('bands', {}):
        try:
            occupations = inputs['bands']['occupations']['$']  # yapf: disable
        except TypeError:  # "string indices must be integers" -- might have attribute 'nspin'
            occupations = inputs['bands']['occupations']
    else:
        occupations = None

    starting_magnetization = []
    magnetization_angle1 = []
    magnetization_angle2 = []

    for specie in outputs['atomic_species']['species']:
        starting_magnetization.append(specie.get('starting_magnetization', 0.0))
        magnetization_angle1.append(specie.get('magnetization_angle1', 0.0))
        magnetization_angle2.append(specie.get('magnetization_angle2', 0.0))

    constraint_mag = 0
    spin_constraints = inputs.get('spin_constraints', {}).get('spin_constraints', None)
    if spin_constraints == 'atomic':
        constraint_mag = 1
    elif spin_constraints == 'atomic direction':
        constraint_mag = 2
    elif spin_constraints == 'total':
        constraint_mag = 3
    elif spin_constraints == 'total direction':
        constraint_mag = 6

    lsda = inputs.get('spin', {}).get('lsda', False)
    spin_orbit_calculation = inputs.get('spin', {}).get('spinorbit', False)
    non_colinear_calculation = outputs['magnetization']['noncolin']
    do_magnetization = outputs['magnetization'].get('do_magnetization', False)

    # Time reversal symmetry of the system
    if non_colinear_calculation and do_magnetization:
        time_reversal = False
    else:
        time_reversal = True

    # If no specific tags are present, the default is 1
    if non_colinear_calculation or spin_orbit_calculation:
        nspin = 4
    elif lsda:
        nspin = 2
    else:
        nspin = 1

    symmetries = []
    lattice_symmetries = []  # note: will only contain lattice symmetries that are NOT crystal symmetries
    inversion_symmetry = False

    # See also PW/src/setup.f90
    nsym = outputs.get('symmetries', {}).get('nsym', None)  # crystal symmetries
    nrot = outputs.get('symmetries', {}).get('nrot', None)  # lattice symmetries

    for symmetry in outputs.get('symmetries', {}).get('symmetry', []):

        # There are two types of symmetries, lattice and crystal. The pure inversion (-I) is always a lattice symmetry,
        # so we don't care. But if the pure inversion is also a crystal symmetry, then then the system as a whole
        # has (by definition) inversion symmetry; so we set the global property inversion_symmetry = True.
        symmetry_type = symmetry['info']['$']
        symmetry_name = symmetry['info']['@name']
        if symmetry_type == 'crystal_symmetry' and symmetry_name.lower() == 'inversion':
            inversion_symmetry = True

        sym = {
            'rotation': [
                symmetry['rotation']['$'][0:3],
                symmetry['rotation']['$'][3:6],
                symmetry['rotation']['$'][6:9],
            ],
            'name': symmetry_name,
        }

        try:
            sym['t_rev'] = '1' if symmetry['info']['@time_reversal'] else '0'
        except KeyError:
            sym['t_rev'] = '0'

        try:
            sym['equivalent_atoms'] = symmetry['equivalent_atoms']['$']
        except KeyError:
            pass

        try:
            sym['fractional_translation'] = symmetry['fractional_translation']
        except KeyError:
            pass

        if symmetry_type == 'crystal_symmetry':
            symmetries.append(sym)
        elif symmetry_type == 'lattice_symmetry':
            lattice_symmetries.append(sym)
        else:
            raise XMLParseError(f'Unexpected type of symmetry: {symmetry_type}')

    if (nsym != len(symmetries)) or (nrot != len(symmetries) + len(lattice_symmetries)):
        logs.warning.append(
            'Inconsistent number of symmetries: nsym={}, nrot={}, len(symmetries)={}, len(lattice_symmetries)={}'.
            format(nsym, nrot, len(symmetries), len(lattice_symmetries))
        )

    xml_data = {
        #'pp_check_flag': True, # Currently not printed in the new format.
        # Signals whether the XML file is complete
        # and can be used for post-processing. Everything should be in the XML now, but in
        # any case, the new XML schema should mostly protect from incomplete files.
        'lkpoint_dir': False,  # Currently not printed in the new format.
        # Signals whether kpt-data are written in sub-directories.
        # Was generally true in the old format, but now all the eigenvalues are
        # in the XML file, under output / band_structure, so this is False.
        'charge_density': './charge-density.dat',  # A file name. Not printed in the new format.
        # The filename and path are considered fixed: <outdir>/<prefix>.save/charge-density.dat
        # TODO: change to .hdf5 if output format is HDF5 (issue #222)
        'rho_cutoff_units': 'eV',
        'wfc_cutoff_units': 'eV',
        'fermi_energy_units': 'eV',
        'k_points_units': '1 / angstrom',
        'symmetries_units': 'crystal',
        'constraint_mag': constraint_mag,
        'magnetization_angle2': magnetization_angle2,
        'magnetization_angle1': magnetization_angle1,
        'starting_magnetization': starting_magnetization,
        'has_electric_field': has_electric_field,
        'has_dipole_correction': has_dipole_correction,
        'lda_plus_u_calculation': 'dftU' in outputs,
        'format_name': xml_dictionary['general_info']['xml_format']['@NAME'],
        'format_version': xml_dictionary['general_info']['xml_format']['@VERSION'],
        # TODO: check that format version: a) matches the XSD schema version; b) is updated as well
        #       See line 43 in Modules/qexsd.f90
        'creator_name': xml_dictionary['general_info']['creator']['@NAME'].lower(),
        'creator_version': xml_dictionary['general_info']['creator']['@VERSION'],
        'non_colinear_calculation': non_colinear_calculation,
        'do_magnetization': do_magnetization,
        'time_reversal_flag': time_reversal,
        'symmetries': symmetries,
        'lattice_symmetries': lattice_symmetries,
        'do_not_use_time_reversal': inputs.get('symmetry_flags', {}).get('noinv', None),
        'spin_orbit_domag': do_magnetization,
        'fft_grid': [value for _, value in sorted(outputs['basis_set']['fft_grid'].items())],
        'lsda': lsda,
        'number_of_spin_components': nspin,
        'no_time_rev_operations': inputs.get('symmetry_flags', {}).get('no_t_rev', None),
        'inversion_symmetry':
        inversion_symmetry,  # the old tag was INVERSION_SYMMETRY and was set to (from the code): "invsym    if true the system has inversion symmetry"
        'number_of_bravais_symmetries': nrot,  # lattice symmetries
        'number_of_symmetries': nsym,  # crystal symmetries
        'wfc_cutoff': inputs.get('basis', {}).get('ecutwfc', -1.0) * CONSTANTS.hartree_to_ev,
        'rho_cutoff': outputs['basis_set']['ecutrho'] * CONSTANTS.hartree_to_ev,  # not always printed in input->basis
        'smooth_fft_grid': [value for _, value in sorted(outputs['basis_set']['fft_smooth'].items())],
        'dft_exchange_correlation': inputs.get('dft', {}).get('functional',
                                                              None),  # TODO: also parse optional elements of 'dft' tag
        # WARNING: this is different between old XML and new XML
        'spin_orbit_calculation': spin_orbit_calculation,
        'q_real_space': outputs['algorithmic_info']['real_space_q'],
    }

    # alat is technically an optional attribute according to the schema,
    # but I don't know what to do if it's missing. atomic_structure is mandatory.
    output_alat_bohr = outputs['atomic_structure']['@alat']
    output_alat_angstrom = output_alat_bohr * CONSTANTS.bohr_to_ang

    # Band structure
    if 'band_structure' in outputs:
        band_structure = outputs['band_structure']

        smearing_xml = None

        if 'smearing' in outputs['band_structure']:
            smearing_xml = outputs['band_structure']['smearing']
        elif 'smearing' in inputs:
            smearing_xml = inputs['smearing']

        if smearing_xml:
            degauss = smearing_xml['@degauss']

            # Versions below 19.03.04 (Quantum ESPRESSO<=6.4.1) incorrectly print degauss in Ry instead of Hartree
            if xml_version < StrictVersion('19.03.04'):
                degauss *= CONSTANTS.ry_to_ev
            else:
                degauss *= CONSTANTS.hartree_to_ev

            xml_data['degauss'] = degauss
            xml_data['smearing_type'] = smearing_xml['$']

        num_k_points = band_structure['nks']
        num_electrons = band_structure['nelec']
        num_atomic_wfc = band_structure['num_of_atomic_wfc']
        num_bands = band_structure.get('nbnd', None)
        num_bands_up = band_structure.get('nbnd_up', None)
        num_bands_down = band_structure.get('nbnd_dw', None)

        if num_bands is None and num_bands_up is None and num_bands_down is None:
            raise XMLParseError('None of `nbnd`, `nbnd_up` or `nbdn_dw` could be parsed.')

        # If both channels are `None` we are dealing with a non spin-polarized or non-collinear calculation
        elif num_bands_up is None and num_bands_down is None:
            spins = False

        # If only one of the channels is `None` we raise, because that is an inconsistent result
        elif num_bands_up is None or num_bands_down is None:
            raise XMLParseError('Only one of `nbnd_up` and `nbnd_dw` could be parsed')

        # Here it is a spin-polarized calculation, where for pw.x the number of bands in each channel should be identical.
        else:
            spins = True
            if num_bands_up != num_bands_down:
                raise XMLParseError(f'different number of bands for spin channels: {num_bands_up} and {num_bands_down}')

            if num_bands is not None and num_bands != num_bands_up + num_bands_down:
                raise XMLParseError(
                    'Inconsistent number of bands: nbnd={}, nbnd_up={}, nbnd_down={}'.format(
                        num_bands, num_bands_up, num_bands_down
                    )
                )

            if num_bands is None:
                num_bands = num_bands_up + num_bands_down  # backwards compatibility;

        k_points = []
        k_points_weights = []
        ks_states = band_structure['ks_energies']
        for ks_state in ks_states:
            k_points.append([kp * 2 * np.pi / output_alat_angstrom for kp in ks_state['k_point']['$']])
            k_points_weights.append(ks_state['k_point']['@weight'])

        if not spins:
            band_eigenvalues = [[]]
            band_occupations = [[]]
            for ks_state in ks_states:
                band_eigenvalues[0].append(ks_state['eigenvalues']['$'])
                band_occupations[0].append(ks_state['occupations']['$'])
        else:
            band_eigenvalues = [[], []]
            band_occupations = [[], []]
            for ks_state in ks_states:
                band_eigenvalues[0].append(ks_state['eigenvalues']['$'][0:num_bands_up])
                band_eigenvalues[1].append(ks_state['eigenvalues']['$'][num_bands_up:num_bands])
                band_occupations[0].append(ks_state['occupations']['$'][0:num_bands_up])
                band_occupations[1].append(ks_state['occupations']['$'][num_bands_up:num_bands])

        band_eigenvalues = np.array(band_eigenvalues) * CONSTANTS.hartree_to_ev
        band_occupations = np.array(band_occupations)

        if not spins:
            parser_assert_equal(
                band_eigenvalues.shape, (1, num_k_points, num_bands), 'Unexpected shape of band_eigenvalues'
            )
            parser_assert_equal(
                band_occupations.shape, (1, num_k_points, num_bands), 'Unexpected shape of band_occupations'
            )
        else:
            parser_assert_equal(
                band_eigenvalues.shape, (2, num_k_points, num_bands_up), 'Unexpected shape of band_eigenvalues'
            )
            parser_assert_equal(
                band_occupations.shape, (2, num_k_points, num_bands_up), 'Unexpected shape of band_occupations'
            )

        if not spins:
            xml_data['number_of_bands'] = num_bands
        else:
            # For collinear spin-polarized calculations `spins=True` and `num_bands` is sum of both channels. To get the
            # actual number of bands, we divide by two using integer division
            xml_data['number_of_bands'] = num_bands // 2

        for key, value in [('number_of_bands_up', num_bands_up), ('number_of_bands_down', num_bands_down)]:
            if value is not None:
                xml_data[key] = value

        if 'fermi_energy' in band_structure:
            xml_data['fermi_energy'] = band_structure['fermi_energy'] * CONSTANTS.hartree_to_ev

        if 'two_fermi_energies' in band_structure:
            xml_data['fermi_energy_up'], xml_data['fermi_energy_down'] = [
                energy * CONSTANTS.hartree_to_ev for energy in band_structure['two_fermi_energies']
            ]

        bands_dict = {
            'occupations': band_occupations,
            'bands': band_eigenvalues,
            'bands_units': 'eV',
        }

        xml_data['number_of_atomic_wfc'] = num_atomic_wfc
        xml_data['number_of_k_points'] = num_k_points
        xml_data['number_of_electrons'] = num_electrons
        xml_data['k_points'] = k_points
        xml_data['k_points_weights'] = k_points_weights
        xml_data['bands'] = bands_dict

    try:
        monkhorst_pack = inputs['k_points_IBZ']['monkhorst_pack']
    except KeyError:
        pass  # not using Monkhorst pack
    else:
        xml_data['monkhorst_pack_grid'] = [monkhorst_pack[attr] for attr in ['@nk1', '@nk2', '@nk3']]
        xml_data['monkhorst_pack_offset'] = [monkhorst_pack[attr] for attr in ['@k1', '@k2', '@k3']]

    if occupations is not None:
        xml_data['occupations'] = occupations

    if 'boundary_conditions' in outputs and 'assume_isolated' in outputs['boundary_conditions']:
        xml_data['assume_isolated'] = outputs['boundary_conditions']['assume_isolated']

    # This is not printed by QE 6.3, but will be re-added before the next version
    if 'real_space_beta' in outputs['algorithmic_info']:
        xml_data['beta_real_space'] = outputs['algorithmic_info']['real_space_beta']

    conv_info = {}
    conv_info_scf = {}
    conv_info_opt = {}
    # NOTE: n_scf_steps refers to the number of SCF steps in the *last* loop only.
    # To get the total number of SCF steps in the run you should sum up the individual steps.
    # TODO: should we parse 'steps' too? Are they already added in the output trajectory?
    for key in ['convergence_achieved', 'n_scf_steps', 'scf_error']:
        try:
            conv_info_scf[key] = outputs['convergence_info']['scf_conv'][key]
        except KeyError:
            pass
    for key in ['convergence_achieved', 'n_opt_steps', 'grad_norm']:
        try:
            conv_info_opt[key] = outputs['convergence_info']['opt_conv'][key]
        except KeyError:
            pass
    if conv_info_scf:
        conv_info['scf_conv'] = conv_info_scf
    if conv_info_opt:
        conv_info['opt_conv'] = conv_info_opt
    if conv_info:
        xml_data['convergence_info'] = conv_info

    if 'status' in xml_dictionary:
        xml_data['exit_status'] = xml_dictionary['status']
        # 0 = convergence reached;
        # -1 = SCF convergence failed;
        # 3 = ionic convergence failed
        # These might be changed in the future. Also see PW/src/run_pwscf.f90

    try:
        berry_phase = outputs['electric_field']['BerryPhase']
    except KeyError:
        pass
    else:
        # This is what I would like to do, but it's not retro-compatible
        # xml_data['berry_phase'] = {}
        # xml_data['berry_phase']['total_phase']         = berry_phase['totalPhase']['$']
        # xml_data['berry_phase']['total_phase_modulus'] = berry_phase['totalPhase']['@modulus']
        # xml_data['berry_phase']['total_ionic_phase']      = berry_phase['totalPhase']['@ionic']
        # xml_data['berry_phase']['total_electronic_phase'] = berry_phase['totalPhase']['@electronic']
        # xml_data['berry_phase']['total_polarization']           = berry_phase['totalPolarization']['polarization']['$']
        # xml_data['berry_phase']['total_polarization_modulus']   = berry_phase['totalPolarization']['modulus']
        # xml_data['berry_phase']['total_polarization_units']     = berry_phase['totalPolarization']['polarization']['@Units']
        # xml_data['berry_phase']['total_polarization_direction'] = berry_phase['totalPolarization']['direction']
        # parser_assert_equal(xml_data['berry_phase']['total_phase_modulus'].lower(), '(mod 2)',
        #                    "Unexpected modulus for total phase")
        # parser_assert_equal(xml_data['berry_phase']['total_polarization_units'].lower(), 'e/bohr^2',
        #                    "Unsupported units for total polarization")
        # Retro-compatible keys:
        polarization = berry_phase['totalPolarization']['polarization']['$']
        polarization_units = berry_phase['totalPolarization']['polarization']['@Units']
        polarization_modulus = berry_phase['totalPolarization']['modulus']
        parser_assert(
            polarization_units in ['e/bohr^2', 'C/m^2'],
            f"Unsupported units '{polarization_units}' of total polarization"
        )
        if polarization_units == 'e/bohr^2':
            polarization *= e_bohr2_to_coulomb_m2
            polarization_modulus *= e_bohr2_to_coulomb_m2

        xml_data['total_phase'] = berry_phase['totalPhase']['$']
        xml_data['total_phase_units'] = '2pi'
        xml_data['ionic_phase'] = berry_phase['totalPhase']['@ionic']
        xml_data['ionic_phase_units'] = '2pi'
        xml_data['electronic_phase'] = berry_phase['totalPhase']['@electronic']
        xml_data['electronic_phase_units'] = '2pi'
        xml_data['polarization'] = polarization
        xml_data['polarization_module'] = polarization_modulus  # should be called "modulus"
        xml_data['polarization_units'] = 'C / m^2'
        xml_data['polarization_direction'] = berry_phase['totalPolarization']['direction']
        # TODO: add conversion for (e/Omega).bohr (requires to know Omega, the volume of the cell)
        # TODO (maybe): Not parsed:
        # - individual ionic phases
        # - individual electronic phases and weights

    # TODO: We should put the `non_periodic_cell_correction` string in (?)
    atoms = [[atom['@name'], [coord * CONSTANTS.bohr_to_ang
                              for coord in atom['$']]]
             for atom in outputs['atomic_structure']['atomic_positions']['atom']]
    species = outputs['atomic_species']['species']
    structure_data = {
        'atomic_positions_units':
        'Angstrom',
        'direct_lattice_vectors_units':
        'Angstrom',
        # ??? 'atoms_if_pos_list': [[1, 1, 1], [1, 1, 1]],
        'number_of_atoms':
        outputs['atomic_structure']['@nat'],
        'lattice_parameter':
        output_alat_angstrom,
        'reciprocal_lattice_vectors': [
            outputs['basis_set']['reciprocal_lattice']['b1'], outputs['basis_set']['reciprocal_lattice']['b2'],
            outputs['basis_set']['reciprocal_lattice']['b3']
        ],
        'atoms':
        atoms,
        'cell': {
            'lattice_vectors': lattice_vectors,
            'volume': cell_volume(*lattice_vectors),
            'atoms': atoms,
        },
        'lattice_parameter_xml':
        output_alat_bohr,
        'number_of_species':
        outputs['atomic_species']['@ntyp'],
        'species': {
            'index': [i + 1 for i, specie in enumerate(species)],
            'pseudo': [specie['pseudo_file'] for specie in species],
            'mass': [specie['mass'] for specie in species],
            'type': [specie['@name'] for specie in species]
        },
    }

    xml_data['structure'] = structure_data

    return xml_data, logs

max_stars_repo_name: lin-cp/aiida-quantumespresso

max_issues_repo_name: lin-cp/aiida-quantumespresso

max_forks_repo_name: lin-cp/aiida-quantumespresso


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

# To set up pyvista use Miniconda and pip e.g.
#
#  wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh
#  chmod +x Miniconda3-latest-MacOSX-x86_64.sh
#  ./Miniconda3-latest-MacOSX-x86_64.sh -b -p miniconda3
# source miniconda3/bin/activate
# conda create --name vtk_env python=3.8
# conda activate vtk_env
# conda install nodejs matplotlib
# pip install jupyter pyvista ipyvtk_simple
#
# can be be run as a script e.g.
#
# python plot-field.py
#
# or in notebook -
#
# jupyter notebook

import vtk
import pyvista as pv
from matplotlib import cm
import numpy as np

fcode='θ'
# fcode='state_debug'

mesh22 = pv.read('Coupler_UnitTest_atmosphere_long_num0004.pvtu')

xc=mesh22.get_array('xc')
yc=mesh22.get_array('yc')
zc=mesh22.get_array('zc')

R=np.sqrt(np.square(xc)+np.square(yc)+np.square(zc))

# Mask to surface and "top" half of sphere.
tol=1e-1
sval=( np.abs(R[0]-R)>tol ) | ( zc>0 )

l=np.array(list(range(len( R )))); l[ np.abs(R[0]-R)<tol ];len(l);len( R[np.abs(R[0]-R)<tol] )
ptnan=l[sval]
ptplot=l[~sval]

mesh22[fcode][ptnan]=np.nan

msurf=mesh22.extract_points(ptplot,adjacent_cells=True,include_cells=True)

plotter = pv.Plotter()
# plotter.add_mesh(msurf,style='points',show_edges=True,scalars='θ',culling=True,cmap='jet',lighting=True,nan_opacity=0,point_size=10)
plotter.add_mesh(msurf,style='surface',show_edges=False,scalars=fcode,culling=True,cmap='jet',lighting=True,nan_opacity=0)
# plotter.add_mesh(msurf,style='points',scalars='θ',culling=True,cmap='jet',lighting=True,nan_opacity=0)
plotter.show()

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
import pyqpanda.pyQPanda as pq
from pyqpanda.Visualization.circuit_draw import *
import numpy as np
import math
import networkx as nx
import matplotlib.pyplot as plt

class InitQMachine:
    def __init__(self, machineType = pq.QMachineType.CPU):
        self.m_machine = pq.init_quantum_machine(machineType)

    def __del__(self):
        pq.destroy_quantum_machine(self.m_machine)

def plot(list):
    import matplotlib.pyplot as plt
    import matplotlib as mpl

    mpl.rcParams["font.sans-serif"] = ["SimHei"]
    mpl.rcParams["axes.unicode_minus"] = False
    x = [key for key in list.keys()]
    y = [val for val in list.values()]
    y1 = [val/sum(y) for val in list.values()]
    plt.bar(x, y1, align = "center", color = "b", alpha = 0.6)
    plt.ylabel("Probabilities")
    plt.grid(True, axis = "y", ls = ":", color = "r", alpha = 0.3)
    plt.show()

def test_layer1():
    init_machine = InitQMachine()
    machine = init_machine.m_machine
    q = machine.qAlloc_many(6)
    c = machine.cAlloc_many(6)
    prog = pq.QProg()
    cir = pq.QCircuit()
    cir.insert(pq.T(q[0])).insert(pq.SWAP(q[1], q[5])).insert(pq.S(q[1])).insert(pq.CNOT(q[1], q[0]))
    cir.insert(pq.CZ(q[0], q[2])).insert(pq.CU(np.pi/3, 3, 4, 5, q[5], q[2]))
    cir.set_control([q[4],q[3]])
    cir.set_dagger(True)

    prog.insert(pq.T(q[0])).insert(pq.SWAP(q[1], q[2])).insert(pq.X1(q[3]))
    prog.insert(pq.X(q[0])).insert(pq.CZ(q[1], q[2])).insert(pq.Y(q[3])).insert(pq.Z(q[3]))
    prog.insert(pq.S(q[3])).insert(pq.CNOT(q[1], q[4]))
    prog.insert(pq.CR(q[0], q[2], 2.3)).insert(pq.CR(q[0], q[1], np.pi/2.3)).insert(pq.Y1(q[4])).insert(pq.Z1(q[5]))
    prog.insert(pq.CR(q[5], q[2], np.pi/3))
    prog.insert(pq.CR(q[2], q[3], np.pi/5))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.iSWAP(q[4], q[1]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(pq.U1(q[0], np.pi/4)).insert(pq.U2(q[1], 1, np.pi/2)).insert(pq.U3(q[2], 1.2, 2, np.pi/3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(pq.U1(q[0], np.pi/4)).insert(pq.U2(q[1], 1, np.pi/2)).insert(pq.U3(q[2], 1.2, 2, np.pi/3))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[5], q[2]))
    prog.insert(pq.T(q[0])).insert(pq.SWAP(q[1], q[2])).insert(pq.X1(q[3]))
    prog.insert(pq.X(q[0])).insert(pq.CZ(q[1], q[2])).insert(pq.Y(q[3])).insert(pq.Z(q[3]))
    prog.insert(pq.S(q[3]))
    prog.insert(cir)
    prog.insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[1]))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[3], q[2]))
    prog.insert(pq.CR(q[0], q[2], 2.3)).insert(pq.CR(q[0], q[1], np.pi/2.3)).insert(pq.Y1(q[4])).insert(pq.Z1(q[5]))
    prog.insert(pq.CR(q[5], q[2], np.pi/3))
    prog.insert(pq.CR(q[2], q[3], np.pi/5))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.iSWAP(q[4], q[1]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(pq.CR(q[0], q[2], 2.3)).insert(pq.CR(q[0], q[1], np.pi/2.3)).insert(pq.Y1(q[4])).insert(pq.Z1(q[5]))
    prog.insert(pq.CR(q[5], q[2], np.pi/3))
    prog.insert(pq.CR(q[2], q[3], np.pi/5))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.iSWAP(q[4], q[1]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(pq.U1(q[0], np.pi/4)).insert(pq.U2(q[1], 1, np.pi/2)).insert(pq.U3(q[2], 1.2, 2, np.pi/3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(pq.U1(q[0], np.pi/4)).insert(pq.U2(q[1], 1, np.pi/2)).insert(pq.U3(q[2], 1.2, 2, np.pi/3))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[5], q[2]))
    prog.insert(pq.T(q[0])).insert(pq.SWAP(q[1], q[2])).insert(pq.X1(q[3]))
    prog.insert(pq.X(q[0])).insert(pq.CZ(q[1], q[2])).insert(pq.Y(q[3])).insert(pq.Z(q[3]))
    prog.insert(pq.S(q[3]))
    prog.insert(cir)
    prog.insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[1]))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[3], q[2]))
    prog.insert(pq.S(q[3]))
    prog.insert(cir)
    prog.insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[1]))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[3], q[2]))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[5], q[2]))
    prog.insert(pq.T(q[0])).insert(pq.SWAP(q[1], q[2])).insert(pq.X1(q[3]))
    prog.insert(pq.X(q[0])).insert(pq.CZ(q[1], q[2])).insert(pq.Y(q[3])).insert(pq.Z(q[3]))
    prog.insert(pq.S(q[3]))
    prog.insert(cir)
    prog.insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[1]))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[3], q[2]))
    prog.insert(pq.S(q[3]))
    prog.insert(cir)
    prog.insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[5]))
    prog.insert(pq.CR(q[0], q[2], 2.3)).insert(pq.CR(q[0], q[1], np.pi/2.3)).insert(pq.Y1(q[4])).insert(pq.Z1(q[5]))
    prog.insert(pq.CR(q[5], q[2], np.pi/3))
    prog.insert(pq.CR(q[2], q[3], np.pi/5))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.iSWAP(q[4], q[1]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(pq.U1(q[0], np.pi/4)).insert(pq.U2(q[1], 1, np.pi/2)).insert(pq.U3(q[2], 1.2, 2, np.pi/3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(cir)
    prog.insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[5]))
    prog.insert(pq.CR(q[0], q[2], 2.3)).insert(pq.CR(q[0], q[1], np.pi/2.3)).insert(pq.Y1(q[4])).insert(pq.Z1(q[5]))
    prog.insert(pq.CR(q[5], q[2], np.pi/3))
    prog.insert(pq.CR(q[2], q[3], np.pi/5))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.iSWAP(q[4], q[1]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[2], np.pi/3))
    prog.insert(pq.U1(q[0], np.pi/4)).insert(pq.U2(q[5], 1, np.pi/2)).insert(pq.U3(q[5], 1.2, 2, np.pi/3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.U4(1.5, 2.2, np.pi/3, np.pi/4, q[3]))
    prog.insert(pq.RX(q[2], np.pi/6.2)).insert(pq.RY(q[2], 6.6)).insert(pq.RZ(q[5], np.pi/3))
    prog.insert(pq.U1(q[0], np.pi/4)).insert(pq.U2(q[1], 1, np.pi/2)).insert(pq.U3(q[2], 1.2, 2, np.pi/3))
    prog.insert(pq.CU(np.pi/3, 3, 4, 5, q[5], q[2]))
    prog.insert(pq.H(q[0])).insert(pq.iSWAP(q[5], q[2], np.pi/2.3)).insert(pq.iSWAP(q[4], q[1]))
    prog.insert(pq.H(q[0])).insert(pq.H(q[1])).insert(pq.H(q[2])).insert(pq.H(q[3])).insert(pq.H(q[4])).insert(pq.H(q[5])).insert(pq.measure_all(q, c))
    prog.insert(pq.H(q[0])).insert(pq.H(q[1])).insert(pq.Measure(q[0], c[0]))

    print(prog)
    draw_circuit_pic(prog, 'D:/test_cir_my_draw.jpg', verbose=True)

if __name__=="__main__":
    test_layer1()
    print("Test over.")

reason: import numpy,import networkx