text:
import numpy as np
from math import *
def get_parameters():
    f=open('Base.txt','r')
    l=f.readline()
    base=[]
    while len(l)>0:
        l=l.replace('\n','')
        l=l.split(',')
        base.append(l)
        l=f.readline()
    f.close()
    return base

def get_ligand(file):
    f=open(file,'r')
    atom=f.readline()
    atoms=[]
    while len(atom)>0:
        if atom[0:4]=='ATOM':
            serial=atom[6:11].replace(' ','')
            name=atom[11:16].replace(' ','')
            x=atom[31:38].replace(' ','')
            y=atom[38:46].replace(' ','')
            z=atom[46:54].replace(' ','')
            o=atom[55:61].replace(' ','')
            t=atom[61:67].replace(' ','')
            p=atom[67:77].replace(' ','')
            at=atom[77:79].replace(' ','')
            atoms.append([serial,name,float(x),float(y),float(z),o,t,p,at])
        atom=f.readline()
    f.close()
    return atoms

def get_root(file):
    f=open(file,'r')
    atom=f.readline()
    atoms=[]
    while (atom[0:4]!='ROOT')&(len(atom)>0):
        atom=f.readline()
    atom=f.readline()
    while (atom[0:7]!='ENDROOT')&(len(atom)>0):
        if atom[0:4]=='ATOM':
            serial=atom[6:11].replace(' ','')
            name=atom[11:16].replace(' ','')
            x=atom[31:39].replace(' ','')
            y=atom[39:47].replace(' ','')
            z=atom[47:55].replace(' ','')
            o=atom[55:61].replace(' ','')
            t=atom[61:67].replace(' ','')
            p=atom[67:77].replace(' ','')
            at=atom[77:79].replace(' ','')
            atoms.append([serial,name,x,y,z,o,t,p,at])
        atom=f.readline()
    f.close()
    return atoms

def get_branch(file):
    f=open(file,'r')
    line=f.readline()
    branches=[]
    while len(line)>0:
        if (line[0:6]=='BRANCH')&(len(line)>0):
            start=line[6:10].replace(' ','')
            end=line[10:14].replace(' ','')
            branches.append([start,end])
        line=f.readline()
    f.close()
    branch_list=[]
    for i in range(0,len(branches)):
        f=open(file,'r')
        lst=[]
        while 1:
            line=f.readline()
            if line[0:4]=='ATOM':
                start=line[6:11].replace(' ','')
                if (branches[i][0]==start):
                    while 1:
                        if line[0:4]=='ATOM':
                            serial=line[6:11].replace(' ','')
                            lst.append(serial)
                        end=branches[i][1]
                        line=f.readline()
                        if line[0:9]=='ENDBRANCH':
                            s=line[9:13].replace(' ','')
                            e=line[13:17].replace(' ','')
                            if (s==start)&(e==end):
                                break
                    branch_list.append(lst)
                    f.close()
                    break
    return branches,branch_list

def select_atoms(R,L):
    from math import sqrt
    SA=[]
    for i in range(0,len(R)):
        x1=float(R[i][2])
        y1=float(R[i][3])
        z1=float(R[i][4])
        for j in range (0,len(L)):
            x=float(L[j][2])
            y=float(L[j][3])
            z=float(L[j][4])
            dis=sqrt(pow(x-x1,2)+pow(y-y1,2)+pow(z-z1,2))
            if dis<4.72:
                SA.append(R[i])
                break
    return SA

def select_closet(L,P):
    from math import sqrt
    l=len(L)
    x=float(P[0])
    y=float(P[1])
    z=float(P[2])
    a=float(L[0][2])
    b=float(L[0][3])
    c=float(L[0][4])
    dis=sqrt(pow(x-a,2)+pow(y-b,2)+pow(z-c,2))
    j=0
    for i in range(1,len(L)):
        a=float(L[i][2])
        b=float(L[i][3])
        c=float(L[i][4])
        d=sqrt(pow(x-a,2)+pow(y-b,2)+pow(z-c,2))
        if d<dis:
            dis=d
            j=i
    return j

def set_ligand(L):
    x=np.random.rand()
    y=np.random.rand()
    z=np.random.rand()
    for i in range(0,len(L)):
        L[i][2]=L[i][2]+x
        L[i][3]=L[i][3]+y
        L[i][4]=L[i][4]+z
    return L

def RMSD(Ligand,L):
    r=0
    for i in range(0,len(Ligand)):
        x=pow(L[i][2]-Ligand[i][2],2)
        y=pow(L[i][3]-Ligand[i][3],2)
        z=pow(L[i][4]-Ligand[i][4],2)
        r=x+y+z
    r=r/len(L)
    return sqrt(r)

max_stars_repo_name: MasoudiYosef/ProteinPeptideDocking

max_issues_repo_name: MasoudiYosef/ProteinPeptideDocking

max_forks_repo_name: MasoudiYosef/ProteinPeptideDocking


text:
# IMPORTS --------------------------------------------------------------

# Standard packages
import unittest
import pdb

# Third party packages
import pandas as pd
import numpy as np

# Local packages
from ccat import config
from ccat.controller.helper import time as t

from ccat.model.database.bucket import Bucket
from ccat.controller.signal.overtraded import Overtraded
import ccat.controller.helper.df_magic as df_magic

# TESTS ----------------------------------------------------------------

class Test_controller_signal_overtraded_real_data(unittest.TestCase):
    '''Test if the overtraded module returns the expexted results by
    exposing it to a real dataset
    '''

    @classmethod
    def setUpClass(cls):

        cls.len_rsi = 20
        cls.high = 60
        cls.low = 40
        cls.col = 'price_close'
        cls.prefix = 'overtraded'

        # Get a bucket object from Bucket
        cls.bucket = Bucket(
            market_id = 1,
            timeframe_id = 1)

        # Get a dataframe with all the data for the market and timeframe
        cls.df_bucket = cls.bucket.read_all()

        # Get an instance of the Overtraded object
        cls.overtraded = Overtraded(
            df_bucket = cls.df_bucket,
            len_rsi = cls.len_rsi,
            high = cls.high,
            low = cls.low ,
            col = cls.col,
            prefix = cls.prefix)

    def test_overtraded_get_current_signal_only(self):
        signal = self.overtraded.signal()
        self.assertIn(signal, [0, 1, -1])

    def test_overtraded_get_custom_rows(self):
        df = self.overtraded.get(rows = 1)
        self.assertEqual(len(df), 1)

    def test_overtraded_get_custom_columns(self):

        cols = [
            'id',
            # f'{self.prefix}_long',
            f'{self.prefix}_short',
            f'{self.prefix}_signal']

        df = self.overtraded.get(cols = cols)

        self.assertCountEqual(cols, df.columns)
        self.assertEqual(len(df.columns), len(cols))
        self.assertEqual(len(df.head(10)), 10)

    def test_overtraded_get_default(self):

        df = self.overtraded.get()

        self.assertEqual(len(df.columns), 4)
        self.assertEqual(len(df.head(10)), 10)

    def test_overtraded_short(self):
        df = self.overtraded.short()
        self.assertEqual(len(df.columns), 2)
        self.assertEqual(len(df.head(10)), 10)

    def test_overtraded_long(self):
        df = self.overtraded.long()
        self.assertEqual(len(df.columns), 2)
        self.assertEqual(len(df.head(10)), 10)

    def test_overtraded_pandas_dataframe(self):
        self.df_overtraded = self.overtraded.get()
        self.assertIsInstance(self.df_overtraded, pd.DataFrame)

    def test_overtraded_arguments_passed_ok(self):

        self.assertCountEqual(
            self.df_bucket.columns,
            self.overtraded.df_bucket.columns)

        self.assertEqual(self.len_rsi, self.overtraded.len_rsi)
        self.assertEqual(self.high, self.overtraded.high)
        self.assertEqual(self.low, self.overtraded.low)
        self.assertEqual(self.col, self.overtraded.col_bucket)
        self.assertEqual(self.prefix, self.overtraded.prefix)

    def test_overtraded_not_none(self):
        self.assertIsNotNone(self.overtraded)
        self.assertIsNotNone(self.overtraded.df_bucket)
        self.assertIsNotNone(self.overtraded.len_rsi)
        self.assertIsNotNone(self.overtraded.high)
        self.assertIsNotNone(self.overtraded.low)
        self.assertIsNotNone(self.overtraded.col_bucket)
        self.assertIsNotNone(self.overtraded.prefix)

class Test_controller_signal_overtraded_fake_data(unittest.TestCase):
    '''Test if the overtraded module returns the expexted results by
    exposing it to a known dataset
    '''

    @classmethod
    def setUpClass(cls):

        # Set params
        cls.len_rsi = 20
        cls.high = 60
        cls.low = 40
        cls.col = 'price_close'
        cls.prefix = 'overtraded'

        # Generate fake data
        x = np.arange(0, 10, 0.1)
        y = (np.sin(x) + 1)/2 *1000

        # Create a dataframe with the fake data
        cls.df_in = pd.DataFrame({'price_close':y})

        # Set id column = index and then as the id column for merging
        cls.df_in['id'] = cls.df_in.index
        cls.df_in = cls.df_in[['id', 'price_close']]

        # Get an instance of the Overtraded object
        cls.overtraded = Overtraded(
            df_bucket = cls.df_in,
            len_rsi = cls.len_rsi,
            high = cls.high,
            low = cls.low,
            col = cls.col,
            prefix = cls.prefix)

    def test_overtraded_get_custom_columns(self):

        cols = [
            'id',
            f'{self.prefix}_long',
            # f'{self.prefix}_short',
            f'{self.prefix}_signal']

        df = self.overtraded.get(cols = cols)
        self.assertCountEqual(cols, df.columns)
        self.assertEqual(len(df.columns), len(cols))
        self.assertEqual(len(df.head(10)), 10)

    def test_overtraded_get_default(self):
        df = self.overtraded.get()
        self.assertEqual(len(df.columns), 4)
        self.assertEqual(len(df.head(10)), 10)

    def test_overtraded_short(self):
        df = self.overtraded.short()
        self.assertEqual(len(df.columns), 2)
        self.assertEqual(len(df.head(10)), 10)

    def test_overtraded_long(self):
        df = self.overtraded.long()
        self.assertEqual(len(df.columns), 2)
        self.assertEqual(len(df.head(10)), 10)

# MAIN -----------------------------------------------------------------

if __name__ == '__main__':

    unittest.main()

ext: py

max_stars_repo_name: bhaveshrawal/Python

max_issues_repo_name: bliiir/ccat

max_forks_repo_name: bliiir/ccat

reason: import numpy


text:
# -*- coding: utf8 -*-
"""
======================================
    Project Name: NLP
    File Name: kg_globalpointer
    Author: czh
    Create Date: 2022/2/10
--------------------------------------
    Change Activity:
======================================
"""
import os
import sys
import json
import codecs
import logging
from tqdm import tqdm

import numpy as np
import torch
from torch.optim import AdamW, Adam, SGD
from torch.utils.data import DataLoader
from transformers import BertConfig, BertTokenizerFast, get_scheduler
sys.path.append("/data/chenzhihao/NLP")

from nlp.models.bert_for_relation_extraction import GlobalPointerForRel
from nlp.callback.optimizers.ema import EMA
from nlp.processors.dataset import MyDataset
from nlp.losses.loss import global_pointer_crossentropy
from nlp.utils.tokenizers import tokenize
from nlp.tools.common import init_wandb_writer

LOGGER = logging.getLogger(__name__)
root = os.path.abspath(os.path.dirname(__file__))

config = {
    'batch_size': 10,
    'max_seq_length': 300,
    'max_epochs': 100,
    'lr_rate': 1e-5,
    'gradient_accumulation_steps': 1,
    'logging_steps': 500,
    'warmup_ratio': 0.1
}
# bert_model_name_or_path = "/Users/czh/Downloads/chinese-roberta-ext"
# data_dir = "/Users/czh/PycharmProjects/NLP/datas/kg"
bert_model_name_or_path = "/data/chenzhihao/chinese-roberta-ext"
data_dir = "/data/chenzhihao/NLP/datas/kg"
train_data_dir = os.path.join(data_dir, 'train_data.json')
dev_data_dir = os.path.join(data_dir, 'dev_data.json')
scheme_data_dir = os.path.join(data_dir, 'all_50_schemes.json')
output_dir = "../output_file_dir/kg"
num_worker = 0
device = torch.device('cuda:3') if torch.cuda.is_available() else torch.device('cpu')

if not os.path.exists(output_dir):
    os.makedirs(output_dir)

WANDB, run = init_wandb_writer(project_name='kg_globalpointer',
                               train_args=config,
                               group_name="NLP",
                               experiment_name="test")

def load_data(filename):
    """加载数据
    单条格式：{'text': text, 'spo_list': [(s, p, o)]}
    """
    samples = []
    with codecs.open(filename, encoding='utf-8') as f:
        for line in tqdm(f, desc=f"load datas from {filename}"):
            line = line.strip()
            if not line:
                continue
            try:
                item = json.loads(line)
            except Exception as e:
                print("error line: ", line)
                raise e
            samples.append({
                'text': item['text'],
                'spo_list': [(spo['subject'], spo['predicate'], spo['object'])
                             for spo in item['spo_list']]
            })
    return samples

def load_schemes(filename):

    with codecs.open(filename) as f:
        alist = json.load(f)

    id2predicate, predicate2id = alist
    id2predicate = {int(k): v for k, v in id2predicate.items()}
    predicate2id = {k: int(v) for k, v in predicate2id.items()}
    return predicate2id, id2predicate

rel2id, id2rel = load_schemes(scheme_data_dir)
print("number of relations: ", len(rel2id))
print("rel2id: ", rel2id)
print("id2rel: ", id2rel)

tokenizer = BertTokenizerFast.from_pretrained(bert_model_name_or_path, do_lower_case=False, add_special_tokens=True)
bert_config = BertConfig.from_pretrained(bert_model_name_or_path)
my_model = GlobalPointerForRel(config=bert_config, entity_types_num=2,
                               encoder_model_path=bert_model_name_or_path,
                               relation_num=len(rel2id), efficient=False)
my_model.to(device)
no_decay = ["bias", "LayerNorm.weight", "LayerNorm.bias"]
bert_param_optimizer = list(my_model.encoder.named_parameters())
optimizer_grouped_parameters = [
    {'params': [p for n, p in bert_param_optimizer if not any(nd in n for nd in no_decay)],
     'weight_decay': 0.01, 'lr': config['lr_rate']},
    {'params': [p for n, p in bert_param_optimizer if any(nd in n for nd in no_decay)], 'weight_decay': 0.0,
     'lr': config['lr_rate']}
]

optimizer = AdamW(params=optimizer_grouped_parameters, lr=config['lr_rate'])
ema = EMA(my_model, 0.999)
ema.register()

class SPO(tuple):
    """用来存三元组的类
    表现跟tuple基本一致，只是重写了 __hash__ 和 __eq__ 方法，
    使得在判断两个三元组是否等价时容错性更好。
    """
    def __init__(self, spo):
        self.spox = (
            tuple(tokenizer.tokenize(spo[0], add_special_tokens=False)),
            spo[1],
            tuple(tokenizer.tokenize(spo[2], add_special_tokens=False)),
        )

    def __hash__(self):
        return self.spox.__hash__()

    def __eq__(self, spo):
        return self.spox == spo.spox

def search(pattern, sequence):
    """从sequence中寻找子串pattern
    如果找到，返回第一个下标；否则返回-1。
    """
    n = len(pattern)
    for i in range(len(sequence)):
        if sequence[i:i + n] == pattern:
            return i
    return -1

def generate_inputs(datasets, max_seq_len):
    all_inputs = []
    error = 0
    for sample_ in tqdm(datasets):
        # tokens = tokenizer.tokenize(sample_['text'])
        # if len(tokens) > max_seq_len-2:
        #     tokens = tokens[:max_seq_len]
        # tokens = ["[CLS]"] + tokens + ["[SEP]"]
        # input_ids = tokenizer.convert_tokens_to_ids(tokens)
        # attention_mask = [1] * len(input_ids)
        # token_type_ids = [0] * len(input_ids)
        outputs = tokenizer.encode_plus(text=sample_['text'], add_special_tokens=True,
                                        max_length=max_seq_len, return_offsets_mapping=True)
        input_ids = outputs['input_ids']
        attention_mask = outputs['attention_mask']
        token_type_ids = outputs['token_type_ids']
        offset_mapping = outputs['offset_mapping']

        spoes = set()
        for s, p, o in sample_['spo_list']:
            s_ids = tokenizer.encode_plus(s, add_special_tokens=False)['input_ids']
            p = rel2id[p]
            o_ids = tokenizer.encode_plus(o, add_special_tokens=False)['input_ids']
            sh = search(s_ids, input_ids)
            oh = search(o_ids, input_ids)
            if sh != -1 and oh != -1:
                spoes.add((sh, sh + len(s_ids) - 1, p, oh, oh + len(o_ids) - 1))
            else:
                print(sample_)
                error += 1
        item = {
            'input_ids': input_ids,
            'attention_mask': attention_mask,
            'token_type_ids': token_type_ids,
            'labels': list(spoes),
            'text': sample_['text'],
            'spo_list': sample_['spo_list'],
            'offset_mapping': offset_mapping
        }
        all_inputs.append(item)
    print(f"总计 {error} 样本无法对齐")
    return all_inputs

def collate_fn(batch):
    bs = len(batch)
    seq_lens = [len(sample['input_ids']) for sample in batch]
    max_len = max(seq_lens)
    label_lens = [len(sample['labels']) for sample in batch]
    max_spo_num = max(label_lens)

    batch_input_ids = torch.zeros(bs, max_len, dtype=torch.long)
    batch_attention_masks = torch.zeros(bs, max_len, dtype=torch.long)
    batch_token_type_ids = torch.zeros(bs, max_len, dtype=torch.long)
    batch_entity_labels = torch.zeros(bs, 2, max_spo_num, 2, dtype=torch.long)
    batch_head_labels = torch.zeros(
        bs, len(rel2id), max_spo_num, 2, dtype=torch.long
    )
    batch_tail_labels = torch.zeros(
        bs, len(rel2id), max_spo_num, 2, dtype=torch.long
    )
    batch_texts = []
    batch_spo_list = []
    batch_offset_mapping = []

    for i, sample in enumerate(batch):
        batch_input_ids[i, :seq_lens[i]] = torch.tensor(sample['input_ids'])
        batch_attention_masks[i, :seq_lens[i]] = torch.tensor(sample['attention_mask'])
        batch_token_type_ids[i, :seq_lens[i]] = torch.tensor(sample['token_type_ids'])
        for j, (sh, st, p, oh, ot) in enumerate(sample['labels']):
            batch_entity_labels[i, 0, j, :] = torch.tensor([sh, st])
            batch_entity_labels[i, 1, j, :] = torch.tensor([oh, ot])
            batch_head_labels[i, p, j, :] = torch.tensor([sh, oh])
            batch_tail_labels[i, p, j, :] = torch.tensor([st, ot])
        batch_texts.append(sample['text'])
        batch_offset_mapping.append(sample['offset_mapping'])
        batch_spo_list.append(sample['spo_list'])

    item = {'input_ids': batch_input_ids.to(device),
            'attention_mask': batch_attention_masks.to(device),
            'token_type_ids': batch_token_type_ids.to(device),
            'entity_labels': batch_entity_labels.to(device),
            'head_labels': batch_head_labels.to(device),
            'tail_labels': batch_tail_labels.to(device),
            'texts': batch_texts,
            'spo_list': batch_spo_list,
            'offset_mapping': batch_offset_mapping
            }
    return item

# 加载数据集
def data_generator(data_type="train"):
    """
    读取数据，生成DataLoader。
    """
    if data_type == "train":
        train_data = load_data(train_data_dir)
        valid_data = load_data(dev_data_dir)
    elif data_type == "valid":
        valid_data = load_data(dev_data_dir)
        train_data = []

    all_data = train_data + valid_data

    # TODO:句子截取
    max_tok_num = 0
    for sample in all_data:
        tokens = tokenizer.tokenize(sample["text"], add_special_tokens=True)
        max_tok_num = max(max_tok_num, len(tokens))
    assert max_tok_num <= config['max_seq_length'], f"数据文本最大token数量{max_tok_num}超过预设{config['max_seq_length']}"
    max_seq_len = min(max_tok_num, config['max_seq_length'])

    if data_type == "train":
        train_inputs = generate_inputs(train_data, max_seq_len)
        valid_inputs = generate_inputs(valid_data, max_seq_len)
        train_dataloader = DataLoader(MyDataset(train_inputs),
                                      batch_size=config['batch_size'],
                                      shuffle=True,
                                      num_workers=num_worker,
                                      drop_last=False,
                                      collate_fn=collate_fn
                                      )
        valid_dataloader = DataLoader(MyDataset(valid_inputs),
                                      batch_size=config['batch_size'],
                                      shuffle=False,
                                      num_workers=num_worker,
                                      drop_last=False,
                                      collate_fn=collate_fn
                                      )
        return train_dataloader, valid_dataloader
    elif data_type == "valid":
        valid_inputs = generate_inputs(valid_data, max_seq_len)
        valid_dataloader = DataLoader(MyDataset(valid_inputs),
                                      batch_size=config['batch_size'],
                                      shuffle=False,
                                      num_workers=num_worker,
                                      drop_last=False,
                                      )
        return valid_dataloader

def extract_spoes(model, eval_dataloader, threshold=0):
    """抽取输入text所包含的三元组
    """
    # output = tokenizer.encode_plus(text, return_offsets_mapping=True, add_special_tokens=True)
    # token_ids = torch.LongTensor([output['input_ids']]).to(device)
    # segment_ids = torch.LongTensor([output['token_type_ids']]).to(device)
    # attention_mask = torch.LongTensor([output['attention_mask']]).to(device)
    # offset_mapping = output['offset_mapping']

    all_predicts = []
    all_texts = []
    all_spo_list = []

    model.eval()
    for batch in tqdm(eval_dataloader, desc="Evaluating"):

        with torch.no_grad():
            outputs = model(batch['input_ids'], batch['attention_mask'], batch['token_type_ids'])
            outputs = [o.detach().cpu().numpy() for _, o in outputs.items()]
        offset_mappings = batch['offset_mapping']
        texts = batch['texts']
        all_texts.extend(texts)
        all_spo_list.extend(batch['spo_list'])
        assert len(outputs[0]) == len(texts)

        for entity_output, head_output, tai_output, offset_mapping, text in zip(
            outputs[0], outputs[1], outputs[2], offset_mappings, texts
        ):

            # 抽取subject和object
            subjects, objects = set(), set()
            entity_output[:, [0, -1]] -= np.inf
            entity_output[:, :, [0, -1]] -= np.inf
            for l_, h, t in zip(*np.where(entity_output > threshold)):
                if l_ == 0:
                    subjects.add((h, t))
                else:
                    objects.add((h, t))
            # 识别对应的predicate
            spoes = set()
            for sh, st in subjects:
                for oh, ot in objects:
                    p1s = np.where(head_output[:, sh, oh] > threshold)[0]
                    p2s = np.where(tai_output[:, st, ot] > threshold)[0]
                    ps = set(p1s) & set(p2s)
                    for p in ps:
                        try:
                            spoes.add((
                                text[offset_mapping[sh][0]:offset_mapping[st][-1]],
                                id2rel[p],
                                text[offset_mapping[oh][0]:offset_mapping[ot][-1]]
                            ))
                        except Exception as e:
                            print(len(id2rel), len(offset_mapping), offset_mapping)
                            print((sh, st), p, (oh, ot))
                            raise e
            all_predicts.append(list(spoes))
    assert len(all_predicts) == len(all_texts), f"{len(all_predicts)}\t{len(all_texts)}"
    assert len(all_predicts) == len(all_spo_list), f"{len(all_predicts)}\t{len(all_spo_list)}"
    return all_predicts, all_texts, all_spo_list

def evaluate(model, eval_dataloader):
    """评估函数，计算f1、precision、recall
    """
    all_predicts, all_texts, all_spo_list = extract_spoes(model, eval_dataloader)

    x, y, z = 0, 0, 0
    f = open('dev_pred.json', 'w', encoding='utf-8')
    pbar = tqdm()
    f1, precision, recall = 1e-10, 1e-10, 1e-10
    for pred, text, golds in zip(all_predicts, all_texts, all_spo_list):
        r = set(pred)
        t = set(golds)
        x += len(r & t)
        y += len(r)
        z += len(t)
        # pbar.update()
        # pbar.set_description(
        #     'f1: %.5f, precision: %.5f, recall: %.5f' % (f1, precision, recall)
        # )
        s = json.dumps({
            'text': text,
            'spo_list': list(t),
            'spo_list_pred': list(r),
            'new': list(r - t),
            'lack': list(t - r),
        },
            ensure_ascii=False,
            indent=4)
        f.write(s + '\n')
    f1 = 2 * x / (y + z) if y and z else 0
    precision = x / y if y else 0
    recall = x / z if z else 0
    # pbar.close()
    f.close()
    return f1, precision, recall

def train(model, train_dataloader, eval_data, scheduler):
    WANDB.watch(model, log="all")

    model.train()
    model.zero_grad()
    best_f1 = 0.0
    best_epoch = 0
    global_steps = 0
    print("Training")
    for epoch in range(config['max_epochs']):
        tr_loss = 0.0
        epoch_steps = 0
        print("**********************training**********************\n")
        pbar = tqdm()
        for step, batch in enumerate(train_dataloader):
            outputs = model(batch['input_ids'], batch['attention_mask'], batch['token_type_ids'])
            entity_logits = outputs['entity_logits']
            head_logits = outputs['head_logits']
            tail_logits = outputs['tail_logits']

            entity_loss = global_pointer_crossentropy(entity_logits, batch['entity_labels'], sparse=True, mask_zero=True)
            head_loss = global_pointer_crossentropy(head_logits, batch['head_labels'], sparse=True, mask_zero=True)
            tail_loss = global_pointer_crossentropy(tail_logits, batch['tail_labels'], sparse=True, mask_zero=True)

            loss = (entity_loss+head_loss+tail_loss) / 3 / config['gradient_accumulation_steps']

            pbar.update()
            pbar.set_description(
                f"Training at epoch {epoch}\ttotal loss: {loss.item()}\tentity loss: {entity_loss.item()}\t"
                f"head loss: {head_loss.item()}\ttail loss: {tail_loss.item()}"
            )

            loss.backward()
            WANDB.log({'Train/total_loss': loss.item(),
                       'Train/entity_loss': entity_loss.item(),
                       'Train/head_loss': head_loss.item(),
                       'Train/tail_loss': tail_loss.item()}, step=global_steps)
            tr_loss += loss.item()
            epoch_steps += 1
            if (step + 1) % config['gradient_accumulation_steps'] == 0:
                optimizer.step()
                # ema.update()
                # scheduler.step()
                optimizer.zero_grad()
                global_steps += 1

        print("***********************evaluating***********************\n")
        f1, precision, recall = evaluate(model, eval_data)
        print(
            f"Evaluating epoch: {epoch}\ttotal loss: {tr_loss / global_steps}\tf1: {f1}\tprecision: "
            f"{precision}\trecall: {recall}\n")
        WANDB.log({'Eval/f1': f1,
                   'Eval/precision': precision, 'Eval/recall': recall}, step=global_steps)
        if f1 > best_f1:
            best_epoch = epoch
            best_f1 = f1
            model_to_save = (
                model.module if hasattr(model, "module") else model
            )  # Take care of distributed/parallel training
            # model_to_save.save_pretrained(output_dir)
            torch.save(model_to_save.state_dict(), output_dir + '/pytorch.bin')
            tokenizer.save_vocabulary(output_dir)
            print(f"save model to {output_dir}")
        pbar.close()
    return best_epoch, best_f1

def main():
    train_dataloader, eval_dataloader = data_generator()
    t_total = len(train_dataloader) // config['gradient_accumulation_steps'] * config['max_epochs']
    scheduler = get_scheduler(name='linear', optimizer=optimizer, num_warmup_steps=0, num_training_steps=t_total)

    best_epoch, best_f1 = train(my_model, train_dataloader, eval_dataloader, scheduler)
    print(f"best epoch: {best_epoch}\tbest f1: {best_f1}")
    print("Evaluating")
    eval_model = GlobalPointerForRel(config=bert_config,
                                     entity_types_num=2, relation_num=len(rel2id))
    eval_model.load_state_dict(torch.load(output_dir+'/pytorch.bin', map_location='cpu'))
    eval_model.to(device)
    f1, precision, recall = evaluate(eval_model, eval_dataloader)
    print(f"evaluate result: f1: {f1}\tprecision: {precision}\trecall: {recall}")

if __name__ == "__main__":
    main()

max_stars_repo_name: zhihao-chen/NLP-experiments

max_issues_repo_name: zhihao-chen/NLP-experiments

max_forks_repo_name: zhihao-chen/NLP-experiments


text:
import numpy as np
import tensorflow as tf
import gym
import logz
import scipy.signal
import os
import time
from multiprocessing import Process
import shutil

class MyArgument(object):
    def __init__(self,
                 exp_name='vpg',
                 env_name='CartPole-v1',
                 n_iter=100,
                 gamma=.99,
                 min_batch_size=1000,
                 max_path_length=None,
                 learning_rate=1e-3,
                 reward_to_go=True,
                 render=False,
                 normalize_advantage=True,
                 nn_baseline=True,
                 seed=1,
                 n_layers=1,
                 size=32,
                 debug=False,
                 max_loss=100):
        self.exp_name = exp_name
        self.env_name = env_name
        self.n_iter = n_iter
        self.gamma = gamma
        self.min_batch_size = min_batch_size
        self.max_path_length = max_path_length
        self.learning_rate = learning_rate
        self.reward_to_go = reward_to_go
        self.render = render
        self.debug = debug
        self.max_loss = max_loss

        self.normalize_advantage = normalize_advantage
        self.nn_baseline = nn_baseline
        self.seed = seed
        self.n_layers = n_layers
        self.size = size

        base_dir = '/tmp/pg/%s' % self.env_name

        self.log_dir = 'bl_' if self.nn_baseline else ''
        self.log_dir += 'rtg_' if self.reward_to_go else ''
        self.log_dir += 'norm_' if self.normalize_advantage else ''
        self.log_dir += 'nn_%d_%d_' % (self.n_layers, self.size)
        self.log_dir += 'lr_%6f_' % self.learning_rate
        self.log_dir += 'batch_%d_' % self.min_batch_size
        self.log_dir += 't_%d' % self.max_loss

        self.log_dir = os.path.join(base_dir, self.log_dir)
        self.log_dir = os.path.join(self.log_dir, 'seed%d' % self.seed)

class PgModel(object):
    def __init__(self, env, n_layers, size, learning_rate, nn_baseline, debug):
        self.observation_dim = env.observation_space.shape[0]
        self.ph_observation = tf.placeholder(shape=[None, self.observation_dim], name="Observation", dtype=tf.float32)
        self.ph_advance = tf.placeholder(shape=[None], name='advance', dtype=tf.float32)
        self.nn_baseline = nn_baseline
        self.ph_q_value = tf.placeholder(shape=[None], name='QValue', dtype=tf.float32)

        self.debug = debug

        if self.debug:
            self.ph_mean_reward = tf.placeholder(name="reward", dtype=tf.float32)
            tf.summary.scalar("MeanReward", self.ph_mean_reward)

        self.predict_action = None
        self.critic_opt = None
        self.actor_opt = None
        self.ph_action = None
        self.predict_baseline = None
        self.merged = None
        self.critic_loss = None

    @staticmethod
    def build_mlp(input_placeholder, output_size, scope,
                  n_layers=2, size=64, activation=tf.nn.relu, output_activation=None):
        with tf.variable_scope(scope):
            x = tf.keras.Input(tensor=input_placeholder)
            for i in range(n_layers):
                x = tf.keras.layers.Dense(units=size, activation=activation)(x)
            x = tf.keras.layers.Dense(units=output_size, activation=output_activation)(x)
        return x

    def get_predict_action(self, sess, observation):
        action = sess.run(self.predict_action, feed_dict={self.ph_observation: observation[None]})
        return action[0]

    def update(self, sess, observations, actions, q, normalize_advance, mean_reward, max_loss):
        if self.nn_baseline:
            # Update Cirtic Network

            baseline = sess.run(self.predict_baseline, feed_dict={self.ph_observation: observations})
            updated_baseline = baseline * .9 + q * .1
            sess.run([self.critic_opt], feed_dict={self.ph_observation: observations,
                                                   self.ph_q_value: updated_baseline})
            advance = q - baseline
        else:
            advance = q.copy()

        if normalize_advance:
            advance = (advance - np.mean(advance)) / (np.std(advance) + 1e-8)

        # Update the Actor network

        if self.debug:
            _, summary = sess.run([self.action_opt, self.merged], feed_dict={self.ph_observation: observations,
                                                                             self.ph_action: actions,
                                                                             self.ph_advance: advance,
                                                                             self.ph_q_value: q,
                                                                             self.ph_mean_reward: mean_reward})
        else:
            sess.run(self.action_opt, feed_dict={self.ph_observation: observations,
                                                 self.ph_action: actions,
                                                 self.ph_advance: advance})
            summary = None

        return summary

class PgModelContinuous(PgModel):
    def __init__(self, env, n_layers, size, learning_rate, nn_baseline, debug):
        super().__init__(env, n_layers, size, learning_rate, nn_baseline, debug)
        self.action_dim = env.action_space.shape[0]
        self.ph_action = tf.placeholder(shape=[None, self.action_dim], name="Action", dtype=tf.float32)

        # Define the Actor Model
        with tf.variable_scope("Actor"):
            # N x action dim
            # Output activation ?
            self.action_mean = self.build_mlp(input_placeholder=self.ph_observation,
                                              output_size=self.action_dim,
                                              scope="Mean_%d_%d" % (n_layers, size),
                                              size=size,
                                              n_layers=n_layers,
                                              activation=tf.nn.relu,
                                              output_activation=None)

            # action dim
            # self.action_sigma = tf.get_variable(name='Sigma', shape=[self.action_dim],
            #                                     dtype=tf.float32, trainable=True, initializer=tf.ones_initializer())

            self.action_sigma = self.build_mlp(input_placeholder=self.ph_observation,
                                               output_size=self.action_dim,
                                               scope="Sigma_%d_%d" % (n_layers, size),
                                               size=32,
                                               n_layers=1,
                                               activation=tf.nn.relu,
                                               output_activation=tf.nn.sigmoid)

            tf.summary.histogram('Mean', self.action_mean)
            tf.summary.histogram('Std', self.action_sigma)

            # Broadcast expected here
            # Get N x action dim distributions
            self.normal_dist = tf.distributions.Normal(self.action_mean, (self.action_sigma + 1e-8),
                                                       name="PredictDistribution")

            # Expected N* action dis distributions
            self.predict_action = self.normal_dist.sample(name="PredictAction")
            # self.predict_action = tf.clip_by_value(normal_dist.sample(), env.action_space.low, env.action_space.high,
            #                                      name="PredictAction")

            with tf.name_scope("Loss"):
                self.action_prob = self.normal_dist.log_prob(self.ph_action, name="Prob")
                self.actor_loss = - tf.reduce_mean(self.action_prob * tf.expand_dims(self.ph_advance, -1))
                tf.summary.scalar('Actor Loss', self.actor_loss)
        with tf.name_scope("Opt/"):
            self.action_opt = tf.train.AdamOptimizer(learning_rate).minimize(self.actor_loss)

        # Define the Critic Model
        if nn_baseline:
            with tf.name_scope("Critic"):
                self.predict_baseline_2d = self.build_mlp(input_placeholder=self.ph_observation,
                                                          output_size=1,
                                                          scope="NN_%d_%d" % (n_layers, size),
                                                          n_layers=n_layers,
                                                          size=size,
                                                          activation=tf.nn.relu)
                self.predict_baseline = tf.squeeze(self.predict_baseline_2d, axis=1, name="PredictBaseline")

                with tf.name_scope("Loss"):
                    self.critic_loss = tf.losses.mean_squared_error(self.ph_q_value, self.predict_baseline)
                    tf.summary.scalar('Critic Loss', self.critic_loss)
            with tf.name_scope("Opt/"):
                self.critic_opt = tf.train.AdamOptimizer(learning_rate).minimize(self.critic_loss)

        weights = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES)
        for w in weights:
            tf.summary.histogram(w.name, w)
        self.merged = tf.summary.merge_all()

class PgModelDiscrete(PgModel):
    def __init__(self, env, n_layers, size, learning_rate, nn_baseline, debug):
        super().__init__(env, n_layers, size, learning_rate, nn_baseline, debug)
        self.action_dim = env.action_space.n
        self.ph_action = tf.placeholder(shape=[None], name="Action", dtype=tf.int32)

        # Define the Actor Model
        with tf.name_scope("Actor"):
            self.action_logist = self.build_mlp(input_placeholder=self.ph_observation,
                                                output_size=self.action_dim, scope="NN_%d_%d" % (n_layers, size),
                                                size=size, n_layers=n_layers)
            self.predict_action_2d = tf.multinomial(self.action_logist, 1)
            self.predict_action = tf.squeeze(self.predict_action_2d, axis=1, name="PredictAction")

            self.batch_size = tf.shape(self.ph_observation)[0]
            with tf.name_scope('Loss'):
                indices = tf.stack([tf.range(self.batch_size), self.ph_action], axis=1)
                action_prob = tf.gather_nd(tf.nn.softmax(self.action_logist), indices)
                self.actor_loss = tf.reduce_mean(-tf.log(action_prob) * self.ph_advance)
                tf.summary.scalar('Actor loss', self.actor_loss)

        with tf.name_scope("Opt/"):
            self.action_opt = tf.train.AdamOptimizer(learning_rate).minimize(self.actor_loss)

        # Define the Critic Model
        if nn_baseline:
            with tf.name_scope("Critic"):
                self.predict_baseline_2d = self.build_mlp(input_placeholder=self.ph_observation,
                                                          output_size=1,
                                                          scope="NN_%d_%d" % (n_layers, size),
                                                          n_layers=n_layers,
                                                          size=size,
                                                          activation=tf.nn.relu)
                self.predict_baseline = tf.squeeze(self.predict_baseline_2d, axis=1, name="PredictBaseline")

                with tf.name_scope("Loss"):
                    self.critic_loss = tf.losses.mean_squared_error(self.ph_q_value, self.predict_baseline)
                    tf.summary.scalar('Critic Loss', self.critic_loss)
            with tf.name_scope("Opt/"):
                self.critic_opt = tf.train.AdamOptimizer(learning_rate).minimize(self.critic_loss)

        weights = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES)
        for w in weights:
            tf.summary.histogram(w.name, w)
        self.merged = tf.summary.merge_all()

def discount_reward(paths, gamma, reward_to_go):
    if reward_to_go:
        discounted_reward = []
        for path in paths:
            path_len = len(path['reward'])
            discount_factor = [1 * (gamma ** i) for i in range(path_len)]
            for i in range(path_len):
                discounted_reward.append(
                    np.sum(np.array(path['reward'][i:]) * np.array(discount_factor[:path_len - i])))
    else:
        discounted_reward = []
        for path in paths:
            ret_tau = 0
            discount_factor = 1
            for reward in path['reward']:
                ret_tau += reward * discount_factor
                discount_factor *= gamma
            discounted_reward.extend([ret_tau for i in range(len(path['reward']))])

    q_n = np.array(discounted_reward, dtype=np.float32)
    return q_n

def verify_model(sess, model, env):
    action_dim = env.action_space.shape[0]
    observation_dim = env.observation_space.shape[0]

    observations = np.random.randn(observation_dim)
    print("Observation shape:", observations.shape)
    print(observations)
    actions = model.get_predict_action(sess, observations)
    print("action shape:", actions.shape)
    print(actions)
    assert actions.shape == (action_dim,)

    N = 10
    observations = np.random.randn(N, observation_dim)
    print("Observation shape:", observations.shape)
    print(observations)
    predict_action = sess.run(model.predict_action, feed_dict={model.ph_observation: observations})
    print("Action shape:", predict_action.shape)
    print(predict_action)
    assert predict_action.shape == (N, action_dim)

    action_prob = sess.run(model.action_prob, feed_dict={model.ph_observation: observations,
                                                         model.ph_action: predict_action})
    print("Prob:", action_prob)
    assert action_prob.shape == (N, action_dim)

    loss = sess.run(model.actor_loss, feed_dict={model.ph_observation: observations,
                                                 model.ph_action: predict_action,
                                                 model.ph_advance: np.ones(N)})
    print("Loss:", loss)

    for i in range(20):
        sess.run(model.action_opt, feed_dict={model.ph_observation: observations,
                                              model.ph_action: predict_action,
                                              model.ph_advance: np.ones(N)})

    action_prob2 = sess.run(model.action_prob, feed_dict={model.ph_observation: observations,
                                                          model.ph_action: predict_action})
    print("Prob2:", action_prob2)
    assert action_prob2.shape == (N, action_dim)

    loss2 = sess.run(model.actor_loss, feed_dict={model.ph_observation: observations,
                                                  model.ph_action: predict_action,
                                                  model.ph_advance: np.ones(N)})
    print("Loss2:", loss2)
    assert loss2 < loss
    print(np.sum(action_prob2 > action_prob))
    assert np.sum(action_prob2 > action_prob) == N * action_dim

def train_pg(args):
    if os.path.exists(args.log_dir):
        shutil.rmtree(args.log_dir)
    logz.configure_output_dir(args.log_dir)
    logz.save_params(args.__dict__)

    # Set random seeds
    tf.set_random_seed(args.seed)
    np.random.seed(args.seed)

    # Make the gym environment
    env = gym.make(args.env_name)

    # Is this env continuous, or discrete?
    discrete = isinstance(env.action_space, gym.spaces.Discrete)

    if discrete:
        model = PgModelDiscrete(env, args.n_layers, args.size, args.learning_rate, args.nn_baseline, debug=args.debug)
    else:
        model = PgModelContinuous(env, args.n_layers, args.size, args.learning_rate, args.nn_baseline, debug=args.debug)

    tf_config = tf.ConfigProto(inter_op_parallelism_threads=1, intra_op_parallelism_threads=1)
    sess = tf.Session(config=tf_config)
    sess.__enter__()  # equivalent to `with sess:`
    tf.global_variables_initializer().run()  # pylint: disable=E1101
    writer = tf.summary.FileWriter(args.log_dir, sess.graph)

    # if not discrete:
    #     verify_model(sess, model, env)

    max_path_length = args.max_path_length or env.spec.max_episode_steps

    start = time.time()
    for itr in range(args.n_iter):
        end = time.time()
        cost = end - start
        start = end
        if itr > 0:
            if itr == 1:
                mean_cost = cost
            else:
                mean_cost = .9 * mean_cost + .1 * cost
            print("Time: %1f " % mean_cost, "Togo:%1f min" % ((args.n_iter - itr) * mean_cost / 60))

        print("********** Iteration %i ************" % itr)
        time_steps_this_batch = 0
        paths = []
        while True:
            observation = env.reset()
            obs, acs, rewards = [], [], []

            render_this_episode = (len(paths) == 0 and (itr % 10 == 0) and args.render)
            path_steps = 0
            while True:
                if render_this_episode:
                    env.render()
                    time.sleep(0.0001)

                obs.append(observation)
                action = model.get_predict_action(sess, observation)
                acs.append(action)
                observation, rew, done, _ = env.step(action)
                rewards.append(rew)

                path_steps += 1
                if done or path_steps > max_path_length:
                    break

            path = {"observation": np.array(obs),
                    "reward": np.array(rewards),
                    "action": np.array(acs)}
            paths.append(path)

            time_steps_this_batch += len(obs)
            if time_steps_this_batch > args.min_batch_size:
                break

        path_reward = [sum(path['reward']) for path in paths]
        mean_reward = sum(path_reward) / len(path_reward)
        print("Average Reward:", mean_reward)

        ob_batch = np.concatenate([path["observation"] for path in paths])
        ac_batch = np.concatenate([path["action"] for path in paths])
        q_batch = discount_reward(paths, args.gamma, args.reward_to_go)
        summary = model.update(sess, observations=ob_batch, actions=ac_batch,
                               q=q_batch, normalize_advance=args.normalize_advantage,
                               mean_reward=mean_reward, max_loss=args.max_loss)

        if args.debug:
            writer.add_summary(summary, itr)
            logz.pickle_tf_vars()

def main():
    for baseline in [True]:
        for normalize in [True]:
            for reward_to_go in [True]:
                for min_batch in [10000]:
                    for lr in [1e-2]:
                        for max_loss in [1e6]:
                            for seed in [23]:
                                # env_name = 'CartPole-v0'
                                # env_name = 'MountainCar-v0'
                                # env_name = 'MountainCarContinuous-v0'
                                env_name = 'InvertedPendulum-v1'
                                # env_name = "Ant-v1"
                                env_name = 'HalfCheetah-v1'
                                args = MyArgument(env_name=env_name,
                                                  seed=seed,
                                                  debug=True,
                                                  n_layers=3,
                                                  size=64,
                                                  reward_to_go=reward_to_go,
                                                  normalize_advantage=normalize,
                                                  nn_baseline=baseline,
                                                  n_iter=1200,
                                                  learning_rate=lr,
                                                  render=False,
                                                  max_path_length=500,
                                                  min_batch_size=min_batch,
                                                  max_loss=max_loss)
                                p = Process(target=train_pg, args=(args,))
                                p.start()
                                p.join()

if __name__ == "__main__":
    main()

reason: import numpy,import scipy


text:
r"""
===============================================================
Apply Imperfections (:mod:`desicos.abaqus.apply_imperfections`)
===============================================================

.. currentmodule:: desicos.abaqus.apply_imperfections

Routines to apply geometric and thickness imperfections into the finite
element model.

"""
from random import sample

import numpy as np
from numpy import cos, sin, tan, arctan2, deg2rad

from desicos.logger import log
from desicos.constants import FLOAT
from desicos.conecylDB.measured_imp_ms import calc_nodal_translations
from desicos.conecylDB.measured_imp_t import calc_elems_t
from desicos.conecylDB.read_write import read_theta_z_imp
from desicos.conecylDB.interpolate import inv_weighted
from desicos.abaqus.utils import vec_calc_elem_cg, index_within_linspace

def calc_translations_ABAQUS(imperfection_file_name,
                             model_name,
                             part_name,
                             H_model,
                             H_measured,
                             R_model,
                             R_best_fit=None,
                             semi_angle=0.,
                             stretch_H=False,
                             z_offset_bot=None,
                             rotatedeg=0.,
                             scaling_factor=1.,
                             r_TOL=1.,
                             num_closest_points=5,
                             power_parameter=2,
                             use_theta_z_format=True,
                             ignore_bot_h=None,
                             ignore_top_h=None,
                             sample_size=None,
                             T=None):
    r"""Reads an imperfection file and calculates the nodal translations

    Parameters
    ----------
    imperfection_file_name : str
        Full path to the imperfection file.
    model_name : str
        Model name.
    part_name : str
        Part name.
    H_model : float
        Total height of the model where the imperfections will be applied to,
        considering also eventual resin rings.
    H_measured : float
        The total height of the measured test specimen, including eventual
        resin rings at the edges.
    R_model : float
        Radius **at the bottom edge** of the model where the imperfections
        will be applied to.
    R_best_fit : float, optional
        Best fit radius obtained with functions :func:`.best_fit_cylinder`
        or :func:`.best_fit_cone`.
    semi_angle : float, optional
        Cone semi-vertex angle in degrees, when applicable.
    stretch_H : bool, optional
        If the measured imperfection data should be stretched to the current
        model (which may happen when ``H_model!=H_measured``.
    z_offset_bot : float, optional
        It is common to have the measured data not covering the whole test
        specimen, and therefore it will be centralized, if a non-centralized
        position is desired this parameter can be used for the adjustment.
    rotatedeg : float, optional
        Rotation angle in degrees telling how much the imperfection pattern
        should be rotated about the `X_3` (or `Z`) axis.
    scaling_factor : float, optional
        A scaling factor that can be used to study the imperfection
        sensitivity.
    r_TOL : float, optional
        Percentage tolerance to ignore noisy data from the measurements.
    num_closest_points : int, optional
        See :func:`the inverse-weighted interpolation algorithm
        <.inv_weighted>`.
    power_parameter : float, optional
        See :func:`the inverse-weighted interpolation algorithm
        <.inv_weighted>`.
    use_theta_z_format : bool, optional
        If the new format `\theta, Z, imp` should be used instead of the old
        `X, Y, Z`.
    ignore_bot_h : None or float, optional
        Used to ignore nodes from the bottom resin ring.
    ignore_top_h : None or float, optional
        Used to ignore nodes from the top resin ring.
    sample_size : int, optional
        If the input file containing the measured data is too large it may be
        required to limit the sample size.
    T : None or np.ndarray, optional
        A transformation matrix (cf. :func:`.transf_matrix`) required when the
        mesh is not in the :ref:`default coordinate system <figure_conecyl>`.

    """
    import abaqus

    import desicos.abaqus.abaqus_functions as abaqus_functions

    mod = abaqus.mdb.models[model_name]
    part = mod.parts[part_name]
    part_nodes = np.array(part.nodes)
    coords = np.array([n.coordinates for n in part_nodes], dtype=FLOAT)

    if T is not None:
        tmp = np.vstack((coords.T, np.ones((1, coords.shape[0]))))
        coords = np.dot(T, tmp).T
        del tmp

    if ignore_bot_h is not None:
        if ignore_bot_h <= 0:
            ignore_bot_h = None
        else:
            mask = coords[:, 2] > ignore_bot_h
            coords = coords[mask]
            part_nodes = part_nodes[mask]

    if ignore_top_h is not None:
        if ignore_top_h <= 0:
            ignore_top_h = None
        else:
            mask = coords[:, 2] < (H_model - ignore_top_h)
            coords = coords[mask]
            part_nodes = part_nodes[mask]

    if use_theta_z_format:
        d, d, data = read_theta_z_imp(path=imperfection_file_name,
                                      H_measured=H_measured,
                                      stretch_H=stretch_H,
                                      z_offset_bot=z_offset_bot)
        if sample_size:
            num = data.shape[0]
            if sample_size < num:
                log('Using sample_size={0}'.format(sample_size), level=1)
                data = data[sample(range(num), int(sample_size)), :]

        if r_TOL:
            max_imp = R_model * r_TOL / 100.
            imp = data[:, 2]
            cond = np.any(np.array((imp > max_imp, imp < (-max_imp))), axis=0)
            log('Skipping {0} points'.format(len(imp[cond])))
            data = data[np.logical_not(cond), :]

        data3D = np.zeros((data.shape[0], 4), dtype=FLOAT)
        if rotatedeg:
            data[:, 0] += deg2rad(rotatedeg)
        z = data[:, 1]
        z *= H_model

        alpharad = deg2rad(semi_angle)
        tana = tan(alpharad)
        def r_local(z):
            return R_model - z*tana
        data3D[:, 0] = r_local(z)*cos(data[:, 0])
        data3D[:, 1] = r_local(z)*sin(data[:, 0])
        data3D[:, 2] = z
        data3D[:, 3] = data[:, 2]

        dist, w0 = inv_weighted(data3D, coords,
                          ncp = num_closest_points,
                          power_parameter = power_parameter)

        thetas = arctan2(coords[:, 1], coords[:, 0])

        trans = np.zeros_like(coords)
        trans[:, 0] = w0*cos(alpharad)*cos(thetas)
        trans[:, 1] = w0*cos(alpharad)*sin(thetas)
        trans[:, 2] = w0*sin(alpharad)

    else:
        #NOTE perhaps remove this in the future, when the imperfection files
        #     are stored as theta, z, amplitude only
        nodes = np.array([[n.coordinates[0],
                           n.coordinates[1],
                           n.coordinates[2],
                           n.label] for n in part_nodes], dtype=FLOAT)

        # calling translate_nodes function
        trans = calc_nodal_translations(
                                imperfection_file_name,
                                nodes = nodes,
                                H_model = H_model,
                                H_measured = H_measured,
                                R_model = R_model,
                                R_best_fit = R_best_fit,
                                semi_angle = semi_angle,
                                stretch_H = stretch_H,
                                z_offset_bot = z_offset_bot,
                                rotatedeg = rotatedeg,
                                r_TOL = r_TOL,
                                num_closest_points = num_closest_points,
                                power_parameter = power_parameter,
                                sample_size = sample_size)
        trans = trans[:, :3]

    return trans

def translate_nodes_ABAQUS(imperfection_file_name,
                           model_name,
                           part_name,
                           H_model,
                           H_measured,
                           R_model,
                           R_best_fit=None,
                           semi_angle=0.,
                           stretch_H=False,
                           z_offset_bot=None,
                           rotatedeg=0.,
                           scaling_factor=1.,
                           r_TOL=1.,
                           num_closest_points=5,
                           power_parameter=2,
                           nodal_translations=None,
                           use_theta_z_format=False,
                           ignore_bot_h=None,
                           ignore_top_h=None,
                           sample_size=None,
                           T=None):
    r"""Translates the nodes in Abaqus based on imperfection data

    The imperfection amplitude for each node is calculated using an inversed
    weight function (see :func:`desicos.conecylDB.interpolate.inv_weighted`).

    Parameters
    ----------
    imperfection_file_name : str
        The full path to the imperfection file, which must be a file with
        three columns containing the ``x, y, z`` coordinates when
        ``use_theta_z_format=False`` or containing ``x, theta, amplitude``
        when ``use_theta_z_format=True``.
    model_name : str
        Must be a valid key in the dictionary ``mdb.models``, in the
        interactive Python inside Abaqus.
    part_name : str
        Must be a valid key in the dictionary
        ``mdb.models[model_name].parts``, in the interactive Python inside
        Abaqus.
    H_model : float
        Total height of the model where the imperfections will be applied to,
        considering also eventual resin rings.
    H_measured : float
        The total height of the measured test specimen, including eventual
        resin rings at the edges.
    R_model : float
        The radius of the current model. In case of cones this should be the
        bottom radius.
    R_best_fit : float, optional
        Best fit radius obtained with functions :func:`.best_fit_cylinder`
        or :func:`.best_fit_cone`.
    semi_angle : float, optional
        The cone semi-vertex angle (a null value indicates that a cylinder is
        beeing analyzed).
    stretch_H : float, optional
        A boolean indicating if the imperfection pattern should be stretched
        when applied to the model. The measurement systems usually cannot
        obtain data for the whole surface, making it an option to stretch the
        data to fit the whole surface. In case ``stretch_H=False`` the
        measured data of the extremities will be extruded up to the end of the
        domain.
    z_offset_bot : float, optional
        This parameter allows the analyst to adjust the height of the measured
        data about the model, when the measured data is not available for the
        whole domain.
    rotatedeg : float, optional
        Rotation angle in degrees telling how much the imperfection pattern
        should be rotated about the `X_3` (or `Z`) axis.
    scaling_factor : float, optional
        The scaling factor that will multiply the calculated imperfection
        amplitude.
    r_TOL : float, optional
        Parameter to ignore noisy data in the imperfection file, the points
        with a radius higher than `r=r_{model} \times (1 + r_{TOL})` will not
        be considered in the interpolation.
    num_closest_points : int, optional
        See :func:`the inverse-weighted interpolation algorithm
        <.inv_weighted>`.
    power_parameter : int, optional
        See :func:`the inverse-weighted interpolation algorithm
        <.inv_weighted>`.
    nodal_translations : None or numpy.ndarray, optional
        An array containing the interpolated traslations, which is passed to
        avoid repeated calls to the interpolation functions.
    use_theta_z_format : bool, optional
        A boolean to indicate whether the imperfection file contains ``x, y,
        z`` positions or ``theta, z, amplitude``.
    ignore_bot_h : None or float, optional
        Used to ignore nodes from the bottom resin ring.
    ignore_top_h : None or float, optional
        Used to ignore nodes from the top resin ring.
    sample_size : int, optional
        If the input file containing the measured data is too large it may be
        required to limit the sample size.
    T : None or np.ndarray, optional
        A transformation matrix (cf. :func:`.transf_matrix`) required when the
        mesh is not in the :ref:`default coordinate system <figure_conecyl>`.

    Returns
    -------
    nodal_translations : numpy.ndarray
        A 2-D array containing the translations ``x, y, z`` for each column.

    Notes
    -----
    Despite the nodal traslations are returned all the nodes belonging to this
    model will already be translated.

    """
    from abaqus import mdb, session

    mod = mdb.models[model_name]
    part = mod.parts[part_name]

    part_nodes = np.array(part.nodes)
    coords = np.array([n.coordinates for n in part_nodes])

    if T is not None:
        tmp = np.vstack((coords.T, np.ones((1, coords.shape[0]))))
        coords = np.dot(T, tmp).T
        del tmp

    if ignore_bot_h is not None:
        if ignore_bot_h <= 0:
            ignore_bot_h = None
        else:
            log('Applying ignore_bot_h: ignoring nodes with z <= {0}'.format(
                ignore_bot_h))
            mask = coords[:, 2] > ignore_bot_h
            coords = coords[mask]
            part_nodes = part_nodes[mask]

    if ignore_top_h is not None:
        if ignore_top_h <= 0:
            ignore_top_h = None
        else:
            log('Applying ignore_top_h: ignoring nodes with z >= {0}'.format(
                H_model - ignore_top_h))
            mask = coords[:, 2] < (H_model - ignore_top_h)
            coords = coords[mask]
            part_nodes = part_nodes[mask]

    if use_theta_z_format:
        if nodal_translations is None:
            trans = calc_translations_ABAQUS(
                        imperfection_file_name = imperfection_file_name,
                        model_name = model_name,
                        part_name = part_name,
                        H_model = H_model,
                        H_measured = H_measured,
                        R_model = R_model,
                        R_best_fit = R_best_fit,
                        semi_angle = semi_angle,
                        stretch_H = stretch_H,
                        z_offset_bot = z_offset_bot,
                        rotatedeg = rotatedeg,
                        scaling_factor = scaling_factor,
                        r_TOL = r_TOL,
                        num_closest_points = num_closest_points,
                        power_parameter = power_parameter,
                        use_theta_z_format = use_theta_z_format,
                        ignore_bot_h = ignore_bot_h,
                        ignore_top_h = ignore_top_h,
                        sample_size = sample_size,
                        T = T)

        else:
            trans = nodal_translations

        # applying translations
        viewport = session.viewports[session.currentViewportName]
        log('Applying new nodal positions in ABAQUS CAE to model {0} ...'.
            format(model_name))
        log('    (using a scaling factor of {0})'.format(scaling_factor))

        new_coords = coords + trans*scaling_factor

        if T is not None:
            Tinv = np.zeros_like(T)
            Tinv[:3, :3] = T[:3, :3].T
            Tinv[:, 3] = -T[:, 3]
            tmp = np.vstack((new_coords.T, np.ones((1, new_coords.shape[0]))))
            new_coords = np.dot(Tinv, tmp).T
            del tmp

        meshNodeArray = part.nodes.sequenceFromLabels(
                            [n.label for n in part_nodes])
        new_coords = np.ascontiguousarray(new_coords)
        part.editNode(nodes=part_nodes.tolist(), coordinates=new_coords)

        log('Application of new nodal positions finished!')

        ra = mod.rootAssembly
        viewport.setValues(displayedObject=ra)
        ra.regenerate()

        return trans

    else:
        nodes = np.array([[n.coordinates[0],
                           n.coordinates[1],
                           n.coordinates[2],
                           n.label] for n in part_nodes], dtype=FLOAT)

        # calling translate_nodes function
        if nodal_translations is None:
            nodal_translations = calc_translations_ABAQUS(
                         imperfection_file_name = imperfection_file_name,
                         model_name = model_name,
                         part_name = part_name,
                         H_model = H_model,
                         H_measured = H_measured,
                         R_model = R_model,
                         R_best_fit = R_best_fit,
                         semi_angle = semi_angle,
                         stretch_H = stretch_H,
                         z_offset_bot = z_offset_bot,
                         rotatedeg = rotatedeg,
                         scaling_factor = scaling_factor,
                         r_TOL = r_TOL,
                         num_closest_points = num_closest_points,
                         power_parameter = power_parameter,
                         use_theta_z_format = use_theta_z_format,
                         ignore_bot_h = ignore_bot_h,
                         ignore_top_h = ignore_top_h,
                         sample_size = sample_size,
                         T = T)

        # applying translations
        viewport = session.viewports[session.currentViewportName]
        log('Applying new nodal positions in ABAQUS CAE for model {0} ...'.
            format(model_name))
        log('    (using a scaling factor of {0})'.format(scaling_factor))

        trans = nodal_translations

        new_coords = coords + trans*scaling_factor

        if T is not None:
            Tinv = np.zeros_like(T)
            Tinv[:3, :3] = T[:3, :3].T
            Tinv[:, 3] = -T[:, 3]
            tmp = np.vstack((new_coords.T, np.ones((1, new_coords.shape[0]))))
            new_coords = np.dot(Tinv, tmp).T
            del tmp

        meshNodeArray = part.nodes.sequenceFromLabels(
                            [n.label for n in part_nodes])
        new_coords = np.ascontiguousarray(new_coords)
        part.editNode(nodes=meshNodeArray, coordinates=new_coords)

        log('Application of new nodal positions finished!')
        # regenerating ra
        ra = mod.rootAssembly
        viewport.setValues(displayedObject=ra)
        ra.regenerate()

        return nodal_translations

def translate_nodes_ABAQUS_c0(m0, n0, c0, funcnum,
                              model_name,
                              part_name,
                              H_model,
                              semi_angle=0.,
                              scaling_factor=1.,
                              fem_meridian_bot2top=True,
                              ignore_bot_h=None,
                              ignore_top_h=None,
                              T=None):
    r"""Translates the nodes in Abaqus based on a Fourier series

    The Fourier Series can be a half-sine, half-cosine or a complete Fourier
    Series as detailed in :func:`desicos.conecylDB.fit_data.calc_c0`.

    Parameters
    ----------
    m0 : int
        Number of terms along the `x` coordinate.
    n0 : int
        Number of terms along the `\theta` coordinate.
    c0 : numpy.ndarray
        The coefficients that will give the imperfection pattern.
    funcnum : int
        The function type, as detailed in
        :func:`desicos.conecylDB.fit_data.calc_c0`.
    model_name : str
        Must be a valid key in the dictionary ``mdb.models``, in the
        interactive Python inside Abaqus.
    part_name : str
        Must be a valid key in the dictionary
        ``mdb.models[model_name].parts``, in the interactive Python inside
        Abaqus.
    H_model : float
        Total height of the model where the imperfections will be applied to,
        considering also eventual resin rings.
    semi_angle : float, optional
        The cone semi-vertex angle (a null value indicates that a cylinder is
        beeing analyzed).
    scaling_factor : float, optional
        The scaling factor that will multiply ``c0`` when applying the
        imperfections.
    fem_meridian_bot2top : bool, optional
        A boolean indicating if the finite element has the `x` axis starting
        at the bottom or at the top.
    ignore_bot_h : None or float, optional
        Used to ignore nodes from the bottom resin ring.
    ignore_top_h : None or float, optional
        Used to ignore nodes from the top resin ring.
    T : None or np.ndarray, optional
        A transformation matrix (cf. :func:`.transf_matrix`) required when the
        mesh is not in the :ref:`default coordinate system <figure_conecyl>`.

    Returns
    -------
    nodal_translations : numpy.ndarray
        A 2-D array containing the translations ``x, y, z`` for each column.

    Notes
    -----
    Despite the nodal traslations are returned all the nodes belonging to this
    model will be already translated.

    """
    from abaqus import mdb, session

    from desicos.conecylDB import fit_data

    log('Calculating imperfection amplitudes for model {0} ...\n\t'.
        format(model_name) +
        '(using a scaling factor of {0})'.format(scaling_factor))
    c0 = c0*scaling_factor
    mod = mdb.models[model_name]
    part = mod.parts[part_name]

    part_nodes = np.array(part.nodes)
    coords = np.array([n.coordinates for n in part_nodes])

    if T is not None:
        tmp = np.vstack((coords.T, np.ones((1, coords.shape[0]))))
        coords = np.dot(T, tmp).T
        del tmp

    if ignore_bot_h is not None:
        if ignore_bot_h <= 0:
            ignore_bot_h = None
        else:
            log('Applying ignore_bot_h: ignoring nodes with z <= {0}'.format(
                ignore_bot_h))
            mask = coords[:, 2] > ignore_bot_h
            coords = coords[mask]
            part_nodes = part_nodes[mask]

    if ignore_top_h is not None:
        if ignore_top_h <= 0:
            ignore_top_h = None
        else:
            log('Applying ignore_top_h: ignoring nodes with z >= {0}'.format(
                H_model - ignore_top_h))
            mask = coords[:, 2] < (H_model - ignore_top_h)
            coords = coords[mask]
            part_nodes = part_nodes[mask]

    if ignore_bot_h is None:
        ignore_bot_h = 0
    if ignore_top_h is None:
        ignore_top_h = 0
    H_eff = H_model - (ignore_bot_h + ignore_top_h)
    if fem_meridian_bot2top:
        xs_norm = (coords[:, 2]-ignore_bot_h)/H_eff
    else:
        xs_norm = (H_eff-(coords[:, 2]-ignore_bot_h))/H_eff

    thetas = arctan2(coords[:, 1], coords[:, 0])

    alpharad = deg2rad(semi_angle)
    w0 = fit_data.fw0(m0, n0, c0, xs_norm, thetas, funcnum)
    nodal_translations = np.zeros_like(coords)
    nodal_translations[:, 0] = w0*cos(alpharad)*cos(thetas)
    nodal_translations[:, 1] = w0*cos(alpharad)*sin(thetas)
    nodal_translations[:, 2] = w0*sin(alpharad)

    log('Calculation of imperfection amplitudes finished!')

    # applying translations
    viewport = session.viewports[session.currentViewportName]
    log('Applying new nodal positions in ABAQUS CAE to model {0} ...'.
        format(model_name))
    log('    (using a scaling factor of {0})'.format(scaling_factor))

    new_coords = coords + nodal_translations

    if T is not None:
        Tinv = np.zeros_like(T)
        Tinv[:3, :3] = T[:3, :3].T
        Tinv[:, 3] = -T[:, 3]
        tmp = np.vstack((new_coords.T, np.ones((1, new_coords.shape[0]))))
        new_coords = np.dot(Tinv, tmp).T
        del tmp

    meshNodeArray = part.nodes.sequenceFromLabels(
                        [n.label for n in part_nodes])
    new_coords = np.ascontiguousarray(new_coords)
    part.editNode(nodes=meshNodeArray, coordinates=new_coords)

    log('Application of new nodal positions finished!')

    ra = mod.rootAssembly
    viewport.setValues(displayedObject=ra)
    ra.regenerate()

    return nodal_translations

def change_thickness_ABAQUS(imperfection_file_name,
                            model_name,
                            part_name,
                            t_model,
                            t_measured,
                            H_model,
                            H_measured,
                            R_model,
                            R_best_fit = None,
                            number_of_sets = None,
                            semi_angle = 0.,
                            stretch_H = False,
                            z_offset_bot = None,
                            scaling_factor = 1.,
                            num_closest_points = 5,
                            power_parameter = 2,
                            elems_t = None,
                            t_set = None,
                            use_theta_z_format = False):
    r"""Applies a given thickness imperfection to the finite element model

    Assumes that a percentage variation of the laminate thickness can be
    represented by the same percentage veriation of each ply, i.e., each
    ply thickness is varied in order to reflect a given measured thickness
    imperfection field.

    Parameters
    ----------

    imperfection_file_name : str
        Full path to the imperfection file.
    model_name : str
        Model name.
    part_name : str
        Part name.
    t_model : float
        The nominal shell thickness of the current model.
    t_measured : float
        The nominal thickness of the measured specimen.
    H_model : float
        Total height of the model where the imperfections will be applied to,
        considering also eventual resin rings.
    H_measured : float
        The total height of the measured test specimen, including eventual
        resin rings at the edges.
    R_model : float
        Radius **at the bottom edge** of the model where the imperfections
        will be applied to.
    R_best_fit : float, optional
        Best fit radius obtained with functions :func:`.best_fit_cylinder`
        or :func:`.best_fit_cone`.
    number_of_sets : int, optional
        Defines in how many levels the thicknesses should be divided. If
        ``None`` it will be based on the input file, and if the threshold
        of ``100`` is exceeded, ``10`` sections are used.
    semi_angle : float, optional
        Cone semi-vertex angle in degrees, when applicable.
    stretch_H : bool, optional
        If the measured imperfection data should be stretched to the current
        model (which may happen when ``H_model!=H_measured``.
    z_offset_bot : float, optional
        It is common to have the measured data not covering the whole test
        specimen, and therefore it will be centralized, if a non-centralized
        position is desired this parameter can be used for the adjustment.
    scaling_factor : float, optional
        A scaling factor that can be used to study the imperfection
        sensitivity.
    num_closest_points : int, optional
        See :func:`the inverse-weighted interpolation algorithm
        <.inv_weighted>`.
    power_parameter : float, optional
        See :func:`the inverse-weighted interpolation algorithm
        <.inv_weighted>`.
    elems_t : np.ndarray, optional
        Interpolated thickness for each element. Can be used to avoid the same
        interpolation to be performed twice.
    t_set : set, optional
        A ``set`` object containing the unique thicknesses that will be used
        to create the new properties.
    use_theta_z_format : bool, optional
        If the new format `\theta, Z, imp` should be used instead of the old
        `X, Y, Z`.

    """
    from abaqus import mdb

    import desicos.abaqus.abaqus_functions as abaqus_functions

    mod = mdb.models[model_name]
    part = mod.parts[part_name]
    part_cyl_csys = part.features['part_cyl_csys']
    part_cyl_csys = part.datums[part_cyl_csys.id]

    if use_theta_z_format:
        if elems_t is None or t_set is None:
            log('Reading coordinates for elements...')
            elements = vec_calc_elem_cg(part.elements)

            log('Coordinates for elements read!')
            d, d, data = read_theta_z_imp(path = imperfection_file_name,
                                          H_measured = H_measured,
                                          stretch_H = stretch_H,
                                          z_offset_bot = z_offset_bot)

            data3D = np.zeros((data.shape[0], 4), dtype=FLOAT)
            z = data[:, 1]
            z *= H_model

            alpharad = deg2rad(semi_angle)
            tana = tan(alpharad)
            def r_local(z):
                return R_model - z*tana
            data3D[:, 0] = r_local(z)*cos(data[:, 0])
            data3D[:, 1] = r_local(z)*sin(data[:, 0])
            data3D[:, 2] = z
            data3D[:, 3] = data[:, 2]

            dist, ans = inv_weighted(data3D, elements[:, :3],
                               ncp = num_closest_points,
                               power_parameter = power_parameter)

            t_set = set(ans)
            t_set.discard(0.) #TODO why does inv_weighted return an array with 0.?
            elems_t = np.zeros((elements.shape[0], 2), dtype=FLOAT)
            elems_t[:, 0] = elements[:, 3]
            elems_t[:, 1] = ans

        else:
            log('Thickness differences already calculated!')

    else:
        if elems_t is None or t_set is None:
            # reading elements data
            log('Reading coordinates for elements...')
            elements = vec_calc_elem_cg(part.elements)
            log('Coordinates for elements read!')
            # calling translate_nodes function
            elems_t, t_set = calc_elems_t(
                                imperfection_file_name,
                                nodes = elements,
                                t_model = t_model,
                                t_measured = t_measured,
                                H_model = H_model,
                                H_measured = H_measured,
                                R_model = R_model,
                                R_best_fit = R_best_fit,
                                semi_angle = semi_angle,
                                stretch_H = stretch_H,
                                z_offset_bot = z_offset_bot,
                                num_closest_points = num_closest_points,
                                power_parameter = power_parameter
                                )
        else:
            log('Thickness differences already calculated!')
    # creating sets
    t_list = []
    max_len_t_set = 100
    if len(t_set) >= max_len_t_set and number_of_sets in (None, 0):
        number_of_sets = 10
        log('More than {0:d} different thicknesses measured!'.format(
            max_len_t_set))
        log('Forcing a number_of_sets = {0:d}'.format(number_of_sets))
    if number_of_sets is None or number_of_sets == 0:
        number_of_sets = len(t_set)
        t_list = list(t_set)
        t_list.sort()
    else:
        t_min = min(t_set)
        t_max = max(t_set)
        t_list = list(np.linspace(t_min, t_max, number_of_sets+1))

    # grouping elements
    sets_ids = [[] for i in range(len(t_list))]
    for entry in elems_t:
        elem_id, t = entry
        index = index_within_linspace(t_list, t)
        sets_ids[index].append(int(elem_id))
    # putting elements in sets
    original_layup = part.compositeLayups['CompositePlate']
    for i, set_ids in enumerate(sets_ids):
        if len(set_ids) == 0:
            # since t_set_norm * t_model <> t_set originally measured
            # there may be empty set_ids at the end
            continue
        elements = part.elements.sequenceFromLabels(labels=set_ids)
        suffix = 'measured_imp_t_{0:03d}'.format(i)
        set_name = 'Set_' + suffix
        log('Creating set ({0: 7d} elements): {1}'.format(
            len(set_ids), set_name))
        part.Set(name = set_name, elements = elements)
        region = part.sets[set_name]
        layup_name = 'CLayup_' + suffix
        t_diff = (float(t_list[i]) - t_model) * scaling_factor
        t_scaling_factor = (t_model + t_diff)/t_model

        def modify_ply(index, kwargs):
            kwargs['thickness'] *= t_scaling_factor
            kwargs['region'] = region
            return kwargs

        layup = part.CompositeLayup(name=layup_name,
                                    objectToCopy=original_layup)
        layup.resume()
        abaqus_functions.modify_composite_layup(part=part,
            layup_name=layup_name, modify_func=modify_ply)
    # suppress needed to put the new properties to the input file
    original_layup.suppress()

    return elems_t, t_set

def change_thickness_per_element_ABAQUS(
                            elems_t,
                            model_name,
                            part_name,
                            scaling_factor=1.,
                            ):
    r"""Applies a given thickness imperfection to the finite element model

    Assumes that a percentage variation of the laminate thickness can be
    represented by the same percentage veriation of each ply, i.e., each
    ply thickness is varied in order to reflect a given measured thickness
    imperfection field.

    Parameters
    ----------

    elems_t : np.ndarray, optional
        Interpolated thickness for each element. Can be used to avoid the same
        interpolation to be performed twice.
    model_name : str
        Model name.
    part_name : str
        Part name.
    scaling_factor : float, optional
        A scaling factor that can be used to study the imperfection
        sensitivity.

    """
    from abaqus import mdb

    import desicos.abaqus.abaqus_functions as abaqus_functions

    mod = mdb.models[model_name]
    part = mod.parts[part_name]
    part_cyl_csys = part.features['part_cyl_csys']
    part_cyl_csys = part.datums[part_cyl_csys.id]

    # creating sets
    t_list = elems_t[:, 1]
    # grouping elements
    sets_ids = [[] for i in range(len(t_list))]
    for i, entry in enumerate(elems_t):
        elem_id, t = entry
        sets_ids[i].append(int(elem_id))
    # finding original layup for each element
    comp_layup_elems = {}
    for name, composite_layup in part.compositeLayups.items():
        comp_layup_elems[name] = set()
        for ply in composite_layup.plies.values():
            for e in part.sets[ply.region[0]].elements:
                comp_layup_elems[name].add(e.label)
    elem_layup = {}
    for eid in elems_t[:, 0].astype(int):
        found = False
        for name, eids in comp_layup_elems.items():
            if eid in eids:
                elem_layup[elem_id] = name
                found = True
                break
        if not found:
            raise RuntimeError('Element ID %d not assigned to any composite layup' % eid)

    # putting elements in sets
    for i, set_ids in enumerate(sets_ids):
        if len(set_ids) == 0:
            # since t_set_norm * t_model <> t_set originally measured
            # there may be empty set_ids at the end
            continue
        elements = part.elements.sequenceFromLabels(labels=set_ids)
        suffix = 'measured_imp_t_{0:07d}'.format(i)
        set_name = 'Set_' + suffix
        log('Creating set ({0: 7d} elements): {1}'.format(
            len(set_ids), set_name))
        part.Set(name = set_name, elements = elements)
        region = part.sets[set_name]
        layup_name = 'CLayup_' + suffix
        compositeLayup = part.compositeLayups[name]
        t_model = sum([ply.thickness for ply in compositeLayup.plies.values()])
        t_diff = (float(t_list[i]) - t_model) * scaling_factor
        t_scaling_factor = (t_model + t_diff)/t_model

        def modify_ply(index, kwargs):
            kwargs['thickness'] *= t_scaling_factor
            kwargs['region'] = region
            return kwargs

        layup = part.CompositeLayup(name=layup_name,
                                    objectToCopy=compositeLayup)
        layup.resume()
        abaqus_functions.modify_composite_layup(part=part,
            layup_name=layup_name, modify_func=modify_ply)
    # suppress needed to put the new properties to the input file
    for name, composite_layup in part.compositeLayups.items():
        if name in comp_layup_elems:
            if len(comp_layup_elems[name]) > 0:
                composite_layup.suppress()

if __name__ == '__main__':
    cc = stds['desicos_study'].ccs[1]
    cc.created_model = False
    cc.create_model()
    part_nodes = mdb.models[cc.model_name].parts[cc.part_name_shell].nodes
    coords = np.array([n.coordinates for n in part_nodes])

    sf = 100

    path = r'C:\clones\desicos\desicos\conecylDB\files\dlr\degenhardt_2010_z20\degenhardt_2010_z20_msi_theta_z_imp.txt'
    H_measured = 510.
    H_model = 510.
    d, d, data = read_theta_z_imp(path=path,
                                  H_measured=H_measured,
                                  stretch_H=False,
                                  z_offset_bot=None)

    log('init sample')
    data = np.array(sample(data, 10000))
    log('end sample')
    data = data[np.argsort(data[:, 0])]
    data = np.vstack((data[-3000:], data[:], data[:3000]))
    data[:, 0] /= data[:, 0].max()

    mesh = np.zeros((coords.shape[0], 2), dtype=FLOAT)
    mesh[:, 0] = arctan2(coords[:, 1], coords[:, 0])
    mesh[:, 1] = coords[:, 2]
    mesh_norm = mesh.copy()
    mesh_norm[:, 0] /= mesh_norm[:, 0].max()
    mesh_norm[:, 1] /= mesh_norm[:, 1].max()

    import desicos.conecylDB.interpolate
    reload(desicos.conecylDB.interpolate)
    dist, ans = desicos.conecylDB.interpolate.inv_weighted(
            data, mesh_norm, ncp=10, power_parameter=1.5)

    alpharad = deg2rad(0.)
    trans = np.zeros_like(coords)
    trans[:, 0] = ans*cos(mesh[:, 0])*cos(alpharad)
    trans[:, 1] = ans*sin(mesh[:, 0])*cos(alpharad)
    trans[:, 2] = ans*sin(alpharad)

    meshNodeArray = part.nodes.sequenceFromLabels(
                        [n.label for n in part_nodes])
    cc.part.editNode(nodes=meshNodeArray,
                  coordinates=(coords + trans*100))

max_stars_repo_name: saullocastro/desicos

max_issues_repo_name: saullocastro/desicos

max_forks_repo_name: saullocastro/desicos

reason: import numpy,from numpy


text:
"""

Generates cost Hamiltonian for the maxcut problem

from qiskit.aqua/translators/ising/maxcut.py

# Copyright 2018 IBM.
#
# 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
from qiskit.quantum_info import Pauli
from qiskit.aqua.operators import WeightedPauliOperator

def get_maxcut_qubitops(weight_matrix):
    """Generate Hamiltonian for the max-cut problem of a graph.
    Args:
        weight_matrix (numpy.ndarray) : adjacency matrix.
    Returns:
        WeightedPauliOperator: operator for the Hamiltonian
        float: a constant shift for the obj function.
    """
    num_nodes = weight_matrix.shape[0]
    pauli_list = []
    shift = 0
    for i in range(num_nodes):
        for j in range(i):
            if weight_matrix[i, j] != 0:
                x_p = np.zeros(num_nodes, dtype=np.bool)
                z_p = np.zeros(num_nodes, dtype=np.bool)
                z_p[i] = True
                z_p[j] = True
                pauli_list.append([0.5 * weight_matrix[i, j], Pauli(z_p, x_p)])
                shift -= 0.5 * weight_matrix[i, j]
    return WeightedPauliOperator(paulis=pauli_list), shift

max_stars_repo_name: rsln-s/ibmqxbackend

max_issues_repo_name: rsln-s/ibmqxbackend

max_forks_repo_name: rsln-s/ibmqxbackend


text:
# This code is part of Mapomatic.
#
# (C) Copyright IBM 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""Test best mappings"""
import numpy as np
from qiskit import transpile, QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.test.mock import FakeBelem, FakeQuito, FakeLima

import mapomatic as mm

def test_best_mapping_ghz_state_full_device_multiple_qregs():
    """Test best mappings with multiple registers"""
    qr_a = QuantumRegister(2)
    qr_b = QuantumRegister(3)
    qc = QuantumCircuit(qr_a, qr_b)
    qc.h(qr_a[0])
    qc.cx(qr_a[0], qr_a[1])
    qc.cx(qr_a[0], qr_b[0])
    qc.cx(qr_a[0], qr_b[1])
    qc.cx(qr_a[0], qr_b[2])
    qc.measure_all()
    trans_qc = transpile(qc, FakeLima(), seed_transpiler=102442)
    backends = [FakeBelem(), FakeQuito(), FakeLima()]
    res = mm.best_overall_layout(trans_qc, backends, successors=True)
    expected_res = [
        ([0, 1, 2, 3, 4], 'fake_lima', 0.2948138379010775),
        ([0, 1, 2, 3, 4], 'fake_belem', 0.3099503939385677),
        ([2, 1, 0, 3, 4], 'fake_quito', 0.5360875795095078)
    ]
    for index, expected in enumerate(expected_res):
        assert res[index][0] == expected[0]
        assert res[index][1] == expected[1]
        assert np.allclose(res[index][2], expected[2])

def test_best_mapping_ghz_state_deflate_multiple_registers():
    """Test best mappings with multiple registers with deflate"""
    qr_a = QuantumRegister(2)
    qr_b = QuantumRegister(2)
    cr_a = ClassicalRegister(2)
    cr_b = ClassicalRegister(2)
    qc = QuantumCircuit(qr_a, qr_b, cr_a, cr_b)
    qc.h(qr_a[0])
    qc.cx(qr_a[0], qr_a[1])
    qc.cx(qr_a[0], qr_b[0])
    qc.cx(qr_a[0], qr_b[1])
    qc.measure(qr_a, cr_b)
    qc.measure(qr_b, cr_a)
    trans_qc = transpile(qc, FakeLima(), seed_transpiler=102442)
    small_circ = mm.deflate_circuit(trans_qc)
    backends = [FakeBelem(), FakeQuito(), FakeLima()]
    res = mm.best_overall_layout(small_circ, backends, successors=True)
    expected_res = [
        ([3, 1, 0, 2], 'fake_lima', 0.1466490604029853),
        ([0, 1, 3, 2], 'fake_belem', 0.18757249682201993),
        ([3, 1, 2, 0], 'fake_quito', 0.3202504720264385)
    ]
    for index, expected in enumerate(expected_res):
        assert res[index][0] == expected[0]
        assert res[index][1] == expected[1]
        assert np.allclose(res[index][2], expected[2])

max_stars_repo_name: nonhermitian/mapomatic

max_issues_repo_name: nonhermitian/mapomatic

max_forks_repo_name: nonhermitian/mapomatic


text:
from numpy import *
import matplotlib.pyplot as plt
def loadSimpData():
    datMat = array([[1.,2.1],
                     [2.,1.1],
                     [1.3,1.],
                     [1.,1.],
                     [2.,1.]])
    classLabels = [1.0,1.0,-1.0,-1.0,1.0]
    return datMat, classLabels
dataMat, classLabels = loadSimpData()

# 最佳单层决策树生成函数
def stumpClassify(dataMatrix, dimen, threshVal, threshIneq):
    retArray = ones((shape(dataMatrix)[0], 1))
    if threshIneq == 'lt':
        retArray[dataMatrix[:, dimen] <= threshVal] = -1.0
    else:
        retArray[dataMatrix[:, dimen] > threshVal] = -1.0
    return retArray

def buildStump(dataArr, classLabels, D):
    dataMatrix = mat(dataArr);labelMat = mat(classLabels).T
    m,n = shape(dataMatrix)
    numSteps = 10.0; bestStump = {}; bestClassEst = mat(zeros((m, 1)))
    minError = inf
    for i in range(n):
        rangeMin = dataMatrix[:, i].min();rangeMax = dataMatrix[:, i].max();
        stepSize = (rangeMax - rangeMin) / numSteps
        for j in range(-1, int(numSteps) + 1):
            for inequal in ['lt', 'gt']:
                threshVal = (rangeMin + float(j) * stepSize)
                predictVals = stumpClassify(dataMatrix, i, threshVal, inequal)
                errArr = mat(ones((m, 1)))
                errArr[predictVals == labelMat] = 0
                weightedError = D.T * errArr
                # print("split: dim %d, thresh %.2f, thresh ineqal: %s, the weighted error is %.3f" % (i, threshVal, inequal, weightedError))
                if weightedError < minError:
                    minError = weightedError
                    bestClassEst = predictVals.copy()
                    bestStump['dim'] = i
                    bestStump['thresh'] = threshVal
                    bestStump['ineq'] = inequal
    return bestStump, minError, bestClassEst

# 测试
# # 初始权重 1/N
# D = mat(ones((5,1)) / 5)
# print(buildStump(dataMat, classLabels, D))

# 基于单层决策树的AdaBoost训练过程
# dataArr 数据集, classLabels 类别标签, numlt = 40 迭代次数（默认40）
def adaBoostTrainDS(dataArr, classLabels, numlt = 40):
    weakClassArr = []
    m = shape(dataArr)[0]
    D = mat(ones((m,1)) / m)
    # aggClassEst记录每个数据点的类别估计累计值
    aggClassEst = mat(zeros((m, 1)))
    for i in range(numlt):
        bestStump, error, classEst = buildStump(dataArr, classLabels, D)
        print("D:", D.T)
        # max(error, 1e-16))用于确保在没有错误时不会发生除0溢出
        alpha = float(0.5 * log((1.0 - error) / max(error, 1e-16))) # 科学计数法 即exp(-16)
        bestStump['alpha'] = alpha
        weakClassArr.append(bestStump)
        print("classEst:", classEst.T)
        expon = multiply(-1 * alpha * mat(classLabels).T, classEst)
        D = multiply(D, exp(expon))
        D = D / D.sum()
        aggClassEst += alpha * classEst
        print("aggClassEst:", aggClassEst.T)
        aggErrors = multiply(sign(aggClassEst) != mat(classLabels).T, ones((m, 1)))
        errorRate = aggErrors.sum() / m
        print("total error:", errorRate, "\n")
        # 训练错误为0时，退出循环
        if errorRate == 0.0: break
    return weakClassArr,aggClassEst
# 测试
# classifierArray = adaBoostTrainDS(dataMat, classLabels, 9)
# print(classifierArray)

# AdaBoost分类函数
def adaClassify(dataToClass, classifierArray):
    dataMatrix = mat(dataToClass)
    m = shape(dataMatrix)[0]
    aggClassEst = mat(zeros((m,1)))
    for i in range(len(classifierArray)):
        classEst = stumpClassify(dataMatrix, classifierArray[i]['dim'], classifierArray[i]['thresh'], classifierArray[i]['ineq'])
        aggClassEst += classifierArray[i]['alpha'] * classEst
        print(aggClassEst)
    return sign(aggClassEst)

# 测试
# print(adaClassify([0,0],classifierArray))

# 马疝病的预测

# 自适应数据加载函数
def loadDataSet(fileName):
    numFeat = len(open(fileName).readline().split('\t'))
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = []
        curLine = line.strip().split('\t')
        for i in range(numFeat - 1):
            lineArr.append(float(curLine[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    return dataMat, labelMat

# 测试
# dataArr, labelArr = loadDataSet('horseColicTraining2.txt')
# classifierArray = adaBoostTrainDS(dataArr, labelArr, 10)
# print(classifierArray)
# print("-------------------------------------------------")
# testArr, testLabelArr = loadDataSet('horseColicTest2.txt')
# prediction10 = adaClassify(testArr, classifierArray)
# errArr = mat(ones((67,1)))
# print(errArr[prediction10 != mat(testLabelArr).T].sum())

# ROC曲线的绘制及AUC计算函数
def plotROC(predStrengths, classLabels):
    cur = (1.0, 1.0)  # cursor
    ySum = 0.0  # variable to calculate AUC
    numPosClas = sum(array(classLabels) == 1.0)
    yStep = 1 / float(numPosClas);
    xStep = 1 / float(len(classLabels) - numPosClas)
    sortedIndicies = predStrengths.argsort()  # get sorted index, it's reverse
    fig = plt.figure()
    fig.clf()
    ax = plt.subplot(111)
    # loop through all the values, drawing a line segment at each point
    for index in sortedIndicies.tolist()[0]:
        if classLabels[index] == 1.0:
            delX = 0;
            delY = yStep;
        else:
            delX = xStep;
            delY = 0;
            ySum += cur[1]
        # draw line from cur to (cur[0]-delX,cur[1]-delY)
        ax.plot([cur[0], cur[0] - delX], [cur[1], cur[1] - delY], c='b')
        cur = (cur[0] - delX, cur[1] - delY)
    ax.plot([0, 1], [0, 1], 'b--')
    plt.xlabel('False positive rate'); plt.ylabel('True positive rate')
    plt.title('ROC curve for AdaBoost horse colic detection system')
    ax.axis([0, 1, 0, 1])
    plt.show()
    print("the Area Under the Curve is: ", ySum * xStep)

# 测试
dataArr, labelArr = loadDataSet('horseColicTraining2.txt')
classifierArray, aggClassEst = adaBoostTrainDS(dataArr, labelArr, 10)
plotROC(aggClassEst.T, labelArr)

reason: from numpy