text:
#from math import *
import numpy as np
import math
try:
    import psyco
    psyco.full()
except ImportError:
    pass

R_earth = 6371000.785

sin = np.sin
cos = np.cos
tan = np.tan
atan = np.arctan
acos = np.arccos
asin = np.arcsin
log = np.log
exp = np.exp
sqrt = np.sqrt
pi = math.pi
radians = np.radians
degrees = np.degrees

def UTMReverse(E, N, lmb0, FN = 0.0 ): #328.):
    # parameters
    M0 = 0. # natural origin north
    FE = 500000. # false eastening
    # constants
    k0 = 0.9996
    # ellipsoid
    a = 6377563.396
    #oneDivF = 299.32496
    eSqr = 0.00667054

    eSSqr = eSqr/(1 - eSqr)
    e1 = (1 - (1 - eSqr)**0.5)/(1 + (1 + eSqr)**0.5)
    M1 = M0 + (N - FN)/k0
    mu1 = M1/(a*(1 - eSqr/4 - 3*(eSqr**2)/64 - 5*(eSqr**3)/256))
    phi1 = mu1 + (3 *(e1**1)/2  - 27*(e1**3)/32) * sin(2*mu1) + \
                 (21*(e1**2)/16 - 55*(e1**4)/32) * sin(4*mu1) + \
                 (151*(e1**3)/96               ) * sin(6*mu1) + \
                 (1097*(e1**4)/512             ) * sin(8*mu1)
    T1 = tan(phi1)**2
    C1 = eSSqr * (cos(phi1)**2)
    nu1 = a / ((1 - eSqr*(sin(phi1)**2))**0.5)
    ro1 = a*(1-eSqr)/((1 - eSqr*(sin(phi1)**2))**1.5)
    D = (E - FE)/(nu1 * k0)

    phi = phi1 - (nu1 * tan(phi1) / ro1) * ((D**2)/2 - \
                  (5 + 3*T1 + 10*C1 - 4*(C1**2) - 9*eSSqr) * (D**4) / 24 + \
                  (61 + 90*T1 + 298*C1 + 45*(T1**2) - 252 * eSSqr - 3 * (C1**2)) * (D**6)/720)
    lmb = lmb0 + (D - (1 + 2*T1 + C1)*(D**3)/6 + \
                  (5 - 2*C1 + 28*T1 - 3*(C1**2) + 8*eSSqr + 24*(T1**2))*(D**5)/120)/cos(phi1)
    return (lmb, phi)

def UTMForward(lmb, phi, lmb0):
    # parameters
    M0 = 0. # natural origin north
    FN = 0. # false northing
    FE = 500000. # false eastening
    # constants
    k0 = 0.9996
    # ellipsoid
    a = 6377563.396
    #oneDivF = 299.32496
    eSqr = 0.00667054

    eSSqr = eSqr/(1 - eSqr)
    T = tan(phi)**2
    C = eSqr * (cos(phi)**2) / (1 - eSqr)
    A = (lmb - lmb0) * cos(phi)
    nu = a / ((1 - eSqr*(sin(phi)**2))**0.5)
    M = a * ((1 - eSqr/4 - 3 * (eSqr**2)/64 - 5 * (eSqr**3)/256) * phi - \
             (3*eSqr/8 + 3*(eSqr**2)/32 + 45*(eSqr**3)/1024) * sin(2*phi) + \
             (15*(eSqr**2)/256 + 45*(eSqr**3)/1024) * sin(4*phi) + \
             (35*(eSqr**3)/3072) * sin(6*phi))
    E = FE + k0*nu*(A + (1 - T + C)*(A**3)/6 + \
                    (5 - 18*T + T**2 + 72*C - 58*eSSqr)*(A**5)/120)
    N = FN + k0*(M - M0 + nu * tan(phi)*((A**2)/2 + \
                    (5 - T + 9*C + 4*(C**2))*(A**4)/24 + \
                    (61 - 58*T + T**2 + 600*C - 330*eSSqr)*(A**6)/720))
    return (E,N)

def TopocentricForward(lmb, phi, lmb0, phi0):
    # ellipsoid
    a = 6377563.396
    #oneDivF = 299.32496
    eSqr = 0.00667054

    h = 0.0 # height above ground, approx.
    # convert to radians
    lmb = radians(lmb)
    phi = radians(phi)
    lmb0 = radians(lmb0)
    phi0 = radians(phi0)
    # precompute sin/cos
    sinlmbMlmb0 = sin(lmb-lmb0)
    coslmbMlmb0 = cos(lmb-lmb0)
    sinphi = sin(phi)
    cosphi = cos(phi)
    sinphi0 = sin(phi0)
    cosphi0 = cos(phi0)

    nu  = a / ((1. - eSqr*(sinphi **2))**0.5)
    nu0 = a / ((1. - eSqr*(sinphi0**2))**0.5)

    U = (nu + h)*cosphi * sinlmbMlmb0
    V = (nu + h)*(sinphi*cosphi0 - cosphi*sinphi0*coslmbMlmb0) + eSqr*(nu0*sinphi0 - nu*sinphi)*cosphi0
    return (U,V)

def TopocentricBackward(U,V,lmb0,phi0):
    # ellipsoid
    a = 6377563.396
    oneDivF = 299.32496
    eSqr = 0.00667054

    W = 0.0
    h0 = 0.0 #?
    # convert to radians
    lmb0 = radians(lmb0)
    phi0 = radians(phi0)
    # precompute sin/cos
    sinphi0 = sin(phi0)
    cosphi0 = cos(phi0)
    sinlmb0 = sin(lmb0)
    coslmb0 = cos(lmb0)

    nu0 = a / ((1. - eSqr*(sinphi0**2))**0.5)

    X0 = (nu0 + h0)*cosphi0*coslmb0
    Y0 = (nu0 + h0)*cosphi0*sinlmb0
    Z0 = ((1-eSqr)*nu0 + h0)*sinphi0

    X = X0 - U*sinlmb0 - V*sinphi0*coslmb0 + W*cosphi0*coslmb0
    Y = Y0 + U*coslmb0 - V*sinphi0*sinlmb0 + W*cosphi0*sinlmb0
    Z = Z0 + V*cosphi0 + W*sinphi0

    eps = eSqr/(1-eSqr)
    b = a*(1 - oneDivF)
    p = (X**2 + Y**2)**0.5
    q = atan((Z*a)/(p*b))
    phi = atan((Z+eps*b*(sin(q)**3))/(p - eps*a*(cos(q)**3)))
    lmb = atan(Y/X)
    return degrees(phi), degrees(lmb)

def UTMZoneForward(lmb, phi, zone):
    lmb0 = radians(zone * 6 - (180 + 6/2))
    return UTMForward(lmb, phi, lmb0)

def UTMZoneReverse(E, N, zone, FN = 0.0):
    lmb0 = radians(zone * 6 - (180 + 6/2))
    return UTMReverse(E, N, lmb0, FN)

def SwissForward(lmb, phi):
    #a = 6377397.155
    EE = 0.006674372230614
    E = sqrt(EE)
    #phi0 = radians(46.0 + 57./60. + (8.66)/3600.)
    lambda0 = radians(7.0 + 26./60. + (22.50)/3600.)
    R =  6378815.90365
    alpha =  1.00072913843038
    b0 = radians(46.0 + 54./60. + (27.83324844/3600.))
    K = 0.0030667323772751
    S = -alpha * log(tan(pi/4. - phi/2)) - alpha*E/2 * log( (1+E*sin(phi))/(1-E*sin(phi)) ) + K
    b = 2*(atan(exp(S)) - pi/4)
    l = alpha * (lmb - lambda0)
    l_ = atan( sin(l) / (sin(b0)*tan(b) + cos(b0)*cos(l)) )
    b_ = asin( cos(b0)*sin(b) - sin(b0)*cos(b)*cos(l) )
    Y = R*l_
    X = R/2 * log( (1+sin(b_))/(1-sin(b_)) )
    return X+200000, Y+600000

# no idea what went wrong in the swiss transformations ...
def SwissReverse(yy,  xx):
    #a = 6377397.155
    #EE = 0.006674372230614
    #E = sqrt(EE)
    #phi0 = radians(46.0 + 57./60. + (8.66)/3600.)
    #lambda0 = radians(7.0 + 26./60. + (22.50)/3600.)
    #R =  6378815.90365
    #alpha =  1.00072913843038
    #b0 = radians(46.0 + 54./60. + (27.83324844/3600.))
    #K = 0.0030667323772751
    Y = (yy-600000.)/1000000.
    X = (xx-200000.)/1000000.

    a1 = (+ 4.72973056
          + 0.7925714  * X
          + 0.132812 * (X**2)
          + 0.02550 * (X**3)
          + 0.0048 * (X**4))
    a3 = (- 0.044270
          - 0.02550 * X
          - 0.0096 * (X**2))
    a5 = + 0.00096
    lmb = 2.67825 + a1*Y + a3*(Y**3) + a5*(Y**5)

    p0 = (  0
          + 3.23864877 * X
          - 0.0025486 * (X**2)
          - 0.013245 * (X**3)
          + 0.000048 * (X**4))
    p2 = (- 0.27135379
          - 0.0450442 * X
          - 0.007553 * (X**2)
          - 0.00146 * (X**3))
    p4 = 0.002442 + 0.00132 * X
    phi = 16.902866 + p0 + p2*(Y**2) + p4*(Y**4)

    phi = phi * 10000./3600. * pi/180
    lmb = lmb * 10000./3600. * pi/180

    return (lmb,phi)

def distInMeter(E1, N1, E2, N2):
    E1 = radians(E1)
    N1 = radians(N1)
    E2 = radians(E2)
    N2 = radians(N2)
    acosarg = sin(N1)*sin(N2) + cos(N1)*cos(N2)*cos(E1-E2)
    acosarg = min(1,acosarg)
    acosarg = max(-1, acosarg)
    return acos(acosarg) * R_earth

if __name__ == "__main__":
    from math import *
    print ("Rigi:", map(degrees, SwissReverse(212273,679520)))

max_stars_repo_name: nooneisperfect/ReadYourMAPFile

max_issues_repo_name: nooneisperfect/ReadYourMAPFile

max_forks_repo_name: nooneisperfect/ReadYourMAPFile


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Calculate the pitch needed to maintain contact between the front wheel and
ground.
"""
from sympy import simplify, symbols
from sympy import Symbol
from sympy.physics.mechanics import ReferenceFrame, Point
from sympy.physics.mechanics import msprint
from sympy.mpmath import findroot
from sympy.utilities import lambdify
from sympy.printing import ccode

## define coordinates
# phi: roll
# theta: pitch
# delta: steer
phi, theta, delta = symbols('φ θ δ')
# rR: rear radius
# rF: front radius
rR, rF = symbols('rR rF')
# cR: distance from rear wheel center to steer axis
# cF: distance from front wheel center to steer axis
# ls: steer axis separation
cR, cF, ls = symbols('cR cF ls')

benchmark_parameters = {
    rR: 0.3,
    rF: 0.35,
    cR: 0.9534570696121847,
    ls: 0.2676445084476887,
    cF: 0.0320714267276193,
}

## define reference frames
# N: inertial frame
# B: rear aseembly frame
# H: front assembly frame
N = ReferenceFrame('N')
B = N.orientnew('B', 'body', [0, phi, theta], 'zxy') # yaw is ignored
H = B.orientnew('H', 'axis', [delta, B.z])

## define points
# rear wheel/ground contact point
pP = Point('P')

# define unit vectors from rear/front wheel centers to ground
# along the wheel plane
R_z = ((B.y ^ N.z) ^ B.y).normalize()
F_z = ((H.y ^ N.z) ^ H.y).normalize()

# define rear wheel center point
pRs = pP.locatenew('R*', -rR*R_z)

# "top" of steer axis, point of SA closest to R*
# orthogonal projection of rear wheel center on steer axis
pRh = pRs.locatenew('R^', cR*B.x)

# orthogonal projection of front wheel center on steer axis
pFh = pRh.locatenew('S^', ls*B.z)

# front wheel center point
pFs = pFh.locatenew('S*', cF*H.x)

# front wheel/ground contact point
pQ = pFs.locatenew('Q', rF*F_z)

# N.z component of vector to pQ from pP
# this is our configuration constraint
f = simplify(pQ.pos_from(pP) & N.z)

print("f = {}\n".format(msprint(f)))

# calculate the derivative of f for use with newton-raphson
df = simplify(f.diff(theta))
print("df/dθ = {}\n".format(msprint(df)))

# constraint function for zero steer/lean configuration and
# using the benchmark parameters
f0 = lambdify(theta, f.subs({phi: 0, delta: 0}).subs(benchmark_parameters))
df0 = lambdify(theta, df.subs({phi: 0, delta: 0}).subs(benchmark_parameters))

print("verifying constraint equations are correct")
print("for zero steer/lean, pitch should be pi/10")
findroot(f0, 0.3, solver="newton", tol=1e-8, verbose=True, df=df0)

c_sym = [Symbol('lean'), Symbol('pitch'), Symbol('steer')]
c_sym_dict = dict(zip([phi, theta, delta], c_sym))

fc = ccode(f.subs(c_sym_dict))
dfc = ccode(df.subs(c_sym_dict))

cs_math = {
    'cos': 'Math.Cos',
    'sin': 'Math.Sin',
    'pow': 'Math.Pow',
    'sqrt': 'Math.Sqrt'
}

fcs = fc
dfcs = dfc

for k, v in cs_math.items():
    fcs = fcs.replace(k, v)
    dfcs = dfcs.replace(k, v)

print(fcs)
print(dfcs)

ext: py

max_stars_repo_name: mickvangelderen/bikesim

max_issues_repo_name: mickvangelderen/bikesim

max_forks_repo_name: mickvangelderen/bikesim

reason: from sympy


text:
import os
import sys
import csv

import numpy as np
import cv2
from sklearn.externals import joblib

class ImageProcessor(object):
    def __init__(self, image_folder, classes_csv=None, classes_lst=None,
                 save_encoding_fname=None, load_encoding_fname=None):
        """
        Initialise image processor object, internal module used for preprocessing scraped/downloaded images
        """
        self.original_folder = os.getcwd()
        os.chdir(image_folder)
        self.image_folder = os.getcwd()
        self.images_filenames = [f for f in os.listdir(self.image_folder)]
        self.unprocessed_images = []
        self.max_img_height = 0
        self.max_img_width = 0
        self.val2idx = None
        self.idx2val = None
        self.classes_filenames_dict = {}

        if not (classes_csv or classes_lst):
            sys.exit('Provided neither csv filename or Python list for image classes')

        for image_filename in self.images_filenames:
            img = cv2.imread(image_filename)
            if img is None:
                self.unprocessed_images.append(image_filename)
            else:
                image_shape = img.shape
                if image_shape[0] > self.max_img_height:
                    self.max_img_height = image_shape[0]
                if image_shape[1] > self.max_img_width:
                    self.max_img_width = image_shape[1]

        self.processable_images = [filename for filename in self.images_filenames if
                                   filename not in self.unprocessed_images]

        if load_encoding_fname:
            encoding_dicts = joblib.load(load_encoding_fname)
            self.val2idx = encoding_dicts['val2idx']
            self.val2idx = encoding_dicts['idx2val']
        else:
            if classes_csv:
                with open(classes_csv) as csvfile:
                    data = csv.reader(csvfile)
                    # skip header
                    next(data, None)
                    self.classes = [row[0] for row in data]
            elif classes_lst:
                self.classes = classes_lst
            self.val2idx, self.idx2val = create_encode_dicts(self.classes)
            if save_encoding_fname:
                joblib.dump({'val2idx': self.val2idx, 'idx2val': self.idx2val}, save_encoding_fname)
            for class_name, index in self.val2idx.items():
                for image_filename in self.processable_images:
                    if class_name in image_filename:
                        if index not in self.classes_filenames_dict:
                            self.classes_filenames_dict[index] = []
                        self.classes_filenames_dict[index].append(image_filename)

    def load_images(self, pad=True, stretch=True, normalise='dataset'):
        if not self.classes_filenames_dict:
            sys.exit('Classes have not been loaded into Image Processor')
        else:
            loaded_images = []
            output_vectors = []
            for class_index, class_files in self.classes_filenames_dict.items():
                print('loading images for:', self.idx2val[class_index])
                for class_file in class_files:
                    output_vec = [0] * len(self.idx2val)
                    output_vec[class_index] = 1
                    file_img = cv2.imread(class_file)
                    final_image = resize_image(file_img, pad=pad, stretch=stretch)
                    if normalise == 'image':
                        # For per image per channel normalisation
                        normalised_image = normalise_image(final_image)
                        loaded_images.append(np.array(normalised_image))
                    loaded_images.append(np.array(final_image))
                    output_vectors.append(np.array(output_vec))
            if normalise == 'dataset':
                loaded_images = normalise_dataset(np.array(loaded_images))
            else:
                loaded_images = np.array(loaded_images)
            output_vectors = np.array(output_vectors)
            return loaded_images, output_vectors

def normalise_image(img, channels=3):
    output_image = None
    for ch_num in range(channels):
        ch_image = img[:, :, ch_num]
        normalised_image = (ch_image - ch_image.mean()) / ch_image.std()
        if ch_num == 0:
            output_image = normalised_image
        else:
            output_image = np.dstack((output_image, normalised_image))
    return output_image

def normalise_dataset(dataset, channels=3):
    ch_dataset = None
    output_dataset = None
    for ch_num in range(channels):
        if channels > 1:
            ch_dataset = np.expand_dims(dataset[:, :, :, ch_num], axis=3)
        elif channels == 1:
            ch_dataset = dataset
        normalised_dataset = (ch_dataset - ch_dataset.mean()) / ch_dataset.std()
        if ch_num == 0:
            output_dataset = normalised_dataset
        else:
            output_dataset = np.concatenate((output_dataset, normalised_dataset), axis=3)
    return output_dataset

def resize_image(img, pad=True, stretch=True, resized_width=200, resized_height=200):
    resized_aspect_ratio = float(resized_width)/resized_height
    (height, width, channel) = img.shape
    if stretch or float(width)/height == resized_aspect_ratio:
        output_image = cv2.resize(img, (resized_width, resized_height))
    else:
        if height > width:
            ratio = height/resized_height
            output_image = cv2.resize(img, (int(width / ratio), resized_height))
            if pad:
                padding_pixel = [0, 0, 0]
                padding_size = resized_width - int(width/ratio)
                left, right = __padding_size_split__(padding_size)
                output_image = cv2.copyMakeBorder(output_image, 0, 0, left, right,
                                                  cv2.BORDER_CONSTANT, value=padding_pixel)
        else:
            ratio = width/resized_width
            output_image = cv2.resize(img, (resized_width, int(height / ratio)))
            if pad:
                padding_pixel = [0, 0, 0]
                padding_size = resized_height - int(height / ratio)
                top, bottom = __padding_size_split__(padding_size)
                output_image = cv2.copyMakeBorder(output_image, top, bottom, 0, 0,
                                                  cv2.BORDER_CONSTANT, value=padding_pixel)
    return output_image

def __padding_size_split__(padding_size):
    if padding_size % 2:
        return int((padding_size/2) - 0.5), int((padding_size/2) + 0.5)
    else:
        return int(padding_size/2), int(padding_size/2)

def create_encode_dicts(iterable):
    val2idx = None
    if isinstance(iterable, list):
        val2idx = {value: index for index, value in enumerate(iterable)}
    elif isinstance(iterable, dict):
        val2idx = {value: index for index, value in enumerate(iterable.keys())}
    idx2val = {index: value for value, index in val2idx.items()}
    return val2idx, idx2val

def encode_val(value, val2idx_dict):
    return val2idx_dict[value]

def decode_idx(index, idx2val_dict):
    return idx2val_dict[index]

def encode_dict_keys(dict2encode, val2idx_dict):
    for old_key in dict2encode.keys():
        dict2encode[val2idx_dict[old_key]] = dict2encode.pop(old_key)
    return dict2encode

if __name__ == '__main__':
    processor = ImageProcessor('images', '../image_scrapper/doggotime_breeds - original.csv')
    processor.load_images(pad=True)

max_stars_repo_name: limjiaxiang/doggo-cv-chatbot

max_issues_repo_name: limjiaxiang/doggo-cv-chatbot

max_forks_repo_name: limjiaxiang/doggo-cv-chatbot

reason: import numpy


text:
def app():

    import os
    import pickle
    import numpy as np
    import pandas as pd
    import streamlit as st
    import lang.diabetes as lng

    st.title(lng.diabetes_diagnosis)

    with st.form(key='input_form'):

        fields = (
            st.slider(lng.age, 1, 100, 25, 1),
            st.radio(lng.gender, (lng.male, lng.female)),
            st.checkbox(lng.polyuria),
            st.checkbox(lng.polydipsia),
            st.checkbox(lng.sudden_weight_loss),
            st.checkbox(lng.weakness),
            st.checkbox(lng.polyphagia),
            st.checkbox(lng.genital_thrush),
            st.checkbox(lng.visual_blurring),
            st.checkbox(lng.itching),
            st.checkbox(lng.irritability),
            st.checkbox(lng.delayed_healing),
            st.checkbox(lng.partial_paresis),
            st.checkbox(lng.muscle_stiffness),
            st.checkbox(lng.alopecia),
            st.checkbox(lng.obesity))

        submit = st.form_submit_button(label=lng.check)

    if submit:
        data = pd.DataFrame([fields])
        data.replace({lng.male: 1, lng.female: 0}, inplace=True)
        data.replace({True: 1, False: 0}, inplace=True)
        app_path = os.path.dirname(__file__)
        model_path = os.path.join(app_path, 'models', 'diabetes.pkl')
        model = pickle.load(open(model_path,'rb'))
        probas = model.predict_proba(data)[:,1]
        predicted = model.predict(data)
        e = '⚠' if predicted[0] else '✅'
        p = lng.high_prob if predicted[0] else lng.low_prob
        st.markdown(f'{p} ({probas[0]:.4f}) {e}')

max_stars_repo_name: vaaliferov/paranormal

max_issues_repo_name: vaaliferov/paranormal

max_forks_repo_name: vaaliferov/paranormal


text:
from agent.TradingAgent import TradingAgent
import pandas as pd
import numpy as np
import os
from contributed_traders.util import get_file

class valdni_agentLinearReg(TradingAgent):
    """
    Simple Trading Agent that compares the past mid-price observations and places a
    buy limit order if the first window mid-price exponential average >= the second window mid-price exponential average or a
    sell limit order if the first window mid-price exponential average < the second window mid-price exponential average
    """

    def __init__(self, id, name, type, symbol, starting_cash,
                 min_size, max_size, wake_up_freq='60s',
                 log_orders=False, random_state=None):

        super().__init__(id, name, type, starting_cash=starting_cash, log_orders=log_orders, random_state=random_state)
        self.symbol = symbol
        self.min_size = min_size  # Minimum order size
        self.max_size = max_size  # Maximum order size
        self.size = self.random_state.randint(self.min_size, self.max_size)
        self.wake_up_freq = wake_up_freq
        self.mid_list, self.avg_win1_list, self.avg_win2_list = [], [], []
        self.log_orders = log_orders
        self.lookback = 50
        self.trades = []
        self.state = "AWAITING_WAKEUP"

    def kernelStarting(self, startTime):
        super().kernelStarting(startTime)

    def wakeup(self, currentTime):
        """ Agent wakeup is determined by self.wake_up_freq """
        can_trade = super().wakeup(currentTime)
        if not can_trade: return
        self.getLastTrade(self.symbol)
        self.state = "AWAITING_LAST_TRADE"

    def dump_shares(self):
        # get rid of any outstanding shares we have
        if self.symbol in self.holdings and len(self.orders) == 0:
            order_size = self.holdings[self.symbol]
            bid, _, ask, _ = self.getKnownBidAsk(self.symbol)
            if bid:
                self.placeLimitOrder(self.symbol, quantity=order_size, is_buy_order=False, limit_price=0)

    def receiveMessage(self, currentTime, msg):
        """ Momentum agent actions are determined after obtaining the best bid and ask in the LOB """
        super().receiveMessage(currentTime, msg)
        if self.state == "AWAITING_LAST_TRADE" and msg.body['msg'] == "QUERY_LAST_TRADE":
            dt = (self.mkt_close - currentTime) / np.timedelta64(1, 'm')
            if dt < 10:
                self.dump_shares()
            last = self.last_trade[self.symbol]
            self.trades = (self.trades + [last])[-self.lookback:]
            if len(self.trades) >= self.lookback:
                m,b = np.polyfit(range(len(self.trades)),self.trades,1)
                pred = (self.lookback + 1) * m + b
                holdings = self.getHoldings(self.symbol)
                # bid, _, ask, _ = self.getKnownBidAsk(self.symbol)
                if pred > last:
                    self.placeLimitOrder(self.symbol, quantity=self.size, is_buy_order=True, limit_price=last + 1)
                else:
                    if self.symbol in self.holdings and self.holdings[self.symbol] > 0:
                        order_size = min(self.size, self.holdings[self.symbol])
                        self.placeLimitOrder(self.symbol, quantity=order_size, is_buy_order=False, limit_price=last - 1)
            self.setWakeup(currentTime + pd.Timedelta ("1m"))
            self.state = "AWAITING_WAKEUP"

    def getWakeFrequency(self):
        return pd.Timedelta(self.wake_up_freq)


text:
import torch
from torch.nn import Linear,Tanh,Sequential,ReLU,Softplus
import torch.nn.functional as F
from torch.autograd import Variable
import torch.utils.data as Data
import torch.nn as nn
import torch.nn.functional as func
import random
from pyDOE import lhs
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import time
from scipy.interpolate import griddata
import re
import os
import os.path
import operator
import sys
import dataset
import configure as config
import HCProjection as HCP
import simulation as simu

seed = config.seed
torch.random.manual_seed(seed)
np.random.seed(seed)
torch.cuda.manual_seed_all(seed)

device = config.device
BATCH_SIZE_hyper = config.BATCH_SIZE_hyper
use_HCP = config.use_HCP
use_ANN = config.use_ANN
num_epoch = config.num_epoch_hyper
Nf_hyper = config.Nf_hyper
N_boun_hyper = config.N_boun_hyper
nx = config.nx     #网格个数
ny = config.ny
nt = config.nt      #时间步数

k = simu.k
x = simu.x
y = simu.y
t = simu.t
logk = simu.logk
h_boun1_0 = simu.h_boun1_0-200
h_boun2_0 = simu.h_boun2_0-200
hh = simu.hh-200
Ss = simu.Ss

H_boun_train = dataset.H_boun_train
H_boun2_train = dataset.H_boun2_train
H_ic_train = dataset.H_ic_train
H_ic2_train = dataset.H_ic2_train
Y = dataset.Y
X = dataset.X
T = dataset.T
TXY_boun_train = dataset.TXY_boun_train
TXY_boun2_train = dataset.TXY_boun2_train
TXY_noflow1 = dataset.TXY_noflow1
TXY_noflow1_ghost = dataset.TXY_noflow1_ghost
TXY_noflow2 = dataset.TXY_noflow2
TXY_noflow2_ghost = dataset.TXY_noflow2_ghost
TXY_noflow = dataset.TXY_noflow
TXY_noflow_ghost = dataset.TXY_noflow_ghost
TXY_ic_train = dataset.TXY_ic_train
TXY_ic2_train = dataset.TXY_ic2_train
n_train0 = dataset.n_train0
n_colloc = dataset.n_colloc
TXY_f_train1 = dataset.TXY_f_train1
TXY_f_train2 = dataset.TXY_f_train2 # 去除边界后的配点，对应所有时间步
TXY_f_train3 = dataset.TXY_f_train3
TXY_f_train4 = dataset.TXY_f_train4
TXY_f_train5 = dataset.TXY_f_train5
TXY_f_train6 = dataset.TXY_f_train6
TXYK_kxky_f_train = dataset.TXYK_kxky_f_train
TXY_train = dataset.TXY_train
H_train = dataset.H_train

#########################################################
#建立模型
#########################################################
#设置保留小数位数
torch.set_printoptions(precision=7, threshold=None, edgeitems=None, linewidth=None, profile=None)
#定义神经网络
Net_pinn=Sequential(
    Linear(3,50),
    #torch.nn.BatchNorm1d(30, momentum=0.5),
    Softplus(beta=1, threshold=20),
    Linear(50,50),
    #torch.nn.BatchNorm1d(30, momentum=0.5),
    Softplus(beta=1, threshold=20),
    Linear(50,50),
    #torch.nn.BatchNorm1d(30, momentum=0.5),
    Softplus(beta=1, threshold=20),
    Linear(50, 50),
    #torch.nn.BatchNorm1d(30, momentum=0.5),
    Softplus(beta=1, threshold=20),
    Linear(50, 50),
    #torch.nn.BatchNorm1d(30, momentum=0.5),
    Softplus(beta=1, threshold=20),
    Linear(50, 50),
    #torch.nn.BatchNorm1d(30, momentum=0.5),
    Softplus(beta=1, threshold=20),
    Linear(50, 50),
    #torch.nn.BatchNorm1d(30, momentum=0.5),
    Softplus(beta=1, threshold=20),
    Linear(50, 1),
)
#网络参数初始化
def init_weights(m):
    if type(m) == Linear:
        torch.nn.init.xavier_uniform(m.weight)
        m.bias.data.fill_(0.01)
#初始化网络
Net_pinn.apply(init_weights)
#定义Loss
loss_pinn_set=[]
f1_set=[]
f2_set=[] # AH
f2updated_set=[] # 映射后的AH_
h_diff_set=[] # (H-H_)/H
diff0_set = [] # a0-pde0
diff1_set = [] # a1-pde1
diff2_set = [] # a2-pde2
pde012_set = []
a012_set = []
f2_pde_set=[]
f3_set=[]
f4_set=[]
f5_set=[]
f6_set=[]
f7_set=[]
f7_hc_set=[]
hc_loss_set=[]
hc_loss_project_set=[]
hc_loss_std_set=[]
hc_loss_project_std_set=[]
gradient_diff_set=[]
start_time = time.time()
Net_pinn=Net_pinn.to(device)
# 定义神经网络优化器
optimizer2=torch.optim.Adam([
    {'params': Net_pinn.parameters()},
])

#########################################################
#训练模型
#########################################################
#训练数据加载与分批
BATCH_SIZE = BATCH_SIZE_hyper     # 批训练的数据个数
torch_dataset = Data.TensorDataset(TXY_train, H_train) # 前t_train个有观测值的时间步的数据
loader = Data.DataLoader(
    dataset=torch_dataset,      # torch TensorDataset format
    batch_size=BATCH_SIZE,      # mini batch size
    shuffle=True,               # 要不要打乱数据 (打乱比较好)
    num_workers=0,              # 多线程来读数据
)
N_batch=math.ceil(n_train0/BATCH_SIZE) # 观测值数目/batch size获得批数目
BATCH_SIZE2=math.ceil(n_colloc/N_batch ) # 根据批数目计算每批中的配点数目
TXY_f_train_set1=[]
TXY_f_train_set2=[]
TXY_f_train_set3=[]
TXY_f_train_set4=[]
TXY_f_train_set5=[]
TXY_f_train_set6=[]
if use_HCP:
    for i_batch in range(int(N_batch)):
        TXY_f_train_set1.append(TXY_f_train1[BATCH_SIZE2*i_batch:BATCH_SIZE2*(i_batch+1),:])
        TXY_f_train_set2.append(TXY_f_train2[BATCH_SIZE2*i_batch:BATCH_SIZE2*(i_batch+1),:])
        TXY_f_train_set3.append(TXY_f_train3[BATCH_SIZE2*i_batch:BATCH_SIZE2*(i_batch+1),:])
        TXY_f_train_set4.append(TXY_f_train4[BATCH_SIZE2*i_batch:BATCH_SIZE2*(i_batch+1),:])
        TXY_f_train_set5.append(TXY_f_train5[BATCH_SIZE2*i_batch:BATCH_SIZE2*(i_batch+1),:])
        TXY_f_train_set6.append(TXY_f_train6[BATCH_SIZE2*i_batch:BATCH_SIZE2*(i_batch+1),:])

TXYK_kxky_f_train_set=[]
for i_batch in range(int(N_batch)):
    TXYK_kxky_f_train_set.append(TXYK_kxky_f_train[BATCH_SIZE2*i_batch:BATCH_SIZE2*(i_batch+1),:])
start_time = time.time()
#迭代训练
for epoch in range(num_epoch):   # 训练所有!整套!数据
    for ite, (batch_x, batch_y) in enumerate(loader):  # 每一步 loader 释放一小批数据用来学习
        batch_x= Variable(batch_x, requires_grad=True).to(device) # 抽取的有观测的数据，包含边界
        batch_y= Variable(batch_y).to(device) # 抽取的有观测的数据，包含边界
        if use_HCP == True:
            batch_xf1=Variable(TXY_f_train_set1[ite],requires_grad=True).to(device) # 配点数据，全部时间步，无边界
            batch_xf2=Variable(TXY_f_train_set2[ite],requires_grad=True).to(device) # 配点数据，全部时间步，无边界
            batch_xf3=Variable(TXY_f_train_set3[ite],requires_grad=True).to(device) # 配点数据，全部时间步，无边界
            batch_xf4=Variable(TXY_f_train_set4[ite],requires_grad=True).to(device) # 配点数据，全部时间步，无边界
            batch_xf5=Variable(TXY_f_train_set5[ite],requires_grad=True).to(device) # 配点数据，全部时间步，无边界
            batch_xf6=Variable(TXY_f_train_set6[ite],requires_grad=True).to(device) # 配点数据，全部时间步，无边界
        batch_xf=Variable(TXYK_kxky_f_train_set[ite],requires_grad=True).to(device)
        batch_x2= batch_xf[:,0:3]
        K=batch_xf[:,-3:-2]
        k_x_=batch_xf[:,-2:-1]
        k_y_=batch_xf[:,-1:]
        if use_HCP == True:
            n_train=len(batch_xf1)
        else:
            n_train=len(batch_xf)
        optimizer2.zero_grad()
        # 预测观测

        def HCP_layer_observation (T_test_col, X_test_col, Y_test_col, k):
        # 此函数只用于观测值处。根据渗透率场和核心点坐标，给出经过映射后的H。
        # 输入：核心点位的坐标T_test_col, X_test_col, Y_test_col。
        # 输出：h_pred_pinn, TXY_all, hc_loss, hc_loss_, Hreal_, H_diff, Hreal
            # 产生含有ghost grid的全场坐标
            TXY_all_1, TXY_all, TXY_all_3, TXY_all_4, TXY_all_5, TXY_all_6 = HCP.generate_HCP_area(T_test_col.cpu().detach().numpy(), X_test_col.cpu().detach().numpy(), Y_test_col.cpu().detach().numpy())
            # 获得含有ghost grid的全场预测
            h_pred_pinn_1 =Net_pinn(TXY_all_1) # t-1
            h_pred_pinn = Net_pinn(TXY_all) # 核心点
            h_pred_pinn_3=Net_pinn(TXY_all_3) # x-1
            h_pred_pinn_4=Net_pinn(TXY_all_4) # x+1
            h_pred_pinn_5=Net_pinn(TXY_all_5) # y-1
            h_pred_pinn_6=Net_pinn(TXY_all_6) # y+1
            # H=torch.cat((h_pred_pinn_1_bc.data,     h_pred_pinn_bc.data,        h_pred_pinn_3_bc.data,\
            #             h_pred_pinn_4_bc.data,     h_pred_pinn_5_bc.data,      h_pred_pinn_6_bc.data),1)
            H=torch.cat((h_pred_pinn_1,     h_pred_pinn,        h_pred_pinn_3,\
                        h_pred_pinn_4,     h_pred_pinn_5,      h_pred_pinn_6),1)
            k_origin = k
            k_hc1 = np.vstack((np.vstack((np.zeros(len(k)), k)), np.zeros(len(k)))) # 增加无流条件对应的上下k=0
            k_hc2 = np.hstack((np.hstack((k_hc1[:,0].reshape(-1,1),k_hc1)), k_hc1[:,-1].reshape(-1,1)))
            # HCProjection
            A = HCP.buildA_batch(k_hc2, (torch.round(TXY_all[:, 1:2]*100/20-1)+1).long().cpu(),\
                    (torch.round(TXY_all[:, 2:3]*100/20-1)+1).long().cpu(), \
                    0.2,0.2,0.2, 1) # 坐标中的最后一个+1是因为此处对k进行了扩充（k_hc1和k_hc2），增加了ghost grid，所以原先下边界y=0处实际变成了y=1，同理，左边界x=0处变成了x=1
            # A = HCP.buildA_batch(k, (torch.round(TXY_all[:, 1:2]*100/20-1)).long().cpu(),\
            #         (torch.round(TXY_all[:, 2:3]*100/20-1)).long().cpu(), \
            #         0.2,0.2,0.2, 1)
            H_, _, _, Hreal_ = HCP.projection_batch(A, H.unsqueeze(2))
            H_diff= torch.abs((H-H_)/H).mean()
            Hreal = H.unsqueeze(2)*(h_boun1_0-h_boun2_0)+h_boun2_0
            h_pred_pinn = torch.mm(H_, torch.tensor([[0],[1],[0],[0],[0],[0]]).type(torch.FloatTensor).to(device)) # 经过映射调整后的配点预测值
            return h_pred_pinn
        if use_HCP == True:
            prediction = HCP_layer_observation (batch_x[:,0].reshape(-1,1), batch_x[:,1].reshape(-1,1), batch_x[:,2].reshape(-1,1), k)
        else:
            prediction=Net_pinn(batch_x)

        # 预测配点或者配区域
        if use_HCP == True:
            prediction_f1=Net_pinn(batch_xf1)
            prediction_f2=Net_pinn(batch_xf2)
            prediction_f3=Net_pinn(batch_xf3)
            prediction_f4=Net_pinn(batch_xf4)
            prediction_f5=Net_pinn(batch_xf5)
            prediction_f6=Net_pinn(batch_xf6)
            # HCProjection
            A = HCP.buildA_batch(k, (torch.round(batch_xf2[:, 1:2]*100/20-1)).long().cpu(),\
                    (torch.round(batch_xf2[:, 2:3]*100/20-1)).long().cpu(), \
                    0.2,0.2,0.2, 1)
            H=torch.cat((prediction_f1.data,prediction_f2,prediction_f3.data,\
                        prediction_f4.data,prediction_f5.data,prediction_f6.data),1)
            H_, hc_loss, hc_loss_, Hreal_ = HCP.projection_batch(A, H.unsqueeze(2))
            hc_loss_set.append(torch.mean(torch.abs(hc_loss)))
            hc_loss_std_set.append(torch.std(torch.abs(hc_loss)))
            hc_loss_project_set.append(torch.mean(torch.abs(hc_loss_)))
            hc_loss_project_std_set.append(torch.std(torch.abs(hc_loss)))
            H_diff= torch.abs((H-H_)/H).mean()
            Hreal = H.unsqueeze(2)*(h_boun1_0-h_boun2_0)+h_boun2_0
            #print('mean: (Hreal)[:,1,:]: ', torch.abs((Hreal)[:,1,:]).mean())
            #print('mean: (Hreal_)[:,1,:]: ', torch.abs((Hreal_)[:,1,:]).mean())
            #print('mean: (Hreal - Hreal_)[:,1,:]: ', torch.abs((Hreal-Hreal_)[:,1,:]).mean())
            #print('max: (Hreal - Hreal_)[/H:,1,:]: ', (torch.abs((Hreal-Hreal_)[:,1,:])/torch.abs(Hreal[:,1,:])).max())
            prediction_f = torch.mm(H_, torch.tensor([[0],[1],[0],[0],[0],[0]]).type(torch.FloatTensor).to(device)) # 经过映射调整后的配点预测值
            #print('mean: (H - H_)/H[:,1,:]:',torch.abs(prediction_f - prediction_f2.permute(1,0)).mean())
            # if (epoch+1)%100 == 0:
            #     if ite == 0:
            #         np.savetxt('result_epoch_%d_H.csv'%(epoch+1),np.array(H.detach().cpu()).reshape(-1,6), delimiter=',')
            #         np.savetxt('result_epoch_%d_H_.csv'%(epoch+1),np.array(H_.detach().cpu()).reshape(-1,6), delimiter=',')
        else:
            prediction_f = Net_pinn(batch_x2) # 直接PINN预测

        # 预测边界
        boun_pred=Net_pinn(TXY_boun_train) # 左边界，x=0.2处H=1
        boun2_pred=Net_pinn(TXY_boun2_train) # 右边界，x=10.2处H=0
        noflow_pred=Net_pinn(TXY_noflow)
        if use_HCP:
            noflow_pred_ghost=Net_pinn(TXY_noflow_ghost)
        def HCP_layer_boundary (boundary_type, direction, pressure, T_test_col, X_test_col, Y_test_col, k):
        # 此函数只用于边界条件处。根据渗透率场和核心点坐标，给出经过映射后的H。
        # 输入：boundary_type边界类型,0定压1无流； direction边界位置； pressure只用于定压边界；核心点位的坐标T_test_col, X_test_col, Y_test_col。
        # 输出：h_pred_pinn, TXY_all, hc_loss, hc_loss_, Hreal_, H_diff, Hreal
            # 产生含有ghost grid的全场坐标
            TXY_all_1, TXY_all, TXY_all_3, TXY_all_4, TXY_all_5, TXY_all_6 = HCP.generate_HCP_area(T_test_col.cpu().detach().numpy(), X_test_col.cpu().detach().numpy(), Y_test_col.cpu().detach().numpy())
            # 获得含有ghost grid的全场预测
            h_pred_pinn_1 =Net_pinn(TXY_all_1) # t-1
            h_pred_pinn = Net_pinn(TXY_all) # 核心点
            h_pred_pinn_3=Net_pinn(TXY_all_3) # x-1
            h_pred_pinn_4=Net_pinn(TXY_all_4) # x+1
            h_pred_pinn_5=Net_pinn(TXY_all_5) # y-1
            h_pred_pinn_6=Net_pinn(TXY_all_6) # y+1
            # 设置边界条件
            # 根据定压和无流条件，生成H场
            if boundary_type == 0: # 定压
                h_pred_pinn_1_bc, h_pred_pinn_bc, h_pred_pinn_3_bc, h_pred_pinn_4_bc, h_pred_pinn_5_bc, h_pred_pinn_6_bc = HCP.H_set_boundary_pressure(direction, pressure, h_pred_pinn_1, h_pred_pinn, h_pred_pinn_3, h_pred_pinn_4, h_pred_pinn_5, h_pred_pinn_6)
            if boundary_type == 1: # 无流
                h_pred_pinn_1_bc, h_pred_pinn_bc, h_pred_pinn_3_bc, h_pred_pinn_4_bc, h_pred_pinn_5_bc, h_pred_pinn_6_bc = HCP.H_set_boundary_noflow(direction, h_pred_pinn_1, h_pred_pinn, h_pred_pinn_3, h_pred_pinn_4, h_pred_pinn_5, h_pred_pinn_6)
            # H=torch.cat((h_pred_pinn_1_bc.data,     h_pred_pinn_bc.data,        h_pred_pinn_3_bc.data,\
            #             h_pred_pinn_4_bc.data,     h_pred_pinn_5_bc.data,      h_pred_pinn_6_bc.data),1)
            H=torch.cat((h_pred_pinn_1_bc,     h_pred_pinn_bc,        h_pred_pinn_3_bc,\
                        h_pred_pinn_4_bc,     h_pred_pinn_5_bc,      h_pred_pinn_6_bc),1)
            # 根据定压和无流条件，生成k场
            k_origin = k
            k_hc1 = np.vstack((np.vstack((np.zeros(len(k)), k)), np.zeros(len(k)))) # 增加无流条件对应的上下k=0
            k_hc2 = np.hstack((np.hstack((k_hc1[:,0].reshape(-1,1),k_hc1)), k_hc1[:,-1].reshape(-1,1)))
            # HCProjection
            A = HCP.buildA_batch(k_hc2, (torch.round(TXY_all[:, 1:2]*100/20-1)+1).long().cpu(),\
                    (torch.round(TXY_all[:, 2:3]*100/20-1)+1).long().cpu(), \
                    0.2,0.2,0.2, 1) # 坐标中的最后一个+1是因为此处对k进行了扩充（k_hc1和k_hc2），增加了ghost grid，所以原先下边界y=0处实际变成了y=1，同理，左边界x=0处变成了x=1
            H_, _, _, Hreal_ = HCP.projection_batch(A, H.unsqueeze(2))
            H_diff= torch.abs((H-H_)/H).mean()
            Hreal = H.unsqueeze(2)*(h_boun1_0-h_boun2_0)+h_boun2_0
            h_pred_pinn = torch.mm(H_, torch.tensor([[0],[1],[0],[0],[0],[0]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
            h_pred_pinn_5 = torch.mm(H_, torch.tensor([[0],[0],[0],[0],[1],[0]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
            h_pred_pinn_6 = torch.mm(H_, torch.tensor([[0],[0],[0],[0],[0],[1]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
            # return h_pred_pinn, h_pred_pinn_5, h_pred_pinn_6, TXY_all, hc_loss, hc_loss_, Hreal_, H_diff, Hreal
            return h_pred_pinn, h_pred_pinn_5, h_pred_pinn_6
        # boun_pred, _, _ = HCP_layer_boundary (0, 0, 1, TXY_boun_train[:,0].reshape(-1,1), TXY_boun_train[:,1].reshape(-1,1), TXY_boun_train[:,2].reshape(-1,1), k)
        # boun2_pred, _, _ = HCP_layer_boundary (0, 1, 0, TXY_boun2_train[:,0].reshape(-1,1), TXY_boun2_train[:,1].reshape(-1,1), TXY_boun2_train[:,2].reshape(-1,1), k)
        # noflow_pred1, noflow_pred1_ghost, _ = HCP_layer_boundary (1, 0, -999, TXY_noflow1[:,0].reshape(-1,1), TXY_noflow1[:,1].reshape(-1,1), TXY_noflow1[:,2].reshape(-1,1), k)
        # noflow_pred2, _, noflow_pred2_ghost = HCP_layer_boundary (1, 1, -999, TXY_noflow2[:,0].reshape(-1,1), TXY_noflow2[:,1].reshape(-1,1), TXY_noflow2[:,2].reshape(-1,1), k)
        # noflow_pred = torch.cat((noflow_pred1, noflow_pred2), 0)
        # noflow_pred_ghost = torch.cat((noflow_pred1_ghost, noflow_pred2_ghost), 0)

        # 初始实现
        ic_pred=Net_pinn(TXY_ic_train)
        ic2_pred=Net_pinn(TXY_ic2_train)

        # 配点处自动微分
        if use_HCP:
            H_grad = torch.autograd.grad(outputs=prediction_f.sum(), inputs=batch_xf2, create_graph=True)[0]
        else:
            H_grad = torch.autograd.grad(outputs=prediction_f.sum(), inputs=batch_x2, create_graph=True)[0]
        Ht = H_grad[:, 0].contiguous().view(n_train, 1)
        Hx = H_grad[:, 1].contiguous().view(n_train, 1)
        Hy = H_grad[:, 2].contiguous().view(n_train, 1)
        if use_HCP:
            Hxx=torch.autograd.grad(outputs=Hx.sum(), inputs=batch_xf2,\
                                    create_graph=True)[0][:,1].contiguous().view(n_train, 1)
            Hyy=torch.autograd.grad(outputs=Hy.sum(), inputs=batch_xf2,\
                                    create_graph=True)[0][:,2].contiguous().view(n_train, 1)
        else:
            Hxx=torch.autograd.grad(outputs=Hx.sum(), inputs=batch_x2,\
                                    create_graph=True)[0][:,1].contiguous().view(n_train, 1)
            Hyy=torch.autograd.grad(outputs=Hy.sum(), inputs=batch_x2,\
                                    create_graph=True)[0][:,2].contiguous().view(n_train, 1)
        # 无流边界微分
        H_noflow_grad = torch.autograd.grad(outputs=noflow_pred.sum(), inputs=TXY_noflow, create_graph=True)[0]
        Hy_noflow=H_noflow_grad[:,2:3].contiguous()
        # 对比差分梯度与自动微分差异
        if use_HCP == True:
            gradient_diff, diff_t_ID, diff_x_ID, diff_y_ID = HCP.compare_gradient(H, Ht, Hx, Hy)
            diff0, diff1, diff2, diff, pde012, a012, hc_loss_3_item = HCP.compare_loss (H, k, batch_xf2, Ht, Hxx, Hx, Hyy, Hy, K, k_x_, k_y_ )
        # 计算Loss
        f1=torch.pow((prediction-batch_y),2).mean()
        if use_HCP == True:
            f2=torch.pow(hc_loss_3_item.type(torch.FloatTensor),2).mean()
            #f2=torch.pow(hc_loss_3,2).mean()  # AH
            f2_=torch.pow(hc_loss_,2).mean()
        f2_pde=torch.pow((Ss*Ht-K*Hxx-K*Hyy-k_x_*Hx-k_y_*Hy)*1,2).mean()
        f3=torch.pow((boun_pred-H_boun_train),2).mean()+torch.pow((boun2_pred-H_boun2_train),2).mean()
        f6=torch.pow((ic_pred-H_ic_train),2).mean()+torch.pow((ic2_pred-H_ic2_train),2).mean()
        if use_HCP == True:
            f7_hc=25*torch.pow(torch.abs(noflow_pred-noflow_pred_ghost)*10,2).mean() #基于差分的无流边界条件
        else:
            f7=torch.pow(Hy_noflow*10,2).mean()
        #f7=torch.pow(Hy_noflow*10,2).mean()
        # 组合训练用loss
        ################################################################################
        ################################################################################
        if use_HCP == True:
            loss=1*f1+100*f2_pde+1*f3+1*f6+1*f7_hc
        else:
            if use_ANN == True:
                loss=1*f1
            else:
                loss=1*f1+100*f2_pde+1*f3+1*f6+1*f7
        ################################################################################
        ################################################################################
        # 更新网络
        loss.backward()
        optimizer2.step()
        # 打印各个loss
        print('f1:', f1)
        if use_HCP == True:
            print('f2:', f2)
            print('h_diff:', H_diff)
        print('f2_pde', f2_pde)
        print('f3:', f3)
        print('f6:', f6)
        if use_HCP:
            print('f7_hc:', f7_hc)
        else:
            print('f7:', f7)
        loss=loss.data
        f1=f1.data
        if use_HCP == True:
            f2=f2.data
            f2_=f2_.data
            H_diff=H_diff.data
        f2_pde=f2_pde.data
        f3=f3.data
        f6=f6.data
        if use_HCP:
            f7_hc=f7_hc.data
        else:
            f7=f7.data
        loss_pinn_set.append(loss)
        f1_set.append(f1)
        if use_HCP == True:
            f2_set.append(f2)
            f2updated_set.append(f2_)
            h_diff_set.append(H_diff)
            diff0_set.append(diff0)
            diff1_set.append(diff1)
            diff2_set.append(diff2)
            pde012_set.append(pde012)
            a012_set.append(a012)
        f2_pde_set.append(f2_pde)
        f3_set.append(f3)
        f6_set.append(f6)
        if use_HCP:
            f7_hc_set.append(f7_hc)
        else:
            f7_set.append(f7)
        if use_HCP == True:
            gradient_diff_set.append(gradient_diff)
        print('Epoch: ', epoch, '| Step: ', ite, '|loss: ',loss)
# np.savetxt('result_PDE.csv',np.array(f2_pde_set).reshape(-1,1), delimiter=',')
# np.savetxt('result_HC.csv',np.array(f2_set).reshape(-1,1), delimiter=',')

elapsed = time.time() - start_time
print('Training time: %.4f' % (elapsed))

def HCP_layer(T_test_col, X_test_col, Y_test_col, k):
# 根据渗透率场和核心点坐标，给出经过映射后，符合边界条件的H。
# 注意！此时核心点坐标的排布顺序是有要求的，不是任意顺序的坐标点都可以（因为HCP.H_set_boundary函数对于坐标点排布顺序有要求）
# 输入：核心点位的坐标T_test_col, X_test_col, Y_test_col。
# 输出：h_pred_pinn, TXY_all, hc_loss, hc_loss_, Hreal_, H_diff, Hreal
    # 产生含有ghost grid的全场坐标
    TXY_all_1, TXY_all, TXY_all_3, TXY_all_4, TXY_all_5, TXY_all_6 = HCP.generate_HCP_area(T_test_col, X_test_col, Y_test_col)
    # 获得含有ghost grid的全场预测
    h_pred_pinn_1 =Net_pinn(TXY_all_1) # t-1
    h_pred_pinn = Net_pinn(TXY_all) # 核心点
    h_pred_pinn_3=Net_pinn(TXY_all_3) # x-1
    h_pred_pinn_4=Net_pinn(TXY_all_4) # x+1
    h_pred_pinn_5=Net_pinn(TXY_all_5) # y-1
    h_pred_pinn_6=Net_pinn(TXY_all_6) # y+1
    # 设置边界条件
    # 根据定压和无流条件，生成H场
    h_pred_pinn_1_bc, h_pred_pinn_bc, h_pred_pinn_3_bc, h_pred_pinn_4_bc, h_pred_pinn_5_bc, h_pred_pinn_6_bc = HCP.H_set_boundary(h_pred_pinn_1, h_pred_pinn, h_pred_pinn_3, h_pred_pinn_4, h_pred_pinn_5, h_pred_pinn_6)
    H=torch.cat((h_pred_pinn_1_bc.data,     h_pred_pinn_bc.data,        h_pred_pinn_3_bc.data,\
                 h_pred_pinn_4_bc.data,     h_pred_pinn_5_bc.data,      h_pred_pinn_6_bc.data),1)
    # 根据定压和无流条件，生成k场
    k_origin = k
    k_hc1 = np.vstack((np.vstack((np.zeros(len(k)), k)), np.zeros(len(k)))) # 增加无流条件对应的上下k=0
    k_hc2 = np.hstack((np.hstack((k_hc1[:,0].reshape(-1,1),k_hc1)), k_hc1[:,-1].reshape(-1,1)))
    # HCProjection
    A = HCP.buildA_batch(k_hc2, (torch.round(TXY_all[:, 1:2]*100/20-1)+1).long().cpu(),\
            (torch.round(TXY_all[:, 2:3]*100/20-1)+1).long().cpu(), \
            0.2,0.2,0.2, 1) # 坐标中的最后一个+1是因为此处对k进行了扩充（k_hc1和k_hc2），增加了ghost grid，所以原先下边界y=0处实际变成了y=1，同理，左边界x=0处变成了x=1
    H_, hc_loss, hc_loss_, Hreal_ = HCP.projection_batch(A, H.unsqueeze(2))
    H_diff= torch.abs((H-H_)/H).mean()
    Hreal = H.unsqueeze(2)*(h_boun1_0-h_boun2_0)+h_boun2_0
    h_pred_pinn_1 = torch.mm(H_, torch.tensor([[1],[0],[0],[0],[0],[0]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
    h_pred_pinn = torch.mm(H_, torch.tensor([[0],[1],[0],[0],[0],[0]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
    h_pred_pinn_3 = torch.mm(H_, torch.tensor([[0],[0],[1],[0],[0],[0]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
    h_pred_pinn_4 = torch.mm(H_, torch.tensor([[0],[0],[0],[1],[0],[0]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
    h_pred_pinn_5 = torch.mm(H_, torch.tensor([[0],[0],[0],[0],[1],[0]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
    h_pred_pinn_6 = torch.mm(H_, torch.tensor([[0],[0],[0],[0],[0],[1]]).type(torch.FloatTensor).to(device)) #此处是经过映射后的对应位置的结果
    # 设置边界条件
    _, h_pred_pinn, _, _, _, _ = HCP.H_set_boundary(h_pred_pinn_1, h_pred_pinn, h_pred_pinn_3, h_pred_pinn_4, h_pred_pinn_5, h_pred_pinn_6)
    return h_pred_pinn, TXY_all, hc_loss, hc_loss_, Hreal_, H_diff, Hreal

# PINN predict
H_all=hh.flatten()[:,None]
# 给出核心点坐标（整个流场所有时间步的）
T_test_col=T.flatten()[:,None]
X_test_col=X.flatten()[:,None]
Y_test_col=Y.flatten()[:,None]
if use_HCP == True:
    h_pred_pinn, TXY_all, hc_loss, hc_loss_, Hreal_, H_diff, Hreal = HCP_layer(T_test_col, X_test_col, Y_test_col, k)
else:
    TXY_all = np.hstack((T_test_col,X_test_col,Y_test_col))
    TXY_all = Variable(torch.from_numpy(TXY_all).type(torch.FloatTensor), requires_grad=True).to(device)
    h_pred_pinn = Net_pinn(TXY_all)
h_pred_pinn = (h_pred_pinn*(h_boun1_0-h_boun2_0)+h_boun2_0).cpu().data.numpy()
TXY_all=TXY_all.cpu()
h_array_pred_pinn=h_pred_pinn.reshape(nt,nx,ny)
#  calculate error
if config.use_noise == True:
    hh_origion = simu.hh_origion-200
    error_l2 = np.linalg.norm(hh_origion.flatten()-h_array_pred_pinn.flatten(),2)/np.linalg.norm(hh_origion.flatten(),2)
    print('Error L2: %e' % (error_l2))
    R2=1-np.sum((hh_origion.flatten()-h_array_pred_pinn.flatten())**2)/np.sum((hh_origion.flatten()-hh_origion.flatten().mean())**2)
    print('coefficient of determination  R2: %e' % (R2))
    # # 绘制交汇图
    # noise = config.noise
    # noise_level = '%d'%(noise*100)
    # hh_all = np.array(dataset.H_train)
    # hh_origion_all = dataset.H_train_origion.flatten()
    # # hh_origion_all = (hh_origion_all)/2
    # # idx2 = np.random.choice(range(hh_all.shape[0]), 1000, replace=False)
    # plt.figure(figsize=(8,8))
    # plt.scatter(hh_origion_all,hh_all,marker='o',edgecolors='b')
    # plt.xlabel('Reference',fontsize=20)
    # plt.ylabel('Noisy observation',fontsize=20)
    # plt.title('Noise level '+noise_level+'%',fontsize=20)
    # plt.xlim(-0.1,1+0.1)
    # plt.ylim(-0.3,1+0.3)
    # plt.savefig('Noise level 555555 1800 normal %d.png'%(noise*100), dpi=600)
    # plt.show()
else:
    error_l2 = np.linalg.norm(hh.flatten()-h_array_pred_pinn.flatten(),2)/np.linalg.norm(hh.flatten(),2)
    print('Error L2: %e' % (error_l2))
    R2=1-np.sum((hh.flatten()-h_array_pred_pinn.flatten())**2)/np.sum((hh.flatten()-hh.flatten().mean())**2)
    print('coefficient of determination  R2: %e' % (R2))

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

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

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

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

##################################################################################################################
#绘图
##################################################################################################################
# 目录
# path = "\\4_no_label_parameter_results\\"#切换工作目录
# now = os.getcwd()# 查看当前工作目录
# os.chdir( now+path )# 修改当前工作目录
# new_path = os.getcwd()# 查看修改后的工作目录
# if use_HCP:
#     np.savetxt('result_lastepoch_H.csv',np.array(H.detach().cpu()).reshape(-1,6), delimiter=',')
#     np.savetxt('result_lastepoch_H_.csv',np.array(H_.detach().cpu()).reshape(-1,6), delimiter=',')
#     np.savetxt('result_AH.csv',np.array(hc_loss_set).reshape(-1,1), delimiter=',')
#     np.savetxt('result_AH_.csv',np.array(hc_loss_project_set).reshape(-1,1), delimiter=',')
#     np.savetxt('result_AH_std.csv',np.array(hc_loss_std_set).reshape(-1,1), delimiter=',')
#     np.savetxt('result_AH_std_.csv',np.array(hc_loss_project_std_set).reshape(-1,1), delimiter=',')
#     np.savetxt('result_gradient.csv',gradient_diff_set, delimiter=',')

# # 画局部取值图，看梯度差分计算是否合理
# def NN_value_surface (t, x, y, n=1, n_unit = 10):
# # 功能：计算某一个区域内差分尺度范围的神经网络取值面，可以对比对应配点附近临域自动微分的梯度和配区域差分的梯度是否相似
# # 输入：txy对应配点坐标，n决定绘图区域的长与宽（2n倍的差分间隔），n_unit对应于绘图的精细度（差分间隔分为n_unit份）
# # 输出：value_surface一个尺寸为（2n/delta， 2n/delta）的矩阵，元素为对应位置的神经网预测值
#     x_min = x-n*0.2
#     x_max = x+n*0.2
#     y_min = y-n*0.2
#     y_max = y+n*0.2
#     len_area = 2*n*n_unit+1
#     coordinate = np.zeros((np.int(len_area*len_area), 3))
#     coordinate[:,0] = t
#     coordinate[:,1] = np.tile(np.arange(x_min,x_max, 0.2/n_unit), np.int(len_area))
#     coordinate[:,2] = np.arange(x_min,x_max,0.2/n_unit).repeat(len_area, axis=0)
#     coordinate = Variable(torch.from_numpy(coordinate).type(torch.FloatTensor),requires_grad=True).to(device)
#     prediction_area = Net_pinn(coordinate)
#     prediction_area_square = prediction_area.cpu().detach().numpy().reshape((len_area, len_area))
#     return prediction_area_square
# if config.plot_value_surface == True:
#     for i in range (3):
#         for j in range (2):
#             t_plot = i*20+4
#             y_plot = j*20+14
#             n_unit=50
#             Z_plot =  NN_value_surface (t_plot, 5, y_plot, n=1, n_unit=n_unit)
#             X_plot=np.arange(5-0.2, 5+0.2001, 0.2/n_unit)
#             Y_plot=np.arange(y_plot-0.2, y_plot+0.2001, 0.2/n_unit)
#             fig = plt.figure()
#             from mpl_toolkits.mplot3d import Axes3D
#             ax = Axes3D(fig)
#             X_plot,Y_plot=np.meshgrid(X_plot,Y_plot)
#             plt.title('Time step %d'%(t_plot+1),fontsize=10)
#             ax.plot_surface(X_plot, Y_plot, Z_plot, rstride = 1, cstride = 1, cmap = plt.get_cmap('rainbow'))

# # #画loss图

# plt.figure()
# plt.plot(range(len(loss_pinn_set)),loss_pinn_set)
# plt.xlabel('Iteration')
# plt.ylabel('Loss')
# loss_result_loss = np.array(loss_pinn_set).reshape(-1,1)
# np.savetxt('result_loss.csv',loss_result_loss, delimiter=',')
# plt.savefig('result_loss.png', dpi=600)

# plt.figure()
# plt.plot(range(len(f1_set)),f1_set)
# plt.xlabel('Iteration')
# plt.ylabel('Observation loss')
# #plt.ylim([0,0.001])
# loss_result_f1=np.array(f1_set).reshape(-1,1)
# loss_result_all=np.hstack((loss_result_loss,loss_result_f1))
# np.savetxt('result_f1.csv',loss_result_f1, delimiter=',')
# plt.savefig('result_f1.png', dpi=600)

# if use_HCP == True:
#     plt.figure()
#     plt.plot(range(len(f2_set)),f2_set)
#     plt.xlabel('Iteration')
#     plt.ylabel('HCP loss')
#     #plt.ylim([0,0.001])
#     loss_result_f2AH=np.array(f2_set).reshape(-1,1)
#     loss_result_all=np.hstack((loss_result_all,loss_result_f2AH))
#     np.savetxt('result_f2AH.csv',loss_result_f2AH, delimiter=',')
#     plt.savefig('result_f2AH.png', dpi=600)

#     plt.figure()
#     plt.plot(range(len(f2updated_set)),f2updated_set)
#     plt.xlabel('Iteration')
#     plt.ylabel('Projected HCP loss')
#     loss_result_f2updatedAH=np.array(f2updated_set).reshape(-1,1)
#     loss_result_all=np.hstack((loss_result_all,loss_result_f2updatedAH))
#     np.savetxt('result_f2updatedAH_.csv',loss_result_f2updatedAH, delimiter=',')
#     plt.savefig('result_f2updatedAH.png', dpi=600)

#     plt.figure()
#     plt.plot(range(len(h_diff_set)),h_diff_set)
#     plt.xlabel('Iteration')
#     plt.ylabel('HCP change ratio')
#     #plt.ylim([0,0.005])
#     np.savetxt('result_h_diff.csv',np.array(h_diff_set).reshape(-1,1), delimiter=',')
#     plt.savefig('result_h_diff.png', dpi=600)

#     plt.figure()
#     plt.plot(range(len(diff0_set)),diff0_set, label='a0-pde0')
#     plt.plot(range(len(diff1_set)),diff1_set, label='a1-pde1')
#     plt.plot(range(len(diff2_set)),diff2_set, label='a2-pde2')
#     plt.xlabel('Iteration')
#     plt.ylabel('diff:(a0-pde0)(a1-pde1)(a2-pde2)')
#     plt.legend()

#     plt.figure()
#     plt.plot(range(len(pde012_set)),pde012_set)
#     plt.xlabel('Iteration')
#     plt.ylabel('pde012')
#     plt.legend()

#     plt.figure()
#     plt.plot(range(len(a012_set)),a012_set)
#     plt.xlabel('Iteration')
#     plt.ylabel('a012')
#     plt.legend()
# plt.figure()
# plt.plot(range(len(f2_pde_set)),f2_pde_set)
# plt.xlabel('Iteration')
# plt.ylabel('PDE loss')
# #plt.ylim([0,0.001])
# loss_result_f2pde=np.array(f2_pde_set).reshape(-1,1)
# loss_result_all=np.hstack((loss_result_all,loss_result_f2pde))
# np.savetxt('result_f2pde.csv',loss_result_f2pde, delimiter=',')
# plt.savefig('result_f2pde.png', dpi=600)
# if use_HCP == True:
#     plt.figure(figsize=(5,5))
#     plt.scatter(np.array(f2_pde_set),np.array(f2_set),marker='o',edgecolors='b')
#     plt.xlabel('pde',fontsize=18)
#     plt.ylabel('hc',fontsize=18)
#     plt.title("PDE v.s. HCP",fontsize=18)

# plt.figure()
# plt.plot(range(len(f3_set)),f3_set)
# plt.xlabel('Iteration')
# plt.ylabel('Constant pressure boundary loss')
# #plt.ylim([0,0.05])
# loss_result_f3=np.array(f3_set).reshape(-1,1)
# loss_result_all=np.hstack((loss_result_all,loss_result_f3))
# np.savetxt('result_f3.csv',loss_result_f3, delimiter=',')
# plt.savefig('result_f3.png', dpi=600)

# '''
# plt.figure()
# plt.plot(range(len(f4_set)),f4_set)
# plt.xlabel('Iteration')
# plt.ylabel('f4_loss')
# plt.ylim([0,0.05])
# plt.figure()
# plt.plot(range(len(f5_set)),f5_set)
# plt.xlabel('Iteration')
# plt.ylabel('f5_loss')
# plt.ylim([0,0.05])
# '''

# plt.figure()
# plt.plot(range(len(f6_set)),f6_set)
# plt.xlabel('Iteration')
# plt.ylabel('Initial condition loss')
# loss_result_f6=np.array(f6_set).reshape(-1,1)
# loss_result_all=np.hstack((loss_result_all,loss_result_f6))
# np.savetxt('result_f6.csv',loss_result_f6, delimiter=',')
# plt.savefig('result_f6.png', dpi=600)
# #plt.ylim([0,0.05])
# # plt.figure()
# # plt.plot(range(len(f7_set)),f7_set)
# # plt.xlabel('Iteration')
# # plt.ylabel('f7_loss')
# #plt.ylim([0,0.05])
# if use_HCP:
#     plt.figure()
#     plt.plot(range(len(f7_hc_set)),f7_hc_set)
#     plt.xlabel('Iteration')
#     plt.ylabel('No flow boundary loss')
#     loss_result_f7=np.array(f7_hc_set).reshape(-1,1)
#     loss_result_all=np.hstack((loss_result_all,loss_result_f7))
#     #plt.ylim([0,0.005])
#     np.savetxt('result_f7_hc.csv',loss_result_f7, delimiter=',')
#     plt.savefig('result_f7_hc.png', dpi=600)
# else:
#     plt.figure()
#     plt.plot(range(len(f7_set)),f7_set)
#     plt.xlabel('Iteration')
#     plt.ylabel('No flow boundary loss')
#     loss_result_f7=np.array(f7_set).reshape(-1,1)
#     loss_result_all=np.hstack((loss_result_all,loss_result_f7))
#     #plt.ylim([0,0.005])
#     np.savetxt('result_f7.csv',loss_result_f7, delimiter=',')
#     plt.savefig('result_f7.png', dpi=600)
# # plt.figure(figsize=(5,5))
# # plt.scatter(np.array(f7_set),np.array(f7_hc_set),marker='o',edgecolors='b')
# # plt.xlabel('pde',fontsize=18)
# # plt.ylabel('hc',fontsize=18)
# # plt.title("pde v.s. hc",fontsize=18)
# np.savetxt('result_loss_all.csv',loss_result_all, delimiter=',')

# # PINN plot
# x=x*100
# y=y*100
# for i_t in range(nt-31,nt,10):
#     # Row 0: h(t,x)
#     plt.figure(figsize=(3,3))
#     mm1=plt.imshow(hh[i_t], interpolation='nearest', cmap='rainbow',
#                   extent=[x.min(), x.max(), y.min(),y.max()],
#                   origin='lower')
#     plt.xlabel('$x$')
#     plt.ylabel('$y$')
#     plt.title('Reference $H$ at time step %d'%(i_t+1),fontsize=10)
#     plt.colorbar(mm1,fraction=0.046, pad=0.04)
#     line = np.linspace(x.min(), x.max(), 2)[:,None]
#     plt.plot(line, y[15]*np.ones((2,1)), 'w-', linewidth = 1)
#     plt.plot(line, y[30]*np.ones((2,1)), 'w-', linewidth = 1)
#     plt.plot(line, y[45]*np.ones((2,1)), 'w-', linewidth = 1)
#     plt.savefig('Reference $H$ at time step %d .png'%(i_t+1), dpi=600)
#     # Row 1: pred h(t,x)
#     plt.figure(figsize=(3,3))
#     mm2=plt.imshow(h_array_pred_pinn[i_t], interpolation='nearest', cmap='rainbow',
#                   extent=[x.min(), x.max(), y.min(),y.max()],
#                   origin='lower')
#     plt.xlabel('$x$')
#     plt.ylabel('$y$')
#     plt.title('Predicted $H$ at time step %d'%(i_t+1),fontsize=10)
#     plt.colorbar(mm2,fraction=0.046, pad=0.04)
#     line = np.linspace(x.min(), x.max(), 2)[:,None]
#     plt.plot(line, y[15]*np.ones((2,1)), 'w-', linewidth = 1)
#     plt.plot(line, y[30]*np.ones((2,1)), 'w-', linewidth = 1)
#     plt.plot(line, y[45]*np.ones((2,1)), 'w-', linewidth = 1)
#     plt.savefig('Predicted $H$ at time step %d .png'%(i_t+1), dpi=600)
#     # Row 2: u(t,x) slices
#     plt.figure(figsize=(10,5))
#     gs1 = gridspec.GridSpec(1, 3)
#     # gs1.update(top=-0.05, bottom=-0.5, left=0.1, right=0.9, wspace=0.6)
#     # 320处
#     ax = plt.subplot(gs1[0, 0])
#     ax.plot(x,hh[i_t,15,:], 'b-', linewidth = 2, label = 'Reference')
#     ax.plot(x,h_array_pred_pinn[i_t,15,:], 'r--', linewidth = 2, label = 'Prediction')
#     ax.set_xlabel('$x$')
#     ax.set_ylabel('$H$')
#     ax.set_title('$y = 320$', fontsize = 10)
#     ax.set_xlim([min(x)-0.1*max(x),max(x)+0.1*max(x)])
#     # 620处
#     ax = plt.subplot(gs1[0, 1])
#     ax.plot(x,hh[i_t,30,:], 'b-', linewidth = 2, label = 'Reference')
#     ax.plot(x,h_array_pred_pinn[i_t,30,:], 'r--', linewidth = 2, label = 'Prediction')
#     ax.set_xlabel('$x$')
#     #ax.set_ylabel('$H$')
#     ax.set_xlim([min(x)-0.1*max(x),max(x)+0.1*max(x)])
#     ax.set_title('$y = 620$', fontsize = 10)
#     ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.35), ncol=5, frameon=False)
#     # 920处
#     ax = plt.subplot(gs1[0, 2])
#     ax.plot(x,hh[i_t,45,:], 'b-', linewidth = 2, label = 'Reference')
#     ax.plot(x,h_array_pred_pinn[i_t,45,:], 'r--', linewidth = 2, label = 'Prediction')
#     ax.set_xlabel('$x$')
#     #ax.set_ylabel('$H$')
#     ax.set_xlim([min(x)-0.1*max(x),max(x)+0.1*max(x)])
#     ax.set_title('$y = 920$', fontsize = 10)
#     plt.savefig('Slices $H$ at time step %d .png'%(i_t+1), dpi=600)
#     # plt.show()
#     slice_h = np.hstack((hh[i_t,15,:].reshape(-1,1), hh[i_t,30,:].reshape(-1,1)))
#     slice_h = np.hstack((slice_h, hh[i_t,45,:].reshape(-1,1)))
#     slice_h = np.hstack((slice_h, h_array_pred_pinn[i_t,15,:].reshape(-1,1)))
#     slice_h = np.hstack((slice_h, h_array_pred_pinn[i_t,30,:].reshape(-1,1)))
#     slice_h = np.hstack((slice_h, h_array_pred_pinn[i_t,45,:].reshape(-1,1)))
#     np.savetxt('result_slice_h_%d.csv'%(i_t+1),slice_h, delimiter=',')

# lastepoch_prediction = h_array_pred_pinn[-1,:,:]
# lastepoch_groundtruth = hh[-1,:,:]
# np.savetxt('result_lastepoch_prediction.csv',lastepoch_prediction, delimiter=',')
# np.savetxt('result_lastepoch_groundtruth.csv',lastepoch_groundtruth, delimiter=',')

# time30_prediction = h_array_pred_pinn[-20,:,:]
# time30_groundtruth = hh[-20,:,:]
# np.savetxt('result_lastepoch_time30_prediction.csv',time30_prediction, delimiter=',')
# np.savetxt('result_lastepoch_time30_groundtruth.csv',time30_groundtruth, delimiter=',')

# # 绘制交汇图
# idx2 = np.random.choice(range(H_all.shape[0]), 300, replace=False)
# plt.figure(figsize=(5,5))
# plt.scatter(H_all[idx2],h_pred_pinn[idx2],marker='o',edgecolors='b')
# plt.xlabel('Reference',fontsize=18)
# plt.ylabel('Prediction',fontsize=18)
# if use_HCP == True:
#     plt.title("Results Comparison for HCP",fontsize=18)
# else:
#     if use_ANN == True:
#         plt.title("Results Comparison for ANN",fontsize=18)
#     else:
#         plt.title("Results Comparison for TgNN",fontsize=18)
# plt.xlim(h_boun2_0-0.2,h_boun1_0+0.2)
# plt.ylim(h_boun2_0-0.2,h_boun1_0+0.2)
# plt.savefig('Results Comparison for TgNN.png', dpi=600)

# # 绘制渗透率分布
# plt.figure(figsize=(3,3))
# mm=plt.imshow(logk,origin='lower',cmap='jet',interpolation='nearest',aspect='equal',
#            extent=[0, 51, 0,51])
# plt.xlabel('$x$')
# plt.ylabel('$y$')
# plt.title("ln$K(x,y)$")
# plt.xlim(0,51)
# plt.ylim(0,51)
# plt.colorbar(mm,fraction=0.046, pad=0.04)
# my_x_ticks = np.arange(0,51, 10)
# my_y_ticks = np.arange(0,51, 10)
# plt.xticks(my_x_ticks)
# plt.yticks(my_y_ticks)
# plt.savefig('K.png', dpi=600)

# plt.show()

print('HCP Mission complete')

reason: import numpy,from scipy


text:
import datetime

import numpy as np
from matplotlib import pyplot as plt
from tqdm import tqdm

def plotter(traces: list, param_intervals: dict = {}, density: int = 0, comparison_set=None, savename=None, timestep=1,
            case=0):
    # TODO add things like axes, title, etc
    fig, axs = plt.subplots(nrows=1, ncols=1, constrained_layout=True, figsize=(16, 9))

    for trace in tqdm(traces, desc='Samples plotted', position=0, leave=True):

        parts = splitter_helper(trace)
        if parts:  # this is for the _constraint tweak.. in case splitter_helper returns none, nothing happens for this trace
            for part in parts:
                axs.plot(part[0], part[1], 'b')
        # axs.scatter(traces[key][0][0], traces[key][1][0])
        # axs[1].plot(traces[(duration, slope, (x,y))][0],traces[(duration, slope, (x,y))][1])

        # plt.plot(traces[key][0],traces[key][1],zorder=0)
    if param_intervals and density:
        param_intervals_string_list = [param + ':' + str(param_intervals[param]) for param in param_intervals.keys()]
        param_intervals_string = ",".join(param_intervals_string_list)

        axs.set_xlabel(param_intervals_string + ', density = ' + str(density))

    if comparison_set:
        # axs.plot(np.linspace(0,100,100), np.zeros(100),'c')
        if case == 1:
            for trace in comparison_set:
                axs.plot(trace[:, 0], trace[:, 1], 'c')

        elif case == 0:
            for trace in comparison_set:
                x_vals = np.arange(len(trace))
                axs.plot(x_vals, trace, 'c')  # this is for xvals [0,1,...]

    axs.set_xlabel('time')
    # axs.set_title('Pulse Specification')

    time = str(datetime.datetime.now())
    timestamp = time.replace(':', '_')  # ':' cannot be in file_path

    if savename:
        plt.savefig(savename)
    # plt.savefig(timestamp + '.jpg') # saves plot
    plt.show()

    return fig

def splitter_helper(trace):
    # this is just to get rid of unnecessary vertical lines between two letters
    # it takes the trace apart at repeating time values and returns a list of smaller traces

    # look for where the duplicates appear. this is where the two atomics are joined together
    # if trace[0] == 'discarded':
    #     return None
    between_letters = [0, ]
    for idx in range(trace[0].size):
        if idx == trace[0].size - 1:
            pass
        else:
            if trace[0][idx] == trace[0][idx + 1]:
                between_letters.append(idx + 1)

    trace_parts = []
    for idx in range(len(between_letters) - 1):
        start_idx = between_letters[idx]
        end_idx = between_letters[idx + 1]  # +1 because of python slicing (last index isnt in slice)
        trace_part = np.asarray((trace[0][start_idx:end_idx], trace[1][start_idx:end_idx]))
        trace_parts.append(trace_part)

    # handle last part separately:
    trace_parts.append(np.asarray((trace[0][between_letters[-1]:], trace[1][between_letters[-1]:])))

    return trace_parts

max_stars_repo_name: figlerg/ShapEx

max_issues_repo_name: figlerg/ShapEx

max_forks_repo_name: figlerg/ShapEx


text:
import numpy as np
import ctypes
from collections import OrderedDict
from multiprocessing import Lock, Value, Array

class SharedParams(object):
    """the parameters of a model as numpy arrays arrays are stored in a dict
    {parameter_name : parameter_values}
    these key/value pairs are stored in order of insertion
    """

    def __init__(self, initial, locked=True, dtype="float32"):
        if not isinstance(initial, OrderedDict):
            raise ValueError("SharedParams requires OrderedDict as initial parameters!")
        self._locked = locked
        self._dtype = dtype
        self._ordered_keys = []
        self._arrays = {}
        self._shapes = {}
        for k, v in initial.items():
            self._ordered_keys.append(k)
            self._arrays[k] = Array(ctypes.c_float, v.flatten())
            self._shapes[k] = v.shape

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

    def __contains__(self, key):
        return key in self._ordered_keys

    def __getitem__(self, key):
        if self._locked:
            with self._arrays[key].get_lock():
                ret = np.frombuffer(self._arrays[key].get_obj(), dtype=self._dtype) \
                    .reshape(self._shapes[key])
        else:
            ret = np.frombuffer(self._arrays[key].get_obj(), dtype=self._dtype) \
                .reshape(self._shapes[key])
        return ret

    def __setitem__(self, key, value):
        """sets the value of the given key to the given value
           currently does not allow adding values after creation
           or altering the dimension of the value
        """
        if key not in self._ordered_keys:
            raise ValueError("Adding additional data is not supported (yet)!")
        if value.shape != self._shapes[key]:
            raise ValueError("new value has shape {}, expected {}".format(value.shape,
                                                                          self._shapes[key]))

        if self._locked:
            with self._arrays[key].get_lock():
                np.frombuffer(self._arrays[key].get_obj(), dtype=self._dtype)[:] = \
                    value.astype(self._dtype).flatten()
        else:
            np.frombuffer(self._arrays[key].get_obj(), dtype=self._dtype)[:] = \
                value.astype(self._dtype).flatten()

    def __delitem__(self, key):
        self._ordered_keys.remove(key)
        self._arrays.__delitem__(key)
        self._shapes.__delitem__(key)

    def __iter__(self):
        return self.keys()

    def __repr__(self):
        return "SharedParams: {} parameters with shapes: {}"\
            .format(len(self), self._shapes)

    def keys(self):
        for k in self._ordered_keys:
            yield k

    def values(self):
        for k in self.keys():
            yield self[k]

    def items(self):
        for k in self.keys():
            yield (k, self[k])

    def as_dict(self):
        """returns an OrderedDict of parameters with the *current* values"""
        d = OrderedDict()
        for k in self:
            d[k] = self[k]
        return d

class SharedFloat(object):

    def __init__(self, val=0.):
        self.val = Value("f", val)
        self.lock = Lock()

    def set_value(self, val):
        with self.lock:
            self.val.value = val
        return self

    @property
    def value(self):
        with self.lock:
            return self.val.value

class SharedCounter(object):

    def __init__(self, val=0):
        self.val = Value("i", val)
        self.lock = Lock()

    def increment(self, val=1):
        with self.lock:
            self.val.value += val
            return self.val.value

    @property
    def value(self):
        with self.lock:
            return self.val.value

max_stars_repo_name: valentindey/async-train

max_issues_repo_name: valentindey/async-train

max_forks_repo_name: valentindey/async-train


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

from pycompss.api.task import task
from pycompss.api.api import compss_barrier

from ddf_library.context import COMPSsContext
from ddf_library.utils import generate_info
from ddf_library.functions.ml.feature import StringIndexer, StandardScaler
from ddf_library.columns import col, udf
from ddf_library.types import DataType

import pandas as pd
import numpy as np
import time
import sys
import math
from io import StringIO

####
#   Data generation
####

@task(returns=2)
def generate_partition(x, multiplier, frag):
    x = x * multiplier
    cols = ['PassengerId', 'Survived', 'Pclass', 'Name', 'Sex', 'Age',
            'SibSp', 'Parch', 'Ticket', 'Fare', 'Cabin', 'Embarked']
    dtypes = {'PassengerId': np.int16, 'Survived': np.int8, 'Pclass': np.int8,
              'Name': np.dtype('O'), 'Sex': np.dtype('O'), 'Age': np.float32,
              'SibSp': np.int8, 'Parch': np.int8, 'Ticket': np.dtype('O'),
              'Fare': np.float32, 'Cabin': np.dtype('O'),
              'Embarked': np.dtype('O')}
    df = pd.read_csv(StringIO(x), sep=',',
                     names=cols, header=None, dtype=dtypes)
    schema = generate_info(df, frag)
    return df, schema

def generate_data(n_repeat, nfrag):
    dfs = [[] for _ in range(nfrag)]
    info = [[] for _ in range(nfrag)]
    block = ''
    with open('base_titanic.csv', 'r') as content_file:
        block = content_file.read()

    for f in range(nfrag):
        dfs[f], info[f] = generate_partition(block, n_repeat, f)
    return dfs, info

def titanic(ddf1):

    def title_checker(name):
        titles = {"Mr.": 1, "Miss": 2, "Mrs.": 3, "Master": 4, "Rare": 5}
        for title in titles:
            if title in name:
                return titles[title]
        return -1

    def age_categorizer(age):
        category = 7
        if age <= 11:
            category = 0
        elif (age > 11) and (age <= 18):
            category = 1
        elif (age > 18) and (age <= 22):
            category = 2
        elif (age > 22) and (age <= 27):
            category = 3
        elif (age > 27) and (age <= 33):
            category = 4
        elif (age > 33) and (age <= 40):
            category = 5
        elif (age > 40) and (age <= 66):
            category = 6

        return category

    def fare_categorizer(fare):
        category = 5
        if fare <= 7.91:
            category = 0
        elif (fare > 7.91) and (fare <= 14.454):
            category = 1
        elif (fare > 14.454) and (fare <= 31):
            category = 2
        elif (fare > 31) and (fare <= 99):
            category = 3
        elif (fare > 99) and (fare <= 250):
            category = 4
        return category

    title_checker_udf = udf(title_checker, DataType.INT, col('Name'))
    age_categorizer_udf = udf(age_categorizer, DataType.INT, col('Age'))
    fare_categorizer_udf = udf(fare_categorizer, DataType.INT, col('Fare'))

    features = ['Pclass', 'Sex', 'Age', 'SibSp', 'Parch', 'Fare', 'Embarked']
    all_columns = features + ['Survived']

    ddf2 = StringIndexer() \
        .fit_transform(ddf1, input_col=['Embarked'], output_col=['Embarked']) \
        .drop(['PassengerId', 'Cabin', 'Ticket']) \
        .dropna(all_columns, how='any') \
        .replace({'male': 1, 'female': 0}, subset=['Sex']) \
        .map(title_checker_udf, 'Name') \
        .map(age_categorizer_udf, 'Age') \
        .map(fare_categorizer_udf, 'Fare')  # .cache()

    features.append('Name')
    ddf3 = StandardScaler(with_mean=True, with_std=True) \
        .fit_transform(ddf2, input_col=features, output_col=features)
    ddf3.save.csv("hdfs://localhost:9000/titanic-out", mode='overwrite')

if __name__ == "__main__":

    print("Titanic workflow implemented to DDF v0.6")

    size_block = int(sys.argv[1]) * 1024 * 1024  # in bytes,
    n_frag = int(sys.argv[2])
    artificial_data = size_block > 0
    cc = COMPSsContext()

    if artificial_data:
        # if working with artificial data:
        size_original = 61113
        n_repeat = int(math.ceil(size_block / size_original))
        ddf1 = cc.import_compss_data(generate_data(n_repeat, n_frag))
    else:

        schema = {'PassengerId': DataType.INT, 'Survived': DataType.INT,
                  'Pclass': DataType.INT, 'Name': DataType.STRING,
                  'Sex': DataType.STRING, 'Age': DataType.DECIMAL,
                  'SibSp': DataType.INT, 'Parch': DataType.INT,
                  'Ticket': DataType.STRING, 'Fare': DataType.DECIMAL,
                  'Cabin': DataType.STRING, 'Embarked': DataType.STRING}
        cc = COMPSsContext()
        ddf1 = cc.read.csv("hdfs://localhost:9000/titanic.csv",
                           schema=schema, sep='\t', num_of_parts=n_frag)

    titanic(ddf1)

    cc.stop()

max_stars_repo_name: eubr-bigsea/Compss-Python

max_issues_repo_name: eubr-bigsea/Compss-Python

max_forks_repo_name: eubr-bigsea/Compss-Python


text:
# tests/tagifai/test_data.py
# Test tagifai/data.py components.

import numpy as np
import pandas as pd
import pytest

from iris import data

def test_train_test_split():
    # Process
    path = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/iris.csv"
    df = pd.read_csv(path, header=None, names=["f1", "f2", "f3", "f4", "class"])
    df = df.sample(frac=1).reset_index(drop=True)
    # df["class"] = LabelEncoder().fit_transform(df["class"])
    train_df = df[: int(len(df) * 0.8)].reset_index(drop=True)
    test_df = df[int(len(df) * 0.8) : int(len(df) * 0.9)].reset_index(drop=True)
    val_df = df[int(len(df) * 0.9) :].reset_index(drop=True)
    X_train, y_train = train_df.values[:, :-1], train_df.values[:, -1]
    X_val, y_val = val_df.values[:, :-1], val_df.values[:, -1]
    X_test, y_test = test_df.values[:, :-1], test_df.values[:, -1]

    assert len(X_train) == len(y_train)
    assert len(X_val) == len(y_val)
    assert len(X_test) == len(y_test)
    assert len(X_train) / float(len(df)) == pytest.approx(0.8, abs=0.05)  # 0.8 ± 0.05
    assert len(X_val) / float(len(df)) == pytest.approx(0.1, abs=0.05)  # 0.1 ± 0.05
    assert len(X_test) / float(len(df)) == pytest.approx(0.1, abs=0.05)  # 0.1 ± 0.05

class TestCSVDataset:
    def setup_method(self):
        """Called before every method."""
        self.X = np.array([[4, 2, 3, 0], [2, 4, 3, 3], [2, 3, 0, 0]])
        self.y = np.array([[1], [2], [0]])
        self.batch_size = 1
        self.dataset = data.CSVDataset(X=self.X, y=self.y)

    def teardown_method(self):
        """Called after every method."""
        del self.dataset

    def test_len(self):
        assert len(self.X) == len(self.dataset)

    def test_str(self):
        assert str(self.dataset) == f"<Dataset(N={len(self.dataset)})>"


text:
import numpy as np
import matplotlib.pyplot as plt
from skimage import io
import imageio
from sklearn.cluster import KMeans
from sklearn.svm import SVC

##18----------------------------------------------------------------------------------------------
classes = {'wall': 0, 'polyp': 1, 'dirt': 2, 'dark': 3, 'hole': 4}
n_classes = len(classes)
palette = np.uint8([[234, 164, 128], [234, 112, 78], [234, 207, 79], [0, 0, 0], [67,52,31]])

img = io.imread('./images/18.jpg')
img2 = img
rows, cols, bands = img.shape

X = img.reshape(rows*cols, bands)
kmeans = KMeans(n_clusters=n_classes, random_state=0).fit(X)

supervised = n_classes*np.ones(shape=(rows, cols), dtype=np.int)
supervised[205:227, 90:200] = classes['wall']
supervised[165:205, 125:210] = classes['wall']
supervised[110:195, 0:82] = classes['polyp']
supervised[190:205, 12:35] = classes['polyp']
supervised[105:149, 50:120] = classes['polyp']
supervised[50:140, 130:185] = classes['dirt']
supervised[40:75, 15:150] = classes['dirt']
#supervised[40:45, 50:55] = classes['dirt']
#supervised[50:60, 15:25] = classes['hole']
supervised[0:12, 0:12] = classes['dark']

y = supervised.ravel()
train = np.flatnonzero(supervised < n_classes)
test = np.flatnonzero(supervised == n_classes)

clf = SVC()
clf.fit(X[train], y[train])
y[test] = clf.predict(X[test])
supervised = y.reshape(rows, cols)
#io.imshow(palette[supervised])
plt.imshow(img2)
plt.show()
plt.imshow(palette[supervised], cmap=plt.cm.binary)
plt.show()

supervised_18 = supervised
imageio.imwrite('./images/18_SVM.jpg',palette[supervised])

max_stars_repo_name: d-kozak/sdu-data-science

max_issues_repo_name: d-kozak/sdu-data-science

max_forks_repo_name: d-kozak/sdu-data-science


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

import numpy as np
from scipy.linalg import norm

from vaspy.iter import XdatCar, AniFile

fixed_coord = [0.00000000, 1.92610195, 10.69285799]
O_index = 1
max_x = 0.5
max_y = 0.5

def pdc_coord(coord):
    ''' Process the coordinate using PDC.
    '''
    new_coord = [0.0]*3
    for i, c in enumerate(coord):
        if (i == 0 and c > max_x) or (i == 1 and c > max_y):
            c -= 1.0
        new_coord[i] = c

    return new_coord

if '__main__' == __name__:
    ani = AniFile()
    xdatcar = XdatCar()

    fixed_coord = xdatcar.dir2cart(xdatcar.bases, pdc_coord(xdatcar.cart2dir(xdatcar.bases, fixed_coord)))

    atoms = [0]*5

    for xyz in ani:
        coord = xyz.data[O_index]
        # Convert cartesian to direct.
        dir_coord = xdatcar.cart2dir(xdatcar.bases, coord)
        dir_coord = pdc_coord(dir_coord)

        # Convert coordinate back to cartesian coordinate.
        cart_coord = xdatcar.dir2cart(xdatcar.bases, dir_coord)

        # Distance.
        distance = norm(cart_coord - np.array(fixed_coord))

        if 0.0 < distance < 0.2:
            atoms[0] += 1
        elif 0.2 < distance < 0.4:
            atoms[1] += 1
        elif 0.4 < distance < 0.6:
            atoms[2] += 1
        elif 0.6 < distance < 0.8:
            atoms[3] += 1
        elif 0.8 < distance < 1.0:
            atoms[4] += 1
        else:
            print('WARNING: {} -> {}, distance: {}'.format(cart_coord, fixed_coord, distance))

    with open('result.txt', 'w') as f:
        f.write('atoms = {}'.format(atoms))

max_stars_repo_name: faradaymahe/VASPpy

max_issues_repo_name: hsulab/VASPy

max_forks_repo_name: hsulab/VASPy


text:
import copy

import numpy as np

from .base import Visuals
from . import color

from .. import util
from .. import caching
from .. import grouping

from .material import SimpleMaterial, PBRMaterial, empty_material  # NOQA

class TextureVisuals(Visuals):
    def __init__(self,
                 uv=None,
                 material=None,
                 image=None,
                 face_materials=None):
        """
        Store a single material and per-vertex UV coordinates
        for a mesh.

        If passed UV coordinates and a single image it will
        create a SimpleMaterial for the image.

        Parameters
        --------------
        uv : (n, 2) float
          UV coordinates for the mesh
        material : Material
          Store images and properties
        image : PIL.Image
          Can be passed to automatically create material
        """

        # store values we care about enough to hash
        self.vertex_attributes = caching.DataStore()
        # cache calculated values
        self._cache = caching.Cache(self.vertex_attributes.fast_hash)

        # should be (n, 2) float
        self.uv = uv

        if material is None:
            if image is None:
                self.material = empty_material()
            else:
                # if an image is passed create a SimpleMaterial
                self.material = SimpleMaterial(image=image)
        else:
            # if passed assign
            self.material = material

        self.face_materials = face_materials

    def _verify_crc(self):
        """
        Dump the cache if anything in self.vertex_attributes has changed.
        """
        self._cache.verify()

    @property
    def kind(self):
        """
        Return the type of visual data stored

        Returns
        ----------
        kind : str
          What type of visuals are defined
        """
        return 'texture'

    @property
    def defined(self):
        """
        Check if any data is stored

        Returns
        ----------
        defined : bool
          Are UV coordinates and images set?
        """
        ok = self.material is not None
        return ok

    def crc(self):
        """
        Get a CRC of the stored data.

        Returns
        --------------
        crc : int
          Hash of items in self.vertex_attributes
        """
        return self.vertex_attributes.crc()

    @property
    def uv(self):
        """
        Get the stored UV coordinates.

        Returns
        ------------
        uv : (n, 2) float
          Pixel position per- vertex
        """
        return self.vertex_attributes.get('uv', None)

    @uv.setter
    def uv(self, values):
        """
        Set the UV coordinates.

        Parameters
        --------------
        values : (n, 2) float
          Pixel locations on a texture per- vertex
        """
        if values is None:
            self.vertex_attributes['uv'] = None
        else:
            self.vertex_attributes['uv'] = np.asanyarray(
                values, dtype=np.float64)

    def copy(self, uv=None):
        """
        Return a copy of the current TextureVisuals object.

        Returns
        ----------
        copied : TextureVisuals
          Contains the same information in a new object
        """
        if uv is None:
            uv = self.uv
        if uv is not None:
            uv = uv.copy()
        copied = TextureVisuals(
            uv=uv,
            material=self.material.copy(),
            face_materials=copy.copy(self.face_materials))

        return copied

    def to_color(self):
        """
        Convert textured visuals to a ColorVisuals with vertex
        color calculated from texture.

        Returns
        -----------
        vis : trimesh.visuals.ColorVisuals
          Contains vertex color from texture
        """
        # find the color at each UV coordinate
        colors = self.material.to_color(self.uv)
        # create ColorVisuals from result
        vis = color.ColorVisuals(vertex_colors=colors)
        return vis

    def face_subset(self, face_index):
        """
        Get a copy of
        """
        if self.uv is not None:
            indices = np.unique(self.mesh.faces[face_index].flatten())
            return self.copy(self.uv[indices])
        else:
            return self.copy()

    def update_vertices(self, mask):
        """
        Apply a mask to remove or duplicate vertex properties.

        Parameters
        ------------
        mask : (len(vertices),) bool or (n,) int
          Mask which can be used like: `vertex_attribute[mask]`
        """
        # collect updated masked values
        updates = {}
        for key, value in self.vertex_attributes.items():
            # DataStore will convert None to zero-length array
            if len(value) == 0:
                continue
            try:
                # store the update
                updates[key] = value[mask]
            except BaseException:
                # usual reason is an incorrect size or index
                util.log.warning('failed to update visual: `{}`'.format(key))
        # clear all values from the vertex attributes
        self.vertex_attributes.clear()
        # apply the updated values
        self.vertex_attributes.update(updates)

    def update_faces(self, mask):
        """
        Apply a mask to remove or duplicate face properties,
        not applicable to texture visuals.
        """
        pass

    def concatenate(self, others):
        """
        Concatenate this TextureVisuals object with others
        and return the result without modifying this visual.

        Parameters
        -----------
        others : (n,) Visuals
          Other visual objects to concatenate

        Returns
        -----------
        concatenated : TextureVisuals
          Concatenated visual objects
        """
        util.log.warning('concatenating texture: may result in visual artifacts')
        from .objects import concatenate
        return concatenate(self, others)

def unmerge_faces(faces, *args, **kwargs):
    """
    Textured meshes can come with faces referencing vertex
    indices (`v`) and an array the same shape which references
    vertex texture indices (`vt`) and sometimes even normal (`vn`).

    Vertex locations with different values of any of these can't
    be considered the "same" vertex, and for our simple data
    model we need to not combine these vertices.

    Parameters
    -------------
    faces : (n, d) int
      References vertex indices
    *args : (n, d) int
      Various references of corresponding values
      This is usually UV coordinates or normal indexes
    maintain_faces : bool
      Do not alter original faces and return no-op masks.

    Returns
    -------------
    new_faces : (m, d) int
      New faces for masked vertices
    mask_v : (p,) int
      A mask to apply to vertices
    mask_* : (p,) int
      A mask to apply to vt array to get matching UV coordinates
      Returns as many of these as args were passed
    """
    # unfortunately Python2 doesn't let us put named kwargs
    # after an `*args` sequence so we have to do this ugly get
    maintain_faces = kwargs.get('maintain_faces', False)

    # don't alter faces
    if maintain_faces:
        # start with not altering faces at all
        result = [faces]
        # find the maximum index referenced by faces
        max_idx = faces.max()
        # add a vertex mask which is just ordered
        result.append(np.arange(max_idx + 1))

        # now given the order is fixed do our best on the rest of the order
        for arg in args:
            # create a mask of the attribute-vertex mapping
            # note that these might conflict since we're not unmerging
            masks = np.zeros((3, max_idx + 1), dtype=np.int64)
            # set the mask using the unmodified face indexes
            for i, f, a in zip(range(3), faces.T, arg.T):
                masks[i][f] = a
            # find the most commonly occurring attribute (i.e. UV coordinate)
            # and use that index note that this is doing a float conversion
            # and then median before converting back to int: could also do this as
            # a column diff and sort but this seemed easier and is fast enough
            result.append(np.median(masks, axis=0).astype(np.int64))

        return result

    # stack into pairs of (vertex index, texture index)
    stackable = [np.asanyarray(faces).reshape(-1)]
    # append multiple args to the correlated stack
    # this is usually UV coordinates (vt) and normals (vn)
    for arg in args:
        stackable.append(np.asanyarray(arg).reshape(-1))

    # unify them into rows of a numpy array
    stack = np.column_stack(stackable)
    # find unique pairs: we're trying to avoid merging
    # vertices that have the same position but different
    # texture coordinates
    unique, inverse = grouping.unique_rows(stack)

    # only take the unique pairs
    pairs = stack[unique]
    # try to maintain original vertex order
    order = pairs[:, 0].argsort()
    # apply the order to the pairs
    pairs = pairs[order]

    # we re-ordered the vertices to try to maintain
    # the original vertex order as much as possible
    # so to reconstruct the faces we need to remap
    remap = np.zeros(len(order), dtype=np.int64)
    remap[order] = np.arange(len(order))

    # the faces are just the inverse with the new order
    new_faces = remap[inverse].reshape((-1, 3))

    # the mask for vertices and masks for other args
    result = [new_faces]
    result.extend(pairs.T)

    return result

def power_resize(image, resample=1, square=False):
    """
    Resize a PIL image so every dimension is a power of two.

    Parameters
    ------------
    image : PIL.Image
      Input image
    resample : int
      Passed to Image.resize
    square : bool
      If True, upsize to a square image

    Returns
    -------------
    resized : PIL.Image
      Input image resized
    """
    # what is the current resolution of the image in pixels
    size = np.array(image.size, dtype=np.int64)
    # what is the resolution of the image upsized to the nearest
    # power of two on each axis: allow rectangular textures
    new_size = (2 ** np.ceil(np.log2(size))).astype(np.int64)

    # make every dimension the largest
    if square:
        new_size = np.ones(2, dtype=np.int64) * new_size.max()

    # if we're not powers of two upsize
    if (size != new_size).any():
        return image.resize(new_size, resample=resample)

    return image.copy()

max_stars_repo_name: friggog/trimesh

max_issues_repo_name: friggog/trimesh

max_forks_repo_name: friggog/trimesh


text:
#   Copyright (c) 2022 PaddlePaddle Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import numpy as np

class MetricRecord(object):
    """
    One record for a single metric at a given execution step.
    """

    def __init__(self, value, step):
        self._value = value
        self._step = step

    @property
    def value(self):
        return self._value

    @value.setter
    def value(self, value):
        self._value = value

    @property
    def step(self):
        return self._step

    @step.setter
    def step(self, step):
        self._step = step

    def mean(self):
        return np.mean(self.value)

    def get_state(self):
        return {"value": self.value, "step": self.step}

    @classmethod
    def from_state(cls, state):
        return cls(**state)

    def __eq__(self, other):
        if not isinstance(other, MetricRecord):
            return False
        return other.value == self.value and other.step == self.step

    def __repr__(self):
        return "MetricRecord(value={}, step={})".format(self.value, self.step)

class MetricRecords(object):
    """
    Records of a single metric across different executions.
    """

    def __init__(self, direction="min"):
        if direction not in {"min", "max"}:
            raise ValueError(
                "direction should be one of {min, max}, but got: {}.".format(
                    direction))
        self._direction = direction
        self._records = {}

    @property
    def records(self):
        return sorted(self._records.values(), key=lambda r: r.step)

    @records.setter
    def records(self, records):
        for r in records:
            self.update(r.value, step=r.step)

    @property
    def direction(self):
        return self._direction

    @direction.setter
    def direction(self, direction):
        self._direction = direction

    def update(self, value, step=0):
        if step in self._records:
            self._records[step].set_value(value)
        else:
            self._records[step] = MetricRecord(value, step=step)

    def get_best_value(self):
        values = list(r.mean() for r in self._records.values())
        if not values:
            return None
        if self._direction == "min":
            return np.nanmin(values)
        return np.nanmax(values)

    def get_best_step(self):
        best_value = self.get_best_value()
        if best_value is None:
            return None
        for r in self._records.values():
            if r.mean() == best_value:
                return r.step

    def get_statistics(self):
        records = self.records
        records_values = [r.mean() for r in records]
        if not len(records_values):
            return {}
        return {
            "min": float(np.nanmin(records_values)),
            "max": float(np.nanmax(records_values)),
            "mean": float(np.nanmean(records_values)),
            "median": float(np.nanmedian(records_values)),
            "var": float(np.nanvar(records_values)),
            "std": float(np.nanstd(records_values)),
        }

    def get_state(self):
        state = {}
        state["direction"] = self._direction
        state["records"] = [r.get_state() for r in self.records]
        return state

    @classmethod
    def from_state(cls, state):
        records = cls(state["direction"])
        records.records = [MetricRecord.from_state(r) for r in state["records"]]
        print("here 1", records.records)
        return records

class MetricsRecorder(object):
    """
    Record the values for all metrics.
    """

    def __init__(self, metrics=None):
        self._records = {}
        self.register_metrics(metrics)

    @property
    def records(self):
        return self._records

    def exists(self, name):
        return name in self._records

    def register_metrics(self, metrics=None):
        metrics = metrics or []
        for metric in metrics:
            self.register(metric.name)

    def register(self, name, direction=None):
        if self.exists(name):
            raise ValueError("Metric {} have been registered.".format(name))
        if direction is None:
            direction = "min"
        self._records[name] = MetricRecords(direction)

    def update(self, name, value, step=0):
        value = float(value)
        if not self.exists(name):
            self.register(name)

        prev_best = self._records[name].get_best_value()
        self._records[name].update(value, step=step)
        new_best = self._records[name].get_best_value()

        improved = new_best != prev_best
        return improved

    def get_records(self, name):
        return self._records[name].records

    def set_records(self, name, records):
        if not self.exists(name):
            self.register(name)
        self._records[name].records = records

    def get_best_value(self, name):
        return self._records[name].get_best_value()

    def get_best_step(self, name):
        return self._records[name].get_best_step()

    def get_statistics(self, name):
        return self._records[name].get_statistics()

    def get_state(self):
        return {
            "metrics": {
                name: metric_records.get_state()
                for name, metric_records in self._records.items()
            }
        }

    @classmethod
    def from_state(cls, state):
        recorder = cls()
        recorder._records = {
            name: MetricRecords.from_state(metric_records)
            for name, metric_records in state["metrics"].items()
        }
        return recorder

max_stars_repo_name: DevilCarp/Paddle

max_issues_repo_name: DevilCarp/Paddle

max_forks_repo_name: DevilCarp/Paddle


text:
"""
A script to build the command line text for downloading fnet data in the website.
    @problem contain more stations so that "no avaliable data for some selections"
"""
import obspy
import numpy as np
import sys
import pyperclip

def get_download_dict(events):
    result = {}
    for event in events[int(sys.argv[1]): int(sys.argv[1]) + 4]:
        starttime = event.origins[0].time - 2 * 60
        endtime = event.origins[0].time + 40 * 60
        starttime = starttime.datetime
        endtime = endtime.datetime
        name = event.resource_id.id.split("/")[2]
        result[name] = (
            starttime.strftime("%Y/%m/%d,%H:%M:%S"),
            endtime.strftime("%Y/%m/%d,%H:%M:%S"),
        )
    return result

def main():
    station_list = np.loadtxt("../data/fnet_stationlist", dtype=np.str)
    command = ""
    events = obspy.read_events("../data/event_info.xml")
    download_dict = get_download_dict(events)
    for key in download_dict:
        for item in station_list:
            command += (
                f"get {item} BHZ {download_dict[key][0]} {download_dict[key][1]}\n"
            )
            command += (
                f"get {item} BHN {download_dict[key][0]} {download_dict[key][1]}\n"
            )
            command += (
                f"get {item} BHE {download_dict[key][0]} {download_dict[key][1]}\n"
            )
    pyperclip.copy(command)

if __name__ == "__main__":
    main()

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts