text:
import os
import shutil
import numpy as np
import multiprocessing

from aria_component import aria_component
import sync_times

import time
import logging

def run_point(network, idx, qdpt, workdir):
    return_code = network.execute_quadrature_point(idx, qdpt, workdir)
    assert(return_code == 0)

class network():

    def __init__(self, component_list, workdir=".", restart=False, merge_timelines=False):
        self.components = component_list
        self.working_directory = os.path.abspath(workdir)
        self.linked_endo_ports = []
        self.linked_exo_ports = []
        self.pce_dim = self.get_num_exo_ports()
        if restart:
            for comp in self.components:
                if comp.pc_model is None:
                    comp.generate_pc_model(self.pce_dim)
            print("Collect Data")
            self.endo_pce_coeffs, self.QoI_pce_coeffs = self.collect_data()
        else:
            print("Initialize PCE")
            self.endo_pce_coeffs = []
            self.QoI_pce_coeffs = []
            self.initialize_PCE()
            print("Create subfolders")
            self.create_folders()
        self.historyfile = None
        self.merge_timelines = merge_timelines

    def get_num_exo_ports(self):
        """
        Calculate total number of exogenous ports in network
        """
        count = 0
        for comp in self.components:
            count += comp.get_num_exogenous_ports()
        return count - len(self.linked_exo_ports)

    def get_num_endo_ports(self):
        """
        Calculate total number of endogenous ports in network
        """
        count = 0
        for comp in self.components:
            count += comp.get_num_endogenous_ports()
        return count

    def get_num_endo_pc_coefs(self):
        """
        Count total number of PCE coefficients in network
        """
        count = 0
        for comp in self.components:
            count += comp.get_num_pc_coeffs()
        return count

    def get_num_QoIs(self):
        """
        Count total number of QoIs in network
        """
        count = 0
        for comp in self.components:
            count += comp.get_num_QoI()
        return count

    def link_endogenous_ports(self, componentA, portA, componentB, portB):
        """
        Link endogenous ports between 2 components
        Useful for only evaluating subset of output BCs on components
        """
        self.linked_endo_ports.append([componentA, portA, componentB, portB])

    def link_exogenous_ports(self, componentA, portA, componentB, portB):
        """
        Link exogenous ports between 2 components
        Useful for shared material properties/BCs between components
        """
        if componentA > componentB:
            self.linked_exo_ports.append([componentB, portB, componentA, portA])
        else:
            self.linked_exo_ports.append([componentA, portA, componentB, portB])
        self.pce_dim -= 1

    def execute_quadrature_point(self, current_comp_idx, qdpt, workdir):
        """
        Execute forward problem for single component, quadrature point
        Evaluates endogenous outputs for other components and exogenous inputs
        """
        launchdir = os.getcwd()
        os.chdir(workdir)

        logging.basicConfig(filename='network.log', level=logging.DEBUG)

        start=time.time()
        # Evaluate PCEs for other components
        if len(self.linked_endo_ports) == 0:
            for idx,othercomp in enumerate(self.components):
                logging.info('Evaluating PCEs for component %d', idx)
                if idx!=current_comp_idx:
                    othercomp.set_all_endogenous_values(self.endo_pce_coeffs[idx], qdpt)
        else:
            for link in self.linked_endo_ports:
                compA = link[0]
                compB = link[2]
                if compA == current_comp_idx:
                    self.components[compB].set_port_endogenous_values(self.endo_pce_coeffs[compB], qdpt, link[3])
                elif compB == current_comp_idx:
                    self.components[compA].set_port_endogenous_values(self.endo_pce_coeffs[compA], qdpt, link[1])
        end=time.time()
        logging.info('Evaluating PCEs for other components took %f s.' % (end-start))

        start=time.time()
        # Convert Quadrature point in \xi_i to equivalent samples of input parameters
        # (taking advantage of the fact that inputs are assumed to be Gaussian)
        # This is equivalent to evaluating 1st order PC expansions for the input parameters.
        exogenous_idx = 0
        for compidx, comp in enumerate(self.components):
            for portidx,port in enumerate(comp.get_exogenous_ports()):
                found = False
                for link in self.linked_exo_ports:
                    if compidx == link[2] and portidx == link[3]:
                        port.value = self.components[link[0]].exogenous_ports[link[1]].value
                        found = True
                if not found:
                    port.value = port.mean + port.std * qdpt[0,exogenous_idx]
                    exogenous_idx += 1
        end=time.time()
        logging.info('Evaluating PCEs for exogenous ports took %f s.' % (end-start))

        start=time.time()
        # Solve forward problem on component
        running_sim = self.components[current_comp_idx].execute()
        end=time.time()
        logging.info('Starting simulation took %f s.' % (end-start))

        start=time.time()
        # Wait for simulation to complete
        running_sim.communicate()
        end=time.time()
        logging.info('Simulation execution took %f s.' % (end-start))

        # Check return code
        if not running_sim.returncode == 0:
            print("simulation failed in %s" % workdir)

        os.chdir(launchdir)
        return running_sim.returncode

    def create_folders(self):
        """
        Creates directory structure for working directory
        allows for simultaneous execution across components/quadrature points
        """
        # TODO: utilize symbolic links where possible to reduce overhead

        # Note that this will clear everything from working_directory.  Be careful to avoid data loss
        if os.path.isdir(self.working_directory):
            shutil.rmtree(self.working_directory)
        os.makedirs(self.working_directory)

        for compidx,comp in enumerate(self.components):
            compfolder=os.path.join(self.working_directory, "component_%d" % compidx)
            os.makedirs(compfolder)
            qdpts, totquat = comp.get_quadpts(self.pce_dim)
            for quad_iter in range(totquat):
                # Create subdirectory for each quadrature point
                #   and copy required files into it
                newfolder=os.path.join(compfolder, "qdpt_%d" % quad_iter)
                if not os.path.isdir(newfolder):
                    os.makedirs(newfolder)
                    for c in self.components:
                        for f in c.get_required_files():
                            shutil.copy2(f,newfolder)

    def consolidate_timesteps(self):
        """
        Navigate directory structure and generate composite timeline for data
        """
        for myidx,comp in enumerate(self.components):
            qdpts, totquat = comp.get_quadpts(self.pce_dim)
            folder_list = []

            for quad_iter in range(totquat):
                targetfolder=os.path.join(self.working_directory, "component_%d" % myidx)
                targetfolder=os.path.join(targetfolder, "qdpt_%d" % quad_iter)
                folder_list.append(targetfolder)

            filename = comp.get_output_filename()
            file_list = []
            for folder in folder_list:
                file_list.append(folder + "/" + filename)

            times = sync_times.get_combined_timeline(file_list)
            comp.set_solution_times(times)

            if self.historyfile:
                os.remove(self.historyfile)

    def collect_data(self):
        """
        Navigate directory structure and extract nodal-temporal values
        Perform projection to generate nodal-temporal PCE coefficients
        """
        launchdir = os.getcwd()

        endo_pce_coeffs = []
        QoI_pce_coeffs = []
        for myidx,comp in enumerate(self.components):
            qdpts, totquat = comp.get_quadpts(self.pce_dim)

            endo_data=np.zeros(( totquat, comp.get_num_endogenous_ports() ))
            QoI_data=np.zeros(( totquat, comp.get_num_QoI()*comp.get_num_timesteps() ))

            endo_pce_coeffs.append(np.zeros(( comp.get_num_endogenous_ports(), comp.get_num_pc_terms() )))
            QoI_pce_coeffs.append(np.zeros(( comp.get_num_QoI()*comp.get_num_timesteps(), comp.get_num_pc_terms() )))

            for quad_iter in range(totquat):
                targetfolder=os.path.join(self.working_directory, "component_%d" % myidx)
                targetfolder=os.path.join(targetfolder, "qdpt_%d" % quad_iter)
                os.chdir(targetfolder)

                logging.basicConfig(filename='network.log', level=logging.DEBUG)
                logging.info('Component %d, QuadPt %d' % (myidx, quad_iter))
                logging.info('endo_data.size: %d,%d' % (endo_data.shape[0], endo_data.shape[1]))
                logging.info('QoI_data.size: %d,%d' % (QoI_data.shape[0], QoI_data.shape[1]))

                # Collect endogenous data
                logging.info('Collect Endogenous Data')
                endo_at_qdpt = comp.get_endogenous_data()
                index = 0
                for nodal_data in endo_at_qdpt:
                    for nid in nodal_data:
                        endo_data[quad_iter, index] = nodal_data[nid]
                        index += 1

                # Collect QoI data
                QoI_at_qdpt = comp.get_QoI_data()

                for i,QoI in enumerate(comp.get_QoIs()):
                    QoI_data[[quad_iter],(i*comp.get_num_timesteps()):((i+1)*comp.get_num_timesteps())] = QoI_at_qdpt[QoI]

            # Perform PC projection via quadrature integration
            endo_ck = np.zeros(( endo_data.shape[1], comp.get_num_pc_terms() ))
            for i in range(endo_data.shape[1]):
                endo_ck[i,...] = comp.GalerkinProjection(endo_data[:,i])
            endo_pce_coeffs[myidx] = endo_ck

            QoI_ck = np.zeros(( QoI_data.shape[1], comp.get_num_pc_terms() ))
            for i in range(QoI_data.shape[1]):
                QoI_ck[i,...] = comp.GalerkinProjection(QoI_data[:,i])
            QoI_pce_coeffs[myidx] = QoI_ck

        os.chdir(launchdir)
        return endo_pce_coeffs, QoI_pce_coeffs

    def save_all_data(self, variable_name, comps=None):
        """
        Output nodal data to exodus
        """
        launchdir = os.getcwd()

        if comps is None:
            comps = range(len(self.components))

        for compidx in comps:
            exodus_filename = self.components[compidx].get_output_filename()
            pce_filename = "%s_pce.e" % exodus_filename.split('.')[:-1][0]

            qdpts, totquat = self.components[compidx].get_quadpts(self.pce_dim)
            all_data = None
            for quad_iter in range(totquat):
                targetfolder=os.path.join(self.working_directory, "component_%d" % compidx)
                targetfolder=os.path.join(targetfolder, "qdpt_%d" % quad_iter)
                os.chdir(targetfolder)

                # Collect nodal data
                data = self.components[compidx].get_all_data(variable_name, exodus_filename)

                if quad_iter == 0:
                    # size numpy array to hold data
                    datasize=0
                    for timestep, nodal_data in enumerate(data):
                        datasize += len(nodal_data)
                    all_data = np.zeros(( totquat, datasize ))

                for timestep, nodal_data in enumerate(data):
                    numnodes = len(nodal_data)
                    for i,node_val in enumerate(nodal_data):
                        index = timestep*numnodes + i
                        all_data[quad_iter, index] = node_val

            # Perform PC projection via quadrature integration
            nodal_ck = np.zeros(( all_data.shape[1], self.components[compidx].get_num_pc_terms() ))
            for i in range(all_data.shape[1]):
                nodal_ck[i,...] = self.components[compidx].GalerkinProjection(all_data[:,i])

            # Save nodal PCE
            os.chdir(launchdir)
            self.components[compidx].save_nodal_pce(nodal_ck, exodus_filename, pce_filename)

    def initialize_PCE(self):
        """
        Initialize all PCE coefficients
        """
        self.endo_pce_coeffs = []
        for comp in self.components:
            if comp.get_pc_model() is None:
                comp.generate_pc_model(self.pce_dim)

            my_endo_pce_coeffs = comp.initialize_PCE()

            self.endo_pce_coeffs.append(my_endo_pce_coeffs)
            self.QoI_pce_coeffs.append(np.zeros(( comp.get_num_QoI()*comp.get_num_timesteps(), comp.get_num_pc_terms() )))

    def compute_relative_residual(self, endo_pce_coeffs_new):
        """
        Compute relative stepsize ||x_k - x_(k-1)|| / ||x_k|| where ||.|| is the infinity norm
        """
        delta = []
        norm_factor = []
        for myidx in range(len(self.endo_pce_coeffs)):
            delta.append(np.linalg.norm(endo_pce_coeffs_new[myidx] - self.endo_pce_coeffs[myidx], np.inf))
            norm_factor.append(np.linalg.norm(endo_pce_coeffs_new[myidx], np.inf))

        return max(delta)/max(norm_factor)

    def anderson_acceleration(self, endo_pce_coeffs_new, mk=5):
        """
        Perform anderson accleration to modify step taken in DDUQ iteration
        """
        if self.historyfile is None:
            self.historyfile = os.path.join(self.working_directory, "PCE_history.csv")

        # serialize data
        compstrings = []
        for comp_coeffs in self.endo_pce_coeffs:
            flat = comp_coeffs.flatten()
            compstrings.append(','.join(['%f' % num for num in flat]))
        line = ','.join(compstrings)
        line += '\n'
        compstrings = []
        for comp_coeffs in endo_pce_coeffs_new:
            flat = comp_coeffs.flatten()
            compstrings.append(','.join(['%f' % num for num in flat]))
        line += ','.join(compstrings) + '\n'

        # Append to history file
        with open(self.historyfile, "a+") as f:
            f.write(line)

        # read in data
        data = np.genfromtxt(self.historyfile, delimiter=',', dtype=float)
        xi = data[::2].transpose()

        # Perform acceleration step
        if xi.shape[1] > 1:
            gi = data[1::2].transpose()
            fi = gi - xi
            dfi = np.diff(fi)
            dgi = np.diff(gi)

            # Solve linear least squares problem
            # TODO: use QR factorization, save off Q,R and update with scipy.linalg.qr_insert
            gamma = np.linalg.lstsq(dfi, fi[:,-1], rcond=None)[0]

            # Take step
            next_x = gi[:,-1] - np.matmul(dgi, gamma)

            # Unpack step
            endo_pce_coeffs_accel = []
            offset = 0
            for comp_coeffs in self.endo_pce_coeffs:
                myshape = comp_coeffs.shape
                mydata = next_x[offset:offset+np.product(myshape)].reshape(myshape)
                endo_pce_coeffs_accel.append(mydata)
                offset += np.product(myshape)

            # Delete first 2 lines from file if too many steps are recorded
            if xi.shape[1] > mk:
                os.system("tail -n +3 %s >> tmp.csv" % self.historyfile)
                os.system("mv tmp.csv %s" % self.historyfile)

            return endo_pce_coeffs_accel
        else:
            return endo_pce_coeffs_new

    def run_simulation(self, max_iter=30, tolerance=0.001):
        """
        Perform DDUQ simulation
        """

        # Perform iterations
        for jacobi_iter in range(max_iter):
            print("Performing iteration %d..." % jacobi_iter)
            args=[]
            for myidx,comp in enumerate(self.components):
                qdpts, totquat = comp.get_quadpts(self.pce_dim)
                for quad_iter in range(totquat):
                    targetfolder=os.path.join(self.working_directory, "component_%d" % myidx)
                    targetfolder=os.path.join(targetfolder, "qdpt_%d" % quad_iter)
                    args.append((self, myidx, qdpts[[quad_iter],:], targetfolder))

            # Begin simulations
            print("Starting aria jobs")
            tasks=[]
            for myargs in args:
                tasks.append(multiprocessing.Process(target=run_point, args=myargs))
                tasks[-1].start()

            # Wait for all simulations to complete
            print("Waiting for aria jobs to complete...")
            unfinished = len(tasks)
            print("%d/%d tasks to go..." % (unfinished,unfinished))
            while True:
                incomplete = [task.is_alive() for task in tasks]
                if incomplete.count(True) < unfinished:
                    unfinished = incomplete.count(True)
                    print("%d/%d tasks to go..." % (unfinished,len(incomplete)))
                if unfinished > 0:
                    time.sleep(2)
                else:
                    break

            for task in tasks:
                task.join()
                assert(task.exitcode==0)

            if jacobi_iter==0 and self.merge_timelines:
                self.consolidate_timesteps()

            # Collect update PCE coefficients
            print("Collect Data")
            endo_pce_coeffs_new, self.QoI_pce_coeffs = self.collect_data()

            if len(self.components)>1:
                print("Perform acceleration step")
                # Perform acceleration step
                if jacobi_iter==0:
                    next_step = endo_pce_coeffs_new
                    self.endo_pce_coeffs = endo_pce_coeffs_new
                else:
                    next_step = self.anderson_acceleration(endo_pce_coeffs_new)

                    print("Check for convergence")
                    # Check for convergence
                    stepsize = self.compute_relative_residual(endo_pce_coeffs_new)
                    self.endo_pce_coeffs = next_step
                    print(stepsize)
                    if stepsize < tolerance:
                        break

        return self.QoI_pce_coeffs, self.endo_pce_coeffs

max_stars_repo_name: jtencer/PyNetUQ

max_issues_repo_name: jtencer/PyNetUQ

max_forks_repo_name: jtencer/PyNetUQ


text:
import numpy as np
from matplotlib import pyplot as plt

import chainer
import chainer.functions as F
import chainer.links as L
import chainercv

"""
forked from https://github.com/leetenki/YOLOtiny_v2_chainer
"""

class YOLOv2_tiny(chainer.Chain):

    def __init__(self, n_classes=20, n_boxes=5):
        super(YOLOv2_tiny, self).__init__()
        with self.init_scope():
            self.conv1 = L.Convolution2D(None, 16, 3, stride=1, pad=1, nobias=True)
            self.bn1 = L.BatchNormalization(16)
            self.conv2 = L.Convolution2D(None, 32, 3, stride=1, pad=1, nobias=True)
            self.bn2 = L.BatchNormalization(32)
            self.conv3 = L.Convolution2D(None, 64, 3, stride=1, pad=1, nobias=True)
            self.bn3 = L.BatchNormalization(64)
            self.conv4 = L.Convolution2D(None, 128, 3, stride=1, pad=1, nobias=True)
            self.bn4 = L.BatchNormalization(128)
            self.conv5 = L.Convolution2D(None, 256, 3, stride=1, pad=1, nobias=True)
            self.bn5 = L.BatchNormalization(256)
            self.conv6 = L.Convolution2D(None, 512, 3, stride=1, pad=1, nobias=True)
            self.bn6 = L.BatchNormalization(512)
            self.conv7 = L.Convolution2D(None, 1024, 3, stride=1, pad=1, nobias=True)
            self.bn7 = L.BatchNormalization(1024)
            self.conv8 = L.Convolution2D(None, 1024, 3, stride=1, pad=1, nobias=True)
            self.bn8 = L.BatchNormalization(1024)
            self.conv9 = L.Convolution2D(None, 125, 1, stride=1, pad=0, nobias=True)

    def __call__(self, x):
        h = F.leaky_relu(self.bn1(self.conv1(x)),slope=0.1)
        h = F.max_pooling_2d(h, ksize=2, stride=2, pad=0)
        h = F.leaky_relu(self.bn2(self.conv2(h)),slope=0.1)
        h = F.max_pooling_2d(h, ksize=2, stride=2, pad=0)
        h = F.leaky_relu(self.bn3(self.conv3(h)),slope=0.1)
        h = F.max_pooling_2d(h, ksize=2, stride=2, pad=0)
        h = F.leaky_relu(self.bn4(self.conv4(h)),slope=0.1)
        h = F.max_pooling_2d(h, ksize=2, stride=2, pad=0)
        h = F.leaky_relu(self.bn5(self.conv5(h)),slope=0.1)
        h = F.max_pooling_2d(h, ksize=2, stride=2, pad=0)
        h = F.leaky_relu(self.bn6(self.conv6(h)),slope=0.1)
        h = F.leaky_relu(self.bn7(self.conv7(h)),slope=0.1)
        h = F.leaky_relu(self.bn8(self.conv8(h)),slope=0.1)
        h = self.conv9(h)
        return h

def load_npz(filename, model):
    with np.load(filename) as f:
        weights = dict(f.items())

    for i in range(1, 9):
        exec('model.conv{0}.W.data = weights["c{0}/c/W"]'.format(i))
        exec('model.bn{0}.avg_mean = weights["c{0}/n/avg_mean"]'.format(i))
        exec('model.bn{0}.avg_var = weights["c{0}/n/avg_var"]'.format(i))
        exec('model.bn{0}.N = weights["c{0}/n/N"]'.format(i))
        exec('model.bn{0}.gamma.data = weights["c{0}/n/gamma"]'.format(i))
        exec('model.bn{0}.beta.data = weights["c{0}/b/b"]'.format(i))
    i = 9
    exec('model.conv{0}.W.data = weights["c{0}/c/W"]'.format(i))

def _sigmoid(x):
    return 1.0 / (np.exp(-x) + 1.0)

def _softmax(x):
    x = np.array([x]) # reshape (len(x),) to (1, len(x))
    return F.softmax(x).data

def _get_detected_boxes(ans, n_grid_x, n_grid_y, n_bbox, n_classes, prob_thresh, img_width, img_height, biases):
    detected_boxes = []
    grid_width = img_width / float(n_grid_x)
    grid_height = img_height / float(n_grid_y)

    boxes = []
    scores = []
    labels = []
    for grid_y in range(n_grid_y):
        for grid_x in range(n_grid_x):
            for i in range(n_bbox):
                box = ans[grid_y, grid_x, i, 0:4] # (4,)
                conf = _sigmoid(ans[grid_y, grid_x, i, 4])
                probs = _softmax(ans[grid_y, grid_x, i, 5:])[0] # (20,)

                p_class = probs * conf # (20,)
                if np.max(p_class) < prob_thresh:
                    continue

                class_id = np.argmax(p_class)
                x = (grid_x + _sigmoid(box[0])) * grid_width
                y = (grid_y + _sigmoid(box[1])) * grid_height
                w = np.exp(box[2]) * biases[i][0] * grid_width
                h = np.exp(box[3]) * biases[i][1] * grid_height

                b = [y - 0.5*h, x - 0.5*w, y + 0.5*h, x+0.5*w]

                boxes.append(b)
                labels.append(class_id)
                scores.append(max(p_class))

    return boxes, labels, scores

def postprocess(ans, im_org_):
    img_width, img_height = 352,352
    org_img_height, org_img_width = im_org_.shape[0:2]
    n_grid_x = img_width//32
    n_grid_y = img_height//32
    n_classes = 20
    n_bbox = 5
    biases = [[1.08, 1.19], [3.42, 4.41], [6.63, 11.38], [9.42, 5.11], [16.62, 10.52]]
    prob_thresh = 0.2
    iou_thresh = 0.04
    #org_img_height, org_img_width = im_org.shape[0:2]
    org_img_height, org_img_width = 352, 352
    classes = ["aeroplane", "bicycle", "bird", "boat", "bottle",
              "bus", "car", "cat", "chair", "cow",
              "diningtable", "dog", "horse", "motorbike", "person",
              "pottedplant", "sheep", "sofa", "train","tvmonitor"]

    ans = ans.transpose(1, 2, 0) # (13, 13, 125)
    ans = ans.reshape(n_grid_y, n_grid_x, n_bbox, (n_classes + 5)) # (13, 13, 5, 25)

    # compute detected boxes
    bboxes, labels, scores = _get_detected_boxes(ans, n_grid_x, n_grid_y, n_bbox, n_classes, prob_thresh, org_img_width, org_img_height, biases)

    bboxes = np.asarray(bboxes, dtype=np.int32)
    labels = np.asarray(labels, dtype=np.int32)
    scores = np.asarray(scores, dtype=np.float32)

    #for box, label, score in zip(bboxes, labels, scores):
    #    print(classes[label], score, box)
    #print()
    idx = chainercv.utils.non_maximum_suppression(bboxes, iou_thresh, score=scores, limit=None)
    bboxes, labels, scores = bboxes[idx], labels[idx], scores[idx]

    for box, label, score in zip(bboxes, labels, scores):
        print(classes[label], score, box)

    chainercv.visualizations.vis_bbox(im_org_, bboxes, label=labels, score=scores, label_names=classes, ax=None)

ext: py

reason: import numpy


text:
# A test configuration dict

from paltas.Sampling import distributions
import numpy as np
from scipy.stats import uniform, norm, loguniform, lognorm, multivariate_normal
from paltas.Substructure.los_dg19 import LOSDG19
from paltas.Substructure.subhalos_dg19 import SubhalosDG19
from paltas.MainDeflector.simple_deflectors import PEMDShear
from paltas.Sources.cosmos import COSMOSCatalog
from scipy.stats import truncnorm
from paltas.Sources.sersic import SingleSersicSource
from paltas.PointSource.single_point_source import SinglePointSource

# Define a multivariate distribution we'll use
mean = np.ones(2)
cov = np.array([[1.0,0.7],[0.7,1.0]])
min_values = np.zeros(2)
tmn = distributions.TruncatedMultivariateNormal(mean,cov,min_values,None)

# Define the numerics kwargs
kwargs_numerics = {'supersampling_factor':1}
numpix = 64
seed = 10

# Define some general image kwargs for the dataset
mask_radius = 0.2
no_noise = True

# Define arguments that will be used multiple times
output_ab_zeropoint = 25.127

# Define the cosmos path
cosmos_folder = './test_data/cosmos/'

config_dict = {
	'subhalo':{
		'class': SubhalosDG19,
		'parameters':{
			'sigma_sub':uniform(loc=0,scale=5e-4).rvs,
			'shmf_plaw_index':-1.83,
			'm_pivot': 1e8,'m_min': 1e9,'m_max': 1e10,
			'c_0':18,'conc_zeta':-0.2,'conc_beta':0.8,
			'conc_m_ref': 1e8,'dex_scatter': 0.1, 'k1':0.88,
			'k2':1.7
		}
	},
	'los':{
		'class': LOSDG19,
		'parameters':{
			'm_min':1e8,'m_max':1e10,'z_min':0.01,
			'dz':0.01,'cone_angle':8.0,'r_min':0.5,'r_max':10.0,
			'c_0':18,'conc_zeta':-0.2,'conc_beta':0.8,'conc_m_ref': 1e8,
			'dex_scatter': 0.1,'delta_los':uniform(loc=0,scale=5e-3).rvs,
			'alpha_dz_factor':20.0
		}
	},
	'main_deflector':{
		'class': PEMDShear,
		'parameters':{
			'M200': loguniform(a=1e11,b=5e13).rvs,
			'z_lens': 0.5,
			'gamma': lognorm(scale=2.01,s=0.1).rvs,
			'theta_E': lognorm(scale=1.1,s=0.05).rvs,
			'e1,e2': multivariate_normal(np.zeros(2),
				np.array([[0.03,0.015],[0.015,0.03]])).rvs,
			'center_x': norm(loc=0.0,scale=0.16).rvs,
			'center_y': norm(loc=0.0,scale=0.16).rvs,
			'gamma1': norm(loc=0.0,scale=0.05).rvs,
			'gamma2': norm(loc=0.0,scale=0.05).rvs,
			'ra_0':0.0, 'dec_0':0.0
		}
	},
	'source':{
		'class': COSMOSCatalog,
		'parameters':{
			'z_source':1.5,'cosmos_folder':cosmos_folder,
			'max_z':None,'minimum_size_in_pixels':None,
			'faintest_apparent_mag':None,
			'smoothing_sigma':0.0,'random_rotation':True,
			'center_x':norm(loc=0.0,scale=0.16).rvs,
			'center_y':norm(loc=0.0,scale=0.16).rvs,
			'min_flux_radius':None,'output_ab_zeropoint':output_ab_zeropoint}
	},
	'lens_light':{
		'class': SingleSersicSource,
		'parameters':{
			'z_source':1.5,
			'magnitude':uniform(loc=19,scale=2).rvs,
			'output_ab_zeropoint':output_ab_zeropoint,
			'R_sersic':truncnorm(-1.333,np.inf,loc=0.8,scale=0.15).rvs,
			'n_sersic':truncnorm(-2.,np.inf,loc=3,scale=0.5).rvs,
			'e1,e2':distributions.EllipticitiesTranslation(
				q_dist=truncnorm(-np.inf,1.,loc=0.85,scale=0.15).rvs,
				phi_dist=uniform(loc=-np.pi/2,scale=np.pi).rvs),
			'center_x':0.0,
			'center_y':0.0}
	},
	'point_source':{
		'class': SinglePointSource,
		'parameters':{
			'x_point_source':0.001,
			'y_point_source':0.001,
			'magnitude':uniform(loc=20,scale=2.5).rvs,
			'output_ab_zeropoint':25.127,
			'mag_pert': distributions.MultipleValues(dist=norm(1,0.1).rvs,num=5),
			'compute_time_delays':True,
			'kappa_ext': distributions.KappaTransformDistribution(
				n_dist=norm(loc=1.,scale=0.025).rvs),
			'time_delay_error': distributions.MultipleValues(
				dist=norm(loc=0.,scale=0.25).rvs, num=5)
		}
	},
	'lens_equation_solver':{
		'parameters':{
			'min_distance':0.05,
			'search_window':numpix*0.08,
			'num_iter_max':100,
			'precision_limit':10**(-10)
		}
	},
	'cosmology':{
		'parameters':{
			'cosmology_name': 'planck18'
		}
	},
	'psf':{
		'parameters':{
			'psf_type':'GAUSSIAN',
			'fwhm': 0.1
		}
	},
	'detector':{
		'parameters':{
			'pixel_scale':0.08,'ccd_gain':2.5,'read_noise':4.0,
			'magnitude_zero_point':output_ab_zeropoint,
			'exposure_time':5400.0,'sky_brightness':22,
			'num_exposures':1, 'background_noise':None
		}
	},
	'cross_object':{
		'parameters':{
			'los:delta_los,subhalo:sigma_sub':tmn
		}
	}
}

max_stars_repo_name: swagnercarena/paltas

max_issues_repo_name: swagnercarena/paltas

max_forks_repo_name: swagnercarena/paltas

reason: import numpy,from scipy


text:
import logging
import sys
import time

from PySide6 import QtCore, QtWidgets
import numpy as np
from PySide6.QtCore import qInstallMessageHandler
from matplotlib.figure import Figure

import faulthandler

faulthandler.enable()

from .matplotlibqml import FigureCanvasQTAgg, DemoViewModel, myMessageOutput

class ApplicationWindow(QtWidgets.QMainWindow):
    def __init__(self, vm: DemoViewModel):
        super().__init__()
        self.vm=vm
        self._main = QtWidgets.QWidget()
        self.setCentralWidget(self._main)
        layout = QtWidgets.QVBoxLayout(self._main)

        toolbar = QtWidgets.QHBoxLayout()

        home_btn=QtWidgets.QPushButton(text='home')
        home_btn.clicked.connect(vm.home)
        toolbar.addWidget(home_btn)

        back_btn=QtWidgets.QPushButton(text='back')
        back_btn.clicked.connect(vm.back)
        toolbar.addWidget(back_btn)

        forward_btn=QtWidgets.QPushButton(text='forward')
        forward_btn.clicked.connect(vm.forward)
        toolbar.addWidget(forward_btn)

        pan_btn=QtWidgets.QPushButton(text='pan')
        pan_btn.clicked.connect(vm.pan)
        toolbar.addWidget(pan_btn)

        zoom_btn=QtWidgets.QPushButton(text='zoom')
        zoom_btn.clicked.connect(vm.zoom)
        toolbar.addWidget(zoom_btn)

        pause_chkbtn=QtWidgets.QPushButton(text='pause')
        pause_chkbtn.setCheckable(True)
        pause_chkbtn.toggled.connect(vm.pauseChanged)
        toolbar.addWidget(pause_chkbtn)

        text=QtWidgets.QLabel()
        text.setFixedWidth(128)
        toolbar.addWidget(text)

        vm.coordinatesChanged.connect(lambda x: text.setText(x))

        layout.addLayout(toolbar)

        dynamic_canvas = FigureCanvasQTAgg(Figure(figsize=(5, 3)))
        layout.addWidget(dynamic_canvas)

        vm.updateWithCanvas(canvas=dynamic_canvas, dynamic=True)

if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="[%(levelname)s] %(message)s")

    qInstallMessageHandler(myMessageOutput)

    # Check whether there is already a running QApplication (e.g., if running
    # from an IDE).
    qapp = QtWidgets.QApplication.instance()
    if not qapp:
        qapp = QtWidgets.QApplication(sys.argv)

    app = ApplicationWindow(vm=DemoViewModel())
    app.show()
    app.activateWindow()
    app.raise_()
    qapp.exec_()

max_stars_repo_name: medlab/matplotlibqml

max_issues_repo_name: medlab/matplotlibqml

max_forks_repo_name: medlab/matplotlibqml


text:
#-*-coding:utf-8-*-
# date:2020-04-11
# Author: xiang li

import os
import shutil
import cv2
import numpy as np
import json
import torch
from config import cfg
from layers.functions.prior_box import PriorBox
from utils.box_utils import decode
def mkdir_(path, flag_rm=False):
    if os.path.exists(path):
        if flag_rm == True:
            shutil.rmtree(path)
            os.mkdir(path)
            print('remove {} done ~ '.format(path))
    else:
        os.mkdir(path)

def plot_box(bbox, img, color=None, label=None, line_thickness=None):
    tl = line_thickness or round(0.002 * max(img.shape[0:2])) + 1
    color = color or [random.randint(0, 255) for _ in range(3)]
    c1, c2 = (int(bbox[0]), int(bbox[1])), (int(bbox[2]), int(bbox[3]))
    cv2.rectangle(img, c1, c2, color, thickness=tl)# 目标的bbox
    if label:
        tf = max(tl - 2, 1)
        t_size = cv2.getTextSize(label, 0, fontScale=tl / 3, thickness=tf)[0] # label size
        c2 = c1[0] + t_size[0], c1[1] - t_size[1] - 3 # 字体的bbox
        cv2.rectangle(img, c1, c2, color, -1)  # label 矩形填充
        # 文本绘制
        cv2.putText(img, label, (c1[0], c1[1] - 2), 0, tl / 4, [225, 255, 255],thickness=tf, lineType=cv2.LINE_AA)

class JSON_Encoder(json.JSONEncoder):
    def default(self, obj):
        if isinstance(obj, np.integer):
            return int(obj)
        elif isinstance(obj, np.floating):
            return float(obj)
        elif isinstance(obj, np.ndarray):
            return obj.tolist()
        else:
            return super(JSON_Encoder, self).default(obj)

def draw_landmarks(img,output,r_bboxes,draw_circle):
    img_width = img.shape[1]
    img_height = img.shape[0]
    dict_landmarks = {}

    # 眼镜、鼻子、嘴角5个关键点
    faceswap_list = []
    faceswap_list.append((output[96*2+0]*float(img_width)+r_bboxes[0],output[96*2+1]*float(img_height)+r_bboxes[1]))
    faceswap_list.append((output[97*2+0]*float(img_width)+r_bboxes[0],output[97*2+1]*float(img_height)+r_bboxes[1]))
    faceswap_list.append((output[54*2+0]*float(img_width)+r_bboxes[0],output[54*2+1]*float(img_height)+r_bboxes[1]))
    faceswap_list.append((output[76*2+0]*float(img_width)+r_bboxes[0],output[76*2+1]*float(img_height)+r_bboxes[1]))
    faceswap_list.append((output[82*2+0]*float(img_width)+r_bboxes[0],output[82*2+1]*float(img_height)+r_bboxes[1]))
    # 98个关键点遍历，方便绘制每个部分的轮廓
    for i in range(int(output.shape[0]/2)):
        x = output[i*2+0]*float(img_width)
        y = output[i*2+1]*float(img_height)
        #
        if 41>= i >=33:
            if 'left_eyebrow' not in dict_landmarks.keys():
                dict_landmarks['left_eyebrow'] = []
            dict_landmarks['left_eyebrow'].append([int(x),int(y),(0,255,0)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (0,255,0),-1)
        elif 50>= i >=42:
            if 'right_eyebrow' not in dict_landmarks.keys():
                dict_landmarks['right_eyebrow'] = []
            dict_landmarks['right_eyebrow'].append([int(x),int(y),(0,255,0)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (0,255,0),-1)
        elif 67>= i >=60:
            if 'left_eye' not in dict_landmarks.keys():
                dict_landmarks['left_eye'] = []
            dict_landmarks['left_eye'].append([int(x),int(y),(255,55,255)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (255,0,255),-1)
        elif 75>= i >=68:
            if 'right_eye' not in dict_landmarks.keys():
                dict_landmarks['right_eye'] = []
            dict_landmarks['right_eye'].append([int(x),int(y),(255,55,255)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (255,0,255),-1)
        elif 97>= i >=96:
            cv2.circle(img, (int(x),int(y)), 2, (0,0,255),-1)
        elif 54>= i >=51:
            if 'bridge_nose' not in dict_landmarks.keys():
                dict_landmarks['bridge_nose'] = []
            dict_landmarks['bridge_nose'].append([int(x),int(y),(0,170,255)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (0,170,255),-1)
        elif 32>= i >=0:
            if 'basin' not in dict_landmarks.keys():
                dict_landmarks['basin'] = []
            dict_landmarks['basin'].append([int(x),int(y),(255,30,30)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (255,30,30),-1)
        elif 59>= i >=55:
            if 'wing_nose' not in dict_landmarks.keys():
                dict_landmarks['wing_nose'] = []
            dict_landmarks['wing_nose'].append([int(x),int(y),(0,255,255)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (0,255,255),-1)
        elif 87>= i >=76:
            if 'out_lip' not in dict_landmarks.keys():
                dict_landmarks['out_lip'] = []
            dict_landmarks['out_lip'].append([int(x),int(y),(255,255,0)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (255,255,0),-1)
        elif 95>= i >=88:
            if 'in_lip' not in dict_landmarks.keys():
                dict_landmarks['in_lip'] = []
            dict_landmarks['in_lip'].append([int(x),int(y),(50,220,255)])
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (50,220,255),-1)
        else:
            if draw_circle:
                cv2.circle(img, (int(x),int(y)), 2, (255,0,255),-1)

    faceswap_list_e = []

    for i in range(5):
        faceswap_list_e.append(faceswap_list[i][0])
    for i in range(5):
        faceswap_list_e.append(faceswap_list[i][1])

    return dict_landmarks,faceswap_list_e

def draw_contour(image,dict,r_bbox):
    x0 = r_bbox[0]# 全图偏置
    y0 = r_bbox[1]

    for key in dict.keys():
        # print(key)
        _,_,color = dict[key][0]

        if 'left_eye' == key:
            eye_x = np.mean([dict[key][i][0]+x0 for i in range(len(dict[key]))])
            eye_y = np.mean([dict[key][i][1]+y0 for i in range(len(dict[key]))])
            cv2.circle(image, (int(eye_x),int(eye_y)), 3, (255,255,55),-1)
        if 'right_eye' == key:
            eye_x = np.mean([dict[key][i][0]+x0 for i in range(len(dict[key]))])
            eye_y = np.mean([dict[key][i][1]+y0 for i in range(len(dict[key]))])
            cv2.circle(image, (int(eye_x),int(eye_y)), 3, (255,215,25),-1)

        if 'basin' == key or 'wing_nose' == key:
            pts = np.array([[dict[key][i][0]+x0,dict[key][i][1]+y0] for i in range(len(dict[key]))],np.int32)
            # print(pts)
            cv2.polylines(image,[pts],False,color,thickness = 2)

        else:
            points_array = np.zeros((1,len(dict[key]),2),dtype = np.int32)
            for i in range(len(dict[key])):
                x,y,_ = dict[key][i]
                points_array[0,i,0] = x+x0
                points_array[0,i,1] = y+y0

            # cv2.fillPoly(image, points_array, color)
            cv2.drawContours(image,points_array,-1,color,thickness=2)
def refine_face_bbox(bbox,img_shape):
    height,width,_ = img_shape

    x1,y1,x2,y2 = bbox

    expand_w = (x2-x1)
    expand_h = (y2-y1)

    x1 -= expand_w*0.06
    y1 -= expand_h*0.15
    x2 += expand_w*0.06
    y2 += expand_h*0.06

    x1,y1,x2,y2 = int(x1),int(y1),int(x2),int(y2)

    x1 = int(max(0,x1))
    y1 = int(max(0,y1))
    x2 = int(min(x2,width-1))
    y2 = int(min(y2,height-1))

    return (x1,y1,x2,y2)
def py_cpu_nms(dets, thresh):
    """Pure Python NMS baseline."""
    x1 = dets[:, 0]
    y1 = dets[:, 1]
    x2 = dets[:, 2]
    y2 = dets[:, 3]
    scores = dets[:, 4]

    areas = (x2 - x1 + 1) * (y2 - y1 + 1)
    order = scores.argsort()[::-1]

    keep = []
    while order.size > 0:
        i = order[0]
        keep.append(i)
        xx1 = np.maximum(x1[i], x1[order[1:]])
        yy1 = np.maximum(y1[i], y1[order[1:]])
        xx2 = np.minimum(x2[i], x2[order[1:]])
        yy2 = np.minimum(y2[i], y2[order[1:]])

        w = np.maximum(0.0, xx2 - xx1 + 1)
        h = np.maximum(0.0, yy2 - yy1 + 1)
        inter = w * h
        ovr = inter / (areas[i] + areas[order[1:]] - inter)

        inds = np.where(ovr <= thresh)[0]
        order = order[inds + 1]

    return keep

def check_keys(model, pretrained_state_dict):
    ckpt_keys = set(pretrained_state_dict.keys())
    model_keys = set(model.state_dict().keys())
    used_pretrained_keys = model_keys & ckpt_keys
    unused_pretrained_keys = ckpt_keys - model_keys
    missing_keys = model_keys - ckpt_keys
    print('Missing keys:{}'.format(len(missing_keys)))
    print('Unused checkpoint keys:{}'.format(len(unused_pretrained_keys)))
    print('Used keys:{}'.format(len(used_pretrained_keys)))
    assert len(used_pretrained_keys) > 0, 'load NONE from pretrained checkpoint'
    return True

def remove_prefix(state_dict, prefix):
    ''' Old style model is stored with all names of parameters sharing common prefix 'module.' '''
    print('remove prefix \'{}\''.format(prefix))
    f = lambda x: x.split(prefix, 1)[-1] if x.startswith(prefix) else x
    return {f(key): value for key, value in state_dict.items()}

def load_model(model, pretrained_path, load_to_cpu):
    print('Loading pretrained model from {}'.format(pretrained_path))
    if load_to_cpu:
        pretrained_dict = torch.load(pretrained_path, map_location=lambda storage, loc: storage)
    else:
        device = torch.cuda.current_device()
        pretrained_dict = torch.load(pretrained_path, map_location=lambda storage, loc: storage.cuda(device))
    if "state_dict" in pretrained_dict.keys():
        pretrained_dict = remove_prefix(pretrained_dict['state_dict'], 'module.')
    else:
        pretrained_dict = remove_prefix(pretrained_dict, 'module.')
    check_keys(model, pretrained_dict)
    model.load_state_dict(pretrained_dict, strict=False)
    return model

def detect_faces(ops,detect_model,img_raw,device):
    resize = 1
    img = np.float32(img_raw)
    if resize != 1:
        img = cv2.resize(img, None, None, fx=resize, fy=resize, interpolation=cv2.INTER_LINEAR)
    im_height, im_width, _ = img.shape
    scale = torch.Tensor([img.shape[1], img.shape[0], img.shape[1], img.shape[0]])
    img -= (104, 117, 123)
    img = img.transpose(2, 0, 1)
    img = torch.from_numpy(img).unsqueeze(0)
    img = img.to(device)
    scale = scale.to(device)

    loc, conf = detect_model(img)  # forward pass

    priorbox = PriorBox(cfg, image_size=(im_height, im_width))
    priors = priorbox.forward()
    priors = priors.to(device)
    prior_data = priors.data
    boxes = decode(loc.data.squeeze(0), prior_data, cfg['variance'])
    boxes = boxes * scale / resize
    boxes = boxes.cpu().numpy()
    scores = conf.squeeze(0).data.cpu().numpy()[:, 1]

    # ignore low scores
    inds = np.where(scores > ops.confidence_threshold)[0]
    boxes = boxes[inds]
    scores = scores[inds]

    # keep top-K before NMS
    order = scores.argsort()[::-1][:ops.top_k]
    boxes = boxes[order]
    scores = scores[order]

    # do NMS
    dets = np.hstack((boxes, scores[:, np.newaxis])).astype(np.float32, copy=False)
    #keep = py_cpu_nms(dets, ops.nms_threshold)
    # keep = nms(dets, ops.nms_threshold,force_cpu=True)
    keep = py_cpu_nms(dets, ops.nms_threshold)
    dets = dets[keep, :]

    # keep top-K faster NMS
    dets = dets[:ops.keep_top_k, :]

    return dets

def get_faces_batch_landmarks(ops,landmarks_model,dets,img_raw,scale_img,use_cuda,draw_bbox = True):
    # 绘制图像
    image_batch = None
    r_bboxes = []
    imgs_crop = []
    for b in dets:
        if b[4] < ops.vis_thres:
            continue
        text = "{:.4f}".format(b[4])
        b = list(map(int, b))

        r_bbox = refine_face_bbox((b[0],b[1],b[2],b[3]),img_raw.shape)
        r_bboxes.append(r_bbox)
        img_crop = img_raw[r_bbox[1]:r_bbox[3],r_bbox[0]:r_bbox[2]]
        imgs_crop.append(img_crop)
        img_ = cv2.resize(img_crop, (ops.landmarks_img_size[1],ops.landmarks_img_size[0]), interpolation = cv2.INTER_LINEAR) # INTER_LINEAR INTER_CUBIC

        img_ = img_.astype(np.float32)
        img_ = (img_-128.)/256.

        img_ = img_.transpose(2, 0, 1)
        img_ = np.expand_dims(img_,0)

        if image_batch is None:
            image_batch = img_
        else:
            image_batch = np.concatenate((image_batch,img_),axis=0)
    for b in dets:
        if b[4] < ops.vis_thres:
            continue
        text = "{:.4f}".format(b[4])
        b = list(map(int, b))
        if draw_bbox:
            cv2.rectangle(img_raw, (b[0], b[1]), (b[2], b[3]), (0, 0, 255), 2)
        cx = b[0]
        cy = b[1] - 3
        cv2.putText(img_raw, text, (cx, cy),cv2.FONT_HERSHEY_DUPLEX, 0.6, (155, 155, 255),3)
        cv2.putText(img_raw, text, (cx, cy),cv2.FONT_HERSHEY_DUPLEX, 0.6, (155, 10, 10),1)

    # 填充最大 关键点 批次数据
    if len(dets) < ops.max_batch_size:
        im_mask = np.zeros([1,3,ops.landmarks_img_size[0],ops.landmarks_img_size[1]], dtype = np.float32)
        for i in range(ops.max_batch_size-len(dets)):
            if image_batch is None:
                image_batch = im_mask
            else:
                image_batch = np.concatenate((image_batch,im_mask),axis=0)

    image_batch = torch.from_numpy(image_batch).float()

    if use_cuda:
        image_batch = image_batch.cuda()  # (bs, 3, h, w)

    pre_ = landmarks_model(image_batch.float())
    # print(pre_.size())
    output = pre_.cpu().detach().numpy()
    # print('output shape : ',output.shape)
    # n_array = np.zeros([ops.landmarks_img_size[0],ops.landmarks_img_size[1],3], dtype = np.float)
    faceswap_landmarks = []
    for i in range(len(dets)):

        dict_landmarks,list_e = draw_landmarks(imgs_crop[i],output[i],r_bboxes[i],draw_circle = False)
        faceswap_landmarks.append(list_e)

        draw_contour(img_raw,dict_landmarks,r_bboxes[i])
    # print(dets)
    #-----------------------------------------------------------------------------------

    fs_dets = []
    fs_landmarks = []

    for  i in range(len(r_bboxes)):
        fs_dets.append((r_bboxes[i][0]/scale_img,r_bboxes[i][1]/scale_img,r_bboxes[i][2]/scale_img,r_bboxes[i][3]/scale_img))
        fs_land = []
        for k in range(10):
            fs_land.append(faceswap_landmarks[i][k]/scale_img)
        fs_landmarks.append(fs_land)

        bbox = r_bboxes[i]
        min_x = int(bbox[0])
        min_y = int(bbox[1])
        max_x = int(bbox[2])
        max_y = int(bbox[3])
        # cv2.rectangle(img_raw, (min_x, min_y), (max_x, max_y), (255, 0, 255), thickness=11)
        for k in range(5):
            x = int(faceswap_landmarks[i][k+0])
            y = int(faceswap_landmarks[i][k+5])
            if draw_bbox:
                cv2.circle(img_raw,(x,y),5+k*2,(0,0,255),-1)

    cv2.putText(img_raw, 'face:'+str(len(dets)), (3,35),cv2.FONT_HERSHEY_DUPLEX, 1.45, (55, 255, 255),5)
    cv2.putText(img_raw, 'face:'+str(len(dets)), (3,35),cv2.FONT_HERSHEY_DUPLEX, 1.45, (135, 135, 5),2)
    return fs_dets,fs_landmarks

max_stars_repo_name: XLEric/FaceSwap

max_issues_repo_name: XLEric/FaceSwap

max_forks_repo_name: XLEric/FaceSwap


text:
import os
import hydra
import torch
import numpy as np
import random
from hydra import utils
from torch.utils.data import DataLoader
from deepke.name_entity_re.multimodal.models.IFA_model import IFANERCRFModel
from deepke.name_entity_re.multimodal.modules.dataset import MMPNERProcessor, MMPNERDataset
from deepke.name_entity_re.multimodal.modules.train import Trainer

import warnings
warnings.filterwarnings("ignore", category=UserWarning)
import logging
logging.basicConfig(format = '%(asctime)s - %(levelname)s - %(name)s -   %(message)s',
                    datefmt = '%m/%d/%Y %H:%M:%S',
                    level = logging.INFO)
logger = logging.getLogger(__name__)

# import wandb
# writer = wandb.init(project="DeepKE_NER_MM")
writer=None

DATA_PATH = {
    'twitter15': {'train': 'data/twitter2015/train.txt',
                'dev': 'data/twitter2015/valid.txt',
                'test': 'data/twitter2015/test.txt',
                'train_auximgs': 'data/twitter2015/twitter2015_train_dict.pth',
                'dev_auximgs': 'data/twitter2015/twitter2015_val_dict.pth',
                'test_auximgs': 'data/twitter2015/twitter2015_test_dict.pth',
                'rcnn_img_path': 'data/twitter2015',
                'img2crop': 'data/twitter2015/twitter2015_img2crop.pth'},

    'twitter17': {'train': 'data/twitter2017/train.txt',
                'dev': 'data/twitter2017/valid.txt',
                'test': 'data/twitter2017/test.txt',
                'train_auximgs': 'data/twitter2017/twitter2017_train_dict.pth',
                'dev_auximgs': 'data/twitter2017/twitter2017_val_dict.pth',
                'test_auximgs': 'data/twitter2017/twitter2017_test_dict.pth',
                'rcnn_img_path': 'data/twitter2017',
                'img2crop': 'data/twitter2017/twitter17_img2crop.pth'}
    }

IMG_PATH = {
    'twitter15': 'data/twitter2015/twitter2015_images',
    'twitter17': 'data/twitter2017/twitter2017_images'
}

AUX_PATH = {
    'twitter15': {'train': 'data/twitter2015/twitter2015_aux_images/train/crops',
                'dev': 'data/twitter2015/twitter2015_aux_images/val/crops',
                'test': 'data/twitter2015/twitter2015_aux_images/test/crops'},

    'twitter17': {'train': 'data/twitter2017/twitter2017_aux_images/train/crops',
                'dev': 'data/twitter2017/twitter2017_aux_images/val/crops',
                'test': 'data/twitter2017/twitter2017_aux_images/test/crops'}
}

LABEL_LIST = ["O", "B-MISC", "I-MISC", "B-PER", "I-PER", "B-ORG", "I-ORG", "B-LOC", "I-LOC", "X", "[CLS]", "[SEP]"]

def set_seed(seed=2021):
    """set random seed"""
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True
    np.random.seed(seed)
    random.seed(seed)

@hydra.main(config_path="conf", config_name='config')
def main(cfg):
    cwd = utils.get_original_cwd()
    cfg.cwd = cwd
    print(cfg)

    set_seed(cfg.seed) # set seed, default is 1
    if cfg.save_path is not None:  # make save_path dir
        if not os.path.exists(cfg.save_path):
            os.makedirs(cfg.save_path, exist_ok=True)
    print(cfg)
    label_mapping = {label:idx for idx, label in enumerate(LABEL_LIST, 1)}
    label_mapping["PAD"] = 0
    data_path, img_path, aux_path = DATA_PATH[cfg.dataset_name], IMG_PATH[cfg.dataset_name], AUX_PATH[cfg.dataset_name]
    rcnn_img_path = DATA_PATH[cfg.dataset_name]['rcnn_img_path']

    processor = MMPNERProcessor(data_path, cfg)
    train_dataset = MMPNERDataset(processor, label_mapping, img_path, aux_path, rcnn_img_path, max_seq=cfg.max_seq, ignore_idx=cfg.ignore_idx, aux_size=cfg.aux_size, rcnn_size=cfg.rcnn_size, mode='train',cwd=cwd)
    train_dataloader = DataLoader(train_dataset, batch_size=cfg.batch_size, shuffle=True, num_workers=4, pin_memory=True)

    dev_dataset = MMPNERDataset(processor, label_mapping, img_path, aux_path, rcnn_img_path, max_seq=cfg.max_seq, ignore_idx=cfg.ignore_idx, aux_size=cfg.aux_size, rcnn_size=cfg.rcnn_size, mode='dev',cwd=cwd)
    dev_dataloader = DataLoader(dev_dataset, batch_size=cfg.batch_size, shuffle=False, num_workers=4, pin_memory=True)

    test_dataset = MMPNERDataset(processor, label_mapping, img_path, aux_path, rcnn_img_path, max_seq=cfg.max_seq, ignore_idx=cfg.ignore_idx, aux_size=cfg.aux_size, rcnn_size=cfg.rcnn_size, mode='test',cwd=cwd)
    test_dataloader = DataLoader(test_dataset, batch_size=cfg.batch_size, shuffle=False, num_workers=4, pin_memory=True)

    model = IFANERCRFModel(LABEL_LIST, cfg)

    trainer = Trainer(train_data=train_dataloader, dev_data=dev_dataloader, test_data=test_dataloader, model=model, label_map=label_mapping, args=cfg, logger=logger, writer=writer)
    trainer.train()

if __name__ == '__main__':
    main()

max_stars_repo_name: zjunlp/deepke-dre

max_issues_repo_name: zjunlp/deepke-dre

max_forks_repo_name: zjunlp/deepke-dre


text:
import IMLearn.learners.regressors.linear_regression
from IMLearn.learners.regressors import PolynomialFitting
from IMLearn.utils import split_train_test

import numpy as np
import pandas as pd
import plotly.express as px
import plotly.io as pio
import matplotlib.pyplot as plt
import plotly.graph_objects as go

pio.templates.default = "simple_white"

def load_data(filename: str) -> pd.DataFrame:
    """
    Load city daily temperature dataset and preprocess data.
    Parameters
    ----------
    filename: str
        Path to house prices dataset

    Returns
    -------
    Design matrix and response vector (Temp)
    """
    df = pd.read_csv(filename, parse_dates=['Date'])
    df = df[df["Temp"] > -70]
    df['DayofYear'] = df['Date'].dt.dayofyear
    df = df.drop("Date", axis=1)
    return df

if __name__ == '__main__':
    np.random.seed(0)
    # Question 1 - Load and preprocessing of city temperature dataset
    temp_df = load_data("C:\\Users\\klain\\IML.HUJI\\datasets\\City_Temperature.csv")

    # Question 2 - Exploring data for specific country
    # part 1: the polynomial degree should be 2, the plot looks like a Parabola
    israel_df = temp_df.loc[temp_df["Country"] == "Israel"]
    plot = px.scatter(data_frame=israel_df, x=israel_df.DayofYear, y=israel_df.Temp, color=israel_df.Year
                      ,title="Temperature change as function of DayofYear in Israel ")
    plot.show()
    # part 2:
    std_months_israel_df = israel_df.groupby(['Month']).agg(np.std).Temp
    plot2 = px.bar(data_frame=std_months_israel_df, x=std_months_israel_df.index, y=std_months_israel_df,
                   title="Standard deviation of temp in each month in israel",
                   labels={'x': "Month", 'y': 'Standard deviation of temp'})
    plot2.show()

    # Question 3 - Exploring differences between countries
    q3_df2 = temp_df.groupby(['Country', 'Month']).Temp.agg(["mean", "std"]).reset_index()
    plot3 = px.line(data_frame=q3_df2, x='Month', y='mean', color='Country', error_y='std',
                    title="Average Temperature in countries by Months")
    plot3.show()

    # Question 4 - Fitting model for different values of `k`
    X = israel_df.DayofYear
    y = israel_df.Temp
    train_X, train_y, test_X, test_y = split_train_test(X, y, train_proportion=0.75)
    loss_for_k = []
    for k in range(1, 11):
        poli = PolynomialFitting(k)
        poli.fit(train_X.to_numpy(), train_y.to_numpy())
        k_loss = poli.loss(test_X.to_numpy(), test_y.to_numpy())
        loss_for_k.append(k_loss)
        print(f"The loss for k={k} is: {k_loss}")
    plot4 = px.bar(x=range(1, 11), y=loss_for_k, title="Loss as function of degree",
                   labels={"x": "Degree", "y": "Loss"})
    plot4.show()

    # Question 5 - Evaluating fitted model on different countries
    loss_per_country = []
    q5poli = PolynomialFitting(5)
    q5poli.fit(X, y)
    no_israel_df = temp_df[temp_df.Country != "Israel"]
    countries = no_israel_df.Country.unique()
    for country in countries:
        country_df = no_israel_df[no_israel_df.Country == country]
        c_loss = q5poli.loss(country_df.DayofYear, country_df.Temp)
        loss_per_country.append(c_loss)
    plot5 = px.bar(x=countries, y=loss_per_country, title="Loss as function of Country",
                   labels={"x": "Country name", "y": "Loss"})
    plot5.show()


text:
from collectives.collectives import Collectives
from context.harp_context import HarpContext
from harp_constants import Type, PartitioningMode
from numpy import ndarray

class JavaCollectives(Collectives):

    def __init__(self, harp_ctx: HarpContext):
        self.hctx = harp_ctx

    def barrier(self, ctx_name: str, op_name: str) -> bool:
        pass

    def broadcast(self, ctx_name: str, op_name: str, data: ndarray, data_type: Type,
                  partition_mode: PartitioningMode,
                  bcast_worker_id: int, use_mst_bcast: bool) -> bool:
        self.hctx.broadcast(ctx_name, op_name, data, bcast_worker_id, use_mst_bcast)

max_stars_repo_name: IU-nanoBio/harp

max_issues_repo_name: IU-nanoBio/harp

max_forks_repo_name: IU-nanoBio/harp

reason: from numpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#  pivoted_chol.py Author "Nathan Wycoff <nathanbrwycoff@gmail.com>" Date 01.14.2020

## A pivoted cholesky function for kernel functions
import numpy as np

def pivoted_chol(get_diag, get_row, M, err_tol = 1e-6):
    """
    A simple python function which computes the Pivoted Cholesky decomposition/approximation of positive semi-definite operator. Only diagonal elements and select rows of that operator's matrix represenation are required.

    get_diag - A function which takes no arguments and returns the diagonal of the matrix when called.
    get_row - A function which takes 1 integer argument and returns the desired row (zero indexed).
    M - The maximum rank of the approximate decomposition; an integer.
    err_tol - The maximum error tolerance, that is difference between the approximate decomposition and true matrix, allowed. Note that this is in the Trace norm, not the spectral or frobenius norm.

    Returns: R, an upper triangular matrix of column dimension equal to the target matrix. It's row dimension will be at most M, but may be less if the termination condition was acceptably low error rather than max iters reached.
    """

    d = np.copy(get_diag())
    N = len(d)

    pi = list(range(N))

    R = np.zeros([M,N])

    err = np.sum(np.abs(d))

    m = 0
    while (m < M) and (err > err_tol):

        i = m + np.argmax([d[pi[j]] for j in range(m,N)])

        tmp = pi[m]
        pi[m] = pi[i]
        pi[i] = tmp

        R[m,pi[m]] = np.sqrt(d[pi[m]])
        Apim = get_row(pi[m])
        for i in range(m+1, N):
            if m > 0:
                ip = np.inner(R[:m,pi[m]], R[:m,pi[i]])
            else:
                ip = 0
            R[m,pi[i]] = (Apim[pi[i]] - ip) / R[m,pi[m]]
            d[pi[i]] -= pow(R[m,pi[i]],2)

        err = np.sum([d[pi[i]] for i in range(m+1,N)])
        m += 1

    R = R[:m,:]

    return(R)

# Example:
N = 10
S = np.random.normal(size=[N,N])
A = np.matmul(S.T, S)

# Important to copy as diag may not be writable depending on np version
get_diag = lambda: np.diag(A).copy()
get_row = lambda i: A[i,:]
R = pivoted_chol(get_diag, get_row, M = N)

# R.T @ R is numerically equal to A.
print(np.linalg.norm(np.matmul(R.T, R) - A))

max_stars_repo_name: NathanWycoff/PivotedCholesky

max_issues_repo_name: NathanWycoff/PivotedCholesky

max_forks_repo_name: NathanWycoff/PivotedCholesky


text:
#!/usr/bin/python
# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*-coding:Utf-8 -*

"""
  syntheticpopulation.py
  Auteur : Alexis Petit, PhD Student, Namur University
  Date : 2017 03 22
  Adaptation du code : Morgane Dumont, PhD Student, Namur University
  Date : 2018 02 26

  This scrip analyses a population of space debris in the GEO region.
"""

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import pandas as pd
from sklearn import linear_model
import numpy as np
from scipy import stats
import sys
import math
import os
from datetime import datetime
from scipy.stats import norm,lognorm, ks_2samp
from random import random

def read_simulation(name_file):
  """ Read the snapshot of a simulation of the space debris environment.

  INPUT
  -----

  name_file: string
    Name of the file containing the population.

  RETURN
  ------

  population: dataframe
    The population of space debris.

  """

  f = open(name_file,'r')
  data = f.readlines()
  f.close()

  yy = int(data[1].strip().split()[0])
  mm = int(data[1].strip().split()[1])
  dd = int(data[1].strip().split()[2])
  hh = int(data[1].strip().split()[3])
  mi = int(data[1].strip().split()[4])
  ss = int(data[1].strip().split()[5])
  date = datetime(yy,mm,dd,hh,mi,ss)

  headers = data[2].strip().split()
  population = []
  for row in data[3:]:
    if 'GEOTEST' not in row:
      population.append(row.strip().split())
  population = pd.DataFrame(population,columns=headers)
  population['A/M[m2/kg]'] = population['BC[m2/kg]'].astype('float64')/2.2

  category = {}
  for i in range(0,len(population)-1,1):
    name = population.loc[i]['NAME']
    if name not in category:
        category[name] = 1
    else:
        category[name] += 1

  print
  print '> Categories in the simulation:'
  print '> ',category
  print

  fig, ax = plt.subplots()
  ax.set_title(date.strftime('%y/%m/%d'),fontsize=14)
  ax.set_xlabel('RAAN [deg]',fontsize=14)
  ax.set_ylabel('i [deg]',fontsize=14)
  ax.spines["top"].set_visible(False)
  ax.spines["right"].set_visible(False)

  for name_selected in category:
    flag = population['NAME']==name_selected
    raan = population['RAAN[deg]'][flag]
    inc = population['i[deg]'][flag]
    ax.plot(raan,inc,'o',markersize=2,label=name_selected.replace('_',' '))
  ax.set_xlim([0,360])
  ax.set_ylim([0,20])
  plt.legend(loc='upper left',fontsize=14)
  plt.savefig('set_simulation_1')
  plt.close()

  return population

def read_controls(name_file,var_list,limits,nb_obj,change_n_obj):
  """ Infer contraints from data

  INPUTS
  ------

  name_file: string
    Name of the files used as contraints
  var_list: list
    Variables used [a,i,RAAN,...]
  limits: list

  RETURN
  ------

  controls: dataframe
    Frequencies used as contraints

  """

  #controls empty
  frequencies = {}
  laws = {}

  #read the files containing the contraints
  f = open(name_file,'r')
  data = f.readlines()
  f.close()

  yy = int(data[1].strip().split()[0])
  mm = int(data[1].strip().split()[1])
  dd = int(data[1].strip().split()[2])
  hh = int(data[1].strip().split()[3])
  mi = int(data[1].strip().split()[4])
  ss = int(data[1].strip().split()[5])
  date = datetime(yy,mm,dd,hh,mi,ss)

  headers = data[2].strip().split()
  population = []
  for row in data[3:]:
    if 'GEOTEST' not in row:
      population.append(row.strip().split())
  population = pd.DataFrame(population,columns=headers)
  population['A/M[m2/kg]'] = population['BC[m2/kg]'].astype('float64')/2.2

  category = {}
  for i in range(0,len(population)-1,1):
    name = population.loc[i]['NAME']
    if name not in category:
        category[name] = 1
    else:
        category[name] += 1

  fig, ax = plt.subplots()
  ax.set_title(date.strftime('%y/%m/%d'),fontsize=14)
  ax.set_xlabel('RAAN [deg]',fontsize=14)
  ax.set_ylabel('i [deg]',fontsize=14)
  ax.spines["top"].set_visible(False)
  ax.spines["right"].set_visible(False)

  for name_selected in category:
    flag = population['NAME']==name_selected
    raan = population['RAAN[deg]'][flag]
    inc = population['i[deg]'][flag]
    ax.plot(raan,inc,'o',markersize=2,label=name_selected.replace('_',' '))
  ax.set_xlim([0,360])
  ax.set_ylim([0,20])
  plt.legend(loc='upper left',fontsize=14)
  plt.savefig('set_simulation_2')
  plt.close()

  print '> Plot the distribution (constraints)'
  fig = plt.figure(figsize=(10,8))

  for k,var in enumerate(var_list):

    #we select only the objects with the tag 'EKRAN_2_DEB'
    flag = population['NAME']=='EKRAN_2_DEB'
    var_data = [float(i) for i in population[var][flag].tolist()]

    ax = fig.add_subplot(2,2,k+1)
    ax.grid(True)
    ax.set_title(var)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    #ax.hist(var_data,bins=25)
    xmin =  limits[var][0]
    xmax =  limits[var][-1]
    x = np.linspace(xmin, xmax, 100)
    if var == 'A/M[m2/kg]':
      #param = lognorm.fit(var_data)
      #laws[var] = param
      #p = lognorm.pdf(x,param[0])
      #mu,std = norm.fit(var_data)
      #laws[var] = [mu,std]
      #p = norm.pdf(x,mu,std)
      shape, loc, scale = lognorm.fit(var_data, floc=0)
      std = np.log(shape)
      mu = np.log(scale)
      p = lognorm.pdf(x,std,loc=0, scale = scale)
    else:
      mu,std = norm.fit(var_data)
      laws[var] = [mu,std]
      p = norm.pdf(x,mu,std)
    ax.plot(x,p,'k',linewidth=2)
    ax.set_xlim([xmin,xmax])

    if change_n_obj:
	# to fit the approximated laws:
    	frequencies[var] = create_controls(var,xmin,xmax,mu,std,nb_obj,limits[var])
    else:
	# to fit exactly the same structure as simu2
	frequencies[var]=create_cut(var,var_data,limits[var])

    #print 'freq ',var,frequencies[var]

  plt.savefig('controls')
  plt.close()
  flag = population['NAME']=='EKRAN_2_DEB'
  population = population[flag]
  return frequencies,laws, population

def comparison_clouds(name1,name2,limit_raan,limit_inc):

  #read the files containing the contraints
  f = open(name1,'r')
  data = f.readlines()
  f.close()

  yy = int(data[1].strip().split()[0])
  mm = int(data[1].strip().split()[1])
  dd = int(data[1].strip().split()[2])
  hh = int(data[1].strip().split()[3])
  mi = int(data[1].strip().split()[4])
  ss = int(data[1].strip().split()[5])
  date = datetime(yy,mm,dd,hh,mi,ss)

  headers = data[2].strip().split()
  population = []
  for row in data[3:]:
    if 'GEOTEST' not in row:
      population.append(row.strip().split())
  population = pd.DataFrame(population,columns=headers)
  population['A/M[m2/kg]'] = population['BC[m2/kg]']

  category = {}
  for i in range(0,len(population)-1,1):
    name = population.loc[i]['NAME']
    if name not in category:
        category[name] = 1
    else:
        category[name] += 1

  fig, ax = plt.subplots()
  ax.set_title(date.strftime('%y/%m/%d'),fontsize=14)
  ax.set_xlabel('RAAN [deg]',fontsize=14)
  ax.set_ylabel('i [deg]',fontsize=14)
  ax.spines["top"].set_visible(False)
  ax.spines["right"].set_visible(False)

  flag = population['NAME']=='EKRAN_2_DEB'
  raan = population['RAAN[deg]'][flag]
  inc = population['i[deg]'][flag]
  ax.plot(raan,inc,'o',markersize=4,label='Simulation 1')

  #read the files containing the contraints
  f = open(name2,'r')
  data = f.readlines()
  f.close()

  headers = data[2].strip().split()
  population = []
  for row in data[3:]:
    if 'GEOTEST' not in row:
      population.append(row.strip().split())
  population = pd.DataFrame(population,columns=headers)

  category = {}
  for i in range(0,len(population)-1,1):
    name = population.loc[i]['NAME']
    if name not in category:
        category[name] = 1
    else:
        category[name] += 1

  flag = population['NAME']=='EKRAN_2_DEB'
  raan = population['RAAN[deg]'][flag]
  inc = population['i[deg]'][flag]
  ax.plot(raan,inc,'o',markersize=4,color='darkorange',label='Simulation 2')
  for j in range(0,n_dim['i[deg]'],1):
    x = limits['i[deg]'][j]
    ax.axhline(x, color='k', linestyle='--')
  x = limits['i[deg]'][j+1]
  ax.axhline(x, color='k', linestyle='--')
  for k in range(0,n_dim['RAAN[deg]'],1):
    x = limits['RAAN[deg]'][k]
    ax.axvline(x, color='k', linestyle='--')
  x = limits['RAAN[deg]'][k+1]
  ax.axvline(x, color='k', linestyle='--')

  ax.set_xlim(limit_raan)
  ax.set_ylim(limit_inc)
  plt.legend(loc='upper left',fontsize=14)
  plt.savefig('comparison_simulation')
  plt.close()

def create_cut(var,var_data,limits):
  """ aggregate var_data within classes defined by limits[var]
  """
  frequencies = np.zeros(len(limits)-1)
  for data in var_data:
	found = False
	index = 0
	if data<=limits[len(limits)-1]:
		while found==False:
			if data <= limits[index+1]:
				frequencies[index]+=1
				found=True
			else:
				index+=1
	else:
		frequencies[len(limits)-2]+=1
  return frequencies.tolist()

def def_bounds(population,var_list,n_dim, quantile_cut):
    fig = plt.figure(figsize=(10,8))
    bounds = {}
    limits = {}
    frequencies = {}
    k=0
    for var in var_list:
      #print var

      population[var] = pd.to_numeric(population[var])
      if quantile_cut == True:
      	bounds[var],limits[var] = pd.qcut(population[var],n_dim[var],retbins=True,labels=False)
      else:
      	bounds[var],limits[var] = pd.cut(population[var],n_dim[var],retbins=True,labels=False)
      #print population[var]
      #print limits[var]
      #print bounds[var]
      #print frequencies[var]
      countsTemp = pd.value_counts(bounds[var])

      print 'count avant', countsTemp
      counts = [0 for i in range(0,n_dim[var])]
      for i in  range(0,len(counts),1):
	if i in countsTemp.index:
		counts[i]= countsTemp.loc[i]
      print 'count apres', counts

      #print counts
      frequencies[var] = counts
      #print frequencies[var]

      ax = fig.add_subplot(2,2,k+1)
      #ax.spines["top"].set_visible(False)
      #ax.spines["right"].set_visible(False)
      data_max = len(population[var])

      if var == 'A/M[m2/kg]':
        n, bins, patches = plt.hist(population[var],50,facecolor='green',alpha=0.75,normed=True)#,weights=weights)
        shape, loc, scale = lognorm.fit(population[var], floc=0)
        xmin =  min(population[var])
        xmax =  max(population[var])
        x = np.linspace(xmin,xmax,50)
        p = lognorm.pdf(x,shape,loc=0, scale = scale)
        plt.plot(x,p,'r--',linewidth=2)
      else:
        mu,sigma = norm.fit(population[var])
        weights = np.ones_like(population[var])/float(len(population[var]))
        n, bins, patches = plt.hist(population[var],50,facecolor='green',alpha=0.75,normed=1)#,weights=weights)
        p = mlab.normpdf(bins,mu,sigma)#*data_max
        plt.plot(bins,p,'r--',linewidth=2)

      for i in range(0,n_dim[var]+1,1):
        x = min(population[var])+((max(population[var])-min(population[var])))*i/n_dim[var]
        plt.axvline(x,color='k',linestyle='--')

      k += 1
      #ax.grid(True)
      ax.set_title(var)

    plt.savefig('initial_frequencies_testMo')
    plt.close()
    return frequencies,limits,bounds

def create_controls(var,vmin,vmax,mu,std,nb_obj,limits):
  """ From statistical data we compute frequencies

  """

  counter = 0
  pop = []
  while(len(pop)<nb_obj):
    if var=='A/M[m2/kg]':
      x = np.random.lognormal(mean=mu,sigma=std)
      #print x
    else:
      x = np.random.normal(loc=mu,scale=std)
    if x>vmin and x<vmax:
      pop.append(x)

  frequencies = np.zeros(len(limits)-1)
  for i in range(0,len(limits)-1,1):
    for var in pop:
      if (var>limits[i]) and (var<limits[i+1]):
        frequencies[i] += 1

  return frequencies.tolist()

def pop_2_cross_table(population,var_list,n_dim,quantile_cut):
    """ Create the initial cross table.

    INPUT
    -----

    population: dataframe
    var_list: list
    n_dim: integer

    RETURN
    ------

    cross_table: float-matrix
      The cross table generated from the initial population.
    frequencies: dataframe
      Count the frequencies for each variables.

    """

    frequencies, limits, bounds = def_bounds(population,var_list,n_dim, quantile_cut)

    #sys.exit()

    cross_table = np.zeros((n_dim['a[m]'],n_dim['i[deg]'],n_dim['RAAN[deg]'],n_dim['A/M[m2/kg]']))
    for i in range(0,len(limits['a[m]'])-1,1):
      for j in range(0,len(limits['i[deg]'])-1,1):
        for k in range(0,len(limits['RAAN[deg]'])-1,1):
          for l in range(0,len(limits['A/M[m2/kg]'])-1,1):
            flag = (bounds['a[m]']==i) & (bounds['i[deg]']==j) & (bounds['RAAN[deg]']==k) & (bounds['A/M[m2/kg]']==l)
            cross_table[i,j,k,l] = len(population[flag])

    return cross_table,frequencies,limits

def delete_zero_cell_problem(cross_table,frequencies, n_dim):
	for i in range(0,len(frequencies['a[m]']),1):
		#print i
		#print frequencies
		for j in range(0,len(frequencies['i[deg]']),1):
			for k in range(0,len(frequencies['RAAN[deg]']),1):
				for l in range(0,len(frequencies['A/M[m2/kg]']),1):
					if cross_table[i,j,k,l] == 0:
						cross_table[i,j,k,l] = 0.001
						frequencies['a[m]'][i] += cross_table[i,j,k,l]
						frequencies['i[deg]'][j] += cross_table[i,j,k,l]
						frequencies['RAAN[deg]'][k] += cross_table[i,j,k,l]
						frequencies['A/M[m2/kg]'][l] += cross_table[i,j,k,l]
	#print cross_table
	return cross_table, frequencies

def ipf_process(cross_table,frequencies,controls, n_dim):
    """ Apply the IPF process.

    INPUTS
    ------

    cross_table: float-matrix
    frequencies: dataframe
    controls: dataframe
    var_list: list

    RETURN
    ------

    new_cross_table: float-matrix

    """

    #cross_table, frequencies = delete_zero_cell_problem(cross_table,frequencies, n_dim)

    epsilon = 1E-13 # 10^-13 because eps machine is 10^-16 for each cell of contingency table (+-1000)
    print n_dim['a[m]']*n_dim['i[deg]']*n_dim['RAAN[deg]']*n_dim['A/M[m2/kg]']
    distance = 1
    list_distance = []
    counter = 0

    while (epsilon<distance) and (counter<200):
      #print "iteration",counter+1
      cross_table_saved = cross_table.copy()

      cross_table, frequencies = update_variable_cross_table(cross_table,frequencies,controls)
      total = frequencies

      distance = 0
      for i in range(0,n_dim['a[m]'],1):
        for j in range(0,n_dim['i[deg]'],1):
          for k in range(0,n_dim['RAAN[deg]'],1):
            for l in range(0,n_dim['A/M[m2/kg]'],1):
              distance += np.abs(cross_table[i,j,k,l]-cross_table_saved[i,j,k,l])

      #print distance
      list_distance.append(distance)
      counter += 1

    #print cross_table
    print 'fin IPF', frequencies

    fig, ax = plt.subplots()
    ax.grid(True)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    #plt.rc('xtick', labelsize=12)
    #plt.rc('ytick', labelsize=12)
    #plt.ticklabel_format(style='sci',useOffset=False)
    ax.set_xlabel('Iteration',fontsize=14)
    ax.set_ylabel('Distance',fontsize=14)
    ax.set_xlim([0,len(list_distance)-1])
    ax.set_yscale('log')
    ax.plot(list_distance)
    plt.savefig('convergence')
    plt.close()

    return cross_table

def update_variable_cross_table(cross_table,frequency,controls):
  """ Fitting of the cross table

  INPUTS
  ------

  cross_table: float-matrix
  frequency: dataframe
  controls: dataframes

  RETURN
  ------

  cross_table: float-matrix
  new_frequency: dataframe

  """

  new_frequency = frequency.copy()
  for item in new_frequency:
    dim = len(new_frequency[item])
    new_frequency[item] = np.zeros(dim).tolist()

  for i in range(0,len(frequency['a[m]']),1):
    for j in range(0,len(frequency['i[deg]']),1):
      for k in range(0,len(frequency['RAAN[deg]']),1):
        for l in range(0,len(frequency['A/M[m2/kg]']),1):
          if frequency['i[deg]'][j]==0:
            cross_table[i,j,k,l] = 0
          else:
            cross_table[i,j,k,l] = cross_table[i,j,k,l]*controls['i[deg]'][j]/frequency['i[deg]'][j]
          new_frequency['a[m]'][i] += cross_table[i,j,k,l]
          new_frequency['i[deg]'][j] += cross_table[i,j,k,l]
          new_frequency['RAAN[deg]'][k] += cross_table[i,j,k,l]
          new_frequency['A/M[m2/kg]'][l] += cross_table[i,j,k,l]
  frequency = new_frequency.copy()

  new_frequency = frequency.copy()
  for item in new_frequency:
    dim = len(new_frequency[item])
    new_frequency[item] = np.zeros(dim).tolist()

  for i in range(0,len(frequency['a[m]']),1):
    for j in range(0,len(frequency['i[deg]']),1):
      for k in range(0,len(frequency['RAAN[deg]']),1):
        for l in range(0,len(frequency['A/M[m2/kg]']),1):
          if frequency['RAAN[deg]'][k]!=0:
			cross_table[i,j,k,l] = cross_table[i,j,k,l]*controls['RAAN[deg]'][k]/frequency['RAAN[deg]'][k]
          new_frequency['a[m]'][i] += cross_table[i,j,k,l]
          new_frequency['i[deg]'][j] += cross_table[i,j,k,l]
          new_frequency['RAAN[deg]'][k] += cross_table[i,j,k,l]
          new_frequency['A/M[m2/kg]'][l] += cross_table[i,j,k,l]
  frequency = new_frequency.copy()

  new_frequency = frequency.copy()
  for item in new_frequency:
    dim = len(new_frequency[item])
    new_frequency[item] = np.zeros(dim).tolist()

  for i in range(0,len(frequency['a[m]']),1):
    for j in range(0,len(frequency['i[deg]']),1):
      for k in range(0,len(frequency['RAAN[deg]']),1):
        for l in range(0,len(frequency['A/M[m2/kg]']),1):
          if frequency['A/M[m2/kg]'][l]==0:
            cross_table[i,j,k,l] = 0
          else:
            cross_table[i,j,k,l] = cross_table[i,j,k,l]*controls['A/M[m2/kg]'][l]/frequency['A/M[m2/kg]'][l]
          new_frequency['a[m]'][i] += cross_table[i,j,k,l]
          new_frequency['i[deg]'][j] += cross_table[i,j,k,l]
          new_frequency['RAAN[deg]'][k] += cross_table[i,j,k,l]
          new_frequency['A/M[m2/kg]'][l] += cross_table[i,j,k,l]
  frequency = new_frequency.copy()

  new_frequency = frequency.copy()
  for item in new_frequency:
    dim = len(new_frequency[item])
    new_frequency[item] = np.zeros(dim).tolist()
  for i in range(0,len(frequency['a[m]']),1):
    for j in range(0,len(frequency['i[deg]']),1):
      for k in range(0,len(frequency['RAAN[deg]']),1):
        for l in range(0,len(frequency['A/M[m2/kg]']),1):
          if frequency['a[m]'][i]==0:
            cross_table[i,j,k,l]=0
          else:
            cross_table[i,j,k,l] = cross_table[i,j,k,l]*controls['a[m]'][i]/frequency['a[m]'][i]
          new_frequency['a[m]'][i] += cross_table[i,j,k,l]
          new_frequency['i[deg]'][j] += cross_table[i,j,k,l]
          new_frequency['RAAN[deg]'][k] += cross_table[i,j,k,l]
          new_frequency['A/M[m2/kg]'][l] += cross_table[i,j,k,l]
  frequency = new_frequency.copy()

  return cross_table,new_frequency

def trs(cross_table,n_dim,nb_obj):

  new_cross_table = np.zeros((n_dim['a[m]'],n_dim['i[deg]'],n_dim['RAAN[deg]'],n_dim['A/M[m2/kg]']))
  weight = np.zeros((n_dim['a[m]'],n_dim['i[deg]'],n_dim['RAAN[deg]'],n_dim['A/M[m2/kg]']))
  sum_weight = 0

  total = 0
  for i in range(0,n_dim['a[m]'],1):
    for j in range(0,n_dim['i[deg]'],1):
      for k in range(0,n_dim['RAAN[deg]'],1):
        for l in range(0,n_dim['A/M[m2/kg]'],1):
          new_cross_table[i,j,k,l] += int(cross_table[i,j,k,l])
          total += new_cross_table[i,j,k,l]
          weight[i,j,k,l] = cross_table[i,j,k,l] - new_cross_table[i,j,k,l]
          sum_weight += weight[i,j,k,l]

  bound = np.zeros((n_dim['a[m]'],n_dim['i[deg]'],n_dim['RAAN[deg]'],n_dim['A/M[m2/kg]']))
  sum_bound = 0
  for i in range(0,n_dim['a[m]'],1):
    for j in range(0,n_dim['i[deg]'],1):
      for k in range(0,n_dim['RAAN[deg]'],1):
        for l in range(0,n_dim['A/M[m2/kg]'],1):
          weight[i,j,k,l] = weight[i,j,k,l]/sum_weight
          sum_bound += weight[i,j,k,l]
          bound[i,j,k,l] = sum_bound

  while (total<nb_obj):
    flag = False
    for i in range(0,n_dim['a[m]'],1):
      for j in range(0,n_dim['i[deg]'],1):
        for k in range(0,n_dim['RAAN[deg]'],1):
          for l in range(0,n_dim['A/M[m2/kg]'],1):
            if (bound[i,j,k,l] > random()) and new_cross_table[i,j,k,l]!=0:
               new_cross_table[i,j,k,l] +=1
               total += 1
               flag = True
            if flag:
              break
          if flag:
            break
        if flag:
          break
      if flag:
        break

  return new_cross_table

def cross_table_2_pop(cross_table,frequencies,limits,laws,var_list):
  """ Convert the cross table to a population of objects

  """
  print 'laws :',laws
  population = []
  for i in range(0,len(frequencies['a[m]']),1):
    for j in range(0,len(frequencies['i[deg]']),1):
      for k in range(0,len(frequencies['RAAN[deg]']),1):
        for l in range(0,len(frequencies['A/M[m2/kg]']),1):
          counter = 0
          while counter<cross_table[i,j,k,l]:
            sma = limits['a[m]'][i] + (limits['a[m]'][i+1]-limits['a[m]'][i])*random()
            while(True):
              inc = np.random.normal(laws['i[deg]'][0],laws['i[deg]'][1],1)
              if (inc>=limits['i[deg]'][j] and inc<limits['i[deg]'][j+1]):
                break
            ecc = 0.1*random()
            while(True):
              raan = np.random.normal(laws['RAAN[deg]'][0],laws['RAAN[deg]'][1],1)
              if (raan>=limits['RAAN[deg]'][k] and raan<limits['RAAN[deg]'][k+1]):
                break
            omega = 360.*random()
            ma = 360.*random()
            bc = limits['A/M[m2/kg]'][l] + (limits['A/M[m2/kg]'][l+1]-limits['A/M[m2/kg]'][l])*random()
            population.append([sma,ecc,inc,raan,omega,ma,bc])
            counter += 1

  #['a[m]', 'e[-]', 'i[deg]', 'RAAN[deg]', 'Omega[deg]', 'MA[deg]', 'A/M[m2/kg]', 'S[m]', 'IDNORAD', 'NAME']
  headers = ['a[m]', 'e[-]', 'i[deg]', 'RAAN[deg]', 'Omega[deg]', 'MA[deg]', 'A/M[m2/kg]']
  population = pd.DataFrame(population,columns=headers)

  return population

if __name__ == "__main__":

  # FIRST STEP: read the population coming from our simulation
  print
  print 'Read the data of the simulation'
  population = read_simulation('simulation1.txt')
  flag = population['NAME']=='EKRAN_2_DEB'
  population = population[flag]
  raan_init = pd.to_numeric(population['RAAN[deg]'])
  inc_init = pd.to_numeric(population['i[deg]'])
  print 'initial correlation between RAAN and i for sim 1 ', np.corrcoef(raan_init,inc_init)
  #print len(population)

  # SECOND STEP: create the cross-table
  print
  print 'Compute the cross-table'
  var_list = ['a[m]','i[deg]','RAAN[deg]','A/M[m2/kg]']
  n_dim = {}
  n_dim['a[m]'] = 1
  n_dim['i[deg]'] = 10
  n_dim['RAAN[deg]'] = 10
  n_dim['A/M[m2/kg]'] = 1#3
  quantile_cut = False
  change_n_obj = False
  nb_obj = 920#*2

  cross_table,frequencies,limits = pop_2_cross_table(population,var_list,n_dim,quantile_cut)
  #print 'Cross-table'
  #print cross_table
  print 'Frequencies'
  print frequencies
  print 'Limits'
  print limits
  print

  # THIRD STEP: calculate the controls
  print
  print 'Compute the contraints'
  controls,laws, population_constraint = read_controls('simulation2.txt',var_list,limits,nb_obj,change_n_obj)
  raan_init_cons = pd.to_numeric(population_constraint['RAAN[deg]'])
  inc_init_cons = pd.to_numeric(population_constraint['i[deg]'])
  print 'initial correlation between RAAN and i for sim 2 ', np.corrcoef(raan_init_cons,inc_init_cons)
  print '> ',controls
  print
  nb_obj = sum(controls['a[m]'])
  print nb_obj

  # FOURTH STEP: apply the IPF process
  new_cross_table = ipf_process(cross_table,frequencies,controls, n_dim)
  print 'New cross-table'
  #print cross_table

  # TRS process
  cross_table = trs(new_cross_table,n_dim,nb_obj)
  #cross_table = new_cross_table

  #check total
  total = 0
  for i in range(0,n_dim['a[m]'],1):
    for j in range(0,n_dim['i[deg]'],1):
      for k in range(0,n_dim['RAAN[deg]'],1):
        for l in range(0,n_dim['A/M[m2/kg]'],1):
          total += cross_table[i,j,k,l]
  print 'New total = ',total

  # FITH STEP: compute the synthetic population
  new_population = cross_table_2_pop(cross_table,frequencies,limits,laws,var_list)
  raan = new_population['RAAN[deg]']
  inc = new_population['i[deg]']

  print type(raan_init)
  print 'new correlation between RAAN and i', np.corrcoef([a[0] for a in raan],[a[0] for a in inc])
  raan_approx = np.random.normal(laws['RAAN[deg]'][0],laws['RAAN[deg]'][1],9600)
  inc_approx = np.random.normal(laws['i[deg]'][0],laws['i[deg]'][1],9600)
  print 'distrib raan - approx :', ks_2samp(raan_approx, [a[0] for a in raan])
  print 'distrib inc - approx :', ks_2samp(inc_approx, [a[0] for a in inc])

  ##############################################################
  # linear reg
  slope_init, intercept_init, r_value, p_value, std_err = stats.linregress(raan_init,inc_init)
  print('Coefficients linear reg sim 1: \n', slope_init, intercept_init)
  print('r_value : \n',r_value)
  slope_init_cons, intercept_init_cons, r_value, p_value, std_err = stats.linregress(raan_init_cons,inc_init_cons)
  print('Coefficients linear reg sim 1: \n', slope_init_cons, intercept_init_cons)
  print('r_value : \n',r_value)
  slope_synthPop, intercept_synthPop, r_value, p_value, std_err = stats.linregress([a[0] for a in raan],[a[0] for a in inc])
  print('Coefficients linear reg synthetic pop: \n', slope_synthPop, intercept_synthPop)
  print('r_value : \n',r_value)

  fig, ax = plt.subplots()
  ax.set_xlabel('RAAN [deg]',fontsize=14)
  ax.set_ylabel('i [deg]',fontsize=14)
  ax.spines["top"].set_visible(False)
  ax.spines["right"].set_visible(False)
  ax.plot(population['RAAN[deg]'],population['i[deg]'],'o',markersize=4,label='Simulation ('+str(len(population))+' fragments)',alpha=0.7)
  ax.plot(raan,inc,'o',markersize=4,color='g',label='Synthetic population ('+str(len(raan))+' fragments)',alpha=0.7)
  for j in range(0,n_dim['i[deg]'],1):
    x = limits['i[deg]'][j]
    ax.axhline(x, color='k', linestyle='--')
  x = limits['i[deg]'][j+1]
  ax.axhline(x, color='k', linestyle='--')
  for k in range(0,n_dim['RAAN[deg]'],1):
    x = limits['RAAN[deg]'][k]
    ax.axvline(x, color='k', linestyle='--')
  x = limits['RAAN[deg]'][k+1]
  ax.axvline(x, color='k', linestyle='--')
  ax.set_xlim([limits['RAAN[deg]'][0],limits['RAAN[deg]'][k+1]])
  ax.set_ylim([limits['i[deg]'][0],limits['i[deg]'][j+1]])
  plt.legend(loc='upper left',fontsize=12)
  plt.savefig('synthetic_population')
  plt.close()

################ figure with everything and regression
  fig, ax = plt.subplots()
  ax.set_xlabel('RAAN [deg]',fontsize=14)
  ax.set_ylabel('i [deg]',fontsize=14)
  ax.spines["top"].set_visible(False)
  ax.spines["right"].set_visible(False)
  ax.plot(population['RAAN[deg]'],population['i[deg]'],'o',markersize=4,color='blue',label='Simulation 1 ('+str(len(population))+' fragments)',alpha=0.7)
  X = np.array([min(raan_init), max(raan_init)])
  plt.plot(X, X*slope_init + intercept_init, 'b')
  ax.plot(population_constraint['RAAN[deg]'],population_constraint['i[deg]'],'o',markersize=4,color='darkorange',label='Simulation 2 ('+str(len(population_constraint))+' fragments)',alpha=0.7)
  plt.plot(X, X*slope_init_cons + intercept_init_cons ,color='darkorange')
  ax.plot(raan,inc,'o',markersize=4,label='Synthetic population ('+str(len(raan))+' fragments)',alpha=0.7,color='g')
  plt.plot(X, X*slope_synthPop + intercept_synthPop ,color='g')
  for j in range(0,n_dim['i[deg]'],1):
    x = limits['i[deg]'][j]
    ax.axhline(x, color='k', linestyle='--')
  x = limits['i[deg]'][j+1]
  ax.axhline(x, color='k', linestyle='--')
  for k in range(0,n_dim['RAAN[deg]'],1):
    x = limits['RAAN[deg]'][k]
    ax.axvline(x, color='k', linestyle='--')
  x = limits['RAAN[deg]'][k+1]
  ax.axvline(x, color='k', linestyle='--')
  ax.set_xlim([limits['RAAN[deg]'][0],limits['RAAN[deg]'][k+1]])
  ax.set_ylim([limits['i[deg]'][0],limits['i[deg]'][j+1]])
  plt.legend(loc='upper left',fontsize=12)
  plt.savefig('synthetic_population_all')
  plt.close()

  plt.subplot(1, 2, 1)
  plt.hist([a[0] for a in raan], histtype='stepfilled',color="orange", alpha=0.4, normed=True, bins=40, label="Synthetic population")
  plt.hist(raan_init, histtype='stepfilled', alpha=0.4, normed=True, bins=40, label="Simulation 1")
  plt.hist(raan_init_cons, histtype='stepfilled', alpha=0.4, normed=True, bins=40, label="Simulation 2")
  plt.xlabel('RAAN [deg]')
  plt.ylabel('Density')
  plt.legend(loc='upper left',fontsize=12)
  plt.ylim((0,1))
  plt.subplot(1, 2, 2)
  plt.hist([a[0] for a in inc], histtype='stepfilled',color="orange", alpha=0.4, normed=True, bins=40, label="Synthetic population")
  plt.hist(inc_init, histtype='stepfilled', alpha=0.4, normed=True, bins=40, label="Simulation 1")
  plt.hist(inc_init_cons, histtype='stepfilled', alpha=0.4, normed=True, bins=40, label="Simulation 2")
  plt.xlabel('i[deg]')
  plt.ylim((0,1.4))
  plt.savefig('Distributions_pop_appl2')
  plt.close()

  limit_raan = [limits['RAAN[deg]'][0],limits['RAAN[deg]'][k+1]]
  limit_inc = [limits['i[deg]'][0],limits['i[deg]'][j+1]]
  comparison_clouds('simulation1.txt','simulation2.txt',limit_raan,limit_inc)