text:
__author__ = 'mike'

import sys

from blocks.bricks.base import application
from blocks.bricks.recurrent import BaseRecurrent, recurrent
from blocks.bricks import Linear, Sigmoid, lazy, Initializable
from blocks.initialization import IsotropicGaussian, Constant
import theano
import theano.tensor as T
import numpy as np
from blocks.graph import ComputationGraph

from utils import test_value
from rbm import Rbm
from rnn import Rnn

sys.setrecursionlimit(10000)

from theano.sandbox.rng_mrg import MRG_RandomStreams as RandomStreams

rng = RandomStreams(seed=np.random.randint(1 << 30))

floatX = theano.config.floatX

class Rnnrbm(BaseRecurrent, Initializable):
    @lazy(allocation=['visible_dim', 'hidden_dim'])
    def __init__(self, visible_dim, hidden_dim, rnn_dimensions=(128, 128), **kwargs):
        super(Rnnrbm, self).__init__(**kwargs)
        self.rnn_dimensions = rnn_dimensions
        self.visible_dim = visible_dim
        self.hidden_dim = hidden_dim

        # self.in_layer = Linear(input_dim=input_dim, output_dim=rnn_dimension * 4,
        # weights_init=IsotropicGaussian(0.01),
        # biases_init=Constant(0.0),
        # use_bias=False,
        # name="in_layer")

        self.rbm = Rbm(visible_dim=visible_dim, hidden_dim=hidden_dim,
                       activation=Sigmoid(), weights_init=IsotropicGaussian(0.1),
                       biases_init=Constant(0.1),
                       name='rbm')

        self.uv = Linear(input_dim=rnn_dimensions[-1], output_dim=visible_dim,
                         weights_init=IsotropicGaussian(0.0001),
                         biases_init=Constant(0.001),
                         use_bias=True, name='uv')

        self.uh = Linear(input_dim=rnn_dimensions[-1], output_dim=hidden_dim,
                         weights_init=IsotropicGaussian(0.0001),
                         biases_init=Constant(0.001),
                         use_bias=True, name='uh')

        self.rnn = Rnn([visible_dim] + list(rnn_dimensions), name='rnn')

        self.children = [self.rbm, self.uv, self.uh, self.rnn] + self.rnn.children._items

    def initialize(self, pretrained_rbm=None, pretrained_rnn=None, **kwargs):
        super(Rnnrbm, self).initialize(**kwargs)
        # if pretrained_rbm is not None:
        # self.rbm.bv.set_value(pretrained_rbm.bv.get_value())
        # self.uv.b.set_value(pretrained_rbm.bv.get_value())
        #     self.rbm.bh.set_value(pretrained_rbm.bh.get_value())
        #     self.uh.b.set_value(pretrained_rbm.bh.get_value())
        #     self.rbm.W.set_value(pretrained_rbm.W.get_value())
        # if pretrained_rnn is not None:
        #     for param, trained_param in zip(itertools.chain(*[child.params for child in self.rnn.children]),
        #                                     itertools.chain(*[child.params for child in pretrained_rnn.children])):
        #         param.set_value(trained_param.get_value())
        self.uv.b.name = 'buv'
        self.uv.W.name = 'Wuv'
        self.uh.W.name = 'Wuh'
        self.uh.b.name = 'buh'
        self.rbm.bh.name = 'buh'
        self.rbm.bv.name = 'buv'
        # self.rnn.input_transform.b.name = 'bvu'
        # self.rnn.input_transform.W.name = 'Wvu'

    # def _allocate(self):
    # Wrbm = shared_floatx_nans((self.visible_dim, self.hidden_dim), name='Wrbm')
    # add_role(Wrbm, WEIGHT)
    # self.params.append(Wrbm)
    # self.add_auxiliary_variable(Wrbm.norm(2), name='Wrbm_norm')
    #
    #
    # def _initialize(self):
    # for param in self.params:
    # if has_roles(param, WEIGHT):
    # self.weights_init.initialize(param, self.rng)
    # elif has_roles(param, BIAS):
    # self.biases_init.initialize(param, self.rng)

    def get_dim(self, name):
        dims = {'visible': self.visible_dim,
                'mask': self.visible_dim,
                'gru_state': self.rnn_dimensions[0],
                'lstm_state': self.rnn_dimensions[1],
                'lstm_cells': self.rnn_dimensions[1]
                }
        return dims.get(name, None) or super(Rnnrbm, self).get_dim(name)

    @recurrent(sequences=[], contexts=['visible'], states=['gru_state', 'lstm_state', 'lstm_cells'],
               outputs=['visible', 'gru_state', 'lstm_state', 'lstm_cells'])
    def generate(self, visible=None, gru_state=None, lstm_state=None, lstm_cells=None, rbm_steps=25):
        bv = self.uv.apply(lstm_state)
        bh = self.uh.apply(lstm_state)
        visible_start = rng.binomial(size=visible.shape, n=1, p=bv, dtype=floatX)
        _, visible = self.rbm.apply(visible=visible_start, bv=bv, bh=bh, n_steps=rbm_steps,
                                    batch_size=visible.shape[0])
        visible = visible[-1]
        gru_state, lstm_state, lstm_cells = self.rnn.rnn_apply(inputs=visible, gru_state=gru_state,
                                                               lstm_state=lstm_state,
                                                               lstm_cells=lstm_cells, )  # iterate=False)
        # input_transform = self.vu.apply(visible)
        # gru_state = self.gru_layer.apply(inputs=input_transform, states=gru_state,
        # iterate=False)
        # lstm_state, lstm_cells = self.lstm_layer.apply(inputs=gru_state, states=lstm_state, cells=lstm_state,
        # iterate=False)

        updates = ComputationGraph(lstm_state).updates
        # output = 1.17 * self.out_layer.apply(hidden_state)
        return [visible, gru_state, lstm_state, lstm_cells], updates

    @recurrent(sequences=['visible', 'mask'], contexts=[],
               states=['gru_state', 'lstm_state', 'lstm_cells'],
               outputs=['gru_state', 'lstm_state', 'lstm_cells', 'bv', 'bh'])
    def training_biases(self, visible, mask=None, gru_state=None, lstm_state=None, lstm_cells=None):
        bv = self.uv.apply(lstm_state)
        bh = self.uh.apply(lstm_state)
        # inputs = rbm.apply(visible=visible, bv=bv, bh=bh, n_steps=25)
        # input_transform = self.vu.apply(visible)
        # gru_state = self.gru_layer.apply(inputs=input_transform, states=gru_state,
        # iterate=False)
        # lstm_state, lstm_cells = self.lstm_layer.apply(inputs=gru_state, states=lstm_state,
        #                                                cells=lstm_state,
        #                                                iterate=False)
        gru_state, lstm_state, lstm_cells = self.rnn.rnn_apply(visible, mask=mask, gru_state=gru_state,
                                                               lstm_state=lstm_state,
                                                               lstm_cells=lstm_cells, )  #iterate=False)
        updates = ComputationGraph(lstm_state).updates
        return [gru_state, lstm_state, lstm_cells, bv, bh], updates

    def cost(self, examples, mask, k=10):
        _, _, _, bv, bh = self.training_biases(visible=examples, mask=mask)
        cost, v_samples = self.rbm.cost(visible=examples, bv=bv, bh=bh, k=k,
                                        batch_size=examples.shape[0])
        return cost.astype(floatX), v_samples.astype(floatX)

    # @application(inputs=['visible', 'k', 'batch_size', 'mask'], outputs=['vsample', 'cost'])
    def rbm_pretrain_cost(self, visible, k, batch_size, mask=None):
        return self.rbm.cost(visible=visible, k=k, batch_size=batch_size, mask=mask)

    # @application(inputs=['x', 'x_mask'], outputs=['output_'])
    def rnn_pretrain_pred(self, x, x_mask):
        return self.rnn.apply(inputs=x, mask=x_mask)

if __name__ == "__main__":
    # v = np.float32(np.random.randn(10, visible_dim))

    x = T.tensor3('features')
    x_mask = T.matrix('features_mask')
    y = T.tensor3('targets')
    y_mask = T.matrix('targets_mask')

    x = test_value(x, np.ones((15, 10, 88), dtype=floatX))
    y = test_value(y, np.ones((15, 10, 88), dtype=floatX))
    x_mask = test_value(x_mask, np.ones((15, 10), dtype=floatX))
    y_mask = test_value(y_mask, np.ones((15, 10), dtype=floatX))

    rnnrbm = Rnnrbm(88, 256, rnn_dimensions=(100, 100))
    rnnrbm.allocate()
    rnnrbm.initialize()

    # # y = rnnrbm.generate(n_steps=5, batch_size=x.shape[1])
    # # rbm = Rbm(visible_dim, hidden_dim, weights_init=IsotropicGaussian(0.1), use_bias=False, name='rbm')
    # # rbm.allocate()
    # # rbm.initialize()
    # # # print rbm.initial_state
    #
    # # a = T.ones(visible_dim, dtype=floatX) * 0.5
    # # b = T.ones(hidden_dim, dtype=floatX) * 0.5
    # # y = rbm.apply(visible=x, bv=a, bh=b,
    # # n_steps=10, batch_size=10)
    # cg = ComputationGraph(y)
    # updates = cg.updates
    # #
    # f = theano.function([x, x_mask], y, updates=updates)
    # for i in range(100):
    # t = time.time()
    # print f(np.ones((15, 10, 93), dtype=floatX), np.ones((15, 10), dtype=floatX))
    #     print "took: " + str(time.time() - t)

reason: import numpy,import theano,from theano


text:
"""
This pipeline is intended to extract LBP features from ADC images.
"""

import os

import numpy as np

from protoclass.data_management import ADCModality
from protoclass.data_management import GTModality

from protoclass.preprocessing import PiecewiseLinearNormalization

from protoclass.extraction import LBPExtraction

# Define the path where all the patients are
path_patients = '/data/prostate/experiments'
# Define the path of the modality to normalize
path_adc = 'ADC_reg_bspline'
# Define the path of the ground for the prostate
path_gt = 'GT_inv/prostate'
# Define the label of the ground-truth which will be provided
label_gt = ['prostate']
# Define the path where the information for the piecewise-linear normalization
path_piecewise = '/data/prostate/pre-processing/mp-mri-prostate/piecewise-linear-adc'
# Define the path to store the Tofts data
path_store = '/data/prostate/extraction/mp-mri-prostate/lbp-adc'

# Generate the different path to be later treated
path_patients_list_adc = []
path_patients_list_gt = []
# Create the generator
id_patient_list = [name for name in os.listdir(path_patients)
                   if os.path.isdir(os.path.join(path_patients, name))]
for id_patient in id_patient_list:
    # Append for the ADC data
    path_patients_list_adc.append(os.path.join(path_patients, id_patient,
                                               path_adc))
    # Append for the GT data - Note that we need a list of gt path
    path_patients_list_gt.append([os.path.join(path_patients, id_patient,
                                               path_gt)])

# List where to store the different minimum
for id_p, (p_adc, p_gt) in enumerate(zip(path_patients_list_adc,
                                         path_patients_list_gt)):

    print 'Processing {}'.format(id_patient_list[id_p])

    # Remove a part of the string to have only the id
    nb_patient = id_patient_list[id_p].replace('Patient ', '')

    # Read the image data
    adc_mod = ADCModality()
    adc_mod.read_data_from_path(p_adc)

    # Read the GT
    gt_mod = GTModality()
    gt_mod.read_data_from_path(label_gt, p_gt)

    # Read the normalization information
    pat_chg = id_patient_list[id_p].lower().replace(' ', '_') + '_norm.p'
    filename = os.path.join(path_piecewise, pat_chg)
    adc_norm = PiecewiseLinearNormalization.load_from_pickles(filename)

    # Normalize the data
    adc_mod = adc_norm.normalize(adc_mod)

    # Create the Haralick extractor
    p_arr = (8, 16)
    r_arr = (1., 2.)
    for p, r in zip(p_arr, r_arr):
        ext = LBPExtraction(adc_mod, p=p, r=r)

        # Fit the data
        print 'Compute the LBP maps'
        ext.fit(adc_mod, ground_truth=gt_mod, cat=label_gt[0])

        # Extract the data
        print 'Extract the LBP using the GT'
        data = ext.transform(adc_mod, ground_truth=gt_mod, cat=label_gt[0])

        # Store the data
        print 'Store the data in the right directory'

        # Check that the path is existing
        if not os.path.exists(path_store):
            os.makedirs(path_store)
        pat_chg = (id_patient_list[id_p].lower().replace(' ', '_') +
                   '_lbp_' + str(p) + '_' + str(int(r)) + '_adc.npy')
        filename = os.path.join(path_store, pat_chg)
        np.save(filename, data)

ext: py

max_stars_repo_name: DoraSzasz/mp-mri-prostate

max_issues_repo_name: DoraSzasz/mp-mri-prostate

reason: import numpy


text:
import os, cv2
import matplotlib.pyplot as plt
import numpy as np
import requests

class Image :

	def __init__(self, image_path, local_feature_extractor, image_filter=None, global_feature_extractor=None, url=False):
		self.name = os.path.splitext(os.path.basename(image_path))[0]
		self.path = image_path
		if url:
			self.image = self.read_image_gray(self.fetch_from_url(image_path))
		else:
			self.image = self.read_image_gray(cv2.imread(image_path))
		self.key_points, self.local_descriptors = local_feature_extractor.kps_descriptors(self.image)
		if image_filter is not None:
			image_filter(self.image)
		if global_feature_extractor is not None:
			self.set_global_descriptor(global_feature_extractor)

	def fetch_from_url(self, image_path):
		resp = requests.get(image_path, stream=True).raw
		image = np.asarray(bytearray(resp.read()), dtype="uint8")
		image = cv2.imdecode(image, cv2.IMREAD_COLOR)
		return image

	def set_global_descriptor(self, global_feature_extractor):
		self.global_descriptor = global_feature_extractor.global_feature_vector(self.local_descriptors)

	def has_global_descriptor(self):
		if not hasattr(self, 'global_descriptor'):
			return False
		return True

	def read_image_gray(self, image_path, margin=30):
		return cv2.cvtColor(image_path,cv2.COLOR_BGR2GRAY)[:,margin:-margin]

max_stars_repo_name: Nagisa-sys/WriterIdentification

max_issues_repo_name: Nagisa-sys/WriterIdentification

max_forks_repo_name: Nagisa-sys/WriterIdentification


text:
#!/usr/bin/env python3

## sys build-in package
import tensorflow
import os
import numpy as np
import logging
# set up logging
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

## keras-related package
import keras
from keras import backend as K
from keras.models import Model
from keras.layers import Activation, Input, Conv2D, Conv3D, MaxPooling2D, UpSampling2D, concatenate, BatchNormalization
from keras.layers import add, multiply, LeakyReLU, PReLU, ThresholdedReLU
from keras.initializers import he_normal
from keras.initializers import Constant
from keras.engine.topology import Layer

def m0(model = None, img_row=512, img_col=512, img_slice = 5, img_channel=1, layers = 32, features = 64,
		act_type ='relu', scale_factor = 0.1, init_seed=15213, batch_norm=False,is_pred=False, prow=512,
		pcol=512, pslice=1, pchannel=1):

	logging.info("Current version is 3D : LR -> HR")
	logging.info("Scaling factor is " + str(scale_factor))
	logging.info("Activation function " + act_type)

	# Initialize the inputs
	if is_pred == False:
		inputs = Input((img_row, img_col, img_slice, img_channel))
	else:
		inputs = Input((prow, pcol, pslice, pchannel))
	# Perform the first convolution and save it for residuals at the end
	x = Conv3D(features, kernel_size = (3,3,3),
				kernel_initializer = he_normal(seed=init_seed),
				padding = 'same',
				name = ('conv0'))(inputs)
	conv_1 = x

	# Perform convolutions for the residual residual_block
	for blk_num in range(layers):
		x = residual_block_3d(x, features, act_type, scale_factor, blk_num, init_seed)

	# Perform the penultimate convolution and add to the first conv
	x = Conv3D(features, kernel_size = (3,3,3),
				kernel_initializer = he_normal(seed=init_seed),
				padding = 'same',
				name = ('conv_penultimate'))(x)
	x = add([x, conv_1])

	# Final conv to create the super-res output
	output = Conv3D(1, kernel_size = (3,3,3),
				kernel_initializer = he_normal(seed=init_seed),
				padding = 'same',
				name = ('conv_final'))(x)

	model = Model(inputs = [inputs], outputs = [output])

	return model

class ConstMultiplierLayer(Layer):
    def __init__(self, val=0.1, **kwargs):
        self.val = val
        super(ConstMultiplierLayer, self).__init__(**kwargs)

    # Create a function to initialize this layer with a specific value
    def build(self, input_shape):
        self.k = self.add_weight(
            name='k',
            shape=(),
            initializer=Constant(value=self.val),
            dtype='float32',
            trainable=True,
        )
        super(ConstMultiplierLayer, self).build(input_shape)

    def call(self, x):
        return K.tf.multiply(self.k, x)

    def compute_output_shape(self, input_shape):
        return input_shape

def residual_block_3d(x, filter_size, act_type, scale_factor, blk_num, init_seed):
	tmp = Conv3D(filter_size, kernel_size = (3,3,3), kernel_initializer = he_normal(seed=init_seed),
				padding = 'same',
				name = ('resblock_%i_conv_1' % blk_num))(x)
	if act_type =="LeakyReLU":
		tmp = LeakyReLU(alpha=0.2)(tmp)
	else:
		tmp = Activation(act_type, name = ('resblock_%i_act' % blk_num))(tmp)

	tmp = Conv3D(filter_size, kernel_size = (3,3,3), kernel_initializer = he_normal(seed=init_seed),
				padding = 'same',
				name = ('resblock_%i_conv_2' % blk_num))(tmp)

	tmp = ConstMultiplierLayer(val = scale_factor)(tmp)
	out = keras.layers.Add()([x, tmp])

	return out

max_stars_repo_name: subtlemedicalinc/fastMRI

max_issues_repo_name: subtlemedicalinc/fastMRI

max_forks_repo_name: subtlemedicalinc/fastMRI


text:
# backend/main.py

from concurrent.futures import process
import uuid

import cv2
import uvicorn
from fastapi import File
from fastapi import FastAPI
from fastapi import UploadFile
import numpy as np
from PIL import Image

import config
import inference

import asyncio
from concurrent.futures import ProcessPoolExecutor
from functools import partial

import time

app = FastAPI()

@app.get("/")
def read_root():
    return {"message": "Welcome from the API"}

@app.post("/{style}")
async def get_image(style: str, file: UploadFile = File(...)):
    image = np.array(Image.open(file.file))
    model = config.STYLES[style]
    start = time.time()
    output, resized = inference.inference(model, image)
    name = f"/storage/{str(uuid.uuid4())}.jpg"
    cv2.imwrite(name, output)
    models = config.STYLES.copy()
    del models[style]
    asyncio.create_task(generate_remaining_models(models, image, name))
    return {"name": name, "time": time.time() - start}

async def generate_remaining_models(models, image, name: str):
    executor = ProcessPoolExecutor()
    event_loop = asyncio.get_event_loop()
    await event_loop.run_in_executor(
        executor, partial(process_image, models, image, name)
    )

def process_image(models, image, name: str):
    for model in models:
        output, resized = inference.inference(models[model], image)
        name = name.split(".")[0]
        name = f"{name.split('_')[0]}_{models[model]}.jpg"
        cv2.imwrite(name, output)

if __name__ == "__main__":
    uvicorn.run("main:app", host="0.0.0.0", port=8080)


text:
from __future__ import print_function
import numpy as np
import osgeo.osr as osr

def is_registered(pointcloud):
    """
    Returns True when a pointcloud is registered; ie coordinates are relative
    to a specific spatial reference system or offset.

    In that case, first transform one pointcloud to the reference system
    of the other, before doing processing on the points:

       set_srs(pcA, same_as=pcB)
    """
    return hasattr(pointcloud, 'srs') or hasattr(pointcloud, 'offset')

def same_srs(pc_one, pc_two):
    """
    True if the two pointclouds have the same coordinate system

    Arguments:
        pc_one : pcl.PointCloud
        pc_two : pc..PointCloud
    """

    is_reg_one = is_registered(pc_one)
    is_reg_two = is_registered(pc_two)

    # both pointcloud are pure pcl: no offset nor SRS
    if not is_reg_one and not is_reg_two:
        return True

    # only one of the pointclouds is registered
    if ((not is_reg_one) and is_reg_two) or ((not is_reg_two) and is_reg_one):
        return False

    srs_one = None
    srs_two = None
    try:
        srs_one = pc_one.srs
        srs_two = pc_two.srs

        if not srs_one.IsSame(srs_two):
            # SRS present, but different
            return False
    except:
        # one of the pointclouds does not have a SRS
        return False

    off_one = None
    off_two = None
    try:
        off_one = pc_one.offset
        off_two = pc_two.offset
    except:
        # one of the pointclouds does not have an offset
        return False

    # absolute(off_one - off_two) <= (atol + rtol * absolute(off_two))
    if not np.allclose(off_one, off_two, rtol=1e-06, atol=1e-08):
        return False

    return True

def set_srs(pc, srs=None, offset=np.array([0, 0, 0], dtype=np.float64),
            same_as=None):
    """Set the spatial reference system (SRS) and offset for a pointcloud.
    This function transforms all the points to the new reference system, and
    updates the metadata accordingly.

    Either give a SRS and offset, or a reference pointcloud

    NOTE: Pointclouds in PCL do not have absolute coordinates, ie.
          latitude / longitude. This function sets metadata to the pointcloud
          describing an absolute frame of reference.
          It is left to the user to make sure pointclouds are in the same
          reference system, before passing them on to PCL functions. This
          can be checked with patty.utils.same_srs().

    NOTE: To add a SRS to a point cloud, or to update incorrect metadata,
          use force_srs().

    Example:

        # set the SRS to lat/lon,
        # don't use an offset, so it defaults to [0,0,0]
        set_srs( pc, srs="EPSG:4326" )

    Arguments:
        pc : pcl.Pointcloud, with pcl.is_registered() == True

        same_as : pcl.PointCloud

        offset : np.array([3], dtype=np.float64 )
            Must be added to the points to get absolute coordinates,
            neccesary to retain precision for LAS pointclouds.

        srs : object or osgeo.osr.SpatialReference
            If it is an SpatialReference, it will be used directly.
            Otherwise it is passed to osr.SpatialReference.SetFromUserInput()

    Returns:
        pc : pcl.PointCloud
            The input pointcloud.
    """
    if not is_registered(pc):
        raise TypeError("Pointcloud is not registered")
        return None

    update_offset = False
    update_srs = False

    if same_as:

        # take offset and SRS from reference pointcloud
        update_offset = True
        update_srs = True
        if is_registered(same_as):
            newsrs = same_as.srs
            newoffset = same_as.offset
        else:
            raise TypeError("Reference pointcloud is not registered")
    else:

        # take offset and SRS from arguments

        # sanitize offset
        if offset is not None:
            update_offset = True

            newoffset = np.array(offset, dtype=np.float64)
            if len(newoffset) != 3:
                raise TypeError("Offset should be an np.array([3])")

        # sanitize SRS
        if srs is not None:

            # argument is a SRS, use it
            if type(srs) == type(osr.SpatialReference()):
                update_srs = True
                newsrs = srs

            # argument is not an SRS, try to convert it to one
            elif isinstance(srs, osr.SpatialReference):
                newsrs = osr.SpatialReference()
                if newsrs.SetFromUserInput(srs) == 0:
                    update_srs = True

            # illegal input
            else:
                raise TypeError(
                    "SRS should be a string or a osr.SpatialReference")
    # Apply

    # add old offset
    data = np.asarray(pc)
    precise_points = np.array(data, dtype=np.float64) + pc.offset

    # do transformation, this resets the offset to 0
    if update_srs and not pc.srs.IsSame(newsrs):
        try:
            transform = osr.CoordinateTransformation(pc.srs, newsrs)
            precise_points = np.array(transform.TransformPoints(precise_points),
                                      dtype=np.float64)
            pc.srs = newsrs.Clone()
            pc.offset = np.array([0, 0, 0], dtype=np.float64)
        except:
            print("WARNING, CAN'T DO COORDINATE TRANSFORMATION")

    # substract new offset
    if update_offset:
        precise_points -= newoffset
        pc.offset = np.array(newoffset, dtype=np.float64)

    # copy the float64 to the pointcloud
    data[...] = np.asarray(precise_points, dtype=np.float32)

    return pc

def force_srs(pc, srs=None, offset=None, same_as=None):
    """
    Set a spatial reference system (SRS) and offset for a pointcloud.
    Either give a SRS and offset, or a reference pointcloud
    This function affects the metadata only.

    This is the recommended way to turn a python-pcl pointcloud to a
    registerd pointcloud with absolute coordiantes.

    NOTE: To change the SRS for an already registered pointcloud, use set_srs()

    Example:

        # set the SRS to lat/lon, leave offset unchanged
        force_srs( pc, srs="EPSG:4326" )

    Arguments:
        pc : pcl.Pointcloud

        same_as : pcl.PointCloud

        offset : np.array([3])
            Must be added to the points to get absolute coordinates,
            neccesary to retain precision for LAS pointclouds.

        srs : object or osgeo.osr.SpatialReference
            If it is an SpatialReference, it will be used directly.
            Otherwise it is passed to osr.SpatialReference.SetFromUserInput()

    Returns:
        pc : pcl.PointCloud
            The input pointcloud.
    """
    if same_as:
        if is_registered(same_as):
            pc.srs = same_as.srs.Clone()
            pc.offset = np.array(same_as.offset, dtype=np.float64)
    else:
        if type(srs) == type(osr.SpatialReference()):
            pc.srs = srs.Clone()
        elif srs is not None:
            pc.srs = osr.SpatialReference()
            pc.srs.SetFromUserInput(srs)
        else:
            pc.srs = osr.SpatialReference()

        if offset is not None:
            offset = np.asarray(offset, dtype=np.float64)
            if len(offset) != 3:
                raise TypeError("Offset should be an np.array([3])")
            else:
                pc.offset = offset

    return pc

max_stars_repo_name: NLeSC/PattyAnalytics

max_issues_repo_name: NLeSC/PattyAnalytics

max_forks_repo_name: NLeSC/PattyAnalytics


text:
import cv2
import pyrealsense2 as rs2
import numpy as np


text:
# -*- coding: utf-8 -*-
"""
This module handles the identification and cleaning of  noise in signals. A 'noisy' signal
is one that has spikes in it or sudden variations in a continuous looking
function. Most of these functions are built to detect and handle sudden
spikes in the frequency profile estimates of a sound.

"""
import numpy as np
from scipy import ndimage, stats
from itsfm.signal_processing import moving_rms_edge_robust
from itsfm.signal_processing import median_filter, resize_by_adding_one_sample
from itsfm.signal_processing import dB
from itsfm.sanity_checks import make_sure_its_positive

def exterpolate_over_anomalies(X, fs, anomalous, **kwargs):
    '''
    Ex(tra)+(in)ter-polates --> Exterpolates over  anomalous regions. Anomalous
    regions are either 'edge' or 'island' types. The 'edge' anomalies are those which are
    at the extreme ends of the signal. The 'island' anomalies are regions with
    non-anomalous regions on the left and right.

    An 'edge' anomalous region is handled by running a linear regression on the
    neighbouring non-anomalous region, and using the slope to extrapolate over
    the edge anomaly.

    An 'island' anomaly is handled by interpolating between the end values of the
    neighbouring non-anomalous regions.

    Parameters
    ----------
    X : np.array
    fs : float>0
        Sampling rate in Hz
    anomalous : np.array
        Boolean array of same size as X
        True indicates an anomalous sample.
    extrap_window : float>0, optional
        The duration of the extrapolation window in seconds.
        Defaults to 0.1ms

    Returns
    -------
    smooth_X : np.array
        Same size as X, with the anomalous regions

    Notes
    -----
    Only extrapolation by linear regression is supported currently. The `extrap_window`
    parameter is important especially if there is a high rate of frequency modulation
    towards the edges of the sound. When there is a high freq. mod. at the edges it
    is better to set the `extrap_window` small. However, setting it too small also
    means that the extrapolation may not be as nice anymore.

    Example
    -------

    `not up to date!!!`

    See Also
    --------
    find_closest_normal_region

    '''
    smooth_X = X.copy()
    extrap_window = kwargs.get('extrap_window', 0.1*10**-3)
    ref_region_length = int(extrap_window*fs)
    try:
        anomalous_broader = ndimage.filters.percentile_filter(anomalous, 100,
                                                              ref_region_length)
    except:
        raise ValueError(f"Unable to percentile filter with kernel size {ref_region_length} and sampling rate {fs}")

    anomalous_labelled, num_regions = ndimage.label(anomalous_broader)
    if num_regions == 0:
        return smooth_X

    anomalous_regions = ndimage.find_objects(anomalous_labelled)

    for each_region in anomalous_regions:
        region_type = anomaly_type(each_region, X)
        if region_type == 'edge':
            smooth_X[each_region] = anomaly_extrapolation(each_region, X,
                                                            ref_region_length)
        elif region_type == 'island':
            smooth_X[each_region] = fix_island_anomaly(X, fs, each_region,
                                                       ref_region_length,
                                                       **kwargs)
    return smooth_X

def fix_island_anomaly(X, fs, anomaly, ref_region_length, **kwargs):
    '''
    First tries to interpolate between the edges of the anomaly at hand.
    If the interpolation leads to a very drastic slope, a 'sensible' extrapolation
    is attempted using parts of the non-anomalous signal.

    Parameters
    ----------
    X : np.array
    fs : float>0
    anomaly : tuple slice
        scipy.ndimage.find_objects output
        (slice(start,stop,None),)
    ref_region_length : int>0
        The number of samples to be used as a reference region in
        case of extrapolation
    max_fmrate : float>0, optional
        The maximum fm rate to be tolerated while interpolating in kHz/ms
        Defaults to 100 kHz/ms.

    Returns
    -------
    interpolated : np.array
        Array of same size as anomaly.
    '''
    max_fmrate = kwargs.get('max_fmrate', 100)
    trial_fix = anomaly_interpolation(anomaly, X)
    fmrate_trialfix = calc_coarse_fmrate(trial_fix, fs)
    if fmrate_trialfix <= max_fmrate:
        return trial_fix
    else:
        return extrapolate_sensibly(X, fs, anomaly, ref_region_length, **kwargs)

def extrapolate_sensibly(X, fs, anomaly, ref_region_length, **kwargs):
    '''
    Function called when `fix_island_anomaly` detects direct interpolation
    will lead to unrealistic slopes. This function is called when there's
    a big difference in values across an anomalous region and an
    extrapolation must be performed which will not alter the signal drastically.

    The method tries out the following:
        #. Look left and right of the anomaly to see which region
           has higher frequency content.
        #. Extrapolate in the high-to-low frequency direction.

    This basically means that if the local inspection window around anomaly has
    a sweep between 20-10kHZ on the left and a 0Hz region on the right - the
    anomaly will be extrapolated with the slope from the sweep region because it
    has higher frequency content.

    Example
    -------
    >>> freq_profile = [np.zeros(10), np.arange(15,30,5)*1000]
    >>> fs = 1.0
    >>> x = np.concatenate(freq_profile)[::-1]
    >>> anom = (slice(2, 5, None),)
    >>>
    >>> plt.plot(x, label='noisy frequency profile')
    >>> anom_x = np.zeros(x.size, dtype='bool')
    >>> anom_x[anom[0]] = True
    >>> plt.plot(anom_x*8000, label='identified anomaly')
    >>> extrap_out = extrapolate_sensibly(x, fs, anom, 4)
    >>> sensibly_extrap = x.copy()
    >>> sensibly_extrap[anom_x] = extrap_out
    >>> plt.plot(sensibly_extrap, label='extrapolated')
    >>> plt.legend()
    '''

    left_and_right_of_X = get_neighbouring_regions(X, anomaly, ref_region_length)
    left_median, right_median = map(np.median, left_and_right_of_X)

    start, stop = anomaly[0].start, anomaly[0].stop
    anom_size = stop-start
    if left_median > right_median:
        relevant_region = anomaly
        extrap_chunk = anomaly_extrapolation(relevant_region, X[:stop],
                                                             ref_region_length,
                                                             **kwargs)

    else:
        relevant_region = (slice(0,anom_size,None),)
        extrap_chunk = anomaly_extrapolation(relevant_region,
                                             X[start:],
                                             ref_region_length, **kwargs)

    return extrap_chunk

def get_neighbouring_regions(X, target, region_size):
    '''
    Takes out samples of `region_size` on either size of the target.

    Parameters
    ----------
    X: np.array
    target : slice
        ndimage.find_objects type slice
    region_size : int >0

    Returns
    -------
    left_and_right : list
    '''
    start = target[0].start
    stop = target[0].stop

    before_start = start-region_size
    if before_start < 0 :
        left_of_X = X[:start]
    else:
        left_of_X = X[before_start:start]

    try:
        right_of_X = X[stop:stop+region_size]
    except:
        right_of_X = X[stop:]
    return [left_of_X, right_of_X]

def calc_coarse_fmrate(X,fs,**kwargs):
    '''
    Calculates slope by subtracting the difference between 1st and
    last sample and dividing it by the length of the array.
    The output is then converted to units of kHz/ms.

    Parameters
    ----------
    X : np.array
        Frequency profile with values in Hz.
    fs : float>0

    '''
    diff = X[-1] - X[0]
    length = X.size/fs
    return  np.abs((diff/length)*10**-6)

def anomaly_extrapolation(region, X, num_samples, **kwargs):
    '''
    Takes X values next to the region and fits a linear regression
    into the region. This is only suitable for cases where the
    anomalous region is at an 'edge' - either one of its samples
    is 0 or the last sample of X.

    Parameters
    ----------
    region : object tuple
        A slice type object which is the output from scipy.ndimage.find_objects
        This is a slice inside a list/tuple.
    X : np.array
        The original array over which the extrapolation is to be performed
    num_samples : int>0
        The number of samples next to the region to be used to fit the data
        for extrapolation into the region.

    Returns
    -------
    extrapolated : np.array
        The values corresponding to the extrapolated region.

    Notes
    ------
    1. This function covers 90% of cases...if there is an anomaly right next
    to an edge anomaly with <num_samples distance -- of course things will
    go whack.

    Warning
    -------
    A mod on this function also allows extrapolation to occur if there
    are < num_samples next to the anomaly - this might make the function
    a bit lax in terms of the extrapolations it produces.

    '''
    start, stop = region[0].start, region[0].stop
    x = np.arange(start,stop)
    if start == 0:
        try:
            ref_x = range(stop, stop+num_samples)
            ref_values = X[ref_x]
        except:
            ref_x = range(stop, X.size)
            ref_values = X[stop:]

    elif stop==X.size:
        try:
            if start-num_samples <0:
                raise ValueError()
            else:
                ref_x = range(start-num_samples, start)
                ref_values = X[ref_x]
        except:
            ref_x = range(start)
            ref_values = X[ref_x]
    else:
        print(start, stop, x, 'anomaly x', X.size)
        raise NotImplementedError('the handling of none-edge case is not yet doen')
    m, c,rv, pv, stderr = stats.linregress(ref_x, ref_values)
    extrapolated = m*x + c
    return extrapolated

def anomaly_interpolation(region, X, **kwargs):
    '''
    Interpolates X values using values of X adjacent to the
    region.

    Parameters
    ----------
    region : object tuple
        Output from scipy.ndimage.find_objects
    X : np.array

    Returns
    -------
    full_span : np.array
        The values of interpolated X, of same size as the
        region length.
    '''
    start, stop = region[0].start, region[0].stop
    left_point = start-1
    full_span = np.linspace(X[left_point],X[stop],stop-left_point+1)
    return full_span[1:-1]

def anomaly_type(region, X):
    start, stop = region[0].start, region[0].stop

    at_left_edge = start==0
    at_right_edge = stop==X.size

    if np.logical_and(at_left_edge, at_right_edge):
        raise ValueError('The anomaly spans the whole array - please check again')
    at_either_edge = np.logical_or(at_left_edge, at_right_edge)

    if at_either_edge:
        return 'edge'
    else:
        return 'island'

def smooth_over_potholes(X, fs, **kwargs):
    '''
    A signal can show drastic changes in its value because of measurement errors.
    These drastic variations in signal are called `potholes <https://en.wikipedia.org/wiki/Pothole>`_
    (uneven parts of a road). This method tries to 'level' out the pothole by re-setting the samples of the
    pothole. A linear interpolation is done from the start of a pothole till its end using the closest
    non-pothole samples.

    A pothole is identified by a region of the signal with drastic changes in slope. A moving window
    calculates N slopes between the focal sample and the Nth sample after it to estimate if
    the Nth sample could be part of a pothole or not.

    Parameters
    ----------
    X : np.array
    fs : float>0
    max_stepsize : float>0, optional
        The maximum absolute difference between adjacent samples.
        Defaults to 50.
    pothole_inspection_window : float>0, optional
        The length of the moving window that's used to discover potholes.
        See identify_pothole_samples for default value.

    Returns
    -------
    pothole_covered
    pothole_regions

    See Also
    --------
    identify_pothole_samples
    pothole_inspection_window

    '''
    kwargs['max_stepsize'] = kwargs.get('max_stepsize', 50)
    potholes = identify_pothole_samples(X, fs, **kwargs)

    abnormal_fm, num_regions = ndimage.label(potholes)
    if num_regions < 1:
        return X, []
    else:
        pothole_covered = X.copy()
        pothole_regions = ndimage.find_objects(abnormal_fm)
        spikeish_indices = get_all_spikeish_indices(pothole_regions)
        for each_region in pothole_regions:
            region = each_region[0]
            start, stop = region.start, region.stop
            region_length = stop-start
            # find the sample values next to the start and stop

            next_to_start = find_non_forbidden_index(start, spikeish_indices, -1, X)
            next_to_stop =  find_non_forbidden_index(stop, spikeish_indices, +1, X)
            interpolated_values = np.linspace(pothole_covered[next_to_start],
                                              pothole_covered[next_to_stop],
                                              region_length)
            pothole_covered[start:stop] = interpolated_values

    return pothole_covered, pothole_regions

def identify_pothole_samples(X, fs,  **kwargs):
    '''Moves a sliding window and checks the values of samples in the sliding window.
    If the jump of values between samples is not linearly propotional to the
    expected max_stepsize, then it is labelled a pothole sample.

    A pothole sample is one which represents a sudden jump in the values - indicating
    a noisy tracking of the frequency. The jump in values in a non-noisy signal is expected
    to be proportional to the distance between the samples.

    For instance, if :

    >>> a = np.array([10, 2, 6, 10, 12])

    If the max step size is 2, then because abs(10-2)>2, it causes a pothole to appear on 2.
    There is no pothole label on the 2nd index because abs(10-6) is not >4. Because 10 and 6
    are two samples apart, the maximum allowed jump in value is max_stepsize*2, which is 4.

    For optimal pothole detection the 'look-ahead' span of the pothole_inspection_window
    should at least the size of the longest expected potholes. Smaller window sizes
    will lead to false negatives.

    Parameters
    ----------
    X : np.array
    fs : float>0
    max_stepsize : float>0
        The max absolute difference between the values of one sample to the next.

    pothole_inspection_window : float>0, optional
        Defaults to 0.25ms

    Returns
    -------
    pothole_candidates : np.array
        Boolean array with same size as X. Sample that are True represent pothole candidates.

    See Also
    --------
    detect_local_potholes
    '''
    # forward pass
    left2right_pothole_candidates, _ = onepass_identify_potholes(X, fs,
                                                                 **kwargs)
    # backward pass
    right2left_pothole_candidates, _ = onepass_identify_potholes(X[::-1], fs,
                                                              **kwargs)
    pothole_candidates = np.logical_and(left2right_pothole_candidates>0,
                                        right2left_pothole_candidates[::-1]>0)

    return pothole_candidates

def onepass_identify_potholes(X, fs, max_stepsize, **kwargs):
    '''
    '''
    window_duration = kwargs.get('pothole_inspection_window', 0.25*10**-3)
    window_size = int(fs*window_duration)
    potholes = []
    consensus = np.zeros(X.size)

    for i,each in enumerate(X):
        candidates = detect_local_potholes(X[i:i+window_size], max_stepsize)
        potholes.append(candidates)
        consensus[i+candidates] += 1
    return consensus, potholes

def detect_local_potholes(X, max_step_size):
    '''accepts a 1D array and checks the absolute difference between
    the first sample and all other samples.

    The samples with difference greater than the linearly expected increase
    from max_step_sizes are labelled candidate potholes.

    Parameters
    ----------
    X : np.array
    max_step_size : float>=0

    Returns
    -------
    candidate_potholes : np.array
        Boolean array of same size as X
    '''
    pothole_depth = np.abs(X-X[0])
    max_allowed = np.arange(X.size)*max_step_size
    candidate_potholes = np.argwhere(pothole_depth > max_allowed).flatten()
    return candidate_potholes

def get_all_spikeish_indices(regions):
    '''
    '''
    indices = np.concatenate(list(map(extract_indices, regions)))
    return indices

def find_non_forbidden_index(candidate, forbidden_indices, search_direction, X):
    '''
    '''

    index_in_spike = candidate in forbidden_indices
    if index_in_spike:
        next_candidate = search_operation[search_direction](candidate)
        return find_non_forbidden_index(next_candidate, forbidden_indices, search_direction, X)
    else:
        candidate_index_within_array = np.logical_and(candidate>=0, candidate<=X.size-1)
        if candidate_index_within_array:
            return candidate
        else:
            search_direction *= -1
            next_candidate = search_operation[search_direction](candidate)
            return find_non_forbidden_index(next_candidate, forbidden_indices, search_direction, X)

def extract_indices(X):
    indices = np.arange(X[0].start, X[0].stop)
    return indices

search_operation = {-1: lambda X: X-1,
                    1: lambda X : X+1}

def remove_bursts(X, fs, **kwargs):
    '''Bursts are brief but large jumps in the signal above zero. Even though they satisfy
    most of the other conditions of beginning above the noise floor and of
    being above 0 value, they still are too short to be relevant signals.

    Parameters
    ----------
    X : np.array
        The noisy signal to be handled
    fs : float>0
        Sampling rate in Hz.
    min_element_length : float>0, optional
        The minimum length a section of the signal must be to be
        kept in seconds. Defaults to 5 inter-sample-intervals.

    Returns
    -------
    X_nonspikey : np.array
        Same size as X, and without very short segments.

    See Also
    --------
    segments_above_min_duration

    Notes
    -----
    An inter-sample-interval is defined as 1/fs

    '''
    inter_sample_durn = 1.0/fs
    min_element_length = kwargs.get('min_element_length', 5*inter_sample_durn) #to 5 samples
    min_element_samples = int(fs*min_element_length)

    if  min_element_length <= inter_sample_durn:
        raise ValueError('Please set the min element length.\
        The current value of:%f is less than 1/sampling rate'%(min_element_length))
    min_element_samples = int(fs*min_element_length)

    non_spikey_regions = segments_above_min_duration(X>0, min_element_samples)

    X_nonspikey = np.zeros(X.size)
    X_nonspikey[non_spikey_regions] = X[non_spikey_regions]
    return X_nonspikey

def segments_above_min_duration(satisfies_condition, min_samples):
    '''Accepts a boolean array and looks for continuous chunks
    that are above a minimum length.

    Parameters
    ----------
    satisfies_condition : np.array
        Boolean array where samples with True satisfy a condition.
    min_samples : int >0
        The minimum number of samples a continuous region of True
        must be to be kept.

    Returns
    -------
    above_min_duration : np.array
        Same size as satisfies_condition, with only the continuous
        chunks that are above min_samples.

    '''
    all_regions, number_regions = ndimage.label(satisfies_condition)
    region_stretches = ndimage.find_objects(all_regions)

    above_min_duration = np.tile(False, satisfies_condition.size)

    for each_stretch in region_stretches:
        if satisfies_condition[each_stretch].size > min_samples:
            above_min_duration[each_stretch] = True
    return above_min_duration

def suppress_background_noise(main_signal, input_audio, **kwargs):
    '''
    '''
    signal_level = kwargs.get('signal_level', -20) # dBrms
    signal_dBrms = dB(moving_rms_edge_robust(input_audio, **kwargs))
    bg_noise_suppressed = suppress_to_zero(main_signal, signal_dBrms,
                                           signal_level, 'below')
    return bg_noise_suppressed

def suppress_frequency_spikes(noisy_profile, input_audio, fs, **kwargs):
    '''
    '''
    max_spike_rate = kwargs.get('max_spike_rate', 3000) # Hz jump/sample

    # median filter to get rid of smaller fluctuations in the noisy profile *not*
    # caused by abrupt transitions in the edges.
    med_filtered = median_filter(noisy_profile, fs, **kwargs)

    raw_fmrate = abs(np.diff(med_filtered))
    delta_profile = resize_by_adding_one_sample(raw_fmrate, input_audio)
    spike_suppressed = suppress_to_zero(noisy_profile, delta_profile, max_spike_rate, 'above')
    return spike_suppressed

def suppress_to_zero(target_signal, basis_signal, threshold, mode='below'):
    '''
    Sets the values of the target signal to zero if the
    samples in the basis_signal are $\geq$ or $\leq$ the threshold

    Parameters
    ----------
    target_signal, basis_signal : np.array
    threshold : float
    mode : ['below', 'above'], str

    Returns
    -------
    cleaned_signal : np.array
        A copy of the target signal with the values that are below/above the threshold
        set to zero

    Example
    --------
    # create a basis signal with a 'weak' left half and a 'loud' right hald
    # we want to suppress the we
    >>> basis = np.concatenate((np.arange(10), np.arange(100,200)))
    >>> target_signal = np.random.normal(0,1,basis.size)
    >>> cleaned_target = suppress_to_zero(basis, target_signal, 100, mode='above')
    '''
    if mode == 'below':
        to_suppress = basis_signal < threshold
    elif mode == 'above':
        to_suppress = basis_signal > threshold
    else:
        raise ValueError('Mode should be either "below" or "above" and not: %s'%(mode))
    cleaned_signal = np.copy(target_signal)
    cleaned_signal[to_suppress.flatten()] = 0
    return cleaned_signal

def clip_tfr(tfr, **kwargs):
    '''
    Parameters
    ----------
    tfr : np.array
        2D array with the time-frequency representation of choice
        (pwvd, fft etc). The tfr must have real-valued non-negative
        values as the clip range is defined in dB.
    tfr_cliprange: float >0, optional
        The maximum dynamic range in dB which will be used to
        track the instantaneous frequency. Defaults to
        None. See `Notes` for more details

    Returns
    -------
    clipped_tfr : np.array
        A 2d array of same shape as `tfr`, with values
        clipped between [max, max x 10^(tfr_range/20)]
    Notes
    -----
    The `tfr_cliprange` is used to remove the presence of
    background noise, faint harmonics or revernberations/echoes
    in the audio. This of course all assumes that the main
    signal itself is sufficiently intense in the first place.

    After the PWVD time-frequency represenation is made,
    values below X dB of the maximum value are 'clipped' to
    the same minimum value. eg. if the pwvd had values of
    [0.1, 0.9, 0.3, 1, 0.001, 0.0006] and the tfr_cliprange is
    set to  6dB, then the output of the clipping will be
    [0.5, 0.9, 0.3, 1, 0.5, 0.5].  This step essentially eliminates
    any variation in the array, thus allowing a clear
    tracking of the highest component in it.
    '''
    tfr_cliprange = kwargs.get('tfr_cliprange')
    if tfr_cliprange is None:
        return tfr
    else:
        make_sure_its_positive(tfr_cliprange, variable='tfr_cliprange')
        max_value = np.max(tfr)
        clip_value = max_value*10**(-tfr_cliprange/20.0)

        clipped_tfr = tfr.copy()
        clipped_tfr[clipped_tfr<clip_value] = clip_value

        return clipped_tfr

def conditionally_set_to(X, conditional, bool_state):
    '''Inverts the samples in X where the conditional is True.
    Parameters
    ----------
    X : np.array
        Boolean
    conditional : np.array
        Boolean
    bool_state : [True, False]

    Returns
    -------
    cond_set_X : np.array
        conditionally set X

    Notes
    -----
    this function is useful if you want to 'suppress' a few samples
    conditionally based ont he values of the same samples
    on another array.

    Example
    -------
    >>> x = np.array([True, True, False, False, True])
    >>> y = np.array([0,0,10,10,10])
    Imagine x is some kind of detection array, while y is the
    signal-to-noise ratio at each of the sample points. Of course,
    you'd like to discard all the predictions from low SNR measurements.
    Let's say you want to keep only those entries in X where y is >1.
    >>> x_cond = conditionally_set_to(x, y<10, False)
    >>> x_cond

    np.array([False, False, False, False, True ])
    '''
    if sum(conditional)==0:
        return X
    else:
        cond_set_X = X.copy()
        cond_set_X[conditional] = bool_state
        return cond_set_X

max_stars_repo_name: thejasvibr/itsfm

reason: import numpy,from scipy


text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst:

"""
Tests for model evaluation.
Compare the results of some models with other programs.
"""
import unittest.mock as mk

import numpy as np
# pylint: disable=invalid-name, no-member
import pytest
from numpy.testing import assert_allclose, assert_equal

import astropy.modeling.tabular as tabular_models
from astropy import units as u
from astropy.modeling import fitting, models
from astropy.modeling.bounding_box import ModelBoundingBox
from astropy.modeling.core import FittableModel, Model, _ModelMeta
from astropy.modeling.models import Gaussian2D
from astropy.modeling.parameters import InputParameterError, Parameter
from astropy.modeling.polynomial import PolynomialBase
from astropy.modeling.powerlaws import (
    BrokenPowerLaw1D, ExponentialCutoffPowerLaw1D, LogParabola1D, PowerLaw1D,
    SmoothlyBrokenPowerLaw1D)
from astropy.modeling.separable import separability_matrix
from astropy.tests.helper import assert_quantity_allclose
from astropy.utils import NumpyRNGContext
from astropy.utils.compat.optional_deps import HAS_SCIPY  # noqa: F401

from .example_models import models_1D, models_2D

fitters = [
    fitting.LevMarLSQFitter,
    fitting.TRFLSQFitter,
    fitting.LMLSQFitter,
    fitting.DogBoxLSQFitter
]

@pytest.mark.skipif('not HAS_SCIPY')
@pytest.mark.parametrize('fitter', fitters)
def test_custom_model(fitter, amplitude=4, frequency=1):
    fitter = fitter()

    def sine_model(x, amplitude=4, frequency=1):
        """
        Model function
        """
        return amplitude * np.sin(2 * np.pi * frequency * x)

    def sine_deriv(x, amplitude=4, frequency=1):
        """
        Jacobian of model function, e.g. derivative of the function with
        respect to the *parameters*
        """
        da = np.sin(2 * np.pi * frequency * x)
        df = 2 * np.pi * x * amplitude * np.cos(2 * np.pi * frequency * x)
        return np.vstack((da, df))

    SineModel = models.custom_model(sine_model, fit_deriv=sine_deriv)

    x = np.linspace(0, 4, 50)
    sin_model = SineModel()

    sin_model.evaluate(x, 5., 2.)
    sin_model.fit_deriv(x, 5., 2.)

    np.random.seed(0)
    data = sin_model(x) + np.random.rand(len(x)) - 0.5
    model = fitter(sin_model, x, data)
    assert np.all((np.array([model.amplitude.value, model.frequency.value]) -
                   np.array([amplitude, frequency])) < 0.001)

def test_custom_model_init():
    @models.custom_model
    def SineModel(x, amplitude=4, frequency=1):
        """Model function"""

        return amplitude * np.sin(2 * np.pi * frequency * x)

    sin_model = SineModel(amplitude=2., frequency=0.5)
    assert sin_model.amplitude == 2.
    assert sin_model.frequency == 0.5

def test_custom_model_defaults():
    @models.custom_model
    def SineModel(x, amplitude=4, frequency=1):
        """Model function"""

        return amplitude * np.sin(2 * np.pi * frequency * x)

    sin_model = SineModel()
    assert SineModel.amplitude.default == 4
    assert SineModel.frequency.default == 1

    assert sin_model.amplitude == 4
    assert sin_model.frequency == 1

def test_inconsistent_input_shapes():
    g = Gaussian2D()
    x = np.arange(-1., 1, .2)
    y = x.copy()
    # check scalar input broadcasting works
    assert np.abs(g(x, 0) - g(x, 0 * x)).sum() == 0
    # but not array broadcasting
    x.shape = (10, 1)
    y.shape = (1, 10)
    result = g(x, y)
    assert result.shape == (10, 10)

def test_custom_model_bounding_box():
    """Test bounding box evaluation for a 3D model"""

    def ellipsoid(x, y, z, x0=13, y0=10, z0=8, a=4, b=3, c=2, amp=1):
        rsq = ((x - x0) / a) ** 2 + ((y - y0) / b) ** 2 + ((z - z0) / c) ** 2
        val = (rsq < 1) * amp
        return val

    class Ellipsoid3D(models.custom_model(ellipsoid)):
        @property
        def bounding_box(self):
            return ((self.z0 - self.c, self.z0 + self.c),
                    (self.y0 - self.b, self.y0 + self.b),
                    (self.x0 - self.a, self.x0 + self.a))

    model = Ellipsoid3D()
    bbox = model.bounding_box

    zlim, ylim, xlim = bbox.bounding_box()
    dz, dy, dx = np.diff(bbox) / 2
    z1, y1, x1 = np.mgrid[slice(zlim[0], zlim[1] + 1),
                          slice(ylim[0], ylim[1] + 1),
                          slice(xlim[0], xlim[1] + 1)]
    z2, y2, x2 = np.mgrid[slice(zlim[0] - dz, zlim[1] + dz + 1),
                          slice(ylim[0] - dy, ylim[1] + dy + 1),
                          slice(xlim[0] - dx, xlim[1] + dx + 1)]

    arr = model(x2, y2, z2, with_bounding_box=True)
    sub_arr = model(x1, y1, z1, with_bounding_box=True)

    # check for flux agreement
    assert abs(np.nansum(arr) - np.nansum(sub_arr)) < np.nansum(arr) * 1e-7

class Fittable2DModelTester:
    """
    Test class for all two dimensional parametric models.

    Test values have to be defined in example_models.py. It currently test the
    model with different input types, evaluates the model at different
    positions and assures that it gives the correct values. And tests if the
    model works with non-linear fitters.

    This can be used as a base class for user defined model testing.
    """

    def setup_class(self):
        self.N = 100
        self.M = 100
        self.eval_error = 0.0001
        self.fit_error = 0.1
        self.x = 5.3
        self.y = 6.7
        self.x1 = np.arange(1, 10, .1)
        self.y1 = np.arange(1, 10, .1)
        self.y2, self.x2 = np.mgrid[:10, :8]

    def test_input2D(self, model_class, test_parameters):
        """Test model with different input types."""

        model = create_model(model_class, test_parameters)
        model(self.x, self.y)
        model(self.x1, self.y1)
        model(self.x2, self.y2)

    def test_eval2D(self, model_class, test_parameters):
        """Test model values add certain given points"""

        model = create_model(model_class, test_parameters)
        x = test_parameters['x_values']
        y = test_parameters['y_values']
        z = test_parameters['z_values']
        assert np.all(np.abs(model(x, y) - z) < self.eval_error)

    def test_bounding_box2D(self, model_class, test_parameters):
        """Test bounding box evaluation"""

        model = create_model(model_class, test_parameters)

        # testing setter
        model.bounding_box = ((-5, 5), (-5, 5))
        assert model.bounding_box == ((-5, 5), (-5, 5))

        model.bounding_box = None
        with pytest.raises(NotImplementedError):
            model.bounding_box

        # test the exception of dimensions don't match
        with pytest.raises(ValueError):
            model.bounding_box = (-5, 5)

        del model.bounding_box

        try:
            bbox = model.bounding_box
        except NotImplementedError:
            return

        ddx = 0.01
        ylim, xlim = bbox
        x1 = np.arange(xlim[0], xlim[1], ddx)
        y1 = np.arange(ylim[0], ylim[1], ddx)

        x2 = np.concatenate(([xlim[0] - idx * ddx for idx in range(10, 0, -1)],
                             x1,
                             [xlim[1] + idx * ddx for idx in range(1, 10)]))
        y2 = np.concatenate(([ylim[0] - idx * ddx for idx in range(10, 0, -1)],
                             y1,
                             [ylim[1] + idx * ddx for idx in range(1, 10)]))

        inside_bbox = model(x1, y1)
        outside_bbox = model(x2, y2, with_bounding_box=True)
        outside_bbox = outside_bbox[~np.isnan(outside_bbox)]

        assert np.all(inside_bbox == outside_bbox)

    def test_bounding_box2D_peak(self, model_class, test_parameters):
        if not test_parameters.pop('bbox_peak', False):
            return

        model = create_model(model_class, test_parameters)
        bbox = model.bounding_box

        ylim, xlim = bbox
        dy, dx = np.diff(bbox)/2
        y1, x1 = np.mgrid[slice(ylim[0], ylim[1] + 1),
                          slice(xlim[0], xlim[1] + 1)]
        y2, x2 = np.mgrid[slice(ylim[0] - dy, ylim[1] + dy + 1),
                          slice(xlim[0] - dx, xlim[1] + dx + 1)]

        arr = model(x2, y2)
        sub_arr = model(x1, y1)

        # check for flux agreement
        assert abs(arr.sum() - sub_arr.sum()) < arr.sum() * 1e-7

    @pytest.mark.skipif('not HAS_SCIPY')
    @pytest.mark.parametrize('fitter', fitters)
    def test_fitter2D(self, model_class, test_parameters, fitter):
        """Test if the parametric model works with the fitter."""
        fitter = fitter()

        x_lim = test_parameters['x_lim']
        y_lim = test_parameters['y_lim']

        parameters = test_parameters['parameters']
        model = create_model(model_class, test_parameters)

        if isinstance(parameters, dict):
            parameters = [parameters[name] for name in model.param_names]

        if "log_fit" in test_parameters:
            if test_parameters['log_fit']:
                x = np.logspace(x_lim[0], x_lim[1], self.N)
                y = np.logspace(y_lim[0], y_lim[1], self.N)
        else:
            x = np.linspace(x_lim[0], x_lim[1], self.N)
            y = np.linspace(y_lim[0], y_lim[1], self.N)
        xv, yv = np.meshgrid(x, y)

        np.random.seed(0)
        # add 10% noise to the amplitude
        noise = np.random.rand(self.N, self.N) - 0.5
        data = model(xv, yv) + 0.1 * parameters[0] * noise
        new_model = fitter(model, xv, yv, data)

        params = [getattr(new_model, name) for name in new_model.param_names]
        fixed = [param.fixed for param in params]
        expected = np.array([val for val, fixed in zip(parameters, fixed)
                             if not fixed])
        fitted = np.array([param.value for param in params
                           if not param.fixed])
        assert_allclose(fitted, expected,
                        atol=self.fit_error)

    @pytest.mark.skipif('not HAS_SCIPY')
    @pytest.mark.parametrize('fitter', fitters)
    def test_deriv_2D(self, model_class, test_parameters, fitter):
        """
        Test the derivative of a model by fitting with an estimated and
        analytical derivative.
        """
        fitter = fitter()

        x_lim = test_parameters['x_lim']
        y_lim = test_parameters['y_lim']

        if model_class.fit_deriv is None or issubclass(model_class, PolynomialBase):
            return

        if "log_fit" in test_parameters:
            if test_parameters['log_fit']:
                x = np.logspace(x_lim[0], x_lim[1], self.N)
                y = np.logspace(y_lim[0], y_lim[1], self.M)
                x_test = np.logspace(x_lim[0], x_lim[1], self.N*10)
                y_test = np.logspace(y_lim[0], y_lim[1], self.M*10)
        else:
            x = np.linspace(x_lim[0], x_lim[1], self.N)
            y = np.linspace(y_lim[0], y_lim[1], self.M)
            x_test = np.linspace(x_lim[0], x_lim[1], self.N*10)
            y_test = np.linspace(y_lim[0], y_lim[1], self.M*10)
        xv, yv = np.meshgrid(x, y)
        xv_test, yv_test = np.meshgrid(x_test, y_test)

        try:
            model_with_deriv = create_model(model_class, test_parameters,
                                            use_constraints=False,
                                            parameter_key='deriv_initial')
            model_no_deriv = create_model(model_class, test_parameters,
                                          use_constraints=False,
                                          parameter_key='deriv_initial')
            model = create_model(model_class, test_parameters,
                                 use_constraints=False,
                                 parameter_key='deriv_initial')
        except KeyError:
            model_with_deriv = create_model(model_class, test_parameters,
                                            use_constraints=False)
            model_no_deriv = create_model(model_class, test_parameters,
                                          use_constraints=False)
            model = create_model(model_class, test_parameters,
                                 use_constraints=False)

        # add 10% noise to the amplitude
        rsn = np.random.default_rng(0)
        amplitude = test_parameters['parameters'][0]
        n = 0.1 * amplitude * (rsn.random((self.M, self.N)) - 0.5)

        data = model(xv, yv) + n
        fitter_with_deriv = fitter
        new_model_with_deriv = fitter_with_deriv(model_with_deriv, xv, yv,
                                                 data)
        fitter_no_deriv = fitter
        new_model_no_deriv = fitter_no_deriv(model_no_deriv, xv, yv, data,
                                             estimate_jacobian=True)
        assert_allclose(new_model_with_deriv(xv_test, yv_test),
                        new_model_no_deriv(xv_test, yv_test),
                        rtol=1e-2)
        if model_class != Gaussian2D:
            assert_allclose(new_model_with_deriv.parameters,
                            new_model_no_deriv.parameters,
                            rtol=0.1)

class Fittable1DModelTester:
    """
    Test class for all one dimensional parametric models.

    Test values have to be defined in example_models.py. It currently test the
    model with different input types, evaluates the model at different
    positions and assures that it gives the correct values. And tests if the
    model works with non-linear fitters.

    This can be used as a base class for user defined model testing.
    """

    # These models will fail fitting test, because built in fitting data
    #   will produce non-finite values
    _non_finite_models = [
        BrokenPowerLaw1D,
        ExponentialCutoffPowerLaw1D,
        LogParabola1D,
        PowerLaw1D,
        SmoothlyBrokenPowerLaw1D
    ]

    def setup_class(self):
        self.N = 100
        self.M = 100
        self.eval_error = 0.0001
        self.fit_error = 0.11
        self.x = 5.3
        self.y = 6.7
        self.x1 = np.arange(1, 10, .1)
        self.y1 = np.arange(1, 10, .1)
        self.y2, self.x2 = np.mgrid[:10, :8]

    @pytest.mark.filterwarnings(r'ignore:.*:RuntimeWarning')
    def test_input1D(self, model_class, test_parameters):
        """Test model with different input types."""

        model = create_model(model_class, test_parameters)
        model(self.x)
        model(self.x1)
        model(self.x2)

    def test_eval1D(self, model_class, test_parameters):
        """
        Test model values at certain given points
        """
        model = create_model(model_class, test_parameters)
        x = test_parameters['x_values']
        y = test_parameters['y_values']
        assert_allclose(model(x), y, atol=self.eval_error)

    def test_bounding_box1D(self, model_class, test_parameters):
        """Test bounding box evaluation"""

        model = create_model(model_class, test_parameters)

        # testing setter
        model.bounding_box = (-5, 5)
        model.bounding_box = None

        with pytest.raises(NotImplementedError):
            model.bounding_box

        del model.bounding_box

        # test exception if dimensions don't match
        with pytest.raises(ValueError):
            model.bounding_box = 5

        try:
            bbox = model.bounding_box.bounding_box()
        except NotImplementedError:
            return

        ddx = 0.01
        x1 = np.arange(bbox[0], bbox[1], ddx)
        x2 = np.concatenate(([bbox[0] - idx * ddx for idx in range(10, 0, -1)],
                             x1,
                             [bbox[1] + idx * ddx for idx in range(1, 10)]))

        inside_bbox = model(x1)
        outside_bbox = model(x2, with_bounding_box=True)
        outside_bbox = outside_bbox[~np.isnan(outside_bbox)]

        assert np.all(inside_bbox == outside_bbox)

    def test_bounding_box1D_peak(self, model_class, test_parameters):
        if not test_parameters.pop('bbox_peak', False):
            return

        model = create_model(model_class, test_parameters)
        bbox = model.bounding_box

        if isinstance(model, models.Lorentz1D) or isinstance(model, models.Drude1D):
            rtol = 0.01  # 1% agreement is enough due to very extended wings
            ddx = 0.1  # Finer sampling to "integrate" flux for narrow peak
        else:
            rtol = 1e-7
            ddx = 1

        if isinstance(bbox, ModelBoundingBox):
            bbox = bbox.bounding_box()

        dx = np.diff(bbox) / 2
        x1 = np.mgrid[slice(bbox[0], bbox[1] + 1, ddx)]
        x2 = np.mgrid[slice(bbox[0] - dx, bbox[1] + dx + 1, ddx)]
        arr = model(x2)
        sub_arr = model(x1)

        # check for flux agreement
        assert abs(arr.sum() - sub_arr.sum()) < arr.sum() * rtol

    @pytest.mark.skipif('not HAS_SCIPY')
    @pytest.mark.parametrize('fitter', fitters)
    def test_fitter1D(self, model_class, test_parameters, fitter):
        """
        Test if the parametric model works with the fitter.
        """
        fitter = fitter()

        x_lim = test_parameters['x_lim']
        parameters = test_parameters['parameters']
        model = create_model(model_class, test_parameters)

        if isinstance(parameters, dict):
            parameters = [parameters[name] for name in model.param_names]

        if "log_fit" in test_parameters:
            if test_parameters['log_fit']:
                x = np.logspace(x_lim[0], x_lim[1], self.N)
        else:
            x = np.linspace(x_lim[0], x_lim[1], self.N)

        np.random.seed(0)
        # add 10% noise to the amplitude
        relative_noise_amplitude = 0.01
        data = ((1 + relative_noise_amplitude * np.random.randn(len(x))) *
                model(x))
        new_model = fitter(model, x, data)

        # Only check parameters that were free in the fit
        params = [getattr(new_model, name) for name in new_model.param_names]
        fixed = [param.fixed for param in params]
        expected = np.array([val for val, fixed in zip(parameters, fixed)
                             if not fixed])
        fitted = np.array([param.value for param in params
                           if not param.fixed])
        if model_class == models.BrokenPowerLaw1D:
            atol = 3.75
        else:
            atol = self.fit_error
        assert_allclose(fitted, expected, atol=atol)

    @pytest.mark.skipif('not HAS_SCIPY')
    @pytest.mark.filterwarnings(r'ignore:.*:RuntimeWarning')
    @pytest.mark.parametrize('fitter', fitters)
    def test_deriv_1D(self, model_class, test_parameters, fitter):
        """
        Test the derivative of a model by comparing results with an estimated
        derivative.
        """
        fitter = fitter()

        if model_class in self._non_finite_models:
            return

        x_lim = test_parameters['x_lim']

        if model_class.fit_deriv is None or issubclass(model_class, PolynomialBase):
            return

        if "log_fit" in test_parameters:
            if test_parameters['log_fit']:
                x = np.logspace(x_lim[0], x_lim[1], self.N)
        else:
            x = np.linspace(x_lim[0], x_lim[1], self.N)

        parameters = test_parameters['parameters']
        model_with_deriv = create_model(model_class, test_parameters,
                                        use_constraints=False)
        model_no_deriv = create_model(model_class, test_parameters,
                                      use_constraints=False)

        # NOTE: PR 10644 replaced deprecated usage of RandomState but could not
        #       find a new seed that did not cause test failure, resorted to hardcoding.
        # add 10% noise to the amplitude
        rsn_rand_1234567890 = np.array([
            0.61879477, 0.59162363, 0.88868359, 0.89165480, 0.45756748,
            0.77818808, 0.26706377, 0.99610621, 0.54009489, 0.53752161,
            0.40099938, 0.70540579, 0.40518559, 0.94999075, 0.03075388,
            0.13602495, 0.08297726, 0.42352224, 0.23449723, 0.74743526,
            0.65177865, 0.68998682, 0.16413419, 0.87642114, 0.44733314,
            0.57871104, 0.52377835, 0.62689056, 0.34869427, 0.26209748,
            0.07498055, 0.17940570, 0.82999425, 0.98759822, 0.11326099,
            0.63846415, 0.73056694, 0.88321124, 0.52721004, 0.66487673,
            0.74209309, 0.94083846, 0.70123128, 0.29534353, 0.76134369,
            0.77593881, 0.36985514, 0.89519067, 0.33082813, 0.86108824,
            0.76897859, 0.61343376, 0.43870907, 0.91913538, 0.76958966,
            0.51063556, 0.04443249, 0.57463611, 0.31382006, 0.41221713,
            0.21531811, 0.03237521, 0.04166386, 0.73109303, 0.74556052,
            0.64716325, 0.77575353, 0.64599254, 0.16885816, 0.48485480,
            0.53844248, 0.99690349, 0.23657074, 0.04119088, 0.46501519,
            0.35739006, 0.23002665, 0.53420791, 0.71639475, 0.81857486,
            0.73994342, 0.07948837, 0.75688276, 0.13240193, 0.48465576,
            0.20624753, 0.02298276, 0.54257873, 0.68123230, 0.35887468,
            0.36296147, 0.67368397, 0.29505730, 0.66558885, 0.93652252,
            0.36755130, 0.91787687, 0.75922703, 0.48668067, 0.45967890])
        n = 0.1 * parameters[0] * (rsn_rand_1234567890 - 0.5)

        data = model_with_deriv(x) + n
        fitter_with_deriv = fitter
        new_model_with_deriv = fitter_with_deriv(model_with_deriv, x, data)
        fitter_no_deriv = fitter
        new_model_no_deriv = fitter_no_deriv(model_no_deriv, x, data,
                                             estimate_jacobian=True)
        assert_allclose(new_model_with_deriv.parameters,
                        new_model_no_deriv.parameters, atol=0.15)

def create_model(model_class, test_parameters, use_constraints=True,
                 parameter_key='parameters'):
    """Create instance of model class."""

    constraints = {}
    if issubclass(model_class, PolynomialBase):
        return model_class(**test_parameters[parameter_key])
    elif issubclass(model_class, FittableModel):
        if "requires_scipy" in test_parameters and not HAS_SCIPY:
            pytest.skip("SciPy not found")
        if use_constraints:
            if 'constraints' in test_parameters:
                constraints = test_parameters['constraints']
        return model_class(*test_parameters[parameter_key], **constraints)

@pytest.mark.filterwarnings(r'ignore:Model is linear in parameters.*')
@pytest.mark.filterwarnings(r'ignore:The fit may be unsuccessful.*')
@pytest.mark.parametrize(('model_class', 'test_parameters'),
                         sorted(models_1D.items(), key=lambda x: str(x[0])))
class TestFittable1DModels(Fittable1DModelTester):
    pass

@pytest.mark.filterwarnings(r'ignore:Model is linear in parameters.*')
@pytest.mark.parametrize(('model_class', 'test_parameters'),
                         sorted(models_2D.items(), key=lambda x: str(x[0])))
class TestFittable2DModels(Fittable2DModelTester):
    pass

def test_ShiftModel():
    # Shift by a scalar
    m = models.Shift(42)
    assert m(0) == 42
    assert_equal(m([1, 2]), [43, 44])

    # Shift by a list
    m = models.Shift([42, 43], n_models=2)
    assert_equal(m(0), [42, 43])
    assert_equal(m([1, 2], model_set_axis=False),
                 [[43, 44], [44, 45]])

def test_ScaleModel():
    # Scale by a scalar
    m = models.Scale(42)
    assert m(0) == 0
    assert_equal(m([1, 2]), [42, 84])

    # Scale by a list
    m = models.Scale([42, 43], n_models=2)
    assert_equal(m(0), [0, 0])
    assert_equal(m([1, 2], model_set_axis=False),
                 [[42, 84], [43, 86]])

def test_voigt_model():
    """
    Currently just tests that the model peaks at its origin.
    Regression test for https://github.com/astropy/astropy/issues/3942
    """

    m = models.Voigt1D(x_0=5, amplitude_L=10, fwhm_L=0.5, fwhm_G=0.9)
    x = np.arange(0, 10, 0.01)
    y = m(x)
    assert y[500] == y.max()  # y[500] is right at the center

def test_model_instance_repr():
    m = models.Gaussian1D(1.5, 2.5, 3.5)
    assert repr(m) == '<Gaussian1D(amplitude=1.5, mean=2.5, stddev=3.5)>'

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_interp_1d():
    """
    Test Tabular1D model.
    """
    points = np.arange(0, 5)
    values = [1., 10, 2, 45, -3]
    LookupTable = models.tabular_model(1)
    model = LookupTable(points=points, lookup_table=values)
    xnew = [0., .7, 1.4, 2.1, 3.9]
    ans1 = [1., 7.3, 6.8, 6.3, 1.8]
    assert_allclose(model(xnew), ans1)
    # Test evaluate without passing `points`.
    model = LookupTable(lookup_table=values)
    assert_allclose(model(xnew), ans1)
    # Test bounds error.
    xextrap = [0., .7, 1.4, 2.1, 3.9, 4.1]
    with pytest.raises(ValueError):
        model(xextrap)
    # test extrapolation and fill value
    model = LookupTable(lookup_table=values, bounds_error=False,
                        fill_value=None)
    assert_allclose(model(xextrap),
                    [1., 7.3, 6.8, 6.3, 1.8, -7.8])

    # Test unit support
    xnew = xnew * u.nm
    ans1 = ans1 * u.nJy
    model = LookupTable(points=points*u.nm, lookup_table=values*u.nJy)
    assert_quantity_allclose(model(xnew), ans1)
    assert_quantity_allclose(model(xnew.to(u.nm)), ans1)
    assert model.bounding_box == (0 * u.nm, 4 * u.nm)

    # Test fill value unit conversion and unitless input on table with unit
    model = LookupTable([1, 2, 3], [10, 20, 30] * u.nJy, bounds_error=False,
                        fill_value=1e-33*(u.W / (u.m * u.m * u.Hz)))
    assert_quantity_allclose(model(np.arange(5)),
                             [100, 10, 20, 30, 100] * u.nJy)

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_interp_2d():
    table = np.array([
        [-0.04614432, -0.02512547, -0.00619557, 0.0144165, 0.0297525],
        [-0.04510594, -0.03183369, -0.01118008, 0.01201388, 0.02496205],
        [-0.05464094, -0.02804499, -0.00960086, 0.01134333, 0.02284104],
        [-0.04879338, -0.02539565, -0.00440462, 0.01795145, 0.02122417],
        [-0.03637372, -0.01630025, -0.00157902, 0.01649774, 0.01952131]])

    points = np.arange(0, 5)
    points = (points, points)

    xnew = np.array([0., .7, 1.4, 2.1, 3.9])
    LookupTable = models.tabular_model(2)
    model = LookupTable(points, table)
    znew = model(xnew, xnew)
    result = np.array(
        [-0.04614432, -0.03450009, -0.02241028, -0.0069727, 0.01938675])
    assert_allclose(znew, result, atol=1e-7)

    # test 2D arrays as input
    a = np.arange(12).reshape((3, 4))
    y, x = np.mgrid[:3, :4]
    t = models.Tabular2D(lookup_table=a)
    r = t(y, x)
    assert_allclose(a, r)

    with pytest.raises(ValueError):
        model = LookupTable(points=([1.2, 2.3], [1.2, 6.7], [3, 4]))
    with pytest.raises(ValueError):
        model = LookupTable(lookup_table=[1, 2, 3])
    with pytest.raises(NotImplementedError):
        model = LookupTable(n_models=2)
    with pytest.raises(ValueError):
        model = LookupTable(([1, 2], [3, 4]), [5, 6])
    with pytest.raises(ValueError):
        model = LookupTable(([1, 2] * u.m, [3, 4]), [[5, 6], [7, 8]])
    with pytest.raises(ValueError):
        model = LookupTable(points, table, bounds_error=False,
                            fill_value=1*u.Jy)

    # Test unit support
    points = points[0] * u.nm
    points = (points, points)
    xnew = xnew * u.nm
    model = LookupTable(points, table * u.nJy)
    result = result * u.nJy
    assert_quantity_allclose(model(xnew, xnew), result, atol=1e-7*u.nJy)
    xnew = xnew.to(u.m)
    assert_quantity_allclose(model(xnew, xnew), result, atol=1e-7*u.nJy)
    bbox = (0 * u.nm, 4 * u.nm)
    bbox = (bbox, bbox)
    assert model.bounding_box == bbox

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_nd():
    a = np.arange(24).reshape((2, 3, 4))
    x, y, z = np.mgrid[:2, :3, :4]
    tab = models.tabular_model(3)
    t = tab(lookup_table=a)
    result = t(x, y, z)
    assert_allclose(a, result)

    with pytest.raises(ValueError):
        models.tabular_model(0)

def test_with_bounding_box():
    """
    Test the option to evaluate a model respecting
    its bunding_box.
    """
    p = models.Polynomial2D(2) & models.Polynomial2D(2)
    m = models.Mapping((0, 1, 0, 1)) | p
    with NumpyRNGContext(1234567):
        m.parameters = np.random.rand(12)

    m.bounding_box = ((3, 9), (1, 8))
    x, y = np.mgrid[:10, :10]
    a, b = m(x, y)
    aw, bw = m(x, y, with_bounding_box=True)
    ind = (~np.isnan(aw)).nonzero()
    assert_allclose(a[ind], aw[ind])
    assert_allclose(b[ind], bw[ind])

    aw, bw = m(x, y, with_bounding_box=True, fill_value=1000)
    ind = (aw != 1000).nonzero()
    assert_allclose(a[ind], aw[ind])
    assert_allclose(b[ind], bw[ind])

    # test the order of bbox is not reversed for 1D models
    p = models.Polynomial1D(1, c0=12, c1=2.3)
    p.bounding_box = (0, 5)
    assert(p(1) == p(1, with_bounding_box=True))

    t3 = models.Shift(10) & models.Scale(2) & models.Shift(-1)
    t3.bounding_box = ((4.3, 6.9), (6, 15), (-1, 10))
    assert_allclose(t3([1, 1], [7, 7], [3, 5], with_bounding_box=True),
                    [[np.nan, 11], [np.nan, 14], [np.nan, 4]])

    trans3 = models.Shift(10) & models.Scale(2) & models.Shift(-1)
    trans3.bounding_box = ((4.3, 6.9), (6, 15), (-1, 10))
    assert_allclose(trans3(1, 7, 5, with_bounding_box=True), [11, 14, 4])

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_with_bounding_box():
    points = np.arange(5)
    values = np.array([1.5, 3.4, 6.7, 7, 32])
    t = models.Tabular1D(points, values)
    result = t(1, with_bounding_box=True)

    assert result == 3.4
    assert t.inverse(result, with_bounding_box=True) == 1.

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_bounding_box_with_units():
    points = np.arange(5)*u.pix
    lt = np.arange(5)*u.AA
    t = models.Tabular1D(points, lt)
    result = t(1*u.pix, with_bounding_box=True)

    assert result == 1.*u.AA
    assert t.inverse(result, with_bounding_box=True) == 1*u.pix

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular1d_inverse():
    """Test that the Tabular1D inverse is defined"""
    points = np.arange(5)
    values = np.array([1.5, 3.4, 6.7, 7, 32])
    t = models.Tabular1D(points, values)
    result = t.inverse((3.4, 6.7))
    assert_allclose(result, np.array((1., 2.)))

    # Check that it works for descending values in lookup_table
    t2 = models.Tabular1D(points, values[::-1])
    assert_allclose(t2.inverse.points[0], t2.lookup_table[::-1])

    result2 = t2.inverse((7, 6.7))
    assert_allclose(result2, np.array((1., 2.)))

    # Check that it errors on double-valued lookup_table
    points = np.arange(5)
    values = np.array([1.5, 3.4, 3.4, 32, 25])
    t = models.Tabular1D(points, values)
    with pytest.raises(NotImplementedError):
        t.inverse((3.4, 7.))

    # Check that Tabular2D.inverse raises an error
    table = np.arange(5*5).reshape(5, 5)
    points = np.arange(0, 5)
    points = (points, points)
    t3 = models.Tabular2D(points=points, lookup_table=table)
    with pytest.raises(NotImplementedError):
        t3.inverse((3, 3))

    # Check that it uses the same kwargs as the original model
    points = np.arange(5)
    values = np.array([1.5, 3.4, 6.7, 7, 32])
    t = models.Tabular1D(points, values)
    with pytest.raises(ValueError):
        t.inverse(100)
    t = models.Tabular1D(points, values, bounds_error=False, fill_value=None)
    result = t.inverse(100)
    assert_allclose(t(result), 100)

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_grid_shape_mismatch_error():
    points = np.arange(5)
    lt = np.mgrid[0:5, 0:5][0]
    with pytest.raises(ValueError) as err:
        models.Tabular2D(points, lt)
    assert str(err.value) == "Expected grid points in 2 directions, got 5."

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_repr():
    points = np.arange(5)
    lt = np.arange(5)
    t = models.Tabular1D(points, lt)
    assert repr(t) == "<Tabular1D(points=(array([0, 1, 2, 3, 4]),), lookup_table=[0 1 2 3 4])>"

    table = np.arange(5*5).reshape(5, 5)
    points = np.arange(0, 5)
    points = (points, points)
    t = models.Tabular2D(points=points, lookup_table=table)
    assert repr(t) == (
        "<Tabular2D(points=(array([0, 1, 2, 3, 4]), array([0, 1, 2, 3, 4])), "
        "lookup_table=[[ 0  1  2  3  4]\n"
        " [ 5  6  7  8  9]\n"
        " [10 11 12 13 14]\n"
        " [15 16 17 18 19]\n"
        " [20 21 22 23 24]])>"
    )

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_str():
    points = np.arange(5)
    lt = np.arange(5)
    t = models.Tabular1D(points, lt)
    assert str(t) == (
        "Model: Tabular1D\n"
        "N_inputs: 1\n"
        "N_outputs: 1\n"
        "Parameters: \n"
        "  points: (array([0, 1, 2, 3, 4]),)\n"
        "  lookup_table: [0 1 2 3 4]\n"
        "  method: linear\n"
        "  fill_value: nan\n"
        "  bounds_error: True"
    )

    table = np.arange(5*5).reshape(5, 5)
    points = np.arange(0, 5)
    points = (points, points)
    t = models.Tabular2D(points=points, lookup_table=table)
    assert str(t) == (
        "Model: Tabular2D\n"
        "N_inputs: 2\n"
        "N_outputs: 1\n"
        "Parameters: \n"
        "  points: (array([0, 1, 2, 3, 4]), array([0, 1, 2, 3, 4]))\n"
        "  lookup_table: [[ 0  1  2  3  4]\n"
        " [ 5  6  7  8  9]\n"
        " [10 11 12 13 14]\n"
        " [15 16 17 18 19]\n"
        " [20 21 22 23 24]]\n"
        "  method: linear\n"
        "  fill_value: nan\n"
        "  bounds_error: True"
    )

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_evaluate():
    points = np.arange(5)
    lt = np.arange(5)[::-1]
    t = models.Tabular1D(points, lt)
    assert (t.evaluate([1, 2, 3]) == [3, 2, 1]).all()
    assert (t.evaluate(np.array([1, 2, 3]) * u.m) == [3, 2, 1]).all()

    t.n_outputs = 2
    value = [np.array([3, 2, 1]), np.array([1, 2, 3])]
    with mk.patch.object(tabular_models, 'interpn', autospec=True, return_value=value) as mkInterpn:
        outputs = t.evaluate([1, 2, 3])
        for index, output in enumerate(outputs):
            assert np.all(value[index] == output)
        assert mkInterpn.call_count == 1

@pytest.mark.skipif("not HAS_SCIPY")
def test_tabular_module_name():
    """
    The module name must be set manually because
    these classes are created dynamically.
    """
    for model in [models.Tabular1D, models.Tabular2D]:
        assert model.__module__ == "astropy.modeling.tabular"

class classmodel(FittableModel):
    f = Parameter(default=1)
    x = Parameter(default=0)
    y = Parameter(default=2)

    def __init__(self, f=f.default, x=x.default, y=y.default):
        super().__init__(f, x, y)

    def evaluate(self):
        pass

class subclassmodel(classmodel):
    f = Parameter(default=3, fixed=True)
    x = Parameter(default=10)
    y = Parameter(default=12)
    h = Parameter(default=5)

    def __init__(self, f=f.default, x=x.default, y=y.default, h=h.default):
        super().__init__(f, x, y)

    def evaluate(self):
        pass

def test_parameter_inheritance():
    b = subclassmodel()
    assert b.param_names == ('f', 'x', 'y', 'h')
    assert b.h == 5
    assert b.f == 3
    assert b.f.fixed == True  # noqa: E712

def test_parameter_description():

    model = models.Gaussian1D(1.5, 2.5, 3.5)
    assert model.amplitude._description == "Amplitude (peak value) of the Gaussian"
    assert model.mean._description == "Position of peak (Gaussian)"

    model = models.Voigt1D(x_0=5, amplitude_L=10, fwhm_L=0.5, fwhm_G=0.9)
    assert model.amplitude_L._description == "The Lorentzian amplitude"
    assert model.fwhm_L._description == "The Lorentzian full width at half maximum"
    assert model.fwhm_G._description == "The Gaussian full width at half maximum"

def test_SmoothlyBrokenPowerLaw1D_validators():
    with pytest.raises(InputParameterError) as err:
        SmoothlyBrokenPowerLaw1D(amplitude=-1)
    assert str(err.value) == "amplitude parameter must be > 0"

    with pytest.raises(InputParameterError) as err:
        SmoothlyBrokenPowerLaw1D(delta=0)
    assert str(err.value) == "delta parameter must be >= 0.001"

@pytest.mark.skipif('not HAS_SCIPY')
@pytest.mark.filterwarnings(r'ignore:.*:RuntimeWarning')
@pytest.mark.filterwarnings(r'ignore:The fit may be unsuccessful.*')
def test_SmoothlyBrokenPowerLaw1D_fit_deriv():
    x_lim = [0.01, 100]

    x = np.logspace(x_lim[0], x_lim[1], 100)

    parameters = {'parameters': [1, 10, -2, 2, 0.5],
                  'constraints': {'fixed': {'x_break': True, 'delta': True}}}
    model_with_deriv = create_model(SmoothlyBrokenPowerLaw1D, parameters,
                                    use_constraints=False)
    model_no_deriv = create_model(SmoothlyBrokenPowerLaw1D, parameters,
                                  use_constraints=False)

    # NOTE: PR 10644 replaced deprecated usage of RandomState but could not
    #       find a new seed that did not cause test failure, resorted to hardcoding.
    # add 10% noise to the amplitude
    rsn_rand_1234567890 = np.array([
        0.61879477, 0.59162363, 0.88868359, 0.89165480, 0.45756748,
        0.77818808, 0.26706377, 0.99610621, 0.54009489, 0.53752161,
        0.40099938, 0.70540579, 0.40518559, 0.94999075, 0.03075388,
        0.13602495, 0.08297726, 0.42352224, 0.23449723, 0.74743526,
        0.65177865, 0.68998682, 0.16413419, 0.87642114, 0.44733314,
        0.57871104, 0.52377835, 0.62689056, 0.34869427, 0.26209748,
        0.07498055, 0.17940570, 0.82999425, 0.98759822, 0.11326099,
        0.63846415, 0.73056694, 0.88321124, 0.52721004, 0.66487673,
        0.74209309, 0.94083846, 0.70123128, 0.29534353, 0.76134369,
        0.77593881, 0.36985514, 0.89519067, 0.33082813, 0.86108824,
        0.76897859, 0.61343376, 0.43870907, 0.91913538, 0.76958966,
        0.51063556, 0.04443249, 0.57463611, 0.31382006, 0.41221713,
        0.21531811, 0.03237521, 0.04166386, 0.73109303, 0.74556052,
        0.64716325, 0.77575353, 0.64599254, 0.16885816, 0.48485480,
        0.53844248, 0.99690349, 0.23657074, 0.04119088, 0.46501519,
        0.35739006, 0.23002665, 0.53420791, 0.71639475, 0.81857486,
        0.73994342, 0.07948837, 0.75688276, 0.13240193, 0.48465576,
        0.20624753, 0.02298276, 0.54257873, 0.68123230, 0.35887468,
        0.36296147, 0.67368397, 0.29505730, 0.66558885, 0.93652252,
        0.36755130, 0.91787687, 0.75922703, 0.48668067, 0.45967890])
    n = 0.1 * parameters['parameters'][0] * (rsn_rand_1234567890 - 0.5)

    data = model_with_deriv(x) + n
    fitter_with_deriv = fitting.LevMarLSQFitter()
    new_model_with_deriv = fitter_with_deriv(model_with_deriv, x, data)
    fitter_no_deriv = fitting.LevMarLSQFitter()
    new_model_no_deriv = fitter_no_deriv(model_no_deriv, x, data,
                                         estimate_jacobian=True)
    assert_allclose(new_model_with_deriv.parameters,
                    new_model_no_deriv.parameters, atol=0.5)

class _ExtendedModelMeta(_ModelMeta):
    @classmethod
    def __prepare__(mcls, name, bases, **kwds):
        # this shows the parent class machinery still applies
        namespace = super().__prepare__(name, bases, **kwds)
        # the custom bit
        namespace.update(kwds)
        return namespace

    model = models.Gaussian1D(1.5, 2.5, 3.5)
    assert model.amplitude._description == "Amplitude (peak value) of the Gaussian"
    assert model.mean._description == "Position of peak (Gaussian)"

def test_metaclass_kwargs():
    """Test can pass kwargs to Models"""
    class ClassModel(FittableModel, flag="flag"):
        def evaluate(self):
            pass

    # Nothing further to test, just making the class is good enough.

def test_submetaclass_kwargs():
    """Test can pass kwargs to Model subclasses."""
    class ClassModel(FittableModel, metaclass=_ExtendedModelMeta, flag="flag"):
        def evaluate(self):
            pass

    assert ClassModel.flag == "flag"

class ModelDefault(Model):
    slope = Parameter()
    intercept = Parameter()
    _separable = False

    @staticmethod
    def evaluate(x, slope, intercept):
        return slope * x + intercept

class ModelCustom(ModelDefault):
    def _calculate_separability_matrix(self):
        return np.array([[0, ]])

def test_custom_separability_matrix():
    original = separability_matrix(ModelDefault(slope=1, intercept=2))
    assert original.all()

    custom = separability_matrix(ModelCustom(slope=1, intercept=2))
    assert not custom.any()

max_stars_repo_name: zabop/astropy

max_issues_repo_name: nabobalis/astropy

max_forks_repo_name: nabobalis/astropy

reason: import numpy,from numpy,import astropy,from astropy


text:
import numpy as np
import numpy.linalg as linalg
from simulationTools import *

nCores = 2
n = 8
hz = [.375, -.375, .375, -.375, .375, -.375, .375, -.375]
H = power_law_ising_hamiltonian(n, hz, 1.05, 1, 5, 6)

vAxis = [0,0]
V1 = construct_local_projector(n, 0, 0, vAxis)
V2 = construct_local_projector(n, 0, 1, vAxis)
xValues = np.arange(-30,30,.1)
pg1 = weak_measurement_couplings(10, .1, .16, xValues)
pg2 = weak_measurement_couplings(10, .1, .16, xValues)
wAxisList = [[0,0], ]*n
Wlist = construct_complete_projector_list(n, wAxisList)
tList = np.linspace(0,10,101)

t_unscramble = 5
Lambda, Q = linalg.eigh(H.todense())
rho_unscrambled = Wlist[0]
U_scramble = get_U(Lambda, Q, t_unscramble)
rho = U_scramble.getH() @ rho_unscrambled @ U_scramble

LHSneu, LHSmm, LHSneuTerms, LHSmmTerms, rhoF, rhoR, f, aOrd, aTerm, qProbs = simulation_strong(V1, V2, pg1, pg2, Wlist, H, rho, tList, nWorks=nCores)

np.save('./strong_data/tList.npy', tList)
np.save('./strong_data/LHSneu.npy', LHSneu)
np.save('./strong_data/LHSmm.npy', LHSmm)
np.save('./strong_data/LHSneuTerms.npy', LHSneuTerms)
np.save('./strong_data/LHSmmTerms.npy', LHSmmTerms)
np.save('./strong_data/rhoF.npy', rhoF)
np.save('./strong_data/rhoR.npy', rhoR)
np.save('./strong_data/f.npy', f)
np.save('./strong_data/aOrd.npy', aOrd)
np.save('./strong_data/aTerm.npy', aTerm)
np.save('./strong_data/qProbs.npy', qProbs)

max_stars_repo_name: abartolo-tb/Entropic-Unc-and-Quasiprobs

max_issues_repo_name: abartolo-tb/Entropic-Unc-and-Quasiprobs

max_forks_repo_name: abartolo-tb/Entropic-Unc-and-Quasiprobs


text:
from random import choice
import pytest
import numpy as np
import pandas as pd
from ..benford import utils as ut
from ..benford.constants import CONFS, REV_DIGS
from ..benford.expected import _get_expected_digits_
from ..benford.stats import _two_dist_ks_

@pytest.fixture
def gen_N():
    return np.random.randint(0, 25000)

@pytest.fixture
def gen_decimals():
    return np.random.randint(0, 8)

@pytest.fixture
def gen_N_lower(gen_N):
    return np.random.randint(0, gen_N)

@pytest.fixture
def gen_array(gen_N):
    num = gen_N
    return np.abs(np.random.rand(num) * np.random.randn(num) *
                  np.random.randint(1, num, num))

@pytest.fixture
def choose_digs_rand():
    return choice([1, 2, 3, 22, -2])

@pytest.fixture
def choose_test():
    return choice(["F1D","F2D","F3D","SD","L2D"])

@pytest.fixture
def choose_confidence():
    return choice(list(CONFS.keys())[1:])

@pytest.fixture
def gen_series(gen_array):
    return pd.Series(gen_array)

@pytest.fixture
def gen_data_frame(gen_array):
    return pd.DataFrame({'seq': gen_array, 'col2': gen_array})

@pytest.fixture
def gen_int_df(gen_data_frame):
    return gen_data_frame.astype(int)

small_arrays_type = [
    (np.array([1, 2, 3, 4, 5.0, 6.3, .17]), float),
    (np.array([1, 2, 3, 4, 5, 6, 7]), int),
    (np.array(['1', '2', '3', '4', '5', '6', '7']), float),
    (pd.Series([1, 2, 3, 4, 5.0, 6.3, .17]), float),
    (pd.Series([1, 2, 3, 4, 5, 6, 7]), int),
    (pd.Series(['1', '2', '3', '4', '5', '6', '7']), float)
]

@pytest.fixture(params=small_arrays_type)
def get_small_arrays(request):
    return request.param

@pytest.fixture
def small_str_foo_array():
    return np.array(['foo', 'baar', 'baz', 'hixks'])

@pytest.fixture
def small_str_foo_series():
    return pd.Series(['foo', 'baar', 'baz', 'hixks'])

@pytest.fixture
def gen_get_digs_df(gen_series, gen_decimals):
    return ut.get_digs(gen_series, decimals=gen_decimals)

@pytest.fixture
def gen_proportions_F1D(gen_get_digs_df):
    return ut.get_found_proportions(gen_get_digs_df.F1D)

@pytest.fixture
def gen_proportions_F2D(gen_get_digs_df):
    return ut.get_found_proportions(gen_get_digs_df.F2D)

@pytest.fixture
def gen_proportions_F3D(gen_get_digs_df):
    return ut.get_found_proportions(gen_get_digs_df.F3D)

@pytest.fixture
def gen_proportions_SD(gen_get_digs_df):
    return ut.get_found_proportions(gen_get_digs_df.SD)

@pytest.fixture
def gen_proportions_L2D(gen_get_digs_df):
    return ut.get_found_proportions(gen_get_digs_df.L2D)

@pytest.fixture
def gen_proportions_random_test(choose_test, gen_get_digs_df):
    dig_str = choose_test
    return ut.get_found_proportions(gen_get_digs_df[dig_str]), REV_DIGS[dig_str]

@pytest.fixture
def gen_join_expect_found_diff_random_test(gen_proportions_random_test):
    rand_test, rand_digs = gen_proportions_random_test
    return ut.join_expect_found_diff(rand_test, rand_digs)

@pytest.fixture
def gen_join_expect_found_diff_F1D(gen_proportions_F1D):
    return ut.join_expect_found_diff(gen_proportions_F1D, 1)

@pytest.fixture
def gen_join_expect_found_diff_F2D(gen_proportions_F2D):
    return ut.join_expect_found_diff(gen_proportions_F2D, 2)

@pytest.fixture
def gen_join_expect_found_diff_F3D(gen_proportions_F3D):
    return ut.join_expect_found_diff(gen_proportions_F3D, 3)

@pytest.fixture
def gen_join_expect_found_diff_SD(gen_proportions_SD):
    return ut.join_expect_found_diff(gen_proportions_SD, 22)

@pytest.fixture
def gen_join_expect_found_diff_L2D(gen_proportions_L2D):
    return ut.join_expect_found_diff(gen_proportions_L2D, -2)

@pytest.fixture
def gen_linspaced_zero_one(cuts:int=1000):
    return np.linspace(0, 1, cuts)

@pytest.fixture
def gen_mantissas_ks_dists(gen_array):
    dist2 = ut.get_mantissas(gen_array)
    dist1 = np.linspace(0, 1, len(dist2), endpoint=False)
    return dist1, dist2

def gen_mantissa_distribution():
    num =  np.random.randint(1500, 5000)
    a = np.random.rand(num)
    b = np.random.randint(1, 999, num)
    c = np.random.randn(num)
    abc = np.abs(a * b * c)
    return ut.get_mantissas(abc)

mant_ks_dists_types = [
    (gen_mantissa_distribution(), np.random.choice([True, False]), np.float_)
        for i in range(10)
]

@pytest.fixture(params=mant_ks_dists_types)
def get_mant_ks_types(request):
    dist2, cummulative, ks_type = request.param
    return np.linspace(0, 1, len(dist2), endpoint=False), \
        dist2, cummulative, ks_type

mant_dists = [
    (gen_mantissa_distribution(), np.random.choice([True, False]), 0)
        for i in range(10)
]

@pytest.fixture(params=mant_dists)
def get_mant_ks_s(request):
    dist2, cummulative, zero = request.param
    return np.linspace(0, 1, len(dist2), endpoint=False), \
        dist2, cummulative, zero

def gen_mantissa_distribution_len():
    mants = gen_mantissa_distribution()
    return mants, len(mants)

mant_ks_confidences = [
    (*gen_mantissa_distribution_len(), conf) for conf in CONFS.keys()
]

@pytest.fixture(params=mant_ks_confidences)
def get_mant_ks_confs_limit_N(request):
    mants, mants_lengths, confidence = request.param
    cap_sample = mants_lengths - np.random.randint(500, 1400)
    sample_size = np.random.choice([mants_lengths, cap_sample])
    return mants, confidence, sample_size

@pytest.fixture
def gen_random_digs_and_proportions(gen_linspaced_zero_one, choose_digs_rand):
    exp = _get_expected_digits_(choose_digs_rand).Expected.values
    rand_prop = np.random.choice(gen_linspaced_zero_one, len(exp))
    return exp, rand_prop / rand_prop.sum()

max_issues_repo_name: milcent/Benford

max_forks_repo_name: milcent/Benford