text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import base64
import time
from line_profiler import LineProfiler
from turtle import width
from difflib import SequenceMatcher
import re
from tqdm.std import tqdm
from urllib.parse import quote, unquote
from datetime import datetime
import pickle
import shutil
import glob
import jsonlines
import json
import csv
import os
import pandas as pd
import traceback
from numpy import random

def readlines(path, portion=1, end: int = None, show_bar=True, encoding="utf-8"):
    """将txt的每一行读出来

    Args:
        txt_path (string): 文件路径
        portion (int, optional): 读出文件比例，比如只需要对一部分看看. Defaults to 1.
        end (int, optional): 指定读到哪一行为止. Defaults to None.
        showBar (bool, optional): 是否显示进度条. Defaults to True.
        encoding (str,optional): 'UTF-8-sig'  "utf-8" 'urt8'
    Returns:
        list: each elements is a row
    """
    data_file = []
    with open(path, "r+", encoding=encoding) as f:
        if end is None:
            num_lines = len([1 for line in open(path, "r", encoding=encoding)])
            num_lines *= portion
        else:
            num_lines = end
        if show_bar:
            for idx, item in enumerate(
                    tqdm(f, total=num_lines, desc=path.split(os.sep)[-1] + " is loading...")):
                if idx >= num_lines:
                    break
                data_file.append(item.replace("\n", ""))
        else:
            for idx, item in enumerate(f):
                if idx >= num_lines:
                    break
                data_file.append(item.replace("\n", ""))
    return data_file

def glob_read(path, read_fun, stop_i=None,  show_bar=False):
    """read all files in path by glob

    Args:
        path (str): absolute path
        read_fun ([type]): different read function based on file type
        stop_i (int, optional): stop read at file i. Defaults to None.
        show_bar (bool, optional): display read process or not . Defaults to False.

    Returns:
        text_list (list): a list of lines in all files
        pathlist (list): a list of absolute file path

    Examples:
        >>> text_list, pathlist = glob_read("/home/directory/test", read_fun=read_txt)
        >>> print(text_list)
        ['line1', 'line2', 'line3']
    """
    text_list = []
    pathlist = glob.glob(path)
    if show_bar:
        w = tqdm(pathlist, desc=u'已加载0个text')
    else:
        w = pathlist
    for i, txt in enumerate(w):
        if stop_i is not None and i > stop_i:
            break
        text_list.extend(read_fun(txt, show_bar=False))
        if show_bar:
            w.set_description(u'已加载%s个text' % str(i+1))
    return text_list, pathlist

def rjsonl(path, show_bar=False):
    """read jsonlines file

    Args:
        path (str): file path
        show_bar (bool, optional): display read process or not . Defaults to False.

    Returns:
        each time yield a json object

    Examples:
        >>> for item in File.rjsonl("test.jsonl"):
        >>>     print(item)
        {'a': 1}
        {'b': 2}
        {'c': 3}
    """
    with open(path, "r+", encoding="utf8") as f:
        lines = [1 for _ in open(path, "r", encoding="utf-8")]
        w = lines if not show_bar else tqdm(
            lines, total=len(lines), desc=u'已加载0个text')
        reader = jsonlines.Reader(f)
        for i, _ in enumerate(w):
            try:
                item = reader.read()
                yield item
            except Exception:
                traceback.print_exc()
            if show_bar:
                w.set_description(u'已加载%s个text' % str(i+1))
    f.close()

def wjsonl(data, path):
    """write jsonlines file

    Args:
        data (list): a list of json object
        path (str): file path

    Returns:
        None

    Examples:
        >>> File.wjsonl([{'a': 1}, {'b': 2}, {'c': 3}], "test.jsonl")
    """
    with jsonlines.open(path, 'w') as writer:
        writer.write_all(data)

def wjsons(jsons, path, mode="w"):
    """save json object to a file

    Args:
        jsons (list): a list of json object
        path (str): file path
        mode (str, optional): file mode. Defaults to "w".

    Returns:
        None

    Examples:
        >>> File.savejsons([{'a': 1}, {'b': 2}, {'c': 3}], "test.json")
    """
    with open(path, mode) as outfile:
        for entry in jsons:
            json.dump(entry, outfile)
            outfile.write('\n')

def mkdir(dir_path):
    """make directory if dir_path not exist

    Args:
        dir_path (string): directory absolute path

    Returns:
        state: if successfully created directory return True,
                nontheless return False
    Examples:
        >>> File.mkdir("/home/directory/test")
        True
    """
    state = True
    try:
        if not os.path.exists(dir_path):
            os.makedirs(dir_path)
    except Exception:
        state = False
    return state

def save2pkl(obj: dict, path: str):
    """save object to pickle file

    Args:
        obj (dict): object to be saved
        path (str): file path
    Returns:
        state: if successfully created directory return True,
                nontheless return False
    Examples:
        >>> save2pkl({"a": 1}, "test.pkl")
        True
    """
    try:
        mkdir(os.path.dirname(path))
        f = open(path, "wb")
        pickle.dump(obj, f)
        f.close()
        return True
    except Exception:
        traceback.print_exc()
        return False

def files_with_extension(path: str, extension: str):
    """get files by filename extension

    Args:
        path ([type]): [description]
        extension ([type]): [description]

    Returns:
        file_list: file name list
    """
    file_list = [
        file for file in sorted(os.listdir(path))
        if file.lower().endswith(extension)
    ]
    return file_list

def read_csv(filepath: str):
    """read csv from a given filepath

    Args:
        filepath (str): absolute path to the csv file

    Returns:
        rows (list): each element is a row of the csv file
        column_titles (list): each element is a column title of the csv file

    Examples:
        >>> rows, column_titles = File.read_csv("test.csv")
        >>> print(rows)
        [['TOM',20],['JERRY',30],['JANE',25]]
        >>> print(column_titles)
        ['name', 'age']
    """
    # initializing the titles and rows list
    column_titles = []
    rows = []
    num_lines = len([1 for line in open(filepath, "r")])
    # reading csv file
    with open(filepath, 'r') as csvfile:
        # creating a csv reader object
        csvreader = csv.reader((line.replace('\0', '')
                                for line in tqdm(csvfile,
                                                 desc=filepath.split(
                                                     os.sep)[-1]+" is filtering ",
                                                 total=num_lines)))
        # extracting field names through first row
        column_titles = next(csvreader)
        # extracting each data row one by one
        try:
            for row in tqdm(csvreader, desc=filepath.split(os.sep)[-1]
                            + " is loading..", total=num_lines):
                rows.append(row)
        except Exception as e:
            traceback.print_exc()
    return rows, column_titles

def delete(filepath: str):
    """delete file or folders if file exist

    Args:
        file_path (str): absolute path to the file
    """
    if os.path.isdir(filepath):
        shutil.rmtree(filepath)
    elif os.path.isfile(filepath):
        os.remove(filepath)
    else:
        return False

def zip(src_path: str, tgt_path: str):
    """zip files in a directory

    Args:
        zip_path (str): absolute path to the directory
        zip_name (str): zip file name
    """
    date = datetime.now()
    zip_name = src_path.split(os.sep)[-1] + str(date)
    tgt_path = os.path.join(tgt_path, zip_name)
    delete(tgt_path)
    shutil.make_archive(tgt_path, 'zip', src_path)
    return tgt_path

def unzip(src_path, tgt_path, archive_format="zip"):
    return shutil.unpack_archive(src_path, tgt_path, archive_format)

def copy(source_path: str, target_path: str):
    """copy files in a directory

    Args:
        source_path (str): absolute path to the directory
        target_path (str): absolute path to the directory

    Returns:
        state: if successfully copied return True,

    Examples:
        >>> File.copy("/home/directory/test", "/home/directory/test_copy")
        True
    """

    mkdir(target_path)
    if os.path.exists(source_path):
        # root :refer to current iteration directory
        # dirs :refer to all direct sub directories
        # files :refer to all files in current iteration directory
        for root, dirs, files in os.walk(source_path):
            for file in files:
                src_file = os.path.join(root, file)
                shutil.copy(src_file, target_path)
    return True

def save_table2csv(data, path, sortby=None, ascending=False):
    """save a table to csv file

    Args:
        data (dict): key is column name and value is list of values belong to the column
        path (str): save path
        sortby ([type], optional): [description]. Defaults to None.
        ascending (bool, optional): [description]. Defaults to False.

    Example:
        nme = ["aparna", "pankaj", "sudhir", "Geeku"]
        deg = ["MBA", "BCA", "M.Tech", "MBA"]
        scr = [90, 40, 80, 98]
        dict = {'name': nme, 'degree': deg, 'score': scr}
        save_table2csv(dict, os.path.join(sys.path[0], "dict.csv"))
    """
    df = pd.DataFrame(data)
    if sortby is not None:
        df.sort_values([sortby], axis=0, ascending=ascending, inplace=True)
    # saving the dataframe
    df.to_csv(path)

def parent_dir(path: str, layers: int = 1):
    """get parent directory

    Args:
        path(str): file path
        layers(int, optional): layers of parent directory. Defaults to 1.

    Returns:
        path: parent directory
    """
    dirname = os.path.dirname
    for _ in range(layers):
        path = dirname(path)
    return path

def newpth(path, subfolder="", filename: str = ""):
    """get new path

    Args:
        path(str): absolute path
        subfolder(str, optional): sub-folder inside this path . Defaults to "".
        filename(str, optional): filename inside sub-folder. Defaults to "".

    Returns:
        str: absolute path / path/subfolder/filename

    Examples:
        >>> newpth("/home/directory/test", "subfolder", "filename")
        /home/directory/test/subfolder/filename
    """
    return os.path.join(parent_dir(path, 1), subfolder, filename)

def get_all_sub_folders(path):
    """get all sub folders by path

    Args:
        path(str): target path

    Returns:
        list: sub folder list

    Examples:
        >>> get_all_sub_folders("/home/directory/test")
        ['test1', 'test2']
    """
    sub_directory = []
    for filename in os.listdir(path):
        file_path = os.path.join(path, filename)
        if os.path.isdir(file_path):
            sub_directory.append(filename)
    return sub_directory

def generate_cleanfolder(path: str):
    """clean/generate a directory by path recursively

    Args:
        path(str): target path

    Returns:
        state: if successfully created directory return True,

    Examples:
        >>> generate_cleanfolder("/home/directory/test")
        True
    """
    folder_paths = []
    if isinstance(path, list):
        folder_paths.extend(path)
    else:
        folder_paths.append(path)
    for floder in folder_paths:
        # print('make folder:',floder)
        if not os.path.isdir(floder):
            os.mkdir(floder)
    for folder in folder_paths:
        if not os.path.isdir(folder):
            print("this folder name doesn't exist...")
            break
        files = os.listdir(folder)
        for filename in files:
            file_path = os.path.join(folder, filename)
            try:
                if os.path.isfile(file_path):
                    os.remove(file_path)
                elif os.path.islink(file_path):
                    os.unlink(file_path)
                elif os.path.isdir(file_path):
                    # print("shutil :",file_path)
                    shutil.rmtree(file_path)
            except Exception as e:
                print('Failed to delete %s. Reason: %s' % (file_path, e))
                return False
    return True

def generate_emptyfolder_bylist(root_path, folders_list):
    """generate empty folder by list

    Args:
        root_path(str): absolute path
        folders_list(list): folder list

    Examples:
        >>> generate_emptyfolder_bylist("/home/directory/test", ["test1", "test2"])
    """
    floder_list = []
    for fol in folders_list:
        floder_list.append(os.path.join(root_path, fol))
    generate_cleanfolder(floder_list)

def file_size(path):
    """give file absolute path, return file size

    Args:
        path(str): absolute path

    Returns:
        int: file size in bytes
        int: file size in kb
        int: file size in mb

    Examples:
        >>> file_size("/home/directory/test.txt")
        (9, 0.009, 0.0009)
    """
    size = os.path.getsize(path)
    return size, round(size / (1024**2),
                       2), round(size / (1024**2),
                                 2), round(size / (1024**3), 2)

def savelist(obj, path):
    """save list to file

    Args:
        obj(list): a list to be saved
        path(str): file path

    Returns:
        state: if successfully created directory return True

    Examples:
        >>> savelist([1, 2, 3], "test.txt")
        True
    """
    try:
        mkdir(parent_dir(path, 1))
        with open(path, 'w') as f:
            for item in tqdm(obj, desc=path.split(os.sep)[-1] + " is saving..."):
                f.write("%s\n" % item)
        return True
    except Exception as ex:
        traceback.print_exc()
        return False

def append_list2file(path, obj, show_bar=False):
    """append list to file

    Args:
        path(str): absolute path
        obj(list): a list to be saved
        show_bar(bool, optional): [description]. Defaults to False.

    Returns:
        state: if successfully created directory return True

    Examples:
        >>> append_list2file("test.txt", [1, 2, 3])
        True
        >>> append_list2file("test.txt", [4, 5, 6], True)
        True
    """
    mkdir(os.path.dirname(path))
    try:
        with open(path, 'a+') as f:
            if show_bar:
                for item in tqdm(obj, desc=path.split(os.sep)[-1] + " is saving..."):
                    f.write("%s\n" % item)
            else:
                for item in obj:
                    f.write("%s\n" % item)
        return True
    except Exception as ex:
        traceback.print_exc()
        return False

def rjson(path, encoding='utf-8'):
    """read json file

    Args:
        path(str): absolute path
        encoding(str, optional): [description]. Defaults to 'utf-8'.

    Returns:
        dict: json data

    Examples:
        >>> rjson("test.json")
        {'a': 1, 'b': 2}
    """
    f = open(path)
    data = json.load(f, encoding=encoding)
    return data

def rpkl(dict_path):
    """read pickle file

    Args:
        dict_path(str): pickle file path

    Returns:
        dict_object(dict): pickle file object

    Examples:
        >>> rpkl("test.pkl")
        {'a': 1, 'b': 2}
    """
    with open(dict_path, "rb") as f:
        dict_object = pickle.load(f)
    return dict_object

def union_dict(dict1, dict2):
    """combine two dicts

    Args:
        dict1(dict): only allow dict which value is int
        dict2(dict): only allow dict which value is int

    Returns:
        dict2: combined dict
    Examples:
        >>> d = Dict()
        >>> d.union_dict({"a": 1}, {"b": 2})
        {'a': 1, 'b': 2}
    """
    for key in dict1.keys():
        if dict2.get(key) != None:
            dict2[key] += dict1[key]
        else:
            dict2[key] = dict1[key]
    return dict2

def split_dict(dictionary: dict, split_nums: int):
    """split dict into several parts

    Args:
        dictionary(dict): a dict to be split
        split_nums(int): split nums

    Returns:
        list: each element is a dict

    Examples:
        >>> d = Dict()
        >>> d.split_dict({"a": 1, "b": 2, "c": 3}, 2)
        [{'a': 1, 'b': 2}, {'c': 3}]
    """
    dict_lengths = len(dictionary)
    batch_size = dict_lengths // split_nums
    batch_dict = []
    for n in range(split_nums + 1):
        cur_idx = batch_size * n
        end_idx = batch_size * (n + 1)
        cur_batch = dict(list(dictionary.items())[cur_idx: end_idx])
        batch_dict.append(cur_batch)
    return batch_dict

def viw_pkl(path, start=0, end=10):
    """view dict in pickle file from start to end

    Args:
        path(str): absolute path
        start(int, optional): start index of dict. Defaults to 0.
        end(int, optional): end index of dict. Defaults to 10.

    Returns:
        result(dict): a small dict

    Examples:
        >>> d = Dict()
        >>> d.viw_pkl("/home/directory/test.pkl", 0, 10)
        {'a': 1, 'b': 2, 'c': 3, 'd': 4, 'e': 5,
            'f': 6, 'g': 7, 'h': 8, 'i': 9, 'j': 10}
    """
    n_pkl = []
    with open(path, "rb") as f:
        dict_object = pickle.load(f)
        result = dict(list(dict_object.items())[start: end])
    return result

def get_keys(slef, val, obj: dict):
    """get keys by value in dict

    Args:
        val([type]): value
        obj(dict): dict

    Returns:
        list: keys\

    Examples:
        >>> d = Dict()
        >>> d.get_keys(1, {"a": 1, "b": 2})
        ['a']
    """
    return [k for k, v in obj.items() if v == val]

def sort_dict_by_value(d, increase=True):
    """sort dict by value

    Args:
        d(dict): dict to be sorted
        increase(bool, optional): increase sort or decrease sort. Defaults to True.

    Returns:
        [type]: [description]

    Examples:
        >>> d = Dict()
        >>> d.sort_dict_by_value({"a": 1, "b": 2, "c": 3}, increase=False)
        [{'c': 3}, {'b': 2}, {'a': 1}]
    """
    return dict(sorted(d.items(), key=lambda x: x[1], reverse=not increase))

def gen_a_random_probability(self):
    """generate a random probability obey uniform distribution

    Returns:
        float: a random probability

    Examples:
        >>> s = String()
        >>> s.gen_a_random_probability()
        0.5
    """
    return random.uniform(0, 1)

def get_window_content(center, whole_context, windows_size):
    """given whole context, center then return window's content in range windows_size by character

    Args:
        center(str): center word
        whole_context(str): whole context
        windows_size(int): window size

    Returns:
        list: window's content

    Examples:
        >>> s = String()
        >>> s.get_window_content("a", "abcdefg", 3)
        ['a', 'b', 'c']
    """
    neighbors = []
    length = len(whole_context)
    # 窗口左边界下标
    left_boundary = max(center - windows_size // 2, 0)
    # 完整窗口坐标list
    window_idx = (
        list(range(left_boundary, center)) + [center] +
        list(range(center + 1, min(center + windows_size // 2 + 1, length))))
    # 完整窗口内容
    window_content = [whole_context[i] for i in window_idx]

    return "".join(window_content)

def get_current_time(hour=True):
    """get current time

    Args:
        hour(bool, optional): show hours-minutes-seconds or not. Defaults to True.

    Returns:
        str: current time

    Examples:
        >>> s = String()
        >>> s.get_current_time()
        '2019-07-24 15:00:00'
        >>> s.get_current_time(False)
        '2019-07-24'
    """
    now = datetime.now()
    hour = "_%H-%M-%S" if hour else ""
    return now.strftime("%Y-%m-%d" + hour)

def unquote_text(text):
    """parse quoted str to unquoted str

    Args:
        test (str): quoted text

    Returns:
        str: unquoted string

    Examples:
        >>> s = String()
        >>> s.url2str("%E6%8A%98%E5%8F%A0%E8%87%AA%E8%A1%8C%E8%BD%A6")
        '折叠自行车'
    """
    unquoted = unquote(text, encoding="utf-8", errors="replace")
    return unquoted

def str2url(text):
    """convert unquoted str to quoted str

    Args:
        text (str): unquoted str

    Returns:
        str: quoted str

    Examples:
        >>> s = String()
        >>> s.str2url("折叠自行车")
        '%E6%8A%98%E5%8F%A0%E8%87%AA%E8%A1%8C%E8%BD%A6'
    """
    quoted = quote(text, encoding="utf-8", errors="replace")
    return quoted

def remove_blank(text):
    """
    Args:
        text (str): input text, contains blank between zh and en, zh and zh, en and en
    Returns:
        str: text without blank between zh and en, zh and zh, but keep en and en

    Examples:
        >>> text = "比如 Convolutional Neural Network，CNN 对应中 文是卷 积神 经网络。"
        >>> remove_blank(text)
        "比如Convolutional Neural Network，CNN对应中文是卷积神经网络。"
    """
    # filter blank space between Chinese characters
    text = re.sub(r'([^a-zA-Z])([\u0020]*)', r'\1', text)
    # remove blank space between English characters and Chinese characters
    text = re.sub(r'([\u0020]*)([^a-zA-Z])', r'\2', text)
    return text

def iszh(text):
    """check if text is Chinese

    Args:
        text (str): text

    Returns:
        bool: True if text is Chinese, False otherwise

    Examples:
        >>> iszh("比如")
        True
        >>> iszh("比如Convolutional Neural Network，CNN对应中文是卷积神经网络。")
        False
    """
    for ch in text:
        if not "\u4e00" <= ch <= "\u9fff":
            return False
    return True

def colourful_text(text, color):
    """add color to text

    Args:
        text (str): [description]
        color (str): red, green, yellow, blue, black, none

    Returns:
        str: coloured text

    Examples:
        >>> s = String()
        >>> s.colourful_text("比如Convolutional Neural Network，CNN对应中文是卷积神经网络。", "red")
        '\x1b[31m比如Convolutional Neural Network，CNN对应中文是卷积神经网络。\x1b[0m'
    """
    colourful = {
        "red": u"\033[1;31;1m%s\033[0m",
        "green": u"\033[1;32;1m%s\033[0m",
        "yellow": u"\033[1;33;1m%s\033[0m",
        "blue": u"\033[1;34;1m%s\033[0m",
        "black": u"\033[1;30;1m%s\033[0m",
    }
    return colourful[color] % text if color != "none" else text

def type_restore(value):
    """restore str go back to original type

    Args:
        value (str): input str

    Returns:
        value (str): original type

    Examples:
        >>> s = String()
        >>> s.type_restore("1")
        1
        >>> s.type_restore("1.0")
        1.0
        >>> s.type_restore("True")
        True
        >>> s.type_restore("False")
        False
        >>> s.type_restore("None")
        None
        >>> s.type_restore("[1, 2, 3]")
        [1, 2, 3]
        >>> s.type_restore("{'a': 1, 'b': 2}")
        {'a': 1, 'b': 2}
    """
    try:
        value = ast.literal_eval(value)
    except Exception:
        value = value
    return value

def extract(text, expression=None, type: str = None):
    """extract text from text by regular expression

    Args:
        text (str): input text
        expression (str, optional): regular expression. Defaults to None.
        type (str, optional): type of extracted text. Defaults to None.

    Returns:
        str: extracted text

    Examples:
        >>> s = String()
        >>> text ="我今天很happy，因为我喜欢的notebook终于卖完了。我是在刚刚发布的时候以12900元的价格买的，现在以12800元的价格卖出了。"
        >>> s.extract(text, type="en")
        (['happy', 'notebook'], True)
        >>> s.extract(text, type="zh")
        (['我今天很', '因为我喜欢的', '终于卖完了', '我是在刚刚发布的时候以', '元的价格买的', '现在以', '元的价格卖出了'], True)
        >>> s.extract(text, type="zh", expression="num")
        (['12900', '12800'], True)
        >>> s.extract(text, type="zh", expression="punctuation")
        (['，', '。', '，', '。'], True)
    """
    if type is not None:
        type2expression = {"en": r"[a-zA-Z]+", "zh": r"[\u4e00-\u9fa5]+", "num": r"\d+",
                           "punctuation": u"[\u3002\uff1b\uff0c\uff1a\u201c\u201d\uff08\uff09\u3001\uff1f\u300a\u300b]"}
        expression = type2expression[type]
    try:
        pattern = re.compile(expression)
        results = re.findall(pattern, text)
        return results, True if len(results) != 0 else False
    except Exception:
        traceback.print_exc()

def half_to_full(text):
    """convert half-width to full-width

    Args:
        text (str): half-width text

    Returns:
        change_text (str): full-width text

    Examples:
        >>> s = String()
        >>> s.half_to_full("abc")
        'ａｂｃ'
    """
    change_text = ""
    for word in text:
        inside_code = ord(word)
        if inside_code == 32:  # 半角空格直接转化
            inside_code = 12288
        elif inside_code >= 32 and inside_code <= 126:  # 半角字符（除空格）根据关系转化
            inside_code += 65248
        change_text += chr(inside_code)
    return change_text

def full_to_half(text: str):  # 输入为一个句子
    """convert full-width to half-width

    Args:
        text (str): full-width text

    Returns:
        change_text (str): half-width text
    """
    change_text = ""
    for word in text:
        inside_code = ord(word)
        # print("before:", inside_code)
        if inside_code == 12288:  # 全角空格直接转换
            inside_code = 32
        # elif inside_code in [ord("“"), ord("”")]:
        #     inside_code = ord('"')
        # elif inside_code == ord("。"):
        #     inside_code = ord(".")
        # elif inside_code in [ord("‘"), ord("’")]:
        #     inside_code = ord("'")
        elif inside_code >= 65281 and inside_code <= 65374:  # 全角字符（除空格）根据关系转化
            inside_code -= 65248
        change_text += chr(inside_code)
        # print("after:", inside_code)
    return change_text

def lc_subsequence(str1: str, str2: str, split_char: str = " "):
    """longest common subsequence

    Args:
        str1 (str): str1
        str2 (str): str2
        split_char (str): split char, when character isn't joined, use split_char to split

    Returns:
        str: longest common subsequence

    Examples:
        >>> s = String()
        >>> s.lc_subsequence("我今天吃了碗饭", "我昨天也吃了一碗饭", "_")
        我_天_吃了_碗饭
        >>> s.lc_subsequence("我今天吃了碗饭", "我昨天也吃了一碗饭", "")
        我天吃了碗饭
    """
    m, n = len(str1), len(str2)
    L = [[0 for x in range(n + 1)] for x in range(m + 1)]
    for i in range(m + 1):
        for j in range(n + 1):
            if i == 0 or j == 0:
                L[i][j] = 0
            elif str1[i - 1] == str2[j - 1]:
                L[i][j] = L[i - 1][j - 1] + 1
            else:
                L[i][j] = max(L[i - 1][j], L[i][j - 1])
    index = L[m][n]
    lcs = [""] * (index + 1)
    lcs[index] = ""
    i = m
    j = n
    join = True
    while i > 0 and j > 0:
        if str1[i - 1] == str2[j - 1]:
            split_code = ""
            if not join:
                split_code = split_char
            lcs[index - 1] = str1[i - 1] + split_code
            i -= 1
            j -= 1
            index -= 1
            join = True
        elif L[i - 1][j] > L[i][j - 1]:
            i -= 1
            join = False
        else:
            j -= 1
            join = False
    return "".join(lcs)

def lc_substring(str1, str2):
    """longest common substring

    Args:
        str1 (str): str1
        str2 (str): str2

    Returns:
        str: longest common substring
        int: length of longest common substring

    Examples:
        >>> s = String()
        >>> s.lc_substring("我今天是吃了碗饭", "我昨天也是吃了一碗饭")
        ('是吃了', 3)
    """
    m = [[0 for i in range(len(str2) + 1)]
         for j in range(len(str1) + 1)]
    mmax = 0  # 最长匹配的长度
    p = 0  # 最长匹配对应在s1中的最后一位
    for i in range(len(str1)):
        for j in range(len(str2)):
            if str1[i] == str2[j]:
                m[i + 1][j + 1] = m[i][j] + 1
                if m[i + 1][j + 1] > mmax:
                    mmax = m[i + 1][j + 1]
                    p = i + 1
    return str1[p - mmax:p], mmax  # 返回最长子串及其长

def index_substring(text, sub_text):
    """given a text, find all the sub-text's index

    Args:
        text (str): text
        sub_text (str): sub-text

    Yields:
        int: sub-text's index

    Examples:
        >>> s = String()
        >>> sub_idxs =s.index_substring("我今天吃了碗饭，我昨天也吃了一碗饭", "碗饭")
        >>> [_ for _ in sub_idxs]
        [5, 15]
    """
    start = 0
    while True:
        start = text.find(sub_text, start)
        if start == -1:
            return
        yield start
        # use start += 1 to find overlapping matches
        start += len(sub_text)

def extact_content_between_marks(text, expression: str = None, type: str = None):
    """extract content in bracket

    Args:
        text (str): text
        expression (str): expression

    Returns:
        str: content in bracket

    Examples:
        >>> s = String()
        >>> text = "如果说【上天】愿意再次给我（一次机会），那么我(希望)保持现状"
        >>> s.extact_content_in_bracket(text=text)
        ['【上天】', '（一次机会）', '(希望)']
    """
    if type is not None:
        type2expression = {
            "bracket": r"(\[|（|\(|【|《).+?(）|\]|\)|】|》)", "quotation": r"(\"|\'|‘|“).+?(’|”|\'|\")"}
        expression = type2expression[type]
    elif expression is None:
        raise ValueError("expression is None")
    pattern = re.compile(expression)
    matchs = re.finditer(pattern, text)
    results = []
    for match in matchs:
        results.append(match.group())
    return results if len(results) != 0 else None

def split_text(text, pattern=r";|。|；|,|，"):
    """split text by pattern

    Args:
        text (str): text
        pattern (regexp, optional): expression. Defaults to r";|。|；|,|，".

    Returns:
        str: text split by pattern

    Examples:
        >>> s = String()
        >>> text = "收快递的时候最怕收不到货，所以购物的时候一定要把地址写清楚，这样才会精准的送到你手里，我告诉大家以后怎么写：“本宇宙-拉尼凯亚超星系团-室女座星系团-本星系群-银河系-猎户臂-太阳系-第三行星-地球-亚洲板块-中国-xxx-xxx-xxx”这样可以保证不会送到其他宇宙去"
        >>> s.split_text(text=text)
        ['收快递的时候最怕收不到货', '所以购物的时候一定要把地址写清楚', '这样才会精准的送到你手里', '我告诉大家以后怎么写：“本宇宙-拉尼凯亚超星系团-室女座星系团-本星系群-银河系-猎户臂-太阳系-第三行星-地球-亚洲板块-中国-xxx-xxx-xxx”这样可以保证不会送到其他宇宙去']
    """
    txts = re.split(pattern, text)
    return txts

def similar(str1, str2):
    """calculate similarity between two strings

    Args:
        str1 (str): str1
        str2 (str): str2

    Returns:
        float: similarity

    Examples:
        >>> s = String()
        >>> s.similar("我今天是吃了碗饭", "我昨天也是吃了一碗饭")
        0.75
    """
    return SequenceMatcher(None, str1, str2).ratio()

def str_num_type(string):
    """check string whether belong to int,float,complex or not

    Args:
        string (str): string

    Returns:
        str: int,float,complex or None

    Examples:
        >>> s = String()
        >>> s.str_num_type("1")
        'int'
        >>> s.str_num_type("1.0")
        'float'
        >>> s.str_num_type("1+1j")
        'complex'
        >>> s.str_num_type("1.0.5.355")
        None

    """
    types = [int, float, complex]
    _string = type_restore(string)

    for type in types:
        if isinstance(_string, type):
            return str(type.__name__)
    return None

def proper_nouns_alignment(original_text, changed_text):
    """this function align the content which surround by bracket or
    quotation marks. Prevent some proper nouns changed by the reducer system.

    Args:
        original_text (string): original text
        changed_text (string): changed text

    Returns:
        [string]: alignment text
    """

    expression_types = ["bracket", "quotation"]
    for regular_exp in expression_types:
        result1 = extact_content_between_marks(
            original_text, type=regular_exp)
        result2 = extact_content_between_marks(
            changed_text, type=regular_exp)
        if len(result1) != len(result2):
            return changed_text, False
        for idx, ele in enumerate(result2):
            origin_ele = result1[idx]  # 按序匹配
            changed_text = changed_text.replace(ele, origin_ele)
    return changed_text, True

def detect_repeat(text, expression=r"([\u4e00-\u9fa5])(\1+)"):
    """detect repeat characters

    Args:
        text (str): text
        expression (regexp, optional): regular expression . Defaults to r"([\u4e00-\u9fa5])(\1+)".

    Returns:
        (list): repeat characters

    Examples:
        >>> s = String()
        >>> text = "锄禾日当当午，汗滴禾下下土"
        >>> s.detect_repeat(text)
        ['当当', '下下']
    """
    pattern = re.compile(expression)
    matchs = re.finditer(pattern, text)
    results = []
    for match in matchs:
        results.append(match.group())
    return results

if __name__ == '__main__':
    lines = readlines(path="/mnt/f/git_repository/notebook/pages/URI.md")
    print(lines)

max_stars_repo_name: TITC/nlptools

max_issues_repo_name: TITC/nlptools

max_forks_repo_name: TITC/nlptools

reason: from numpy


text:
import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.optimize import curve_fit

#pdf for particle tracking
def p(x,mu,norm):
    f=norm*mu*np.exp(-mu*x)
    return f

#inverse cumulative pdf
def c_inv(x,mu):
    c=-(np.log(x))/(mu)
    return c

N=10**6             #number of entries for Monte Carlo
n_bins=5*10**2           #number of bins

# compute the individual step lengths for process 1 and 2, keep the minimum, count how many times  s1 is the minimum
#Geant4-like sampling
s=[]
counter=0
for i in range(0,N):
    x=random.random()
    y=random.random()
    s_1=c_inv(x,1)
    s_2=c_inv(y,2)
    if s_1 < s_2:
        s.append(s_1)
        counter+=1
    else:
        s.append(s_2)
s=np.array(s)

#fit the histogram with the exponential function used by Penelope for sampling
hist, bin_edges=np.histogram(s, bins=n_bins, density=False)
x=np.linspace(s.min(),s.max(),n_bins)
fit, cov =curve_fit(p, xdata=x, ydata=hist, p0=(3,10000))
mu_fit, norm_fit = fit
mu_err, norm_err =  np.sqrt(np.diag(cov)[0]) , np.sqrt(np.diag(cov)[1])

#plot the results
plt.figure()
plt.hist(s, color='orange', bins=n_bins, histtype='step')
plt.plot(x, p(x,mu_fit, norm_fit), 'k')
plt.xlabel('min($s_1, s_2$)', fontsize=14)
plt.ylabel('$Counts$', fontsize=14)
plt.yscale('log')
plt.xlim(0, 3, 0.5)
plt.xticks(np.arange(0,3.5, 0.5), fontsize=10)
plt.yticks(fontsize=10)
props = dict(boxstyle='round', facecolor='orange', alpha=0.5)
plt.text(2, 10000, '$\mu = %.3f \pm %.3f $' %(mu_fit, mu_err), fontsize=10, bbox=props)

print('probability that s_1 is the minimum: %.2f \n probability that s_2 is the minimum: %.2f' %(counter/N, (N-counter)/N))

ext: py

reason: import numpy,from scipy


text:
import sympy as sp
x, y, z = sp.symbols('x y z')
sp.init_printing()

# Here we simply give the input needed as The Following:
# ( Xi ) is The Initial Point
# ( Fx ) is the equation of the function
# ( n ) is the number of Iterations needed

Xi = 0

Fx = x**3 -0.2589*x**2+0.2262*x - 0.001122

n = 5

print("\n_____________________________________________________")
print(" Iteration |   Xi    |   F(Xi)  |   F(x)\'  |  Error "
      "\n-----------------------------------------------------")
Error = 100 ; XiPlus1 = 0
for i in range (n):
    FxVal = float(Fx.subs(x, Xi).evalf())
    drFx = float(sp.diff(Fx, x).doit().subs(x, Xi).evalf())
    print("\t{0:1d}  | {1:.5f} | {2:.4f}  | {3:.4f}  | {4:.4f}%\n"
          .format(i+1, Xi,FxVal ,drFx , Error)) ; XiPlus1 = Xi - FxVal/drFx
    Error = abs(((XiPlus1 - Xi)/XiPlus1) *100) ; Xi = XiPlus1

# In The End The Program Prints The Output as well as The Percentage of error between Iterations

reason: import sympy


text:
## @ingroup Methods-Geometry-Two_Dimensional-Cross_Section-Airfoil
# compute_airfoil_polars.py
#
# Created:  Mar 2019, M. Clarke
# Modified: Mar 2020, M. Clarke
#           Jan 2021, E. Botero
#           Jan 2021, R. Erhard
#           Nov 2021, R. Erhard

# ----------------------------------------------------------------------
#  Imports
# ----------------------------------------------------------------------

import SUAVE
from SUAVE.Core               import Data , Units
from SUAVE.Methods.Aerodynamics.AERODAS.pre_stall_coefficients import pre_stall_coefficients
from SUAVE.Methods.Aerodynamics.AERODAS.post_stall_coefficients import post_stall_coefficients
from .import_airfoil_geometry import import_airfoil_geometry
from .import_airfoil_polars   import import_airfoil_polars
from scipy.interpolate        import RectBivariateSpline
import numpy as np

## @ingroup Methods-Geometry-Two_Dimensional-Cross_Section-Airfoil
def compute_airfoil_polars(a_geo,a_polar,use_pre_stall_data=True):
    """This computes the lift and drag coefficients of an airfoil in stall regimes using pre-stall
    characterstics and AERODAS formation for post stall characteristics. This is useful for
    obtaining a more accurate prediction of wing and blade loading. Pre stall characteristics
    are obtained in the form of a text file of airfoil polar data obtained from airfoiltools.com

    Assumptions:
    Uses AERODAS formulation for post stall characteristics

    Source:
    Models of Lift and Drag Coefficients of Stalled and Unstalled Airfoils in Wind Turbines and Wind Tunnels
    by D Spera, 2008

    Inputs:
    a_geo                  <string>
    a_polar                <string>
    use_pre_stall_data     [Boolean]

    Outputs:
    airfoil_data.
        cl_polars          [unitless]
        cd_polars          [unitless]
        aoa_sweep          [unitless]

    Properties Used:
    N/A
    """

    num_airfoils = len(a_geo)

    # check number of polars per airfoil in batch
    num_polars   = 0
    for i in range(num_airfoils):
        n_p = len(a_polar[i])
        if n_p < 3:
            raise AttributeError('Provide three or more airfoil polars to compute surrogate')

        num_polars = max(num_polars, n_p)

    # read airfoil geometry
    airfoil_data = import_airfoil_geometry(a_geo)

    # Get all of the coefficients for AERODAS wings
    AoA_sweep_deg = np.linspace(-14,90,105)
    AoA_sweep_radians = AoA_sweep_deg*Units.degrees
    CL = np.zeros((num_airfoils,num_polars,len(AoA_sweep_deg)))
    CD = np.zeros((num_airfoils,num_polars,len(AoA_sweep_deg)))
    aoa0 = np.zeros((num_airfoils,num_polars))
    cl0 = np.zeros((num_airfoils,num_polars))

    CL_surs = Data()
    CD_surs = Data()
    aoa_from_polars = []

    # Create an infinite aspect ratio wing
    geometry              = SUAVE.Components.Wings.Wing()
    geometry.aspect_ratio = np.inf
    geometry.section      = Data()

    # Create dummy settings and state
    settings = Data()
    state    = Data()
    state.conditions = Data()
    state.conditions.aerodynamics = Data()
    state.conditions.aerodynamics.pre_stall_coefficients = Data()
    state.conditions.aerodynamics.post_stall_coefficients = Data()

    # read airfoil polars
    airfoil_polar_data = import_airfoil_polars(a_polar)

    # AERODAS
    for i in range(num_airfoils):

        # Modify the "wing" slightly:
        geometry.thickness_to_chord = airfoil_data.thickness_to_chord[i]

        for j in range(len(a_polar[i])):
            # Extract from polars
            airfoil_cl         = airfoil_polar_data.lift_coefficients[i,j]
            airfoil_cd         = airfoil_polar_data.drag_coefficients[i,j]
            airfoil_aoa        = airfoil_polar_data.angle_of_attacks

            # computing approximate zero lift aoa
            airfoil_cl_plus = airfoil_cl[airfoil_cl>0]
            idx_zero_lift = np.where(airfoil_cl == min(airfoil_cl_plus))[0][0]
            airfoil_cl_crossing = airfoil_cl[idx_zero_lift-1:idx_zero_lift+1]
            airfoil_aoa_crossing = airfoil_aoa[idx_zero_lift-1:idx_zero_lift+1]
            try:
                A0  = np.interp(0,airfoil_cl_crossing, airfoil_aoa_crossing)* Units.deg
            except:
                A0 = airfoil_aoa[idx_zero_lift] * Units.deg

            # max lift coefficent and associated aoa
            CL1max = np.max(airfoil_cl)
            idx_aoa_max_prestall_cl = np.where(airfoil_cl == CL1max)[0][0]
            ACL1 = airfoil_aoa[idx_aoa_max_prestall_cl] * Units.degrees

            # computing approximate lift curve slope
            linear_idxs = [int(np.where(airfoil_aoa==0)[0]),int(np.where(airfoil_aoa==4)[0])]
            cl_range = airfoil_cl[linear_idxs]
            aoa_range = airfoil_aoa[linear_idxs] * Units.degrees
            S1 = (cl_range[1]-cl_range[0])/(aoa_range[1]-aoa_range[0])

            # max drag coefficent and associated aoa
            CD1max  = np.max(airfoil_cd)
            idx_aoa_max_prestall_cd = np.where(airfoil_cd == CD1max)[0][0]
            ACD1   = airfoil_aoa[idx_aoa_max_prestall_cd] * Units.degrees

            # Find the point of lowest drag and the CD
            idx_CD_min = np.where(airfoil_cd==min(airfoil_cd))[0][0]
            ACDmin     = airfoil_aoa[idx_CD_min] * Units.degrees
            CDmin      = airfoil_cd[idx_CD_min]

            # Setup data structures for this run
            ones = np.ones_like(AoA_sweep_radians)
            settings.section_zero_lift_angle_of_attack                = A0
            state.conditions.aerodynamics.angle_of_attack             = AoA_sweep_radians * ones
            geometry.section.angle_attack_max_prestall_lift           = ACL1 * ones
            geometry.pre_stall_maximum_drag_coefficient_angle         = ACD1 * ones
            geometry.pre_stall_maximum_lift_coefficient               = CL1max * ones
            geometry.pre_stall_maximum_lift_drag_coefficient          = CD1max * ones
            geometry.section.minimum_drag_coefficient                 = CDmin * ones
            geometry.section.minimum_drag_coefficient_angle_of_attack = ACDmin
            geometry.pre_stall_lift_curve_slope                       = S1

            # Get prestall coefficients
            CL1, CD1 = pre_stall_coefficients(state,settings,geometry)

            # Get poststall coefficents
            CL2, CD2 = post_stall_coefficients(state,settings,geometry)

            # Take the maxes
            CL_ij = np.fmax(CL1,CL2)
            CL_ij[AoA_sweep_radians<=A0] = np.fmin(CL1[AoA_sweep_radians<=A0],CL2[AoA_sweep_radians<=A0])

            CD_ij = np.fmax(CD1,CD2)

            # Pack this loop
            CL[i,j,:] = CL_ij
            CD[i,j,:] = CD_ij
            aoa0[i,j] = A0
            cl0[i,j]  = np.interp(0,airfoil_aoa,airfoil_cl)

            if use_pre_stall_data == True:
                CL[i,j,:], CD[i,j,:] = apply_pre_stall_data(AoA_sweep_deg, airfoil_aoa, airfoil_cl, airfoil_cd, CL[i,j,:], CD[i,j,:])

        # remove placeholder values (for airfoils that have different number of polars)
        n_p = len(a_polar[i])
        CL_sur = RectBivariateSpline(airfoil_polar_data.reynolds_number[i][0:n_p],AoA_sweep_radians, CL[i,0:n_p,:])
        CD_sur = RectBivariateSpline(airfoil_polar_data.reynolds_number[i][0:n_p],AoA_sweep_radians, CD[i,0:n_p,:])

        CL_surs[a_geo[i]]  = CL_sur
        CD_surs[a_geo[i]]  = CD_sur
        aoa_from_polars.append(airfoil_aoa)

    airfoil_data.angle_of_attacks              = AoA_sweep_radians
    airfoil_data.lift_coefficient_surrogates   = CL_surs
    airfoil_data.drag_coefficient_surrogates   = CD_surs

    airfoil_data.lift_coefficients_from_polar  = airfoil_polar_data.lift_coefficients
    airfoil_data.drag_coefficients_from_polar  = airfoil_polar_data.drag_coefficients
    airfoil_data.re_from_polar  = airfoil_polar_data.reynolds_number
    airfoil_data.aoa_from_polar = aoa_from_polars #airfoil_polar_data.angle_of_attacks

    return airfoil_data

def apply_pre_stall_data(AoA_sweep_deg, airfoil_aoa, airfoil_cl, airfoil_cd, CL, CD):
    # Coefficients in pre-stall regime taken from experimental data:
    aoa_locs = (AoA_sweep_deg>=airfoil_aoa[0]) * (AoA_sweep_deg<=airfoil_aoa[-1])
    aoa_in_data = AoA_sweep_deg[aoa_locs]

    # if the data is within experimental use it, if not keep the surrogate values
    CL[aoa_locs] = airfoil_cl[abs(aoa_in_data[:,None] - airfoil_aoa[None,:]).argmin(axis=-1)]
    CD[aoa_locs] = airfoil_cd[abs(aoa_in_data[:,None] - airfoil_aoa[None,:]).argmin(axis=-1)]

    # remove kinks/overlap between pre- and post-stall
    data_lb = np.where(CD == airfoil_cd[0])[0][0]
    data_ub = np.where(CD == airfoil_cd[-1])[0][-1]
    CD[0:data_lb] = np.maximum(CD[0:data_lb], CD[data_lb]*np.ones_like(CD[0:data_lb]))
    CD[data_ub:]  = np.maximum(CD[data_ub:],  CD[data_ub]*np.ones_like(CD[data_ub:]))

    return CL, CD


text:
import taichi as ti
import numpy as np
ti.init()

N = 128

_et = np.array(
    [
        [[-1, -1], [-1, -1]],  #
        [[0, 1], [-1, -1]],  #a
        [[0, 2], [-1, -1]],  #b
        [[1, 2], [-1, -1]],  #ab
        [[1, 3], [-1, -1]],  #c
        [[0, 3], [-1, -1]],  #ca
        [[2, 3], [0, 1]],  #cb
        [[2, 3], [-1, -1]],  #cab
        [[2, 3], [-1, -1]],  #d
        [[2, 3], [0, 1]],  #da
        [[0, 3], [-1, -1]],  #db
        [[1, 3], [-1, -1]],  #dab
        [[1, 2], [-1, -1]],  #dc
        [[0, 2], [-1, -1]],  #dca
        [[0, 1], [-1, -1]],  #dcb
        [[-1, -1], [-1, -1]],  #dcab
    ],
    np.int32)

m = ti.field(float, (N, N))  # field value of voxels

r = ti.Vector.field(2, float, (N**2, 2))  # generated edges
et = ti.Vector.field(2, int, _et.shape[:2])  # edge table (constant)
et.from_numpy(_et)

@ti.func
def gauss(x):
    return ti.exp(-6 * x**2)

@ti.kernel
def touch(mx: float, my: float):
    for i, j in m:
        p = ti.Vector([i, j]) / N
        a = gauss((p - 0.5).norm() / 0.25)
        p.x -= mx - 0.5 / N
        p.y -= my - 0.5 / N
        b = gauss(p.norm() / 0.25)
        m[i, j] = (a + b) * 3

@ti.func
def list_subscript(a,
                   i):  # magic method to subscript a list with dynamic index
    ret = sum(a) * 0
    for j in ti.static(range(len(a))):
        if i == j:
            ret = a[j]
    return ret

@ti.kernel
def march() -> int:
    r_n = 0

    for i, j in ti.ndrange(N - 1, N - 1):
        id = 0
        if m[i, j] > 1: id |= 1
        if m[i + 1, j] > 1: id |= 2
        if m[i, j + 1] > 1: id |= 4
        if m[i + 1, j + 1] > 1: id |= 8

        E = [ti.Vector(_) + .5 for _ in [(.5, 0), (0, .5), (1, .5), (.5, 1)]]

        for k in range(2):
            if et[id, k][0] == -1:
                break

            n = ti.atomic_add(r_n, 1)
            for l in ti.static(range(2)):
                e = et[id, k][l]
                r[n, l] = ti.Vector([i, j]) + list_subscript(E, e)

    return r_n

gui = ti.GUI('Marching rect')
while gui.running and not gui.get_event(gui.ESCAPE):
    touch(*gui.get_cursor_pos())
    ret_len = march()
    ret = r.to_numpy()[:ret_len] / N
    gui.set_image(ti.imresize(m, *gui.res))
    gui.lines(ret[:, 0], ret[:, 1], color=0xff66cc, radius=1.5)
    gui.show()

max_stars_repo_name: GuangSunshine/taichi

max_issues_repo_name: zxgtz/taichi

max_forks_repo_name: zxgtz/taichi

reason: import numpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 21 11:59:11 2020

@author: seanselzer
"""

import pandas as pd
import numpy as np
import copy
import time
from pathlib import Path
from xs_calculator import xs_calculator

def create_summary_dataframe(lvl2_main_df, log_df, sample_variable, sub_sample_variable, output_directory, linear_data_type):

    half_life223 = 11.43 #days
    half_life224 = 3.65 #days
    half_life_th228 = 1.9*365 #days
    half_life_ra228 = 5.7*365 #days

    ra223_lambda_days = (np.log(2)/(half_life223))
    ra224_lambda_days = (np.log(2)/(half_life224))
    th228_lambda_days = (np.log(2)/(half_life_th228))
    ra228_lambda_days = (np.log(2)/(half_life_ra228))

    summary_df = copy.deepcopy(log_df)
    if linear_data_type == True:
        None_list = ['None' for i in range (len( summary_df))]
        summary_df['subsample_dummy_column'] = None_list

    xs224_list = []
    xs224_err_list = []
    xs223_list = []
    xs223_err_list = []
    xs223_thstdonly_list = []
    xs223_thstdonly_err_list = []
    xs224_t0_list = []
    xs223_t0_list = []
    xs223_thstdonly_t0_list = []
    ra228_list = []
    ra228_err_list = []
    th228_list = []
    th228_err_list = []
    ac227_list = []
    ac227_err_list = []
    ac227_thstdonly_list = []
    ac227_thstdonly_err_list = []
    fraction_decayed_224_list = []
    fraction_decayed_223_list = []
    fraction_decayed_223_thstdonly_list = []
    error_list = []
    read_error_list = []
    ra226_list = []
    ra226_err_list = []

#   Index list for xs_calculator ##################################################################################################################
    xs = 0
    xs_err = 1
    fraction_decayed = 2
    xs_t0 = 3
    parent_activity = 4
    parent_activity_err = 5
    row_errors = 6
    read_errors = 7

    error_flag = -999

    for index, row in summary_df.iterrows():
        row_specific_errors = []
        read_specific_errors = {}

        row_sample_variable = row[sample_variable]
        if linear_data_type == False:
            row_sub_sample_variable = str(row[sub_sample_variable])
            read_number_set = set(lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable) &
                                 (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)), 'read_number'])
        else:
            row_sub_sample_variable = str(row['subsample_dummy_column'])
            read_number_set = set(lvl2_main_df.loc[(lvl2_main_df[sample_variable]==row_sample_variable), 'read_number'])

#        print (row_sample_variable, row_sub_sample_variable)

###########################################################################################################################################
#       224xs calculations
###########################################################################################################################################

#####   Read1 - Read3    ##################################################################################################################

        if 1 in read_number_set and 3 in read_number_set:

            xs_calc_results = xs_calculator (lvl2_main_df, sample_variable, sub_sample_variable, row_sample_variable, row_sub_sample_variable,
                   isotope = 'Ra-224', isotope_column_string = 'vdpm224 (dpm/m^3)', isotope_column_string_err = 'vdpm224_err (dpm/m^3)',
                   isotope_lambda_days = ra224_lambda_days,
                   read_a = 1, read_b = 3,
                   read_number_set = read_number_set)

            xs224_list.append(xs_calc_results[xs])
            xs224_err_list.append(xs_calc_results[xs_err])
            fraction_decayed_224_list.append(xs_calc_results[fraction_decayed])
            xs224_t0_list.append(xs_calc_results[xs_t0])
            th228_list.append(xs_calc_results[parent_activity])
            th228_err_list.append(xs_calc_results[parent_activity_err])
            row_specific_errors.append(xs_calc_results[row_errors])
            read_specific_errors.update(xs_calc_results[read_errors])

#####   Read1 - Read4    ##################################################################################################################
        elif 1 in read_number_set and 4 in read_number_set:
            row_specific_errors.append('(224xs, 228Th) using read 4 instead of 3')

            xs_calc_results = xs_calculator (lvl2_main_df, sample_variable, sub_sample_variable, row_sample_variable, row_sub_sample_variable,
                   isotope = 'Ra-224', isotope_column_string = 'vdpm224 (dpm/m^3)', isotope_column_string_err = 'vdpm224_err (dpm/m^3)',
                   isotope_lambda_days = ra224_lambda_days,
                   read_a = 1, read_b = 4,
                   read_number_set = read_number_set)

            xs224_list.append(xs_calc_results[xs])
            xs224_err_list.append(xs_calc_results[xs_err])
            fraction_decayed_224_list.append(xs_calc_results[fraction_decayed])
            xs224_t0_list.append(xs_calc_results[xs_t0])
            th228_list.append(xs_calc_results[parent_activity])
            th228_err_list.append(xs_calc_results[parent_activity_err])
            row_specific_errors.append(xs_calc_results[row_errors])
            read_specific_errors.update(xs_calc_results[read_errors])

#####   Reads missing      ##################################################################################################################
        else:
            row_specific_errors.append('(224xs, 228Th) Required reads not available')
            read_specific_errors.update({'(224xs, 228Th)': 'Required reads not available'})
            xs224_list.append(error_flag)
            xs224_err_list.append(error_flag)
            xs224_t0_list.append(error_flag)
            th228_list.append(error_flag)
            th228_err_list.append(error_flag)
            fraction_decayed_224_list.append(error_flag)

###########################################################################################################################################
#       223xs calculations
###########################################################################################################################################

#####   Read2 - Read4   ##################################################################################################################
        if 2 in read_number_set and 4 in read_number_set:

            xs_calc_results = xs_calculator (lvl2_main_df, sample_variable, sub_sample_variable, row_sample_variable, row_sub_sample_variable,
                   isotope = 'Ra-223', isotope_column_string = 'vdpm223 (dpm/m^3)', isotope_column_string_err = 'vdpm223_err (dpm/m^3)',
                   isotope_lambda_days = ra223_lambda_days,
                   read_a = 2, read_b = 4,
                   read_number_set = read_number_set)

            xs223_list.append(xs_calc_results[xs])
            xs223_err_list.append(xs_calc_results[xs_err])
            fraction_decayed_223_list.append(xs_calc_results[fraction_decayed])
            xs223_t0_list.append(xs_calc_results[xs_t0])
            ac227_list.append(xs_calc_results[parent_activity])
            ac227_err_list.append(xs_calc_results[parent_activity_err])
            row_specific_errors.append(xs_calc_results[row_errors])
            read_specific_errors.update(xs_calc_results[read_errors])

#####   Read1 - Read4      #################################################################################################################
        elif 1 in read_number_set and 4 in read_number_set:
#            xs223_t0_list.append('(223xs) using read 1 instead of 2')
            row_specific_errors.append('(223xs) using read 1 instead of 2')
            xs_calc_results = xs_calculator (lvl2_main_df, sample_variable, sub_sample_variable, row_sample_variable, row_sub_sample_variable,
                   isotope = 'Ra-223', isotope_column_string = 'vdpm223 (dpm/m^3)', isotope_column_string_err = 'vdpm223_err (dpm/m^3)',
                   isotope_lambda_days = ra223_lambda_days,
                   read_a = 1, read_b = 4,
                   read_number_set = read_number_set)

            xs223_list.append(xs_calc_results[xs])
            xs223_err_list.append(xs_calc_results[xs_err])
            fraction_decayed_223_list.append(xs_calc_results[fraction_decayed])
            xs223_t0_list.append(xs_calc_results[xs_t0])
            ac227_list.append(xs_calc_results[parent_activity])
            ac227_err_list.append(xs_calc_results[parent_activity_err])
            row_specific_errors.append(xs_calc_results[row_errors])
            read_specific_errors.update(xs_calc_results[read_errors])

#####   Reads missing       ##################################################################################################################
        else:
            row_specific_errors.append('(223xs,227Ac) Required reads not available')
            read_specific_errors.update({'(223xs,227Ac)': 'Required reads not available'})
            ac227_list.append(error_flag)
            ac227_err_list.append(error_flag)
            xs223_list.append(error_flag)
            xs223_err_list.append(error_flag)
            xs223_t0_list.append(error_flag)
            fraction_decayed_223_list.append(error_flag)

###########################################################################################################################################
#       223xs_from_thstd calculations
###########################################################################################################################################

#####   Read2 - Read4   ##################################################################################################################
        if 2 in read_number_set and 4 in read_number_set:

            xs_calc_results = xs_calculator (lvl2_main_df, sample_variable, sub_sample_variable, row_sample_variable, row_sub_sample_variable,
                   isotope = 'Ra-223', isotope_column_string = 'vdpm223_thstdonly (dpm/m^3)', isotope_column_string_err = 'vdpm223_thstdonly_err (dpm/m^3)',
                   isotope_lambda_days = ra223_lambda_days,
                   read_a = 2, read_b = 4,
                   read_number_set = read_number_set)

            xs223_thstdonly_list.append(xs_calc_results[xs])
            xs223_thstdonly_err_list.append(xs_calc_results[xs_err])
            fraction_decayed_223_thstdonly_list.append(xs_calc_results[fraction_decayed])
            xs223_thstdonly_t0_list.append(xs_calc_results[xs_t0])
            ac227_thstdonly_list.append(xs_calc_results[parent_activity])
            ac227_thstdonly_err_list.append(xs_calc_results[parent_activity_err])
            row_specific_errors.append(xs_calc_results[row_errors])
            read_specific_errors.update(xs_calc_results[read_errors])

#####   Read1 - Read4      #################################################################################################################
        elif 1 in read_number_set and 4 in read_number_set:
#            xs223_t0_list.append('(223xs) using read 1 instead of 2')
            row_specific_errors.append('(223xs) using read 1 instead of 2')
            xs_calc_results = xs_calculator (lvl2_main_df, sample_variable, sub_sample_variable, row_sample_variable, row_sub_sample_variable,
                   isotope = 'Ra-223', isotope_column_string = 'vdpm223_thstdonly (dpm/m^3)', isotope_column_string_err = 'vdpm223_thstdonly_err (dpm/m^3)',
                   isotope_lambda_days = ra223_lambda_days,
                   read_a = 1, read_b = 4,
                   read_number_set = read_number_set)

            xs223_thstdonly_list.append(xs_calc_results[xs])
            xs223_thstdonly_err_list.append(xs_calc_results[xs_err])
            fraction_decayed_223_thstdonly_list.append(xs_calc_results[fraction_decayed])
            xs223_thstdonly_t0_list.append(xs_calc_results[xs_t0])
            ac227_thstdonly_list.append(xs_calc_results[parent_activity])
            ac227_thstdonly_err_list.append(xs_calc_results[parent_activity_err])
            row_specific_errors.append(xs_calc_results[row_errors])
            read_specific_errors.update(xs_calc_results[read_errors])

#####   Reads missing       ##################################################################################################################
        else:
            row_specific_errors.append('(223xs_thstdonly,227Ac_thstdonly) Required reads not available')
            read_specific_errors.update({'(223xs_thstdonly, 227Ac_thstdonly)': 'Required reads not available'})
            ac227_thstdonly_list.append(error_flag)
            ac227_thstdonly_err_list.append(error_flag)
            xs223_thstdonly_list.append(error_flag)
            xs223_thstdonly_err_list.append(error_flag)
            xs223_thstdonly_t0_list.append(error_flag)
            fraction_decayed_223_thstdonly_list.append(error_flag)

###########################################################################################################################################
#       228Ra calculations
###########################################################################################################################################

        if 4 in read_number_set and 5 in read_number_set:

            read4_vdpm = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==4))
                    , 'vdpm224 (dpm/m^3)']
            read4_vdpm_err = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==4))
                    , 'vdpm224_err (dpm/m^3)']
            read4_datetime = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==4))
                    , 'Mid_Read_Datetime']

            read5_vdpm = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==5))
                    , 'vdpm224 (dpm/m^3)']
            read5_vdpm_err = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==5))
                    , 'vdpm224_err (dpm/m^3)']
            read5_datetime = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==5))
                    , 'Mid_Read_Datetime']

            read4_errors = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==4))
                    , 'Error_List']
            read5_errors = lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable)
                    & (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)
                    & (lvl2_main_df['read_number']==5))
                    , 'Error_List']

            if len(read4_vdpm)>1 and len(read5_vdpm)==1:
                # print ('>=')
                timedelta_read4_to_read5_list = []
                time_read4_to_read5_list = []
                for i in range(len(read4_datetime)):
                    timedelta_read4_to_read5_list.append(read5_datetime.iloc[0]-read4_datetime.iloc[i])
                    time_read4_to_read5_list.append(timedelta_read4_to_read5_list[-1].days + (timedelta_read4_to_read5_list[-1].seconds/(24*60*60)))
                time_read4_to_read5 = np.average(time_read4_to_read5_list)
                # print(time_read4_to_read5)
                ra228_list.append(
                                    read5_vdpm.iloc[0] - (np.average(read4_vdpm) - np.exp(-th228_lambda_days* (time_read4_to_read5)/(1.499* (np.exp(-ra228_lambda_days*time_read4_to_read5)-np.exp(-th228_lambda_days*time_read4_to_read5)))))
                                    )
                ra228_err_list.append(
                    np.sqrt(read5_vdpm_err.iloc[0]**2 + np.average(read4_vdpm_err)**2)
                )
                row_specific_errors.append('Multiple 4th reads Averaged')
                # read_specific_errors.update(xs_calc_results[read_errors])

            if len(read4_vdpm)==1 and len(read5_vdpm)>1:
                # print ('=>')

                timedelta_read4_to_read5_list = []
                time_read4_to_read5_list = []
                for i in range(len(read5_datetime)):
                    timedelta_read4_to_read5_list.append(read5_datetime.iloc[i]-read4_datetime.iloc[0])
                    time_read4_to_read5_list.append(timedelta_read4_to_read5_list[-1].days + (timedelta_read4_to_read5_list[-1].seconds/(24*60*60)))
                time_read4_to_read5 = np.average(time_read4_to_read5_list)
                # print(time_read4_to_read5)
                ra228_list.append(
                                    np.average(read5_vdpm) - (read4_vdpm.iloc[0] - np.exp(-th228_lambda_days* (time_read4_to_read5)/(1.499* np.exp((-ra228_lambda_days*time_read4_to_read5)-(-th228_lambda_days*time_read4_to_read5)))))
                                    )
                ra228_err_list.append(
                    np.sqrt(np.average(read5_vdpm_err)**2 + read4_vdpm_err.iloc[0]**2)
                )
                row_specific_errors.append('Multiple 5th reads Averaged')
                # read_specific_errors.update(xs_calc_results[read_errors])

            if len(read4_vdpm)>1 and len(read5_vdpm)>1:
                # print ('>>')
                timedelta_read4_to_read5_list = []
                time_read4_to_read5_list = []
                for i in range(len(read5_datetime)):
                    for j in range(len(read4_datetime)):
                        timedelta_read4_to_read5_list.append(read5_datetime.iloc[i]-read4_datetime.iloc[j])
                        time_read4_to_read5_list.append(timedelta_read4_to_read5_list[-1].days + (timedelta_read4_to_read5_list[-1].seconds/(24*60*60)))

                time_read4_to_read5 = np.average(time_read4_to_read5_list)

                # print(time_read4_to_read5)
                ra228_list.append(
                                    np.average(read5_vdpm) - (np.average(read4_vdpm) - np.exp(-th228_lambda_days* (time_read4_to_read5)/(1.499* np.exp((-ra228_lambda_days*time_read4_to_read5)-(-th228_lambda_days*time_read4_to_read5)))))
                                    )

                ra228_err_list.append(
                    np.sqrt(np.average(read5_vdpm_err)**2 + np.average(read4_vdpm_err)**2)
                )
                row_specific_errors.append('Multiple 4th and 5th reads averaged')
                # read_specific_errors.update(xs_calc_results[read_errors])

            if len(read4_vdpm)==1 and len(read5_vdpm)==1:
                # print ('==')
                timedelta_read4_to_read5 = read5_datetime.iloc[0]-read4_datetime.iloc[0]
                time_read4_to_read5 = timedelta_read4_to_read5.days + (timedelta_read4_to_read5.seconds/(24*60*60))
                # print(time_read4_to_read5)
                ra228_list.append(
                                    read5_vdpm.iloc[0] - (read4_vdpm.iloc[0] - np.exp(-th228_lambda_days* (time_read4_to_read5)/(1.499* np.exp((-ra228_lambda_days*time_read4_to_read5)-(-th228_lambda_days*time_read4_to_read5)))))
                                    )
                ra228_err_list.append(
                    np.sqrt(read5_vdpm_err.iloc[0]**2 + read4_vdpm_err.iloc[0]**2)
                )
            ###############################################################################################################################
            # Bring Diego-Feliu Errors through from results dataframe to summary dataframe.
            ###############################################################################################################################
            read_specific_errors.update({'Read4' : list(read4_errors.iloc[0].keys()),'Read5' :list(read5_errors.iloc[0].keys())})

        else:
            ra228_list.append('(228Ra) Required reads not available')
            ra228_err_list.append('(228Ra) Required reads not available')
            # read_specific_errors.update({'(228Ra)': 'Required reads not available'})

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

###########################################################################################################################################
#       226Ra calculations
###########################################################################################################################################
        ra226_row_list = list(lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable) &
                                 (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)), 'vdpm226 (dpm/m^3)'])
        ra226_err_row_list = list(lvl2_main_df.loc[((lvl2_main_df[sample_variable]==row_sample_variable) &
                                 (lvl2_main_df[sub_sample_variable]==row_sub_sample_variable)), 'vdpm226_err (dpm/m^3)'])

        upper_limit_226_err = np.average(ra226_err_row_list)+np.std(ra226_err_row_list)
        lower_limit_226_err = np.average(ra226_err_row_list)-np.std(ra226_err_row_list)

        if len(ra226_row_list)==0 :
            ra226_list.append(error_flag)
            row_specific_errors.append('No_226_Value calculated')

        else:
            ra226_row_list = [i for i in ra226_row_list if i != -999]
            upper_limit_226 = np.average(ra226_row_list)+np.std(ra226_row_list)
            lower_limit_226 = np.average(ra226_row_list)-np.std(ra226_row_list)
            ra226_outliers_removed_list = [i for i in ra226_row_list if i < upper_limit_226 and i > lower_limit_226]
            # ra226_outliers_removed_list = [i for i in ra226_row_list if i != -999]

            if len(ra226_outliers_removed_list)>1:
                ra226_list.append(np.average(ra226_outliers_removed_list))
            if len(ra226_outliers_removed_list)==1:
                ra226_list.append(np.average(ra226_outliers_removed_list))
                row_specific_errors.append('Only one 226-Ra value averaged')

            if len(ra226_outliers_removed_list)<1:
                ra226_list.append(np.average(error_flag))
                row_specific_errors.append('226-Ra values non-concordant')

        if len(ra226_err_row_list)==0 :
            ra226_err_list.append(error_flag)
            row_specific_errors.append('No_226_err_Value calculated')

        else:
            ra226_err_row_list = [i for i in ra226_err_row_list if i != -999]
            upper_limit_226_err = np.average(ra226_err_row_list)+np.std(ra226_err_row_list)
            lower_limit_226_err = np.average(ra226_err_row_list)-np.std(ra226_err_row_list)
            ra226_err_outliers_removed_list = [i for i in ra226_err_row_list if i < upper_limit_226_err and i > lower_limit_226_err]
            # ra226_err_outliers_removed_list = [i for i in ra226_err_row_list if i != -999]
            #print (ra226_err_outliers_removed_list)

            if len(ra226_err_outliers_removed_list)>1:
                ra226_err_list.append(np.average(ra226_err_outliers_removed_list))
            if len(ra226_err_outliers_removed_list)==1:
                ra226_err_list.append(np.average(ra226_err_outliers_removed_list))
                row_specific_errors.append('Only one 226-Ra_err value averaged')

            if len(ra226_err_outliers_removed_list)<1:
                ra226_err_list.append(error_flag)
                row_specific_errors.append('226-Ra_err values non-concordant')

###########################################################################################################################################
        error_list.append(row_specific_errors)
        read_error_list.append(read_specific_errors)

        # print (len(error_list), len(read_error_list), len(xs224_list))

    summary_df['224xs'] = xs224_list
    summary_df['224xs_err'] = xs224_err_list
    summary_df['224xs_t0'] = xs224_t0_list
    summary_df['224xs_t0_err'] = summary_df['224xs_err']*(summary_df['224xs_t0']/summary_df['224xs'])
    summary_df['Fraction_of_original_224_remaining'] = fraction_decayed_224_list
    summary_df['228Th'] = th228_list
    summary_df['228Th_err'] = th228_err_list

    summary_df['223xs'] = xs223_list
    summary_df['223xs_err'] = xs223_err_list
    summary_df['223xs_t0'] = xs223_t0_list
    summary_df['223xs_t0_err'] = summary_df['223xs_err']*(summary_df['223xs_t0']/summary_df['223xs'])
    summary_df['Fraction_of_original_223_remaining'] = fraction_decayed_223_list
    summary_df['227Ac'] = ac227_list
    summary_df['227Ac_err'] = ac227_err_list

    summary_df['223xs_thstdonly'] = xs223_thstdonly_list
    summary_df['223xs_thstdonly_err'] = xs223_thstdonly_err_list
    summary_df['223xs_thstdonly_t0'] = xs223_thstdonly_t0_list
    summary_df['223xs_thstdonly_t0_err'] = summary_df['223xs_thstdonly_err']*(summary_df['223xs_thstdonly_t0']/summary_df['223xs_thstdonly'])
    summary_df['Fraction_of_original_223_thstdonly_remaining'] = fraction_decayed_223_thstdonly_list
    summary_df['227Ac_thstdonly'] = ac227_thstdonly_list
    summary_df['227Ac_thstdonly_err'] = ac227_thstdonly_err_list

    summary_df['228Ra'] = ra228_list
    summary_df['228Ra_err'] = ra228_err_list
    summary_df['xs_calc_errors'] = error_list
    summary_df['read_errors'] = read_error_list

    summary_df['226Ra'] = ra226_list
    summary_df['226Ra_err'] = ra226_err_list

    cols = list(log_df.columns) + ['224xs','224xs_err', '224xs_t0', '224xs_t0_err',	'Fraction_of_original_224_remaining', '228Th',	 '228Th_err',
                                   '223xs','223xs_err', '223xs_t0', '223xs_t0_err', 'Fraction_of_original_223_remaining', '227Ac',	'227Ac_err',
                                   '223xs_thstdonly','223xs_thstdonly_err', '223xs_thstdonly_t0', '223xs_thstdonly_t0_err', 'Fraction_of_original_223_thstdonly_remaining', '227Ac_thstdonly',	'227Ac_thstdonly_err',
                                   '226Ra', '226Ra_err', '228Ra', '228Ra_err', 'xs_calc_errors', 'read_errors' ]
    summary_df = summary_df[cols]

    summary_df.to_csv(output_directory/Path('Dataframes/Summary_Sample_Results_Dataframe'+time.strftime("%Y-%m-%d_%H%M%S")+'.csv'))

    return(summary_df)

# print (create_summary_dataframe(lvl2_main_df, log_df, sample_variable, sub_sample_variable,output_directory))


text:
""" Base of the object-oriented model for modelpacks """
#***************************************************************************************************
# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
# in this software.
# 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 or in the LICENSE file in the root pyGSTi directory.
#***************************************************************************************************

from abc import ABC as _ABC, abstractmethod as _abstractmethod
from pathlib import Path as _Path
import numpy as _np
import gzip as _gzip
import pickle as _pickle

from ..objects.polynomial import bulk_load_compact_polys as _bulk_load_compact_polys
from ..objects.circuit import Circuit as _Circuit
from ..construction.circuitconstruction import circuit_list as _circuit_list
from ..construction.modelconstruction import build_explicit_model as _build_explicit_model
from ..construction.stdlists import make_lsgst_structs as _make_lsgst_structs
from ..protocols import gst as _gst

class ModelPack(_ABC):
    """ ABC of all derived modelpack types"""
    description = None
    gates = None
    _sslbls = None

    def __init__(self):
        self._gscache = {}

    @_abstractmethod
    def _target_model(self, sslbls):
        pass

    def _build_explicit_target_model(self, sslbls, gate_names, gate_expressions, **kwargs):
        """
        A helper function for derived classes which create explicit models.
        Updates gate names and expressions with a given set of state-space labels.
        """
        def update_gatename(gn):
            return gn[0:1] + tuple([sslbls[i] for i in gn[1:]])
        full_sslbls = [sslbls]  # put all sslbls in single tensor product block
        updated_gatenames = [update_gatename(gn) for gn in gate_names]
        updated_gateexps = [gexp.format(*sslbls) for gexp in gate_expressions]
        return _build_explicit_model(full_sslbls, updated_gatenames, updated_gateexps, **kwargs)

    def target_model(self, parameterization_type="full", sim_type="auto", qubit_labels=None):
        """
        Returns a copy of the target model in the given parameterization.

        Parameters
        ----------
        parameterization_type : {"TP", "CPTP", "H+S", "S", ... }
            The gate and SPAM vector parameterization type. See
            :function:`Model.set_all_parameterizations` for all allowed values.

        sim_type : {"auto", "matrix", "map", "termorder" }
            The simulator type to be used for model calculations (leave as
            "auto" if you're not sure what this is).

        qubit_labels : tuple, optional
            A tuple of qubit labels, e.g. ('Q0', 'Q1') or (0, 1).  The default
            are the integers starting at 0.

        Returns
        -------
        Model
        """
        qubit_labels = self._sslbls if (qubit_labels is None) else tuple(qubit_labels)
        assert(len(qubit_labels) == len(self._sslbls)), \
            "Expected %d qubit labels and got: %s!" % (len(self._sslbls), str(qubit_labels))

        if (parameterization_type, sim_type, qubit_labels) not in self._gscache:
            # cache miss
            mdl = self._target_model(qubit_labels)
            mdl.set_all_parameterizations(parameterization_type)  # automatically sets sim_type
            if parameterization_type == "H+S terms":
                assert (sim_type == "auto" or sim_type in ("termorder", "termgap", "termdirect")), \
                    "Invalid `sim_type` argument for H+S terms: %s!" % sim_type
                if sim_type == "auto":
                    simt = "termorder"
                    simt_kwargs = {'max_order': 1}
                else:
                    simt = sim_type  # don't update sim_type b/c gscache
                    simt_kwargs = {'max_order': 1}  # TODO: update so user can specify other args?

                key_path, val_path = self._get_cachefile_names(parameterization_type, simt)
                if key_path.exists() and val_path.exists():
                    simt_kwargs['cache'] = _load_calccache(key_path, val_path)

                mdl.set_simtype(simt, **simt_kwargs)
            else:
                if sim_type != "auto":
                    mdl.set_simtype(sim_type)

            # finally cache result
            self._gscache[(parameterization_type, sim_type, qubit_labels)] = mdl

        return self._gscache[(parameterization_type, sim_type, qubit_labels)].copy()

    def _get_cachefile_names(self, param_type, sim_type):
        """ Get the standard cache file names for a modelpack """

        if param_type == "H+S terms":
            cachePath = _Path(__file__).absolute().parent / "caches"

            assert (sim_type == "auto" or sim_type.startswith("termorder:")), "Invalid `sim_type` argument!"
            termOrder = 1 if sim_type == "auto" else int(sim_type.split(":")[1])
            fn = ("cacheHS%d." % termOrder) + self.__module__
            return cachePath / (fn + "_keys.pkz"), cachePath / (fn + "_vals.npz")
        else:
            raise ValueError("No cache files used for param-type=%s" % param_type)

class GSTModelPack(ModelPack):
    """ ABC for modelpacks with GST information"""
    _germs = None
    _germs_lite = None
    _fiducials = None
    _prepfiducials = None
    _measfiducials = None

    global_fidPairs = None
    global_fidPairs_lite = None
    _pergerm_fidPairsDict = None
    _pergerm_fidPairsDict_lite = None

    def __init__(self):
        self._gscache = {}

    def _indexed_circuits(self, prototype, index):
        if index is None: index = self._sslbls
        assert(len(index) == len(self._sslbls)), "Wrong number of labels in: %s" % str(index)
        if prototype is not None:
            return _circuit_list(_transform_circuittup_list(prototype, index), index)

    def _indexed_circuitdict(self, prototype, index):
        if index is None: index = self._sslbls
        assert(len(index) == len(self._sslbls)), "Wrong number of labels in: %s" % str(index)
        if prototype is not None:
            return {_Circuit(_transform_circuit_tup(k, index), line_labels=index): val for k, val in prototype.items()}

    def germs(self, qubit_labels=None, lite=True):
        if lite and self._germs_lite is not None:
            return self._indexed_circuits(self._germs_lite, qubit_labels)
        else:
            return self._indexed_circuits(self._germs, qubit_labels)

    def fiducials(self, qubit_labels=None):
        return self._indexed_circuits(self._fiducials, qubit_labels)

    def prep_fiducials(self, qubit_labels=None):
        return self._indexed_circuits(self._prepfiducials, qubit_labels)

    def meas_fiducials(self, qubit_labels=None):
        return self._indexed_circuits(self._measfiducials, qubit_labels)

    def pergerm_fidpair_dict(self, qubit_labels=None):
        return self._indexed_circuitdict(self._pergerm_fidPairsDict, qubit_labels)

    def pergerm_fidpair_dict_lite(self, qubit_labels=None):
        return self._indexed_circuitdict(self._pergerm_fidPairsDict_lite, qubit_labels)

    def get_gst_experiment_design(self, max_max_length, qubit_labels=None, fpr=False, lite=True, **kwargs):
        """ Construct a :class:`protocols.gst.StandardGSTDesign` from this modelpack

        Parameters
        ----------
        max_max_length : number
            The greatest maximum-length to use. Equivalent to
            constructing a :class:`StandardGSTDesign` with a
            `maxLengths` list of powers of two less than or equal to
            the given value.

        qubit_labels : tuple, optional
            A tuple of qubit labels.  None means the integers starting at 0.

        fpr : bool, optional
            Whether to reduce the number of sequences using fiducial
            pair reduction (FPR).

        lite : bool, optional
            Whether to use a smaller "lite" list of germs. Unless you know
            you have a need to use the more pessimistic "full" set of germs,
            leave this set to True.

        **kwargs :
            Additional arguments to pass to :class:`StandardGSTDesign`

        Returns
        -------
        :class:`StandardGSTDesign`
        """
        for k in kwargs.keys():
            if k not in ('germLengthLimits', 'keepFraction', 'keepSeed', 'includeLGST', 'nest', 'sequenceRules',
                         'opLabelAliases', 'dscheck', 'actionIfMissing', 'verbosity', 'add_default_protocol'):
                raise ValueError("Invalid argument '%s' to StandardGSTDesign constructor" % k)

        if qubit_labels is None: qubit_labels = self._sslbls
        assert(len(qubit_labels) == len(self._sslbls)), \
            "Expected %d qubit labels and got: %s!" % (len(self._sslbls), str(qubit_labels))

        if fpr:
            fidpairs = self.pergerm_fidpair_dict_lite(qubit_labels) if lite else \
                self.pergerm_fidpair_dict(qubit_labels)
            if fidpairs is None:
                raise ValueError("No FPR information for lite=%s" % lite)
        else:
            fidpairs = None

        return _gst.StandardGSTDesign(
            self._target_model(qubit_labels),
            self.prep_fiducials(qubit_labels),
            self.meas_fiducials(qubit_labels),
            self.germs(qubit_labels, lite),
            list(_gen_max_length(max_max_length)),
            kwargs.get('germLengthLimits', None),
            fidpairs,
            kwargs.get('keepFraction', 1),
            kwargs.get('keepSeed', None),
            kwargs.get('includeLGST', True),
            kwargs.get('nest', True),
            kwargs.get('sequenceRules', None),
            kwargs.get('opLabelAliases', None),
            kwargs.get('dscheck', None),
            kwargs.get('actionIfMissing', None),
            qubit_labels,
            kwargs.get('verbosity', 0),
            kwargs.get('add_default_protocol', False),
        )

    def get_gst_circuits_struct(self, max_max_length, qubit_labels=None, fpr=False, lite=True, **kwargs):
        """ Construct a :class:`pygsti.objects.LsGermsStructure` from this modelpack.

        Parameters
        ----------
        max_max_length : number
            The greatest maximum-length to use. Equivalent to
            constructing a cicuit struct with a `maxLengths`
            list of powers of two less than or equal to
            the given value.

        qubit_labels : tuple, optional
            A tuple of qubit labels.  None means the integers starting at 0.

        fpr : bool, optional
            Whether to reduce the number of sequences using fiducial
            pair reduction (FPR).

        lite : bool, optional
            Whether to use a smaller "lite" list of germs. Unless you know
            you have a need to use the more pessimistic "full" set of germs,
            leave this set to True.

        **kwargs :
            Additional arguments to pass to :function:`make_lsgst_structs`

        Returns
        -------
        :class:`pygsti.objects.LsGermsStructure`
        """
        if fpr:
            fidpairs = self.pergerm_fidpair_dict_lite(qubit_labels) if lite else \
                self.pergerm_fidpair_dict(qubit_labels)
            if fidpairs is None:
                raise ValueError("No FPR information for lite=%s" % lite)
        else:
            fidpairs = None

        qubit_labels = self._sslbls if (qubit_labels is None) else tuple(qubit_labels)
        assert(len(qubit_labels) == len(self._sslbls)), \
            "Expected %d qubit labels and got: %s!" % (len(self._sslbls), str(qubit_labels))

        structs = _make_lsgst_structs(self._target_model(qubit_labels),  # Note: only need gate names here
                                      self.prep_fiducials(qubit_labels),
                                      self.meas_fiducials(qubit_labels),
                                      self.germs(qubit_labels, lite),
                                      list(_gen_max_length(max_max_length)),
                                      fidpairs,
                                      **kwargs)
        return structs[-1]  # just return final struct (for longest sequences)

class RBModelPack(ModelPack):
    _clifford_compilation = None

    def clifford_compilation(self, qubit_labels=None):
        if qubit_labels is None: qubit_labels = self._sslbls
        assert(len(qubit_labels) == len(self._sslbls)), "Wrong number of labels in: %s" % str(qubit_labels)
        return {clifford_name: _Circuit(_transform_circuit_tup(circuittup_of_native_gates,
                                                               qubit_labels), line_labels=qubit_labels)
                for clifford_name, circuittup_of_native_gates in self._clifford_compilation.items()}

# Helper functions
def _load_calccache(key_path, val_path):
    """
    The complement to _write_calccache, this function loads a cache
    dictionary from key and value filenames.

    Parameters
    ----------
    key_path, val_path : str
        key and value filenames.

    Returns
    -------
    dict
        The cache of calculated (compact) polynomials.
    """
    #print("Loading cache..."); t0 = _time.time()
    with _gzip.open(key_path, "rb") as f:
        keys = _pickle.load(f)
    npfile = _np.load(val_path)
    vals = _bulk_load_compact_polys(npfile['vtape'], npfile['ctape'], keep_compact=True)
    calc_cache = {k: v for k, v in zip(keys, vals)}
    #print("Done in %.1fs" % (_time.time()-t0))
    return calc_cache

def _transform_layer_tup(layer_tup, index_tup):
    if len(layer_tup) > 1:
        # This assumes simple layer labels that == a gate label
        # We could use the Label object here in future to support more complex labels
        lbl, *idx = layer_tup
        return (lbl, *[index_tup[i] for i in idx])
    else:
        return layer_tup

def _transform_circuit_tup(circuit_tup, index_tup):
    return tuple(_transform_layer_tup(layer_lbl, index_tup) for layer_lbl in circuit_tup)

def _transform_circuittup_list(circuittup_list, index_tup):
    """ Transform the indices of the tuples in a circuit list with the given index factory """
    for circuit_tup in circuittup_list:
        yield _transform_circuit_tup(circuit_tup, index_tup)

def _gen_max_length(max_max_length):
    i = 1
    while i <= max_max_length:
        yield i
        i *= 2

max_stars_repo_name: drewrisinger/pyGSTi

max_issues_repo_name: drewrisinger/pyGSTi

max_forks_repo_name: drewrisinger/pyGSTi


text:
"""
This recipe trains and evaluates a deep clustering model
on the clean data from the WHAM dataset with 8k. It's divided into
three big chunks: data preparation, training, and evaluation.
Final output of this script:

┌───────────────────┬────────────────────┬────────────────────┐
│                   │ OVERALL (N = 6000) │                    │
╞═══════════════════╪════════════════════╪════════════════════╡
│        SAR        │        SDR         │        SIR         │
├───────────────────┼────────────────────┼────────────────────┤
│ 11.07829052874508 │ 10.737156798640111 │ 23.704177123014816 │
└───────────────────┴────────────────────┴────────────────────┘

Last run on 3/20/20.
"""
import os
import multiprocessing
import logging
import shutil
import json
import glob
from concurrent.futures import ThreadPoolExecutor

import torch
from torch import optim
import tqdm
import numpy as np
import termtables

import nussl
from nussl import ml, datasets, utils, separation, evaluation

# ----------------------------------------------------
# ------------------- SETTING UP ---------------------
# ----------------------------------------------------

# seed this recipe for reproducibility
utils.seed(0)

# set up logging
logging.basicConfig(
    format='%(asctime)s,%(msecs)d %(levelname)-8s [%(filename)s:%(lineno)d] %(message)s',
    datefmt='%Y-%m-%d:%H:%M:%S', level=logging.INFO)

# make sure this is set to WHAM root directory
WHAM_ROOT = os.getenv("WHAM_ROOT")
CACHE_ROOT = os.getenv("CACHE_ROOT")
NUM_WORKERS = multiprocessing.cpu_count() // 2
OUTPUT_DIR = os.path.expanduser('~/.nussl/recipes/wham_dpcl/run3_1e-2')
RESULTS_DIR = os.path.join(OUTPUT_DIR, 'results')
MODEL_PATH = os.path.join(OUTPUT_DIR, 'checkpoints', 'best.model.pth')
DEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'
BATCH_SIZE = 25
MAX_EPOCHS = 100
CACHE_POPULATED = True
LEARNING_RATE = 1e-3
PATIENCE = 5
GRAD_NORM = 1e-2

shutil.rmtree(os.path.join(RESULTS_DIR), ignore_errors=True)
os.makedirs(RESULTS_DIR, exist_ok=True)
shutil.rmtree(os.path.join(OUTPUT_DIR, 'tensorboard'), ignore_errors=True)

def construct_transforms(cache_location):
    # stft will be 32ms wlen, 8ms hop, sqrt-hann, at 8khz sample rate by default
    tfm = datasets.transforms.Compose([
        datasets.transforms.MagnitudeSpectrumApproximation(),  # take stfts and get ibm
        datasets.transforms.MagnitudeWeights(),  # get magnitude weights
        datasets.transforms.ToSeparationModel(),  # convert to tensors
        datasets.transforms.Cache(cache_location),  # up to here gets cached
        datasets.transforms.GetExcerpt(400)  # get 400 frame excerpts (3.2 seconds)
    ])
    return tfm

def cache_dataset(_dataset):
    cache_dataloader = torch.utils.data.DataLoader(
        _dataset, num_workers=NUM_WORKERS, batch_size=BATCH_SIZE)
    ml.train.cache_dataset(cache_dataloader)
    _dataset.cache_populated = True

tfm = construct_transforms(os.path.join(CACHE_ROOT, 'tr'))
dataset = datasets.WHAM(WHAM_ROOT, split='tr', transform=tfm,
                        cache_populated=CACHE_POPULATED)

tfm = construct_transforms(os.path.join(CACHE_ROOT, 'cv'))
val_dataset = datasets.WHAM(WHAM_ROOT, split='cv', transform=tfm,
                            cache_populated=CACHE_POPULATED)

if not CACHE_POPULATED:
    # cache datasets for speed
    cache_dataset(dataset)
    cache_dataset(val_dataset)

# ----------------------------------------------------
# -------------------- TRAINING ----------------------
# ----------------------------------------------------

# reload after caching
train_sampler = torch.utils.data.sampler.RandomSampler(dataset)
val_sampler = torch.utils.data.sampler.RandomSampler(val_dataset)

dataloader = torch.utils.data.DataLoader(dataset, num_workers=NUM_WORKERS,
                                         batch_size=BATCH_SIZE, sampler=train_sampler)
val_dataloader = torch.utils.data.DataLoader(val_dataset, num_workers=NUM_WORKERS,
                                             batch_size=BATCH_SIZE, sampler=val_sampler)

n_features = dataset[0]['mix_magnitude'].shape[1]

# builds a baseline model with 4 recurrent layers, 600 hidden units, bidirectional
# and 20 dimensional embedding
config = ml.networks.builders.build_recurrent_dpcl(
    n_features, 600, 4, True, 0.3, 20, ['sigmoid', 'unit_norm'],
    normalization_class='BatchNorm')
model = ml.SeparationModel(config).to(DEVICE)
logging.info(model)

optimizer = optim.Adam(model.parameters(), lr=LEARNING_RATE)
scheduler = optim.lr_scheduler.ReduceLROnPlateau(
    optimizer, factor=0.5, patience=PATIENCE)

# set up the loss function
loss_dictionary = {
    'DeepClusteringLoss': {'weight': 1.0}
}

# set up closures for the forward and backward pass on one batch
train_closure = ml.train.closures.TrainClosure(loss_dictionary, optimizer, model)
val_closure = ml.train.closures.ValidationClosure(loss_dictionary, model)

# set up engines for training and validation
trainer, validator = ml.train.create_train_and_validation_engines(train_closure, val_closure,
                                                                  device=DEVICE)

# attach handlers for visualizing output and saving the model
ml.train.add_stdout_handler(trainer, validator)
ml.train.add_validate_and_checkpoint(OUTPUT_DIR, model, optimizer, dataset,
                                     trainer, val_data=val_dataloader, validator=validator)
ml.train.add_tensorboard_handler(OUTPUT_DIR, trainer)

# add a handler to set up patience
@trainer.on(ml.train.ValidationEvents.VALIDATION_COMPLETED)
def step_scheduler(trainer):
    val_loss = trainer.state.epoch_history['validation/loss'][-1]
    scheduler.step(val_loss)

# add a handler to set up gradient clipping
@trainer.on(ml.train.BackwardsEvents.BACKWARDS_COMPLETED)
def clip_gradient(trainer):
    torch.nn.utils.clip_grad_norm_(model.parameters(), GRAD_NORM)

# train the model
trainer.run(dataloader, max_epochs=MAX_EPOCHS)

# ----------------------------------------------------
# ------------------- EVALUATION ---------------------
# ----------------------------------------------------

test_dataset = datasets.WHAM(WHAM_ROOT, sample_rate=8000, split='tt')
# make a deep clustering separator with an empty audio signal initially
# this one will live on gpu and be used in a threadpool for speed
dpcl = separation.deep.DeepClustering(
    nussl.AudioSignal(), num_sources=2, model_path=MODEL_PATH, device='cuda')

def forward_on_gpu(audio_signal):
    # set the audio signal of the object to this item's mix
    dpcl.audio_signal = audio_signal
    features = dpcl.extract_features()
    return features

def separate_and_evaluate(item, features):
    separator = separation.deep.DeepClustering(item['mix'], num_sources=2)
    estimates = separator(features)

    evaluator = evaluation.BSSEvalScale(
        list(item['sources'].values()), estimates, compute_permutation=True)
    scores = evaluator.evaluate()
    output_path = os.path.join(RESULTS_DIR, f"{item['mix'].file_name}.json")
    with open(output_path, 'w') as f:
        json.dump(scores, f)

pool = ThreadPoolExecutor(max_workers=NUM_WORKERS)
for i, item in enumerate(tqdm.tqdm(test_dataset)):
    features = forward_on_gpu(item['mix'])
    if i == 0:
        separate_and_evaluate(item, features)
    else:
        pool.submit(separate_and_evaluate, item, features)
pool.shutdown(wait=True)

json_files = glob.glob(f"{RESULTS_DIR}/*.json")
df = evaluation.aggregate_score_files(json_files)

overall = df.mean()
headers = ["", f"OVERALL (N = {df.shape[0]})", ""]
metrics = ["SAR", "SDR", "SIR"]
data = np.array(df.mean()).T

data = [metrics, data]
termtables.print(data, header=headers, padding=(0, 1), alignment="ccc")

max_issues_repo_name: titocaco/nussl

max_forks_repo_name: titocaco/nussl


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Nov  7 20:21:46 2016

@author: smrak
"""

import numpy as np
import datetime
import plotting
folder = '/home/smrak/Documents/TheMahali/imf/'

#def openIMF():
fname = 'omniweb20151007.csv'
Y, DinY, H, M, Bx, By, Bz, AE = np.loadtxt(folder+fname,dtype='S',
                                           unpack=True,delimiter=',')
Bx = Bx.astype(float)
By = By.astype(float)
Bz = Bz.astype(float)
AE = AE.astype(int)
H = H.astype(int)
M = M.astype(int)
t = []
for i in range(len(AE)):
    t.append(datetime.datetime(2015,10,7, H[i], M[i], 0))
t = np.array(t)

start = datetime.datetime(2015,10,7, 6, 15, 0)
stop = datetime.datetime(2015,10,7, 6, 40, 0)
idx = np.where( (t>=start) & (t<=stop) )[0]

obstimes = t[idx]

plotting.plotimf(t, Bx, By, Bz, AE/1000, ylim1=[-20,20], ylim2=[0,1.8],
                 ytick1=[-20,-15,-10,-5,0,5,10,15,20], ytick2=[0,0.4,0.8,1.2,1.6],
                 obstimes=obstimes, centerline=True, lw=2,
                 xlim=[datetime.datetime(2015,10,7,2,0,0),
                      datetime.datetime(2015,10,7,11,0,0)])

#mahaliPlot.plotImfField(np.array(t), Bx, By, Bz, AE)