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

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

#

import numpy as np

from kernel.base.instance import Instance
from kernel.base.statics import MultivariateStatisticalSummary
from kernel.utils import consts

class LeastSquaredErrorLoss(object):
    @staticmethod
    def initialize(y):
        y_inst = y.mapValues(lambda label: Instance(features=np.asarray([label])))
        y_inst.schema = {"header": ["label"]}
        statistics = MultivariateStatisticalSummary(y_inst, -1)
        mean = statistics.get_mean()["label"]
        return y.mapValues(lambda x: np.asarray([mean])), np.asarray([mean])

    @staticmethod
    def predict(value):
        return value

    @staticmethod
    def compute_loss(y, y_pred):
        lse_loss = y.join(y_pred, lambda y, yp: ((y - yp) * (y - yp), 1))
        lse_sum, sample_num = lse_loss.reduce(lambda tuple1, tuple2: (tuple1[0] + tuple2[0], tuple1[1] + tuple2[1]))
        return lse_sum / sample_num

    @staticmethod
    def compute_grad(y, y_pred):
        return 2 * (y_pred - y)

    @staticmethod
    def compute_hess(y, y_pred):
        if type(y).__name__ == "ndarray" or type(y_pred).__name__ == "ndarray":
            shape = (y - y_pred).shape
            return np.full(shape, 2)
        else:
            return 2

class LeastAbsoluteErrorLoss(object):
    @staticmethod
    def initialize(y):
        y_inst = y.mapValues(lambda label: Instance(features=np.asarray([label])))
        y_inst.schema = {"header": ["label"]}
        statistics = MultivariateStatisticalSummary(y_inst, -1)
        median = statistics.get_median()["label"]
        return y.mapValues(lambda x: np.asarray([median])), np.asarray([median])

    @staticmethod
    def predict(value):
        return value

    @staticmethod
    def compute_loss(y, y_pred):
        lae = y.join(y_pred, lambda y, yp: (np.abs(y - yp), 1))
        lae_sum, sample_num = lae.reduce(lambda tuple1, tuple2: (tuple1[0] + tuple2[0], tuple1[1] + tuple2[1]))
        return lae_sum / sample_num

    @staticmethod
    def compute_grad(y, y_pred):
        if type(y).__name__ == "ndarray" or type(y_pred).__name__ == "ndarray":
            diff = y_pred - y
            diff[diff > consts.FLOAT_ZERO] = 1
            diff[diff < consts.FLOAT_ZERO] = -1
            diff[np.abs(diff) <= consts.FLOAT_ZERO] = 0

            return diff
        else:
            diff = y_pred - y
            if diff > consts.FLOAT_ZERO:
                return 1
            elif diff < consts.FLOAT_ZERO:
                return -1
            else:
                return 0

    @staticmethod
    def compute_hess(y, y_pred):
        if type(y).__name__ == "ndarray" or type(y_pred).__name__ == "ndarray":
            shape = (y - y_pred).shape
            return np.full(shape, 1)
        else:
            return 1

class HuberLoss(object):
    @staticmethod
    def initialize(y):
        y_inst = y.mapValues(lambda label: Instance(features=np.asarray([label])))
        y_inst.schema = {"header": ["label"]}
        statistics = MultivariateStatisticalSummary(y_inst, -1)
        mean = statistics.get_mean()["label"]
        return y.mapValues(lambda x: np.asarray([mean])), np.asarray([mean])

    def __init__(self, delta):
        if delta is None:
            self.delta = consts.FLOAT_ZERO
        else:
            self.delta = delta

        if np.abs(self.delta) < consts.FLOAT_ZERO:
            self.delta = consts.FLOAT_ZERO

    def compute_loss(self, y, y_pred):
        huber_loss = y.join(y_pred, lambda y, yp:
        (self.delta ** 2 * (np.sqrt(1 + ((yp - y) / self.delta) ** 2) - 1), 1))
        huber_sum, sample_num = huber_loss.reduce(lambda tuple1, tuple2: (tuple1[0] + tuple2[0], tuple1[1] + tuple2[1]))
        return huber_sum / sample_num

    @staticmethod
    def predict(value):
        return value

    def compute_grad(self, y, y_pred):
        diff = y_pred - y
        return diff / np.sqrt(1.0 + diff * diff / (self.delta ** 2))

    def compute_hess(self, y, y_pred):
        diff = y_pred - y
        return 1.0 / (1.0 + diff * diff / (self.delta ** 2)) ** 1.5

class FairLoss(object):
    @staticmethod
    def initialize(y):
        y_inst = y.mapValues(lambda label: Instance(features=np.asarray([label])))
        y_inst.schema = {"header": ["label"]}
        statistics = MultivariateStatisticalSummary(y_inst, -1)
        mean = statistics.get_mean()["label"]
        return y.mapValues(lambda x: np.asarray([mean])), np.asarray([mean])

    def __init__(self, c):
        if c is None:
            self.c = consts.FLOAT_ZERO
        else:
            self.c = c

        if np.abs(self.c) < consts.FLOAT_ZERO:
            self.c = consts.FLOAT_ZERO

    @staticmethod
    def predict(value):
        return value

    def compute_loss(self, y, y_pred):
        fair_loss = y.join(y_pred, lambda y, yp:
        (self.c * np.abs(yp - y) - self.c ** 2 * np.log(np.abs(yp - y) / self.c + 1), 1))
        fair_loss_sum, sample_num = fair_loss.reduce(
            lambda tuple1, tuple2: (tuple1[0] + tuple2[0], tuple1[1] + tuple2[1]))
        return fair_loss_sum / sample_num

    def compute_grad(self, y, y_pred):
        diff = y_pred - y
        return self.c * diff / (np.abs(diff) + self.c)

    def compute_hess(self, y, y_pred):
        diff = y_pred - y
        return self.c ** 2 / (np.abs(diff) + self.c) ** 2

class LogCoshLoss(object):
    @staticmethod
    def initialize(y):
        y_inst = y.mapValues(lambda label: Instance(features=np.asarray([label])))
        y_inst.schema = {"header": ["label"]}
        statistics = MultivariateStatisticalSummary(y_inst, -1)
        mean = statistics.get_mean()
        return y.mapValues(lambda x: np.asarray([mean])), np.asarray([mean])

    @staticmethod
    def predict(value):
        return value

    @staticmethod
    def compute_loss(y, y_pred):
        log_cosh_loss = y.join(y_pred, lambda y, yp: (np.log(np.cosh(yp - y)), 1))
        log_cosh_sum, sample_num = log_cosh_loss.reduce(
            lambda tuple1, tuple2: (tuple1[0] + tuple2[0], tuple1[1] + tuple2[1]))
        return log_cosh_sum / sample_num

    @staticmethod
    def compute_grad(y, y_pred):
        return np.tanh(y_pred - y)

    @staticmethod
    def compute_hess(y, y_pred):
        return 1 - np.tanh(y_pred - y) ** 2

class TweedieLoss(object):
    @staticmethod
    def initialize(y):
        y_inst = y.mapValues(lambda label: Instance(features=np.asarray([label])))
        y_inst.schema = {"header": ["label"]}
        statistics = MultivariateStatisticalSummary(y_inst, -1)
        mean = statistics.get_mean()["label"]
        return y.mapValues(lambda x: np.asarray([mean])), np.asarray([mean])

    def __init__(self, rho=None):
        if rho is None:
            self.rho = consts.FLOAT_ZERO
        else:
            self.rho = rho

    @staticmethod
    def predict(value):
        return value

    def compute_loss(self, y, y_pred):
        tweedie_loss = y.join(y_pred,
                              lambda y, yp:
                              (-y * np.exp(1 - self.rho) * np.log(max(yp, consts.FLOAT_ZERO)) / (1 - self.rho) +
                               np.exp(2 - self.rho) * np.log(max(consts.FLOAT_ZERO, yp)) / (2 - self.rho), 1))
        tweedie_loss_sum, sample_num = tweedie_loss.reduce(lambda tuple1, tuple2:
                                                           (tuple1[0] + tuple2[0], tuple1[1] + tuple2[1]))
        return tweedie_loss_sum / sample_num

    def compute_grad(self, y, y_pred):
        return -y * np.exp(1 - self.rho) * y_pred + np.exp(2 - self.rho) * y_pred

    def compute_hess(self, y, y_pred):
        return -y * (1 - self.rho) * np.exp(1 - self.rho) * y_pred + (2 - self.rho) * np.exp(2 - self.rho) * y_pred

max_stars_repo_name: rinceyuan/WeFe

max_issues_repo_name: rinceyuan/WeFe

max_forks_repo_name: rinceyuan/WeFe


text:
import ctypes

import numpy as np
import jpype
import scyjava

from imglyb import util

def _java_setup():
    """
    Lazy initialization function for Java-dependent data structures.
    Do not call this directly; use scyjava.start_jvm() instead.
    """
    Intervals = scyjava.jimport('net.imglib2.util.Intervals')

scyjava.when_jvm_starts(_java_setup)

dtype_selector = {
    'FloatType': np.dtype('float32'),
    'DoubleType': np.dtype('float64'),
    'ByteType': np.dtype('int8'),
    'UnsignedByteType': np.dtype('uint8'),
    'ShortType': np.dtype('int16'),
    'UnsignedShortType': np.dtype('uint16'),
    'IntType': np.dtype('int32'),
    'UnsignedIntType': np.dtype('uint32'),
    'LongType': np.dtype('int64'),
    'UnsignedLongType': np.dtype('uint64'),
}

ctype_conversions_imglib = {
    'FloatType': ctypes.c_float,
    'DoubleType': ctypes.c_double,
    'ByteType': ctypes.c_int8,
    'UnsignedByteType': ctypes.c_uint8,
    'ShortType': ctypes.c_int16,
    'UnsignedShortType': ctypes.c_uint16,
    'IntType': ctypes.c_int32,
    'UnsignedIntType': ctypes.c_uint32,
    'LongType': ctypes.c_int64,
    'UnsignedLongType': ctypes.c_uint64
}

def get_address(rai):
    scyjava.start_jvm()
    class_name = util.Helpers.classNameSimple(rai)
    class_name_full = util.Helpers.className(rai)
    if class_name in ('ArrayImg', 'UnsafeImg'):
        access = jpype.JObject(class_name_full, rai).update(None)
        access_type = util.Helpers.className(access)

        if 'basictypelongaccess.unsafe' in access_type:
            return jpype.JObject(access_type, access).getAddress()
        else:
            raise ValueError("Excpected unsafe access but got {}".format(access_type))

    else:
        raise ValueError("Excpected ArrayImg or UnsafeImg but got {}".format(class_name))

class ImgLibReferenceGuard(np.ndarray):
    def __new__(cls, rai):
        access = rai.randomAccess()
        rai.min(access)
        imglib_type = util.Helpers.classNameSimple(access.get())
        address = get_address(rai)

        shape = tuple(Intervals.dimensionsAsLongArray(rai))[::-1]
        dtype = dtype_selector[imglib_type]
        pointer = ctypes.cast(address, ctypes.POINTER(ctype_conversions_imglib[imglib_type]))
        order = 'C'

        obj = np.ndarray.__new__(cls, buffer=np.ctypeslib.as_array(pointer, shape=shape), shape=shape, dtype=dtype)
        obj.setflags(write=True)
        obj.rai = rai
        return obj

    def __array_finalize__(self, obj):
        if obj is None:
            return
        self.rai = obj.rai

if __name__ == "__main__":
    ArrayImgs = scyjava.jimport('net.imglib2.img.array.ArrayImgs')
    UnsafeUtil = scyjava.jimport('net.imglib2.img.basictypelongaccess.unsafe.UnsafeUtil')
    Arrays = scyjava.jimport('java.util.Arrays')
    OwningFloatUnsafe = scyjava.jimport('net.imglib2.img.basictypelongaccess.unsafe.owning.OwningFloatUnsafe')
    Fraction = scyjava.jimport('net.imglib2.util.Fraction')
    LongStream = scyjava.jimport('java.util.stream.LongStream')

    shape = (2, 3, 4)
    n_elements = int(np.prod(shape))
    data_store = OwningFloatUnsafe(n_elements)
    dim_array = LongStream.of(*shape).toArray()
    print(Arrays.toString(dim_array))
    rai = util.Helpers.toArrayImg(jpype.JObject(util.Helpers.className(data_store), data_store), dim_array)
    # rai = ArrayImgs.floats( *shape )
    c = rai.cursor()
    count = 23
    while c.hasNext():
        c.next().setReal(count)
        count += 1
    print(util.Helpers.className(rai.randomAccess().get()))
    print(util.Helpers.classNameSimple(rai.randomAccess().get()))
    arr = ImgLibReferenceGuard(rai)
    print(arr, arr.mean())
    c = rai.cursor()
    c.fwd()
    c.next().setReal(0)
    print(arr, arr.mean())

ext: py

max_stars_repo_name: hinerm/imglyb

max_issues_repo_name: hinerm/imglyb

max_forks_repo_name: hinerm/imglyb

reason: import numpy


text:
#!/usr/bin/env python
"""
lambdata - a collection of Data Science helper functions
"""
import pandas as pd
import numpy as np
from . import example_module

Y = example_module.increment(example_module.x)
TEST = pd.DataFrame(np.ones(10))

""" For Fun """
class Dog:
  """Testing Woof"""
    def __init__(self, name):
      self.name = name
      self.tricks = []

    def add_trick(self, trick):
      self.tricks.append(trick)

class DataCleaner:
  """A Class for cleaning up your data"""
    def __init__(self, df):
       self.df = df
    def remove_nulls(self):
       self.df


text:
import numpy as np

from ..utils.magic import MagicSquare

class Searcher():

    def __init__(self, magic_square, succ_num):
        self.magic_square = magic_square
        self.max_violation_number = (self.magic_square['n'] + 1) * 2
        self.iter = None
        self.all_violations = None
        self.succ_num = succ_num

    def find(self, iterations):
        pass

    def print_solution(self):
        if self.magic_square is None:
            return
        print('============================================================\n'
              'Search type: {}\n'
              'iterations: {}\n'
              'violation number: {} of {}\n'
              'solution: \n{}'.format(
                    self.type, self.iter, self['viol num'],
                    self.max_violation_number, self['sol']
        ))

    def _init_start_state(self, iterations):
        self.magic_square.init_random()
        self.iter = iterations
        self.all_violations = [self['viol num']]

    def _should_break(self, curr_it):
        new_violation_num = self.magic_square.violation_number()
        self.all_violations.append(new_violation_num)
        if new_violation_num == 0:
            self.iter = curr_it + 1
            return True
        return False

    def __getitem__(self, key):
        if key == 'sol':
            return self.magic_square['matrix']
        if key == 'viol num':
            return self.magic_square.violation_number()
        if key == 'max viol num':
            return self.max_violation_number
        if key == 'type':
            return self.type
        if key == 'iter':
            return self.iter
        if key == 'viol. through iter.':
            return self.all_violations
        return None

max_stars_repo_name: IvanovicM/ai

max_issues_repo_name: IvanovicM/ai

max_forks_repo_name: IvanovicM/ai


text:
# This script uses MDAnalysis to count the number of C-alpha atoms at the interface of two proteins A and B.
# author @dmachalz
import MDAnalysis
import os
import numpy as np

# Variables
# The working directory containing the PDB files. Proteins A and B should be stored in one PDB file.
foldername='REPLACEME'
# Outputfile in CSV format
outcsv='REPLACEME.csv'
# Project required renaming
replacer='1jma_prep_corr_orient_'

outlist=[]
# Change to working directory
os.chdir(foldername)
pdblist = [p for p in os.listdir(foldername) if p.endswith('.pdb')]
for pdbfile in pdblist:
    # Getting the System props - highly project specific
    sysname = pdbfile.replace('.pdb','')
    sysnumber = sysname.replace(replacer, '')
    # Create MDAnalysis universe
    u = MDAnalysis.Universe(pdbfile)
    # Get atom selection
    asel = u.select_atoms("((name CA and segid A) and around 4.5 segid B) or ((name CA and segid B) and around 4.5 segid A)")
    # Get number of Calpha atoms (acount) in selection
    acount = len(asel)
    # Write pdbname, and acount to csv file
    #print(sysname, acount)
    outlist.append([int(sysnumber), acount])

# Converting output list to numpy array for easier writing
outlist = np.array(outlist)
# Actually writing out the list using 'integer' formatting and delimiter ',' to OUTCSV!
np.savetxt(outcsv, outlist, fmt='%i', delimiter=',')

max_stars_repo_name: dmachalz/mdanalysis

max_issues_repo_name: dmachalz/mdanalysis

max_forks_repo_name: dmachalz/mdanalysis


text:
import tensorflow as tf
import numpy as np

def preprocess_image(image_path, input_size):
  """Preprocess the input image to feed to the TFLite model"""
  img = tf.io.read_file(image_path)
  img = tf.io.decode_image(img, channels=3)
  img = tf.image.convert_image_dtype(img, tf.uint8)
  original_image = img
  resized_img = tf.image.resize(img, input_size)
  resized_img = resized_img[tf.newaxis, :]
  resized_img = tf.cast(resized_img, dtype=tf.uint8)
  return resized_img, original_image

def detect_objects(interpreter, image, threshold):
  """Returns a list of detection results, each a dictionary of object info."""
  # Get the Model Runner
  signature_fn = interpreter.get_signature_runner()
  # Feed the input image to the model and evaluate
  output = signature_fn(images=image)
  # Get all outputs from the model
  count = int(np.squeeze(output['output_0']))
  scores = np.squeeze(output['output_1'])
  classes = np.squeeze(output['output_2'])
  boxes = np.squeeze(output['output_3'])

  results = []
  for i in range(count):
    if scores[i] >= threshold:
      result = {
        'bounding_box': boxes[i],
        'class_id': classes[i],
        'score': scores[i]
      }
      results.append(result)
  return results


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 14 19:14:40 2021

Falling body estimation. Example 14.2 in Dan Simon's book "Optimal State Estimation"

@author: halvorak
"""

import numpy as np
import scipy.stats
import scipy.integrate
import matplotlib.pyplot as plt
import pathlib
import os
import scipy.linalg
import matplotlib.patches as plt_patches
import matplotlib.transforms as plt_transforms
import copy
import pickle
# Did some modification to these packages
from myFilter import UKF
from myFilter import sigma_points as ukf_sp
# from myFilter import UKF_constrained

#Self-written modules
import sigma_points_classes as spc
import unscented_transformation as ut
import utils_falling_body as utils_fb

N_LHS_list = np.logspace(1,3,num = 20, dtype = int)
# N_LHS_list = np.array([10, 15, 20]) #for test
n_filters = 5
mean_cost = {"x1": np.zeros((n_filters, N_LHS_list.shape[0])),
             "x2": np.zeros((n_filters, N_LHS_list.shape[0])),
             "x3": np.zeros((n_filters, N_LHS_list.shape[0]))
             }
sigma_cost = {"x1": np.zeros((n_filters, N_LHS_list.shape[0])),
             "x2": np.zeros((n_filters, N_LHS_list.shape[0])),
             "x3": np.zeros((n_filters, N_LHS_list.shape[0]))
             }
for q in range(len(N_LHS_list)):
    #%% For running the sim N times
    N = 100 # number of times every simulation is repeated
    N_lhs_dist = N_LHS_list[q]
    N_mc_dist = N_LHS_list[q]
    dim_x = 3
    j_valappil = np.zeros((dim_x, N))
    j_valappil_lhs = np.zeros((dim_x, N))
    j_valappil_mc = np.zeros((dim_x, N))
    j_valappil_mcm = np.zeros((dim_x, N))
    j_valappil_qf = np.zeros((dim_x, N))
    Ni = 0
    rand_seed = 1234
    while Ni < N:
        try:
            np.random.seed(rand_seed) #to get reproducible results. rand_seed updated in every iteration
            #%% Import the distributions of the parameters for the fx equations (states)
            #Modes of the dists are used for the true system, and mean of the dists are the parameters for the UKF

            # utils_fb.uncertainty_venturi2()
            _, _, Q_nom, R_nom = utils_fb.get_literature_values()
            par_dist_fx, par_det_fx, fig_p_fx, ax_p_fx = utils_fb.get_param_ukf_case1(
                par_fx = True,
                std_dev_prct = 0.05,
                plot_dist = False)

            mode_dist = np.zeros(len(par_dist_fx)) #for the mode of the distributions
            par_true_fx = par_det_fx.copy()
            par_kf_fx = par_det_fx.copy()
            del par_det_fx
            for key, dist in par_dist_fx.items():
                #Need to find the mode numericlly
                mode = scipy.optimize.minimize(lambda theta: -dist.pdf(theta),
                                               dist.mean(), #use mean as theta0
                                               tol = 1e-10)
                par_true_fx[key] = mode.x #true system uses the mode values
                # par_true_fx[key] = dist.mean()+dist.std() #true system uses the mean - std_dev
                par_kf_fx[key] = dist.mean() #the ukf uses mean values reported in the literaure

                # print(f"{key}_true/{key}_KF: {par_true_fx[key]/par_kf_fx[key]}")
            # par_kf_fx = par_true_fx
            #%% Import the distributions of the parameters for the hx equations (measurements)
            #Modes of the dists are used for the true system, and mean of the dists are the parameters for the UKF

            par_dist_hx, par_det_hx, fig_p_hx, ax_p_hx = utils_fb.get_param_ukf_case1(
                par_fx = False, #return based on par_true_hx in stead
                std_dev_prct = 0.05,
                plot_dist = False)
            mode_dist = np.zeros(len(par_dist_hx)) #for the mode of the distributions
            par_true_hx = par_det_hx.copy()
            par_kf_hx = par_det_hx.copy()
            del par_det_hx
            for key, dist in par_dist_hx.items():
                #Need to find the mode numericlly
                mode = scipy.optimize.minimize(lambda theta: -dist.pdf(theta),
                                               dist.mean(), #use mean as theta0
                                               tol = 1e-10)
                par_true_hx[key] = mode.x #true system uses the mode values
                # par_true_hx[key] = dist.mean() - dist.std() #true system uses the mean - std_dev
                par_kf_hx[key] = dist.mean() #the ukf uses mean values reported in the literaure

                # print(f"{key}_true/{key}_KF: {par_true_hx[key]/par_kf_hx[key]}")
            # par_kf_hx = par_true_hx
            #Measurement equation also need y_repeatability
            y_rep = scipy.stats.norm(loc = 0.,
                                     scale = np.sqrt(R_nom[0])) #this is just 1D, so ok
            par_true_hx["y_rep"] = y_rep.rvs()
            par_kf_hx["y_rep"] = y_rep.mean()
            par_dist_hx["y_rep"] = y_rep

            #%% Define samples for LHS
            # N_lhs_dist = 10
            labels_fx = par_true_fx.copy()
            labels_fx.pop("g")
            par_lhs_fx, samples_lhs_fx, fig_lhs, ax_lhs = utils_fb.get_lhs_points(par_dist_fx,
                                                                          N_lhs_dist = N_lhs_dist,
                                                                          plot_mc_samples=False,
                                                                          labels = list(labels_fx.keys()))
            par_lhs_fx["g"] = np.ones(N_lhs_dist)*par_kf_fx["g"] #append with constant g
            par_lhs_hx, samples_lhs_hx, fig_lhs, ax_lhs = utils_fb.get_lhs_points(par_dist_hx,
                                                                                N_lhs_dist = N_lhs_dist,
                                                                                plot_mc_samples=False,
                                                                                labels = list(par_true_hx.keys())
                                                                                )
            #%% Define samples for MC, random sampling
            # N_mc_dist = int(1e2)

            par_mc_fx, fig_mc, ax_mc = utils_fb.get_mc_points(par_dist_fx,
                                                                N_mc_dist = N_mc_dist,
                                                                plot_mc_samples=False,
                                                                labels = labels_fx
                                                                 )
            par_mc_fx["g"] = np.ones(N_mc_dist)*par_kf_fx["g"] #append with constant g
            par_mc_hx, fig_mc, ax_mc = utils_fb.get_mc_points(par_dist_hx,
                                                                N_mc_dist = N_mc_dist,
                                                                plot_mc_samples=False,
                                                                labels = list(par_true_hx.keys())
                                                                 )
            #%% Define samples for MCm, random sampling
            N_mcm_dist = int(5e2)

            par_mcm_fx, fig_mcm, ax_mcm = utils_fb.get_mc_points(par_dist_fx,
                                                                N_mc_dist = N_mcm_dist,
                                                                plot_mc_samples=False,
                                                                labels = labels_fx
                                                                 )
            par_mcm_fx["g"] = np.ones(N_mcm_dist)*par_kf_fx["g"] #append with constant g
            par_mcm_hx, fig_mcm, ax_mcm = utils_fb.get_mc_points(par_dist_hx,
                                                                N_mc_dist = N_mcm_dist,
                                                                plot_mc_samples=False,
                                                                labels = list(par_true_hx.keys())
                                                                 )

            #%% Define dimensions and initialize arrays
            # sigmas, w = utils_fb.get_sigmapoints_and_weights(par_dist)

            x0 = utils_fb.get_x0_literature()
            P0 = np.diag([3e8,#[ft^2], altitute, initial covariance matrix for UKF
                        4e6, # [ft^2], horizontal range
                        1e-6 # [?] ballistic coefficient
                        ])
            # print(x0+np.sqrt(np.diag(P0)))

            x0_kf = copy.deepcopy(x0) + np.sqrt(np.diag(P0)) #+ the standard deviation
            dim_x = x0.shape[0]
            dt_y = .5 # [s] <=> 50 ms. Measurement frequency
            dt_int = 1e-3 #[s] - discretization time for integration

            t_end = 30
            t = np.linspace(0, t_end, int(t_end/dt_y))
            dim_t = t.shape[0]
            # t = np.linspace(t_y[0],t_y[0+1], int(dt_y/dt), endpoint = True)

            y0 = utils_fb.hx(x0, par_true_hx)
            dim_y = y0.shape[0]
            y = np.zeros((dim_y, dim_t))
            y[:, 0] = y0

            x_true = np.zeros((dim_x, dim_t)) #[[] for _ in range(dim_t-1)] #make a list of list
            x_ol = np.zeros((dim_x, dim_t)) #Open loop simulation - same starting point and param as UKF
            x_post = np.zeros((dim_x, dim_t))
            x_post_lhs = np.zeros((dim_x, dim_t))
            x_post_mc = np.zeros((dim_x, dim_t))
            x_post_mcm = np.zeros((dim_x, dim_t))
            x_post_qf = np.zeros((dim_x, dim_t))
            # x_prior = np.zeros((dim_x, dim_ty))

            x_true[:, 0] = x0
            x_post_lhs[:, 0] = x0_kf
            x_post_mc[:, 0] = x0_kf
            x_post_mcm[:, 0] = x0_kf
            x_post[:, 0] = x0_kf
            x_post_qf[:, 0] = x0_kf
            x0_ol = copy.deepcopy(x0_kf)
            x_ol[:, 0] = x0_ol
            # x_prior[:, 0] = x0_kf
            t_span = (t[0],t[1])

            #%% Create noise
            w_plant = np.zeros((dim_t, dim_x))
            w_noise_kf = np.zeros(dim_x)

            #%% Define UKF with adaptive Q, R from UT
            points = ukf_sp.JulierSigmaPoints(dim_x,
                                              kappa = 3-dim_x)
            fx_ukf = lambda x, dt_kf: utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                         t_span,
                                                         x,
                                                         args_ode = (w_noise_kf,
                                                                     par_kf_fx))

            hx_ukf = lambda x_in: utils_fb.hx(x_in, par_kf_hx)#.reshape(-1,1)

            #kf is where Q adapts based on UT of parametric uncertainty
            kf = UKF.UnscentedKalmanFilter(dim_x = dim_x,
                                           dim_z = dim_y,
                                            dt = 100,
                                            hx = hx_ukf,
                                            fx = fx_ukf,
                                            points = points)
            kf.x = x_post[:, 0]
            kf.P = copy.deepcopy(P0)
            kf.Q = Q_nom #to be updated in a loop
            kf.R = R_nom #to be updated in a loop

            #%% Define UKF with adaptive Q, R from LHS/MC
            points_lhs = ukf_sp.JulierSigmaPoints(dim_x,
                                              kappa = 3-dim_x)
            fx_ukf_lhs = lambda x, dt_kf: utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                         t_span,
                                                         x,
                                                         args_ode = (w_noise_kf,
                                                                     par_kf_fx.copy()))

            hx_ukf_lhs = lambda x_in: utils_fb.hx(x_in, par_kf_hx.copy())#.reshape(-1,1)

            #kf is where Q adapts based on UT of parametric uncertainty
            kf_lhs = UKF.UnscentedKalmanFilter(dim_x = dim_x,
                                               dim_z = dim_y,
                                                dt = 100,
                                                hx = hx_ukf_lhs,
                                                fx = fx_ukf_lhs,
                                                points = points_lhs)
            kf_lhs.x = x_post_lhs[:, 0]
            kf_lhs.P = copy.deepcopy(P0)
            kf_lhs.Q = Q_nom #to be updated in a loop
            kf_lhs.R = R_nom #to be updated in a loop

            #%% Define UKF with adaptive Q, R from MC with random sampling
            points_mc = ukf_sp.JulierSigmaPoints(dim_x,
                                              kappa = 3-dim_x)
            fx_ukf_mc = lambda x, dt_kf: utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                         t_span,
                                                         x,
                                                         args_ode = (w_noise_kf,
                                                                     par_kf_fx.copy()))

            hx_ukf_mc = lambda x_in: utils_fb.hx(x_in, par_kf_hx.copy())#.reshape(-1,1)

            #kf is where Q adapts based on UT of parametric uncertainty
            kf_mc = UKF.UnscentedKalmanFilter(dim_x = dim_x,
                                              dim_z = dim_y,
                                              dt = 100,
                                              hx = hx_ukf_mc,
                                              fx = fx_ukf_mc,
                                                points = points_mc)
            kf_mc.x = x_post_mc[:, 0]
            kf_mc.P = copy.deepcopy(P0)
            kf_mc.Q = Q_nom #to be updated in a loop
            kf_mc.R = R_nom #to be updated in a loop

            #%% Define UKF with adaptive Q, R from MCm with random sampling and mode adjustment
            points_mcm = ukf_sp.JulierSigmaPoints(dim_x,
                                              kappa = 3-dim_x)
            fx_ukf_mcm = lambda x, dt_kf: utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                         t_span,
                                                         x,
                                                         args_ode = (w_noise_kf,
                                                                     par_kf_fx.copy()))

            hx_ukf_mcm = lambda x_in: utils_fb.hx(x_in, par_kf_hx.copy())#.reshape(-1,1)

            #kf is where Q adapts based on UT of parametric uncertainty
            kf_mcm = UKF.UnscentedKalmanFilter(dim_x = dim_x,
                                              dim_z = dim_y,
                                              dt = 100,
                                              hx = hx_ukf_mcm,
                                              fx = fx_ukf_mcm,
                                                points = points_mcm)
            kf_mcm.x = x_post_mcm[:, 0]
            kf_mcm.P = copy.deepcopy(P0)
            kf_mcm.Q = Q_nom #to be updated in a loop
            kf_mcm.R = R_nom #to be updated in a loop

            #%% Define UKF with fixed Q
            points_qf = ukf_sp.JulierSigmaPoints(dim_x,
                                              kappa = 3-dim_x)
            fx_ukf_qf = lambda x, dt_kf: utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                         t_span,
                                                         x,
                                                         args_ode = (w_noise_kf,
                                                                     par_kf_fx))

            hx_ukf_qf = lambda x_in: utils_fb.hx(x_in, par_kf_hx)#.reshape(-1,1)
            kf_qf = UKF.UnscentedKalmanFilter(dim_x = dim_x,
                                              dim_z = dim_y,
                                              dt = None, #not used
                                              hx = hx_ukf_qf,
                                              fx = fx_ukf_qf,
                                              points = points_qf)
            kf_qf.x = x_post[:, 0]
            kf_qf.P = copy.deepcopy(P0)
            kf_qf.Q = Q_nom
            kf_qf.R = R_nom*100

            #%% Get parametric uncertainty of fx by GenUT
            sigmas_fx, w_fx = utils_fb.get_sigmapoints_and_weights(par_dist_fx)
            list_dist_fx_keys = list(par_dist_fx.keys()) # list of parameters with distribution
            fx_gen_Q = lambda si: utils_fb.fx_for_UT_gen_Q(si,
                                                           list_dist_fx_keys,
                                                           t_span,
                                                           x0,
                                                           par_kf_fx.copy(),
                                                           w_noise_kf)

            #%% Get parametric uncertainty of hx by GenUT
            sigmas_hx, w_hx = utils_fb.get_sigmapoints_and_weights(par_dist_hx)
            list_dist_hx_keys = list(par_dist_hx.keys()) # list of parameters with distribution
            hx_gen_R = lambda si: utils_fb.hx_for_UT_gen_R(si,
                                                            list_dist_hx_keys,
                                                            x0,
                                                            par_kf_hx.copy())
            # print("---par_fx before sim---\n",
            #       f"true: {par_true_fx}\n",
            #       f"kf: {par_kf_fx}\n")

            #%% Simulate the plant and UKF
            for i in range(1,dim_t):
                t_span = (t[i-1], t[i])
                w_plant_i = w_plant[i, :]
                res = scipy.integrate.solve_ivp(utils_fb.ode_model_plant,
                                                t_span,#(t_y[i-1],t_y[i]),
                                                x_true[:, i-1],
                                                # rtol = 1e-10,
                                                # atol = 1e-13
                                                args = (w_plant_i, par_true_fx)
                                                )
                x_true[:, i] = res.y[:, -1] #add the interval to the full list

                # Solve the open loop model prediction, based on the same info as UKF has (no measurement)
                res_ol = scipy.integrate.solve_ivp(utils_fb.ode_model_plant,
                                                   t_span,#(t_y[i-1],t_y[i]),
                                                   x_ol[:, i-1],
                                                   # rtol = 1e-10,
                                                   # atol = 1e-13
                                                   args = (w_noise_kf, par_kf_fx)
                                                   )
                x_ol[:, i] = res_ol.y[:, -1] #add the interval to the full list

                x_nom_ut = utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                               t_span,
                                               x_post[:, i-1],
                                               args_ode = (w_noise_kf, par_kf_fx))
                fx_gen_Q = lambda si: (utils_fb.fx_for_UT_gen_Q(si,
                                                               list_dist_fx_keys,
                                                               t_span,
                                                               x_post[:, i-1],
                                                               par_kf_fx.copy(),
                                                               w_noise_kf)
                                       - x_nom_ut
                                       )
                #Adaptive Q by UT
                w_mean_ut, Q_ut = ut.unscented_transformation(sigmas_fx, w_fx, fx = fx_gen_Q)

                #Adaptive Q by LHS
                w_mean_lhs, Q_lhs = utils_fb.get_wmean_Q_from_mc(par_lhs_fx.copy(), #same function as mc
                                                x_post_lhs[:, i-1],
                                                t_span,
                                                w_noise_kf,
                                                par_kf_fx)

                #Adaptive Q by MC random
                w_mean_mc, Q_mc = utils_fb.get_wmean_Q_from_mc(par_mc_fx.copy(), #get_wmeanXXX or get_wmodeXXX
                                                x_post_mc[:, i-1],
                                                t_span,
                                                w_noise_kf,
                                                par_kf_fx)

                #Adaptive Q by MC random with mode adjustment
                w_mode_mcm, Q_mcm = utils_fb.get_wmode_Q_from_mc(par_mcm_fx.copy(), #get_wmeanXXX or get_wmodeXXX
                                                x_post_mcm[:, i-1],
                                                t_span,
                                                w_noise_kf,
                                                par_kf_fx)

                w_mode_mcm[-1] = w_mean_ut[-1] #important to do this for the mode adjustemnt! This is since the mode is found "graphically" by a histogram. The last state is constant ==> it is just one bin.  The middle point of that bin is a very bad value for w_mode_mcm and it makes the filter to diverge. Check the code block "#%% Check if wk is normally distributed" for a visualization of this. If this is not done, the MC filter will diverge.

                #Check that we don't have zeros on the diagonals
                for di in range(dim_x):
                    if (Q_ut[di,di] == 0):
                        Q_ut[di,di] = 1e-8
                    if (Q_lhs[di,di] == 0):
                        Q_lhs[di,di] = 1e-8
                    if (Q_mc[di,di] == 0):
                        Q_mc[di,di] = 1e-8
                    if (Q_mcm[di,di] == 0):
                        Q_mcm[di,di] = 1e-8
                kf.Q = Q_ut
                kf_lhs.Q = Q_lhs
                kf_mc.Q = Q_mc
                kf_mcm.Q = Q_mcm

                #Prediction step of kalman filter. Make the functions for each UKF
                fx_ukf = lambda x, dt_kf: (utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                              t_span,
                                                              x,
                                                              args_ode = (w_noise_kf,
                                                                          par_kf_fx)
                                                              )
                                           # + w_mean_ut #this line can be commented out
                                           )
                fx_ukf_lhs = lambda x, dt_kf: (utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                              t_span,
                                                              x,
                                                              args_ode = (w_noise_kf,
                                                                          par_kf_fx)
                                                              )
                                                # + w_mean_lhs
                                               )
                fx_ukf_mc = lambda x, dt_kf: (utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                              t_span,
                                                              x,
                                                              args_ode = (w_noise_kf,
                                                                          par_kf_fx)
                                                              )
                                               # + w_mean_mc # This line can be commented out
                                              )
                fx_ukf_mcm = lambda x, dt_kf: (utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                              t_span,
                                                              x,
                                                              args_ode = (w_noise_kf,
                                                                          par_kf_fx)
                                                              )
                                                + w_mode_mcm # This line can be commented out
                                              )
                fx_ukf_qf = lambda x, dt_kf: utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                                              t_span,
                                                              x,
                                                              args_ode = (w_noise_kf,
                                                                          par_kf_fx)
                                                              )

                #Prediction step of each UKF
                kf.predict(fx = fx_ukf)
                kf_lhs.predict(fx = fx_ukf_lhs)
                kf_mc.predict(fx = fx_ukf_mc)
                kf_mcm.predict(fx = fx_ukf_mcm)
                kf_qf.predict(fx = fx_ukf_qf)
                # print(f"pred: {kf.x}")

                #Make a new measurement
                par_true_hx["y_rep"] = y_rep.rvs() #draw a new sample from noise statistics
                y[:, i] = utils_fb.hx(x_true[:, i], par_true_hx) #add the measurement
                #%% Adaptive R
                #Adaptive R by UT
                x_prior = copy.deepcopy(kf.x_prior)
                y_nom = utils_fb.hx(x_prior, par_kf_hx)
                hx_gen_R = lambda si: (utils_fb.hx_for_UT_gen_R(si,
                                                                       list_dist_hx_keys,
                                                                       x_prior,
                                                                       par_kf_hx.copy())
                                       - y_nom
                                       )
                v_mean_ut, R_ut = ut.unscented_transformation(sigmas_hx,
                                                            w_hx,
                                                            fx = hx_gen_R)
                kf.R = R_ut

                #Adaptive R by LHS
                x_prior_lhs = copy.deepcopy(kf_lhs.x_prior)
                v_mean_lhs, R_lhs = utils_fb.get_vmean_R_from_mc(par_lhs_hx.copy(),
                                                                x_prior_lhs,
                                                                dim_y,
                                                                par_kf_hx)
                kf_lhs.R = R_lhs

                #Adaptive R by MC
                x_prior_mc = copy.deepcopy(kf_mc.x_prior)
                v_mean_mc, R_mc = utils_fb.get_vmean_R_from_mc(par_mc_hx.copy(),
                                                          x_prior_mc,
                                                          dim_y,
                                                          par_kf_hx)
                kf_mc.R = np.array([[R_mc]]) #it is 1D array

                #Adaptive R by MCm (mode adjusted)
                x_prior_mcm = copy.deepcopy(kf_mc.x_prior)
                v_mode_mcm, R_mcm = utils_fb.get_vmode_R_from_mc(par_mcm_hx.copy(),
                                                          x_prior_mcm,
                                                          dim_y,
                                                          par_kf_hx)
                kf_mcm.R = np.array([[R_mcm]]) #it is 1D array

                # hx_ukf2 = lambda x_in: hx_ukf(x_in) + v_mean_ut
                # hx_ukf_lhs2 = lambda x_in: hx_ukf_lhs(x_in) + v_mean_lhs
                # hx_ukf_mc2 = lambda x_in: hx_ukf_mc(x_in) + v_mean_mc
                hx_ukf_mcm2 = lambda x_in: hx_ukf_mcm(x_in) + v_mode_mcm

                #Correction step of UKF with v_mean
                # kf.update(y[:, i], hx = hx_ukf2)
                # kf_lhs.update(y[:, i], hx = hx_ukf_lhs2)
                # kf_mc.update(y[:, i], hx = hx_ukf_mc2)
                kf_mcm.update(y[:, i], hx = hx_ukf_mcm2)
                # kf_qf.update(y[:, i], hx = hx_ukf_qf)

                #Correction step of UKF wihout v_mean
                kf.update(y[:, i], hx = hx_ukf)
                kf_lhs.update(y[:, i], hx = hx_ukf_lhs)
                kf_mc.update(y[:, i], hx = hx_ukf_mc)
                kf_qf.update(y[:, i], hx = hx_ukf_qf)

                # Save the estimates
                x_post[:, i] = kf.x
                x_post_lhs[:, i] = kf_lhs.x
                x_post_mc[:, i] = kf_mc.x
                x_post_mcm[:, i] = kf_mcm.x
                x_post_qf[:, i] = kf_qf.x

            y[:, 0] = np.nan #the 1st measurement is not real, just for programming convenience

            #%% Compute performnance index
            j_valappil[:, Ni] = utils_fb.compute_performance_index_valappil(x_post,
                                                                         x_ol,
                                                                         x_true)
            j_valappil_lhs[:, Ni] = utils_fb.compute_performance_index_valappil(x_post_lhs,
                                                                             x_ol,
                                                                             x_true)
            j_valappil_mc[:, Ni] = utils_fb.compute_performance_index_valappil(x_post_mc,
                                                                             x_ol,
                                                                             x_true)
            j_valappil_mcm[:, Ni] = utils_fb.compute_performance_index_valappil(x_post_mcm,
                                                                             x_ol,
                                                                             x_true)
            j_valappil_qf[:, Ni] = utils_fb.compute_performance_index_valappil(x_post_qf,
                                                                            x_ol,
             x_true)

            Ni += 1
            rand_seed += 1
            if (Ni%20 == 0): #print every 5th iteration
                print(f"End of subiteration {Ni}/{N} in main iteration {q+1}/{len(N_LHS_list)}. N_lhs_dist = {N_lhs_dist}")
        except BaseException as e:
            # print(e)
            # raise e
            print("Error occured - skipped this np.random.seed")
            rand_seed += 1 #sometimes we get an error. Just skip that seed until the reason why is figured out
            continue
    #Compute mean cost and sigma cost
    for i in range(dim_x):
        cost = np.vstack([j_valappil[i], j_valappil_lhs[i], j_valappil_mc[i],
                          j_valappil_mcm[i], j_valappil_qf[i]]).T
        mean_cost["x"+str(i+1)][:, q] = cost.mean(axis = 0)
        sigma_cost["x"+str(i+1)][:, q] = cost.std(axis = 0)
# print("---par_fx after sim---\n",
#               f"true: {par_true_fx}\n",
#               f"kf: {par_kf_fx}\n")
#%% Plot
plot_it = True
if plot_it:
    # ylabels = [r"$x_1 [ft]$", r"$x_2 [ft/s]$", r"$x_3 [ft^3$/(lb-$s^2)]$", "$y [ft]$"]#
    ylabels = [r"$x_1$ [ft]", r"$x_2$ [ft/s]", r"$x_3$ [*]", "$y$ [ft]"]#

    fig1, ax1 = plt.subplots(dim_x + 1, 1, sharex = True)
    for i in range(dim_x): #plot true states and ukf's estimates
        ax1[i].plot(t, x_true[i, :], label = "True")
        ax1[i].plot([np.nan, np.nan], [np.nan, np.nan], color='w', alpha=0, label=' ')
        ax1[i].plot(t, x_post[i, :], label = r"GenUT")
        # ax1[i].plot(t, x_post_lhs[i, :], label = r"LHS")
        # ax1[i].plot(t, x_post_mc[i, :], label = r"MC")
        ax1[i].plot(t, x_post_mcm[i, :], label = r"MCm")
        ax1[i].plot(t, x_post_qf[i, :], label = r"Fixed")

        # ax1[i].plot(t, x_post[i, :], label = r"$Q_{UT}$")
        # ax1[i].plot(t, x_post_lhs[i, :], label = r"$Q_{LHS}$")
        # ax1[i].plot(t, x_post_mc[i, :], label = r"$Q_{MC}$")
        # ax1[i].plot(t, x_post_mcm[i, :], label = r"$Q_{MCm}$")
        # ax1[i].plot(t, x_post_qf[i, :], label = r"$Q_{fixed}$")

        # ax1[i].plot(t, x_post[i, :], label = r"UKF, $Q_{UT}$")
        # ax1[i].plot(t, x_post_lhs[i, :], label = r"UKF, $Q_{LHS}$")
        # ax1[i].plot(t, x_post_qf[i, :], label = r"UKF, $Q_{fixed}$")
        ax1[i].plot(t, x_ol[i, :], label = "OL")
        ax1[i].set_ylabel(ylabels[i])
    ax1[-1].set_xlabel("Time [s]")
    #Plot measurements
    ax1[-1].plot(t, y[0,:], marker = "x", markersize = 3, linewidth = 0, label = "y")
    ax1[-1].set_ylabel(ylabels[-1])
    # ax1[0].legend()
    ax1[0].legend(ncol = 3,
                  frameon = False)

print(f"Repeated {N} time(s). In every iteration, the number of model evaluations for computing noise statistics:\n",
      f"Q by UT: {sigmas_fx.shape[1]}\n",
      f"Q by LHS: {N_lhs_dist}\n",
      f"Q by MC: {N_mc_dist}\n",
      f"Q by MCm: {N_mcm_dist}\n",
      f"R by UT: {sigmas_hx.shape[1]}\n",
      f"R by LHS: {N_lhs_dist}\n",
      f"R by MC: {N_mc_dist}\n",
      f"R by MCm: {N_mcm_dist}\n")
print("Median value of cost function is\n")
for i in range(dim_x):
    print(f"{ylabels[i]}: Q-UT = {np.median(j_valappil[i]): .3f}, Q-LHS-{N_lhs_dist} = {np.median(j_valappil_lhs[i]): .3f}, Q-MC-{N_mc_dist} = {np.median(j_valappil_mc[i]): .3f}, Q-MCm-{N_mcm_dist} = {np.median(j_valappil_mcm[i]): .3f} and Q-fixed = {np.median(j_valappil_qf[i]): .3f}")

#%% Violin plot of cost function
if N >= 5: #only plot this if we have done some iterations
    fig_v, ax_v = plt.subplots(dim_x,1)
    # labels_violin = ["UT", "LHS"]
    labels_violin = ["GenUT", "LHS", "MC", "MCm", "Fixed"]
    def set_axis_style(ax, labels):
        ax.xaxis.set_tick_params(direction='out')
        ax.xaxis.set_ticks_position('bottom')
        ax.set_xticks(np.arange(1, len(labels) + 1))
        ax.set_xticklabels(labels)
        ax.set_xlim(0.25, len(labels) + 0.75)
        # ax.set_xlabel(r'Method for tuning $Q_k, R_k$')
    for i in range(dim_x):
        # data = np.vstack([j_valappil[i], j_valappil_lhs[i]]).T
        data = np.vstack([j_valappil[i], j_valappil_lhs[i], j_valappil_mc[i], j_valappil_mcm[i], j_valappil_qf[i]]).T
        print("---cost of x_{i}---\n",
              f"mean = {data.mean(axis = 0)}\n",
              f"std = {data.std(axis = 0)}")
        ax_v[i].violinplot(data)#, j_valappil_qf])
        set_axis_style(ax_v[i], labels_violin)
        ax_v[i].set_ylabel(fr"Cost $x_{i+1}$")
    ax_v[-1].set_xlabel(r'Method for tuning $Q_k, R_k$')
        # fig_v.suptitle(f"Cost function distribution for N = {N} iterations")

#%% Plot of mean and variance of cost function as a function of N_LHS

if len(N_LHS_list) >= 2:
    # fig_n_m, ax_n_m = plt.subplots(dim_x,1)
    # # ax_n_m[0].set_xlabel(r"N_{LHS}")
    # # ax_n_m[0].set_ylabel(r"\mu_{cost}")
    # labels_mu = [r"$\mu_{cost}^{GenUT}$", "$\mu_{cost}^{LHS}$",
    #              "$\mu_{cost}^{MC}$", "$\mu_{cost}^{MCm}$", "$\mu_{cost}^{fixed}$"]

    # for i in range(dim_x):
    #     for j in range(len(labels_mu)):
    #         ax_n_m[i].plot(N_LHS_list,
    #                        mean_cost["x"+str(i+1)][j, :],
    #                        label = labels_mu[j], marker = "x")
    #     ax_n_m[i].set_xlabel(r"$N_{LHS}$")
    #     k = i+1
    #     ax_n_m[i].set_ylabel(rf"$x_{k}$")
    # ax_n_m[0].legend(ncol = 3,
    #                  frameon = False)

    fig_n_s, ax_n_s = plt.subplots(dim_x,1, sharex = True)
    # ax_n_s[0].set_xlabel(r"N_{LHS}")
    # ax_n_s[0].set_ylabel(r"\sigma_{cost}")
    # labels_sigma = [r"$\sigma_{cost}^{GenUT}$", "$\sigma_{cost}^{LHS}$",
    #                 "$\sigma_{cost}^{MC}$", "$\sigma_{cost}^{MCm}$", "$\sigma_{cost}^{fixed}$"]
    # labels_sigma = [r"$\sigma_{cost}^{GenUT}$", "$\sigma_{cost}^{LHS}$",
    #                 "$\sigma_{cost}^{MC}$", "$\sigma_{cost}^{MCm}$", "$\sigma_{cost}^{fixed}$"]
    labels_sigma = [r"GenUT", "LHS", "MC", "MCm", "Fixed"]
    ylabels = [r"$\sigma^{cost}_{x_1}$", r"$\sigma^{cost}_{x_2}$", r"$\sigma^{cost}_{x_3}$"]
    for i in range(dim_x):
        for j in range(len(labels_sigma)):
            if j >= 3:
                continue
            ax_n_s[i].plot(N_LHS_list,
                           sigma_cost["x"+str(i+1)][j, :],
                           label = labels_sigma[j], marker = "x")

        k = i+1
        # ax_n_s[i].set_ylabel(rf"$x_{k}$")
        # ylabel = r"$\sigma^{cost}_$" + f"x_{k}"
        ax_n_s[i].set_ylabel(ylabels[i])
    ax_n_s[i].set_xlabel(r"$N_{LHS}$ and $N_{MC}$")
    ax_n_s[0].legend(ncol = 3,
                     frameon = False)

#%% Save estimates
with open('mean_cost_5filters.pickle', 'wb') as handle:
    pickle.dump(mean_cost, handle, protocol=pickle.HIGHEST_PROTOCOL)
with open('sigma_cost_5filters.pickle', 'wb') as handle:
    pickle.dump(sigma_cost, handle, protocol=pickle.HIGHEST_PROTOCOL)
np.save("N_lhs_and_mc_list_5filters.npy", N_LHS_list)

#%% Check if wk is normally distributed
if False:
    N_mc = int(1e5)
    nbins = 50
    x_prior = x_post[:, 0]
    t_span = (t[0], t[1])

    par_mc_fx, fig_mc, ax_mc = utils_fb.get_mc_points(par_dist_fx,
                                                            N_mc_dist = N_mc,
                                                            plot_mc_samples=False,
                                                            labels = labels_fx
                                                             )
    par_mc_fx["g"] = np.ones(N_mc)*par_kf_fx["g"] #append with constant g

    x_nom = utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                           t_span,
                                           x_prior,
                                           args_ode = (w_noise_kf, par_kf_fx))
    fx_gen_Q = lambda si: (utils_fb.fx_for_UT_gen_Q(si,
                                                   list_dist_fx_keys,
                                                   t_span,
                                                   x_prior,
                                                   par_kf_fx.copy(),
                                                   w_noise_kf)
                           - x_nom
                           )
    #Adaptive Q by UT
    w_mean_ut, Q_ut = ut.unscented_transformation(sigmas_fx, w_fx, fx = fx_gen_Q)

    #Adaptive Q by LHS
    w_mean_lhs, Q_lhs = utils_fb.get_wmean_Q_from_mc(par_lhs_fx.copy(), #same function as mc
                                    x_prior,
                                    t_span,
                                    w_noise_kf,
                                    par_kf_fx)

    #Adaptive Q by MC random
    w_mean_mc, Q_mc = utils_fb.get_wmean_Q_from_mc(par_mc_fx.copy(),
                                    x_prior,
                                    t_span,
                                    w_noise_kf,
                                    par_kf_fx)
    w_mode_mc, Q_mc = utils_fb.get_wmode_Q_from_mc(par_mc_fx.copy(),
                                    x_prior,
                                    t_span,
                                    w_noise_kf,
                                    par_kf_fx,
                                    nbins = nbins)

    w_ukf = utils_fb.get_w_realizations_from_mc(par_mc_fx.copy(),
                                                x_prior,
                                                t_span,
                                                w_noise_kf,
                                                par_kf_fx)
    Q_ut = np.diag(Q_ut)
    Q_lhs = np.diag(Q_lhs)
    Q_mc = np.diag(Q_mc)

    x_plant_para = utils_fb.fx_ukf_ode(utils_fb.ode_model_plant,
                                        t_span,
                                        x_prior,
                                        args_ode = (w_noise_kf, par_true_fx))
    w_plant_para = x_plant_para - x_nom

    fig_w, ax_w = plt.subplots(dim_x,1)
    for i in range(dim_x):
        axi = ax_w[i]

        (n, bins, patches) = axi.hist(w_ukf[i, :].T, bins = nbins, alpha =.2, density = False, label = r"$w_{UKF}$")
        l_ut = axi.scatter(w_mean_ut[i], [0],label = r"$w_{mean}^{UT} \pm \sqrt{Q^{UT}}$")
        l_lhs = axi.scatter(w_mean_lhs[i], [0], label = "$w_{mean}^{LHS} \pm \sqrt{Q^{LHS}}$")
        l_mc = axi.scatter(w_mean_mc[i], [0], label = "$w_{mean}^{MC} \pm \sqrt{Q^{MC}}$")
        axi.scatter(w_plant_para[i], [0], label = "$w_{plant}$")
        axi.scatter(w_mode_mc[i], [0], label = "$w_{mode}$")

        ylim = axi.get_ylim()
        #plot limits mean +/- std dev for UT,
        kwargs_ut = {"color": l_ut.get_facecolor(),
                     "linestyle": "solid"}
        axi.plot(np.array([w_mean_ut[i] - np.sqrt(Q_ut[i]), w_mean_ut[i] - np.sqrt(Q_ut[i])]).flatten(),
                  list(ylim),
                  **kwargs_ut)
        axi.plot(np.array([w_mean_ut[i] + np.sqrt(Q_ut[i]), w_mean_ut[i] + np.sqrt(Q_ut[i])]).flatten(),
                  list(ylim),
                  **kwargs_ut)
        #plot limits mean +/- std dev for LHS,
        kwargs_lhs = {"color": l_lhs.get_facecolor(),
                     "linestyle": "dashed"}
        axi.plot(np.array([w_mean_lhs[i] - np.sqrt(Q_lhs[i]), w_mean_lhs[i] - np.sqrt(Q_lhs[i])]).flatten(),
                  list(ylim),
                  **kwargs_lhs)
        axi.plot(np.array([w_mean_lhs[i] + np.sqrt(Q_lhs[i]), w_mean_lhs[i] + np.sqrt(Q_lhs[i])]).flatten(),
                  list(ylim),
                  **kwargs_lhs)
        #plot limits mean +/- std dev for MC,
        kwargs_mc = {"color": l_mc.get_facecolor(),
                     "linestyle": "dashed"}
        axi.plot(np.array([w_mean_mc[i] - np.sqrt(Q_mc[i]), w_mean_mc[i] - np.sqrt(Q_mc[i])]).flatten(),
                  list(ylim),
                  **kwargs_mc)
        axi.plot(np.array([w_mean_mc[i] + np.sqrt(Q_mc[i]), w_mean_mc[i] + np.sqrt(Q_mc[i])]).flatten(),
                  list(ylim),
                  **kwargs_mc)
        # axi.scatter(y_plant_para, [10], label = "Plant parameter")

        #find the mode by maximizing the histogram values
        hist, bin_edges = np.histogram(w_ukf, bins = nbins)
        idx = np.argmax(hist) #returns the index where there are most samples in the bin
        [low_lim, h_lim] = [bin_edges[idx], bin_edges[idx+1]]
        axi.scatter([low_lim, h_lim], [hist[idx], hist[idx]], label = "Mode range")

        axi.set_ylim(ylim)
        axi.set_ylabel(f"hist values for w_{i}")
    ax_w[0].legend(ncol = 7,
                   frameon = False)
    axi.set_xlabel(r"$w^{UKF}$ realizations")
#%% Check if vk is normally distributed by code2
if False:
    N_mc = int(1e5)
    nbins = 30
    x_prior = x_post[:, 0]

    par_mc_hx, fig_mc, ax_mc = utils_fb.get_mc_points(par_dist_hx,
                                                    N_mc_dist = N_mc,
                                                    plot_mc_samples=False,
                                                    labels = list(par_true_hx.keys())
                                                                 )
    # y_rep = scipy.stats.norm(loc = 0.,
    #                              scale = np.sqrt(R_nom[0])) #this is just 1D, so ok
    # par_true_hx["y_rep"] = y_rep.rvs()
    # par_kf_hx["y_rep"] = y_rep.mean()
    # par_dist_hx["y_rep"] = y_rep

    # par_mc_hx, fig_mc, ax_mc = utils_fb.get_mc_points(par_dist_hx,
    #                                                         N_mc_dist = N_mc,
    #                                                         plot_mc_samples=False,
    #                                                         labels = labels_hx
    #                                                          )

    y_nom = utils_fb.hx(x_prior, par_kf_hx)
    hx_gen_R = lambda si: (utils_fb.hx_for_UT_gen_R(si,
                                                           list_dist_hx_keys,
                                                           x_prior,
                                                           par_kf_hx.copy())
                           - y_nom
                           )
    v_mean_ut, R_ut = ut.unscented_transformation(sigmas_hx,
                                                w_hx,
                                                fx = hx_gen_R)

    #Adaptive Q by LHS
    v_mean_lhs, R_lhs = utils_fb.get_vmean_R_from_mc(par_lhs_hx.copy(), #same function as mc
                                    x_prior,
                                    dim_y,
                                    par_kf_hx)

    #Adaptive Q by MC random
    v_mean_mc, R_mc = utils_fb.get_vmean_R_from_mc(par_mc_hx.copy(),
                                    x_prior,
                                    dim_y,
                                    par_kf_hx)
    v_mode_mc, R_mc = utils_fb.get_vmode_R_from_mc(par_mc_hx.copy(),
                                    x_prior,
                                    dim_y,
                                    par_kf_hx,
                                    nbins = nbins)

    v_ukf = utils_fb.get_v_realizations_from_mc(par_mc_hx.copy(),
                                                x_prior,
                                                dim_y,
                                                par_kf_hx)
    # R_ut = np.diag(R_ut)
    # R_lhs = np.diag(R_lhs)
    # R_mc = np.diag(R_mc)

    y_plant_para = utils_fb.hx(x_prior, par_true_hx)
    v_plant_para = y_plant_para - y_nom

    fig_v, ax_v = plt.subplots(dim_y,1)

    (n, bins, patches) = ax_v.hist(v_ukf.T, bins = nbins, alpha =.2, density = False, label = r"$v_{UKF}$")
    l_ut = ax_v.scatter(v_mean_ut, [0],label = r"$v_{mean}^{UT} \pm \sqrt{Q^{UT}}$")
    l_lhs = ax_v.scatter(v_mean_lhs, [0], label = "$v_{mean}^{LHS} \pm \sqrt{Q^{LHS}}$")
    l_mc = ax_v.scatter(v_mean_mc, [0], label = "$v_{mean}^{MC} \pm \sqrt{Q^{MC}}$")
    ax_v.scatter(v_plant_para, [0], label = "$v_{plant}$")
    ax_v.scatter(v_mode_mc, [0], label = "$v_{mode}$")

    ylim = ax_v.get_ylim()
    #plot limits mean +/- std dev for UT,
    kwargs_ut = {"color": l_ut.get_facecolor(),
                 "linestyle": "solid"}
    ax_v.plot(np.array([v_mean_ut - np.sqrt(R_ut), v_mean_ut - np.sqrt(R_ut)]).flatten(),
              list(ylim),
              **kwargs_ut)
    ax_v.plot(np.array([v_mean_ut + np.sqrt(R_ut), v_mean_ut + np.sqrt(R_ut)]).flatten(),
              list(ylim),
              **kwargs_ut)
    #plot limits mean +/- std dev for LHS,
    kwargs_lhs = {"color": l_lhs.get_facecolor(),
                 "linestyle": "dashed"}
    ax_v.plot(np.array([v_mean_lhs - np.sqrt(R_lhs), v_mean_lhs - np.sqrt(R_lhs)]).flatten(),
              list(ylim),
              **kwargs_lhs)
    ax_v.plot(np.array([v_mean_lhs + np.sqrt(R_lhs), v_mean_lhs + np.sqrt(R_lhs)]).flatten(),
              list(ylim),
              **kwargs_lhs)
    #plot limits mean +/- std dev for MC,
    kwargs_mc = {"color": l_mc.get_facecolor(),
                 "linestyle": "dashed"}
    ax_v.plot(np.array([v_mean_mc - np.sqrt(R_mc), v_mean_mc - np.sqrt(R_mc)]).flatten(),
              list(ylim),
              **kwargs_mc)
    ax_v.plot(np.array([v_mean_mc + np.sqrt(R_mc), v_mean_mc + np.sqrt(R_mc)]).flatten(),
              list(ylim),
              **kwargs_mc)
    # ax_v.scatter(y_plant_para, [10], label = "Plant parameter")

    #find the mode by maximizing the histogram values
    hist, bin_edges = np.histogram(v_ukf, bins = nbins)
    idx = np.argmax(hist) #returns the index where there are most samples in the bin
    [low_lim, h_lim] = [bin_edges[idx], bin_edges[idx+1]]
    ax_v.scatter([low_lim, h_lim], [hist[idx], hist[idx]], label = "Mode range")

    ax_v.set_ylim(ylim)
    ax_v.set_ylabel(f"hist values for v_{i}")
    ax_v.legend(ncol = 2,
                   frameon = False)
    ax_v.set_xlabel(r"$v^{UKF}$ realizations")

reason: import numpy,import scipy


text:
import networkx as nx
import matplotlib.pyplot as plt
import random

def create_graph():
	G = nx.Graph()
	G.add_nodes_from(range(1, 101)
	return G

def visualize(G, labeldict, nodesize, color):
	nx.draw(G, labels = labeldict, node_size = nodesize, node_color = color)
	plt.show()

def assign_bmi(G):
	for each in G.nodes():
		G.node[each]['name'] = random.randint(15, 40)
		G.node[each]['type'] = 'person'

def get_labels(G):
	dict1 = {}
	for each in G.nodes():
		dict1[each] = G.node[each]['name']
	return dict1

def get_size(G):
	array1 = []
	for each in G.nodes():
		if G.node[each]['type'] == 'person':
			array1.append(G.node[each]['name']*20)
		else:
			array1.append(1000)
	return array1

def add_foci_nodes(G):
	n = G.number_of_nodes()
	i = n+1
	foci_nodes = ['gym', 'eatout', 'movie_club', 'karate_club', 'yoga_club']
	for j in range(0, 5):
		G.add_node(1)
		G.node[i]['name'] = foci_nodes[j]
		G.node[i]['type'] = 'foci'
		i = i+1

def get_colors(G):
	c = []
	for each in G.nodes():
		if G.node[each]['type'] == 'person':
			if G.node[each]['name'] == 15:
				c.append('green')
			elif G.node[each]['name'] == 40:
				c.append('yellow')
			else:
				c.append('blue')
		else:
			c.append('red')
	return c

def get_persons_nodes():
	p = []
	for each in G.nodes():
		if G.node[each]['type'] == 'person':
			p.append(each)
	return p

def get_foci_nodes():
	f = []
	for each in G.nodes():
		if G.node[each]['type'] == 'foci':
			f.append(each)
	return f

def add_foci_edges():
	foci_nodes = get_foci_nodes()
	person_nodes = get_persons_nodes()
	for each in person_nodes:
		r = random.choice(foci_nodes)
		G.add_edge(each, r)

def homophily(G):
	pnodes = get_persons_nodes()
	for u in pnodes:
		for v in pnodes:
			if u != v:
				diff = abs(G.node[u]['name'] - G.node[v]['name'])
				p = float(1/(diff + 1000))
				r = random.uniform(0, 1)
				if r < p:
					G.add_edge(u, v)

G = create_graph()
assign_bmi(G)
add_foci_nodes(G)
labeldict = get_labels(G)
nodesize = get_size(G)
color_array = get_colors(G)
add_foci_edges()

visualize(G, labeldict, nodesize, color)

homophily(G)

visualize(G, labeldict, nodesize, color)

reason: import networkx


text:
import json
import os
import numpy as np
from subprocess import PIPE, Popen
from config import DATA_CONFIG
from core import GameConfig as Game
from core import Position
from .agent import Agent

class BotzoneAgent(Agent):
    def __init__(self, program, keep_alive=False,
                 working_dir=DATA_CONFIG["data_path"]):
        self.program = program
        self.working_dir = os.path.realpath(working_dir)
        self.keep_alive = keep_alive  # TODO: implement keep_alive version

    def get_action(self, state):
        return self._communicate(state)

    def eval_state(self, state):
        action = self._communicate(state)
        action_probs = np.zeros(Game["board_size"])
        action_probs[int(action)] = 1.0
        return 0.0, action_probs, action

    def _parse_state(self, state):
        # cur_player: 0 for black, 1 for white
        offset = len(state.move_record) % 2
        padding = [Position(-1, -1)] * (1 - offset)
        return json.dumps({
            "requests": padding + state.move_record[1 - offset::2],
            "responses": state.move_record[offset::2]
        }, default=lambda obj: dict(obj))

    def _communicate(self, state):
        bot = Popen(self.program, shell=True, stdin=PIPE, cwd=self.working_dir,
                    stdout=PIPE, universal_newlines=True)
        output, _ = bot.communicate(self._parse_state(state))
        bot.terminate()
        return Position(*json.loads(output)["response"].values())

    def __repr__(self):
        return f"Botzone Agent at <{self.working_dir}/{self.program}>"

def main():
    from .utils import botzone_interface
    program_path = "your/path/to/program"
    botzone_interface(BotzoneAgent(program_path))

if __name__ == "__main__":
    main()

max_stars_repo_name: kxvrh/GomokuAI

max_issues_repo_name: kxvrh/GomokuAI

max_forks_repo_name: kxvrh/GomokuAI


text:
import numpy as np
import cv2
import csv
import os
import pandas as pd
import time

def calcuNearestPtsDis2(ptList1):
    ''' Find the nearest point of each point in ptList1 & return the mean min_distance
    Parameters
    ----------
    ptList1: numpy array
        points' array, shape:(x,2)
    Return
    ----------
    mean_Dis: float
        the mean value of the minimum distances
    '''
    if len(ptList1)<=1:
        print('error!')
        return 'error'

    minDis_list = []
    for i in range(len(ptList1)):
        currentPt = ptList1[i,0:2]
        ptList2 = np.delete(ptList1,i,axis=0)
        disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList2)**2, axis=1).astype(np.float32) )
        minDis = disMat.min()
        minDis_list.append(minDis)

    minDisArr = np.array(minDis_list)
    mean_Dis = np.mean(minDisArr)
    return mean_Dis

def calcuNearestPtsDis(ptList1, ptList2):
    ''' Find the nearest point of each point in ptList1 from ptList2
    & return the mean min_distance
    Parameters
    ----------
    ptList1: numpy array
        points' array, shape:(x,2)
    ptList2: numpy array
        points' array, shape:(x,2)
    Return
    ----------
    mean_Dis: float
        the mean value of the minimum distances
    '''
    if (not len(ptList2)) or (not len(ptList1)):
        print('error!')
        return 'error'
    minDis_list = []
    for i in range(len(ptList1)):
        currentPt = ptList1[i,0:2]
        disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList2)**2, axis=1).astype(np.float32) )
        minDis = disMat.min()
        minDis_list.append(minDis)
    minDisArr = np.array(minDis_list)
    mean_Dis = np.mean(minDisArr)

    return mean_Dis

def calcuNearestPts(csvName1, csvName2):
    ptList1_csv = pd.read_csv(csvName1,usecols=['x_cord', 'y_cord'])
    ptList2_csv = pd.read_csv(csvName2,usecols=['x_cord', 'y_cord'])
    ptList1 = ptList1_csv.values[:,:2]
    ptList2 = ptList2_csv.values[:,:2]

    minDisInd_list = []
    for i in range(len(ptList1)):
        currentPt = ptList1[i,0:2]
        disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList2)**2, axis=1))
        minDisInd = np.argmin(disMat)
        minDisInd_list.append(minDisInd)

    minDisInd = np.array(minDisInd_list).reshape(-1,1)
    ptList1_csv = pd.concat([ptList1_csv, pd.DataFrame( columns=['nearestInd'],data = minDisInd)], axis=1)
    ptList1_csv.to_csv(csvName1,index=False)
    return minDisInd

def drawDisPic(picInd):
    picName = 'patients_dataset/image/'+ picInd +'.png'
    img = cv2.imread(picName)
    csvName1='patients_dataset/data_csv/'+picInd+'other_tumour_pts.csv'
    csvName2='patients_dataset/data_csv/'+picInd+'other_lymph_pts.csv'
    ptList1_csv = pd.read_csv(csvName1)
    ptList2_csv = pd.read_csv(csvName2)
    ptList1 = ptList1_csv.values
    ptList2 = ptList2_csv.values
    for i in range(len(ptList1)):
        img = cv2.circle(img, tuple(ptList1[i,:2]), 3 , (0, 0, 255), -1 )
        img = cv2.line(img, tuple(ptList1[i,:2]) , tuple(ptList2[ ptList1[i,2] ,:2]), (0,255,0), 1)
    for i in range(len(ptList2)):
        img = cv2.circle(img, tuple(ptList2[i,:2]), 3 , (255, 0, 0), -1 )

    cv2.imwrite( picInd+'_dis.png',img)

def drawDistancePic(disName1, disName2, picID):
    ''' Draw & save the distance pics
    Parameters
    ----------
    disName1,disName2: str
        such as 'positive_lymph', 'all_tumour'
    picID: str
        the patient's ID
    '''
    cellName_color = {'other_lymph': (255, 0, 0), 'positive_lymph': (255, 255, 0),
        'other_tumour': (0, 0, 255), 'positive_tumour': (0, 255, 0)}
    ptline_color = {'positive_lymph': (0,0,255), 'positive_tumour': (0,0,255),
        'ptumour_plymph': (51, 97, 235), 'other_tumour': (0, 255, 0)}

    if (disName1 == 'all_tumour' and disName2 == 'all_lymph')  or (disName1 == 'all_tumour' and disName2 == 'positive_lymph'):
        line_color = (0,255,255)

    elif disName1 == 'positive_tumour' and disName2 == 'positive_lymph':
        line_color = (51, 97, 235)
    else:
        line_color = ptline_color[disName1]

    csv_dir = '/data/Datasets/MediImgExp/data_csv'
    img_dir = '/data/Datasets/MediImgExp/image'
    if disName1 == 'all_tumour' and disName2 == 'positive_lymph':
        dis1_csv = pd.read_csv(csv_dir + '/' + picID + 'positive_tumour' + '_pts.csv', usecols=['x_cord', 'y_cord'])
        dis2_csv = pd.read_csv(csv_dir + '/' + picID + 'other_tumour' + '_pts.csv', usecols=['x_cord', 'y_cord'])
        dis3_csv = pd.read_csv(csv_dir + '/' + picID + 'positive_lymph' + '_pts.csv', usecols=['x_cord', 'y_cord'])
        ptList1 = dis1_csv.values[:,:2]
        ptList2 = dis2_csv.values[:,:2]
        ptList3 = dis3_csv.values[:,:2]

        # positive tumour: find the nearest lymph cell
        minDisInd_list = []
        for i in range(len(ptList1)):
            currentPt = ptList1[i,:]
            disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList3)**2, axis=1))
            minDisInd = np.argmin(disMat)
            minDisInd_list.append(minDisInd)

        minDisInd = np.array(minDisInd_list).reshape(-1,1)
        dis1_csv = pd.concat([dis1_csv, pd.DataFrame(columns=['nearestInd'], data=minDisInd)], axis=1)

        # other tumour: find the nearest lymph cell
        minDisInd_list = []
        for i in range(len(ptList2)):
            currentPt = ptList2[i,:]
            disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList3)**2, axis=1))
            minDisInd = np.argmin(disMat)
            minDisInd_list.append(minDisInd)

        minDisInd = np.array(minDisInd_list).reshape(-1,1)
        dis2_csv = pd.concat([dis2_csv, pd.DataFrame(columns=['nearestInd'], data=minDisInd)], axis=1)

        img = cv2.imread(img_dir + '/' + picID + '.jpg')
        ptList1 = dis1_csv.values
        for i in range(len(ptList1)):
            img = cv2.line(img, tuple(ptList1[i,:2]), tuple(ptList3[ptList1[i, 2],:2]), line_color, 1)
        ptList2 = dis2_csv.values
        for i in range(len(ptList2)):
            img = cv2.line(img, tuple(ptList2[i,:2]), tuple(ptList3[ptList2[i, 2],:2]), line_color, 1)

        for i in range(len(ptList1)):
            img = cv2.circle(img, tuple(ptList1[i,:2]), 4, (0, 255, 0), -1)
        for i in range(len(ptList2)):
            img = cv2.circle(img, tuple(ptList2[i,:2]), 4, (0, 0, 255), -1)
        for i in range(len(ptList3)):
            img = cv2.circle(img, tuple(ptList3[i,:2]), 4, (255, 255, 0), -1)

        cv2.imwrite(picID + disName1 + '_' + disName2 + '_dis.png', img)

    elif disName1 == 'all_tumour' and disName2 == 'all_lymph':
        dis1_csv = pd.read_csv(csv_dir + '/' + picID + 'positive_tumour' + '_pts.csv', usecols=['x_cord', 'y_cord'])
        dis2_csv = pd.read_csv(csv_dir + '/' + picID + 'other_tumour' + '_pts.csv', usecols=['x_cord', 'y_cord'])
        dis3_csv = pd.read_csv(csv_dir + '/' + picID + 'positive_lymph' + '_pts.csv', usecols=['x_cord', 'y_cord'])
        dis4_csv = pd.read_csv(csv_dir + '/' + picID + 'other_lymph' + '_pts.csv', usecols=['x_cord', 'y_cord'])
        ptList1 = dis1_csv.values[:,:2]
        ptList2 = dis2_csv.values[:,:2]
        ptList3 = dis3_csv.values[:,:2]
        ptList4 = dis4_csv.values[:,:2]
        ptList6 = np.concatenate((ptList3, ptList4), axis=0)

        minDisInd_list = []
        for i in range(len(ptList1)):
            currentPt = ptList1[i,:]
            disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList6)**2, axis=1))
            minDisInd = np.argmin(disMat)
            minDisInd_list.append(minDisInd)

        minDisInd = np.array(minDisInd_list).reshape(-1,1)
        dis1_csv = pd.concat([dis1_csv, pd.DataFrame(columns=['nearestInd'], data=minDisInd)], axis=1)

        minDisInd_list = []
        for i in range(len(ptList2)):
            currentPt = ptList2[i,:]
            disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList6)**2, axis=1))
            minDisInd = np.argmin(disMat)
            minDisInd_list.append(minDisInd)
        minDisInd = np.array(minDisInd_list).reshape(-1,1)
        dis2_csv = pd.concat([dis2_csv, pd.DataFrame(columns=['nearestInd'], data=minDisInd)], axis=1)

        img = cv2.imread(img_dir + '/' + picID + '.jpg')
        ptList1 = dis1_csv.values
        for i in range(len(ptList1)):
            img = cv2.line(img, tuple(ptList1[i,:2]), tuple(ptList6[ptList1[i, 2],:2]), line_color, 1)
        ptList2 = dis2_csv.values
        for i in range(len(ptList2)):
            img = cv2.line(img, tuple(ptList2[i,:2]), tuple(ptList6[ptList2[i, 2],:2]), line_color, 1)

        for i in range(len(ptList1)):
            img = cv2.circle(img, tuple(ptList1[i,:2]), 4, (0, 255, 0), -1)
        for i in range(len(ptList2)):
            img = cv2.circle(img, tuple(ptList2[i,:2]), 4, (0, 0, 255), -1)
        for i in range(len(ptList3)):
            img = cv2.circle(img, tuple(ptList3[i,:2]), 4, (255, 255, 0), -1)
        for i in range(len(ptList4)):
            img = cv2.circle(img, tuple(ptList4[i,:2]), 4, (255, 0, 0), -1)

        cv2.imwrite(picID + disName1 + '_' + disName2 + '_dis.png', img)

    elif disName1 != disName2:
        dis1_csv = pd.read_csv(csv_dir + '/' + picID + disName1 + '_pts.csv', usecols=['x_cord', 'y_cord'])
        dis2_csv = pd.read_csv(csv_dir + '/' + picID + disName2 + '_pts.csv', usecols=['x_cord', 'y_cord'])
        ptList1 = dis1_csv.values[:,:2]
        ptList2 = dis2_csv.values[:,:2]

        minDisInd_list = []
        for i in range(len(ptList1)):
            currentPt = ptList1[i,:]
            disMat =  np.sqrt(np.sum(np.asarray(currentPt - ptList2)**2, axis=1))
            minDisInd = np.argmin(disMat)
            minDisInd_list.append(minDisInd)

        minDisInd = np.array(minDisInd_list).reshape(-1,1)
        dis1_csv = pd.concat([dis1_csv, pd.DataFrame( columns=['nearestInd'],data = minDisInd)], axis=1)

        img = cv2.imread(img_dir + '/' + picID + '.jpg')
        img[:,:, 0] = 255
        img[:,:, 1] = 255
        img[:,:, 2] = 255
        ptList1 = dis1_csv.values
        for i in range(len(ptList1)):
            img = cv2.line(img, tuple(ptList1[i,:2]) , tuple(ptList2[ ptList1[i,2] ,:2]), line_color, 1)
        for i in range(len(ptList1)):
            img = cv2.circle(img, tuple(ptList1[i,:2]), 5, cellName_color[disName1], -1)
        for i in range(len(ptList2)):
            img = cv2.circle(img, tuple(ptList2[i,:2]), 5, cellName_color[disName2], -1)
        cv2.imwrite(picID + disName1 + '_' + disName2 + '_dis.png', img)

    elif disName1 == disName2:
        dis1_csv = pd.read_csv(csv_dir + '/' + picID + disName1 + '_pts.csv', usecols=['x_cord', 'y_cord'])
        ptList1 = dis1_csv.values[:,:2]

        minDisInd_list = []
        for i in range(len(ptList1)):
            currentPt = ptList1[i, :2]
            disMat = np.sqrt(np.sum(np.asarray(currentPt - ptList1)** 2, axis=1).astype(np.float32))
            minDisInd = np.argmin(disMat)
            disMat[minDisInd] = 1000.0
            minDisInd = np.argmin(disMat)
            minDisInd_list.append(minDisInd)

        minDisInd = np.array(minDisInd_list).reshape(-1,1)
        dis1_csv = pd.concat([dis1_csv, pd.DataFrame( columns=['nearestInd'],data = minDisInd)], axis=1)

        img = cv2.imread(img_dir + '/' + picID + '.jpg')
        img[:,:, 0] = 255
        img[:,:, 1] = 255
        img[:,:, 2] = 255
        ptList1 = dis1_csv.values
        for i in range(len(ptList1)):
            img = cv2.line(img, tuple(ptList1[i,:2]), tuple(ptList1[ptList1[i, 2],:2]), line_color, 1)
        for i in range(len(ptList1)):
            img = cv2.circle(img, tuple(ptList1[i,:2]), 5, cellName_color[disName1], -1)

        cv2.imwrite(picID + disName1 + '_dis.png', img)

def getAllPicsDisCSV():
    '''
    Get all distance data from the saved csv files (get from the above functions)
    '''
    base_dir = '/data/Datasets/MediImgExp'
    f = open( base_dir + '/' + 'AllDisData.csv','w',encoding='utf-8',newline="")
    csv_writer = csv.writer(f)
    csv_writer.writerow([ 'Ind','PosiTumourRatio','PosiLymphRatio',
        'DisTumourLymph','DisPosiTumour','DisPosiLymph',
        'DisPosiTumourPosiLymph','DisTumourPosiLymph'])

    process_dir = base_dir + '/process'
    csv_dir = base_dir + '/data_csv'
    pic_name = os.listdir(process_dir)
    picIDList = []
    for pic_name_ in pic_name:
        picIDList.append( pic_name_.split('_')[0] )
    for picID in picIDList:
        list_data = []
        list_data.append(picID)

        # PosiTumourRatio
        PosiTumourCsv = pd.read_csv( csv_dir+'/'+ picID +'positive_tumour_pts.csv')
        OtherTumourCsv = pd.read_csv( csv_dir+'/'+ picID +'other_tumour_pts.csv')
        Num_PosiTumour = PosiTumourCsv.shape[0]
        Num_OtherTumour = OtherTumourCsv.shape[0]
        if (Num_PosiTumour + Num_OtherTumour)!=0 :
            PosiTumourRatio =  Num_PosiTumour / (Num_PosiTumour + Num_OtherTumour)
        else:
            PosiTumourRatio = 'error'
        list_data.append(PosiTumourRatio)

        # PosiLymphRatio
        PosiLymphCsv = pd.read_csv( csv_dir+'/'+ picID +'positive_lymph_pts.csv')
        OtherLymphCsv = pd.read_csv( csv_dir+'/'+ picID +'other_lymph_pts.csv')
        Num_PosiLymph = PosiLymphCsv.shape[0]
        Num_OtherLymph = OtherLymphCsv.shape[0]
        if (Num_PosiLymph + Num_OtherLymph)!=0 :
            PosiLymphRatio =  Num_PosiLymph / (Num_PosiLymph + Num_OtherLymph)
        else:
            PosiLymphRatio = 'error'
        list_data.append(PosiLymphRatio)

        # DisTumourLymph
        ptList1_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_tumour_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList2_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_lymph_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList1 = ptList1_csv.values[:,:2]
        ptList2 = ptList2_csv.values[:,:2]
        ptList3_csv = pd.read_csv(csv_dir+'/'+ picID +'other_tumour_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList4_csv = pd.read_csv(csv_dir+'/'+ picID +'other_lymph_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList3 = ptList3_csv.values[:,:2]
        ptList4 = ptList4_csv.values[:,:2]
        ptList1 = np.concatenate((ptList1,ptList3), axis=0)
        ptList2 = np.concatenate((ptList2,ptList4), axis=0)
        DisTumourLymph = calcuNearestPtsDis(ptList1, ptList2)
        list_data.append(DisTumourLymph)

        # DisPosiTumour
        ptList1_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_tumour_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList1 = ptList1_csv.values[:,:2]
        DisPosiTumour = calcuNearestPtsDis2(ptList1)
        list_data.append(DisPosiTumour)

        # DisPosiLymph
        ptList1_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_lymph_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList1 = ptList1_csv.values[:,:2]
        DisPosiLymph = calcuNearestPtsDis2(ptList1)
        list_data.append(DisPosiLymph)

        # DisPosiTumourPosiLymph
        ptList1_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_tumour_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList2_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_lymph_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList1 = ptList1_csv.values[:,:2]
        ptList2 = ptList2_csv.values[:,:2]
        DisPosiTumourPosiLymph = calcuNearestPtsDis(ptList1, ptList2)
        list_data.append(DisPosiTumourPosiLymph)

        # DisTumourPosiLymph
        ptList1_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_tumour_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList2_csv = pd.read_csv(csv_dir+'/'+ picID +'positive_lymph_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList1 = ptList1_csv.values[:,:2]
        ptList2 = ptList2_csv.values[:,:2]
        ptList3_csv = pd.read_csv(csv_dir+'/'+ picID +'other_tumour_pts.csv',usecols=['x_cord', 'y_cord'])
        ptList3 = ptList3_csv.values[:,:2]
        ptList1 = np.concatenate((ptList1,ptList3), axis=0)
        DisTumourPosiLymph = calcuNearestPtsDis(ptList1, ptList2)
        list_data.append(DisTumourPosiLymph)

        csv_writer.writerow(list_data)

def adjustToMultiCSV():
    '''
    Divide the AllDisData.csv into  6+1=7 csv
    '''
    base_dir = '/data/Datasets/MediImgExp'
    alldata = pd.read_csv( base_dir + '/' + 'AllDisData.csv' )
    IndData = alldata['Ind'].values

    patient_Ind = []
    for IndName in IndData:
        patient_Ind.append(IndName.split('-')[0])
    patient_Ind = np.unique(patient_Ind)
    patient_Ind = sorted( list(map(int,patient_Ind)) )
    column_name = ['Ind','2D','3D','GAL9','LAG3','MHC','OX40','OX40L','PD1','PDL1','TIM3']

    # stage 1 calculate the 6 csv (10 cols for each csv)
    DisPosiTumour = pd.DataFrame(columns=column_name,index= patient_Ind)
    DisPosiTumour['Ind'] = patient_Ind
    patient_Id = patient_Ind
    column_names = ['2D','3D','GAL9','LAG3','MHC','OX40','OX40L','PD1','PDL1','TIM3']
    for patient in patient_Id:
        for column in column_names:
            combine_name = str(patient) + '-' + column
            exist_flag = (alldata['Ind'].str[0:len(combine_name)]== combine_name).any()
            if not exist_flag:
                continue
            valid_slice = alldata[ alldata['Ind'].str[0:len(combine_name)]== combine_name ]
            arr = valid_slice['DisTumourPosiLymph'].values
            if  arr.__contains__('error'):
                arr = np.setdiff1d(arr, ['error'])
                if not arr.shape[0]:
                    continue
            valid_slice_mean = np.mean( arr.astype(np.float32))
            DisPosiTumour.loc[ patient ,column ] = valid_slice_mean
    DisPosiTumour.to_csv( base_dir + '/' + 'DisTumourPosiLymph.csv',index=False )

    # stage 2 add the outputs (4 cols)
    all_data_name = base_dir + '/' + 'alldata2.csv'
    all_data = pd.read_csv(all_data_name)
    all_data.index = all_data['Ind']
    valid_columns = ['RELAPSE','RFS','DEATH','OS']
    valid_slice = all_data.loc[ patient_Ind, valid_columns ]
    DisPosiTumour = pd.read_csv( base_dir + '/' + 'PosiTumourRatio.csv',index_col=0)
    DisPosiTumour = pd.concat([DisPosiTumour,valid_slice],axis = 1)
    DisPosiTumour.to_csv( base_dir + '/' + 'PosiTumourRatio.csv' )

    # stage 3 calculate DisTumourLymph (use all markers' mean values)
    DisTumourLymph = pd.DataFrame(columns=['mean_10markers'],index= patient_Ind)
    patient_Id = patient_Ind
    column_names = [ 'mean_10markers']
    for patient in patient_Id:
        for column in column_names:
            combine_name = str(patient) + '-'
            exist_flag = (alldata['Ind'].str[0:len(combine_name)]== combine_name).any()
            if not exist_flag:
                continue
            valid_slice = alldata[ alldata['Ind'].str[0:len(combine_name)]== combine_name ]
            arr = valid_slice['DisTumourLymph'].values
            if  arr.__contains__('error'):
                arr = np.setdiff1d(arr, ['error'])
                if not arr.shape[0]:
                    continue
            valid_slice_mean = np.mean( arr.astype(np.float32))
            DisTumourLymph.loc[ patient ,column ] = valid_slice_mean

    DisTumourLymph.to_csv( base_dir + '/' + 'DisTumourLymph.csv' )
    all_data_name = base_dir + '/' + 'alldata2.csv'
    all_data = pd.read_csv(all_data_name)
    all_data.index = all_data['Ind']
    valid_columns = ['RELAPSE','RFS','DEATH','OS']
    valid_slice = all_data.loc[ patient_Ind, valid_columns]
    DisTumourLymph = pd.concat([DisTumourLymph,valid_slice],axis = 1)
    DisTumourLymph.to_csv( base_dir + '/' + 'DisTumourLymph.csv')

def getAllFeatureCSV():
    base_dir = '/data/Datasets/MediImgExp/csv'
    alldata = pd.read_csv( base_dir + '/' + 'AllDisData.csv' )
    oridata = pd.read_csv( base_dir + '/' + 'alldata2.csv',index_col=0 )
    ori_columns = oridata.columns.values
    ori_columns = ori_columns[4:-4] # original 40 feature names
    csv_name = [ 'DisPosiLymph','DisPosiTumour', 'DisPosiTumourPosiLymph',
        'DisTumourPosiLymph','PosiLymphRatio','PosiTumourRatio', 'DisTumourLymph']

    meaningful_feature_OS = { 'DisPosiLymph_PD1':24,'DisPosiLymph_OX40L':48.14,'DisPosiLymph_OX40':98.33,'DisPosiLymph_3D':13.44,
        'DisPosiTumour_TIM3':546.86, 'DisPosiTumour_GAL9':85.97,'DisPosiTumour_2D':24.22,
        'DisPosiTumourPosiLymph_3D':20.88, 'DisPosiTumourPosiLymph_MHC':18.68,'DisPosiTumourPosiLymph_OX40L':40.02,
        'DisTumourPosiLymph_OX40L':173.56, 'DisTumourPosiLymph_2D':223.71,'DisTumourPosiLymph_3D':21.19,'DisTumourPosiLymph_OX40':445.89,
        'PosiLymphRatio_3D':0.97, 'PosiLymphRatio_OX40':0.44,'PosiLymphRatio_GAL9':0.23,'PosiLymphRatio_2D':0.37,
        'PosiTumourRatio_MHC':0.93,'PosiTumourRatio_GAL9':0.41,'PosiTumourRatio_3D':1.0 }

    IndData = alldata['Ind'].values
    patient_Ind = []
    for IndName in IndData:
        patient_Ind.append(IndName.split('-')[0])
    patient_Ind = np.unique(patient_Ind)
    patient_Ind = sorted( list(map(int,patient_Ind)) )

    # create a super_csv including all features 61+40=101 features / 4 outputs
    super_csv = pd.DataFrame(index= patient_Ind)
    super_csv_copy = pd.DataFrame(index= patient_Ind)

    super_csv['Ind'] = patient_Ind
    for column in ori_columns: # 40 features
        super_csv[column] = oridata.loc[ patient_Ind,column]
    for i in range(0,7):
        csvName = os.path.join(base_dir,csv_name[i] + '.csv')
        csvData = pd.read_csv(csvName,index_col=0)
        if i == 6: # 1 feature
            column_name = csvData.columns.values[:1]
            super_csv[ csv_name[i]+'_'+column_name[0] ] = csvData.loc[ patient_Ind, column_name[0] ]
            output_name = csvData.columns.values[-4:]
            for k in range(0,4):
                super_csv[ output_name[k] ] = csvData.loc[ patient_Ind,output_name[k] ]
            break
        column_name = csvData.columns.values[:10]
        for j in range(0,10): # 6*10 = 60 features
            super_csv[ csv_name[i]+'_'+column_name[j] ] = csvData.loc[ patient_Ind, column_name[j] ]

    for key, values in meaningful_feature_OS.items():
        super_csv_copy[key] = super_csv.loc[patient_Ind, key ]
        super_csv_copy.loc[ super_csv_copy[key] <= values, key] = 0
        super_csv_copy.loc[ super_csv_copy[key] > values, key ] = 1

    for k in range(0,4):
        super_csv_copy[ output_name[k] ] = csvData.loc[ patient_Ind,output_name[k] ]
    super_csv_copy = super_csv_copy.dropna(axis=0,how='any')

    #super_csv.to_csv(base_dir+'/'+ 'super1.csv',index =False)
    #super_csv.to_csv(base_dir+'/'+ 'super2.csv',index =False)
    #super_csv_copy.to_csv(base_dir+'/'+ 'super3.csv',index =True)
    #super_csv_copy.to_csv(base_dir+'/'+ 'super4.csv',index =True)

def getAllFeatureCSV2():
    base_dir = '/data/Datasets/MediImgExp/csv'
    alldata = pd.read_csv( base_dir + '/' + 'AllDisData.csv' )
    oridata = pd.read_csv( base_dir + '/' + 'alldata2.csv',index_col=0 )
    ori_columns = oridata.columns.values
    ori_columns = np.concatenate(( ori_columns[19:32],ori_columns[44:48] ) ,axis=0 ) # original 40 feature names
    csv_name = [ 'DisPosiLymph','DisPosiTumour', 'DisPosiTumourPosiLymph',
        'DisTumourPosiLymph','PosiLymphRatio','PosiTumourRatio', 'DisTumourLymph']

    meaningful_feature_OS = { 'DisPosiLymph_OX40L':61.55,
        'DisPosiTumour_GAL9':119.40,
        'DisPosiTumourPosiLymph_3D':20.28, 'DisPosiTumourPosiLymph_OX40L':40.02,
        'DisTumourPosiLymph_3D':21.19, 'DisTumourPosiLymph_MHC':135.56,
        'PosiLymphRatio_3D':0.97, 'PosiLymphRatio_OX40':0.52,'PosiLymphRatio_GAL9':0.2 }

    IndData = alldata['Ind'].values
    patient_Ind = []
    for IndName in IndData:
        patient_Ind.append(IndName.split('-')[0])
    patient_Ind = np.unique(patient_Ind)
    patient_Ind = sorted( list(map(int,patient_Ind)) )

    # create a super_csv including all features 61+xx features / 4 outputs
    super_csv = pd.DataFrame(index= patient_Ind)
    super_csv_copy = pd.DataFrame(index= patient_Ind)

    super_csv['Ind'] = patient_Ind
    for column in ori_columns: # 40 features
        super_csv[column] = oridata.loc[ patient_Ind,column]
    for i in range(0,7):
        csvName = os.path.join(base_dir,csv_name[i] + '.csv')
        csvData = pd.read_csv(csvName,index_col=0)
        if i == 6: # 1 feature
            column_name = csvData.columns.values[:1]
            super_csv[ csv_name[i]+'_'+column_name[0] ] = csvData.loc[ patient_Ind, column_name[0] ]
            output_name = csvData.columns.values[-4:]
            for k in range(0,4):
                super_csv[ output_name[k] ] = csvData.loc[ patient_Ind,output_name[k] ]
            break
        column_name = csvData.columns.values[:10]
        for j in range(0,10): # 6*10 = 60 features
            super_csv[ csv_name[i]+'_'+column_name[j] ] = csvData.loc[ patient_Ind, column_name[j] ]

    for key, values in meaningful_feature_OS.items():
        super_csv_copy[key] = super_csv.loc[patient_Ind, key ]
        super_csv_copy.loc[ super_csv_copy[key] <= values, key] = 0
        super_csv_copy.loc[ super_csv_copy[key] > values, key ] = 1

    for k in range(0,4):
        super_csv_copy[ output_name[k] ] = csvData.loc[ patient_Ind,output_name[k] ]
    #super_csv_copy = super_csv_copy.dropna(axis=0,how='any')
    #super_csv.to_csv(base_dir+'/rfs_csv'+'/'+ 'super1.csv',index =False)
    #super_csv.to_csv(base_dir+'/'+ 'super2.csv',index =False)
    super_csv_copy.to_csv(base_dir+'/rfs_csv'+'/'+ 'super3.csv',index =True)
    #super_csv_copy.to_csv(base_dir+'/'+ 'super4.csv',index =True)

if __name__ == '__main__':

    getAllPicsDisCSV()
    adjustToMultiCSV()
    getAllFeatureCSV2()

    # draw the distance pics
    drawDistancePic(disName1='all_tumour', disName2='positive_lymph', picID='0-GAL9-1')

max_stars_repo_name: dddtqshmpmz/IHC

max_issues_repo_name: dddtqshmpmz/IHC

max_forks_repo_name: dddtqshmpmz/IHC