text:
import numpy as np

class RandomMachine(object):
    """
    List of numpy random state
    """

    MAX_SEED = int(10**9)

    def __init__(self, basis_seed=0, call_length=0):
        self.rand_state = np.random.RandomState(basis_seed)
        self.length = int(call_length)
        self.rand_seed_seq = None
        self.update_seed_seq()
        self.list_random_state = self.gen_random_state_list()

    def renew(self):
        update_nums = self.rand_state.randint(self.MAX_SEED, size=(len(self), ))
        for i in range(len(self.list_random_state)):
            self.list_random_state[i].randint(update_nums[i])

    def gen_random_state_list(self):
        return [np.random.RandomState(self.rand_seed_seq[i])
                for i in range(len(self))]

    def update_seed_seq(self):
        self.rand_seed_seq = self.rand_state.randint(self.MAX_SEED, size=(len(self), ))

    def get_fix_machine(self, ind):
        new_seed = self.list_random_state[ind].randint(self.MAX_SEED)
        return np.random.RandomState(new_seed)

    def __len__(self):
        return self.length


text:
"""
This file takes the raw JSON output from the GT and Google API and
creates a nice paired-sentence JSON file which contains both the ground
truth and transcribed phrases, as well as some metadata.

This standardized format is used for all downstream metrics computation.

Example entry:
 '556': {
	 'hash': '0dd01803e57a3b95fcfab584bfb09aa604d80c0452ce5cd90f02669b0f9b9b5e',
	 'ts': 102.0,
	 'speaker': 'P',
	 'gt': 'hello there how are you doing',
	 'pred': 'hello their how arent you doing'
	}
"""
import re
import os
import sys
import json
import argparse
import numpy as np
from tqdm import tqdm
from typing import List, Dict
import preproc.util
import evaluation.config as config

def main(args):
    if not os.path.exists(args.machine_dir):
        print(f"Path does not exist: {args.machine_dir}")
        sys.exit(0)
    if not os.path.exists(args.gt_dir):
        print(f"Path does not exist: {args.gt_dir}")
        sys.exit(0)

    # Creates a single dictionary, indexed by GID, and contains
    # values: gt sentence, pred sentence, etc.
    paired: Dict[int, Dict] = {}

    # Loop over all sessions.
    ls = sorted(os.listdir(args.machine_dir))
    gid = 0  # Global sentence ID.
    for i in tqdm(range(len(ls)), desc="Processing Files"):
        filename = ls[i]
        hash_ = filename.replace(".json", "")
        # Load and standardize the ground truth.
        machine = load_json(os.path.join(args.machine_dir, filename))
        gt = load_json(os.path.join(args.gt_dir, filename))

        # Get segment timestamps.
        seg_ts = []
        seg_tag = []
        current_tag = None
        for j, tag in enumerate(gt["speaker_tags"]):
            # If first element, start the sequence OR
            # If we encounter a different tag, start a new sequence.
            if j == 0 or tag != current_tag:
                current_tag = tag
                seg_ts.append(gt["timestamps"][j])
                seg_tag.append(tag)

        # Create segments for GT and machine using timestamps.
        machine_segments = create_segments(seg_ts, machine)
        gt_segments = create_segments(seg_ts, gt)

        # Create the final dictionary.
        for j in range(len(gt_segments)):
            value = {
                "hash": hash_,
                "ts": seg_ts[j],
                "speaker": seg_tag[j],
                "gt": gt_segments[j],
                "pred": machine_segments[j],
            }
            paired[gid] = value
            gid += 1

    # Write output json file.
    with open(config.PAIRED_FQN, "w") as f:
        json.dump(paired, f, indent=2)

def is_repeated(arr: List[str], repeat: int) -> bool:
    """
	Checks whether a segment array of strings is doubled. That is,
	the first half contains the same elements as the second half.
	:param arr: List of strings.
	:param repeat: Check if `arr` is repeated `repeat` times.
	:return: True if array is doubled, False otherwise.
	"""
    if len(arr) % repeat != 0:
        return False

    pointers = np.arange(0, repeat) * int(len(arr) / repeat)

    while pointers[-1] < len(arr):
        unique = set([arr[idx] for idx in pointers])
        if len(unique) > 1:
            return False
        pointers += 1

    return True

def is_doubled(arr: List[str]) -> bool:
    """
	Checks whether a segment array of strings is doubled. That is,
	the first half contains the same elements as the second half.
	:param arr: List of strings.
	:return: True if array is doubled, False otherwise.
	"""
    if len(arr) % 2 != 0:
        return False

    first = 0
    second = int(len(arr) / 2)
    while second < len(arr):
        if arr[first] != arr[second]:
            return False
        first += 1
        second += 1
    return True

def create_segments(seg_ts: List[float], data: Dict) -> List[str]:
    """
	Takes a transcription dictionary and returns a list of segments.
	Each segment are the words occurring between two timestamps.

	:param seg_ts: List of timestamps denoting the start time of the segment.
	:param data: Dictionary with keys: timestamps, speakerTags, words.
	:return: List of segments, where each segment is a string.
	"""
    segments: List[str] = []

    data_ts = np.asarray(data["timestamps"])
    # Process everything except the last segment/bucket.
    for i in range(0, len(seg_ts) - 1):
        start = seg_ts[i]
        end = seg_ts[i + 1]
        # Find all words with this segment's timestamps.
        idx = np.where(np.logical_and(start <= data_ts, data_ts < end))[0]
        # Grab the actual words.
        if len(idx) > 0:
            buffer = []
            for j in idx:
                buffer.append(data["words"][j])
            if is_repeated(buffer, 2):
                clean = buffer[: int(len(buffer) / 2)]
            elif is_repeated(buffer, 3):
                clean = buffer[: int(len(buffer) / 3)]
            else:
                clean = buffer
            sentence = " ".join(clean)
            sentence = re.sub("\s+", " ", sentence).strip()
            segments.append(sentence)
        # If we have no words, append an empty segment.
        else:
            segments.append("")

    # Ignore the last segment because we don't know when it ends.
    # That is, the ground truth can end at 20 minutes, but the audio goes
    # until 30 minutes. As a result, we will have 10 minutes of machine
    # transcript that will be considered wrong.
    # Process the last segment/bucket.
    # idx = np.where(seg_ts[-1] <= data_ts)[0]
    # if len(idx) == 0:
    # 	segments.append('')
    # else:
    # 	buffer = []
    # 	for j in idx:
    # 		buffer.append(data['words'][j])
    # 	seg = ' '.join(buffer)
    # 	segments.append(seg)

    segments = [preproc.util.canonicalize_sentence(x) for x in segments]
    return segments

def get_phrases(buckets: List, data: Dict) -> Dict[float, List[str]]:
    """
	Returns a dictionary of (timestamp, words) where words contains the phrases spoken at `timestamp`.
	:param buckets: Python list of bucket timestamps (from GT).
	:param data: Dictionary of timestamps and words.
	:return: Dictionary of phrases.
	"""
    idx = 0  # Bucket index.
    phrases = {x: [] for x in buckets}
    for i in range(len(data["timestamps"])):
        ts, word = data["timestamps"][i], data["words"][i]
        if idx + 1 < len(buckets) and ts >= buckets[idx + 1]:
            idx += 1
        phrases[buckets[idx]].append(word)
    return phrases

def load_json(fqn: str):
    """
	Loads a json transcription file.

	A = json.load(fqn)
	A['results'] is a Python list of dictionaries.
	B = A['results'][0]
	B['alternatives'] is a Python list of dictionaries.
	C = B['alternatives'][0] is a dictionary of transcription results.
		C Keys: transcript, confidence, words
	C['words'] is a Python list of dictionaries.
	D = C['words'][0] contains the transcription for a single word.
		D Keys: startTime, endTime, word, confidence, speakerTag

	:param fqn: Path to the machine-generated transcription.
	:return:
	"""
    with open(fqn, "r") as f:
        A = json.load(f)

    timestamps = []
    words = []
    speaker_tags = []
    # For each word, add it to our list.
    for B in A["results"]:
        for C in B["alternatives"]:
            # Sometimes the 'words' key is not present.
            if "words" not in C:
                continue
            for D in C["words"]:
                # Get the core content.
                start_time = float(D["startTime"].replace("s", ""))
                timestamps.append(start_time)
                word = D["word"]
                word = preproc.util.canonicalize_word(word)
                words.append(word)

                # Add the speaker information.
                if "speakerTag" in D:
                    speaker_tags.append(D["speakerTag"])
                else:
                    speaker_tags.append("")

    result = {
        "timestamps": timestamps,
        "words": words,
        "speaker_tags": speaker_tags,
    }
    return result

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "machine_dir",
        type=str,
        help="Location of the predicted json transcriptions.",
    )
    parser.add_argument(
        "gt_dir",
        type=str,
        help="Location of the ground truth json transcription.",
    )
    main(parser.parse_args())

ext: py

max_stars_repo_name: som-shahlab/psych-audio

max_issues_repo_name: som-shahlab/psych-audio

max_forks_repo_name: som-shahlab/psych-audio

reason: import numpy


text:
#!/bin/python3
# Author: Dino Bollinger
# Licensed under BSD 3-Clause License, see included LICENSE file
"""
This script outputs a pickled numpy array indicating which indices are non-empty tiles.
"""

import os
import re
import sys
import pickle
import numpy as np

outfilename = "./non_empty.pkl"

MAXTILES:int = -1
if len(sys.argv) >= 2:
    MAXTILES = int(sys.argv[1])
else:
    print("Must specify maxtiles a first argument!", file=sys.stderr)
    exit(1)

indicator = np.zeros(MAXTILES)

tile_dir = "./tiles"
filename_prefix = "tile"
filename_suffix = ".PNG"

image_filenames = os.listdir(tile_dir)

for i in image_filenames:
    new_i = re.sub(filename_prefix, "", i)
    new_i = re.sub(filename_suffix, "", new_i)
    tile_num = int(new_i)
    indicator[tile_num] = 1

print(f"Total number of tiles in ART: {np.count_nonzero(indicator)}")
with open(outfilename, "wb") as fd:
    pickle.dump(indicator, fd, pickle.HIGHEST_PROTOCOL)
    print(f"Indicator array written to '{outfilename}'")


text:
import pickle
import re
import copy

import numpy as np
import pandas as pd
from torch import nn
import torch

from counterfit.core import config
from counterfit.core.state import TextTarget, Target

class MovieReviewsSentimentLSTM(nn.Module):
    """pre-trained LSTM model on 25 epochs for building sentiment analysis model on IMDB movies review dataset.
    """
    def __init__(self, no_layers, vocab_size, hidden_dim, embedding_dim, output_dim, drop_prob=0.5):
        # embedding_dim: number of expected features in the input `x`
        # hidden_dim: number of features in the hidden state `h`

        super(MovieReviewsSentimentLSTM, self).__init__()

        self.no_layers = no_layers # number of recurrent layers
        self.vocab_size = vocab_size
        self.hidden_dim = hidden_dim # The number of features in the hidden state h
        # embedding and LSTM layers
        self.embedding = nn.Embedding(vocab_size, embedding_dim)
        self.proj_size = 0
        self.output_dim = output_dim # The size of the output you desire from your RNN
        # dropout layer
        self.dropout = nn.Dropout(drop_prob) # a Dropout layer on the outputs of each LSTM layer except the last layer, with dropout probability equal to dropout

        #lstm
        self.lstm = nn.LSTM(input_size=embedding_dim, hidden_size=self.hidden_dim,
                           num_layers=no_layers, batch_first=True, proj_size=self.proj_size)

        # final fully connected linear and sigmoid layer
        self.fc = nn.Linear(self.hidden_dim, self.output_dim)
        self.sig = nn.Sigmoid()

    def forward(self, x, hidden):
        """Forward process of LSTM model

        Args:
            x ([tensor]): training data/batch_first

        Returns:
            Last sigmoid output and hidden state
        """
        batch_size = x.size(0)
        # embeddings and lstm_out
        embeds = self.embedding(x)  # shape: Batch x Sequence x Feature   since batch_first = True
        lstm_out, hidden = self.lstm(embeds, hidden)

        lstm_out = lstm_out.contiguous().view(-1, self.hidden_dim)

        # dropout and fully connected layer
        out = self.dropout(lstm_out)
        out = self.fc(out)

        # sigmoid function
        sig_out = self.sig(out)

        # reshape to be batch_size first
        sig_out = sig_out.view(batch_size, -1)

        sig_out = sig_out[:, -1] # get last batch of labels

        return sig_out, hidden

    def init_hidden(self, batch_size, device='cpu'):
        # initialize hidden weights (h,c) to 0
        weights = next(self.parameters()).data
        h = (weights.new(self.no_layers, batch_size, self.hidden_dim).zero_().to(device),
             weights.new(self.no_layers, batch_size, self.hidden_dim).zero_().to(device))

        return h

class MovieReviewsTarget(TextTarget):
    """Defining movie reviews target which is responsible for predicting the scores for a given input and convert scores to labels."""
    model_name = "moviereviews"
    model_data_type = "text"
    model_location = "local"
    model_endpoint = f"{config.targets_path}/{model_name}/movie_reviews_sentiment_analysis.pt"
    model_input_shape = (1,)
    model_output_classes = [0, 1]  # TextAttack current requires these to be integers that can be used as an index
    sample_input_path = f"{config.targets_path}/{model_name}/movie-reviews-scores-full.csv"
    vocab_file = f"{config.targets_path}/{model_name}/movie-reviews-vocab.pkl"

    X = []  # we'll populate this in the constructor

    def __init__(self):
        self.data = pd.read_csv(self.sample_input_path)
        print(f"[+] Total Movie Reviews: {len(self.data)}")
        self._load_x()
        self.vocab = self._load_vocab()
        self.model = self._load_model()

    def _load_x(self):
        # Append input reviews to X list
        for idx in range(len(self.data)):
            self.X.append(self.data['review'][idx])

    def _load_vocab(self):
        # Load vocabulary file; 1000 most occurence words
        with open(self.vocab_file, 'rb') as fp:
            vocab = pickle.load(fp)
        return vocab

    def preprocess_string(self, s):
        # Remove all non-word characters (everything except numbers and letters)
        s = re.sub(r"[^\w\s]", '', s)
        # Replace all runs of whitespaces with no space
        s = re.sub(r"\s+", '', s)
        # replace digits with no space
        s = re.sub(r"\d", '', s)

        return s

    def _load_model(self):
        # Load the LST model that's already trained
        no_layers = 2
        vocab_size = len(self.vocab) + 1 #extra 1 for padding purpose
        embedding_dim = 64
        output_dim=1
        hidden_dim = 256
        model = MovieReviewsSentimentLSTM(no_layers, vocab_size, hidden_dim, embedding_dim, output_dim, drop_prob=0.5)
        model.load_state_dict(copy.deepcopy(torch.load(self.model_endpoint, 'cpu')))
        model.eval()
        return model

    def padding_(self, sentences, seq_len):
        # Padding with zeros if sentence is less than required seq length
        features = np.zeros((len(sentences), seq_len), dtype=int)
        for ii, review in enumerate(sentences):
            if len(review) != 0:
                features[ii, -len(review):] = np.array(review)[:seq_len]
        return features

    def __call__(self, x):
        """This function takes list of input texts. For example., ["how are you?"]

        Args:
            x (list): [input_text]

        Returns:
            final_prob_scores: [[0.98, 0.02]] 0.98 probability score represents the sentence tone is positive and 0.02 score represents
        """
        final_prob_scores = []
        for text in x:
            word_seq = np.array([self.vocab[self.preprocess_string(word)] for word in text.split()
                            if self.preprocess_string(word) in self.vocab.keys()])
            word_seq = np.expand_dims(word_seq, axis=0)
            pad = torch.from_numpy(self.padding_(word_seq, 500))
            inputs = pad.to('cpu')
            batch_size = 1
            h = self.model.init_hidden(batch_size)
            h = tuple([each.data for each in h])
            output, h = self.model(inputs, h)
            probability = output.item()
            final_prob_scores.append([probability, 1.0-probability])
        return final_prob_scores # this must produce a list of class probabilities


text:
import numpy as np

a = np.array([1, 2, 3])
print(a)
print(a.dtype)
# [1 2 3]
# int64

a_float = a.astype(np.float32)
print(a_float)
print(a_float.dtype)
# [1. 2. 3.]
# float32

print(a)
print(a.dtype)
# [1 2 3]
# int64

a_float = a.astype(float)
print(a_float)
print(a_float.dtype)
# [1. 2. 3.]
# float64

a_str = a.astype('str')
print(a_str)
print(a_str.dtype)
# ['1' '2' '3']
# <U21

a_int = a.astype('int32')
print(a_int)
print(a_int.dtype)
# [1 2 3]
# int32

a = np.arange(50).reshape((5, 10)) / 10 - 2
print(a)
print(a.dtype)
# [[-2.  -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1]
#  [-1.  -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1]
#  [ 0.   0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9]
#  [ 1.   1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9]
#  [ 2.   2.1  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9]]
# float64

a_int = a.astype('int64')
print(a_int)
print(a_int.dtype)
# [[-2 -1 -1 -1 -1 -1 -1 -1 -1 -1]
#  [-1  0  0  0  0  0  0  0  0  0]
#  [ 0  0  0  0  0  0  0  0  0  0]
#  [ 1  1  1  1  1  1  1  1  1  1]
#  [ 2  2  2  2  2  2  2  2  2  2]]
# int64

print(np.round(a).astype(int))
# [[-2 -2 -2 -2 -2 -2 -1 -1 -1 -1]
#  [-1 -1 -1 -1 -1  0  0  0  0  0]
#  [ 0  0  0  0  0  0  1  1  1  1]
#  [ 1  1  1  1  1  2  2  2  2  2]
#  [ 2  2  2  2  2  2  3  3  3  3]]

my_round_int = lambda x: np.round((x * 2 + 1) // 2)

print(my_round_int(a).astype(int))
# [[-2 -2 -2 -2 -2 -1 -1 -1 -1 -1]
#  [-1 -1 -1 -1 -1  0  0  0  0  0]
#  [ 0  0  0  0  0  1  1  1  1  1]
#  [ 1  1  1  1  1  2  2  2  2  2]
#  [ 2  2  2  2  2  3  3  3  3  3]]

def my_round(x, digit=0):
    p = 10 ** digit
    s = np.copysign(1, x)
    return (s * x * p * 2 + 1) // 2 / p * s

print(my_round(a).astype(int))
# [[-2 -2 -2 -2 -2 -2 -1 -1 -1 -1]
#  [-1 -1 -1 -1 -1 -1  0  0  0  0]
#  [ 0  0  0  0  0  1  1  1  1  1]
#  [ 1  1  1  1  1  2  2  2  2  2]
#  [ 2  2  2  2  2  3  3  3  3  3]]


text:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import ConvexHull
import pandas as pd
import math
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import sympy as sym
from mpl_toolkits.mplot3d import Axes3D
import plotter

import pyenvelope

import homogenous_matrixes

xy_to_uv_file =  "C:/Users/Samarth/OneDrive - HKUST Connect/Thesis Project/IEDA/march/8/coarse_sheet.xls"
# xy_to_uv_file =  "C:/Users/Samarth/OneDrive - HKUST Connect/Thesis Project/IEDA/march/8/bunny.xls"
b = pd.read_excel(xy_to_uv_file)
uv_conv = pd.DataFrame(b)
uv_conv = uv_conv.drop(columns="u")
uv_conv = uv_conv.drop(columns="v")
uv_conv = uv_conv.drop(columns="i")
uv_conv = uv_conv.drop(columns="j")
pts = array_from_dataframe = uv_conv.to_numpy()

hull = ConvexHull(pts)

# plotter.threeD_plot(pts)

filler_points = []
for i in range(len(hull.points)):
    filler_points.append([0])

# print(filler_points)

points_matrix = np.append(hull.points, filler_points, axis=1)
# print(points_matrix)

face_eq = []

for i in range(len(hull.simplices)):
    point_1 = hull.simplices[i][0]
    point_2 = hull.simplices[i][1]
    point_3 = hull.simplices[i][2]
    # print(point_1, point_2, point_3)

    cord_point1 = hull.points[point_1]
    cord_point2 = hull.points[point_2]
    cord_point3 = hull.points[point_3]

    # print(cord_point1, cord_point2, cord_point3)

    v1 = cord_point3 - cord_point1
    v2 = cord_point2 - cord_point1
    cp = np.cross(v1, v2)
    coff_a , coff_b , coff_c = cp
    coff_d = (-1) * (np.dot(cp, cord_point3))

    # face equations in the form of ax +by + cz + d = 0 (abcd and coffecients )

    current_face = [coff_a, coff_b, coff_c, coff_d]
    # print(current_face)
    face_eq.append(current_face)

# print(face_eq)

#  get the distance of all the points from a given face and get the maximum of all

face_point_distance_array = []
max_face_point = []
face = j = 0
for j in range(len(face_eq)):

    current_distance_array = []
    for i in range(len(hull.points)):

        current_point = hull.points[i]
        p_1 = hull.points[i][0]
        p_2 = hull.points[i][1]
        p_3 = hull.points[i][2]

        # print(p_1, p_2, p_3)
        current_distance_num = abs((p_1 * face_eq[j][0]) + (p_2 * face_eq[j][1]) + (p_3 * face_eq[j][2]) + (face_eq[j][3]))
        current_distance_den =math.sqrt((face_eq[j][0] ** 2) + ((face_eq[j][1] ** 2 )+ (face_eq[j][2] ** 2) ))
        current_distance =  current_distance_num / current_distance_den
        # print(current_distance)
        current_distance_array.append(current_distance)
        face_point_distance_array_current = [j, i, current_distance]
        # print(face_point_distance_array_current)
        face_point_distance_array.append(face_point_distance_array_current)
    # print(current_distance_array)
    # print(max(current_distance_array))
    # print('\n')
    max_face_point.append(max(current_distance_array))
# print(max_face_point)
# print(len(max_face_point))

shortest_len = min(max_face_point)
print(shortest_len)
# print(len(face_point_distance_array))

for i in range(len(face_point_distance_array)):
    if face_point_distance_array[i][2] == shortest_len:
        # print(face_point_distance_array[i])
        edge_face = face_point_distance_array[i][0]
        edge_point = face_point_distance_array[i][1]
        edge_distance = face_point_distance_array[i][2]

print('edge face')
print(edge_face)
print('edge face points ')
print(hull.simplices[edge_face])
print('edge point')
print(edge_point)

print('base face coordinates')

edge_face_p1 = hull.simplices[edge_face][0]
edge_face_p2 = hull.simplices[edge_face][1]
edge_face_p3 = hull.simplices[edge_face][2]

print(hull.points[edge_face_p1], hull.points[edge_face_p2], hull.points[edge_face_p3])
print(edge_face_p1, edge_face_p2, edge_face_p3)

print('farthest point is = ')

print(hull.points[edge_point])

# print('edge face equation ')
# print(face_eq[edge_face])

    # getting the projection of the points on this plane

# this plane we know via hull simplices. z is the normal to this plane. x is the line joining the first two points of the hull simplices (hull.simplices[0][0,1])

#  standard point transformation. original coordinates are aligned with 0,0,0 and cartesian planes
# new co-ordinate system has origin at [point 0 of simplices 0],
# z prependicular to face above, x axis [joining simplices[0][point 0 nd 1] and y is normal to given plane and line or ( cross of x and line joining simplices[0][0] and highest point)

#  fist transform to now origin, simplices[0][0] and then rotate about it

new_org  = hull.points[hull.simplices[edge_face][0]]

x_v = hull.points[hull.simplices[edge_face][1]] - hull.points[hull.simplices[edge_face][0]]
x_v_unit =( x_v/np.linalg.norm(x_v)).tolist()
# print(x_v_unit)

z_v = np.cross((hull.points[hull.simplices[edge_face][2]] - hull.points[hull.simplices[edge_face][0]])  , (hull.points[hull.simplices[edge_face][1]] - hull.points[hull.simplices[edge_face][0]]))
# print(z_v)
#
#  making the z normal from plaane to point (so that point is on the top)

v_face_to_p = hull.points[edge_point] - hull.points[hull.simplices[edge_face][0]]

if np.dot(z_v, v_face_to_p) >=0:
    z_v = z_v
    # print('same directions')
else:
    z_v = -1 * (z_v)
    # print('changing directions')

z_v_unit = (z_v/np.linalg.norm(z_v)).tolist()

# print(z_v_unit)
y_v = np.cross(z_v, x_v)
y_v_unit = (np.cross(z_v_unit, x_v_unit)).tolist()

# print(y_v_unit)

new_frame = [x_v_unit, y_v_unit, z_v_unit]

old_x = [1,0,0]
old_y = [0,1,0]
old_z = [0,0,1]
old_frame = [old_x, old_y, old_z]

old_org = [0,0,0]
new_origin = hull.points[hull.simplices[edge_face][0]]

front_rot = homogenous_matrixes.homogenous_transformation_matrix(old_frame, new_frame)

# print('front rot')
# print(front_rot)
# print('front rot transpose')
# print(front_rot.transpose())
# print('front rot inverse')
# print(np.linalg.inv(front_rot))
# print(np.matmul(front_rot, np.linalg.inv(front_rot)))
# print(old_frame)
# print(new_frame)
# #
# print(np.dot(x_v, y_v))
# print(np.dot(y_v, z_v))
# print(np.dot(z_v, x_v))
#

# rotating all the points to new frame

rot_points = []
for i in range(len(hull.points)):
    curr = []
    curr = np.matmul(front_rot, hull.points[i])
    rot_points.append(curr)
#
# print(len(hull.points))
# print(len(rot_points))
# translating  all the points in the new origin [old origin to new origin]

point_rot_trans = []
for i in range(len(rot_points)):
    curr = []
    curr = homogenous_matrixes.translation(old_org, new_org, rot_points[i] )
    # curr = curr.tolist()
    point_rot_trans.append(curr)

# print(len(point_rot_trans))

# # to check the distance of all the points from the plaene if they are alomst 0 (for my given sheet )

# for i in range(len(point_rot_trans)):
#     distaa_n = abs((face_eq[edge_face][0] * point_rot_trans[i][0]) * (face_eq[edge_face][1] * point_rot_trans[i][1]) * (face_eq[edge_face][2] * point_rot_trans[i][2]) + (face_eq[edge_face][3]))
#     dista_d = ((face_eq[edge_face][0]**2) + (face_eq[edge_face][1]**2) + (face_eq[edge_face][2]**2))**0.5
#     d = distaa_n/dista_d
#     print(d)

#
# #  cutting the z of all the points

point_proj_2d = []
for i in range(len(point_rot_trans)):
    curr_xy = []
    x = point_rot_trans[i][0]
    y = point_rot_trans[i][1]
    curr_xy = [x, y]
    point_proj_2d.append(curr_xy)
# print(point_proj_2d)

xy_proj_1 = np.array(point_proj_2d)
xy_proj_1_tuple = tuple([tuple(row) for row in xy_proj_1])
# print(xy_proj_1)

#
a  = pyenvelope.get_minimum_bounding_rectangle(xy_proj_1)
#
print(a)

hull_2d = ConvexHull(xy_proj_1)
#
# # a = pyenvelope.get_minimum_bounding_rectangle(xy_proj_1)
#
rect_points  = a

#
rect_points_3d = []
for i in range(len(rect_points)):
    curr = []
    curr = [rect_points[i][0],rect_points[i][1], 0]
    rect_points_3d.append(curr)

rect_points_3d = list(rect_points_3d)
#
print(rect_points_3d)

print(len(rect_points_3d))
# tranlate back to origin

rect_points_3d_trans = []
#
for i in range(len(rect_points_3d)):
    curr = []
    curr = homogenous_matrixes.translation(new_org, old_org, rect_points_3d[i])

    rect_points_3d_trans.append(curr)
# print(len(rect_points_3d_trans))
print(rect_points_3d_trans)

# rotate all these points back to original frame

# # rotate back there 3d rectangle points
back_rot = front_rot.transpose()

final_rect_points= []
for i in range(len(rect_points_3d_trans)):
    curr = []
    curr = np.matmul(back_rot, rect_points_3d_trans[i])

    final_rect_points.append(curr)

print(final_rect_points)

original_poiints = pts
rotated_points = np.array(rot_points)
translated_points = np.array(point_rot_trans)
proj_rectangle = np.array(rect_points_3d)
proj_rectangle_trans = np.array(rect_points_3d_trans)
proj_rect_final = np.array(final_rect_points)

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

# origins and xyz vectors

o_1 = np.array(new_frame)
o_2 = np.array(old_frame)

#
# plotting origins
#
# ax.plot(o_1.T[0], o_1.T[1], o_1.T[2], "m")
# ax.plot(o_1.T[0], o_1.T[1], o_1.T[2], "m")

# Plot defining model points
ax.plot(original_poiints.T[0], original_poiints.T[1], original_poiints.T[2], "r")
ax.plot(rotated_points.T[0], rotated_points.T[1], rotated_points.T[2], "b")
ax.plot(translated_points.T[0], translated_points.T[1], translated_points.T[2], "y")

ax.plot(proj_rectangle.T[0], proj_rectangle.T[1], proj_rectangle.T[2], "y")
ax.plot(proj_rectangle_trans.T[0], proj_rectangle_trans.T[1], proj_rectangle_trans.T[2], "y")
ax.plot(proj_rect_final.T[0], proj_rect_final.T[1], proj_rect_final.T[2], "r")

# 12 = 2 * 6 faces are the simplices (2 simplices per square face)
# for s in hull.simplices:
#     s = np.append(s, s[0])  # Here we cycle back to the first coordinate
#     ax.plot(pts[s, 0], pts[s, 1], pts[s, 2], "r-")

plt.plot(xy_proj_1[:,0], xy_proj_1[:,1], 'o')
for simplex in hull.simplices:
    plt.plot(xy_proj_1[simplex,0], xy_proj_1[simplex,1], 'k-')

# plt.plot(xy_proj_1[hull_2d.vertices,0], xy_proj_1[hull_2d.vertices,1], 'r--', lw=2)
# plt.plot(xy_proj_1[hull_2d.vertices[0],0], xy_proj_1[hull_2d.vertices[0],1], 'ro')

# Make axis label
for i in ["x", "y", "z"]:
    eval("ax.set_{:s}label('{:s}')".format(i, i))

# plots making the rectangles

# ax.plot(final_box.T[0], final_box.T[1], final_box.T[2], "y")

plt.show()

max_stars_repo_path: face to point  hull.py

max_issues_repo_path: face to point  hull.py

max_forks_repo_path: face to point  hull.py

reason: import numpy,from scipy,import sympy


text:
# Copyright (c) 2020, Fabio Muratore, Honda Research Institute Europe GmbH, and
# Technical University of Darmstadt.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
#    notice, this list of conditions and the following disclaimer in the
#    documentation and/or other materials provided with the distribution.
# 3. Neither the name of Fabio Muratore, Honda Research Institute Europe GmbH,
#    or Technical University of Darmstadt, nor the names of its contributors may
#    be used to endorse or promote products derived from this software without
#    specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL FABIO MURATORE, HONDA RESEARCH INSTITUTE EUROPE GMBH,
# OR TECHNICAL UNIVERSITY OF DARMSTADT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

"""
Plot trajectories of the directly controlled joints of the Barrett WAM.
The plot always includes the trajectories from the MuJoCo simulation as baseline and optionally
    1) a trajectory recorded on the real Barrett WAM and
    2) the desired trajectories specified by the policy.
"""
import os.path as osp

import matplotlib.pyplot as plt
import numpy as np

import pyrado
from pyrado.environment_wrappers.domain_randomization import remove_all_dr_wrappers
from pyrado.logger.experiment import ask_for_experiment
from pyrado.sampling.rollout import rollout
from pyrado.utils.argparser import get_argparser
from pyrado.utils.data_types import RenderMode
from pyrado.utils.experiments import load_experiment
from pyrado.utils.input_output import print_cbt

if __name__ == "__main__":
    # Parse command line arguments
    args = get_argparser().parse_args()

    # Get the experiment's directory to load from if not given as command line argument
    ex_dir = ask_for_experiment(hparam_list=args.show_hparams) if args.dir is None else args.dir

    # Load real trajectories
    mode = ""
    real_data_exists = True
    try:
        while mode not in ["ep", "sb"]:
            mode = input("Pass ep for episodic and sb for step-based control mode: ").lower()
        qpos_real = np.load(osp.join(ex_dir, f"qpos_real_{mode}.npy"))
        qvel_real = np.load(osp.join(ex_dir, f"qvel_real_{mode}.npy"))
    except FileNotFoundError:
        real_data_exists = False
        print_cbt(
            f"Did not find a recorded real trajectory (qpos_real_{mode} and qvel_real_{mode}) for this policy. "
            f"Run deployment/run_policy_wam.py to get real-world trajectories.",
            "y",
            bright=True,
        )

    # Load the policy and the environment
    env, policy, _ = load_experiment(ex_dir, args)

    # Get nominal environment
    env = remove_all_dr_wrappers(env)
    env.domain_param = env.get_nominal_domain_param()
    env.stop_on_collision = False

    # Fix seed for reproducibility
    pyrado.set_seed(args.seed)

    # Use the recorded initial state from the real system
    init_state = env.init_space.sample_uniform()
    if real_data_exists:
        if input("Use the recorded initial state from the real system? [y] / n ").lower() == "" or "y":
            init_state[: env.num_dof] = qpos_real[0, :]

    # Define indices of actuated joints
    act_idcs = [1, 3, 5] if env.num_dof == 7 else [1, 3]

    # Do rollout in simulation
    ro = rollout(env, policy, eval=True, render_mode=RenderMode(video=False), reset_kwargs=dict(init_state=init_state))
    t = ro.env_infos["t"]
    qpos_sim, qvel_sim = ro.env_infos["qpos"], ro.env_infos["qvel"]
    qpos_des, qvel_des = ro.env_infos["qpos_des"], ro.env_infos["qvel_des"]

    plot_des = False
    if input("Plot the desired joint states and velocities (i.e. the policy features)? [y] / n ") == "" or "y":
        plot_des = True

    # Plot trajectories of the directly controlled joints and their corresponding desired trajectories
    fig, ax = plt.subplots(nrows=len(act_idcs), ncols=2, figsize=(12, 8), sharex="all", constrained_layout=True)
    fig.canvas.manager.set_window_title("Trajectory Comparison")

    # Compute the RMSE (root mean squared error) in degree
    if real_data_exists:
        err_deg = 180 / np.pi * (qpos_real - qpos_sim)
        rmse = np.sqrt(np.mean(np.power(err_deg, 2), axis=0)).round(2)
        suptitle = f"RMSE: q_1 = {rmse[1]} deg, q_3 = {rmse[3]} deg"
        if env.num_dof == 7:
            suptitle = suptitle + f", q_5 {rmse[5]} deg"
        fig.suptitle(suptitle)

    for i, idx in enumerate(act_idcs):
        ax[i, 0].plot(t, 180 / np.pi * qpos_sim[:, idx], label="sim")
        ax[i, 1].plot(t, 180 / np.pi * qvel_sim[:, idx], label="sim")

        if real_data_exists:
            ax[i, 0].plot(t, 180 / np.pi * qpos_real[:, idx], label="real")
            ax[i, 1].plot(t, 180 / np.pi * qvel_real[:, idx], label="real")

        if plot_des:
            ax[i, 0].plot(t, 180 / np.pi * qpos_des[:, idx], label="des", ls="--")
            ax[i, 1].plot(t, 180 / np.pi * qvel_des[:, idx], label="des", ls="--")

        ax[i, 0].set_ylabel(f"q_{idx} [deg]")
        ax[i, 1].set_ylabel(f"qd_{idx} [deg/s]")

    ax[0, 0].legend()
    ax[0, 1].legend()

    ax[0, 0].set_title("joint position")
    ax[0, 1].set_title("joint velocity")

    ax[len(act_idcs) - 1, 0].set_xlabel("t [s]")
    ax[len(act_idcs) - 1, 1].set_xlabel("t [s]")

    plt.show()

max_stars_repo_name: theogruner/SimuRLacra

max_stars_repo_licenses: DOC Zlib BSD-3-Clause

max_issues_repo_name: theogruner/SimuRLacra

max_issues_repo_licenses: DOC Zlib BSD-3-Clause

max_forks_repo_name: theogruner/SimuRLacra

max_forks_repo_licenses: DOC Zlib BSD-3-Clause


text:
import torch
import torch.nn as nn
import torch.nn.functional as F
import numpy as np

def divide_safe(num, den):
    eps = 1e-8
    den += eps * torch.eq(den, 0).type(torch.float)
    return num / den

def bilinear_wrapper(imgs, coords):
    '''Wrapper around bilinear sampling function

    Args:
        imgs: images to sample [..., height_s, width_s, #channels]
        coords: pixel location to sample from [..., height_t, width_t, 2]
    Returns:
        Sampled images [..., height_t, width_t, #channels]
    '''
    init_dims = imgs.shape[:-3:]
    end_dims_img = imgs.shape[-3::]
    end_dims_coords = coords.shape[-3::]
    prod_init_dims = init_dims[0]
    for i in range(1, len(init_dims)):
        prod_init_dims *= init_dims[i]
    src = imgs.contiguous().view((prod_init_dims,) + end_dims_img)
    coords = coords.view((prod_init_dims,) + end_dims_coords)
    src = src.permute([0, 3, 1, 2])
    tgt = F.grid_sample(src, coords, align_corners=False)
    tgt = tgt.permute([0, 2, 3, 1])
    tgt = tgt.view(init_dims + tgt.shape[-3::])
    return tgt

def inv_homography(k_s, k_t, rot, t, n_hat, a):
    '''Compute inverse homography matrix between two cameras via a plane.

    Args:
        k_s: source camera intrinsics [..., 3, 3]
        k_t: target camera intrinsics [..., 3, 3]
        rot: relative roation [..., 3, 3]
        t: translation from source to target camera [..., 3, 1]
        n_hat: plane normal w.r.t source camera frame [..., 1, 3]
        a: plane equation displacement [..., 1, 1]
    Returns:
        Inverse homography matrices (mapping from target to source)
        [..., 3, 3]
    '''
    rot_t = rot.transpose(-2, -1)
    k_t_inv = torch.inverse(k_t)

    denom = a - torch.matmul(torch.matmul(n_hat, rot_t), t)
    numerator = torch.matmul(torch.matmul(torch.matmul(rot_t, t),
                                          n_hat),
                             rot_t)
    inv_hom = torch.matmul(
            torch.matmul(k_s, rot_t + divide_safe(numerator, denom)),
            k_t_inv)
    return inv_hom

def transform_points(points, homography):
    '''Transforms input points according to the homography.

    Args:
        points: pixel coordinates [..., height, width, 3]
        homography: transformation [..., 3, 3]
    Returns:
        Transformed coordinates [..., height, width, 3]
    '''
    orig_shape = points.shape
    points_reshaped = points.view(orig_shape[:-3] +
                                  (-1,) +
                                  (orig_shape[-1:]))
    dim0 = len(homography.shape) - 2
    dim1 = dim0 + 1
    transformed_points = torch.matmul(points_reshaped,
                                      homography.transpose(dim0, dim1))
    transformed_points = transformed_points.view(orig_shape)
    return transformed_points

def normalize_homogenous(points):
    '''Converts homogenous coordinates to euclidean coordinates.

    Args:
        points: points in homogenous coordinates [..., #dimensions + 1]
    Returns:
        Points in standard coordinates after dividing by the last entry
        [..., #dimensions]
    '''
    uv = points[..., :-1]
    w = points[..., -1].unsqueeze(-1)
    return divide_safe(uv, w)

def transform_plane_imgs(imgs, pixel_coords, k_s, k_t, rot, t, n_hat, a):
    '''Transforms input images via homographies for corresponding planes.

    Args:
        imgs: input images [..., height_s, width_t, #channels]
        pixel_coords: pixel coordinates [..., height_t, width_t, 3]
        k_s: source camera intrinsics [..., 3, 3]
        k_t: target camera intrinsics [..., 3, 3]
        rot: relative rotation [..., 3, 3]
        t: translation from source to target camera [..., 3, 1]
        n_hat: plane normal w.r.t source camera frame [..., 1, 3]
        a: plane equation displacement [..., 1, 1]
    Returns:
        Images after bilinear sampling from the input.
    '''
    tgt2src = inv_homography(k_s, k_t, rot, t, n_hat, a)
    pixel_coords_t2s = transform_points(pixel_coords, tgt2src)
    pixel_coords_t2s = normalize_homogenous(pixel_coords_t2s)
    pixel_coords_t2s[..., 0] = pixel_coords_t2s[..., 0] /\
                               pixel_coords_t2s.shape[-2] * 2 - 1
    pixel_coords_t2s[..., 1] = pixel_coords_t2s[..., 1] /\
                               pixel_coords_t2s.shape[-3] * 2 - 1
    imgs_s2t = bilinear_wrapper(imgs, pixel_coords_t2s)

    return imgs_s2t

def planar_transform(src_imgs, pixel_coords, k_s, k_t, rot, t, n_hat, a):
    '''Transforms images, masks and depth maps according to
       planar transformation.

    Args:
        src_imgs: input images [layer, batch, height_s, width_s, #channels]
        pixel_coords: coordinates of target image pixels
                      [batch, height_t, width_t, 3]
        k_s: source camera intrinsics [batch, 3, 3]
        k_t: target camera intrinsics [batch, 3, 3]
        rot: relative rotation [batch, 3, 3]
        t: translation from source to target camera [batch, 3, 1]
        n_hat: plane normal w.r.t source camera frame [layer, batch, 1, 3]
        a: plane equation displacement [layer, batch, 1, 1]
    Returns:
        Images projected to target frame [layer, height, width, #channels]
    '''
    layer = src_imgs.shape[0]
    rot_rep_dims = [layer]
    rot_rep_dims += [1 for _ in range(len(k_s.shape))]

    cds_rep_dims = [layer]
    cds_rep_dims += [1 for _ in range(len(pixel_coords.shape))]

    k_s = k_s.repeat(rot_rep_dims)
    k_t = k_t.repeat(rot_rep_dims)
    t = t.repeat(rot_rep_dims)
    rot = rot.repeat(rot_rep_dims)
    pixel_coords = pixel_coords.repeat(cds_rep_dims)

    tgt_imgs = transform_plane_imgs(
            src_imgs, pixel_coords, k_s, k_t, rot, t, n_hat, a)

    return tgt_imgs


text:
from collections import namedtuple

import numpy
import pytest

from thyme.model.fvcom import node_to_centroid
from thyme.model.fvcom import vertical_interpolation

VerticalValues = namedtuple(
    'VerticalValues',
    ['u',
     'v',
     'h',
     'zeta',
     'siglay_centroid',
     'num_nele',
     'num_siglay',
     'time_index',
     'target_depth_default',
     'target_depth_surface',
     'target_depth_deep',
     'expected_u_target_depth_default',
     'expected_v_target_depth_default',
     'expected_u_target_depth_surface',
     'expected_v_target_depth_surface',
     'expected_u_target_depth_deep',
     'expected_v_target_depth_deep'])

@pytest.fixture
def vertical_values():

    time_index = 0
    num_nele = 7
    num_siglay = 20

    target_depth_default = 4.5
    target_depth_surface = 0
    target_depth_deep = 15

    u = numpy.array(
        [[

            [-0.7951693, 0.1987104, 0.20282207, 0.22636837, 0.22338444, 0.22267905, 0.25851583],
            [-0.79512405, 0.17984135, 0.18718982, 0.20845266, 0.20881361, 0.2097686, 0.24099022],
            [-0.79499704, 0.11546519, 0.11235027, 0.12845743, 0.14677778, 0.17036788, 0.14718369],
            [-0.7947852, -0.04345991, -0.0362343, -0.01943639, 0.03145187, 0.05482485, 0.01756552],
            [-0.794485, -0.11175916, -0.10817654, -0.07648413, -0.05387694, -0.02363495, -0.07206099],
            [-0.79409266, -0.15002064, -0.1458673, -0.11296615, -0.09673776, -0.07272214, -0.11500881],
            [-0.7936041, -0.18279588, -0.17216434, -0.15122096, -0.13439055, -0.12422991, -0.15305918],
            [-0.7930144, -0.22257489, -0.21475703, -0.19943172, -0.19317457, -0.18883257, -0.19185503],
            [-0.7923182, -0.26589555, -0.26381484, -0.2563959, -0.2623669, -0.2418433, -0.22976391],
            [-0.7915093, -0.29618508, -0.3025776, -0.30066454, -0.30951524, -0.2785821, -0.2603346],
            [-0.62858844, -0.31844422, -0.3287656, -0.32818496, -0.34010696, -0.2991507, -0.28074494],
            [-0.59509104, -0.3200005, -0.3282557, -0.33371437, -0.3424295, -0.3079574, -0.2808675],
            [-0.5648844, -0.28674564, -0.29246253, -0.3045539, -0.3127165, -0.2805328, -0.26253417],
            [-0.53362864, -0.2308568, -0.23583505, -0.25023454, -0.25649628, -0.2430087, -0.21772045],
            [-0.49928024, -0.1515698, -0.17267239, -0.1813446, -0.19186072, -0.19516239, -0.16393802],
            [-0.45988485, -0.13119768, -0.10612836, -0.11343177, -0.12980227, -0.1302748, -0.12355448],
            [-0.41265842, -0.11103977, -0.09538705, -0.09509068, -0.11300715, -0.10742138, -0.10589793],
            [-0.35249037, -0.09441119, -0.0877323, -0.08854523, -0.10422827, -0.09832109, -0.098252],
            [-0.26714894, -0.08462795, -0.0790456, -0.07945603, -0.09364957, -0.08832049, -0.08927543],
            [-0.10820818, -0.06875627, -0.06421904, -0.06397585, -0.07394486, -0.07045165, -0.07252091]

        ]])

    v = numpy.array(
        [[
            [-0.6447422, 0.04272371, 0.03626542, 0.03629172, 0.02365562, 0.04447182, 0.02940587],
            [-0.64464253, 0.04550566, 0.03864311, 0.03937801, 0.02540983, 0.04566612, 0.03169739],
            [-0.6444702, 0.05464792, 0.05194416, 0.05265207, 0.03223358, 0.04813668, 0.0406685],
            [-0.6442224, 0.05445727, 0.05883322, 0.04983669, 0.0348991, 0.03746793, 0.04016836],
            [-0.64389604, 0.00415152, -0.00149629, -0.00077556, 0.01316911, 0.00495433, 0.01126673],
            [-0.6434879, -0.04056849, -0.0399819, -0.02377908, -0.01824965, -0.0271148, -0.02322751],
            [-0.64299417, -0.05177833, -0.05065423, -0.04175215, -0.03915307, -0.05457062, -0.0499136],
            [-0.6424107, -0.0646683, -0.0626291, -0.05882293, -0.06086947, -0.07759117, -0.06596889],
            [-0.64173275, -0.06379095, -0.06311549, -0.06392502, -0.07486805, -0.08554274, -0.06285658],
            [-0.6409548, -0.05030349, -0.05021876, -0.05602735, -0.06358989, -0.07058521, -0.04798502],
            [-0.48587024, -0.02668872, -0.02443161, -0.03338933, -0.03849523, -0.04692423, -0.03141288],
            [-0.45516193, -0.0020759, 0.00556405, -0.00915105, -0.01343074, -0.02884087, -0.01940906],
            [-0.4284155, 0.0134382, 0.01965301, 0.00118551, -0.00354341, -0.0235792, -0.01681828],
            [-0.40160605, 0.00625094, 0.0087254, -0.00298712, -0.01015626, -0.03078334, -0.02404925],
            [-0.3729973, -0.02794126, -0.02092287, -0.02809138, -0.03298125, -0.04707557, -0.04252515],
            [-0.34107724, -0.02949184, -0.03453156, -0.03266618, -0.04048001, -0.05651411, -0.03880842],
            [-0.30380374, -0.02897151, -0.0317935, -0.02676859, -0.03367185, -0.05099317, -0.03234524],
            [-0.2575078, -0.02654107, -0.02930111, -0.0243634, -0.03028261, -0.04743995, -0.02939028],
            [-0.19346614, -0.02417713, -0.0263748, -0.02134923, -0.02650077, -0.04303462, -0.02615478],
            [-0.07716001, -0.01982844, -0.02149794, -0.01682694, -0.02036272, -0.03471823, -0.02087618]

        ]])

    h = numpy.array(
        [

            1.372000018755595, 1.91399347, 4.6028045, 3.29994912, 30.32372402, 9.3237240, 12.19253056

        ])

    zeta = numpy.array(
        [
            -1.1558676163355508, -0.1443424, -0.14408446, -0.15902775, -0.16805193, -0.16938841, -0.16457143,

        ])

    siglay_centroid = numpy.array(
        [

            [-0.025, -0.02500004, -0.02500004, -0.02500004, -0.02500004, -0.02500004, -0.02500004],
            [-0.075, -0.07500005, -0.07500005, -0.07500005, -0.07500005, -0.07500005, -0.07500005],
            [-0.125, -0.12500003, -0.12500003, -0.12500003, -0.12500003, -0.12500003, -0.12500003],
            [-0.17500001, -0.17500001, -0.17500001, -0.17500001, -0.17500001, -0.17500001, -0.17500001],
            [-0.225, -0.22500002, -0.22500002, -0.22500002, -0.22500002, -0.22500002, -0.22500002],
            [-0.275, -0.27500007, -0.27500007, -0.27500007, -0.27500007, -0.27500007, -0.27500007],
            [-0.325, -0.32500005, -0.32500005, -0.32500005, -0.32500005, -0.32500005, -0.32500005],
            [-0.375, -0.37500003, -0.37500003, -0.37500003, -0.37500003, -0.37500003, -0.37500003],
            [-0.425, -0.42500004, -0.42500004, -0.42500004, -0.42500004, -0.42500004, -0.42500004],
            [-0.475, -0.47500002, -0.47500002, -0.47500002, -0.47500002, -0.47500002, -0.47500002],
            [-0.525, -0.52500004, -0.52500004, -0.52500004, -0.52500004, -0.52500004, -0.52500004],
            [-0.57500005, -0.57500005, -0.57500005, -0.57500005, -0.57500005, -0.57500005, -0.57500005],
            [-0.625, -0.6250001, -0.6250001, -0.6250001, -0.6250001, -0.6250001, -0.6250001],
            [-0.67499995, -0.6750001, -0.6750001, -0.6750001, -0.6750001, -0.6750001, -0.6750001],
            [-0.725, -0.725, -0.725, -0.725, -0.725, -0.725, -0.725],
            [-0.775, -0.7750001, -0.7750001, -0.7750001, -0.7750001, -0.7750001, -0.7750001],
            [-0.82500005, -0.82500005, -0.82500005, -0.82500005, -0.82500005, -0.82500005,  -0.82500005],
            [-0.875, -0.87500006, -0.87500006, -0.87500006, -0.87500006, -0.87500006, -0.87500006],
            [-0.92499995, -0.9250001, -0.9250001, -0.9250001, -0.9250001, -0.9250001,  -0.9250001],
            [-0.975, -0.975, -0.975, -0.975, -0.975, -0.975, -0.975]

        ])

    expected_u_target_depth_default = numpy.array(
        [
            -0.7100488057968346, -0.30731464, -0.31567158, -0.31442473, 0.09090091, -0.28539867, -0.19117865
        ])

    expected_v_target_depth_default = numpy.array(
        [
            -0.5634124710903746, -0.03849612, -0.0373252, -0.04470835,  0.03352506, -0.06274381, -0.06568897
        ])

    expected_u_target_depth_surface = numpy.array(
        [
            -0.79519192, 0.20814494, 0.21063821, 0.23532624, 0.23066987, 0.22913428, 0.26727865
        ])

    expected_v_target_depth_surface = numpy.array(
        [
            -0.64479204, 0.04133273, 0.03507657, 0.03474857, 0.02277851, 0.04387467, 0.02826011
        ])

    expected_u_target_depth_deep = numpy.array(
        [
            -0.71004887, -0.30731464, -0.31567158, -0.31442473, -0.32323185, -0.28886639, -0.27053976
        ])

    expected_v_target_depth_deep = numpy.array(
        [
            -0.56341252, -0.03849612, -0.0373252, -0.04470835, -0.05233803, -0.05875473, -0.03969896
        ])

    return VerticalValues(u, v, h, zeta, siglay_centroid, num_nele, num_siglay, time_index, target_depth_default,
                          target_depth_surface, target_depth_deep, expected_u_target_depth_default,
                          expected_v_target_depth_default, expected_u_target_depth_surface,
                          expected_v_target_depth_surface, expected_u_target_depth_deep, expected_v_target_depth_deep)

def test_vertical_interpolation(vertical_values):
    """Test vertical interpolation."""
    u_target_depth, v_target_depth = vertical_interpolation(vertical_values.u, vertical_values.v, vertical_values.h,
                                                            vertical_values.zeta, vertical_values.siglay_centroid,
                                                            vertical_values.num_nele, vertical_values.num_siglay,
                                                            vertical_values.time_index,
                                                            vertical_values.target_depth_default)

    # print(f"u_target_depth: {u_target_depth}")
    # print(f"v_target_depth: {v_target_depth}")
    assert numpy.allclose(u_target_depth, vertical_values.expected_u_target_depth_default)
    assert numpy.allclose(v_target_depth, vertical_values.expected_v_target_depth_default)

def test_vertical_interpolation_at_surface(vertical_values):
    """Test vertical interpolation."""
    u_target_depth, v_target_depth = vertical_interpolation(vertical_values.u, vertical_values.v, vertical_values.h,
                                                            vertical_values.zeta, vertical_values.siglay_centroid,
                                                            vertical_values.num_nele, vertical_values.num_siglay,
                                                            vertical_values.time_index,
                                                            vertical_values.target_depth_surface)

    # print(f"u_target_depth_surface: {u_target_depth}")
    # print(f"v_target_depth_surface: {v_target_depth}")
    assert numpy.allclose(u_target_depth, vertical_values.expected_u_target_depth_surface)
    assert numpy.allclose(v_target_depth, vertical_values.expected_v_target_depth_surface)

def test_vertical_interpolation_deep(vertical_values):
    """Test vertical interpolation."""
    u_target_depth, v_target_depth = vertical_interpolation(vertical_values.u, vertical_values.v, vertical_values.h,
                                                            vertical_values.zeta, vertical_values.siglay_centroid,
                                                            vertical_values.num_nele, vertical_values.num_siglay,
                                                            vertical_values.time_index,
                                                            vertical_values.target_depth_deep)

    # print(f"u_target_depth_deep: {u_target_depth}")
    # print(f"v_target_depth_deep: {v_target_depth}")
    assert numpy.allclose(u_target_depth, vertical_values.expected_u_target_depth_deep)
    assert numpy.allclose(v_target_depth, vertical_values.expected_v_target_depth_deep)

NodeToCentroidValues = namedtuple(
    'NodeToCentroidValues',
    ['zeta',
     'h',
     'lat_node',
     'lon_node',
     'lat_centroid',
     'lon_centroid',
     'time_index',
     'expected_h_centroid',
     'expected_zeta_centroid'])

@pytest.fixture
def node_to_centroid_values():
    zeta = numpy.array(
        [
            [-0.146551, -0.146383, -0.145782, -0.146226, -0.145688, -0.145071, -0.146047, -0.145581, -0.145099,
             -0.144439, -0.145729, -0.145431, -0.144983, -0.144627, -0.144069, -0.145415, -0.145071, -0.144579,
             -0.144391, -0.143806, -0.144057,

             -0.169269, -0.169959, -0.170395, -0.170738, -0.166229, -0.166729, -0.167032, -0.162118, -0.162876,
             -0.163448, -0.163812, -0.158493, -0.159403, -0.159949, -0.160389
             ]
        ])

    h = numpy.array(
        [
            1.6438, 1.8724, 1.875, 2.0302, 1.9946, 2.0446, 2.1169, 2.1307, 2.1457, 2.2134, 2.1768, 2.178, 2.2341,
            2.2623, 2.331, 2.1949, 2.2047, 2.2817, 2.3537, 2.434, 2.404,

            31.461, 31.874, 32.37, 32.361, 30.565, 31.089, 31.196, 29.231, 29.033, 28.687, 28.761, 26.053, 25.788,
            26.267, 26.835
        ])
    lat_centroid = numpy.array(
        [
            30.260946, 30.262287, 30.260996, 30.259796, 30.26362, 30.262392, 30.261309, 30.260002, 30.258938, 30.264635,
            30.263693, 30.262672, 30.261517, 30.260561, 30.259289, 30.258574, 30.264799, 30.264191, 30.263008,
            30.262444, 30.261324, 30.260527, 30.259325, 30.259058, 30.260839, 30.259773,

            29.969425, 29.963093, 29.967611, 29.961212, 29.96579, 29.959322, 29.980326, 29.985025,
            29.978645, 29.983408, 29.976954, 29.981771, 29.996161, 30.001049, 29.994678, 29.999624,
            29.993172, 29.998173
        ])

    lon_centroid = numpy.array(
        [
            -88.126221, -88.124329, -88.124817, -88.123871, -88.122406, -88.122925, -88.12204, -88.122498, -88.12149,
            -88.120361, -88.120972, -88.120026, -88.120575, -88.119843, -88.1203, -88.119263, -88.118958, -88.118195,
            -88.118652, -88.117645, -88.117676, -88.118591, -88.118286, -88.117096, -88.116516, -88.116272,

            -88.264221, -88.25528, -88.245392, -88.236328, -88.226379, -88.217194, -88.265015,
            -88.256958, -88.246277, -88.238129, -88.227295, -88.219116, -88.260498, -88.253448,
            -88.241791, -88.23468, -88.222839, -88.215698
        ])

    lat_node = numpy.array(
        [
            30.260969, 30.262211, 30.259659, 30.263531, 30.26112, 30.258608, 30.264805, 30.262526, 30.260281, 30.257925,
            30.265354, 30.263748, 30.261744, 30.259661, 30.258137, 30.265295, 30.263531, 30.262054, 30.260176, 30.25886,
            30.260286,

            29.958637, 29.956696, 29.954746, 29.952793, 29.973944, 29.972195, 29.970428, 29.991339, 29.989796,
            29.988235, 29.986652, 30.007345, 30.006006, 30.004629, 30.003235
        ])

    lon_node = numpy.array(
        [
            -88.127625, -88.125763, -88.125275, -88.12384, -88.123413, -88.122894, -88.121857, -88.12149, -88.121185,
            -88.120422, -88.119629, -88.119568, -88.119049, -88.119293, -88.118073, -88.117676, -88.11734, -88.116516,
            -88.117462, -88.115784, -88.11557,

            -88.265167, -88.246277, -88.227234, -88.207947, -88.254395, -88.235504, -88.2164, -88.267578,
            -88.248901, -88.230011, -88.210937, -88.265045, -88.246429, -88.2276, -88.208527
        ])

    time_index = 0

    expected_h_centroid = numpy.array(
        [
            1.79706003, 1.96634453, 1.91399347, 1.97068808, 2.09215354, 2.05125885, 2.08959921, 2.06159611, 2.13525152,
            2.15697391, 2.14185806, 2.1812977, 2.17015265, 2.2140166,  2.20713333, 2.26888473, 2.18323025, 2.19253056,
            2.20559465, 2.23989205, 2.28985895, 2.28394075, 2.31515434, 2.37337259, 2.3464489, 2.39725858,

            30.74697778, 31.29994912, 31.1760438, 31.77788929, 31.55164158, 31.9756924, 30.00374208, 29.60978718,
            30.22881187, 29.6028045, 30.32372402, 29.54809818, 28.10520577, 26.95770733, 27.83656214, 26.91390017,
            27.90491454, 27.28733678

        ])

    expected_zeta_centroid = numpy.array(
        [
           -0.1462386, -0.14609691, -0.14595087, -0.1455159, -0.14595329, -0.14583367, -0.14545806, -0.14528592,
           -0.14486764, -0.14573735, -0.1456864, -0.14532972, -0.14522106, -0.14490293, -0.14472167, -0.14437832,
           -0.14552505, -0.14530566, -0.14516168, -0.14487977, -0.14465093, -0.14466511, -0.1443643, -0.14408686,
           -0.1443424, -0.14408446,

           -0.16695492, -0.16848552, -0.16763914, -0.16902775, -0.16805193, -0.16938841, -0.16457143, -0.16374129,
           -0.16527793, -0.16435127, -0.16573603, -0.16476413, -0.16116232, -0.16025768, -0.16190908, -0.16093319,
           -0.16240301, -0.16138305
        ])

    return NodeToCentroidValues(zeta, h, lat_node, lon_node, lat_centroid, lon_centroid, time_index,
                                expected_h_centroid, expected_zeta_centroid)

def test_node_to_centroid(node_to_centroid_values):
    """Test horizontal interpolation from nodes to centroids(elements)."""
    h_centroid, zeta_centroid = node_to_centroid(node_to_centroid_values.zeta, node_to_centroid_values.h,
                                                 node_to_centroid_values.lon_node, node_to_centroid_values.lat_node,
                                                 node_to_centroid_values.lon_centroid,
                                                 node_to_centroid_values.lat_centroid,
                                                 node_to_centroid_values.time_index)

    # print(f"h_centroid: {h_centroid}")
    # print(f"zeta_centroid: {zeta_centroid}")
    assert numpy.allclose(h_centroid, node_to_centroid_values.expected_h_centroid)
    assert numpy.allclose(zeta_centroid, node_to_centroid_values.expected_zeta_centroid)

max_stars_repo_name: noaa-ocs-modeling/thyme

max_issues_repo_name: noaa-ocs-modeling/thyme

max_forks_repo_name: noaa-ocs-modeling/thyme


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

Governing equations: Double porosity/permeability
Formulation: mixed RT0 formulation
Problem: 2D convergence - Square domain
Run as:
  python <>.py <rtol> <show_monitor> <solver_type> <amg_type> <nx>

  where:
  <rtol> = relative solver tolerance
  <show_monitor> = Turn off KSP residual monitor (0), turn on (1)
  <solver_type> = naive (1), split by fields ("fields"), split by scales ("scales")
  <amg_type> = algebraic multigrid type, either gamg, hypre, or ml
  <nx> = number of cells in each spatial direction
  <element> = T3/Q4

  I recommend starting off with:
  <rtol> = 1e-5
  <show_monitor> = 1
  <solver_type> = fields or scales
  <amg_type> = hypre
  <nx> = 10
  <element> = T3

  Play around with these as necessary. Few notes:
  - Turn off monitor to improve timing metrics
  - Penaties eta_u and eta_p affect performances
  - Solver count will fluctuate with number of MPI processes
  - If solving 2D problems, comment out the hypre boomeramg parameters
  - Preliminary runs suggest gamg is the worst amg_type
  - Preliminary runs suggest that splitting by scales is the best
"""
from firedrake import *
import numpy as np
from pyop2.profiling import timed_region
from mpi4py import MPI
import math,sys,time
rank = MPI.COMM_WORLD.Get_rank()
size = MPI.COMM_WORLD.Get_size()
rtol = float(sys.argv[1])
show_monitor = int(sys.argv[2])
solvestrategy = sys.argv[3]
amg_type = sys.argv[4]
nx = int(sys.argv[5])
element = sys.argv[6]
#=====================================;
#  Create mesh and identify boundary  ;
#=====================================;
if element == 'T3':
	mesh = UnitSquareMesh(nx,nx)
else:
	mesh = RectangleMesh(nx,nx,1,1,quadrilateral = True)
#====================================================;
#  Define function spaces and mixed (product) space  ;
#====================================================;
#velocity
if element == 'T3':
	vSpace = FunctionSpace(mesh,"RT",1)
else:
	vSpace = FunctionSpace(mesh, "RTCF", 1)
#pressure
pSpace = FunctionSpace(mesh,"DG",0)
wSpace = MixedFunctionSpace([vSpace,pSpace,vSpace,pSpace])

#===================================;
#  Define trial and test functions  ;
#===================================;
(v1,p1,v2,p2) = TrialFunctions(wSpace)
(w1,q1,w2,q2) = TestFunctions(wSpace)

#======================;
#  Define body forces  ;
#======================;
rhob1, rhob2 = Constant((0.0,0.0)), Constant((0.0,0.0))

#============================;
#  Define medium properties  ;
#============================;
mu = Constant(1.0)
beta = Constant(1.0)
fact = 1.0  #  fact = beta / mu
k1, k2 = Constant(1.0), Constant(0.1)
alpha1, alpha2 = Constant(mu/k1), Constant(mu/k2)

eta = Constant(3.316625)

#----------------;
#  Pressure BCs  ;
#----------------;
x, y = SpatialCoordinate(mesh)
p1_left = Function(pSpace).interpolate((1/pi)*sin(pi*y) - exp(3.316625*y))
p1_right = Function(pSpace).interpolate((1/pi)*exp(pi)*sin(pi*y) - exp(3.316625*y))
p1_bottom = Constant(-1.0)
p1_top = Constant(-27.5671484)
p2_left = Function(pSpace).interpolate((1/pi)*sin(pi*y) + 10.0 *exp(3.316625*y))
p2_right = Function(pSpace).interpolate((1/pi)*exp(pi)*sin(pi*y) + 10.0*exp(3.316625*y))
p2_bottom = Constant(10.0)
p2_top = Constant(275.671484)

bcs = []
#=====================================;
# Define normal vector and mesh size  ;
#=====================================;
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2

#===========================;
#  Define variational form  ;
#===========================;
a = dot(w1, alpha1*v1)*dx + dot(w2, alpha2*v2)*dx \
    - div(w1) * p1 *dx - div(w2) * p2 * dx \
    + q1 * div(v1) * dx + q2 * div(v2) * dx +\
    q1 * fact * (p1 - p2) * dx -\
    q2 * fact * (p1 - p2) * dx

L = dot(w1,rhob1)*dx +\
    dot(w2,rhob2)*dx -\
    dot(w1,n) * p1_left * ds(1) -\
    dot(w2,n) * p2_left * ds(1) -\
    dot(w1,n) * p1_right * ds(2) -\
    dot(w2,n) * p2_right * ds(2) -\
    dot(w1,n) * p1_bottom * ds(3) -\
    dot(w2,n) * p2_bottom * ds(3) -\
    dot(w1,n) * p1_top * ds(4) -\
        dot(w2,n) * p2_top * ds(4)
#==================;
#  Solver options  ;
#==================;
if show_monitor:
  monitor_view = True
else:
  monitor_view = False
parameters_standard = {
  "ksp_type": "gmres",
  "pc_type": "bjacobi",
  "ksp_monitor_true_residual": show_monitor,
  "ksp_converged_reason": True,
  "ksp_rtol": rtol
}
parameters_twofields = {
  "ksp_type": "gmres",
  #"mat_type": "aij",
  "pc_type": "fieldsplit",
  # first split contains velocities
  # second split contains pressures
  "pc_fieldsplit_0_fields": "0,2",
  "pc_fieldsplit_1_fields": "1,3",
  # Use schur complement
  "pc_fieldsplit_type": "schur",
  "pc_fieldsplit_schur_fact_type": "full",
  "pc_fieldsplit_schur_precondition": "selfp",
  # Velocities
  "fieldsplit_0_ksp_type": "preonly",
  "fieldsplit_0_pc_type": "bjacobi",
  # Pressures - fieldsplit into scales
  "fieldsplit_1_ksp_type": "preonly",
  "fieldsplit_1_pc_type": "fieldsplit",
  "fieldsplit_1_pc_fieldsplit_type": "additive",
  # Macro scale pressure
  "fieldsplit_1_fieldsplit_0_ksp_type": "preonly",
  "fieldsplit_1_fieldsplit_0_pc_type": amg_type,
  # Only needed if using hypre for 3D problems
  # Comment out for 2D problems
  #"fieldsplit_1_fieldsplit_0_pc_hypre_boomeramg_strong_threshold": 0.75,
  #"fieldsplit_1_fieldsplit_0_pc_hypre_boomeramg_agg_nl": 2,
  # Micro scale pressure
  "fieldsplit_1_fieldsplit_1_ksp_type": "preonly",
  "fieldsplit_1_fieldsplit_1_pc_type": amg_type,
  # Only needed if using hypre for 3D problems
  # Comment out for 2D problems
  #"fieldsplit_1_fieldsplit_1_pc_hypre_boomeramg_strong_threshold": 0.75,
  #"fieldsplit_1_fieldsplit_1_pc_hypre_boomeramg_agg_nl": 2,
  "ksp_monitor_true_residual": monitor_view,
  "ksp_converged_reason": True,
  "ksp_rtol": rtol
}
parameters_twoscales = {
  "ksp_type": "gmres",
  #"mat_type": "aij",
  "pc_type": "fieldsplit",
  # First split contains macro scale
  # Second split contains micro scale
  "pc_fieldsplit_0_fields": "0,1",
  "pc_fieldsplit_1_fields": "2,3",
  "pc_fieldsplit_type": "additive", # Can change this to multiplicative
  # Macro scale Darcy
  "fieldsplit_0_ksp_type": "preonly",
  "fieldsplit_0_pc_type": "fieldsplit",
  "fieldsplit_0_pc_fieldsplit_type": "schur",
  "fieldsplit_0_pc_fieldsplit_schur_fact_type": "full",
  "fieldsplit_0_pc_fieldsplit_schur_precondition": "selfp",
  # Fieldsplitting of macro scale
  "fieldsplit_0_fieldsplit_0_ksp_type": "preonly",
  "fieldsplit_0_fieldsplit_0_pc_type": "bjacobi",
  "fieldsplit_0_fieldsplit_1_ksp_type": "preonly",
  "fieldsplit_0_fieldsplit_1_pc_type": amg_type,
  # Only needed if using hypre for 3D problems
  # Comment out for 2D problems
  #"fieldsplit_0_fieldsplit_1_pc_hypre_boomeramg_strong_threshold": 0.75,
  #"fieldsplit_0_fieldsplit_1_pc_hypre_boomeramg_agg_nl": 2,
  # Micro scale Darcy
  "fieldsplit_1_ksp_type": "preonly",
  "fieldsplit_1_pc_type": "fieldsplit",
  "fieldsplit_1_pc_fieldsplit_type": "schur",
  "fieldsplit_1_pc_fieldsplit_schur_fact_type": "full",
  "fieldsplit_1_pc_fieldsplit_schur_precondition": "selfp",
  # Fieldsplitting of micro scale
  "fieldsplit_1_fieldsplit_0_ksp_type": "preonly",
  "fieldsplit_1_fieldsplit_0_pc_type": "bjacobi",
  "fieldsplit_1_fieldsplit_1_ksp_type": "preonly",
  "fieldsplit_1_fieldsplit_1_pc_type": amg_type,
  # Only needed if using hypre for 3D problems
  # Comment out for 2D problems
  #"fieldsplit_1_fieldsplit_1_pc_hypre_boomeramg_strong_threshold": 0.75,
  #"fieldsplit_1_fieldsplit_1_pc_hypre_boomeramg_agg_nl": 2,
  "ksp_monitor_true_residual": monitor_view,
  "ksp_converged_reason": True,
  "ksp_rtol": rtol
}

#=================;
#  Solve problem  ;
#=================;
solution = Function(wSpace)
with timed_region('solve_firedrake'):
  initialtime = time.time()
  A = assemble(a, bcs=bcs, mat_type='aij')
  b = assemble(L)
  assembletime = time.time()
  # Option 1: default solver
  if solvestrategy == "naive":
    solver = LinearSolver(A,P=None,options_prefix="standard_",solver_parameters=parameters_standard)
  # Option 2: split by fields
  elif solvestrategy == "fields":
    solver = LinearSolver(A,P=None,options_prefix="twofields_",solver_parameters=parameters_twofields)
  # Option 3: split by scales
  elif solvestrategy == "scales":
    solver = LinearSolver(A,P=None,options_prefix="twoscales_",solver_parameters=parameters_twoscales)
  solver.solve(solution,b)
  solvetime = time.time()

#=======================;
#  Performance metrics  ;
#=======================;
with solution.dat.vec as solution_vec:
  DoF = solution_vec.getSize()
totaltime = solvetime - initialtime
dofsec = DoF/totaltime
dofsecLOG = np.log10(dofsec)
totaltimeLOG = np.log10(totaltime)

if rank == 0:
  print("================================")
  print("%d MPI processes" % size)
  print("AMG and solver type: %s %s" % (amg_type,solvestrategy))
  print("Relative residual: %1.1e" % (rtol))
  print("Nx = %d" % nx)
  print("Dof = %d\n" % DoF)
  print("Assembly time = %1.3e seconds" % (assembletime - initialtime))
  print("Solve time = %1.3e seconds" % (solvetime - assembletime))
  print("Totaltime = %1.3e seconds" % totaltime)
  print("DofSec = %1.3e" % dofsec)
  print("DofSecLOG = %1.3e" % dofsecLOG) # used for static scaling
  print("TotaltimeLOG = %1.3e log(seconds)" % totaltimeLOG) # used for static scaling , True static scaling, and DoE plots
  #print("L2 error = %1.3e" % L2error) # True static-scaling
  #print("TrueDof = %1.3e" % truedof) # True static-scaling
  #print("DoE = %1.3e" % DoE) # Digits of efficacy
  print("================================")
#=======================================;
#  Dump solution to file in VTK format  ;
#=======================================;
v1sol,p1sol,v2sol,p2sol = solution.split()

file = File('Output/v1_2D_RT0.pvd')
file.write(v1sol)

file = File('Output/v2_2D_RT0.pvd')
file.write(v2sol)

file = File('Output/p1_2D_RT0.pvd')
file.write(p1sol)

file = File('Output/p2_2D_RT0.pvd')
file.write(p2sol)

#==========================;
#  Define exact solutions  ;
#==========================;
p1_ex = Function(pSpace)
p2_ex = Function(pSpace)
v1_ex = Function(vSpace)
v2_ex = Function(vSpace)
p1_exact = Expression("(1/pi)*exp(pi*x[0])*sin(pi*x[1]) - (1/(1.0*1.0))*exp(3.316625*x[1])", degree = 5)
p2_exact = Expression("(1/pi)*exp(pi*x[0])*sin(pi*x[1]) + (1/(1.0*0.1))*exp(3.316625*x[1])", degree = 5)
v1_exact = Expression(("-1*exp(pi*x[0])*sin(pi*x[1])","-1*exp(pi*x[0])*cos(pi*x[1]) + (3.316625)*exp(3.316625*x[1])"), degree = 5)
v2_exact = Expression(("-0.1*exp(pi*x[0])*sin(pi*x[1])","-0.1*exp(pi*x[0])*cos(pi*x[1]) - (3.316625)*exp(3.316625*x[1])"), degree = 5)
# interpolate and project are the same.
p1_ex = project(p1_exact, pSpace)
v1_ex = project(v1_exact, vSpace)
p2_ex = project(p2_exact, pSpace)
v2_ex = project(v2_exact, vSpace)

#file = File('Output/v1ex_2D.pvd')
#file.write(v1_ex)

#file = File('Output/v2ex_2D.pvd')
#file.write(v2_ex)

#file = File('Output/p1ex_2D.pvd')
#file.write(p1_ex)

#file = File('Output/p2ex_2D.pvd')
#file.write(p2_ex)

##plot(p2_ex,title = "Exact solution_microP")
L2_p1 = errornorm(p1_ex,p1sol,norm_type='L2',degree_rise= 3)
#H1_p1 = errornorm(p1_ex,p1sol,norm_type='H1',degree_rise=3)
L2_v1 = errornorm(v1_ex,v1sol,norm_type='L2',degree_rise= 3)
#print ("H1 error in p1 = ", H1_p1)
L2_p2 = errornorm(p2_ex,p2sol,norm_type='L2',degree_rise= 3)
#H1_p2 = errornorm(p2_ex,p2sol,norm_type='H1',degree_rise=3)
L2_v2 = errornorm(v2_ex,v2sol,norm_type='L2',degree_rise= 3)
#print ("H1 error in p2 = ", H1_p2)
#L2error = ...
#-------------------;
# CONVERGENCE PLOTS ;
#-------------------;
#DoA = -np.log10(L2error)
DoA_p1 = -np.log10(L2_p1)
DoA_p2 = -np.log10(L2_p2)
DoA_v1 = -np.log10(L2_v1)
DoA_v2 = -np.log10(L2_v2)
DoS = np.log10(DoF)
#--------------;
# TRUE SCALING ;
#--------------;
#truedof = DoA/DoS*dofsec
truedof_p1 = DoA_p1/DoS*dofsec
truedof_p2 = DoA_p2/DoS*dofsec
truedof_v1 = DoA_v1/DoS*dofsec
truedof_v2 = DoA_v2/DoS*dofsec
#-----;
# DoE ;
#-----;
#DoE = -np.log10(L2error*totaltime)
DoE_p1 = -np.log10(L2_p1*totaltime)
DoE_p2 = -np.log10(L2_p2*totaltime)
DoE_v1 = -np.log10(L2_v1*totaltime)
DoE_v2 = -np.log10(L2_v2*totaltime)

if rank == 0:
	print ("L2 error in p1 =", L2_p1)
	print ("L2 error in v1 =", L2_v1)
	print ("L2 error in p2 =", L2_p2)
	print ("L2 error in v2 =", L2_v2)
	print ("DoA_p1 =" , DoA_p1)
	print ("DoA_p2 =" , DoA_p2)
	print ("DoA_v1 =" , DoA_v1)
	print ("DoA_v2 =" , DoA_v2)
	print ("DoS =" , DoS)
	print ("truedof_p1 =" , truedof_p1)
	print ("truedof_p2 =" , truedof_p2)
	print ("truedof_v1 =" , truedof_v1)
	print ("truedof_v2 =" , truedof_v2)
	print ("DoE_p1 =" , DoE_p1)
	print ("DoE_p2 =" , DoE_p2)
	print ("DoE_v1 =" , DoE_v1)
	print ("DoE_v2 =" , DoE_v2)


text:
from numpy import *
import operator

def createDataSet():
    group = array( [[1.0, 1.1], [1.0, 1.0], [0,0], [0, 0.1] ])
    labels = ['A', 'A', 'B', 'B']
    return group, labels

def classify0(inX, dataSet, lables, k):
   dataSetSize = dataSet.shape[0]
   diffMat = tile(inX, (dataSetSize,1)) - dataSet
   sqDiffMat = diffMat ** 2
   sqDistances = sqDiffMat.sum(axis=1)
   distances = sqDistances ** 0.5
   sortedDistIndicies = distances.argsort()
   classCount = {}
   for i in range(k):
      voteIlabel = labels[sortedDistIndicies[i]]
      classCount[voteIlable] = classCount.get(voteIlabel,0) + 1
   sortedClassCount = sorted(classCount.iteritems(), key=operator.itergetter(1),reverse=True)
   return sortedClassCount[0][0]

def file2Matrix(filename):
    fr = open(filename)
    arrayOLines = fr.readlines()
    lineSize = len(arrayOLines)
    returnMat = zeros((lineSize, 3))
    classLabelVector = []
    index = 0
    for line in arrayOLines:
        line = line.strip()
        listFromLine = line.split('\t')
        returnMat[index,:] = listFromLine[0:3]
        classLabelVector.append(int(listFromLine[-1] ))
        index += 1
    return returnMat, classLabelVector

reason: from numpy


text:
import argparse
import cv2
import glob
import math
import numpy as np
import os
import re
import librosa
import tensorflow as tf
from collections import namedtuple
from datetime import datetime
from scipy.io import wavfile
from traceback import print_exc

Image = namedtuple('Image', 'rows cols depth data')
Audio = namedtuple('Audio', 'mics samples data')

_NUMBER_OF_SAMPLES = 1024
_FRAMES_PER_SECOND = 12

'''Create tfrecords for VGG sound dataset'''

def one_microphone(audio_data):
    """Python function to build a waveform audio from audio samples."""
    # choose index
    mic_id = 0
    # consider audio of one microphone
    audio_data_mic = audio_data[mic_id, :]
    return audio_data_mic
def _read_wav_audio_data(filename, video_time):
    print('{} - Reading {}'.format(datetime.now(), filename))
    fs, data = wavfile.read(filename)
    samples = fs*video_time
    data = data[0:samples]
    data = np.int32(data)
    return data

def _read_raw_audio_data(audio_data_sample):
    audio_serialized = audio_data_sample.tostring()

    return Audio(mics=1, samples=_NUMBER_OF_SAMPLES, data=audio_serialized)

def str2dir(dir_name):
    if not os.path.isdir(dir_name):
        raise argparse.ArgumentTypeError('{} is not a directory!'.format(dir_name))
    # elif os.access(dirname, os.R_OK):
    #     return argparse.ArgumentTypeError('{} is not a readable directory!'.format(dirname))
    else:
        return os.path.abspath(os.path.expanduser(dir_name))

def _aspect_preserving_resize(image, smallest_side):
    """Resize images preserving the original aspect ratio.
    Args:
      image: A 3-D image.
      smallest_side: A python integer or scalar indicating the size of
        the smallest side after resize.
    Returns:
      resized_image: A 3-D resized image.
    """
    shape = np.shape(image)
    height = shape[0]
    width = shape[1]
    new_height, new_width = _smallest_size_at_least(height, width, smallest_side)
    resized_image = cv2.resize(image, (new_width, new_height), interpolation=cv2.INTER_CUBIC)
    return resized_image

def _smallest_size_at_least(height, width, smallest_side):
    """Computes new shape with the smallest side equal to `smallest_side`.
    Computes new shape with the smallest side equal to `smallest_side` while
    preserving the original aspect ratio.
    Args:
      height: an int32 scalar indicating the current height.
      width: an int32 scalar indicating the current width.
      smallest_side: A python integer or scalar indicating the size of
        the smallest side after resize.
    Returns:
      new_height: an int32 scalar indicating the new height.
      new_width: and int32 scalar indicating the new width.
    """
    height = float(height)
    width = float(width)
    smallest_side = float(smallest_side)

    if height > width:
        scale = smallest_side / width
    else:
        scale = smallest_side / height
    new_height = int(height * scale)
    new_width = int(width * scale)
    return new_height, new_width

def _crop(image, offset_height, offset_width, crop_height, crop_width):
    """Crops the given image using the provided offsets and sizes.
    Note that the method doesn't assume we know the input image size but it does
    assume we know the input image rank.
    Args:
      image: an image of shape [height, width, channels].
      offset_height: a scalar indicating the height offset.
      offset_width: a scalar indicating the width offset.
      crop_height: the height of the cropped image.
      crop_width: the width of the cropped image.
    Returns:
      the cropped (and resized) image.
    Raises:
      InvalidArgumentError: if the rank is not 3 or if the image dimensions are
        less than the crop size.
    """
    original_shape = np.shape(image)

    # Use tf.slice instead of crop_to_bounding box as it accepts tensors to
    # define the crop size.
    image = image[offset_height:crop_height + offset_height, offset_width:crop_width + offset_width, :]
    return image

def _central_crop(image, crop_height, crop_width):
    """Performs central crops of the given image.
    Args:
      image: image.
      crop_height: the height of the image following the crop.
      crop_width: the width of the image following the crop.
    Returns:
      cropped image.
    """
    outputs = []
    image_height = np.shape(image)[0]
    image_width = np.shape(image)[1]

    offset_height = (image_height - crop_height) // 2
    offset_width = (image_width - crop_width) // 2

    image = _crop(image, offset_height, offset_width,
                  crop_height, crop_width)
    return image

def _read_video_frame(filename):
    print('{} - Reading {}'.format(datetime.now(), filename))

    image_raw = cv2.imread(filename)

    # rows = image_raw.shape[0]
    # cols = image_raw.shape[1]
    # depth = image_raw.shape[2]
    # image rescaled to give in input image aligned with acoustic image
    image = _aspect_preserving_resize(image_raw, 224)
    image = _central_crop(image, 224, 298)
    rows = image.shape[0]
    cols = image.shape[1]
    depth = image.shape[2]
    # don't take crop
    # image = _central_crop(image, _IMAGE_SIZE, _IMAGE_SIZE)
    image_serialized = image.tostring()
    return Image(rows=rows, cols=cols, depth=depth, data=image_serialized)

def _int64_feature(value):
    return tf.train.Feature(int64_list=tf.train.Int64List(value=[value]))

def _bytes_feature(value):
    return tf.train.Feature(bytes_list=tf.train.BytesList(value=[value]))

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('root_raw_dir', help='Synchronized raw files data set root directory', type=str2dir)
    parser.add_argument('out_dir', help='Directory where to store the converted data', type=str2dir)
    parser.add_argument('--modalities', help='Modalities to consider. 0: Audio images. 1: Audio data. 2: Video data.',
                        nargs='*', type=int)
    parsed_args = parser.parse_args()
    root_raw_dir = parsed_args.root_raw_dir
    out_dir = parsed_args.out_dir
    modalities = parsed_args.modalities
    include_audio_data = modalities is None or 1 in modalities
    include_video_data = modalities is None or 2 in modalities

    data_dirs = sorted(glob.glob('{}/*/*/video/'.format(root_raw_dir)))

    for data_mat_dir in data_dirs:
        # / media / vsanguineti / D90A - E56D / dualcam_actions_dataset / sync / Location_1 / 02
        # _Nuno / data_001
        # / media / vsanguineti / D90A - E56D / dualcam_actions_dataset / sync / classes/ data_001
        splitted_data_dir = data_mat_dir.split('/')
        data_dir_file = str.join('/', splitted_data_dir[:-2])
        video_time_filename = 'video_time.txt'
        try:
            fd = open(data_dir_file + '/' + video_time_filename, 'r')
            info_time = fd.readline().split(':')[1]
            fd.close()
            video_time = int(info_time.strip())

        except:
            print('error during reading video_time.txt')
            print_exc()
        if video_time>0 :
            classes = int(next(filter(re.compile('class_.*').match, splitted_data_dir)).split('_')[1])
            # subject = filter(re.compile('0[1-9]_.*').match, splitted_data_dir)[0]
            # subject_idx = int(subject.split('_')[0])
            location = int(next(filter(re.compile('data_.*').match, splitted_data_dir)).split('_')[1])

            # data_raw_audio_dir = data_mat_dir.replace(root_dir, root_raw_dir).replace('Multispectral_Acoustic_Image', 'audio')
            # data_raw_video_dir = data_mat_dir.replace(root_dir, root_raw_dir).replace('Multispectral_Acoustic_Image', 'video')
            data_raw_audio_dir = data_mat_dir.replace(root_raw_dir, root_raw_dir).replace('video', 'audio')
            data_raw_video_dir = data_mat_dir

            num_raw_video_files = len([name for name in os.listdir(data_raw_video_dir) if name.endswith('.bmp')])

            # compute number of videos (number ot total video files)

            frames_per_video = _FRAMES_PER_SECOND * video_time

            num_samples = _FRAMES_PER_SECOND
            filename = '{}/output_audio2.wav'.format(data_raw_audio_dir)
            audio_data_tot = _read_wav_audio_data(filename,  video_time)
            # changed
            for idx in range(video_time):
                start_index = idx * num_samples

                if include_audio_data:

                    audio_data_raw = [audio_data_tot[index*_NUMBER_OF_SAMPLES:(index+1)*_NUMBER_OF_SAMPLES] for
                        index in range(start_index, start_index + num_samples)]
                    audio_data = [_read_raw_audio_data(audio_d) for audio_d in audio_data_raw]
                else:
                    audio_data = None

                if include_video_data:
                    raw_video_files = ['{}/I_{:06d}.bmp'.format(data_raw_video_dir, index + 1) for
                                       index in range(start_index, start_index + num_samples)]
                    video_images = [_read_video_frame(filename) for filename in raw_video_files]
                else:
                    video_images = None

                out_data_dir = '{}/class_{}/data_{:0>3d}/'.format(out_dir, classes, location)
                out_filename = '{}/Data_{:0>3d}.tfrecord'.format(out_data_dir, idx + 1)

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

                print('{} - Writing {}'.format(datetime.now(), out_filename))

                with tf.python_io.TFRecordWriter(out_filename, options=tf.python_io.TFRecordOptions(
                        compression_type=tf.python_io.TFRecordCompressionType.GZIP)) as writer:
                    # Store audio and video data properties as context features, assuming all sequences are the same size
                    feature = {
                        'classes': _int64_feature(classes),
                        'location': _int64_feature(location)
                    }
                    if include_audio_data:
                        feature.update({
                            'audio_data/mics': _int64_feature(audio_data[0].mics),
                            'audio_data/samples': _int64_feature(audio_data[0].samples)
                        })
                    if include_video_data:
                        feature.update({
                            'video/height': _int64_feature(video_images[0].rows),
                            'video/width': _int64_feature(video_images[0].cols),
                            'video/depth': _int64_feature(video_images[0].depth),
                        })
                    feature_list = {}
                    if include_audio_data:
                        feature_list.update({
                            'audio/data': tf.train.FeatureList(
                                feature=[_bytes_feature(audio_sample.data) for audio_sample in audio_data])
                        })
                    if include_video_data:
                        feature_list.update({
                            'video/image': tf.train.FeatureList(
                                feature=[_bytes_feature(video_image.data) for video_image in video_images])
                        })
                    context = tf.train.Features(feature=feature)
                    feature_lists = tf.train.FeatureLists(feature_list=feature_list)
                    sequence_example = tf.train.SequenceExample(context=context, feature_lists=feature_lists)
                    writer.write(sequence_example.SerializeToString())

max_stars_repo_name: IIT-PAVIS/Acoustic-Image-Generation

max_issues_repo_name: IIT-PAVIS/Acoustic-Image-Generation

max_forks_repo_name: IIT-PAVIS/Acoustic-Image-Generation

reason: import numpy,from scipy


text:
#!/usr/bin/python

import sys
import numpy as np

start_time = time.time()

marker = 1
node = {}
local_skyline = {}
candidate_skyline = []
global_skyline = []
shadow_skyline = {}
virtual_point = {}
why_not_points = {}
n_updated_flag = {}
data_length = 0
t = 5
rsl = []
safe_region = []
number_of_preference = 0
customer_skyline = {}
customer_index = 0
query_point = []
list_customer = []
ct = []
ct_has_skyline = True
product_list = "TESTING_FC_D3_N100.txt"
user_preference = "TESTING_USER_D3_N10.txt"
intersection = []
ct_cost = []
q_cost = []
jumlah_rsl = 0
list_rsl = []

def insert_local_skyline(current_specs, bitmap):
	global local_skyline
	global shadow_skyline
	global virtual_point

	for i in range(0, len(local_skyline[bitmap])):
		dominating_local = False
		dominated_by_local = False
		for j in range(1, len(current_specs)):
			if(current_specs[j] != 'null' and local_skyline[bitmap][i][j] != 'null'):
				if(current_specs[j] < local_skyline[bitmap][i][j]):
					dominating_local = True
				elif(current_specs[j] > local_skyline[bitmap][i][j]):
					dominated_by_local = True
		if(dominating_local == True and dominated_by_local == False):
			local_skyline[bitmap][i][-1] = 'delete'
		elif(dominating_local == False and dominated_by_local == True):
			for k in range(0, i+1):
				local_skyline[bitmap][k][-1] = 'ok'
			return False
	dominated = 0
	for i in range(0, len(virtual_point[bitmap])):
		dominating_virtual = False
		dominated_by_virtual = False
		for j in range(1, len(current_specs)):
			if(current_specs[j] != 'null' and virtual_point[bitmap][i][j] != 'null'):
				if(current_specs[j] < virtual_point[bitmap][i][j]):
					dominating_virtual = True
				elif(current_specs[j] > virtual_point[bitmap][i][j]):
					dominated_by_virtual = True
		if(dominating_virtual == False and dominated_by_virtual == True):
			dominated = 1
			break
	if(dominated == 0):
		content = list(current_specs)
		content.append('ok')
		local_skyline[bitmap].append(content)
		for i in sorted(local_skyline[bitmap], reverse=True):
			if (i[-1] == 'delete'):
				local_skyline[bitmap].remove(i)
		for i in range(0, len(shadow_skyline[bitmap])):
			dominating_shadow = False
			dominated_by_shadow = False
			for j in range(1, len(current_specs)):
				if(current_specs[j] != 'null' and shadow_skyline[bitmap][i][j] != 'null'):
					if(current_specs[j] < shadow_skyline[bitmap][i][j]):
						dominating_shadow = True
					elif(current_specs[j] > shadow_skyline[bitmap][i][j]):
						dominated_by_shadow = True
			if(dominating_shadow == True and dominated_by_shadow == False):
				shadow_skyline[bitmap][i][-1] = 'delete'
		for i in sorted(shadow_skyline[bitmap], reverse=True):
			if (i[-1] == 'delete'):
				shadow_skyline[bitmap].remove(i)
		return True
	elif(dominated == 1):
		n_updated_flag[bitmap] = True
		for i in range(0, len(shadow_skyline[bitmap])):
			dominating_shadow = False
			dominated_by_shadow = False
			for j in range(1, len(current_specs)):
				if(current_specs[j] != 'null' and shadow_skyline[bitmap][i][j] != 'null'):
					if(current_specs[j] < shadow_skyline[bitmap][i][j]):
						dominating_shadow = True
					elif(current_specs[j] > shadow_skyline[bitmap][i][j]):
						dominated_by_shadow = True
			if(dominating_shadow == True and dominated_by_shadow == False):
				shadow_skyline[bitmap][i][-1] = 'delete'
		for i in sorted(shadow_skyline[bitmap], reverse=True):
			if (i[-1] == 'delete'):
				shadow_skyline[bitmap].remove(i)
		content = list(current_specs)
		content.append('ok');
		shadow_skyline[bitmap].append(content)
	return False

def insert_candidate_skyline(current_specs, bitmap):
	global candidate_skyline
	list_bit_inserted = []
	dominated = 0

	for i in range(0, len(candidate_skyline)):
		dominating_candidate = False
		dominated_by_candidate = False
		for j in range(1, len(current_specs)):
			if(current_specs[j] != 'null' and candidate_skyline[i][j] != 'null'):
				if(current_specs[j] < candidate_skyline[i][j]):
					dominating_candidate = True
				elif(current_specs[j] > candidate_skyline[i][j]):
					dominated_by_candidate = True
		if(dominating_candidate == True and dominated_by_candidate == False):
			candidate_skyline[i][-1] = 'delete'
			if(candidate_skyline[i][-2] not in list_bit_inserted):
				insert_virtual_point(current_specs, candidate_skyline[i][-2])
				list_bit_inserted.append(candidate_skyline[i][-2])
		elif(dominating_candidate == False and dominated_by_candidate == True):
			content = list(candidate_skyline[i][:-2])
			insert_virtual_point(content, bitmap)
			dominated = 1
	candidate_skyline = [i for i in candidate_skyline if i[-1] == 'ok']
	if(dominated == 0):
		content = list(current_specs)
		content.append(bitmap)
		content.append('ok')
		candidate_skyline.append(content)

def insert_virtual_point(current_specs, bitmap):
	global local_skyline
	global virtual_point
	global shadow_skyline
	#MEMBANDINGKAN DENGAN LOCAL SKYLINE
	#Move all dominated local_skyline N to shadow_skyline
	for i in range(0, len(local_skyline[bitmap])):
		dominating_local = False
		dominated_by_local = False
		for j in range(1, len(current_specs)):
			if(current_specs[j] != 'null' and local_skyline[bitmap][i][j] != 'null'):
				if(current_specs[j] < local_skyline[bitmap][i][j]):
					dominating_local = True
				elif(current_specs[j] > local_skyline[bitmap][i][j]):
					dominated_by_local = True
		if(dominating_local == True and dominated_by_local == False):
			local_skyline[bitmap][i][-1] = 'delete'
	for i in reversed(local_skyline[bitmap]):
		if (i[-1] == 'delete'):
			shadow_skyline[bitmap].append(i)
			local_skyline[bitmap].remove(i)
			shadow_skyline[bitmap][-1][-1] = 'ok'

	for i in range(0, len(virtual_point[bitmap])):
		dominating_virtual = False
		dominated_by_virtual = False
		superset_check = 0

		for j in range(1, len(current_specs)):
			if(current_specs[j] != 'null' and virtual_point[bitmap][i][j] != 'null'):
				if(current_specs[j] < virtual_point[bitmap][i][j]):
					dominating_virtual = True
				elif(current_specs[j] > virtual_point[bitmap][i][j]):
					dominated_by_virtual = True
			if(current_specs[j] != 'null'):
				superset_check += 1
			elif(current_specs[j] == 'null' and virtual_point[bitmap][i][j] == 'null'):
				superset_check += 1
		if(dominating_virtual == True and dominated_by_virtual == False and superset_check == len(current_specs)):
			virtual_point[bitmap][i][-1] = 'delete'

	for i in reversed(virtual_point[bitmap]):
		if(i[-1] == 'delete'):
			virtual_point[bitmap].remove(i)
	content = list(current_specs)
	content.append('ok')
	virtual_point[bitmap].append(content)

def update_global_skyline():
	global global_skyline
	global candidate_skyline
	global shadow_skyline
	global data_length
	for c in range(0, len(candidate_skyline)):
		for g in range(0, len(global_skyline)):
			dominating_global = False
			dominating_candidate = False
			for i in range(1, data_length):
				if(candidate_skyline[c][i] != 'null' and global_skyline[g][i] != 'null'):
					if(candidate_skyline[c][i] < global_skyline[g][i]):
						dominating_global = True
					elif(candidate_skyline[c][i] > global_skyline[g][i]):
						dominating_candidate = True
			if(dominating_global == True and dominating_candidate == False):
				global_skyline[g][-1] = 'delete'
			elif(dominating_global == False and dominating_candidate == True):
				candidate_skyline[c][-1] = 'delete'

	for i in reversed(global_skyline):
		if(i[-1] == 'delete'):
			global_skyline.remove(i)
	for i in reversed(candidate_skyline):
		if(i[-1] == 'delete'):
			candidate_skyline.remove(i)
	for g in range(0, len(global_skyline)):
		for i in n_updated_flag:
			if (n_updated_flag[i] == True and i in shadow_skyline):
				for s in range(0, len(shadow_skyline[i])):
					dominating_global = False
					dominating_shadow = False
					for j in range(1, data_length):
						if(global_skyline[g][j] != 'null' and shadow_skyline[i][s][j] != 'null' ):
							if(global_skyline[g][j] < shadow_skyline[i][s][j]):
								dominating_shadow = True
							elif(global_skyline[g][j] > shadow_skyline[i][s][j]):
								dominating_global = True
					if(dominating_global == True and dominating_shadow == False):
						global_skyline[g][-1] == 'delete'

	for c in range(0, len(candidate_skyline)):
		for i in node:
			for s in range(0, len(shadow_skyline[i])):
				dominating_candidate = False
				dominating_shadow = False
				for j in range(1, data_length):
					if(candidate_skyline[c][j] != 'null' and shadow_skyline[i][s][j] != 'null'):
						if(candidate_skyline[c][j] < shadow_skyline[i][s][j]):
							dominating_shadow = True
						elif(candidate_skyline[c][j] > shadow_skyline[i][s][j]):
							dominating_candidate = True
				if(dominating_candidate == True and dominating_shadow == False):
					candidate_skyline[c][-1] = 'delete'

	for i in reversed(global_skyline):
		if(i[-1] == 'delete'):
			global_skyline.remove(i)
	for i in reversed(candidate_skyline):
		if(i[-1] == 'delete'):
			candidate_skyline.remove(i)
	for i in candidate_skyline:
		global_skyline.append(i)

	for i in n_updated_flag:
		n_updated_flag[i] == False

def calculate_rsl_q(customer_skyline, query_point):
	global data_length
	for dict_index in customer_skyline:
		transformed_query_point = []
		for q in range(0, data_length):
			transformed_value = abs(float(query_point[q+1]) - customer_skyline[dict_index][-2][q])
			transformed_query_point.append(transformed_value)
		for data_index in range(0, len(customer_skyline[dict_index]) - 2):
			dominating_q = False
			dominating_customer = False
			for i in range(0, data_length):
				if(customer_skyline[dict_index][data_index][i+1] != 'null'):
					if(customer_skyline[dict_index][data_index][i+1] < transformed_query_point[i]):
						dominating_q = True
					elif(customer_skyline[dict_index][data_index][i+1] > transformed_query_point[i]):
						dominating_customer = True
			if(dominating_q == True and dominating_customer == False):
				customer_skyline[dict_index][-1] = 'delete'
			elif(dominating_q == False and dominating_customer == True):
				customer_skyline[dict_index][data_index][-1] = 'delete'
		if(customer_skyline[dict_index][-1] == 'ok'):
			for i in range(len(customer_skyline[dict_index]) - 3, -1, -1):
				if(customer_skyline[dict_index][i][-1] == 'delete'):
					customer_skyline[dict_index].remove(customer_skyline[dict_index][i])
		if(len(customer_skyline[dict_index]) <= 2):
			customer_skyline[dict_index][-1] = 'delete'
	return customer_skyline

def Prepare_Data(line, customer):
	global bitmap
	global node
	global local_skyline
	global candidate_skyline
	global global_skyline
	global shadow_skyline
	global virtual_point
	global n_updated_flag
	global data_length

	current_spec = line.split()
	data = []
	transformed_data = []
	transformed_data.append(current_spec[0])
	data.append(current_spec[0])
	for i in range(1, data_length+1):
		if(current_spec[i] == "null"):
			bitmap += "0"
			data.append(current_spec[i])
			transformed_data.append(current_spec[i])
		else:
			bitmap += "1"
			difference = abs(float(current_spec[i]) - customer[i-1])
			data.append(float(current_spec[i]))
			transformed_data.append(float(difference))
	if bitmap not in node:
		node[bitmap] = []
		local_skyline[bitmap] = []
		shadow_skyline[bitmap] = []
		virtual_point[bitmap] = []
		n_updated_flag[bitmap] = False
	return transformed_data

data_length = len(ct)
fu = open(user_preference)
for list_user in fu:
	temp = [float(x) for x in list_user.split()]
	list_customer.append(temp)

for x in range(0, len(list_customer)):
	fp = open(product_list)
	node.clear()
	local_skyline.clear()
	candidate_skyline.clear()
	global_skyline.clear()
	shadow_skyline.clear()
	virtual_point.clear()
	number_of_preference += 1
	for line in fp:
		bitmap = ""
		transformed_data = Prepare_Data(line, list_customer[x])
		is_skyline = insert_local_skyline(transformed_data, bitmap)
		if is_skyline == True:
			insert_candidate_skyline(transformed_data, bitmap)
			if(len(candidate_skyline) > t):
				update_global_skyline()
				candidate_skyline.clear()
	fp.close()
	update_global_skyline()
	if(len(global_skyline) > 0):
		customer_skyline[str(customer_index)] = list(global_skyline)
		customer_skyline[str(customer_index)].append(list(list_customer[x]))
		customer_skyline[str(customer_index)].append("ok")
		customer_index += 1
fu.close()
ddr_prime_ct = generate_ddr_prime_ct(ct)
safe_region = generate_safe_region_q()

if(ct_has_skyline == True):
	intersection_status = check_intersection(safe_region, ddr_prime_ct)
else:
	intersection_status = False
if(intersection_status == True):
	move_query_point()
else:
	Move_Why_Not_And_Query_Point()