text:
import numpy as np

g = 9.817

def position(t, omega, u0, v0, x0):
    """Finds position of puck on paraboloid. The equations come from
    Cushman-Roisin Chapter 2, in particular, 2.17a, and 2.17b.

    Parameters
    ----------
    t : float
        time of given moment [s]
    omega : float
        effective rotation [1/s]
    u0 : float
        initial x-component of the velocity [cm/s]
    v0 : float
        initial y-component of the velocity [cm/s]
    x0 : float
        initial x-component of the position [cm]

    Returns
    -------
    x, y, z : floats
        x, y, and z componets of the position [cm]
    """
    ot = omega * t

    x = x0 * np.cos(ot) + u0 * np.sin(ot) / omega
    y = v0 * np.sin(ot) / omega
    z = omega**2 / (2*g) * (x**2 + y**2)

    return x, y, z

def circle(r):
    """Creates circle that represents edge of the paraboloid.

    Parameters
    ----------
    r : float
        radius of the circle [cm]

    Returns
    -------
    xc, yc : array_like
        x and y values of the circle
    """
    twopi = 2 * np.pi

    k = np.arange(1, 102)
    xc = r * np.cos(k * twopi / 100)
    yc = r * np.sin(k * twopi / 100)

    return xc, yc

def parabola(r, omega):
    """Creates parabola that represents the sides of the paraboloid along y=0.

    Parameters
    ----------
    r : float
        radius of the paraboloid [cm]
    omega : float
        effective rotation [1/s]

    Returns:
    x, z : array_like
        x and z values of the parabola
    """

    x = np.linspace(start=-r,stop=r)
    z = omega**2 / (2 * g) * x**2

    return x, z

def check_edge(x, y, r):
    """Checks if puck is off the edge of the paraboloid.

    Parameters
    ----------
    x, y :  floats
        x and y positions of the puck [cm]
    r : float
        radius of the circle

    Returns
    -------
    bool
        True if puck crossed the radius of the paraboloid
    """

    if(r > np.sqrt(x**2 + y**2)):
        return True
    return False

max_stars_repo_name: DIYnamics/digipyro

max_issues_repo_name: DIYnamics/digipyro

max_forks_repo_name: DIYnamics/digipyro


text:
"""Run ensemble model on Annual 2011 Data.

Please open a data directory according to the directory
specified in _DATA_ADDR_PREFIX,
and place both training data and validation data there.

"""
import os
import time

from importlib import reload

import pickle as pk
import pandas as pd

import numpy as np
from scipy.stats import norm as norm_dist

from sklearn.cluster import KMeans
from sklearn.isotonic import IsotonicRegression
from sklearn.model_selection import KFold

import tensorflow as tf
import tensorflow_probability as tfp
from tensorflow_probability import edward2 as ed

# sys.path.extend([os.getcwd()])

from calibre.model import gaussian_process as gp
from calibre.model import tailfree_process as tail_free
from calibre.model import gp_regression_monotone as gpr_mono
from calibre.model import adaptive_ensemble

from calibre.inference import mcmc

from calibre.calibration import score

import calibre.util.misc as misc_util
import calibre.util.metric as metric_util
import calibre.util.visual as visual_util
import calibre.util.matrix as matrix_util
import calibre.util.ensemble as ensemble_util
import calibre.util.calibration as calib_util

import calibre.util.experiment_pred as pred_util

from calibre.util.inference import make_value_setter

import matplotlib.pyplot as plt
import seaborn as sns

os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'

_DATA_ADDR_PREFIX = "./example/data"
_SAVE_ADDR_PREFIX = "./result/calibre_2d_annual_pm25_example"

_MODEL_DICTIONARY = {"root": ["IK", "QD", "AV"]}

DEFAULT_LOG_LS_WEIGHT = np.log(0.35).astype(np.float32)
DEFAULT_LOG_LS_RESID = np.log(0.1).astype(np.float32)

"""""""""""""""""""""""""""""""""
# 0. Prepare data
"""""""""""""""""""""""""""""""""
# TODO(jereliu): re-verify the correspondence between prediction observation
os.makedirs(os.path.join(_SAVE_ADDR_PREFIX, 'base'), exist_ok=True)
if not os.path.isdir(_DATA_ADDR_PREFIX):
    raise ValueError("Data diretory {} doesn't exist!".format(_DATA_ADDR_PREFIX))

""" 0. prepare training data dictionary """
y_obs_2011 = pd.read_csv("{}/training_data_2011.csv".format(_DATA_ADDR_PREFIX))

X_train = np.asarray(y_obs_2011[["lon", "lat"]].values.tolist()).astype(np.float32)
y_train = np.asarray(y_obs_2011["pm25_obs"].tolist()).astype(np.float32)

base_train_feat = dict()
base_train_pred = dict()
for model_name in tail_free.get_leaf_model_names(_MODEL_DICTIONARY):
    base_train_feat[model_name] = X_train
    base_train_pred[model_name] = y_obs_2011["pred_{}".format(model_name)].astype(np.float32)

""" 1. prepare prediction data dictionary """
base_valid_feat = dict()
base_valid_pred = dict()
for model_name in tail_free.get_leaf_model_names(_MODEL_DICTIONARY):
    data_pd = pd.read_csv("{}/{}_2011_align.csv".format(
        _DATA_ADDR_PREFIX, model_name))
    base_valid_feat[model_name] = np.asarray(data_pd[["lon", "lat"]].values.tolist()).astype(np.float32)
    base_valid_pred[model_name] = np.asarray(data_pd["pm25"].tolist()).astype(np.float32)

X_valid = base_valid_feat[model_name]
N_pred = X_valid.shape[0]

with open(os.path.join(_SAVE_ADDR_PREFIX, 'base/base_train_feat.pkl'), 'wb') as file:
    pk.dump(base_train_feat, file, protocol=pk.HIGHEST_PROTOCOL)
with open(os.path.join(_SAVE_ADDR_PREFIX, 'base/base_train_pred.pkl'), 'wb') as file:
    pk.dump(base_train_pred, file, protocol=pk.HIGHEST_PROTOCOL)

with open(os.path.join(_SAVE_ADDR_PREFIX, 'base/base_valid_feat.pkl'), 'wb') as file:
    pk.dump(base_valid_feat, file, protocol=pk.HIGHEST_PROTOCOL)
with open(os.path.join(_SAVE_ADDR_PREFIX, 'base/base_valid_pred.pkl'), 'wb') as file:
    pk.dump(base_valid_pred, file, protocol=pk.HIGHEST_PROTOCOL)

""" 3. standardize data """
# standardize
X_centr = np.mean(X_valid, axis=0)
X_scale = np.max(X_valid, axis=0) - np.min(X_valid, axis=0)

X_valid = (X_valid - X_centr) / X_scale
X_train = (X_train - X_centr) / X_scale

"""""""""""""""""""""""""""""""""
# 1. Model Estimation using MCMC
"""""""""""""""""""""""""""""""""
family_name = "hmc"
family_name_full = "Hamilton MC"

os.makedirs("{}/{}".format(_SAVE_ADDR_PREFIX, family_name),
            exist_ok=True)

with open(os.path.join(_SAVE_ADDR_PREFIX, 'base/base_train_pred.pkl'), 'rb') as file:
    base_train_pred = pk.load(file)

with open(os.path.join(_SAVE_ADDR_PREFIX, 'base/base_valid_pred.pkl'), 'rb') as file:
    base_valid_pred = pk.load(file)

"""2.1. sampler basic config"""
family_tree_dict = _MODEL_DICTIONARY

num_mcmc_steps = 5000
num_burnin_steps = 1000

"""2.2. run mcmc estimation"""

# define mcmc computation graph
(mcmc_graph, init_op,
 parameter_samples, is_accepted) = (
    mcmc.make_inference_graph_tailfree(
        X_train, y_train,
        base_pred=base_train_pred,
        family_tree=family_tree_dict,
        default_log_ls_weight=DEFAULT_LOG_LS_WEIGHT,
        default_log_ls_resid=DEFAULT_LOG_LS_RESID,
        num_mcmc_samples=num_mcmc_steps,
        num_burnin_steps=num_burnin_steps))

""" 2.2. execute sampling"""
# this will take some time
parameter_samples_val = mcmc.run_sampling(mcmc_graph,
                                          init_op,
                                          parameter_samples,
                                          is_accepted)

with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_train_parameter_samples_dict.pkl'.format(family_name)), 'wb') as file:
    pk.dump(parameter_samples_val, file, protocol=pk.HIGHEST_PROTOCOL)

"""""""""""""""""""""""""""""""""
# 2. Perform Model Prediction
"""""""""""""""""""""""""""""""""
# load mcmc posterior samples
with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_train_parameter_samples_dict.pkl'.format(family_name)), 'rb') as file:
    parameter_samples_val = pk.load(file)

# extract parameters
sigma_sample_val = parameter_samples_val["sigma_sample"]
resid_sample_val = parameter_samples_val["ensemble_resid_sample"]
temp_sample_val = parameter_samples_val["temp_sample"]
weight_sample_val = parameter_samples_val["weight_sample"]

# since validation data is very large, perform prediction by data into batch,
kf = KFold(n_splits=20)

# prepare output container
ensemble_sample_val = np.zeros(shape=(X_valid.shape[0], num_mcmc_steps))
ensemble_mean_val = np.zeros(shape=(X_valid.shape[0], num_mcmc_steps))

for fold_id, (_, pred_index) in enumerate(kf.split(X_valid)):
    print("Running fold {} out of {}".format(fold_id + 1, kf.n_splits))

    # prepare X_pred and base_pred_dict for each batch
    X_pred_fold = X_valid[pred_index]
    base_pred_dict_fold = {
        model_name: model_pred_val[pred_index]
        for (model_name, model_pred_val) in base_valid_pred.items()}

    # run prediction routine
    (ensemble_sample_fold, ensemble_mean_fold,
     _, _, _) = (
        pred_util.prediction_tailfree(X_pred=X_pred_fold,
                                      base_pred_dict=base_pred_dict_fold,
                                      X_train=X_train,
                                      family_tree=family_tree_dict,
                                      weight_sample_list=weight_sample_val,
                                      resid_sample=resid_sample_val,
                                      temp_sample=temp_sample_val,
                                      default_log_ls_weight=DEFAULT_LOG_LS_WEIGHT,
                                      default_log_ls_resid=DEFAULT_LOG_LS_RESID, )
    )

    # save to output container
    ensemble_sample_val[pred_index] = ensemble_sample_fold.T
    ensemble_mean_val[pred_index] = ensemble_mean_fold.T

print("Estimated ls_weight {:.4f}, ls_resid {:.4f}".format(
    np.exp(DEFAULT_LOG_LS_WEIGHT), np.exp(DEFAULT_LOG_LS_RESID)
))

with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_pred_dist_sample.pkl'.format(family_name)), 'wb') as file:
    pk.dump(ensemble_sample_val, file, protocol=pk.HIGHEST_PROTOCOL)
with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_pred_mean_sample.pkl'.format(family_name)), 'wb') as file:
    pk.dump(ensemble_mean_val, file, protocol=pk.HIGHEST_PROTOCOL)
with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_sigma_sample.pkl'.format(family_name)), 'wb') as file:
    pk.dump(sigma_sample_val, file, protocol=pk.HIGHEST_PROTOCOL)

"""""""""""""""""""""""""""""""""
# 3. Visualization
"""""""""""""""""""""""""""""""""
""" 3.1. prep: load data, compute posterior mean/sd, color config """
with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_pred_dist_sample.pkl'.format(family_name)), 'rb') as file:
    ensemble_sample_val = pk.load(file)
with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_pred_mean_sample.pkl'.format(family_name)), 'rb') as file:
    ensemble_mean_val = pk.load(file)
with open(os.path.join(_SAVE_ADDR_PREFIX,
                       '{}/ensemble_posterior_sigma_sample.pkl'.format(family_name)), 'rb') as file:
    sigma_sample_val = pk.load(file)

post_uncn_dict = {
    "overall": np.var(ensemble_sample_val, axis=1) + np.mean(np.exp(2 * sigma_sample_val)),
    "mean": np.var(ensemble_mean_val, axis=1),
    "resid": np.var(ensemble_sample_val - ensemble_mean_val, axis=1),
    "noise": np.mean(np.exp(2 * sigma_sample_val)) * np.ones(shape=(ensemble_sample_val.shape[0]))
}

post_mean_dict = {
    "overall": np.mean(ensemble_sample_val, axis=1),
    "mean": np.mean(ensemble_mean_val, axis=1),
    "resid": np.mean(ensemble_sample_val - ensemble_mean_val, axis=1)
}

# prepare color norms for plt.scatter
color_norm_unc = visual_util.make_color_norm(
    list(post_uncn_dict.values())[:1],  # use "overall" and "mean" for pal
    method="percentile")
color_norm_ratio = visual_util.make_color_norm(
    post_uncn_dict["noise"] / post_uncn_dict["overall"],
    method="percentile")
color_norm_pred = visual_util.make_color_norm(
    list(post_mean_dict.values())[:2],  # exclude "resid" vales from pal
    method="percentile")

""" 3.1. posterior predictive uncertainty """
for unc_name, unc_value in post_uncn_dict.items():
    save_name = os.path.join(_SAVE_ADDR_PREFIX,
                             '{}/ensemble_posterior_uncn_{}.png'.format(
                                 family_name, unc_name))

    color_norm = visual_util.posterior_heatmap_2d(unc_value,
                                                  X=X_valid, X_monitor=X_train,
                                                  cmap='inferno_r',
                                                  norm=color_norm_unc,
                                                  norm_method="percentile",
                                                  save_addr=save_name)

""" 3.2. posterior predictive mean """
for mean_name, mean_value in post_mean_dict.items():
    save_name = os.path.join(_SAVE_ADDR_PREFIX,
                             '{}/ensemble_posterior_mean_{}.png'.format(
                                 family_name, mean_name))
    color_norm = visual_util.posterior_heatmap_2d(mean_value,
                                                  X=X_valid, X_monitor=X_train,
                                                  cmap='RdYlGn_r',
                                                  norm=color_norm_pred,
                                                  norm_method="percentile",
                                                  save_addr=save_name)

ext: py

max_stars_repo_name: jereliu/aden-tf

max_issues_repo_name: jereliu/aden-tf

max_forks_repo_name: jereliu/aden-tf

reason: import numpy,from scipy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jul  1 14:40:46 2020

@author: arsi

Calculate STL decomposition for given time series and plot the components.

"""
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import STL
from Source.Utils.argument_loader import setup_np, setup_pd
import pytest
from Source.Utils.plot_decorator import plot_decorator

#TODO! = fix the xlabels

@plot_decorator
def _plot_decomposition(Result,
                        title,
                        savepath = False,
                        savename = False,
                        ylabel = "Value",
                        xlabel  = "Date",
                        dates = False,
                        test = False,
                        ):
    """
    Plot the decomposed time series.

    Parameters
    ----------
    Result : statsmodels.tsa.seasonal.DecomposeResult object
        Object containing the decomposition results
    title : str
        Figure title.
    savepath : Path object, optional
        Figure save path The default is False.
    savename : str, optional
        Figure save name. The default is False.
    ylabel : str, optional
        Figure ylabel. The default is "Battery Level (%)".
    xlabel : str, optional
        Figure xlabel. The default is "Date".
    dates : array, optional
        List of daytes to be highlighted in the figure. The default is False.

    Raises
    ------
    Exception
        - savepath does not exist
        - savename or path was not given in correct format

    Returns
    -------
    None.

    """

    fig1 = plt.figure(figsize=(9.3,9.3))

    plt.suptitle(title,fontsize=22,y=1)

    plt.subplot(4,1,1)
    plt.plot(Result.observed)
    plt.title('Observations',fontsize=18)
    if type(dates) != bool:
        for d in dates:
            plt.axvline(x=d,linestyle =":", color ='black')
    plt.ylabel(ylabel,fontsize=14)
    plt.xlabel(xlabel,fontsize=14)

    plt.subplot(4,1,2)
    plt.plot(Result.trend)
    plt.title('Trend',fontsize=18)
    if type(dates) != bool:
        for d in dates:
            plt.axvline(x=d,linestyle =":", color ='black')
    plt.ylabel(ylabel,fontsize=14)
    plt.xlabel(xlabel,fontsize=14)

    plt.subplot(4,1,3)
    plt.plot(Result.seasonal)
    if type(dates) != bool:
        for d in dates:
            plt.axvline(x=d,linestyle =":", color ='black')
    plt.title('Seasonal',fontsize=18)
    plt.ylabel(ylabel,fontsize=14)
    plt.xlabel(xlabel,fontsize=14)

    plt.subplot(4,1,4)
    plt.plot(Result.resid)
    if type(dates) != bool:
        for d in dates:
            plt.axvline(x=d,linestyle =":", color ='black')
    plt.title('Residuals',fontsize=18)
    plt.ylabel(ylabel,fontsize=14)
    plt.xlabel(xlabel,fontsize=14)

    fig1.tight_layout(pad=2)

    return fig1

def STL_decomposition(series,
                      title,
                      test = False,
                      savepath = False,
                      savename = False,
                      ylabel = "Battery Level (%)",
                      xlabel  = "Date",
                      dates = False,
                      ):
    """
    STL Decompose timeseries into Model, Trend, Seasonal and Residual parts.
    Plot the components and their distributions. Optionally save the figure.

    Parameters
    ----------
    series : Numpy ndarray
        Time series to be decomposed
    title : str
        Figure title.
    savepath : Path object, optional
        Figure save path The default is False.
    savename : str, optional
        Figure save name. The default is False.
    ylabel : str, optional
        Figure ylabel. The default is "Battery Level (%)".
    xlabel : str, optional
        Figure xlabel. The default is "Date".
    dates : array, optional
        List of daytes to be highlighted in the figure. The default is False.

    Raises
    ------
    Exception
        - given series is not a numpy array.

    Returns
    -------
    Result : statsmodels.tsa.seasonal.DecomposeResult object
        Object containing the decomposition results

    """
    assert isinstance(series, np.ndarray), "Series is not a numpy array."

    Result = STL(series,
                 period=24,
                 seasonal=7,
                 trend=None,
                 low_pass=None,
                 seasonal_deg=0,
                 trend_deg=0,
                 low_pass_deg=0,
                 robust=False,
                 seasonal_jump=1,
                 trend_jump=1,
                 low_pass_jump=1).fit()

    if test == False:
        _plot_decomposition(Result,
                          title,
                          savepath = savepath,
                          savename = savename,
                          ylabel = "Value",
                          xlabel  = "Date",
                          dates = False,
                          )

    return Result

def test_STL():
    """
    Test STL_decomposition function. Test passes with proper arguments and
    raises an AssertionError if the input time series is not numpy array.

    Returns
    -------
    None.

    """
    res = STL_decomposition(setup_np(),'Test title', test=True)
    assert res is not None
    assert res.observed.all() is not None
    assert res.trend.all() is not None
    assert res.seasonal.all() is not None
    assert res.resid.all() is not None

    test_argument = setup_pd()
    # Store information about raised ValueError in exc_info
    with pytest.raises(AssertionError) as exc_info:
        STL_decomposition(test_argument,'Test title', test=True)
    expected_error_msg = "Series is not a numpy array."
    # Check if the raised ValueError contains the correct message
    assert exc_info.match(expected_error_msg)

max_stars_repo_name: rantahar/tscfat

max_issues_repo_name: rantahar/tscfat

max_forks_repo_name: rantahar/tscfat

reason: import numpy,from statsmodels


text:
import os
import numpy
from subprocess import Popen, PIPE
import pymesh
import tempfile

from default_config.global_vars import apbs_bin, pdb2pqr_bin, multivalue_bin
import random

"""
Modified from:
computeAPBS.py: Wrapper function to compute the Poisson Boltzmann electrostatics for a surface using APBS.
Pablo Gainza - LPDI STI EPFL 2019
"""

def computeAPBS(vertices, pdb_file, tmp_file_base):
    """
        Calls APBS, pdb2pqr, and multivalue and returns the charges per vertex
    """
    pdb2pqr = pdb2pqr_bin + " --ff=parse --whitespace --noopt --apbs-input %s %s"# + tempfile.mktemp()
    make_pqr = pdb2pqr % (pdb_file, tmp_file_base)
    os.system(make_pqr)

    apbs = apbs_bin + " %s"
    make_apbs = apbs % (tmp_file_base+".in")
    os.system(make_apbs)

    vertfile = open(tmp_file_base + ".csv", "w")
    for vert in vertices:
        vertfile.write("{},{},{}\n".format(vert[0], vert[1], vert[2]))
    vertfile.close()

    multivalue = multivalue_bin + " %s %s %s"
    make_multivalue = multivalue % (tmp_file_base+".csv", tmp_file_base+".dx", tmp_file_base+"_out.csv")
    os.system(make_multivalue)

    # Read the charge file
    chargefile = open(tmp_file_base + "_out.csv")
    charges = numpy.array([0.0] * len(vertices))
    for ix, line in enumerate(chargefile.readlines()):
        charges[ix] = float(line.split(",")[3])
    os.system("rm " + tmp_file_base + "*")
    os.system("rm io.mc")

    return charges

"""  ORIGINAL FUNCTION
'''
computeAPBS.py: Wrapper function to compute the Poisson Boltzmann electrostatics for a surface using APBS.
Pablo Gainza - LPDI STI EPFL 2019
This file is part of MaSIF.
Released under an Apache License 2.0
'''

def computeAPBS(vertices, pdb_file, tmp_file_base = tempfile.mktemp()):

        #Calls APBS, pdb2pqr, and multivalue and returns the charges per vertex

    #fields = tmp_file_base.split("/")[0:-1]
    #directory = "/".join(fields) + "/"
    fields = tmp_file_base
    directory = str(fields) + "/"
    filename_base = tmp_file_base.split("/")[-1]
    pdbname = pdb_file.split("/")[-1]
    args = [
        pdb2pqr_bin,
        "--ff=parse",
        "--whitespace",
        "--noopt",
        "--apbs-input",
        pdbname,
        filename_base,
    ]
    p2 = Popen(args, stdout=PIPE, stderr=PIPE, cwd=directory)
    stdout, stderr = p2.communicate()

    args = [apbs_bin, filename_base + ".in"]
    p2 = Popen(args, stdout=PIPE, stderr=PIPE, cwd=directory)
    stdout, stderr = p2.communicate()

    vertfile = open(directory + "/" + filename_base + ".csv", "w")
    for vert in vertices:
        vertfile.write("{},{},{}\n".format(vert[0], vert[1], vert[2]))
    vertfile.close()

    args = [
        multivalue_bin,
        filename_base + ".csv",
        filename_base + ".dx",
        filename_base + "_out.csv",
    ]
    p2 = Popen(args, stdout=PIPE, stderr=PIPE, cwd=directory)
    stdout, stderr = p2.communicate()

    # Read the charge file
    chargefile = open(tmp_file_base + "_out.csv")
    charges = numpy.array([0.0] * len(vertices))
    for ix, line in enumerate(chargefile.readlines()):
        charges[ix] = float(line.split(",")[3])

    return charges
"""

max_stars_repo_name: ricomnl/deepdock

max_issues_repo_name: ricomnl/deepdock

max_forks_repo_name: ricomnl/deepdock

reason: import numpy


text:
# Copyright 2020 Huawei Technologies Co., Ltd
#
# 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
import pytest

import mindspore.context as context
from mindspore.common.tensor import Tensor
from mindspore.ops import operations as P

@pytest.mark.level0
@pytest.mark.platform_x86_gpu_training
@pytest.mark.env_onecard
def test_broadcast():
    context.set_context(mode=context.GRAPH_MODE, device_target='GPU')

    shape = (3, 4, 5, 6)
    x_np = np.random.rand(3, 1, 5, 1).astype(np.float32)
    output = P.BroadcastTo(shape)(Tensor(x_np))
    expect = np.broadcast_to(x_np, shape)
    assert np.allclose(output.asnumpy(), expect)

    x1_np = np.random.rand(3, 1, 5, 1).astype(np.float16)
    output = P.BroadcastTo(shape)(Tensor(x1_np))
    expect = np.broadcast_to(x1_np, shape)
    assert np.allclose(output.asnumpy(), expect)

    shape = (2, 3, 4, 5)
    x1_np = np.random.rand(4, 5).astype(np.float32)
    output = P.BroadcastTo(shape)(Tensor(x1_np))
    expect = np.broadcast_to(x1_np, shape)
    assert np.allclose(output.asnumpy(), expect)

@pytest.mark.level0
@pytest.mark.platform_x86_gpu_training
@pytest.mark.env_onecard
def test_broadcast_dyn_init():
    """
    Test running the op with -1's in the init shape to support varied inputs.
    """
    context.set_context(mode=context.GRAPH_MODE, device_target='GPU')

    ms_shape = (-1, 4, 5, 6)
    np_shape = (3, 4, 5, 6)
    x_np = np.random.rand(3, 1, 5, 1).astype(np.float32)
    output = P.BroadcastTo(ms_shape)(Tensor(x_np))
    expect = np.broadcast_to(x_np, np_shape)
    assert np.allclose(output.asnumpy(), expect)

    x1_np = np.random.rand(3, 1, 5, 1).astype(np.float16)
    output = P.BroadcastTo(ms_shape)(Tensor(x1_np))
    expect = np.broadcast_to(x1_np, np_shape)
    assert np.allclose(output.asnumpy(), expect)

    ms_shape = (2, 3, -1, 5)
    np_shape = (2, 3, 4, 5)
    x1_np = np.random.rand(4, 5).astype(np.float32)
    output = P.BroadcastTo(ms_shape)(Tensor(x1_np))
    expect = np.broadcast_to(x1_np, np_shape)
    assert np.allclose(output.asnumpy(), expect)

@pytest.mark.level0
@pytest.mark.platform_x86_gpu_training
@pytest.mark.env_onecard
def test_broadcast_dyn_invalid_init():
    """
    Test running the op with -1's in the init shape in incorrect positions.
    Expected to fail.
    """
    context.set_context(mode=context.GRAPH_MODE, device_target='GPU')
    ms_shape = (2, -1, 4, 5)
    x_np = np.random.rand(4, 5).astype(np.float32)
    with pytest.raises(ValueError):
        P.BroadcastTo(ms_shape)(Tensor(x_np))

max_stars_repo_name: GuoSuiming/mindspore


text:
"""
Copyright 2017 Steven Diamond

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 cvxpy.settings as s
from cvxpy.constraints.bool_constr import BoolConstr

class IntConstr(BoolConstr):
    """
    An integer constraint:
        X_{ij} in Z for all i,j.

    Attributes:
        noncvx_var: A variable constrained to be elementwise integral.
        lin_op: The linear operator equal to the noncvx_var.
    """
    CONSTR_TYPE = s.INT_IDS

    def __str__(self):
        return "IntConstr(%s)" % self.lin_op

max_stars_repo_name: quantopian/cvxpy

max_issues_repo_name: quantopian/cvxpy

max_forks_repo_name: quantopian/cvxpy

reason: import cvxpy,from cvxpy


text:
# Copyright (C) 2020-2021 Intel Corporation
# SPDX-License-Identifier: Apache-2.0

"""FederatedDataset module."""

import numpy as np

from openfl.utilities.data_splitters import EqualNumPyDataSplitter
from openfl.utilities.data_splitters import NumPyDataSplitter
from .loader_pt import PyTorchDataLoader

class FederatedDataSet(PyTorchDataLoader):
    """
    Data Loader for in memory Numpy data.

    Args:
        X_train: np.array
            Training Features
        y_train: np.array
            Training labels
        X_val: np.array
            Validation features
        y_val: np.array
            Validation labels
        batch_size : int
            The batch size for the data loader
        num_classes : int
            The number of classes the model will be trained on
        **kwargs: Additional arguments to pass to the function

    """

    train_splitter: NumPyDataSplitter
    valid_splitter: NumPyDataSplitter

    def __init__(self, X_train, y_train, X_valid, y_valid,
                 batch_size=1, num_classes=None, train_splitter=None, valid_splitter=None):
        """
        Initialize.

        Args:
            X_train: np.array
                Training Features
            y_train: np.array
                Training labels
            X_val: np.array
                Validation features
            y_val: np.array
                Validation labels
            batch_size : int
                The batch size for the data loader
            num_classes : int
                The number of classes the model will be trained on
            train_splitter: NumPyDataSplitter
                Data splitter for train dataset.
            valid_splitter: NumPyDataSplitter
                Data splitter for validation dataset.
            **kwargs: Additional arguments to pass to the function

        """
        super().__init__(batch_size)

        self.X_train = X_train
        self.y_train = y_train
        self.X_valid = X_valid
        self.y_valid = y_valid

        if num_classes is None:
            num_classes = np.unique(self.y_train).shape[0]
            print(f'Inferred {num_classes} classes from the provided labels...')
        self.num_classes = num_classes
        self.train_splitter = self._get_splitter_or_default(train_splitter)
        self.valid_splitter = self._get_splitter_or_default(valid_splitter)

    @staticmethod
    def _get_splitter_or_default(value):
        if value is None:
            return EqualNumPyDataSplitter()
        if isinstance(value, NumPyDataSplitter):
            return value
        else:
            raise NotImplementedError(f'Data splitter {value} is not supported')

    def split(self, num_collaborators):
        """Create a Federated Dataset for each of the collaborators.

        Args:
            num_collaborators: int
                Collaborators to split the dataset between
            shuffle: boolean
                Should the dataset be randomized?
            equally: boolean
                Should each collaborator get the same amount of data?

        Returns:
            list[FederatedDataSets]
                A dataset slice for each collaborator
        """
        train_idx = self.train_splitter.split(self.y_train, num_collaborators)
        valid_idx = self.valid_splitter.split(self.y_valid, num_collaborators)

        return [
            FederatedDataSet(
                self.X_train[train_idx[i]],
                self.y_train[train_idx[i]],
                self.X_valid[valid_idx[i]],
                self.y_valid[valid_idx[i]],
                batch_size=self.batch_size,
                num_classes=self.num_classes,
                train_splitter=self.train_splitter,
                valid_splitter=self.valid_splitter
            ) for i in range(num_collaborators)
        ]

max_stars_repo_name: brandon-edwards/openfl

max_issues_repo_name: brandon-edwards/openfl

max_forks_repo_name: brandon-edwards/openfl


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import wave
import numpy as np
import pylab

def flatten(l):
    return [value for sublist in l for value in sublist]

def chunks(l, n):
    """Yield successive n-sized chunks from l."""
    for i in xrange(0, len(l), n):
        yield l[i:i+n]

class Normalize(object):
    """docstring for Normalize"""
    def __init__(self, db_level = -1.0):
        super(Normalize, self).__init__()

        self.maximum = -999999
        self.minimum = 999999

        self.db_level = db_level
        self.linear_level = self._db_to_linear(self._bounded_level(db_level))

        # self.current_mean = None
        # self.values_n = 0

        self.av = 0.0
        self.n = 0.0

    def _calc_mean(self, chunk):
        """
            This formula cal help avoid broblem with variable overflow
            (nx + m) can be very very large number.
            But after decomposite this work greatly

             nx + m           n         m
            --------  =  x ------- + -------
              n + 1         n + 1     n + 1

        """

        # for sample in chunk:
        #     self.n += 1.0
        #     self.av = (self.n - 1) / self.n * self.av + sample / self.n

        #   nx + sum(m)            n            sum(m)
        #  -------------- = x ------------ + ------------
        #    n + len(m)        n + len(m)      n + len(m)

        self.n += float(len(chunk))
        self.av = (self.n - len(chunk)) / self.n * self.av + sum(chunk) / self.n

        return self.av

    def _bounded_level(self, level):
        level_min = -145.0
        level_max = 0.0

        current_level = level

        if current_level > level_max:
            current_level = level_max

        if current_level < level_min:
            current_level = level_min

        return current_level

    def _db_to_linear(self, db_level):
        return pow(10.0, float(db_level) / 20.0)

    # def _normalize(self, chunk):
    #     return map(lambda value: value / float(np.iinfo(np.int16).max), chunk)

    # def _de_normalize(self, chunk):
    #     return map(lambda value: int(value * np.iinfo(np.int16).max), chunk)

    def process(self, chunk):

        current_maximim, current_minimum = self._get_max_min(chunk)

        self.maximum = max(self.maximum, current_maximim)
        self.minimum = min(self.minimum, current_minimum)

        offset = self._calc_mean(chunk)

        ratio = self.linear_level * float(np.iinfo(np.int16).max)

        extent = max(abs(self.maximum), abs(self.minimum))

        mult = ratio / float(extent)

        processed = map(lambda sample: int((sample - offset) * mult), chunk)

        return processed

        # return self._de_normalize(processed)

    def _get_max_min(self, chunk):
        maximum = max(chunk)
        minimum = min(chunk)

        return (maximum, minimum)

def main():

    # normalize = Normalize()

    # print normalize._calc_mean([1, 2, 3, 4]) # 2.5
    # print normalize._calc_mean([5, 6, 7, 8]) # 4.5

    # return

    file_name = '1453874060289.wav'

    # data = numpy.memmap(file_name, dtype='h', mode='r')

    with open(file_name, 'rb') as f:
        content = np.fromstring(f.read(), np.int16)

    # print content
    # print content.tostring()

    # pylab.plot(content)
    # pylab.plot(content)

    # maximum = np.max(content)
    # minimum = np.min(content)

    # processed = map(lambda sample: sample * 1.5, content)

    normalize = Normalize()

    # print normalize._calc_mean([1, 2, 3, 4])
    # print normalize._calc_mean([1, 2, 3, 4])

    # return

    # print normalize.process([40, 40])

    result = flatten(map(lambda chunk: normalize.process(chunk), chunks(content, 320)))

    # # print len(content)
    # # print len(result)

    # # result = map(lambda sample: sample * 1.5, content)

    # print result

    # pylab.plot(result)
    # pylab.show()

    q = np.array(result).astype(np.int16)

    with open('res.pcm', 'wb') as f:
        # f.write((np.array(result, dtype='h')).tostring())
        f.write(q.tostring())

    # processed = map(lambda sample: sample * 1.5, data)

    # print processed

if __name__ == '__main__':
    main()


text:
import numpy as np
import tensorflow as tf
import skimage
import skimage.segmentation
import copy
import sklearn.metrics
from sklearn import linear_model

from .mc_drop import get_entropy
from .latent_optimisation import (
    get_latent_representation,
    get_in_class_fiducial,
    get_clue_counterfactual,
)
from .ig import (
    ig_entropy,
    ig_latent_curve,
)

# OUR PROPOSED METHOD
def get_importances(x: tf.Tensor, model, encoder, decoder, known_class=None,
                    max_num_trials=3, entropy_reduction=0.90, jacobian_splits=1):
    """IG attributions with integration path in latent space.
    If getting OOM error, increase the jacobian_splits parameter.
    Make sure that the image size is divisible by jacobian_splits."""

    if x.shape[0] % jacobian_splits:
        raise ValueError(f"image size ({x.shape[0]}) is not divisible by "
                         f"`jacobian_splits` ({jacobian_splits}).")

    z = get_latent_representation(x, encoder, decoder)

    # If the fiducial doesn't reach low enough entropy,
    # run the optimisation algorithm again. Do it at most
    # `max_num_trials` times.
    z_fid = None
    x_ent = get_entropy(x[None, ...], model)
    for _ in range(max_num_trials):
        z_fid = get_in_class_fiducial(x, model, encoder, decoder, known_class=known_class)
        fid_ent = get_entropy(decoder(z_fid), model)
        if (fid_ent <= (1 - entropy_reduction) * x_ent):
            break

    return ig_latent_curve(z_fid, z, x, decoder, model, jacobian_splits=jacobian_splits)

def get_importances_vanilla(x: tf.Tensor, model):
    # Set fiducial to a black image
    fiducial = tf.zeros_like(x)

    return ig_entropy(fiducial, x, model)

def get_importances_CLUE(x: tf.Tensor, model, encoder, decoder,
                         max_num_trials=5, entropy_reduction=0.90):
    # If the fiducial doesn't reach low enough entropy,
    # run the optimisation algorithm again. Do it at most
    # `max_num_trials` times.
    z_fid = None
    x_ent = get_entropy(x[None, ...], model)
    for _ in range(max_num_trials):
        z_fid = get_clue_counterfactual(x, model, encoder, decoder)
        fid_ent = get_entropy(decoder(z_fid), model)
        if (fid_ent <= (1 - entropy_reduction) * x_ent):
            break

    return tf.reduce_sum(x - decoder(z_fid)[0], axis=2).numpy()

def get_importances_LIME(x: tf.Tensor, model, mc_sample=100, num_perturb=None,
                         alpha=None, bs=32, kernel_size=1, max_dist=5, ratio=0.2):
    # Generate superpixels
    channels = x.shape[2]
    if channels == 1:
        x_rgb = skimage.color.gray2rgb(x[:, :, 0]).astype('double')
    else:
        x_rgb = x.astype('double')
    superpixels = skimage.segmentation.quickshift(x_rgb, kernel_size=kernel_size,
                                                  max_dist=max_dist, ratio=ratio)
    num_superpixels = np.unique(superpixels).shape[0]

    if num_perturb is None:
        num_perturb = 10 * num_superpixels + 100

    # Generate perturbations
    perturbations = np.random.binomial(1, 0.8, size=(num_perturb,
                                                     num_superpixels))

    # Generate predictions for perturbations
    perturbed_images = []
    for pert in perturbations:
        perturbed_images.append(
            perturb_image(x_rgb, pert, superpixels)[None, :, :, :channels]
        )
    perturbed_images = tf.concat(perturbed_images, axis=0)
    predictions = get_entropy(perturbed_images, model, mc_sample, bs=bs)

    # Regression weights
    reference = np.ones(num_superpixels)[None, :]
    w = sklearn.metrics.pairwise_distances(perturbations,
                                           reference,
                                           metric='cosine').ravel()
    w = np.exp(-w)

    # Run Lasso regression
    if alpha is not None:
        lf = linear_model.Lasso(alpha=alpha)
    else:
        lf = linear_model.LassoCV(cv=5)
    lf.fit(X=perturbations, y=predictions, sample_weight=w)

    # Get importances in original image shape
    return lf.coef_[superpixels]

def get_importances_SHAP(x: tf.Tensor, model, x_train, mc_sample=100,
                         num_perturb=None, alpha=None, bs=32):
    """Implementation of the kernelShap attribution method."""
    if num_perturb is None:
        num_perturb = 3 * x.shape[0] ** 2 + 2 ** 11
    # Generate perturbations
    perc = 0.05
    feat_c = x.shape[0]*x.shape[1]
    perturbations = np.repeat(x[None, :], num_perturb, axis=0)
    mask = np.ones(shape=perturbations.shape[:3])

    # Sample from marginal x_train
    deac_0 = np.repeat(range(num_perturb), int(feat_c*perc))
    deac_1 = np.random.choice(x.shape[0], deac_0.shape[0], replace=True)
    deac_2 = np.random.choice(x.shape[0], deac_0.shape[0], replace=True)
    idx_train = np.random.choice(x_train.shape[0],
                                 deac_0.shape[0], replace=True)

    perturbations[deac_0, deac_1, deac_2] = \
        x_train[idx_train, deac_1, deac_2].copy()
    mask[deac_0, deac_1, deac_2] = 0

    # Score the perturbations
    predictions = get_entropy(perturbations, model, mc_sample, bs=bs)

    # Run regression
    if alpha is not None:
        lf = linear_model.Lasso(alpha=alpha)
    else:
        lf = linear_model.LassoLarsIC(criterion='aic')
    lf.fit(X=mask.reshape(num_perturb, -1),
           y=predictions)

    # Get importances in original image shape
    return lf.coef_.reshape(x.shape[:2])

def perturb_image(img, perturbation, segments):
    active_pixels = np.where(perturbation == 1)[0]
    mask = np.zeros(segments.shape)
    for active in active_pixels:
        mask[segments == active] = 1

    perturbed_image = copy.deepcopy(img)
    perturbed_image = perturbed_image * mask[:, :, None]

    return perturbed_image

max_stars_repo_name: Featurespace/uncertainty-attribution

max_issues_repo_name: Featurespace/uncertainty-attribution

max_forks_repo_name: Featurespace/uncertainty-attribution


text:
#Copyright (c) 2008 Erik Tollerud (etolleru@uci.edu)
import numpy as np

class Pca(object):
    """
    A basic class for Principal Component Analysis (PCA).

    p is the number of dimensions, while N is the number of data points
    """
    _colors=('r','g','b','c','y','m','k') #defaults

    def __calc(self):
        A = self.A
        M=A-np.mean(A,axis=0)
        N=M/np.std(M,axis=0)

        self.M = M
        self.N = N
        self._eig = None

    def __init__(self,data,names=None):
        """
        p X N matrix input
        """
        A = np.array(data).T
        n,p = A.shape
        self.n,self.p = n,p
        if p > n:
            from warnings import warn
            warn('p > n - intentional?', RuntimeWarning)
        self.A = A
        self._origA=A.copy()

        self.__calc()

        self._colors= np.tile(self._colors,int((p-1)/len(self._colors))+1)[:p]
        if names is not None and len(names) != p:
            raise ValueError('names must match data dimension')
        self.names = None if names is None else tuple([str(x) for x in names])

    def getCovarianceMatrix(self):
        """
        returns the covariance matrix for the dataset
        """
        return np.cov(self.N.T)

    def getEigensystem(self):
        """
        returns a tuple of (eigenvalues,eigenvectors) for the data set.
        """
        if self._eig is None:
            res = np.linalg.eig(self.getCovarianceMatrix())
            sorti=np.argsort(res[0])[::-1]
            res=(res[0][sorti],res[1][:,sorti])
            self._eig=res
        return self._eig

    def getEigenvalues(self):
        return self.getEigensystem()[0]

    def getEigenvectors(self):
        return self.getEigensystem()[1]

    def getEnergies(self):
        """
        "energies" are just normalized eigenvectors
        """
        v=self.getEigenvalues()
        return v/np.sum(v)

    def plot2d(self,ix=0,iy=1,clf=True):
        """
        Generates a 2-dimensional plot of the data set and principle components
        using matplotlib.

        ix specifies which p-dimension to put on the x-axis of the plot
        and iy specifies which to put on the y-axis (0-indexed)
        """
        import matplotlib.pyplot as plt
        x,y=self.N[:,ix],self.N[:,iy]
        if clf:
            plt.clf()
        plt.scatter(x,y)
        vals,evs=self.getEigensystem()
        #evx,evy=evs[:,ix],evs[:,iy]
        xl,xu=plt.xlim()
        yl,yu=plt.ylim()
        dx,dy=(xu-xl),(yu-yl)
        for val,vec,c in zip(vals,evs.T,self._colors):
            plt.arrow(0,0,val*vec[ix],val*vec[iy],head_width=0.05*(dx*dy/4)**0.5,fc=c,ec=c)
        #plt.arrow(0,0,vals[ix]*evs[ix,ix],vals[ix]*evs[iy,ix],head_width=0.05*(dx*dy/4)**0.5,fc='g',ec='g')
        #plt.arrow(0,0,vals[iy]*evs[ix,iy],vals[iy]*evs[iy,iy],head_width=0.05*(dx*dy/4)**0.5,fc='r',ec='r')
        if self.names is not None:
            plt.xlabel('$'+self.names[ix]+'/\\sigma$')
            plt.ylabel('$'+self.names[iy]+'/\\sigma$')

    def plot3d(self,ix=0,iy=1,iz=2,clf=True):
        """
        Generates a 3-dimensional plot of the data set and principle components
        using mayavi.

        ix, iy, and iz specify which of the input p-dimensions to place on each of
        the x,y,z axes, respectively (0-indexed).
        """
        import enthought.mayavi.mlab as M
        if clf:
            M.clf()
        z3=np.zeros(3)
        v=(self.getEigenvectors()*self.getEigenvalues())
        M.quiver3d(z3,z3,z3,v[ix],v[iy],v[iz],scale_factor=5)
        M.points3d(self.N[:,ix],self.N[:,iy],self.N[:,iz],scale_factor=0.3)
        if self.names:
            M.axes(xlabel=self.names[ix]+'/sigma',ylabel=self.names[iy]+'/sigma',zlabel=self.names[iz]+'/sigma')
        else:
            M.axes()

    def sigclip(self,sigs):
        """
        clips out all data points that are more than a certain number
        of standard deviations from the mean.

        sigs can be either a single value or a length-p sequence that
        specifies the number of standard deviations along each of the
        p dimensions.
        """
        if np.isscalar(sigs):
            sigs=sigs*np.ones(self.N.shape[1])
        sigs = sigs*np.std(self.N,axis=1)
        n = self.N.shape[0]
        m = np.all(np.abs(self.N) < sigs,axis=1)
        self.A=self.A[m]
        self.__calc()
        return n-sum(m)

    def reset(self):
        self.A = self._origA.copy()
        self.__calc()

    def project(self,vals=None,enthresh=None,nPCs=None,cumen=None):
        """
        projects the normalized values onto the components

        enthresh, nPCs, and cumen determine how many PCs to use

        if vals is None, the normalized data vectors are the values to project.
        Otherwise, it should be convertable to a p x N array

        returns n,p(>threshold) dimension array
        """
        nonnones = sum([e is not None for e in (enthresh, nPCs, cumen)])
        if nonnones == 0:
            m = slice(None)
        elif nonnones > 1:
            raise ValueError("cannot specify more than one threshold")
        else:
            if enthresh is not None:
                m = self.energies() > enthresh
            elif nPCs is not None:
                m = slice(None,nPCs)
            elif cumen is not None:
                m = np.cumsum(self.energies()) <  cumen
            else:
                raise RuntimeError('Should be unreachable')

        if vals is None:
            vals = self.N.T
        else:
            vals = np.array(vals,copy=False)
            if self.N.T.shape[0] != vals.shape[0]:
                raise ValueError("shape for vals does not match")
        proj = np.matrix(self.getEigenvectors()).T*vals
        return proj[m].T

    def deproject(self,A,normed=True):
        """
        input is an n X q array, where q <= p

        output is p X n
        """
        A=np.atleast_2d(A)
        n,q = A.shape
        p = self.A.shape[1]
        if q > p :
            raise ValueError("q > p")

        evinv=np.linalg.inv(np.matrix(self.getEigenvectors()).T)

        zs = np.zeros((n,p))
        zs[:,:q]=A

        proj = evinv*zs.T

        if normed:
            return np.array(proj.T).T
        else:
            mns=np.mean(self.A,axis=0)
            sds=np.std(self.M,axis=0)
            return (np.array(proj.T)*sds+mns).T

    def subtractPC(self,pc,vals=None):
        """
        pc can be a scalar or any sequence of pc indecies

        if vals is None, the source data is self.A, else whatever is in vals
        (which must be p x m)
        """
        if vals is None:
            vals = self.A
        else:
            vals = vals.T
            if vals.shape[1]!= self.A.shape[1]:
                raise ValueError("vals do not have the correct number of components")

        pcs=self.project()
        zpcs=np.zeros_like(pcs)
        zpcs[:,pc]=pcs[:,pc]
        upc=self.deproject(zpcs,False)

        A = vals.T-upc
        B = A.T*np.std(self.M,axis=0)
        return B+np.mean(self.A,axis=0)

max_stars_repo_name: EkremBayar/bayar

max_issues_repo_name: EkremBayar/bayar

max_forks_repo_name: EkremBayar/bayar


text:
#%% Importing libraries

import os
import glob
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib import rc
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import erfc as erfc
from scipy.optimize import curve_fit

#%% Setting a "LaTeX-like" plotting style

plt.style.use('default')
rc('font', **{'family': 'DejaVu Sans', 'serif': ['Computer Modern']})
rc('text', usetex=True)

#%% For loop to reconstruct beam profile

# Initializing data sets

sets = glob.glob('Data*.txt')
sets.sort()
j    = 1 # custom index for saving plots without overwriting them

# Initializing for loop

for i in sets:
    Data = np.loadtxt(i, comments = '#')

# Defining variables and functions

    z    = Data[:, 0] # [mm]
    x    = Data[:, 1] # [um]
    P    = Data[:, 2] # [W]
    dx   = Data[:, 3] # [um]
    dP   = Data[:, 4] # [W]
    P0   = np.amin(P)
    Pmax = np.amax(P)

    def fit(x, x0, w):
        return P0+Pmax/2*erfc(np.sqrt(2)*(x-x0)/w) # [W]

    xfit = np.linspace(np.amin(x), np.amax(x), 300) # [um]
    pi   = [(np.amax(x)-np.amin(x))/2, 1800]

# Fitting

    popt, pcov = curve_fit(fit, x, P, pi)
    sigma = np.sqrt([pcov[0, 0], pcov[1, 1]])
    fit_values = fit(xfit, popt[0], popt[1])

    w  = popt[1]  # [um]
    dw = sigma[1] # [um]

# Calculating derivative

    def prime(f, a, h):
        return (f(a+h)-f(a-h))/(2*h)

    def fun(x):
        return fit(x, popt[0], w)

    g = -prime(fun, x, 0.001)*1000

# Plotting 2D gaussian profile

    fig = plt.figure()
    grid = plt.GridSpec(2, 8)

    ax1 = fig.add_subplot(grid[1, :4])
    ax2 = fig.add_subplot(grid[1, 4:5])

    ax1.plot(x, g, '-r')
    ax1.set_xlim(2000, 11000)
    ax1.set_ylim(0, 8)
    ax1.set_xlabel('$x$ ($\mu$m)')
    ax1.set_ylabel('$\partial P$ (u.a.)')

    ax2.axhline(z[j-1], color = 'r')
    ax2.set_xlim(0, 1)
    ax2.set_ylim(0, 450)
    ax2.set_ylabel('$z$ (mm)')
    ax2.set_xticks([])
    ax2.yaxis.set_label_position('right')
    ax2.yaxis.set_ticks_position('right')
    plt.savefig('./frames/%d.png' % j, dpi = 300)

# Manually incrementing j

    j = j+1

#%% Converting plots into a gif

gif_name = 'Profile_Evolution'
frames = glob.glob('./frames/*.png')
list.sort(frames, key = lambda x: int(x.split('./frames/')[1].split('.png')[0]))

with open('frames.txt', 'w') as file:
    for i in frames:
        file.write("%s\n" % i)

os.system('convert -delay 60 @frames.txt {}.gif'.format(gif_name))

max_stars_repo_path: Beam Profile (2D).py

max_stars_repo_name: ConteFm/LaBPro

max_issues_repo_path: Beam Profile (2D).py

max_issues_repo_name: ConteFm/LaBPro

max_forks_repo_path: Beam Profile (2D).py

max_forks_repo_name: ConteFm/LaBPro


text:
from builtins import object

import pytest
import numpy as np

from sporco.fista import ccmod
import sporco.linalg as sl
import sporco.cnvrep as cr

class TestSet01(object):

    def setup_method(self, method):
        np.random.seed(12345)

    def test_01(self):
        N = 32
        M = 4
        Nd = 5
        D0 = cr.normalise(cr.zeromean(
            np.random.randn(Nd, Nd, M), (Nd, Nd, M), dimN=2), dimN=2)
        X = np.zeros((N, N, M))
        xr = np.random.randn(N, N, M)
        xp = np.abs(xr) > 3
        X[xp] = np.random.randn(X[xp].size)
        S = np.sum(sl.ifftn(sl.fftn(D0, (N, N), (0,1)) *
                   sl.fftn(X, None, (0,1)), None, (0,1)).real, axis=2)
        L = 3e1
        opt = ccmod.ConvCnstrMOD.Options({'Verbose': False,
                    'MaxMainIter': 3000,
                    'ZeroMean': True, 'RelStopTol': 1e-6, 'L': L,
                    'BackTrack': {'Enabled': False}})
        Xr = X.reshape(X.shape[0:2] + (1,1,) + X.shape[2:])
        Sr = S.reshape(S.shape + (1,))
        c = ccmod.ConvCnstrMOD(Xr, Sr, D0.shape, opt)
        c.solve()
        D1 = cr.bcrop(c.X, D0.shape).squeeze()
        assert(sl.rrs(D0, D1) < 1e-4)
        assert(np.array(c.getitstat().Rsdl)[-1] < 1e-5)

    def test_02(self):
        N = 32
        M = 4
        Nd = 5
        D0 = cr.normalise(cr.zeromean(
            np.random.randn(Nd, Nd, M), (Nd, Nd, M), dimN=2), dimN=2)
        X = np.zeros((N, N, M))
        xr = np.random.randn(N, N, M)
        xp = np.abs(xr) > 3
        X[xp] = np.random.randn(X[xp].size)
        S = np.sum(sl.ifftn(sl.fftn(D0, (N, N), (0,1)) *
                   sl.fftn(X, None, (0,1)), None, (0,1)).real, axis=2)
        L = 1e1
        opt = ccmod.ConvCnstrMOD.Options({'Verbose': False,
                    'MaxMainIter': 3000,
                    'ZeroMean': True, 'RelStopTol': 1e-6, 'L': L,
                    'BackTrack': {'Enabled': True}})
        Xr = X.reshape(X.shape[0:2] + (1,1,) + X.shape[2:])
        Sr = S.reshape(S.shape + (1,))
        c = ccmod.ConvCnstrMOD(Xr, Sr, D0.shape, opt)
        c.solve()
        D1 = cr.bcrop(c.X, D0.shape).squeeze()
        assert(sl.rrs(D0, D1) < 1e-4)
        assert(np.array(c.getitstat().Rsdl)[-1] < 1e-5)

    def test_03(self):
        N = 64
        M = 4
        Nd = 8
        D0 = cr.normalise(cr.zeromean(
            np.random.randn(Nd, Nd, M), (Nd, Nd, M), dimN=2), dimN=2)
        X = np.zeros((N, N, M))
        xr = np.random.randn(N, N, M)
        xp = np.abs(xr) > 3
        X[xp] = np.random.randn(X[xp].size)
        S = np.sum(sl.ifftn(sl.fftn(D0, (N, N), (0,1)) *
                   sl.fftn(X, None, (0,1)), None, (0,1)).real, axis=2)
        L = 1e1
        opt = ccmod.ConvCnstrMOD.Options({'Verbose': False,
                    'MaxMainIter': 3000,
                    'ZeroMean': True, 'RelStopTol': 1e-6, 'L': L,
                    'BackTrack': {'Enabled': True}})
        Xr = X.reshape(X.shape[0:2] + (1,1,) + X.shape[2:])
        Sr = S.reshape(S.shape + (1,))
        c = ccmod.ConvCnstrMOD(Xr, Sr, D0.shape, opt)
        c.solve()
        D1 = cr.bcrop(c.X, D0.shape).squeeze()
        assert(sl.rrs(D0,D1) < 1e-4)
        assert(np.array(c.getitstat().Rsdl)[-1] < 1e-5)

    def test_04(self):
        N = 16
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, 1, M)
        S = np.random.randn(N, N, 1)
        try:
            c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, M))
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_05(self):
        N = 16
        M = 8
        X = np.random.randn(N, N, 1, 1, M)
        S = np.random.randn(N, N, 1)
        try:
            c = ccmod.ConvCnstrMOD(X, S, ((4, 4, 4),(8, 8, 4)))
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_06(self):
        N = 16
        M = 4
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, 1, M)
        S = np.random.randn(N, N, Nc)
        L = 2e3
        try:
            opt = ccmod.ConvCnstrMODOptions({'Verbose': False,
                            'MaxMainIter': 100, 'L' : L})
            c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, 1, M), opt=opt, dimK=0)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)
        assert(np.array(c.getitstat().Rsdl)[-1] < 5e-3)

    def test_07(self):
        N = 16
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, 1, M)
        S = np.random.randn(N, N, 1)
        dt = np.float32
        opt = ccmod.ConvCnstrMODOptions(
            {'Verbose': False, 'MaxMainIter': 20,
             'BackTrack': {'Enabled': True},
             'DataType': dt})
        c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, M), opt=opt)
        c.solve()
        assert(c.X.dtype == dt)
        assert(c.Xf.dtype == sl.complex_dtype(dt))
        assert(c.Yf.dtype == sl.complex_dtype(dt))

    def test_08(self):
        N = 16
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, 1, M)
        S = np.random.randn(N, N, 1)
        dt = np.float64
        opt = ccmod.ConvCnstrMODOptions(
            {'Verbose': False, 'MaxMainIter': 20,
             'BackTrack': {'Enabled': True},
             'DataType': dt})
        c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, M), opt=opt)
        c.solve()
        assert(c.X.dtype == dt)
        assert(c.Xf.dtype == sl.complex_dtype(dt))
        assert(c.Yf.dtype == sl.complex_dtype(dt))

    def test_09(self):
        N = 16
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, 1, M)
        S = np.random.randn(N, N)
        try:
            opt = ccmod.ConvCnstrMODOptions(
            {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, 1, M),
                opt=opt, dimK=0)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_10(self):
        N = 16
        K = 3
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, K, M)
        S = np.random.randn(N, N, K)
        try:
            opt = ccmod.ConvCnstrMODOptions(
            {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, 1, M),
                opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_11(self):
        N = 16
        M = 4
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, 1, M)
        S = np.random.randn(N, N, Nc)
        try:
            opt = ccmod.ConvCnstrMODOptions(
            {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, 1, M),
                opt=opt, dimK=0)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_12(self):
        N = 16
        M = 4
        K = 2
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, K, M)
        S = np.random.randn(N, N, Nc, K)
        try:
            opt = ccmod.ConvCnstrMODOptions(
            {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, 1, M),
                opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_13(self):
        N = 16
        M = 4
        K = 2
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, K, M)
        S = np.random.randn(N, N, Nc, K)
        try:
            opt = ccmod.ConvCnstrMODOptions(
            {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMOD(X, S, (Nd, Nd, Nc, M),
                opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_14(self):
        N = 16
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, 1, M)
        S = np.random.randn(N, N)
        W = np.array([1.0])
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt, dimK=0)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_15(self):
        N = 16
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, 1, M)
        S = np.random.randn(N, N)
        W = np.random.randn(N, N)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt, dimK=0)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_16(self):
        N = 16
        K = 3
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, K, M)
        S = np.random.randn(N, N, K)
        W = np.random.randn(N, N)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_17(self):
        N = 16
        K = 3
        M = 4
        Nd = 8
        X = np.random.randn(N, N, 1, K, M)
        S = np.random.randn(N, N, K)
        W = np.random.randn(N, N, K)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_18(self):
        N = 16
        M = 4
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, 1, M)
        S = np.random.randn(N, N, Nc)
        W = np.random.randn(N, N, Nc)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt, dimK=0)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_19(self):
        N = 16
        M = 4
        K = 2
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, K, M)
        S = np.random.randn(N, N, Nc, K)
        W = np.random.randn(N, N, Nc, K)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_20(self):
        N = 16
        M = 4
        K = 2
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, K, M)
        S = np.random.randn(N, N, Nc, K)
        W = np.random.randn(N, N, Nc)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_21(self):
        N = 16
        M = 4
        K = 2
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, K, M)
        S = np.random.randn(N, N, Nc, K)
        W = np.random.randn(N, N, 1, K)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, 1, M), opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_22(self):
        N = 16
        M = 4
        K = 2
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, K, M)
        S = np.random.randn(N, N, Nc, K)
        W = np.random.randn(N, N, Nc, K)
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options(
                           {'Verbose': False, 'MaxMainIter': 20})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W,
                    (Nd, Nd, Nc, M), opt=opt)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)

    def test_23(self):
        N = 16
        M = 4
        Nc = 3
        Nd = 8
        X = np.random.randn(N, N, Nc, 1, M)
        S = np.random.randn(N, N, Nc)
        W = np.random.randn(N, N)
        L = 5e3
        try:
            opt = ccmod.ConvCnstrMODMaskDcpl.Options({'Verbose': False,
                            'MaxMainIter': 200, 'L' : L})
            c = ccmod.ConvCnstrMODMaskDcpl(X, S, W, (Nd, Nd, 1, M),
                                             opt=opt, dimK=0)
            c.solve()
        except Exception as e:
            print(e)
            assert(0)
        assert(np.array(c.getitstat().Rsdl)[-1] < 5e-3)

max_stars_repo_name: sophiaas/alphacsc

max_issues_repo_name: sophiaas/alphacsc

max_forks_repo_name: sophiaas/alphacsc


text:
# qmpy/analysis/thermodynamics/convex_hull.pyx

import numpy as np
#cimport numpy as np
#DTYPE = np.float64
#ctypedef np.float64_t DTYPE_t

import numpy.linalg as la
import itertools
import random
import logging

logger = logging.getLogger(__name__)

class Simplex(object):
    """
    A Simplex is a geometric object with N vertices in N dimensions.
    """

    def __init__(self, phases):
        space = set.union(*[p.space for p in phases ])
        cc = [ [p.unit_comp[e] for e in self.space ] for p in phases ]

        self.space = space
        self.comps = np.array(cc)
        self.energies = np.array([ p.energy for p in phases ])

        self.bounds = np.array(range(len(space)))
        self.ref_energies = np.array([ self.energies[i] for i in self.bounds ])
        self.ref_comps = np.array([ self.comps[i] for i in self.bounds ])

        self.find_bounds()

    def contains(self, comp):
        """
        Put the point in barycentric coordinates and check that it is in the
        simplex.
        """
        coord = self.bcc(comp)
        return (all([ a >= 0 for a in coord ]) and sum(coord) <= 1)

    def bcc(self, point):
        r = point - self.comps[self.bounds[0]]
        return np.dot(self.inv, r)

    def get_relative_energy(self, ind):
        return self.energies[ind] - self.comps[ind].dot(
                                           self.ref_energies[ind])

max_stars_repo_name: WalterjhShen/qmpy

max_issues_repo_name: WalterjhShen/qmpy

max_forks_repo_name: WalterjhShen/qmpy


text:
#! /usr/bin/env python

""" Wrapper to run lacosmicx.py
"""

import os
import shutil
import inspect
import glob
import yaml

import lacosmicx
from astropy.io import fits

# -----------------------------------------------------------------------------
# Functions (alphabetical)
# -----------------------------------------------------------------------------

def get_file_list(search_string):
    """ Generates a list of files matching the given glob string

    Parameters:
        search_string: string
            glob-formatted search pattern

    Returns:
        file_list: list
            list of files matching pattern

    Outputs: nothing
    """

    file_list = glob.glob(search_string)
    return file_list

# -----------------------------------------------------------------------------

def get_cosmicx_params(header_data):
    """ Accesses the cosmicx configure file and returns the parameters.

    Parameters:
        instrument: string
            The intrument for the image. Can be 'WFPC2', 'WFC3', 'ACS'
        header_data: dictionary
                Dictionary containing information from FITS header

    Returns:
        cosmicx_params: dict
            A list of dictionaries, with one dictionary for each FITS
            extension. Each dictionary has the cosmicx parameters appropriate
            for that extension. If the extension is not a science extension,
            the list has None, instead of a dictionary, for that extension.

    """

    detector_config = { 'WFPC2' : "WFPC2_params.yaml",
                        'HRC' : "HRC_params.yaml",
                        'SBC' : "SBC_params.yaml",
                        'WFC' : "WFC_params.yaml",
                        'UVIS' : "UVIS_params.yaml",
                        'IR' : "IR_params.yaml" }

    detector = header_data["detector"]
    readnoise = header_data["readnoise"]
    gain = header_data["gain"]

    # Load the configuration file for the proper detector:
    config_file = detector_config[detector]

    # os.path.dirname is called three times to move up three directories
    # from this file, to reach to where the cosmicx_cfg directory is located.
    param_path = os.path.dirname(os.path.dirname(os.path.dirname(
                 os.path.abspath(inspect.getfile(inspect.currentframe())))))
    param_path = os.path.join(param_path,'cosmicx_cfg')
    param_path = os.path.join(param_path, config_file)

    cosmicx_params = yaml.load(open(param_path))

    # If we have readnoise and gain from the FITS header, override the
    # manually specified values from the cfg file:
    if readnoise:
        for extension_key in cosmicx_params:
            cosmicx_params[extension_key]["readnoise"] = readnoise
    if gain:
        for extension_key in cosmicx_params:
            cosmicx_params[extension_key]["gain"] = gain

    return cosmicx_params

# -----------------------------------------------------------------------------

def run_cosmicx(filename, output, cosmicx_params, detector):
    """ Driver to run lacosmicx on multi-extension FITS files.

    An equivalent to the run_cosmics function in run_cosmiscs.py,
    borrowing much of its code. Performs cosmic-ray rejection on each
    chip of a FITS file.

    Parameters:
        filename: string
            filename of input FITS file
        output: string
            desired filename of output, a cosmic ray rejected FITS file.
        cosmicx_params: dictionary
            a dictionary of dictionaries, with each entry dictionary
            containing the cosmicx settings appropriate for a
            particular FITS science extension.
        detector: string
            the detector used to take the FITS image. If SBC or IR,
            cr rejection is not run.

    Returns: nothing

    Outputs:
        A cosmic ray rejected FITS file with the filename given by
        'output'
        A symbolic link to the original c1m files that matches the
        output filename, for Astrodrizzle.

    """

    # Assert the input file exists
    error = filename + ' input for run_cosmicx in '
    error += 'run_cosmics.py does not exist.'
    assert os.access(filename, os.F_OK), error

    # If the file is taken by SBC or IR, do not cr reject, copy
    # the file with the output name, and return
    if detector == 'SBC' or detector == 'IR':
        shutil.copyfile(filename,output)
        return

    # Define the output name and delete if exists
    query = os.access(output, os.F_OK)
    if query == True:
        os.remove(output)

    with fits.open(filename, mode='readonly') as HDUlist:

        for key in cosmicx_params:

            params = cosmicx_params[key]

            # It's possible there may not be the same number of SCI extensions
            # in the image as we specify in the cfg files (WFC)
            try:
                HDU = HDUlist[key]
            except IndexError:
                continue

            # It's possible some of the SCI extensions might not have data,
            # as subsets of the CCD are sometimes used, and these extensions
            # are empty (UVIS)
            try:
                array = HDU.data
            except AttributeError:
                continue

            # For a cleaner cr rejection, set all very negative pixels
            # to 0. No significant science data should be lost, as these
            # are already bad pixels
            array[array < -10.0] = 0.0

            cleanarray = lacosmicx.run(array,
                       inmask=params['inmask'],
                       outmaskfile=params['outmaskfile'],
                       pssl=params['pssl'],
                       gain=params['gain'],
                       readnoise=params['readnoise'],
                       sigclip=params['sigclip'],
                       sigfrac=params['sigfrac'],
                       objlim=params['objlim'],
                       satlevel=params['satlevel'],
                       robust=params['robust'],
                       verbose=params['verbose'],
                       niter=params['niter'] )

            HDU.data = cleanarray

        HDUlist.writeto(output)

    # Create the symbolic link
    make_c1m_link(os.path.abspath(filename),
                  os.path.abspath(output))

# -----------------------------------------------------------------------------

def make_c1m_link(filename, output):
    """ Create a link to a c1m.fits that matches the cosmic ray rejected
    naming scheme. This is only done for WFPC2 data (which ends in _c0m.fits)

    Parameters:
        filename: string
            filename of the file to be cleaned.

    Returns: nothing

    Outputs:
        A symbolic link to 'filename'
    """

    if filename[-8:] == 'c0m.fits':
        src = filename.replace('_c0m.fits', '_c1m.fits')
        dst = output.replace('_c0m.fits', '_c1m.fits')
        query = os.path.islink(dst)

        if query == True:
            os.remove(dst)
        os.symlink(src, dst)

max_stars_repo_name: STScI-Citizen-Science/MTPipeline

max_issues_repo_name: STScI-Citizen-Science/MTPipeline

max_forks_repo_name: STScI-Citizen-Science/MTPipeline

reason: from astropy