text:
import streamlit as st
import pandas as pd
import numpy as np
import plotly.express as px
import plotly.graph_objects as go

class st_functions:

    def __init__(self, DataFrame):
        self.df = DataFrame

    def clean_data(self):
        # Column to count how many types of class there is
        self.df['count'] = pd.Series(np.ones(self.df.shape[0])).astype(int)

        # Turning floats into integers
        self.df['year'] = self.df['year'].fillna(0).apply(lambda x: int(x))
        # Turning mass into kg
        self.df['mass (g)'] = self.df['mass (g)'].apply(lambda x: round(x/1000, 3))
        # Mass == 0 turn to NaN
        self.df['mass (g)'] = self.df['mass (g)'].replace(0, np.nan)

        # To numeric
        self.df['reclat'].replace(0, np.nan, inplace=True)
        self.df['reclong'].replace(0, np.nan, inplace=True)
        # df[['reclong', 'reclat']].dropna(inplace=True)
        # pd.to_numeric(df['reclong'])
        # pd.to_numeric(df['reclat'])
        # Rename columns
        self.df.rename(
            columns={'mass (g)': 'mass (kg)', 'reclat': 'lat', 'reclong': 'lon'},
            inplace=True
        )

    def intro(self):
        st.title('EDA METEORITE LANDINGS')
        st.text('By Marcos Alvarez Ots')

        st.write("""
        I'm a current Data Science student at [The Bridge](https://www.thebridge.tech/), challenged to build an EDA (Exploratory Data Analysis).
        Therefore, I have decided to analyse data about meteorite landings on Earth.

        Database provided and created by the [NASA](https://www.nasa.gov/) with license: [CC0 1.0 Universal Public Domain Dedication](https://creativecommons.org/publicdomain/zero/1.0/).

        This EDA has no correlation at all, which means it's meant to be explored by anyone whose intrigued by the matter and want to build own conclusions.
        """)

        # Image
        st.header('Feeding curiosity')
        st.write("Here there's a little example of a meteorite. Yet, not that important until we dive into the meteorites knowledge but, things are going to start make sense once we know more.")
        st.image('img/namibia-hoba.png')
        st.caption('Namibia, Africa. The Hoba meteorite is the largest known intact meteorite (as a single piece). It is also the most massive naturally occurring piece of iron with an estimated mass at more than 60 tonnes, measuring 2.7×2.7×0.9 m. The name "Hoba" comes from Khoekhoegowab word meaning "gift".')

        # Video
        st.header('Short introduction into the matter')
        st.write("""In order to provide a great introduction to meteorites it's important to understand where they come from. This short video gives a great and funny understanding of it.
        A meteorite is a solid object found on Earth's surface after it enters earth's atmosphere from interplanetary space. Meteorites are either observed "falls", or were found some time after they fell, i.e., "finds".""")
        st.video('vid/space_rocks.mp4', start_time=0)
        st.caption('Royal Museums Greenwich, (2016, Jan 27). Space Rocks | What\'s the difference between an asteroid, a meteoroid, a meteorite, and a comet?. Royal Museums Greenwich in UK. [Video]. Youtube. https://www.youtube.com/watch?v=rVdniqc_G_Q')

        # Classification
        st.subheader('Meteorite classification')
        st.write("""Once knowing the procedence here's a classification.
        Meteorites are classified into different types based on their mineral composition, textures, the presence or absence of chondrules, and other criteria. Chondritic meteorites, or \"chondrites\" are from parent bodies (asteroids) that never formed metal-rich cores. They contain many rock-forming elements in the same ratios as observed in the Sun's photosphere. \"Achondrites\" include irons (many from ancient planetary cores), rocks from the Moon, Mars, and Vesta, and meteorites that are not chondrites.""")
        st.image('img/meteoritefamilytree.gif')
        st.caption('Image obtained from [Nature History Museum](https://www.nhm.ac.uk/).')
        st.write('For further information about classification, visit [The Catalogue of Meteorites](https://www.nhm.ac.uk/our-science/data/metcat/search/bgmettypes.dsml) of Nature History Museum.')

    def data(self):
        # VIOLINPLOT
        st.subheader('Understanding data collection')
        st.write('''
        A few entries in the original databse contain date information that was incorrectly parsed into the NASA database.
        As a spot check: any date that is before 860 CE or after 2016 are incorrect; these should actually be BCE years.

        So in order to focus in relevant information I have proceeded to remove them for this plot.
        ''')
        df_0 = self.df[(self.df['year']!=0) & (self.df['year']<2100)]
        fig = px.violin(df_0, x='year', y='fall', color='fall')

        st.plotly_chart(fig, use_container_width=True)
        st.write("Most of the data belong to the last decades due to NASA's involvement to find meteorites, and for the same reason most of the found meteorites are from the USA.")

        # MAPPLOT
        df_1 = self.df.dropna(how='any')
        st.write("Here is where it's confirmed where they landed, showing it's mass by bubble size and showing the class pointing the mouse on each bubble.")
        fig = px.scatter_geo(df_1, color="fall",
            lat='lat',
            lon='lon',
            hover_name='recclass',
            size='mass (kg)',
            projection="natural earth")

        st.plotly_chart(fig, use_container_width=True)

        # BARPLOT
        st.subheader('Classifying meteorites by it\'s mass')
        st.write("This visualization is meant to bring you up to a conclusion build by yourself.")

        switch = st.radio(
            "Choose and option to display in y axis:",
            ('Mean', 'Standard deviation', 'Max'),
            )
        if switch == 'Mean':
            c = 'mean'
        elif switch == 'Standard deviation':
            c = 'std'
        else:
            c = 'max'

        sort = st.radio(
            "Sort by:",
            ('Count', 'Mean', 'Standard deviation', 'Max'),
            )
        if sort == 'Count':
            s = 'count'
        elif sort == 'Mean':
            s = 'mean'
        elif sort == 'Standard deviation':
            s = 'std'
        else:
            s = 'max'

        # recclass slider
        start, end = st.select_slider(
        'Select a range of recclass length:',
        options=list(np.arange(0, 388, 10)),
        value=(0, 100))

        df_2 = self.df.dropna().groupby('recclass')[['mass (kg)']].describe()
        df_2 = df_2['mass (kg)'].sort_values(by=s, ascending=False)
        df_2 = df_2[start:end]

        fig = px.bar(
            df_2,
            x=df_2.index,
            y=c,
            hover_data=['mean', 'std'],
            color='count',
            color_continuous_scale='Blugrn',
            title='Amount rank of each meteorite property classified by a given input',
            height=600,
            text='count'
        )
        fig.update_traces(textfont_size=12, textangle=0, textposition="outside", cliponaxis=False)
        st.plotly_chart(fig, use_container_width=True)
        st.write("""After the effort finding some insights we can conclude that meteorites are difficult to classify by it's mass. There is many groups of them and within them, more variations, doing harder this seeking. So here the effort is not given in classifying them but identifying them.
        So here an example. The most common meteorite are H, L and LL groups which belong to the "Ordinary" class. It's class is giving a valuous hint. They are called ordinary because they conform the 87% of the meteorites landed on Earth.
        Like this we can still find more interesting things and I encourage you to find them.""")


text:
from logging import getLogger
from json import dump
from os import environ
from numpy import linspace, r_, tile
from pathlib import Path
from shutil import move, copyfile, copyfileobj
from subprocess import run, DEVNULL
from tempfile import mkstemp, gettempdir
import gzip

from nibabel import save as nisave
from nibabel import load as niload
from bidso.utils import replace_extension

from .io.parrec import convert_parrec_nibabel
from .utils import rename_task, make_bids_name, find_one_file, make_taskdescription, set_notnone

lg = getLogger(__name__)

# DIRECTION might depend on the way that the data is stored in Nifti file.
# Hopefully nibabel is consistent in how it converts the data but we need to check (also the sign of qform)
DIRECTION = {
    'RL': 'i',
    'LR': 'i-',
    'PA': 'j',
    'AP': 'j-',
    'IS': 'k',
    'SI': 'k-',
    }

def convert_mri(run, rec, dest_path, name, deface=True):
    """Return base name for this run"""
    output_nii = dest_path / f'{make_bids_name(name)}_{rec.modality}.nii.gz'

    file = find_one_file(rec, ('parrec', ))
    if file is not None:
        input_nii, PAR = convert_parrec_nibabel(file.path, MagneticFieldStrength=run.session.MagneticFieldStrength)
        move(input_nii, output_nii)

    else:
        file = find_one_file(rec, ('nifti', ))
        PAR = None
        if file is None:
            return None

        else:
            input_nii = file.path

            if input_nii.name.endswith('.nii.gz'):
                copyfile(file.path, output_nii)
            elif input_nii.name.endswith('.nii'):
                gz(file.path, output_nii)
            else:
                lg.warning(f'Unknown extension for nifti for {input_nii}')
                return None

    if run.task_name == 'MP2RAGE':
        lg.info('Keeping only the first volume for MP2RAGE')
        select(output_nii, 'first')

    nii_shape = _fix_tr(output_nii, rec)

    if PAR is not None and 'phase' in PAR['image_types']:
        phase_file = dest_path / f'{make_bids_name(name)}_phase.nii.gz'
        lg.info(f'Splitting phase info to {phase_file.name}')
        phase_nii = select(output_nii, 'split')
        phase_nii.to_filename(phase_file)

    if deface and rec.modality in ('T1w', 'T2w', 'T2star', 'PD', 'FLAIR'):
        run_deface(output_nii)

    sidecar = _convert_sidecar(run, rec, PAR, nii_shape)
    sidecar_file = replace_extension(output_nii, '.json')

    with sidecar_file.open('w') as f:
        dump(sidecar, f, indent=2)

    return output_nii

def select(nii, slicing):
    """If slicing is "split", it returns the nifti of the second half"""
    img = niload(nii)
    half = int(img.shape[3] / 2)
    secondhalf = None

    if slicing == 'first':
        img = img.slicer[:, :, :, 0]

    elif slicing == 'firsthalf':
        img = img.slicer[:, :, :, :half]

    elif slicing == 'secondhalf':
        img = img.slicer[:, :, :, half:]

    elif slicing == 'split':
        secondhalf = img.slicer[:, :, :, half:]
        img = img.slicer[:, :, :, :half]

    img.to_filename(nii)

    return secondhalf

def _fix_tr(nii, rec):
    """
    Returns
    -------
    tuple
        shape of the nifti file
    """
    img = niload(str(nii))

    # this seems a bug in nibabel. It stores time in sec, not in msec
    img.header.set_xyzt_units('mm', 'sec')

    if rec.modality in ('bold', 'epi') and rec.RepetitionTime is not None:
        img.header['pixdim'][4] = rec.RepetitionTime

    nisave(img, str(nii))

    return img.shape

def _convert_sidecar(run, rec, hdr=None, shape=None):
    D = {
        'InstitutionName': 'University Medical Center Utrecht',
        'InstitutionAddress': 'Heidelberglaan 100, 3584 CX Utrecht, the Netherlands',
        }

    if run.session.MagneticFieldStrength is not None:
        D['MagneticFieldStrength'] = float(run.session.MagneticFieldStrength[:-1])

    if rec.PhaseEncodingDirection is not None:
        D['PhaseEncodingDirection'] = DIRECTION[rec.PhaseEncodingDirection]
    if rec.SliceEncodingDirection is not None:
        D['SliceEncodingDirection'] = DIRECTION[rec.SliceEncodingDirection]

    for field in 'EchoTime', 'EffectiveEchoSpacing', 'FlipAngle':
        set_notnone(D, hdr, field)

    for field in 'Sequence', 'MultibandAccelerationFactor':
        set_notnone(D, rec, field)

    if rec.modality in ('bold', 'epi'):
        set_notnone(D, hdr, 'RepetitionTime')  # first get RepetitionTime from Header,
        set_notnone(D, rec, 'RepetitionTime')  # then, use the one specified in rec
        D['TaskName'] = rename_task(run.task_name)
        D['TaskDescription'] = make_taskdescription(run)

        if hdr is not None:
            add_slicetiming(D, hdr, rec)
        if shape is not None and 'EffectiveEchoSpacing' in D and 'PhaseEncodingDirection' in D:
            NIFTI_INDEX = {'i': 0, 'j': 1, 'k': 2}
            ReconMatrixPE = shape[NIFTI_INDEX[D['PhaseEncodingDirection'][0]]]
            D['TotalReadoutTime'] = D['EffectiveEchoSpacing'] * (ReconMatrixPE - 1)

    return D

def add_slicetiming(D, hdr, rec):
    """

    TODO
    ----
    Is there a SliceTiming when the SliceOrder is 3D?
    """
    if rec.SliceOrder is None:
        lg.warning(f'Please specify SliceOrder for Recording(db, id={rec.id})')
        return

    elif rec.SliceOrder == '3D':
        return

    # get TR from SQL, but if it's not specified used PAR/REC
    TR = D.get('RepetitionTime', hdr['RepetitionTime'])

    n_slices = hdr['n_slices']

    multiband = D.get('MultibandAccelerationFactor', 1)
    n_slices = int(n_slices / multiband)

    SliceTiming = linspace(0, TR, n_slices + 1)[:-1]

    if rec.SliceOrder == 'Interleaved':
        SliceTiming = r_[SliceTiming[::2], SliceTiming[1::2]]

    SliceTiming = tile(SliceTiming, multiband)
    D['SliceTiming'] = SliceTiming.tolist()

def run_deface(nii):
    lg.info(f'Defacing {nii.name}, it might take a while')
    path_avg = Path(environ['FREESURFER_HOME']) / 'average'

    # generate a unique file name in the same folder
    if nii.name.endswith('.nii.gz'):
        suffix = '.nii.gz'
    elif nii.name.endswith('.nii'):
        suffix = '.nii'
    nii_tmp = Path(mkstemp(dir=nii.parent, suffix=suffix)[1])
    nii_tmp.unlink()

    run([
        'mri_deface',  # from freesurfer
        nii,
        path_avg / 'talairach_mixed_with_skull.gca',
        path_avg / 'face.gca',
        nii_tmp],
        stdout=DEVNULL, stderr=DEVNULL,
        cwd=gettempdir(),
        )

    nii_tmp.rename(nii)

def gz(input_file, output_file):
    with input_file.open('rb') as f_in:
        with gzip.open(output_file, 'wb') as f_out:
            copyfileobj(f_in, f_out)

ext: py

reason: from numpy


text:
"""CloneEstimator based on hdbscan."""
from .clones_estimator import ClonesEstimator
from ..utils import merge_dicts
import pandas as pd
import hdbscan
import logging
from scipy.spatial.distance import pdist, squareform

logger = logging.getLogger('chimaera.clustering.hdbscan')

default_parameters = {
    'min_cluster_size': 5
}

class Hdbscan(ClonesEstimator):
    def __init__(self):
        pass

    def estimate(self, mutations_data, **kwargs):
        try:
            precomputed = squareform(
                pdist(mutations_data, 'cityblock') / mutations_data.shape[1]
            )
        except BaseException:
            logger.exception(
                'Problem with mutations_data: {}'.format(mutations_data)
            )
            raise RuntimeError('Problem in computing distances for hdbscan.')
        kwargs.pop('metric', None)
        model = hdbscan.HDBSCAN(
            metric='precomputed',
            **merge_dicts(default_parameters, kwargs)
        )
        model.fit(precomputed)
        return pd.Series(model.labels_ + 1, index=mutations_data.index)

max_stars_repo_name: drugilsberg/chimaera

max_issues_repo_name: drugilsberg/chimaera

max_forks_repo_name: drugilsberg/chimaera

reason: from scipy


text:
from gym.benchmarks import registration
from gym.monitoring.tests import helpers
from gym.envs import rl2
from gym import monitoring
import gym

import numpy as np

def test_bandits():
    env = rl2.BernoulliBanditEnv(n_arms=5, n_episodes=10)
    env.seed(0)
    obs = env.reset()
    assert isinstance(obs, tuple)
    assert len(obs) == 4
    assert obs[-1] == 1
    # run it twice to make sure behavior is consistent
    for _ in range(2):
        path = random_rollout(env, n_steps=10)
        assert path["dones"][-1]
        assert not np.any(path["dones"][:-1])
        raw_obs, last_actions, last_rewards, last_dones = [list(x) for x in zip(*path["observations"])]
        # test that prev actions and actions taken agree
        np.testing.assert_array_equal(
            last_actions[1:],
            path["actions"][:-1]
        )
        # test that prev rewards and rewards taken agree
        np.testing.assert_array_equal(
            last_rewards[1:],
            path["rewards"][:-1]
        )
        # test that it always terminates
        assert np.all(last_dones)

def test_random_tabular_mdps():
    env = rl2.RandomTabularMDPEnv(n_states=10, n_actions=5, n_episodes=10, episode_length=10)
    env.seed(0)
    obs = env.reset()
    assert isinstance(obs, tuple)
    assert len(obs) == 4
    assert obs[-1] == 1
    # run it twice to make sure behavior is consistent
    for _ in range(2):
        path = random_rollout(env, n_steps=100)
        assert path["dones"][-1]
        assert not np.any(path["dones"][:-1])
        raw_obs, last_actions, last_rewards, last_dones = [list(x) for x in zip(*path["observations"])]
        # test that prev actions and actions taken agree
        np.testing.assert_array_equal(
            last_actions[1:],
            path["actions"][:-1]
        )
        # test that prev rewards and rewards taken agree
        np.testing.assert_array_equal(
            last_rewards[1:],
            path["rewards"][:-1]
        )
        last_dones = np.asarray(last_dones)
        # test that it terminates on the right times
        assert np.all(last_dones[::10])
        last_dones[::10] = 0
        assert not np.any(last_dones)
        # test that it always starts episodes on state 0
        assert np.all(np.asarray(raw_obs[::10]) == 0)

def test_benchmarks():
    for benchmark_id in ['BernoulliBandit-v0', 'RandomTabularMDP-v0']:

        benchmark = registration.benchmark_spec(benchmark_id)

        for env_id in benchmark.env_ids:

            with helpers.tempdir() as temp:
                env = gym.make(env_id)
                env.seed(0)
                env.monitor.start(temp, video_callable=False)

                env.monitor.configure(mode='evaluation')
                rollout(env)

                env.monitor.configure(mode='training')
                for i in range(2):
                    rollout(env)

                env.monitor.configure(mode='evaluation')
                rollout(env, good=True)

                env.monitor.close()
                results = monitoring.load_results(temp)
                evaluation_score = benchmark.score_evaluation(env_id, results['data_sources'],
                                                              results['initial_reset_timestamps'],
                                                              results['episode_lengths'],
                                                              results['episode_rewards'], results['episode_types'],
                                                              results['timestamps'])
                benchmark.score_benchmark({
                    env_id: evaluation_score['scores'],
                })

def rollout(env, good=False):
    env.reset()

    action = 0
    d = False
    while not d:
        if good:
            action = 1 - action
        o, r, d, i = env.step(action)

def random_rollout(env, n_steps):
    observations = []
    actions = []
    rewards = []
    dones = []
    obs = env.reset()
    for _ in range(n_steps):
        action = env.action_space.sample()
        next_obs, reward, done, _ = env.step(action)
        observations.append(obs)
        actions.append(action)
        rewards.append(reward)
        dones.append(done)
        obs = next_obs
    return dict(observations=observations, actions=actions, rewards=rewards, dones=dones)

test_benchmarks()

max_stars_repo_name: youngwoon/DnC-RL-Tensorflow

max_issues_repo_name: youngwoon/DnC-RL-Tensorflow

max_forks_repo_name: youngwoon/DnC-RL-Tensorflow

reason: import numpy


text:
# Plots the ssfr vs mass of the galaxies using UltraVISTA data.

'''
Usage: run muzzing_ssfr_mass_completeness dmass dssfr
e.g. run muzzing_ssfr_mass_completeness 0.2 0.5
dmass - size of log(mass) bin - default 0.1
dssfr - size of ssfr bin - default 0.5
'''

# muzzin_mass, muzzin_ssfr
from find_muzzin_mass_ssfr import unique_objs
import sys
import os
import string
import numpy as np
import pandas as pd
from astropy.io import ascii
from itertools import product
from tabulate import tabulate
import matplotlib.pyplot as plt

'''
unique_obs = dataframe of unique objects in our dataset, already filtered for stars. Has columns:
OBJID - Object ID
LMASS - Log Stellar Mass
sSFR - Specific Star Formation Rate
Measured - 1 if we have a measured redshift, 0 if not
'''

try:
    dmass = float(sys.argv[1])
    dssfr = float(sys.argv[2])
except:
    print('No input detected - using default dmass = 0.1, dssfr = 0.5')
    dmass, dssfr = (0.1, 0.5)

# Location to store figure output
figout = '/Users/galaxies-air/COSMOS/Paper_Images/'

def get_bins(dmass, dssfr):
    # Makes bins in mass and ssfr of size dmass and dssfr
    mass_bins = np.arange(9, 10+dmass, dmass)
    ssfr_bins = np.arange(-12, -8+dssfr, dssfr)
    return mass_bins, ssfr_bins

# WORKING HERE - NEED TO SEPARATE ALL AND UNIQUE WHEN BINNING, THEN DIVIDE FOR COMPLETENESS

def compute_completeness(mass_bins, ssfr_bins, unique_objs, measured_objs):
    # Computes the completeness in each box
    unique_binned = np.histogram2d(unique_objs['LMASS'], unique_objs['sSFR'], bins=[
        mass_bins, ssfr_bins])
    measured_binned = np.histogram2d(measured_objs['LMASS'], measured_objs['sSFR'], bins=[
        mass_bins, ssfr_bins])
    # Need to transpose it to get the axes to line up
    # The flip to account for ssfr being backwards. Compared by eye to the ssfr_mass plot

    measured_binned = np.flip(np.transpose(measured_binned[0]), axis=0)
    unique_binned = np.flip(np.transpose(unique_binned[0]), axis=0)
    completeness_arr = measured_binned/unique_binned
    return completeness_arr, unique_binned, measured_binned

measured_objs = unique_objs[unique_objs['Measured'] == 1]

mass_bins, ssfr_bins = get_bins(dmass, dssfr)
completeness_arr, unique_binned, measured_binned = compute_completeness(
    mass_bins, ssfr_bins, unique_objs, measured_objs)

# Fontsizes for plotting
axisfont = 24
ticksize = 18
ticks = 8
titlefont = 24
legendfont = 16
textfont = 16

# Figure setup
fig, ax = plt.subplots(figsize=(9, 7))

mark = '.'
xlab = 'log(Stellar Mass) (M$_\odot$)'
ylab = 'sSFR (yr$^{-1}$)'
xlim = (mass_bins[0], mass_bins[-1])
ylim = (ssfr_bins[0], ssfr_bins[-1])

# Extent sets where to do the imshow in data coords
extent = [mass_bins[0], mass_bins[-1], ssfr_bins[0], ssfr_bins[-1]]
# aspect stretches the pixes to make them boxes
image = ax.imshow(completeness_arr, extent=extent,
                  interpolation='nearest', aspect='auto')
cbar = plt.colorbar(image)
#cbar.ax.set_ylabel('Fraction Measured', rotation=270, fontsize=axisfont-8)
fig.text(0.955, 0.55, 'Fraction Measured', fontsize=axisfont,
         rotation=270, verticalalignment='center')
cbar.ax.tick_params(labelsize=ticksize)

# Plot all the objects
plt.scatter(unique_objs['LMASS'], unique_objs['sSFR'], color='black',
            marker=mark, label='Not measured')
# Overplot the objects that are good, which will cover the bad ones
plt.scatter(measured_objs['LMASS'], measured_objs['sSFR'], color='white',
            marker=mark, label='Measured')

'''
for (x, y) in product(np.arange(len(mass_bins)-1), np.arange(len(ssfr_bins)-1)):
    try:
        # Percentage rounded to nearest two digits, XX%
        plot_text = f'{int(round(completeness_arr[y,x],2)*100)}%'
        # Alternatively, fraction measured in raw counts:
        plot_text = f'{int(measured_binned[y,x])} / {int(unique_binned[y,x])}'
        # Aternatively, number of objects in that bin
        plot_text = f'{int(unique_binned[y,x])}'
        ax.text(mass_bins[x]+dmass/2, ssfr_bins[::-1][y+1]+dssfr/2, plot_text,
                horizontalalignment='center', verticalalignment='center', fontsize=12)
    except ValueError:
        # This error is thrown when the value is nan, so we will not write in that square
        pass
'''

# Set the axis labels
ax.set_xlabel(xlab, fontsize=axisfont)
ax.set_ylabel(ylab, fontsize=axisfont)

# Set the limits
ax.set_xlim(xlim)
ax.set_ylim(ylim)

# Set the tick size
ax.tick_params(labelsize=ticksize, size=ticks)
# ax.legend(fontsize=axisfont-6)

plt.tight_layout()

# Save the figure
fig.savefig(figout+'ssfr_mass_density.pdf')
plt.close('all')

reason: import numpy,from astropy


text:
import cv2
import numpy as np
import matplotlib.pyplot as plt

img=cv2.imread('gsw1.jpg',cv2.IMREAD_GRAYSCALE)
#IMREAD_COLOR
#IMREAD_UNCHANGED

##cv2.imshow('image',img)
##cv2.waitKey(0)
##cv2.destroyAllWindows()

##plt.imshow(img,cmap='gray',interpolation='bicubic')
##plt.plot([50,100],[80,100])
##plt.show()

cv2.imwrite('gsw1.jpg',img)#to write images

max_stars_repo_path: 1intro and loading images.py

max_issues_repo_path: 1intro and loading images.py

max_forks_repo_path: 1intro and loading images.py


text:
"""
BROS
Copyright 2022-present NAVER Corp.
Apache License v2.0
"""

import numpy as np
import torch
import torch.utils.data
from overrides import overrides

from lightning_modules.bros_module import BROSModule

class BROSSPADERELModule(BROSModule):
    def __init__(self, cfg):
        super().__init__(cfg)

        self.eval_kwargs = {"dummy_idx": self.cfg.train.max_seq_length}

    @overrides
    def training_step(self, batch, batch_idx, *args):
        _, loss = self.net(batch)

        log_dict_input = {"train_loss": loss}
        self.log_dict(log_dict_input, sync_dist=True)
        return loss

    @torch.no_grad()
    @overrides
    def validation_step(self, batch, batch_idx, *args):
        head_outputs, loss = self.net(batch)
        step_out = do_eval_step(batch, head_outputs, loss, self.eval_kwargs)
        return step_out

    @torch.no_grad()
    @overrides
    def validation_epoch_end(self, validation_step_outputs):
        scores = do_eval_epoch_end(validation_step_outputs)
        self.print(
            f"precision: {scores['precision']:.4f}, recall: {scores['recall']:.4f}, f1: {scores['f1']:.4f}"
        )

def do_eval_step(batch, head_outputs, loss, eval_kwargs):
    dummy_idx = eval_kwargs["dummy_idx"]

    el_outputs = head_outputs["el_outputs"]

    pr_el_labels = torch.argmax(el_outputs, -1)

    n_batch_gt_rel, n_batch_pr_rel, n_batch_correct_rel = eval_el_spade_batch(
        pr_el_labels,
        batch["el_labels"],
        batch["are_box_first_tokens"],
        dummy_idx,
    )

    step_out = {
        "loss": loss,
        "n_batch_gt_rel": n_batch_gt_rel,
        "n_batch_pr_rel": n_batch_pr_rel,
        "n_batch_correct_rel": n_batch_correct_rel,
    }

    return step_out

def eval_el_spade_batch(
    pr_el_labels,
    gt_el_labels,
    are_box_first_tokens,
    dummy_idx,
):
    n_batch_gt_rel, n_batch_pr_rel, n_batch_correct_rel = 0, 0, 0

    bsz = pr_el_labels.shape[0]
    for example_idx in range(bsz):
        n_gt_rel, n_pr_rel, n_correct_rel = eval_el_spade_example(
            pr_el_labels[example_idx],
            gt_el_labels[example_idx],
            are_box_first_tokens[example_idx],
            dummy_idx,
        )

        n_batch_gt_rel += n_gt_rel
        n_batch_pr_rel += n_pr_rel
        n_batch_correct_rel += n_correct_rel

    return n_batch_gt_rel, n_batch_pr_rel, n_batch_correct_rel

def eval_el_spade_example(pr_el_label, gt_el_label, box_first_token_mask, dummy_idx):
    gt_relations = parse_relations(gt_el_label, box_first_token_mask, dummy_idx)
    pr_relations = parse_relations(pr_el_label, box_first_token_mask, dummy_idx)

    n_gt_rel = len(gt_relations)
    n_pr_rel = len(pr_relations)
    n_correct_rel = len(gt_relations & pr_relations)

    return n_gt_rel, n_pr_rel, n_correct_rel

def parse_relations(el_label, box_first_token_mask, dummy_idx):
    valid_el_labels = el_label * box_first_token_mask
    valid_el_labels = valid_el_labels.cpu().numpy()
    el_label_np = el_label.cpu().numpy()

    valid_token_indices = np.where(
        ((valid_el_labels != dummy_idx) * (valid_el_labels != 0))
    )
    link_map_tuples = []
    for token_idx in valid_token_indices[0]:
        link_map_tuples.append((el_label_np[token_idx], token_idx))

    return set(link_map_tuples)

def do_eval_epoch_end(step_outputs):
    n_total_gt_rel, n_total_pred_rel, n_total_correct_rel = 0, 0, 0

    for step_out in step_outputs:
        n_total_gt_rel += step_out["n_batch_gt_rel"]
        n_total_pred_rel += step_out["n_batch_pr_rel"]
        n_total_correct_rel += step_out["n_batch_correct_rel"]

    precision = 0.0 if n_total_pred_rel == 0 else n_total_correct_rel / n_total_pred_rel
    recall = 0.0 if n_total_gt_rel == 0 else n_total_correct_rel / n_total_gt_rel
    f1 = (
        0.0
        if recall * precision == 0
        else 2.0 * recall * precision / (recall + precision)
    )

    scores = {
        "precision": precision,
        "recall": recall,
        "f1": f1,
    }

    return scores

max_stars_repo_name: clovaai/bros

max_issues_repo_name: clovaai/bros

max_forks_repo_name: clovaai/bros


text:
from tkinter import *
from PIL import ImageTk,Image
from tkinter import filedialog
import os
from datetime import datetime
from flask import Flask, render_template, request, url_for
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import cv2
import imageio
import scipy.ndimage as ndi
import argparse
import torch
import torch.nn as nn
import torch.nn.functional as F
from torchvision.models.vgg import vgg19
import pywt
import pywt.data
from skimage.morphology import extrema
from skimage.morphology import watershed as skwater

global my_label_1,my_label_2,ct_x_label,ct_y_label,mri_x_label,mri_y_label,ct, mri_registered, mri_registered_cv, mri_registered_label, mri_registered_image_label, ct_registered_label, ct_registered_image_label, fusion_button, fusion_img_display_label, fused_image_label

# Root window
root=Tk()
root.title('Multi-modal medical image fusion to detect brain tumors')
# root.geometry('900x500')
scroll_bar = Scrollbar(root)

# Helper functions
def procrustes(X, Y, scaling=True, reflection='best'):
    n,m = X.shape
    ny,my = Y.shape

    muX = X.mean(0)
    muY = Y.mean(0)

    X0 = X - muX
    Y0 = Y - muY

    ssX = (X0**2.).sum()
    ssY = (Y0**2.).sum()

    # centred Frobenius norm
    normX = np.sqrt(ssX)
    normY = np.sqrt(ssY)

    # scale to equal (unit) norm
    X0 /= normX
    Y0 /= normY

    if my < m:
        Y0 = np.concatenate((Y0, np.zeros(n, m-my)),0)

    # optimum rotation matrix of Y
    A = np.dot(X0.T, Y0)
    U,s,Vt = np.linalg.svd(A,full_matrices=False)
    V = Vt.T
    T = np.dot(V, U.T)

    if reflection != 'best':

        # does the current solution use a reflection?
        have_reflection = np.linalg.det(T) < 0

        # if that's not what was specified, force another reflection
        if reflection != have_reflection:
            V[:,-1] *= -1
            s[-1] *= -1
            T = np.dot(V, U.T)

    traceTA = s.sum()

    if scaling:

        # optimum scaling of Y
        b = traceTA * normX / normY

        # standarised distance between X and b*Y*T + c
        d = 1 - traceTA**2
        # transformed coords
        Z = normX*traceTA*np.dot(Y0, T) + muX

    else:
        b = 1
        d = 1 + ssY/ssX - 2 * traceTA * normY / normX
        Z = normY*np.dot(Y0, T) + muX

    # transformation matrix
    if my < m:
        T = T[:my,:]
    c = muX - b*np.dot(muY, T)
    #rot =1
    #scale=2
    #translate=3
    #transformation values
    tform = {'rotation':T, 'scale':b, 'translation':c}

    return d, Z, tform

# VGG19 CNN For Fusion
class VGG19(torch.nn.Module):
    def __init__(self, device='cpu'):
        super(VGG19, self).__init__()
        features = list(vgg19(pretrained=True).features)
        if device == "cuda":
            self.features = nn.ModuleList(features).cuda().eval()
        else:
            self.features = nn.ModuleList(features).eval()

    def forward(self, x):
        feature_maps = []
        for idx, layer in enumerate(self.features):
            x = layer(x)
            if idx == 3:
                feature_maps.append(x)
        return feature_maps

class Fusion:
    def __init__(self, input):
        """
        Class Fusion constructor

        Instance Variables:
            self.images: input images
            self.model: CNN model, default=vgg19
            self.device: either 'cuda' or 'cpu'
        """
        self.input_images = input
        self.device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
        self.model = VGG19(self.device)

    def fuse(self):
        """
        A top level method which fuse self.images
        """
        # Convert all images to YCbCr format
        self.normalized_images = [-1 for img in self.input_images]
        self.YCbCr_images = [-1 for img in self.input_images]
        for idx, img in enumerate(self.input_images):
            if not self._is_gray(img):
                self.YCbCr_images[idx] = self._RGB_to_YCbCr(img)
                self.normalized_images[idx] = self.YCbCr_images[idx][:, :, 0]
            else:
                self.normalized_images[idx] = img / 255.
        # Transfer all images to PyTorch tensors
        self._tranfer_to_tensor()
        # Perform fuse strategy
        fused_img = self._fuse()[:, :, 0]
        # Reconstruct fused image given rgb input images
        for idx, img in enumerate(self.input_images):
            if not self._is_gray(img):
                self.YCbCr_images[idx][:, :, 0] = fused_img
                fused_img = self._YCbCr_to_RGB(self.YCbCr_images[idx])
                fused_img = np.clip(fused_img, 0, 1)

        return (fused_img * 255).astype(np.uint8)
        # return fused_img

    def _fuse(self):
        """
        Perform fusion algorithm
        """
        with torch.no_grad():

            imgs_sum_maps = [-1 for tensor_img in self.images_to_tensors]
            for idx, tensor_img in enumerate(self.images_to_tensors):
                imgs_sum_maps[idx] = []
                feature_maps = self.model(tensor_img)
                for feature_map in feature_maps:
                    sum_map = torch.sum(feature_map, dim=1, keepdim=True)
                    imgs_sum_maps[idx].append(sum_map)

            max_fusion = None
            for sum_maps in zip(*imgs_sum_maps):
                features = torch.cat(sum_maps, dim=1)
                weights = self._softmax(F.interpolate(features,
                                        size=self.images_to_tensors[0].shape[2:]))
                weights = F.interpolate(weights,
                                        size=self.images_to_tensors[0].shape[2:])
                current_fusion = torch.zeros(self.images_to_tensors[0].shape)
                for idx, tensor_img in enumerate(self.images_to_tensors):
                    current_fusion += tensor_img * weights[:,idx]
                if max_fusion is None:
                    max_fusion = current_fusion
                else:
                    max_fusion = torch.max(max_fusion, current_fusion)

            output = np.squeeze(max_fusion.cpu().numpy())
            if output.ndim == 3:
                output = np.transpose(output, (1, 2, 0))
            return output

    def _RGB_to_YCbCr(self, img_RGB):
        """
        A private method which converts an RGB image to YCrCb format
        """
        img_RGB = img_RGB.astype(np.float32) / 255.
        return cv2.cvtColor(img_RGB, cv2.COLOR_RGB2YCrCb)

    def _YCbCr_to_RGB(self, img_YCbCr):
        """
        A private method which converts a YCrCb image to RGB format
        """
        img_YCbCr = img_YCbCr.astype(np.float32)
        return cv2.cvtColor(img_YCbCr, cv2.COLOR_YCrCb2RGB)

    def _is_gray(self, img):
        """
        A private method which returns True if image is gray, otherwise False
        """
        if len(img.shape) < 3:
            return True
        if img.shape[2] == 1:
            return True
        b, g, r = img[:,:,0], img[:,:,1], img[:,:,2]
        if (b == g).all() and (b == r).all():
            return True
        return False

    def _softmax(self, tensor):
        """
        A private method which compute softmax ouput of a given tensor
        """
        tensor = torch.exp(tensor)
        tensor = tensor / tensor.sum(dim=1, keepdim=True)
        return tensor

    def _tranfer_to_tensor(self):
        """
        A private method to transfer all input images to PyTorch tensors
        """
        self.images_to_tensors = []
        for image in self.normalized_images:
            np_input = image.astype(np.float32)
            if np_input.ndim == 2:
                np_input = np.repeat(np_input[None, None], 3, axis=1)
            else:
                np_input = np.transpose(np_input, (2, 0, 1))[None]
            if self.device == "cuda":
                self.images_to_tensors.append(torch.from_numpy(np_input).cuda())
            else:
                self.images_to_tensors.append(torch.from_numpy(np_input))

#Segmentation
def segment_image():
    segmentation_button.destroy()
    fused_image_label.destroy()
    fusion_img_display_label.destroy()

    img = cv2.imread("static/fusion.jpg")
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    ret, thresh = cv2.threshold(gray,0,255,cv2.THRESH_OTSU)

    ret, markers = cv2.connectedComponents(thresh)

    marker_area = [np.sum(markers==m) for m in range(np.max(markers)) if m!=0]
    largest_component = np.argmax(marker_area)+1 #Add 1 since we dropped zero above
    brain_mask = markers==largest_component

    brain_out = img.copy()
    brain_out[brain_mask==False] = (0,0,0)

    img = cv2.imread("static/fusion.jpg")
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    ret, thresh = cv2.threshold(gray,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)

    # noise removal
    kernel = np.ones((3,3),np.uint8)
    opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2)

    # sure background area
    sure_bg = cv2.dilate(opening,kernel,iterations=3)

    # Finding sure foreground area
    dist_transform = cv2.distanceTransform(opening,cv2.DIST_L2,5)
    ret, sure_fg = cv2.threshold(dist_transform,0.7*dist_transform.max(),255,0)

    # Finding unknown region
    sure_fg = np.uint8(sure_fg)
    unknown = cv2.subtract(sure_bg,sure_fg)

    # Marker labelling
    ret, markers = cv2.connectedComponents(sure_fg)

    # Add one to all labels so that sure background is not 0, but 1
    markers = markers+1

    # Now, mark the region of unknown with zero
    markers[unknown==255] = 0
    markers = cv2.watershed(img,markers)
    img[markers == -1] = [255,0,0]

    im1 = cv2.cvtColor(img,cv2.COLOR_HSV2RGB)
    cv2.imwrite("static/segmented.jpg", im1)

    seg_label = Label(root,text="Segmented Image")
    seg_label.grid(row=0,column=0)

    segmented_image_display = Image.fromarray(im1)
    segmented_image_1 = ImageTk.PhotoImage(image=segmented_image_display)
    segmented_image_display_label = Label(image=segmented_image_1)
    segmented_image_display_label.grid(row=1,column=0)

    segmented_image_display_label.image = segmented_image_1

def fuse_image():
    mri_registered_image_label.destroy()
    ct_registered_image_label.destroy()
    fusion_button.destroy()
    mri_registered_label.destroy()
    ct_registered_label.destroy()

    global fused_image_label, fusion_img_display_label, segmentation_button

    # Wavelet transform of image, and plot approximation and details
    coeffs1 = pywt.dwt2(mri_registered_cv, 'haar')
    LL_mri, (LH_mri, HL_mri, HH_mri) = coeffs1

    for i, a in enumerate([LL_mri, LH_mri, HL_mri, HH_mri]):
        path='static/mri_'+str(i)+'.jpg'
        cv2.imwrite(path,a)

    # Wavelet transform of image, and plot approximation and details
    coeffs2 = pywt.dwt2(ct_cv, 'haar')
    LL_ct, (LH_ct, HL_ct, HH_ct) = coeffs2

    for i, a in enumerate([LL_ct, LH_ct, HL_ct, HH_ct]):
        path='static/ct_'+str(i)+'.jpg'
        cv2.imwrite(path,a)

    #Performing fusion
    input_images=[]
    mri = cv2.imread('static/mri_0.jpg')
    mri = cv2.cvtColor(mri, cv2.COLOR_BGR2GRAY)

    ct = cv2.imread('static/ct_0.jpg')
    ct = cv2.cvtColor(ct, cv2.COLOR_BGR2GRAY)
    input_images.append(mri)
    input_images.append(ct)
    FU = Fusion(input_images)
    fusion_1 = FU.fuse()
    cv2.imwrite('static/fusion_1.jpg', fusion_1)

    input_images=[]
    mri = cv2.imread('static/mri_1.jpg')
    mri = cv2.cvtColor(mri, cv2.COLOR_BGR2GRAY)

    ct = cv2.imread('static/ct_1.jpg')
    ct = cv2.cvtColor(ct, cv2.COLOR_BGR2GRAY)
    input_images.append(mri)
    input_images.append(ct)
    FU = Fusion(input_images)
    fusion_2 = FU.fuse()
    cv2.imwrite('static/fusion_2.jpg', fusion_2)

    input_images=[]
    mri = cv2.imread('static/mri_2.jpg')
    mri = cv2.cvtColor(mri, cv2.COLOR_BGR2GRAY)

    ct = cv2.imread('static/ct_2.jpg')
    ct = cv2.cvtColor(ct, cv2.COLOR_BGR2GRAY)
    input_images.append(mri)
    input_images.append(ct)
    FU = Fusion(input_images)
    fusion_3 = FU.fuse()
    cv2.imwrite('static/fusion_3.jpg', fusion_3)

    input_images=[]
    mri = cv2.imread('static/mri_3.jpg')
    mri = cv2.cvtColor(mri, cv2.COLOR_BGR2GRAY)

    ct = cv2.imread('static/ct_3.jpg')
    ct = cv2.cvtColor(ct, cv2.COLOR_BGR2GRAY)
    input_images.append(mri)
    input_images.append(ct)
    FU = Fusion(input_images)
    fusion_4 = FU.fuse()
    cv2.imwrite('static/fusion_4.jpg', fusion_4)

    coeffs=(fusion_1,(fusion_2,fusion_3,fusion_4))
    fusion=pywt.idwt2(coeffs,'haar')

    cv2.imwrite('static/fusion.jpg',fusion)

    fused_image_label=Label(root,text="Fused Image")
    fused_image_label.grid(row=0,column=0)

    fusion_img_display = Image.fromarray(fusion)
    fusion_img_1 = ImageTk.PhotoImage(image=fusion_img_display)
    fusion_img_display_label = Label(image=fusion_img_1)
    fusion_img_display_label.grid(row=1,column=0)

    fusion_img_display_label.image = fusion_img_1

    segmentation_button = Button(root,text="Perform Segmentation",command=segment_image)
    segmentation_button.grid(row=2, column=0)

# Registration
def register():
    frame_file.destroy()
    canvas_mri.destroy()
    canvas_ct.destroy()
    registration_button.destroy()
    my_label_1.destroy()
    my_label_2.destroy()
    # mri_x_label.destroy()
    # mri_y_label.destroy()
    # ct_x_label.destroy()
    # ct_y_label.destroy()

    global mri_registered, mri_registered_cv, mri_registered_label, mri_registered_image_label, ct_registered_label, ct_registered_image_label, fusion_button

    X_pts = np.asarray(ct_points)
    Y_pts = np.asarray(mri_points)
    d,Z_pts,Tform = procrustes(X_pts,Y_pts)
    R = np.eye(3)
    R[0:2,0:2] = Tform['rotation']

    S = np.eye(3) * Tform['scale']
    S[2,2] = 1
    t = np.eye(3)
    t[0:2,2] = Tform['translation']
    M = np.dot(np.dot(R,S),t.T).T
    h=ct_cv.shape[0]
    w=ct_cv.shape[1]
    mri_registered_cv = cv2.warpAffine(mri_cv,M[0:2,:],(h,w))

    mri_registered_image = Image.fromarray(mri_registered_cv)
    mri_registered = ImageTk.PhotoImage(image=mri_registered_image)
    mri_registered_image_label = Label(image=mri_registered)
    mri_registered_image_label.grid(row=1,column=0)

    mri_registered_label=Label(root,text="MRI Registered Image")
    mri_registered_label.grid(row=0,column=0)
    ct_registered_label=Label(root,text="CT Registered Image")
    ct_registered_label.grid(row=0,column=1)

    # ct_registered_image = Image.fromarray(ct)
    # ct_registered_image_1 = ImageTk.PhotoImage(image=ct_registered_image)

    ct_registered_image_1 = ct_image
    ct_2=ct_registered_image_1
    ct_registered_image_label=Label(image=ct_registered_image_1)
    ct_registered_image_label.grid(row=1,column=1)

    ct_registered_image_label.image=ct_registered_image_1

    fusion_button = Button(root,text="Perform Fusion",command=fuse_image)
    fusion_button.grid(row=2, column=0)

# Upload Files frame
frame_file=LabelFrame(root, text="Select files:",pady=20)
frame_file.grid(row=0,column=0)

# Display Uploaded Images
canvas_mri = Canvas(root,width=512, height=512)
canvas_mri.grid(row=2,column=0)

canvas_ct = Canvas(root,width=512, height=512)
canvas_ct.grid(row=2,column=1)

# Points for registration
mri_points=[]
ct_points=[]

def openMRI():
    global my_image_1,mri_image,my_label_1, mri_cv

    path=filedialog.askopenfilename(title="Select MRI Image")

    mri_cv = cv2.imread(path)
    mri_cv = cv2.cvtColor(mri_cv, cv2.COLOR_BGR2GRAY)

    my_label_1=Label(root,text="MRI Image")
    my_label_1.grid(row=1,column=0)
    my_image_1=ImageTk.PhotoImage(Image.open(path))
    mri_image=my_image_1

    canvas_mri.create_image(0, 0, image=my_image_1, anchor="nw")
    canvas_mri.config(scrollregion=canvas_mri.bbox(ALL))

    def printcoordsMRI(event):
        global mri_x_label,mri_y_label

        mri_x_label=Label(root, text="MRI X:"+str(event.x))
        mri_x_label.grid(row=3,column=0)
        mri_y_label=Label(root, text="MRI Y:"+str(event.y))
        mri_y_label.grid(row=4,column=0)
        # print (event.x,event.y)
        mri_points.append([event.x,event.y])

        mri_x_label.after(3000, mri_x_label.destroy)
        mri_y_label.after(3000, mri_y_label.destroy)
    #mouseclick event
    canvas_mri.bind("<Button 1>",printcoordsMRI)

def openCT():
    global my_image_2,ct_image,registration_button,my_label_2,ct_cv

    path=filedialog.askopenfilename(title="Select CT Image")

    ct_cv = cv2.imread(path)
    ct_cv = cv2.cvtColor(ct_cv, cv2.COLOR_BGR2GRAY)

    my_label_2 = Label(root,text="CT Image")
    my_label_2.grid(row=1,column=1)
    my_image_2=ImageTk.PhotoImage(Image.open(path))
    ct_image=my_image_2

    canvas_ct.create_image(0, 0, image=my_image_2, anchor="nw")
    canvas_ct.config(scrollregion=canvas_ct.bbox(ALL))

    # submit registration points button
    registration_button = Button(root,text="Submit Points",command=register)
    registration_button.grid(row=5,column=0)

    def printcoordsCT(event):
        global ct_x_label,ct_y_label

        ct_x_label=Label(root, text="CT X:"+str(event.x))
        ct_x_label.grid(row=3,column=1)
        ct_y_label=Label(root, text="CT Y:"+str(event.y))
        ct_y_label.grid(row=4,column=1)
        # print (event.x,event.y)
        ct_points.append([event.x,event.y])

        ct_x_label.after(3000, ct_x_label.destroy)
        ct_y_label.after(3000, ct_y_label.destroy)
    #mouseclick event
    canvas_ct.bind("<Button 1>",printcoordsCT)

mri_button=Button(frame_file,text="Select MRI File",command=openMRI)
mri_button.grid(row=0,column=0,pady=10,padx=10)

ct_button=Button(frame_file,text="Select CT File",command=openCT)
ct_button.grid(row=0,column=1,pady=10,padx=10)

root.mainloop()

reason: import numpy,import scipy


text:
import numpy as np
import skimage.morphology
import jtmri.fit

def fit(t, data, min_threshold=20):
    selem = np.ones((2,2))

    mask = data < min_threshold
    for idx in np.ndindex(*mask.shape[2:]):
        slc = (slice(None,None), slice(None,None)) + idx
        mask[slc] = skimage.morphology.binary_dilation(mask[slc], selem=selem)
    mask = np.tile(mask.any(axis=-1)[...,np.newaxis], mask.shape[-1])
    return jtmri.fit.fit_r2star_fast(t, np.ma.array(data, mask=mask))

max_stars_repo_name: jthacker/jtmri

max_issues_repo_name: jthacker/jtmri

max_forks_repo_name: jthacker/jtmri


text:
#!/usr/bin/env python

import numpy as np
import os
import math
from collate import readstruc
from glob import glob

def gap_creator(labelend, inter_r = [], rho_deltas = [], temp_deltas = [],
    epsbig_deltas = [], eps_deltas = []):
    '''
    Function that modifies the structure of a disk model saved in
    the fort14.irr file. The modified fort14 is saved into a fort14.mod file.

    INPUTS:
    - labelend: labelend of the model, used to find the correct fort14.irr
    - inter_r: list of radii corresponding to the outer edges of the intervals in
    radius that want to be modified.
    - rho_deltas: list of deltas by which the density will be multiplied in the
    radius intervals provided.
    - temp_deltas: list of deltas for temperature.
    - epsbig_deltas: list of deltas for epsilon big.
    - eps_deltas: list of deltas for epsilon.

    OUTPUTS:
    A modified fort14.mod file with the new structure of the disk.

    EXAMPLE:
    - inter_r = [10.]; rho_deltas = [0.1]:
    Between 0 and 10 au the density will be multiplied by 0.1.
    - inter_r = [5.,20.]; rho_deltas =  [1.0,0.1]:
    Between 0 and 5 au the density will not be modified (multiplied by 1).
    Between 5 and 20 au the density will be multiplied by 0.1.
    '''
    if len(inter_r) == 0:
        print('GAP_CREATOR: No modification provided.')
        return
    if rho_deltas:
        if len(inter_r) != len(rho_deltas):
            raise IOError('GAP_CREATOR: Length of inter_r should be the same as rho_deltas')
    if temp_deltas:
        if len(inter_r) != len(temp_deltas):
            raise IOError('GAP_CREATOR: Length of inter_r should be the same as temp_deltas')
    if epsbig_deltas:
        if len(inter_r) != len(epsbig_deltas):
            raise IOError('GAP_CREATOR: Length of inter_r should be the same as epsbig_deltas')
    if eps_deltas:
        if len(inter_r) != len(eps_deltas):
            raise IOError('GAP_CREATOR: Length of inter_r should be the same as eps_deltas')

    files = glob('fort14*irr.'+labelend+'*')
    if len(files) == 0:
        raise IOError('GAP_CREATOR: fort14.irr file not found for '+labelend)
    elif len(files) > 1:
        raise IOError('GAP_CREATOR: More than one fort14.irr file found for '+labelend)

    fort14 = files[0]
    radii, z, p, t, rho, epsbig, eps = readstruc(fort14) # arrays as [nrad,nz]
    radii = np.meshgrid(radii,z[0,:])[0].T

    # We modify the density according to the given intervals and deltas
    r_0 = 0.0
    for i,r in enumerate(inter_r):
        if rho_deltas:
            rho[(r_0 < radii) & (radii <= r)] *= rho_deltas[i]
        if temp_deltas:
            t[(r_0 < radii) & (radii <= r)] *= temp_deltas[i]
        if epsbig_deltas:
            epsbig[(r_0 < radii) & (radii <= r)] *= epsbig_deltas[i]
        if eps_deltas:
            eps[(r_0 < radii) & (radii <= r)] *= eps_deltas[i]
        r_0 = r

    # We read the file to find where each radius begins
    f = open(fort14,'r')
    # We define the string to parse the file
    string = "irad,err,error="
    lines = f.readlines()
    indexes = [] # position of starting line for each radius
    for i,line in enumerate(lines):
        if string in line:
            indexes.append(i) # We save the line where the string was found
    f.close()
    rstar = float(lines[indexes[0]+1].split()[4])/1.5e13 #in au
    nrad = len(indexes)
    nz = indexes[1] - indexes[0] - 4 # there are four "dummy" lines for each radii

    # We save the new structure
    f = open(fort14.replace('irr','mod'),'w')
    # First dummy line
    f.write(lines[0])
    for i,ind in enumerate(indexes):
        # Each radius starts with four lines
        f.write(lines[ind])
        f.write(lines[ind+1])
        f.write(lines[ind+2])
        f.write(lines[ind+3])
        for j in range(nz):
            f.write("%f %g %g %g %g %f \n" %(z[i,j]/rstar,p[i,j],t[i,j],rho[i,j],epsbig[i,j],eps[i,j]))
    f.close()

    return

max_stars_repo_name: zihuaxin/EDGE

max_issues_repo_name: zihuaxin/EDGE

max_forks_repo_name: zihuaxin/EDGE


text:
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import sympy as sp

# Calculando puntos de intersección
#~ x = sp.Symbol('x')
#~ x = sp.solve(x**2 - x - 2)

x = np.linspace(-1,2)
y1 = x + 3
y2 = x**2 + 1

#~ plt.plot(x,y1)
#~ plt.plot(x,y2)

#~ plt.show()

step = 0.04
maxval = 1.0

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

# create supporting points in polar coordinates
r = np.linspace(-1, 2, 20)
p = np.linspace(0, 2*np.pi, 20)
R, P = np.meshgrid(r, p)

# transform them to cartesian system
X, Y = R*np.cos(P), R*np.sin(P)

Z = R**2
ax.plot_surface(X, Y, Z, rstride=1, cstride=1)

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

reason: import numpy,import sympy


text:
"""Sobol sequence generator

This is a reimplementation of a C++ algorithm by
`Stephen Joe and Frances Y. Kuo <http://web.maths.unsw.edu.au/~fkuo/sobol/>`_.
Directions are based on :file:`new-joe-kuo-6.21201` from the URL above.
"""
import lzma
from typing import Iterable, Tuple, Union, cast

import dask.array as da
import numpy as np
import pkg_resources
from dask.array.core import normalize_chunks
from numba import jit

from .typing import Chunks2D, NormalizedChunks2D

DIRECTIONS = "new-joe-kuo-6.21201.txt.xz"

v_cache = None

def _load_v() -> np.ndarray:
    """Load V from the original author's file. This function is executed
    automatically the first time you call the :func:`sobol` function.
    When using a dask backend, the V array is only loaded when
    actually needed by the kernel; this results in smaller pickle files.
    When using dask distributed, V is loaded locally on the workers instead of
    being transferred over the network.
    """
    global v_cache
    if v_cache is None:
        with pkg_resources.resource_stream("pyscenarios", DIRECTIONS) as fh:
            with lzma.open(fh, "rt") as zfh:
                directions = _load_directions(zfh)  # type: ignore
        v_cache = _calc_v_kernel(directions)
    return v_cache

def _load_directions(fh: Iterable[str]) -> np.ndarray:
    """Load input file containing direction numbers.
    The file must one of those available on the website of the
    original author, or formatted like one.

    :returns:
        contents of the input file as a 2D numpy array.
        Every row contains the information for the matching dimension.
        Column 0 contains the a values, while columns 1+ contain the m values.
        The m values are padded on the right with zeros.
    """
    rows = [row.split() for row in fh]

    # Add padding at end of rows
    # Drop first 2 columns
    # Replace header with element for d=1
    rowlen = len(rows[-1])
    for row in rows:
        row[:] = row[2:] + ["0"] * (rowlen - len(row))
    rows[0] = ["0"] + ["1"] * (rowlen - 3)
    return np.array(rows, dtype="uint32")

@jit("uint32[:,:](uint32[:,:])", nopython=True, nogil=True, cache=True)
def _calc_v_kernel(directions: np.ndarray) -> np.ndarray:
    """Numba kernel for :func:`calc_v`"""
    # Initialise temp array of direction numbers
    v = np.empty((directions.shape[0], 32), dtype=np.uint32)

    for j in range(directions.shape[0]):
        s = 0
        # Compute direction numbers
        for s in range(directions.shape[1] - 1):
            if directions[j, s + 1] == 0:
                break
            v[j, s] = directions[j, s + 1] * 2 ** (31 - s)
        else:
            # need a C-style for loop
            # for(s=0; s<m.size; s++)
            # where at the end of the loop s == m.size
            s += 1

        for t in range(s, 32):
            v[j, t] = v[j, t - s] ^ (v[j, t - s] // 2 ** s)
            for k in range(1, s):
                v[j, t] ^= ((directions[j, 0] // 2 ** (s - 1 - k)) & 1) * v[j, t - k]

    return v

def _sobol_kernel(samples: int, dimensions: int, s0: int, d0: int) -> np.ndarray:
    """Numba kernel for :func:`sobol`

    :returns:
        points 2D array.

        points[i, j] = the jth component of the ith point
        with i indexed from 0 to N-1 and j indexed from 0 to D-1
    """
    output = np.empty((samples, dimensions), order="F")
    _sobol_kernel_jit(samples, dimensions, s0, d0, _load_v(), output)
    return output

@jit(
    "void(uint32, uint32, uint32, uint32, uint32[:, :], float64[:, :])",
    nopython=True,
    nogil=True,
    cache=True,
)
def _sobol_kernel_jit(
    samples: int, dimensions: int, s0: int, d0: int, V: np.ndarray, output: np.ndarray
) -> None:
    """Jit-compiled core of :func:`_sobol_kernel

    When running in dask and there are multiple chunks on the
    samples, calculate and then discard the first s0 samples.
    This is inefficient but preferable to transferring a state
    vector across the graph, which would introduce cross-chunks dependencies.
    """
    for j in range(dimensions):
        state = 0
        for i in range(s0 + samples):
            # c = index from the right of the first zero bit of i
            c = 0
            mask = 1
            while i & mask:
                mask *= 2
                c += 1

            state ^= V[j + d0, c]
            if i >= s0:
                output[i - s0, j] = np.double(state) / np.double(2 ** 32)

def sobol(
    size: Union[int, Tuple[int, int]], d0: int = 0, chunks: Chunks2D = None
) -> Union[np.ndarray, da.Array]:
    """Sobol points generator based on Gray code order

    :param size:
        number of samples (cannot be greater than :math:`2^{32}`) to extract
        from a single dimension, or tuple (samples, dimensions).
        To guarantee uniform distribution, the number of samples should
        always be :math:`2^{n} - 1`.
    :param int d0:
        first dimension. This can be used as a functional equivalent of a
        a random seed. dimensions + d0 can't be greater than
        :func:`max_sobol_dimensions()` - 1.
    :param chunks:
        If None, return a numpy array.

        If set, return a dask array with the given chunk size.
        It can be anything accepted by dask (a positive integer, a
        tuple of two ints, or a tuple of two tuples of ints) for the output
        shape (see result below). e.g. either ``(16384, 50)`` or
        ``((16384, 16383),  (50, 50, 50))`` could be used together with
        ``size=(32767, 150)``.

        .. note::
           The algorithm is not efficient if there are multiple chunks on axis
           0. However, if you do need them, it is typically better to require
           them here than re-chunking afterwards, particularly if (most of) the
           subsequent algorithm is embarassingly parallel.
    :returns:
        If size is an int, a 1-dimensional array of samples.
        If size is a tuple, a 2-dimensional array POINTS, where
        ``POINTS[i, j]`` is the ith sample of the jth dimension.
        Each dimension is a uniform (0, 1) distribution.
    :rtype:
        If chunks is not None, :class:`dask.array.Array`; else
        :class:`numpy.ndarray`
    """
    if isinstance(size, int):
        samples = size
        dimensions = 1
    else:
        samples, dimensions = size

    if not 0 < samples < 2 ** 32:
        raise ValueError("samples must be between 1 and 2^32")
    if not 0 < dimensions + d0 <= max_sobol_dimensions():
        raise ValueError(
            "(dimensions + d0) must be between 1 and %d" % max_sobol_dimensions()
        )

    if chunks is None:
        res = _sobol_kernel(samples, dimensions, 0, d0)
        if isinstance(size, int):
            res = res[:, 0]
        return res

    # dask-specific code
    chunks = cast(
        NormalizedChunks2D, normalize_chunks(chunks, shape=(samples, dimensions))
    )
    name = "sobol-%d-%d-%d" % (samples, dimensions, d0)
    dsk = {}

    offset_i = 0
    for i, size_i in enumerate(chunks[0]):
        offset_j = 0
        for j, size_j in enumerate(chunks[1]):
            dsk[name, i, j] = (_sobol_kernel, size_i, size_j, offset_i, d0 + offset_j)
            offset_j += size_j
        offset_i += size_i

    res = da.Array(dsk, name=name, dtype=float, chunks=chunks)
    if isinstance(size, int):
        res = res[:, 0]
    return res

def max_sobol_dimensions() -> int:
    """Return number of dimensions available. When invoking :func:`sobol`,
    ``size[1] + d0`` must be smaller than this.
    """
    return _load_v().shape[0]

def max_dimensions() -> int:
    import warnings

    warnings.warn(
        "pyscenarios.sobol.max_dimensions has been moved to "
        "pyscenarios.max_sobol_dimensions",
        DeprecationWarning,
    )
    return max_sobol_dimensions()

max_stars_repo_name: crusaderky/pyscenarios

max_issues_repo_name: crusaderky/pyscenarios

max_forks_repo_name: crusaderky/pyscenarios

reason: import numpy,from numba


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

import logging
import tempfile
from pathlib import Path

import langdetect
import numpy as np
from laservec.LASER.source.embed import *
from laservec.LASER.source.lib.text_processing import Token, BPEfastApply

VECTOR_LEN = 1024

class LaserEncoder(object):
    """
    A wrapper for Facebook's LASER (Language-Agnostic SEntence Representations) toolkit.

    See: https://github.com/facebookresearch/LASER
    """
    def __init__(self):
        self.model_dir = Path(__file__).parent / "LASER" / "models"
        self.encoder_path = self.model_dir / "bilstm.93langs.2018-12-26.pt"
        self.bpe_codes_path = self.model_dir / "93langs.fcodes"

        logging.info(f"Loading LASER encoder from: {self.encoder_path}")
        self.encoder = SentenceEncoder(self.encoder_path,
                                       max_sentences=None,
                                       max_tokens=12000,
                                       sort_kind='mergesort',
                                       cpu=True)

    def vectorize(self, text, lang=None):
        """
        Converts the given text into a 1024-dimensional vector representation
        that is agnostic to the source language. The computed vector is "universal"
        in that it lives in an embedding space that is shared among all the languages
        supported by LASER.

        :param text: The input text, consisting of one or more sentences. This parameter is mandatory.

        :param lang: The ISO 639-1 code of the source language. See: https://github.com/facebookresearch/LASER#supported-languages
                     This parameter is optional; if not specified, it will auto-detected.

        :return: A tuple (embedding, lang) containing the 1024-dimensional vector representation
                of the input text (as a numpy array), and the ISO 639-1 language code of the input text.
                Note: if the input text consists of more than one sentence, the returned
                vector is the mean of the individual sentence vectors.
        """
        if not lang:
            lang = langdetect.detect(text)
        if not lang:
            lang = "en"

        with tempfile.TemporaryDirectory() as tmp:
            temp_dir = Path(tmp)
            input_fname = temp_dir / "input.txt"
            bpe_encoded_fname = temp_dir / 'input.bpe'
            raw_output_fname = temp_dir / 'output.raw'

            with input_fname.open("w") as f:
                f.write(text)

            # Tokenize the text
            if lang != '--':
                tokenized_fname = temp_dir / "input.tok"
                Token(str(input_fname),
                      str(tokenized_fname),
                      lang=lang,
                      romanize=True if lang == 'el' else False,
                      lower_case=True,
                      gzip=False,
                      verbose=True,
                      over_write=False)
                input_fname = tokenized_fname

            # Encode the text using BPE
            BPEfastApply(str(input_fname),
                         str(bpe_encoded_fname),
                         str(self.bpe_codes_path),
                         verbose=True, over_write=False)

            # Finally, compute the embedding vector.
            input_fname = bpe_encoded_fname
            EncodeFile(self.encoder,
                       str(input_fname),
                       str(raw_output_fname),
                       verbose=True,
                       over_write=False,
                       buffer_size=10000)

            raw_embeddings = np.fromfile(str(raw_output_fname), dtype=np.float32, count=-1)
            row_count = raw_embeddings.shape[0] // VECTOR_LEN
            col_count = VECTOR_LEN
            raw_embeddings.resize(row_count, col_count)

            # If the input text contains more than one sentence, a vector is computed for each sentence.
            # Combine the individual sentence vectors by computing their mean.
            embedding = raw_embeddings.mean(axis=0)
            return embedding, lang

max_stars_repo_name: cris-almodovar/laservec-docker

max_issues_repo_name: cris-almodovar/laservec-docker

max_forks_repo_name: cris-almodovar/laservec-docker


text:
"""
Copyright (c) 2020 CRISP

Convolutional Orthogonal Matching Pursuit

:author: Andrew H. Song
"""

import numpy as np
import sys
import pickle
import os
import scipy.signal

# Append path
PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..")
sys.path.append(PATH)

from . import BaseCSC
from src.helpers.convolution import sc_cholesky, sc_cholesky_efficient

class COMP(BaseCSC):
	def __init__(self, dlen, error_tol, sparsity_tol, pflag):
		BaseCSC.__init__(self, dlen, error_tol, sparsity_tol, pflag)

	def terminate(self, numOfiter, numOfmaxcoeffs, err_residual, err_bound):
		return (err_residual< err_bound) or (numOfiter >= numOfmaxcoeffs)

	def extractCode_seg(self, y_seg, dictionary, sparsity=None, err=None, boundary=1):

		"""
		Given data segment, extract convolutional codes

		Inputs
		======
		y_seg: A segment of data
		boundary: boolean (default=1)
			If 1, accounts for truncated templates as well (clen = slen + dlen - 1)
			If 0, accounts for whole templates only (cldn = slen - dlen + 1)

		"""
		assert(len(np.where(abs(np.linalg.norm(dictionary,axis=0)-1)>1e-6)[0]) == 0), "Not normalized"

		numOfsamples, numOfelements = dictionary.shape

		slen = len(y_seg)
		clen = slen + numOfsamples - 1
		coeffs = np.zeros(clen*numOfelements)

		numOfmaxcoeffs = self.sparsity_tol if sparsity is None else sparsity
		err_bound = self.error_tol if err is None else err

		chosen_vals = np.zeros(numOfelements)
		chosen_idx = np.zeros(numOfelements, dtype=np.int)

		residual = y_seg
		err_residual = np.linalg.norm(residual)/np.sqrt(np.size(residual))

		# Dictionary to collect expanding set of dictionary
		temp_idx = np.zeros(numOfmaxcoeffs, dtype=np.int)
		dictionary_active = np.zeros((slen, numOfmaxcoeffs))

		iternum = 0
		lower_mat = [1]

		while not self.terminate(iternum, numOfmaxcoeffs, err_residual, err_bound):
			#######################
			# Selection step
			#
			# This step can be fruther sped up with careful book-keeping of residuals
			#######################
			for idx in np.arange(numOfelements):
				d = dictionary[:, idx]
				cross = abs(scipy.signal.correlate(residual, d, mode='full'))

				if boundary:
					t_start = int(np.floor(numOfsamples/2))
					t_end = slen + t_start
				else:
					t_start = numOfsamples-1
					t_end = slen - t_start

				cross = cross/self.computeNorm(d, slen)
				m = np.argmax(cross[t_start:t_end]) + t_start

				chosen_idx[idx] = m
				chosen_vals[idx] = cross[m]

			filter_idx = np.argmax(chosen_vals) # Returns the filter with the highest inner product
			coeff_idx = chosen_idx[filter_idx] # index within the chosen filter

			#######################
			# Projection step
			#######################

			# placeholder for coefficients
			temp_idx[iternum] = filter_idx*clen + coeff_idx

			if coeff_idx < numOfsamples-1:	# Boundary condition
				offset = coeff_idx + 1
				elem = dictionary[-offset:,filter_idx]/np.linalg.norm(dictionary[-offset:,filter_idx])
				elem = np.pad(elem, (0, slen-len(elem)), 'constant')
			elif coeff_idx > slen - 1:	# Boundary condition
				offset = numOfsamples-(coeff_idx-(slen-1))
				elem = dictionary[:offset,filter_idx]/np.linalg.norm(dictionary[:offset,filter_idx])
				elem = np.pad(elem, (slen-len(elem), 0), 'constant')
			else:	# Valid correlation. Entire support of the dictionary lies within the signal
				start_idx = coeff_idx - (numOfsamples-1)
				elem = dictionary[:, filter_idx]/np.linalg.norm(dictionary[:, filter_idx])
				elem = np.pad(elem, (start_idx, slen - numOfsamples - start_idx), 'constant')

			dictionary_active[:, iternum] = elem
			[lower_mat, sparse_code] = sc_cholesky(lower_mat, dictionary_active[:, :iternum], elem, y_seg)

			residual = y_seg - np.matmul(dictionary_active[:, :iternum+1], sparse_code)
			coeffs[temp_idx[:iternum+1]] = sparse_code

			iternum += 1

		err_residual = np.linalg.norm(residual)
		return coeffs, err_residual

	def extractCode_seg_eff(self, y_seg, dictionary, sparsity=None, err=None, boundary=1):

		"""
		Given data segment, extract convolutional codes
		Uses efficient implementation of cholesky decomposition, which should be much faster for long signals

		Inputs
		======
		y_seg: A segment of data
		boundary: boolean (default=1)
			If 1, accounts for truncated templates as well (clen = slen + dlen - 1)
			If 0, accounts for whole templates only (cldn = slen - dlen + 1)
		"""
		assert(len(np.where(abs(np.linalg.norm(dictionary,axis=0)-1)>1e-6)[0]) == 0), "Not normalized"

		numOfsamples, numOfelements = dictionary.shape

		slen = len(y_seg)
		clen = slen + numOfsamples - 1
		coeffs = np.zeros(clen*numOfelements)

		numOfmaxcoeffs = self.sparsity_tol if sparsity is None else sparsity
		err_bound = self.error_tol if err is None else err

		chosen_vals = np.zeros(numOfelements)
		chosen_idx = np.zeros(numOfelements, dtype=np.int)

		residual = y_seg
		err_residual = np.linalg.norm(residual)/np.sqrt(np.size(residual))

		# Dictionary to collect expanding set of dictionary
		temp_idx = np.zeros(self.sparsity_tol, dtype=np.int)

		# Active index sets
		filter_indices = []
		start_indices = []

		iternum = 0
		lower_mat = [1]

		while not self.terminate(iternum, numOfmaxcoeffs, err_residual, err_bound):
			#######################
			# Selection step
			#
			# This step can be fruther sped up with careful book-keeping of residuals
			#######################
			for idx in np.arange(numOfelements):
				d = dictionary[:, idx]
				cross = abs(scipy.signal.correlate(residual, d, mode='full'))

				if boundary:
					t_start =int(np.floor(numOfsamples/2))
					t_end = slen + t_start
				else:
					t_start = numOfsamples-1
					t_end = slen - t_start

				cross = cross/self.computeNorm(d, slen)

				m = np.argmax(cross[t_start:t_end]) + t_start

				chosen_idx[idx] = m
				chosen_vals[idx] = cross[m]

			filter_idx = np.argmax(chosen_vals) # Returns the filter with the highest inner product
			filter_indices.append(filter_idx)
			coeff_idx = chosen_idx[filter_idx] # index within the chosen filter

			#######################
			# Projection step
			#######################

			# placeholder for coefficients
			temp_idx[iternum] = filter_idx*clen + coeff_idx

			# Pad the dictionary element to match the signal length
			start_idx = coeff_idx - numOfsamples + 1
			start_indices.append(start_idx)

			[lower_mat, sparse_code] = sc_cholesky_efficient(lower_mat, dictionary, start_indices, filter_indices, y_seg)

			residual = np.copy(y_seg)
			for i, filter_idx in enumerate(filter_indices):
				elem = dictionary[:, filter_idx]
				start_idx = start_indices[i]
				coeff = sparse_code[i]

				if start_idx < 0:
					residual[:numOfsamples+start_idx] -= coeff * elem[-(numOfsamples+start_idx):]
				elif start_idx > slen - numOfsamples:
					residual[start_idx:] -= coeff * elem[:slen - start_idx]
				else:
					residual[start_idx:start_idx + numOfsamples] -= coeff*elem

			coeffs[temp_idx[:iternum+1]] = sparse_code

			iternum += 1

		err_residual = np.linalg.norm(residual)
		return coeffs, err_residual

max_stars_repo_name: DominicTanzillo/srcdl

max_issues_repo_name: DominicTanzillo/srcdl

max_forks_repo_name: DominicTanzillo/srcdl


text:
import numpy as np
from GPy.util import choleskies
import GPy
from ..util.config import config
import unittest

try:
    from ..util import linalg_cython
    from ...configs_that_do_not_run.GPy_0_8_8.util import choleskies_cython

    config.set('cython', 'working', 'True')
except ImportError:
    config.set('cython', 'working', 'False')

"""
These tests make sure that the pure python and cython codes work the same
"""

@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
class CythonTestChols(np.testing.TestCase):
    def setUp(self):
        self.flat = np.random.randn(45,5)
        self.triang = np.array([np.eye(20) for i in range(3)])
    def test_flat_to_triang(self):
        L1 = choleskies._flat_to_triang_pure(self.flat)
        L2 = choleskies._flat_to_triang_cython(self.flat)
        np.testing.assert_allclose(L1, L2)
    def test_triang_to_flat(self):
        A1 = choleskies._triang_to_flat_pure(self.triang)
        A2 = choleskies._triang_to_flat_cython(self.triang)
        np.testing.assert_allclose(A1, A2)

@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
class test_stationary(np.testing.TestCase):
    def setUp(self):
        self.k = GPy.kern.RBF(10)
        self.X = np.random.randn(300,10)
        self.Z = np.random.randn(20,10)
        self.dKxx = np.random.randn(300,300)
        self.dKzz = np.random.randn(20,20)
        self.dKxz = np.random.randn(300,20)

    def test_square_gradX(self):
        g1 = self.k._gradients_X_cython(self.dKxx, self.X)
        g2 = self.k._gradients_X_pure(self.dKxx, self.X)
        np.testing.assert_allclose(g1, g2)

    def test_rect_gradx(self):
        g1 = self.k._gradients_X_cython(self.dKxz, self.X, self.Z)
        g2 = self.k._gradients_X_pure(self.dKxz, self.X, self.Z)
        np.testing.assert_allclose(g1, g2)

    def test_square_lengthscales(self):
        g1 = self.k._lengthscale_grads_pure(self.dKxx, self.X, self.X)
        g2 = self.k._lengthscale_grads_cython(self.dKxx, self.X, self.X)
        np.testing.assert_allclose(g1, g2)

    def test_rect_lengthscales(self):
        g1 = self.k._lengthscale_grads_pure(self.dKxz, self.X, self.Z)
        g2 = self.k._lengthscale_grads_cython(self.dKxz, self.X, self.Z)
        np.testing.assert_allclose(g1, g2)

@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
class test_choleskies_backprop(np.testing.TestCase):
    def setUp(self):
        a =np.random.randn(10,12)
        A = a.dot(a.T)
        self.L = GPy.util.linalg.jitchol(A)
        self.dL = np.random.randn(10,10)
    def test(self):
        r1 = choleskies._backprop_gradient_pure(self.dL, self.L)
        r2 = choleskies_cython.backprop_gradient(self.dL, self.L)
        r3 = choleskies_cython.backprop_gradient_par_c(self.dL, self.L)
        np.testing.assert_allclose(r1, r2)
        np.testing.assert_allclose(r1, r3)