text:
from itertools import product

import numpy as np
from typing import Tuple
from IMLearn.learners.classifiers import DecisionStump
from IMLearn.metalearners import AdaBoost
from utils import *
import plotly.graph_objects as go
from plotly.subplots import make_subplots

pio.renderers.default = "browser"

def generate_data(n: int, noise_ratio: float) -> Tuple[np.ndarray, np.ndarray]:
    """
    Generate a dataset in R^2 of specified size

    Parameters
    ----------
    n: int
        Number of samples to generate

    noise_ratio: float
        Ratio of labels to invert

    Returns
    -------
    X: np.ndarray of shape (n_samples,2)
        Design matrix of samples

    y: np.ndarray of shape (n_samples,)
        Labels of samples
    """
    '''
    generate samples X with shape: (num_samples, 2) and labels y with shape (num_samples).
    num_samples: the number of samples to generate
    noise_ratio: invert the label for this ratio of the samples
    '''
    X, y = np.random.rand(n, 2) * 2 - 1, np.ones(n)
    y[np.sum(X ** 2, axis=1) < 0.5 ** 2] = -1
    y[np.random.choice(n, int(noise_ratio * n))] *= -1
    return X, y

def add_partial_decision_boundary(fig, X, y, t, learner, lims, row=None, col=None):
    """
    Plot the decision boundary of ensemble with t estimators
    """
    # symbols = np.array(["circle", "x"])[((y + 1) / 2).astype(int)]
    predict = lambda X_: learner.partial_predict(X_, t)
    accuracy = 1 - learner.partial_loss(X, y, t)

    fig.add_trace(decision_surface(predict, lims[0], lims[1], showscale=False),
                  row=row, col=col)

    class0 = y == -1
    fig.add_trace(go.Scatter(x=X[class0][:, 0], y=X[class0][:, 1], mode="markers",
                             name="Class -1", legendgroup='Class -1', showlegend=False,
                             marker=dict(color="red", symbol="circle", line=dict(color="black", width=1))),
                  row=row, col=col)

    class1 = y == 1
    fig.add_trace(go.Scatter(x=X[class1][:, 0], y=X[class1][:, 1], mode="markers",
                             name="Class 1", legendgroup='Class 1', showlegend=False,
                             marker=dict(color="blue", symbol="x", line=dict(color="black", width=1))),
                  row=row, col=col)

    fig.update_xaxes(title_text="x", row=row, col=col)
    fig.update_yaxes(title_text="y", row=row, col=col)
    if row is None:
        fig.update_layout(title_text=f"Decision boundary of ensemble with {t} estimators, Accuracy: {accuracy:.3f}")
    else:
        fig.layout.annotations[2*(row-1)+col-1].update(text=f"Using {t} estimators, Accuracy: {accuracy: .2f}")

    return fig

def fit_and_evaluate_adaboost(noise, n_learners=250, train_size=5000, test_size=500):
    (train_X, train_y), (test_X, test_y) = generate_data(train_size, noise), generate_data(test_size, noise)

    # Question 1: Train- and test errors of AdaBoost in noiseless case
    learner = AdaBoost(DecisionStump, n_learners)
    learner.fit(train_X, train_y)

    # Plot the training- and test errors as a function of the number of fitted learners
    num_of_learners = np.arange(1, n_learners + 1)
    train_errors = [learner.partial_loss(train_X, train_y, t) for t in num_of_learners]
    test_errors = [learner.partial_loss(test_X, test_y, t) for t in num_of_learners]
    fig = go.Figure(data=[go.Scatter(x=num_of_learners, y=train_errors, name="Training Error"),
                          go.Scatter(x=num_of_learners, y=test_errors, name="Test Error")],
                    layout=go.Layout(title='Training and test errors as a function of the number of fitted learners',
                                     xaxis_title='Number of learners',
                                     yaxis_title='Error'))
    fig.show()

    # Question 2: Plotting decision surfaces
    T = [[5, 50], [100, 250]]
    lims = np.array([np.r_[train_X, test_X].min(axis=0),
                     np.r_[train_X, test_X].max(axis=0)]).T + np.array([-.1, .1])

    fig = make_subplots(rows=2, cols=2, specs=2 * [2 * [{"type": "scatter"}]],
                        subplot_titles=4*["Decision surface"],
                        vertical_spacing=0.15,
                        horizontal_spacing=0.10)

    fig.update_layout(title_text=f"Decision boundary of the ensemble",
                      xaxis_title="x", yaxis_title="y", legend_title_text='Test Set',
                      margin_t=50)

    for row, col in product(range(2), range(2)):
        add_partial_decision_boundary(fig, test_X, test_y, T[row][col], learner, lims, row=row+1, col=col+1)

    fig.data[1].showlegend = True
    fig.data[2].showlegend = True
    fig.show()

    # Question 3: Decision surface of best performing ensemble
    best_T = num_of_learners[np.argmin(test_errors)]
    fig = go.Figure(layout=go.Layout(legend_title_text='Test Set'))
    add_partial_decision_boundary(fig, test_X, test_y, best_T, learner, lims)
    fig.data[1].showlegend = True
    fig.data[2].showlegend = True
    fig.show()

    # Question 4: Decision surface with weighted samples
    fig = go.Figure(layout=go.Layout(title=f'Decision boundary of fitted model, with train set',
                                     xaxis=dict(title='x'),
                                     yaxis=dict(title='y')))

    fig.add_trace(decision_surface(learner.predict, lims[0], lims[1], showscale=False))
    max_bubble_size = 50 if noise == 0 else 5
    sizeref = 2. * max(learner.D_) / (max_bubble_size ** 2)
    fig.add_trace(go.Scatter(x=train_X[:, 0], y=train_X[:, 1], mode="markers",
                             marker=dict(color=train_y,
                                         colorscale=[custom[0], custom[-1]],
                                         size=learner.D_,
                                         sizemode='area',
                                         sizeref=sizeref,
                                         sizemin=0.5
                                         )
                             ))
    fig.show()

if __name__ == '__main__':
    np.random.seed(0)
    fit_and_evaluate_adaboost(0)
    fit_and_evaluate_adaboost(0.4)


text:
import numpy as np
import parameter
# Define a class to receive the characteristics of each line detection

class Line():
    def __init__(self):
        # was the line detected in the last iteration?
        self.detected = False
        # x values of the last n fits of the line
        self.recent_xfitted = []
        # average x values of the fitted line over the last n iterations
        self.bestx = None

        # polynomial coefficients averaged over the last n iterations
        self.best_fit = None
        # polynomial coefficients for the most recent fit
        self.current_fit = []

        # difference in fit coefficients between last and new fits
        self.diffs = np.array([0, 0, 0], dtype='float')
        # x values for detected line pixels
        self.allx = None
        # y values for detected line pixels
        self.ally = None

        # radius of curvature of the line in some units
        self.radius_of_curvature = None
        # distance in meters of vehicle center from the line
        self.line_base_pos = None
        self.not_detected_counter = 0

    def check_plausibility(self, lane_fit):
        # check
        return True

    def calc_curvature(self):
        ploty = np.linspace(0, 719, 720)
        y_eval = np.max(ploty)

        # Calc new polynomial coefficients in world space
        if self.best_fit is not None:
            fitA = (parameter.xm_per_pix / parameter.ym_per_pix ** 2) * self.best_fit[0]
            fitB = (parameter.xm_per_pix / parameter.ym_per_pix) * self.best_fit[1]

        # Calculate radius in meter
        self.radius_of_curvature = (1+(2*fitA*y_eval*parameter.ym_per_pix+fitB)**2)**(3/2)/np.absolute(2*fitA)

    def update_best_fit(self, lane_fit):
        if not self.check_plausibility(lane_fit):
            self.not_detected_counter = self.not_detected_counter + 1

        elif lane_fit is not None:
            if self.best_fit is not None:
                self.diffs = abs(lane_fit - self.best_fit)

            if self.diffs[0] > 0.001 or self.diffs[1] > 1.0 or self.diffs[2] > 200.:
                print(self.diffs)
                print("difference too much")
                self.not_detected_counter = self.not_detected_counter + 1
            else:
                self.detected = True

                self.not_detected_counter = 0
                self.current_fit.append(lane_fit)
                if len(self.current_fit) > 5:
                    self.current_fit = self.current_fit[len(self.current_fit)-5:]
                self.best_fit = np.average(self.current_fit, axis=0)
                self.calc_curvature()
        else:
            self.not_detected_counter = self.not_detected_counter + 1

        if self.not_detected_counter >= 2:
            self.detected = False
            self.best_fit = None
            self.diffs = np.array([0, 0, 0], dtype='float')
            self.current_fit = []
            self.not_detected_counter = 0
            print("reset")

ext: py

reason: import numpy


text:
# -*- coding: utf-8 -*-
import re
import io
import time
import numpy
import pandas
import sqlite3
import datetime
import requests
import lxml.html as LH
import xml.etree.ElementTree as xmltree
import Lily.ctao.hostmetadata as chmd
import Lily.ctao.database as cdb

def cwb_melt0():
    import lzma
    #db = cdb.database('data_crawler_cwb_earthquake_list.sqlite')
    #df = db.to_dataframe('data_crawler_cwb_earthquake_list')
    #df = df.reindex ( columns=['id', 'html', 'lzma_html'], fill_value='' )

    #for ind, row in df.iterrows():
    #    df.at[ind,'lzma_html'] = lzma.compress(df.at[ind,'html'])
    #    print (df.at[ind,'id'] )

    #df = df.drop(columns=['html'])
    #df.to_sql('data_crawler_cwb_sensible_earthquake_download', db.connect, if_exists='replace', index=False)
    return

def cwb_melt1():
    import lzma
    db = cdb.database('data_crawler.sqlite')

    sql = '''
        select Id, routeId,  nameZh, nameEn, seqNo, pgp, longitude, showLon, showLat, vector from  group by id
    '''
    df = pandas.read_sql(sql, db.connect , index_col=['rowid'])

    df = df.reindex ( columns=[ 'a','b','c', 'd', 'e', 'f', 'g', 'h','i' ], fill_value=''  )

    for ind, row in df.iterrows():
#        print ('melt', row[0])
        json_tables = pandas.read_html( lzma.decompress( sqlite3.Binary(row['lzma_html']) ), encoding='utf-8')
        arg2 = json_tables[2]
        df.at[ind,'routeId'] = arg2.iat[0,1]
        df.at[ind,'nameZh'] = arg2.iat[1,1]
        df.at[ind,'nameEn'] = arg2.iat[2,1]
        df.at[ind,'seqNo'] = arg2.iat[3,1]
        df.at[ind,'pgp'] = arg2.iat[4,1]
        df.at[ind,'longitude'] = arg2.iat[5,1]

        for ind2, row2 in json_tables[3].iterrows():
            if isinstance(row2, pandas.core.series.Series):
                for elem in row2:
                    if isinstance(elem,str):
                        df.at[ind, 'i'] = df.at[ind, 'i'] + ',' + elem
            else:
                if isinstance(elem,str):
                    df.at[ind, 'i'] = df.at[ind, 'i'] + ',' + row2

    #df = df.drop(columns=['lzma_html'])
    df.to_sql('data_rdset_pylily', db.connect, if_exists='append', index=False)

    db.connect.execute('''delete from {0} where rowid not in
                         (select max (rowid) from {0} group by id)'''.format('data_rdset_pylily') )

    db.connect.commit()
    return

def cwb_melt2():
    db = cdb.database('data_crawler.sqlite')
    df = db.to_dataframe('data_rdset_pylily_cwb_sensible_earthquake')

    station = []
    for ind, row in df.iterrows():

        for st in df.at[ind, 'Stations'].split(';'):
            if u'''地區最大震度''' not in st and st != '':
                rdset = [df.at[ind, 'id'],
                         df.at[ind, 'time'],
                         float(df.at[ind, 'px'][4:-2]),
                         float(df.at[ind, 'py'][4:-2]),
                         float(df.at[ind, 'depth'][:-3]),
                         float(df.at[ind, 'ML']),
                         df.at[ind, 'Location'],
                         ''.join(st.split('\u3000')[:-1]) , float(st.split('\u3000')[-1:][0]) ]
                station.append(rdset)
    df2 = pandas.DataFrame(station, columns=['a','b','c','d', 'e','f','g', 'h', 'i'])
    df2.to_sql('data_rdset_pylily_cwb_sensible_earthquake_LocalSeismicIntensity', db.connect, if_exists='replace', index=False)
    return

if __name__ == '__console__' or __name__ == '__main__':
    import os
    thost = chmd.hostmetadata()
    os.chdir (thost.warehouse)
    cwb_crawler()
    cwb_melt1()
    cwb_melt2()


text:
import numpy as np
import os

import subprocess
import base64
import glob
import io
import json

from gfootball.env.script_helpers import ScriptHelpers as sh

def get_gif_html(videos_path, title, subtitle_eps=None, max_n_videos=4):
    videos = np.array(glob.glob(videos_path))
    dumps = np.array(glob.glob(videos_path))
    if len(videos) == 0:
        return

    n_videos = max(1, min(max_n_videos, len(videos)))
    idxs = np.linspace(0, len(videos) - 1, n_videos).astype(int) if n_videos > 1 else [-1,]
    videos = videos[idxs,...]
    vids = {}
    for video_path in videos:
        basename = os.path.splitext(video_path)[0]
        gif_path = basename + '.gif'
        dump_path = basename + '.dump'

        if not os.path.exists(gif_path):
            ps = subprocess.Popen(
                ('ffmpeg',
                 '-i', video_path,
                 '-r', '10',
                 '-f', 'image2pipe',
                 '-vcodec', 'ppm',
                 '-crf', '20',
                 '-vf', 'scale=300:-1',
                 '-'),
                stdout=subprocess.PIPE)
            output = subprocess.check_output(
                ('convert',
                 '-coalesce',
                 '-delay', '7',
                 '-loop', '0',
                 '-fuzz', '2%',
                 '+dither',
                 '-deconstruct',
                 '-layers', 'Optimize',
                 '-', gif_path),
                stdin=ps.stdout)
            ps.wait()

        gif = io.open(gif_path, 'r+b').read()
        encoded = base64.b64encode(gif)

        html_tag = """

{0}

        """
        prefix = 'Trial ' if subtitle_eps is None else 'Episode '
        d = sh().load_dump(dump_path)
        episode_id = d[0]['debug']['config']['episode_number']
        sufix = str(episode_id if subtitle_eps is None \
                    else subtitle_eps[episode_id])
        vids[episode_id] = html_tag.format(prefix + sufix, encoded.decode('ascii'))

    strm = '
{}
'.format(title)
    for _, v in sorted(vids.items()):
        strm += v
    return strm

def get_gif_html_oai(env_videos, title, subtitle_eps=None, max_n_videos=4):
    videos = np.array(env_videos)
    if len(videos) == 0:
        return

    n_videos = max(1, min(max_n_videos, len(videos)))
    idxs = np.linspace(0, len(videos) - 1, n_videos).astype(int) if n_videos > 1 else [-1,]
    videos = videos[idxs,...]

    strm = '
{}
'.format(title)
    for video_path, meta_path in videos:
        basename = os.path.splitext(video_path)[0]
        gif_path = basename + '.gif'
        if not os.path.exists(gif_path):
            ps = subprocess.Popen(
                ('ffmpeg',
                 '-i', video_path,
                 '-r', '7',
                 '-f', 'image2pipe',
                 '-vcodec', 'ppm',
                 '-crf', '20',
                 '-vf', 'scale=512:-1',
                 '-'),
                stdout=subprocess.PIPE)
            output = subprocess.check_output(
                ('convert',
                 '-coalesce',
                 '-delay', '7',
                 '-loop', '0',
                 '-fuzz', '2%',
                 '+dither',
                 '-deconstruct',
                 '-layers', 'Optimize',
                 '-', gif_path),
                stdin=ps.stdout)
            ps.wait()

        gif = io.open(gif_path, 'r+b').read()
        encoded = base64.b64encode(gif)

        with open(meta_path) as data_file:
            meta = json.load(data_file)

        html_tag = """

{0}

        """
        prefix = 'Trial ' if subtitle_eps is None else 'Episode '
        sufix = str(meta['episode_id'] if subtitle_eps is None \
                    else subtitle_eps[meta['episode_id']])
        strm += html_tag.format(prefix + sufix, encoded.decode('ascii'))
    return strm

max_stars_repo_name: mimoralea/gfootball-intro

max_issues_repo_name: mimoralea/gfootball-intro

max_forks_repo_name: mimoralea/gfootball-intro


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import logging
import os
import h5py
import json_tricks as json
from tqdm import tqdm
from collections import OrderedDict

import torch
import numpy as np
from scipy.io import loadmat, savemat
from poseutils.constants import *

from dataset.JointsDataset import JointsDataset

logger = logging.getLogger(__name__)

'''
keypoints: {
    0: 'Hip',
    1: 'RHip',
    2: 'RKnee',
    3: 'RAnkle',
    4: 'LHip',
    5: 'LKnee',
    6: 'LAnkle',
    7: 'Neck',
    8: 'LUpperArm',
    9: 'LElbow',
    10: 'LWrist',
    11: 'RUpperArm',
    12: 'RElbow',
    13: 'RWrist'
}
edges: [[0, 1], [0, 4], [0, 7], [1, 2], [2, 3], [4, 5], [5, 6], [7, 8], [8, 9], [9, 10], [7, 11], [11, 12], [12, 13]]
'''

class H36MDataset(JointsDataset):
    def __init__(self, cfg, root, image_set, is_train, transform=None):
        super().__init__(cfg, root, image_set, is_train, transform)

        self.num_joints = 14
        self.flip_pairs = [[1, 4],[2, 5],[3, 6],[11, 8],[12, 9],[13, 10]]
        self.parent_ids = None
        self.aspect_ratio = 1.0
        self.pixel_std = 200

        self.upper_body_ids = (7, 8, 9, 10, 11, 12, 13)
        self.lower_body_ids = (0, 1, 2, 3, 4, 5, 6)

        self.db = self._get_db()

        if is_train and cfg.DATASET.SELECT_DATA:
            self.db = self.select_data(self.db)

        logger.info('=> load {} samples'.format(len(self.db)))

    def _get_db(self):
        # create train/val split

        file_name = os.path.join(self.root, f"meta_{self.image_set}.h5")
        image_dir = os.path.join(self.root, "images", self.image_set)

        gt_db = []

        with h5py.File(file_name) as file_:

            key_count = len(list(file_.keys()))

            for _, key in tqdm(enumerate(file_), total=key_count):
                entry_2d = file_[key]["cropped_2d"][:]
                entry_imgs = file_[key]["image_file_names"][:]
                entry_bins = file_[key]["bins"][:]
                total_frames = len(entry_imgs)

                for i_frame in range(total_frames):

                    center, scale = self._box2cs((0, 0, 256, 256))

                    gt_db.append({
                        "image": os.path.join(image_dir, entry_imgs[i_frame].decode("utf-8")),
                        "joints_3d": entry_2d[i_frame],
                        "center": center,
                        "scale": scale,
                        "joints_3d_vis": np.ones_like(entry_2d[i_frame]),
                        "bins": entry_bins[i_frame]
                    })

            file_.close()

        print(f"{len(gt_db)} data points loaded.")

        return gt_db

    def _box2cs(self, box):
        x, y, w, h = box[:4]
        return self._xywh2cs(x, y, w, h)

    def _xywh2cs(self, x, y, w, h):
        center = np.zeros((2), dtype=np.float32)
        center[0] = x + w * 0.5
        center[1] = y + h * 0.5

        if w > self.aspect_ratio * h:
            h = w * 1.0 / self.aspect_ratio
        elif w < self.aspect_ratio * h:
            w = h * self.aspect_ratio
        scale = np.array(
            [w * 1.0 / self.pixel_std, h * 1.0 / self.pixel_std],
            dtype=np.float32)
        if center[0] != -1:
            scale = scale * 1.25

        return center, scale

    def __getitem__(self, idx):
        input, target, target_weight, meta = super().__getitem__(idx)

        return input, target, target_weight, meta, torch.from_numpy(self.db[idx]["bins"])

    def evaluate(self, cfg, preds, output_dir, *args, **kwargs):
        # convert 0-based index to 1-based index
        preds = preds[:, :, 0:2] + 1.0

        if output_dir:
            pred_file = os.path.join(output_dir, 'pred.mat')
            savemat(pred_file, mdict={'preds': preds})

        if 'test' in cfg.DATASET.TEST_SET:
            return {'Null': 0.0}, 0.0

        SC_BIAS = 0.6
        threshold = 0.5

        jnt_missing = np.zeros((self.num_joints, preds.shape[0]))
        pos_gt_src = np.vstack([self.db[i]["joints_3d"].reshape(-1, self.num_joints, 2) for i in range(len(self.db))])
        pos_gt_src = pos_gt_src.transpose(1, 2, 0)

        pos_pred_src = np.transpose(preds, [1, 2, 0])

        head = 7
        lsho = 8
        lelb = 9
        lwri = 10
        lhip = 4
        lkne = 5
        lank = 6

        rsho = 11
        relb = 12
        rwri = 13
        rkne = 2
        rank = 3
        rhip = 0

        jnt_visible = 1 - jnt_missing
        uv_error = pos_pred_src - pos_gt_src
        uv_err = np.linalg.norm(uv_error, axis=1)
        scale = np.ones((len(uv_err), 1))
        scaled_uv_err = np.divide(uv_err, scale)
        scaled_uv_err = np.multiply(scaled_uv_err, jnt_visible)
        jnt_count = np.sum(jnt_visible, axis=1)
        less_than_threshold = np.multiply((scaled_uv_err <= threshold),
                                          jnt_visible)
        PCKh = np.divide(100.*np.sum(less_than_threshold, axis=1), jnt_count)

        # save
        rng = np.arange(0, 0.5+0.01, 0.01)
        pckAll = np.zeros((len(rng), self.num_joints))

        for r in range(len(rng)):
            threshold = rng[r]
            less_than_threshold = np.multiply(scaled_uv_err <= threshold,
                                              jnt_visible)
            pckAll[r, :] = np.divide(100.*np.sum(less_than_threshold, axis=1),
                                     jnt_count)

        PCKh = np.ma.array(PCKh, mask=False)
        PCKh.mask[6:8] = True

        jnt_count = np.ma.array(jnt_count, mask=False)
        jnt_count.mask[6:8] = True
        jnt_ratio = jnt_count / np.sum(jnt_count).astype(np.float64)

        name_value = [
            ('Head', PCKh[head]),
            ('Shoulder', 0.5 * (PCKh[lsho] + PCKh[rsho])),
            ('Elbow', 0.5 * (PCKh[lelb] + PCKh[relb])),
            ('Wrist', 0.5 * (PCKh[lwri] + PCKh[rwri])),
            ('Hip', 0.5 * (PCKh[lhip] + PCKh[rhip])),
            ('Knee', 0.5 * (PCKh[lkne] + PCKh[rkne])),
            ('Ankle', 0.5 * (PCKh[lank] + PCKh[rank])),
            ('Mean', np.sum(PCKh * jnt_ratio)),
            ('Mean@0.1', np.sum(pckAll[11, :] * jnt_ratio))
        ]
        name_value = OrderedDict(name_value)

        return name_value, name_value['Mean']

max_stars_repo_name: SaadManzur/HRNet-Human-Pose-Estimation

max_issues_repo_name: SaadManzur/HRNet-Human-Pose-Estimation

max_forks_repo_name: SaadManzur/HRNet-Human-Pose-Estimation

reason: import numpy,from scipy


text:
from __future__ import print_function
import datetime
import time
import keras
import scipy.io as io
import numpy as np
import os
import scipy.io as sio
import matplotlib.pyplot as plt
from keras.regularizers import l2
from keras.layers import Dense, Dropout, Activation, Flatten, Subtract, Input, BatchNormalization
from keras import backend as K
from keras.models import Sequential, Model
from keras.models import load_model
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.layers.core import Activation
from keras.layers.core import Flatten
from keras.layers.core import Dense
from keras.utils import np_utils
from keras.utils import plot_model
from keras.optimizers import SGD, RMSprop, Adam
import mat73

# loss function
def NMSE_IRS(y_true, y_pred):
    r_true = y_true[:,:,:,0] #shape: (batchsize, 8, 8)
    i_true = y_true[:,:,:,1]
    r_pred = y_pred[:,:,:,0]
    i_pred = y_pred[:,:,:,1]
    mse_r_sum = K.sum(K.sum(K.square(r_pred - r_true), -1), -1)
    r_sum = K.sum(K.sum(K.square(r_true), -1), -1)
    mse_i_sum = K.sum(K.sum(K.square(i_pred - i_true), -1), -1)
    i_sum = K.sum(K.sum(K.square(i_true), -1), -1)
    num = mse_r_sum + mse_i_sum
    den = r_sum + i_sum
    return num/den

class LossHistory(keras.callbacks.Callback):
    def on_train_begin(self, logs={}):
        self.losses = {'batch':[], 'epoch':[]}
        self.accuracy = {'batch':[], 'epoch':[]}
        self.val_loss = {'batch':[], 'epoch':[]}
        self.val_acc = {'batch':[], 'epoch':[]}

    def on_batch_end(self, batch, logs={}):
        self.losses['batch'].append(logs.get('loss'))
        self.accuracy['batch'].append(logs.get('NMSE'))
        self.val_loss['batch'].append(logs.get('val_loss'))
        self.val_acc['batch'].append(logs.get('val_NMSE'))

    def on_epoch_end(self, batch, logs={}):
        self.losses['epoch'].append(logs.get('loss'))
        self.accuracy['epoch'].append(logs.get('NMSE'))
        self.val_loss['epoch'].append(logs.get('val_loss'))
        self.val_acc['epoch'].append(logs.get('val_NMSE'))

    def loss_plot(self, loss_type):
        iters = range(len(self.losses[loss_type]))
        plt.figure()
        plt.plot(iters, self.losses[loss_type], 'g', label='train loss')
        if loss_type == 'epoch':
            plt.plot(iters, self.val_loss[loss_type], 'k', label='val loss')
        plt.grid(True)
        plt.xlabel(loss_type)
        plt.ylabel('NMSE-loss')
        plt.legend(loc="upper right")
        plt.title('Training and Validation loss')
        plt.savefig('loss curve.png')
        plt.show()

def train_model(model, train, test, epochs, batch_size):
    # X-data
    xx_train = train[0]
    xx_train = xx_train.astype('float32')
    X_train = xx_train

    xx_test = test[0]
    xx_test = xx_test.astype('float32')
    X_test = xx_test

    # Y-data
    yy_train = train[1]
    yy_train = yy_train.astype('float32')
    Y_train = yy_train

    yy_test = test[1]
    yy_test = yy_test.astype('float32')
    Y_test = yy_test

    print('x_train shape:', X_train.shape)
    print('y_train shape:', Y_train.shape)
    print('x_test shape:', X_test.shape)
    print('y_test shape:', Y_test.shape)
    print(X_train.shape[0], 'input samples')
    print(Y_train.shape[0], 'label samples')

    model.compile(loss=NMSE_IRS,
                  optimizer=Adam(0.001),
                  metrics=[NMSE_IRS])

    history = LossHistory()

    callbacks_list = [
        # this callback will stop the training when there is no improvement in the validation loss
        keras.callbacks.EarlyStopping(
            monitor='val_loss', min_delta=0, patience=5, verbose=1, mode='min'),
        # save checkpoint
        keras.callbacks.ModelCheckpoint(
            filepath='Model_210106.h5',  # 文件路径
            verbose=1,
            save_best_only=True,
            mode='min'),
        history]

    start_time = time.clock()

    model.fit(X_train, Y_train,
              batch_size=batch_size,
              epochs=epochs,
              verbose=2,
              validation_data=(X_test, Y_test),
              callbacks=callbacks_list)

    # print training time cost
    end_time = time.clock()
    elapsed_time_train = end_time - start_time
    print(elapsed_time_train)

    # serialize model to JSON
    model_json = model.to_json()
    with open("Model_210106.json", "w") as json_file:
        json_file.write(model_json)
    # serialize weights to HDF5
    model.load_weights('Model_210106.h5')
    model_best = model
    print("Saved best model to disk")

    # print evaluating time cost
    start_time = time.clock()
    score = model_best.evaluate(X_test, Y_test, verbose=1)
    end_time = time.clock()
    elapsed_time = end_time - start_time
    print('Inference time for each sample:', elapsed_time / 20000)

    print('Test best epoch NMSE:', score[1])
    history.loss_plot('epoch')

    # Prediction. the output is saved as .mat files.
    residual_predict = []
    for test_sample in X_test: #X_test shape: (20000, 8, 8, 2)  test_sample shape:(8, 8, 2)
        test_sample = test_sample[np.newaxis, :] #shape: (1,8,8,2)
        y = model_best.predict(test_sample)  #output shape: (1,8,8,2)
        residual = test_sample - y
        residual_ = np.squeeze(residual) #output shape: (8,8,2)
        residual_predict.append(residual_)  #output shape: (20000,8,8,2)

    predict_np = np.array(residual_predict)
    io.savemat('Model_CDRN_20dB_Residual_Prediction210106.mat', {'data': predict_np})
    print('prediction saved to .mat file')

def DnCNN_MultiBlock(block, depth, image_channels, filters=64,  use_bnorm=True):

    input = Input(shape=(None, None, image_channels))
    input_ = input
    x = input_

    B = block

    while B:

        input_ = x
        #For depth-1 layers, Conv+BN+relu
        for i in range(depth-1):

            x = Conv2D(filters=filters, kernel_size=(3, 3), strides=(1, 1), kernel_initializer='Orthogonal', \
                       padding='same', use_bias=False, data_format='channels_last')(input_)

            if use_bnorm:
                x = BatchNormalization(axis=3, momentum=0.0, epsilon=0.0001)(x)

            x = Activation('relu')(x)

        #For last layer, Conv
        x = Conv2D(filters=image_channels, kernel_size=(3, 3), strides=(1, 1), kernel_initializer='Orthogonal', \
                   padding='same', use_bias=False, data_format='channels_last')(x)

        x = Subtract()([input_, x])  # input - noise

        B -= 1

    model = Model(inputs=input, outputs=x)
    return model

def main():

    batch_size = 64
    epochs = 400

    model = DnCNN_MultiBlock(block=3, depth=16, image_channels=2, use_bnorm=True)
    model.summary()

    #load .mat data
    data_xtr = sio.loadmat('x_train_Rician_CSCG_K10dB_60000_M8_N32_5dB.mat')
    data_ytr = sio.loadmat('y_train_Rician_CSCG_K10dB_60000_M8_N32_5dB.mat')
    data_xtest = sio.loadmat('x_test_Rician_CSCG_K10dB_20000_M8_N32_5dB.mat')
    data_ytest = sio.loadmat('y_test_Rician_CSCG_K10dB_20000_M8_N32_5dB.mat')

    #rename data
    xa_train = data_xtr['x_train']
    ya_train = data_ytr['y_train']

    xa_test = data_xtest['x_test']
    ya_test = data_ytest['y_test']

    train_model(model,
                (xa_train, ya_train),
                (xa_test, ya_test),
                epochs,
                batch_size)

if __name__ == '__main__':
    main()

max_stars_repo_name: skydvn/CDRN-channel-estimation-IRS

reason: import numpy,import scipy


text:
from input_data import InputData
import numpy
from model import SkipGramModel
from torch.autograd import Variable
import torch
import torch.nn as nn
import torch.optim as optim
from tqdm import tqdm
import sys

class Word2Vec:
    def __init__(self,
                 input_file_name,
                 output_file_name,
                 emb_dimension=100,
                 batch_size=50,
                 window_size=5,
                 iteration=1,
                 initial_lr=0.025,
                 min_count=5):
        """Initilize class parameters.

        Args:
            input_file_name: Name of a text data from file. Each line is a sentence splited with space.
            output_file_name: Name of the final embedding file.
            emb_dimention: Embedding dimention, typically from 50 to 500.
            batch_size: The count of word pairs for one forward.
            window_size: Max skip length between words.
            iteration: Control the multiple training iterations.
            initial_lr: Initial learning rate.
            min_count: The minimal word frequency, words with lower frequency will be filtered.

        Returns:
            None.
        """
        self.data = InputData(input_file_name, min_count)
        self.output_file_name = output_file_name
        self.emb_size = len(self.data.word2id)
        self.emb_dimension = emb_dimension
        self.batch_size = batch_size
        self.window_size = window_size
        self.iteration = iteration
        self.initial_lr = initial_lr
        self.skip_gram_model = SkipGramModel(self.emb_size, self.emb_dimension)
        self.use_cuda = torch.cuda.is_available()
        if self.use_cuda:
            self.skip_gram_model.cuda()
        self.optimizer = optim.SGD(
            self.skip_gram_model.parameters(), lr=self.initial_lr)

    def train(self):
        """Multiple training.

        Returns:
            None.
        """
        pair_count = self.data.evaluate_pair_count(self.window_size)
        batch_count = self.iteration * pair_count / self.batch_size
        process_bar = tqdm(range(int(batch_count)))
        # self.skip_gram_model.save_embedding(
        #     self.data.id2word, 'begin_embedding.txt', self.use_cuda)
        for i in process_bar:
            pos_pairs = self.data.get_batch_pairs(self.batch_size,
                                                  self.window_size)
            neg_v = self.data.get_neg_v_neg_sampling(pos_pairs, 5)
            pos_u = [pair[0] for pair in pos_pairs]
            pos_v = [pair[1] for pair in pos_pairs]

            pos_u = Variable(torch.LongTensor(pos_u))
            pos_v = Variable(torch.LongTensor(pos_v))
            neg_v = Variable(torch.LongTensor(neg_v))
            if self.use_cuda:
                pos_u = pos_u.cuda()
                pos_v = pos_v.cuda()
                neg_v = neg_v.cuda()

            self.optimizer.zero_grad()
            loss = self.skip_gram_model.forward(pos_u, pos_v, neg_v)
            loss.backward()
            self.optimizer.step()

            process_bar.set_description("Loss: %0.8f, lr: %0.6f" %
                                        (loss,
                                         self.optimizer.param_groups[0]['lr']))
            if i * self.batch_size % 100000 == 0:
                lr = self.initial_lr * (1.0 - 1.0 * i / batch_count)
                for param_group in self.optimizer.param_groups:
                    param_group['lr'] = lr
        self.skip_gram_model.save_embedding(
            self.data.id2word, self.output_file_name, self.use_cuda)

if __name__ == '__main__':
    w2v = Word2Vec(input_file_name=sys.argv[1], output_file_name=sys.argv[2])
    w2v.train()


text:
"""Script to nod between two fixed points in pattern ABBA.

History:
2008-04-23 ROwen
2008-07-24 ROwen    Fixed PR 852: end did not always restore the original boresight.
2011-08-02 ROwen    Added variable offset and test for minimum offset.
2012-11-13 ROwen    Stop using Checkbutton indicatoron=False because it is no longer supported on MacOS X.
"""
import math
import numpy
import Tkinter
import RO.Wdg
import RO.PhysConst
import TUI.TCC.TCCModel
import TUI.Inst.ExposeModel
import TUI.Inst.ExposeStatusWdg
import TUI.Inst.ExposeInputWdg

# constants
InstName = "TSpec"
# Offset from A to B in x, y pixels
# (85, 1) is from Matt Nelson 2008-04
ABOffsetPix = (85, 1)
MaxCycles = 9999
 # Instrument scale in unbinned pixels/degree on the sky
# as measured by APO 2008-03-21 but averaged x and y
InstScale = (-14352, 14342)
ABOffsetArcSec = tuple(ABOffsetPix[i] * 3600.0 / InstScale[i] for i in range(2))
MinOffsetArcSec = 0.1 # minimum offset in arcsec
HelpURL = "Scripts/BuiltInScripts/TSpecNod.html"

class ScriptClass(object):
    def __init__(self, sr):
        """Set up widgets to set input exposure time,
        drift amount and drift speed.
        """
        # if True, run in debug-only mode (which doesn't DO anything, it just pretends)
        sr.debug = False
        self.sr = sr

        self.begOffset = numpy.array((numpy.nan, numpy.nan))
        self.currOffset = self.begOffset[:]

        self.tccModel = TUI.TCC.TCCModel.getModel()
        self.expModel = TUI.Inst.ExposeModel.getModel(InstName)

        row=0

        # standard exposure status widget
        expStatusWdg = TUI.Inst.ExposeStatusWdg.ExposeStatusWdg(
            master = sr.master,
            instName = InstName,
            helpURL = HelpURL,
        )
        expStatusWdg.grid(row=row, column=0, sticky="news")
        row += 1

        # separator
        Tkinter.Frame(sr.master,
            bg = "black",
        ).grid(row=row, column=0, pady=2, sticky="ew")
        row += 1

        # standard exposure input widget
        self.expWdg = TUI.Inst.ExposeInputWdg.ExposeInputWdg(
            master = sr.master,
            instName = InstName,
            expTypes = "object",
            helpURL = HelpURL,
        )
        self.expWdg.numExpWdg.helpText = "# of exposures at each point"
        self.expWdg.grid(row=row, column=0, sticky="news")
        row += 1

        # add some controls to the exposure input widget

        self.offsetEnableWdg = RO.Wdg.Checkbutton (
            master = self.expWdg,
            text = "Offset",
            defValue = False,
            helpText = "Enable/disable manual control of offset",
            helpURL = HelpURL,
        )
        self.offsetArcSecWdgSet = []
        offsetFrame = Tkinter.Frame(self.expWdg)
        for i, axis in enumerate(("x", "y")):
            offsetWdg = RO.Wdg.FloatEntry(
                master = offsetFrame,
                callFunc = self._doOffsetWdg,
                helpText = "Amount of %s offset in arcsec" % (axis,),
                helpURL = HelpURL,
            )
            offsetWdg.pack(side="left")
            self.offsetArcSecWdgSet.append(offsetWdg)
        RO.Wdg.Label(master=offsetFrame, text="arcsec").pack(side="left")
        self.offsetPixWdg = RO.Wdg.Label(
            master = offsetFrame,
            helpText = "Amount of x,y offset in pixels",
            helpURL = HelpURL,
        )
        self.offsetPixWdg.pack(side="left")
        self.expWdg.gridder.gridWdg(self.offsetEnableWdg, offsetFrame, colSpan=5)

        # number of cycles
        self.numCyclesWdg = RO.Wdg.IntEntry (
            master = self.expWdg,
            minValue = 1,
            maxValue = MaxCycles,
#            width = 10,
            helpText = "Number of ABBA cycles",
            helpURL = HelpURL,
        )
        self.expWdg.gridder.gridWdg("Cycles", self.numCyclesWdg)

        self.expWdg.gridder.allGridded()

        if sr.debug:
            # set useful debug defaults
            self.expWdg.timeWdg.set("1.0")
            self.expWdg.numExpWdg.set(2)
            self.expWdg.fileNameWdg.set("debug")
            self.numCyclesWdg.set(2)

        self.offsetEnableWdg.addCallback(self._doOffsetEnable, callNow=True)

    def _doOffsetEnable(self, wdg=None):
        """Handle toggle of offsetEnableWdg
        """
        if self.offsetEnableWdg.getBool():
            for wdg in self.offsetArcSecWdgSet:
                wdg.setEnable(True)
        else:
            for i, wdg in enumerate(self.offsetArcSecWdgSet):
                wdg.set(ABOffsetArcSec[i])
                wdg.setEnable(False)

    def _doOffsetWdg(self, wdg=None):
        """Handle changes in offset in arcsec by updating offset in pixels
        """
        if self.isOffsetOK():
            offsetArcSec = self.getOffsetArcSec()
            offsetPix = tuple(offsetArcSec[i] * InstScale[i] / 3600.0 for i in range(2))
            self.offsetPixWdg.set("= %0.1f, %0.1f pix" % offsetPix, severity=RO.Constants.sevNormal)
        else:
            self.offsetPixWdg.set("offset too small", severity=RO.Constants.sevError)

    def getOffsetArcSec(self):
        """Return current x,y offset in arcsec; 0 if field is empty
        """
        return [wdg.getNum() for wdg in self.offsetArcSecWdgSet]

    def isOffsetOK(self):
        """Is offset acceptably large? The test right now is 0
        """
        return math.hypot(*self.getOffsetArcSec()) >= MinOffsetArcSec

    def end(self, sr):
        """If telescope offset, restore original position.
        """
        # restore original boresight position, if changed
        if self.needMove(self.begOffset):
            tccCmdStr = "offset boresight %.7f, %.7f/pabs/vabs/computed" % tuple(self.begOffset)
            #print "sending tcc command %r" % tccCmdStr
            sr.startCmd(
                actor = "tcc",
                cmdStr = tccCmdStr,
            )

    def needMove(self, desOffset):
        """Return True if telescope not at desired offset"""
        if numpy.any(numpy.isnan(self.currOffset)):
            return False
        return not numpy.allclose(self.currOffset, desOffset)

    def run(self, sr):
        """Take one or more exposures while moving the object
        in the +X direction along the slit.
        """
        # make sure the current instrument matches the desired instrument
        if not sr.debug:
            currInst = sr.getKeyVar(self.tccModel.instName)
            if InstName.lower() != currInst.lower():
                raise sr.ScriptError("%s is not the current instrument!" % InstName)

        # record the current boresight position
        begBorePVTs = sr.getKeyVar(self.tccModel.boresight, ind=None)
        if not sr.debug:
            begOffset = [pvt.getPos() for pvt in begBorePVTs]
            if None in begOffset:
                raise sr.ScriptError("Current boresight position unknown")
            self.begOffset = numpy.array(begOffset, dtype=float)
        else:
            self.begOffset = numpy.zeros(2, dtype=float)
        self.currOffset = self.begOffset[:]
        #print "self.begOffset=%r" % self.begOffset

        numCycles = self.numCyclesWdg.getNum()
        if not numCycles:
            raise sr.ScriptError("Must specify number of cycles")

        # exposure command without startNum and totNumExp
        # get it now so that it will not change if the user messes
        # with the controls while the script is running
        numExpPerNode = self.expWdg.numExpWdg.getNumOrNone()
        if numExpPerNode is None:
            raise sr.ScriptError("must specify #Exp")

        nodeNames = ("A1", "B1", "B2", "A2")
        numNodes = len(nodeNames)
        totNumExp = numCycles * numNodes * numExpPerNode
        expCmdPrefix = self.expWdg.getString(totNum = totNumExp)
        if expCmdPrefix is None:
            raise sr.ScriptError("missing inputs")

        NodeOffsetDict = dict (
            A = numpy.zeros(2, dtype=float),
            B = numpy.array(self.getOffsetArcSec(), dtype=float) / 3600.0,
        )

        if not self.isOffsetOK():
            raise sr.ScriptError("offset too small")

        numExpTaken = 0
        for cycle in range(numCycles):
            for nodeName in nodeNames:
                iterName = "Cycle %s of %s, Pos %s" % (cycle + 1, numCycles, nodeName)
                nodeOffset = NodeOffsetDict[nodeName[0]]
                desOffset = self.begOffset + nodeOffset
                if self.needMove(desOffset):
                    sr.showMsg("%s: Offsetting" % (iterName,))
                    yield self.waitOffset(desOffset)

                # expose
                sr.showMsg("%s: Exposing" % (iterName,))
                startNum = numExpTaken + 1
                expCmdStr = "%s startNum=%s" % (expCmdPrefix, startNum)
                #print "sending %s command %r" % (InstName, expCmdStr)
                yield sr.waitCmd(
                    actor = self.expModel.actor,
                    cmdStr = expCmdStr,
                    abortCmdStr = "abort",
                )
                numExpTaken += numExpPerNode

    def waitOffset(self, desOffset):
        """Offset the telescope"""
        tccCmdStr = "offset boresight %.7f, %.7f/pabs/vabs/computed" % tuple(desOffset)
        self.currOffset = desOffset[:]
        yield self.sr.waitCmd(
            actor = "tcc",
            cmdStr = tccCmdStr,
        )

max_stars_repo_name: ApachePointObservatory/TUI

max_issues_repo_name: ApachePointObservatory/TUI

max_forks_repo_name: r-owen/TUI


text:
#%% Imports
import numpy as np
import os
from sys import stdout
import time
import json
import pandas as pd
import random
import pickle
import matplotlib.pyplot as plt
from tqdm import tqdm
import glob

#import seaborn as sns
import umap

import tensorflow as tf
import tensorflow_probability as tfp

#import betavae as cv
#import pcvae as cv
import kcvae as kcv

import utils as ut
import logger
import configs as cf

import importlib

AUTOTUNE = tf.data.experimental.AUTOTUNE

JUPYTER_NOTEBOOK = True

# if JUPYTER_NOTEBOOK:
# %reload_ext autoreload
# %autoreload 2

#%% Setup generic values...
#######
cf_img_size = cf.IMG_SIZE
cf_latent_dim = cf.LATENT_DIM
cf_batch_size = cf.BATCH_SIZE #32
cf_learning_rate = cf.IMGRUN_LR #4e-4
cf_limits = [cf_img_size, cf_img_size]
cf_pixel_dim = cf_img_size*cf_img_size*3
#( *-*) ( *-*)>⌐■-■ ( ⌐■-■)
#
cf_kl_weight = cf.KL_WEIGHT
cf_beta = cf_kl_weight
cf_num_epochs = cf.N_IMGRUN_EPOCH
#dfmeta = ut.read_meta()
cf_val_frac = cf.VALIDATION_FRAC
#%%  are we GPU-ed?
tf.config.experimental.list_physical_devices('GPU')
#%% helpers

# # %%
# #%%
# sns.pairplot(digits_df, hue='digit', palette='Spectral')

# reducer = umap.UMAP(random_state=42)
# reducer.fit(digits_df.to_numpy())

# #%% Run model on all data to get latent vects and loss. Used for streamlit app and other places.
# #preds,losses = ut.dumpReconstruct( model, train_dataset, test_dataset )
# ds = ut.load_and_dump(pix_dim, cf.IMAGE_FILEPATH)
# #or _samples, _labels in ds.take(1) : pass
# # remake this to simply go through all the data and calculate the embedding and loss... new functions probably...
# %%
# HELPERS:
def get_n_samples(ds):
    n_samples = 0
    for _ in ds :
        n_samples += 1
    return n_samples

def get_ds(pix_dim):
    return ut.load_and_dump(pix_dim, cf.IMAGE_FILEPATH)

def create_imgStats(ds,n_samples):
    db = {}
    for sample, label in tqdm(ds,
                                unit_scale=True,
                                desc="Saving data stats: ",
                                unit=" samples",
                                total=n_samples ) :

        key = label.numpy().decode().split('/')[-3]
        if key in db.keys():
            imgs = db[key]
            db[key] = np.concatenate((imgs,sample[None,...]),axis=0)
        else:
            db[key] = sample[None,...]

    return db

def create_snk2dicts(ds,n_samples,vae):
    snk2loss = {}
    snk2vec = {}
    for sample, label in tqdm(ds,
                                unit_scale=True,
                                desc="Saving shape 2 vec: ",
                                unit=" encodes",
                                total=n_samples ) :
        key = label.numpy()  # maybe should have changed this to a (string... but byte is good...
        snk2vec[key] = vae.encoder(sample[None,...]).sample().numpy()[0] #assume that we don't have a batch
        snk2loss[key] = vae.partial_vae_loss(sample[None,...]).numpy()

    return snk2loss, snk2vec

def load_and_prime_model(model,model_name,params,data_dir,epochs):
    filename = f"overtrain-{model_name}-kl_weight{params['kl_weight']:03d}.pkl"
    hist,params = ut.load_pickle(os.path.join(data_dir,filename))

    vae = model(dim_z=params['z_dim'], dim_x=params['x_dim'],
                learning_rate=0.0001, kl_weight=params['kl_weight'])

    vae.compile(optimizer=vae.optimizer, loss = vae.partial_vae_loss)
    sv_path = os.path.join(data_dir,f"kl_weight{params['kl_weight']:03d}")
    vae.load_model(sv_path, epochs)

    return vae,hist

def dump_pickle(filepath, item_to_save):
    f = open(filepath, "wb")
    pickle.dump(item_to_save, f)
    f.close()

def load_pickle(filepath):
    infile = open(filepath, "rb")
    item = pickle.load(infile)
    infile.close()
    return item

def load_snk2pickles(data_dir,kl_weight):
    load_dir = os.path.join(data_dir,f"kl_weight{kl_weight:03d}")
    snk2vec = load_pickle(os.path.join(load_dir,"snk2vec.pkl"))
    snk2loss = load_pickle(os.path.join(load_dir,"snk2loss.pkl"))
    return snk2loss, snk2vec

def dump_snk2pickles(snk2loss,snk2vec,data_dir,kl_weight):
    load_dir = os.path.join(data_dir,f"kl_weight{kl_weight:03d}")
    dump_pickle(os.path.join(load_dir,"snk2vec.pkl"), snk2vec)
    dump_pickle(os.path.join(load_dir,"snk2loss.pkl"), snk2loss)

def load_snk2umap(data_dir,kl_weight):
    load_dir = os.path.join(data_dir,f"kl_weight{kl_weight:03d}")
    umap = load_pickle(os.path.join(load_dir,"snk2umap.pkl"))
    return snk2loss, snk2vec

def dump_snk2umap(snk2umap,data_dir,kl_weight):
    load_dir = os.path.join(data_dir,f"kl_weight{kl_weight:03d}")
    print(load_dir)
    dump_pickle(os.path.join(load_dir,"snk2umap.pkl"), snk2umap)

def get_umap_embedding(snk2vec):
    latents = np.array(list(snk2vec.values()))
    reducer = umap.UMAP(random_state=42)
    reducer.fit(latents)
    embedding = reducer.transform(latents)
    assert(np.all(embedding == reducer.embedding_))
    snk2umap = dict(zip(np.array(list(snk2vec.keys())),embedding))

    return snk2umap

#%% dump the vectors to a dictionary

# ds = get_ds(pix_dim)
# n_samples = get_n_samples(ds)

# model = kcv.K_PCVAE
# model_name = "K_PCVAE"
# epochs = 400

# data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"
# vae,hist = load_and_prime_model(model,model_name,params,data_dir,epochs)

# make_dicts = False
# if make_dicts:
#     snk2loss,snk2vec,svecs = create_snk2dicts(ds,n_samples,vae, dump_dicts=True)
# else:
#     snk2loss,snk2vec = load_snk2pickles(data_dir)

# %%  define function to create snk2vec, snk2loss, and UMAP objects for each network / latents

# %%
#############################
#
#  comparison 1:  KL_Reg model vs KL_loss model @x160
#
###############################
latent_dim = 40
pix_dim = 160
batch_size = 64
kl_weight = 5

p_names = ['z_dim','x_dim','kl_weight','batch_size']
p_vals = [latent_dim, (pix_dim,pix_dim,3), kl_weight, batch_size]
params = dict(zip(p_names,p_vals))

ds = get_ds(pix_dim)
n_samples = get_n_samples(ds)

#loop over models
models = [kcv.K_PCVAE, kcv.K_PCVAE_KL_Reg]
model_names = ["K_PCVAE", "K_PCVAE_KL_Reg"]
data_dirs = [f"{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}" for model_name in model_names]

klws = [1,2,3,5,10]
latents = [32,40,64]
epochs = 400

p_names = ['z_dim','x_dim','kl_weight','batch_size']
p_vals = [latent_dim, (pix_dim,pix_dim,3), kl_weight, batch_size]
params = dict(zip(p_names,p_vals))

for model,model_name in zip(models,model_names):

    for kl in klws:
        params['kl_weight'] = kl
        for l in latents:

            params['z_dim'] = l
            data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"
            print(data_dir)

            vae, hist = load_and_prime_model(model,model_name,params,data_dir,epochs)
            snk2loss,snk2vec = create_snk2dicts(ds,n_samples,vae)

            dump_snk2pickles(snk2loss,snk2vec,data_dir,params['kl_weight'])

#

#############################
#
#  comparison 2:  KL_loss model across X_di, = @x124, @x160, @x192, @x224
#
###############################

model = models[0]
model_name = model_names[0]

# model = kcv.K_PCVAE
# model_name = "K_PCVAE"
klws = [1,2,3,5,10]
latent_dims = [32,40,64]
pix_dims = [128,160,192]
epochs = 400

# finished everything but the biggest KL_WEIHT

# loop over pix_dim
for pix_dim in pix_dims:
    ds = get_ds(pix_dim)
    n_samples = get_n_samples(ds)
    params['x_dim']=(pix_dim,pix_dim,3)
    for kl in klws:
        params['kl_weight'] = kl
        for l in latent_dims:

            params['z_dim'] = l
            data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"
            print(data_dir)

            vae,hist = load_and_prime_model(model,model_name,params,data_dir,epochs)
            snk2loss,snk2vec = create_snk2dicts(ds,n_samples,vae)

            dump_snk2pickles(snk2loss,snk2vec,data_dir,params['kl_weight'])

#%%
ds = get_ds(pix_dim)
n_samples = get_n_samples(ds)

#%%
latent_dim = 64
pix_dim = 160
batch_size = 64
kl_weight = 5

# model = kcv.K_PCVAE
# model_name = "K_PCVAE"
klws = [1,2,3,5,10]
latent_dims = [32,40,64]
pix_dims = [128,160,192]
epochs = 400

p_names = ['z_dim','x_dim','kl_weight','batch_size']
p_vals = [latent_dim, (pix_dim,pix_dim,3), kl_weight, batch_size]
params = dict(zip(p_names,p_vals))

#loop over models
models = [kcv.K_PCVAE, kcv.K_PCVAE_KL_Reg]
model_names = ["K_PCVAE", "K_PCVAE_KL_Reg"]
model = models[1]
model_name = model_names[0]

pix_dim = pix_dims[2]
params['x_dim']=(pix_dim,pix_dim,3)
kl = klws[4]
params['kl_weight'] = kl
l = latent_dims[2]
params['z_dim'] = l
data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"
print(data_dir)
#vae,hist = load_and_prime_model(model,model_name,params,data_dir,epochs)

#%%
snk2loss, snk2vec = load_snk2pickles(data_dir,params['kl_weight'])

desc_df = pd.read_pickle(f"{cf.DESCRIPTIONS}.pkl")
meta_df = pd.read_pickle(f"{cf.META_DATA}.pkl")

# fix the root directory....
#%%
#############################
#
#  make comparison... for each network compare the UMAP...
#
###############################
latents = np.array(list(snk2vec.values()))
losses = np.array(list(snk2loss.values()))
labels = np.array(list(snk2vec.keys()))
reducer = umap.UMAP(random_state=42)
reducer.fit(latents)

embedding = reducer.transform(latents)

assert(np.all(embedding == reducer.embedding_))
embedding.shape

plt.scatter(embedding[:, 0], embedding[:, 1], cmap='Spectral', s=5)
#plt.scatter(embedding[:, 0], embedding[:, 1], c=digits.target, cmap='Spectral', s=5)
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('UMAP projection of VAE ', fontsize=24)
plt.savefig("media/umap_small.png")

latents = np.array(list(snk2vec.values()))
losses = np.array(list(snk2loss.values()))
labels = np.array(list(snk2vec.keys()))
reducer = umap.UMAP(random_state=42)
reducer.fit(latents)

embedding = reducer.transform(latents)

assert(np.all(embedding == reducer.embedding_))
embedding.shape

# %%
#############################
#
#  comparison 1:  KL_Reg model vs KL_loss model @x160
#
###############################
latent_dim = 40
pix_dim = 160
batch_size = 64
kl_weight = 5

p_names = ['z_dim','x_dim','kl_weight','batch_size']
p_vals = [latent_dim, (pix_dim,pix_dim,3), kl_weight, batch_size]
params = dict(zip(p_names,p_vals))

ds = get_ds(pix_dim)
n_samples = get_n_samples(ds)

#loop over models
models = [kcv.K_PCVAE, kcv.K_PCVAE_KL_Reg]
model_names = ["K_PCVAE", "K_PCVAE_KL_Reg"]
data_dirs = [f"{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}" for model_name in model_names]

klws = [1,2,3,5,10]
latents = [32,40,64]
epochs = 400

p_names = ['z_dim','x_dim','kl_weight','batch_size']
p_vals = [latent_dim, (pix_dim,pix_dim,3), kl_weight, batch_size]
params = dict(zip(p_names,p_vals))

for model,model_name in zip(models,model_names):
    for kl in klws:
        params['kl_weight'] = kl
        for l in latents:
            params['z_dim'] = l
            data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"
            #print(data_dir)
            # vae, hist = load_and_prime_model(model,model_name,params,data_dir,epochs)
            # snk2loss,snk2vec = create_snk2dicts(ds,n_samples,vae)
            # dump_snk2pickles(snk2loss,snk2vec,data_dir,params['kl_weight'])
            snk2loss,snk2vec = load_snk2pickles(data_dir,params['kl_weight'])
            snk2umap = get_umap_embedding(snk2vec)
            dump_snk2umap(snk2umap,data_dir,kl)
#

#############################
#
#  comparison 2:  KL_loss model across X_di, = @x124, @x160, @x192, @x224
#
###############################

model = models[0]
model_name = model_names[0]
klws = [1,2,3,5,10]
latent_dims = [32,40,64]
pix_dims = [128,160,192]
epochs = 400

# loop over pix_dim
for pix_dim in pix_dims:
    ds = get_ds(pix_dim)
    n_samples = get_n_samples(ds)
    params['x_dim']=(pix_dim,pix_dim,3)
    for kl in klws:
        params['kl_weight'] = kl
        for l in latent_dims:
            params['z_dim'] = l
            data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"
            print(data_dir)
            # vae, hist = load_and_prime_model(model,model_name,params,data_dir,epochs)
            # snk2loss,snk2vec = create_snk2dicts(ds,n_samples,vae)
            # dump_snk2pickles(snk2loss,snk2vec,data_dir,params['kl_weight'])
            snk2loss,snk2vec = load_snk2pickles(data_dir,params['kl_weight'])
            snk2umap = get_umap_embedding(snk2vec)
            dump_snk2umap(snk2umap,data_dir,kl)

#%%

# %%
from io import BytesIO
from PIL import Image
import base64
def embeddable_image(label):
    # img_data = 255 - 15 * data.astype(np.uint8)
    # image = Image.fromarray(img_data, mode='L').resize((64, 64), Image.BICUBIC)
    # buffer = BytesIO()
    # image.save(buffer, format='jpg')
    # for_encoding = buffer.getvalue()
    # return 'data:image/png;base64,' + base64.b64encode(for_encoding).decode()

    # loads the image if string input

    #image = get_thumbnail(label)
    return image_formatter(label)
    # return 'data:image/jpeg;base64,' + base64.b64encode(for_encoding).decode()

def get_thumbnail(path):
    i = Image.open(path)
    i.thumbnail((64, 64), Image.LANCZOS)
    return i

def image_base64(im):
    if isinstance(im, str):
        im = get_thumbnail(im)
    with BytesIO() as buffer:
        im.save(buffer, 'png')
        return base64.b64encode(buffer.getvalue()).decode()

def image_formatter(im):
    return f"data:image/png;base64,{image_base64(im)}"
    #return f''

#%%

from bokeh.plotting import figure, show, output_notebook
from bokeh.models import HoverTool, ColumnDataSource, CategoricalColorMapper
from bokeh.palettes import Spectral10

output_notebook()

digits_df = pd.DataFrame(embedding, columns=('x', 'y'))
digits_df['digit'] = [str(x.decode()) for x in labels]
digits_df['image'] = digits_df.digit.map(lambda f: embeddable_image(f))
digits_df['fname'] = digits_df.digit.map(lambda x: f"{x.split('/')[-3]} {x.split('/')[-1]}")
digits_df['db'] = digits_df.digit.map(lambda x: f"{x.split('/')[-3]}")
digits_df['loss'] = [f"{x:.1f}" for x in losses]

datasource = ColumnDataSource(digits_df)
color_mapping = CategoricalColorMapper(factors=["sns","goat"],
                                       palette=Spectral10)

plot_figure = figure(
    title='UMAP projection VAE latent',
    plot_width=1000,
    plot_height=1000,
    tools=('pan, wheel_zoom, reset')
)

plot_figure.add_tools(HoverTool(tooltips="""

        @fname
        @loss

"""))

plot_figure.circle(
    'x',
    'y',
    source=datasource,
    color=dict(field='db', transform=color_mapping),
    line_alpha=0.6,
    fill_alpha=0.6,
    size=4
)

show(plot_figure)

#%%%
vecs = np.array(list(snk2vec.values()))

labs = np.array(list(snk2vec.keys()))

reducer = umap.UMAP(random_state = 488)

embedding = reducer.transform(vecs)

reducer.fit(vecs)

plt.scatter(embedding[:, 0], embedding[:, 1], cmap='Spectral', s=5)
#plt.scatter(embedding[:, 0], embedding[:, 1], c=digits.target, cmap='Spectral', s=5)
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('UMAP projection of the Digits dataset', fontsize=24)

#%%

snk2loss, snk2vec = load_snk2pickles(data_dir)

dfdesc = pd.read_pickle(f"{cf.DESCRIPTIONS}.pkl")
df_meta = pd.read_pickle(f"{cf.META_DATA}.pkl")
# fix the root directory....

digits_df = pd.DataFrame.from_dict(snk2vec, orient='index')

#

pix_dim = 160
batch_size = 64
epochs = 400
kl_weight = 5
latent_dim = 40

models = [kcv.K_PCVAE, kcv.K_PCVAE_KL_Reg]
model_names = ["K_PCVAE", "K_PCVAE_KL_Reg"]

data_dirs = [f"{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}" for model_name in model_names]

params['kl_weight'] = kl_weight
params['z_dim'] = latent_dim

print(f"training beta={kl_weight} z={latent_dim}")
data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"
print(data_dir)

model = kcv.K_PCVAE
model_name = "K_PCVAE"
data_dir = f"data/{model_name}-X{params['x_dim'][0]}-Z{params['z_dim']}"

filename = f"overtrain-{model_name}-kl_weight{params['kl_weight']:03d}.pkl"
history,params = ut.load_pickle(os.path.join(data_dir,filename))

vae = model(dim_z=params['z_dim'], dim_x=params['x_dim'],
                learning_rate=0.0001, kl_weight=params['kl_weight'])

vae.compile(optimizer=vae.optimizer, loss = vae.partial_vae_loss)
sv_path = os.path.join(data_dir,f"kl_weight{params['kl_weight']:03d}")
vae.load_model(sv_path, epochs)

# # %%
# x
# # %%
# x.shape
# # %%
# xhat = vae(x)
# # %%
# tf.math.squared_difference(x,xhat).numpy().mean()
# # %%
# xhat.numpy().min()
# # %%
# plt.imshow(xhat[0,].numpy().squeeze())
# # %%

#%%

df_meta.head()

# UMAP(a=None, angular_rp_forest=False, b=None,
#      force_approximation_algorithm=False, init='spectral', learning_rate=1.0,
#      local_connectivity=1.0, low_memory=False, metric='euclidean',
#      metric_kwds=None, min_dist=0.1, n_components=2, n_epochs=None,
#      n_neighbors=15, negative_sample_rate=5, output_metric='euclidean',
#      output_metric_kwds=None, random_state=42, repulsion_strength=1.0,
#      set_op_mix_ratio=1.0, spread=1.0, target_metric='categorical',
#      target_metric_kwds=None, target_n_neighbors=-1, target_weight=0.5,
#      transform_queue_size=4.0, transform_seed=42, unique=False, verbose=False)
            # convert to arrays

latents = np.array(list(snk2vec.values()))
losses = np.array(list(snk2loss.values()))
labels = np.array(list(snk2vec.keys()))
reducer = umap.UMAP(random_state=42)
reducer.fit(latents)

embedding = reducer.transform(digits_df.to_numpy())
# Verify that the result of calling transform is
# idenitical to accessing the embedding_ attribute
assert(np.all(embedding == reducer.embedding_))
embedding.shape

plt.scatter(embedding[:, 0], embedding[:, 1], cmap='Spectral', s=5)
#plt.scatter(embedding[:, 0], embedding[:, 1], c=digits.target, cmap='Spectral', s=5)
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('UMAP projection of the Digits dataset', fontsize=24)

#%%
# Importing Image module from PIL package
from PIL import Image
import PIL

    # if tf.io.is_jpeg(filename):
    #     image_string = tf.io.read_file(filename)
    #     im = tf.image.decode_jpeg(image_string, channels=3)
    #     im = tf.cast(im, tf.float32) / 255.0
    # else:
    #     #print("not a jpeg")
    #     image_string = tf.io.read_file(filename)
    #     im = tf.image.decode_png(image_string, channels=0)  #
    #     mask = tf.cast(tf.expand_dims(im[:,:,3],2),tf.float32) / 255.
    #     mask2 = tf.concat((mask,mask,mask),axis=2)
    #     im = tf.cast(im, tf.float32) / 255.0
    #     im = (1.-mask2) + mask2*im[:,:,0:3]

    # label = filename
    # return im, label

# %%
# this is to fix the PNGs which convert to black backgrounds...
# imoved all the files into "im" and then copied back to the root dir...
files_n = glob.glob(os.path.join(cf.IMAGE_FILEPATH, "goat/im/*"))
out_root = '/home/ergonyc/Projects/DATABASE/SnkrScrpr/data/goat/img'
skip_list = []
its = 0
for f in files:
    #im = load_and_convert_to_dump(f)
    im = Image.open(f)
    its += 1
    outf = os.path.join(out_root,f.split('/')[-1])
    if im.format != 'JPEG':
        print(f"not JPEG>>>>  {im.format}")
        skip_list.append(f)
        #im.thumbnail(im.size, Image.LANCZOS)
        bg = Image.new("RGB", im.size, (255,255,255))
        bg.paste(im,box=(0,0),mask=im.convert('RGBA').split()[-1])
        im = bg

    im.save(outf, "JPEG")

    # PIL.Image.fromarray(im.numpy()).save(outf)

# save a image using extension
ut.dump_pickle("goat_png_fix.pkl",skip_list)

# %%
# this is to fix big border on the edges compared to the goat database
# imoved all the files into "im" and then copied back to the root dir...
files = glob.glob(os.path.join(cf.IMAGE_FILEPATH, "sns/im/*"))
out_root = '/home/ergonyc/Projects/DATABASE/SnkrScrpr/data/sns/img'

skip_list = []
its = 0
for f in files:
    #im = load_and_convert_to_dump(f)
    im = Image.open(f)
    its += 1
    outf = os.path.join(out_root,f.split('/')[-1])
    if im.format != 'JPEG':
        print(f"not JPEG>>>>  {im.format}")
        skip_list.append(f)

        # #im.thumbnail(im.size, Image.LANCZOS)
        # bg = Image.new("RGB", im.size, (255,255,255))
        # bg.paste(im,box=(0,0),mask=im.convert('RGBA').split()[-1])
        # im = bg
    # lets take 10% pixels off
    dp = 100//12
    wd, ht = im.size
    dy, dx = wd//dp, ht//dp
    crop_bx = (dx, dx, wd-dx, ht-dx)
    im_crop = im.crop(crop_bx)
    im_crop.save(outf, "JPEG")

    # PIL.Image.fromarray(im.numpy()).save(outf)

# save a image using extension

# %%

ds = get_ds(pix_dim)
n_samples = get_n_samples(ds)

db = create_imgStats(ds,n_samples)

sn = db['sns']
gt = db['goat']
plt.imshow(sn.std(axis=0).squeeze())
plt.colorbar()

plt.imshow(gt.std(axis=0).squeeze())
plt.colorbar()

# %%

max_stars_repo_name: ergonyc/SneakerGen

max_issues_repo_name: ergonyc/SneakerGen

max_forks_repo_name: ergonyc/SneakerGen


text:
# Problem: https://www.hackerrank.com/challenges/s10-basic-statistics/problem
# Score: 30

import numpy as np
from scipy import stats

n = int(input())
arr = list(map(int, input().split()))
print(np.mean(arr))
print(np.median(arr))
print(stats.mode(arr)[0][0])

max_stars_repo_path: Tutorials/01. 10 Days of Statistics/001. Day 0 - Mean, Median, and Mode.py


text:
#Context module. Use a pretrain model to extract the penultimate layer of the model for surrounding trees.
import tensorflow as tf
import rasterio
import cv2
import numpy as np
import pandas as pd

from DeepTreeAttention.utils.paths import find_sensor_path, elevation_from_tile
from DeepTreeAttention.utils.image import image_normalize, resize, crop_image
from sklearn.neighbors import BallTree

def get_nearest(src_points, candidates, k_neighbors=1, distance_threshold=None):
    """Find nearest neighbors for all source points from a set of candidate points
    Args:
    src_points: an pandas row with a geometry column
    candidates: pandas df
    k_neighbors: number of neighbors
    distance_threshold = minimum distance in meters
    """

    # Create tree from the candidate points
    coordinates = np.vstack(candidates.geometry.centroid.apply(lambda geom: (geom.x,geom.y)))
    tree = BallTree(coordinates, leaf_size=15, metric='haversine')

    # Find closest points and distances
    #src_points = src_points.reset_index()
    src_x = src_points.geometry.centroid.x
    src_y = src_points.geometry.centroid.y

    src_points = np.array([src_x, src_y]).reshape(-1,2)

    #If there are not enough neighbors, reduce K, then pad to original
    if k_neighbors > candidates.shape[0]:
        effective_neighbors = candidates.shape[0]
    else:
        effective_neighbors = k_neighbors

    distances, indices = tree.query(src_points, k=effective_neighbors)

    neighbor_geoms = candidates[candidates.index.isin(indices[0])]
    neighbor_geoms = neighbor_geoms.loc[indices[0]]

    #order by the indices
    neighbor_geoms["distance"] = distances[0]

    if distance_threshold:
        neighbor_geoms = neighbor_geoms[neighbor_geoms.distance > distance_threshold]

    # Return indices and distances
    return neighbor_geoms

def predict_neighbors(target, HSI_size, neighbor_pool, metadata, raster, model, k_neighbors=5):
    """Get features of surrounding n trees
    Args:
        target: geometry object of the target point
    neighbor_pool: geopandas dataframe with points
    metadata: The metadata layer for each of the points, assumed to be identical for all neighbors
    n: Number of neighbors
    model: A model object to predict features
    Returns:
    n * m feature matrix, where n is number of neighbors and m is length of the penultimate model layer
    """

    #Find neighbors
    neighbor_geoms = get_nearest(target, candidates = neighbor_pool , k_neighbors=k_neighbors)

    #Always put itself first
    neighbor_geoms = neighbor_geoms.sort_values("distance")

    #extract crop for each neighbor
    features = [ ]
    distances = [ ]
    for index, row in neighbor_geoms.iterrows():
        try:
            crop = crop_image(src=raster, box=row["geometry"])
        except:
            continue

        #reorder to channels last
        crop = resize(crop, HSI_size, HSI_size)
        crop = image_normalize(crop)
        crop = np.expand_dims(crop, 0)

        #create batch
        elevation = np.expand_dims(metadata[0],axis=0)
        site = np.expand_dims(metadata[1],axis=0)
        domain = np.expand_dims(metadata[2],axis=0)
        batch  = [crop,elevation,site,domain]

        feature = model.predict(batch)
        features.append(feature)
        distances.append(row["distance"])

    #if there are fewer than k_neighbors, pad with 0's and large distances (?)
    if len(features) < k_neighbors:
        for x in np.arange(k_neighbors - len(features)):
            features.append(np.zeros(feature.shape))
            distances.append(9999)

    features = np.vstack(features)

    return features, distances

def extract_features(df, x, model_class, hyperspectral_pool, site_label_dict, domain_label_dict, HSI_size=20, k_neighbors=5):
    """Generate features
    Args:
    df: a geopandas dataframe
    x: individual id to use a target
    model_class: A deeptreeattention model class to extract layer features
    hyperspectral_pool: glob dir to search for sensor files
    HSI_size: size of HSI crop
    site_label_dict: dictionary of numeric site labels
    domain_label_dict: dictionary of numeric domain labels
    k_neighbors: number of neighbors to extract
    Returns:
    feature_array: a feature matrix of encoded bottleneck layer
    """
    #Due to resampling, there will be multiple rows of the same point, all are identical.
    #Always pick itself as neighbor 1
    target  =  df[df.individual == x].head(1)
    target = target.reset_index(drop=True)
    sensor_path = find_sensor_path(bounds=target.total_bounds, lookup_pool=hyperspectral_pool)

    #Encode metadata
    site = target.siteID.values[0]
    numeric_site = site_label_dict[site]
    one_hot_sites = tf.one_hot(numeric_site, model_class.sites)

    domain = target.domainID.values[0]
    numeric_domain = domain_label_dict[domain]
    one_hot_domains = tf.one_hot(numeric_domain, model_class.domains)

    #for tests, dummy elevation variable
    try:
        elevation = elevation_from_tile(sensor_path)/1000
    except:
        print("Dummy variable for elevation debug")
        elevation = 100/1000

    metadata = [elevation, one_hot_sites, one_hot_domains]

    neighbor_pool = df

    #If there are no neighbors, return 0's
    if neighbor_pool.empty:
        feature_array = np.zeros((k_neighbors, model_class.ensemble_model.output.shape[1]))
        distances = np.repeat(9999, k_neighbors)
    else:
        raster = rasterio.open(sensor_path)
        feature_array, distances = predict_neighbors(target, metadata=metadata, HSI_size=HSI_size, raster=raster, neighbor_pool=neighbor_pool, model=model_class.ensemble_model, k_neighbors=k_neighbors)

    #enforce dtype
    return feature_array, distances

def predict_dataframe(df, model_class, hyperspectral_pool, site_label_dict, domain_label_dict, HSI_size=20, k_neighbors=5):
    """Iterate through a geopandas dataframe and get neighbors for each tree.
    Args:
    df: a geopandas dataframe
    model_class: A deeptreeattention model class to extract layer features
    hyperspectral_pool: glob dir to search for sensor files
    HSI_size: size of HSI crop
    site_label_dict: dictionary of numeric site labels
    domain_label_dict: dictionary of numeric domain labels
    k_neighbors: number of neighbors to extract
    Returns:
    feature_array: a feature matrix of encoded bottleneck layer
    """

    #for each target in a dataframe, lookup the correct tile
    neighbor_features = {}
    for index, row in df.iterrows():
        row = pd.DataFrame(row).transpose()
        neighbor_features[index] = extract_features(
            df=df,
            x=row["individual"].values[0],
            model_class=model_class,
            hyperspectral_pool=hyperspectral_pool,
            site_label_dict=site_label_dict,
            domain_label_dict=domain_label_dict,
            k_neighbors=k_neighbors
        )

    return neighbor_features

max_stars_repo_name: zoeyingz/DeepTreeAttention

max_issues_repo_name: zoeyingz/DeepTreeAttention

max_forks_repo_name: zoeyingz/DeepTreeAttention


text:
from __main__ import vtk, qt, ctk, slicer
import string
import numpy
import math
import operator
import collections
import time

# from decimal import *

class TextureGLCM:
    def __init__(self, grayLevels, numGrayLevels, parameterMatrix, parameterMatrixCoordinates, parameterValues,
                 allKeys, checkStopProcessFunction):
        self.textureFeaturesGLCM = collections.OrderedDict()
        self.textureFeaturesGLCMTiming = collections.OrderedDict()

        self.textureFeaturesGLCM["Autocorrelation"] = "self.autocorrelationGLCM(self.P_glcm, self.prodMatrix)"
        self.textureFeaturesGLCM[
            "Cluster Prominence"] = "self.clusterProminenceGLCM(self.P_glcm, self.sumMatrix, self.ux, self.uy)"
        self.textureFeaturesGLCM[
            "Cluster Shade"] = "self.clusterShadeGLCM(self.P_glcm, self.sumMatrix, self.ux, self.uy)"
        self.textureFeaturesGLCM[
            "Cluster Tendency"] = "self.clusterTendencyGLCM(self.P_glcm, self.sumMatrix, self.ux, self.uy)"
        self.textureFeaturesGLCM["Contrast"] = "self.contrastGLCM(self.P_glcm, self.diffMatrix)"
        self.textureFeaturesGLCM[
            "Correlation"] = "self.correlationGLCM(self.P_glcm, self.prodMatrix, self.ux, self.uy, self.sigx, self.sigy)"
        self.textureFeaturesGLCM["Difference Entropy"] = "self.differenceEntropyGLCM(self.pxSuby, self.eps)"
        self.textureFeaturesGLCM["Dissimilarity"] = "self.dissimilarityGLCM(self.P_glcm, self.diffMatrix)"
        self.textureFeaturesGLCM["Energy (GLCM)"] = "self.energyGLCM(self.P_glcm)"
        self.textureFeaturesGLCM["Entropy(GLCM)"] = "self.entropyGLCM(self.P_glcm, self.pxy, self.eps)"
        self.textureFeaturesGLCM["Homogeneity 1"] = "self.homogeneity1GLCM(self.P_glcm, self.diffMatrix)"
        self.textureFeaturesGLCM["Homogeneity 2"] = "self.homogeneity2GLCM(self.P_glcm, self.diffMatrix)"
        self.textureFeaturesGLCM["IMC1"] = "self.imc1GLCM(self.HXY, self.HXY1, self.HX, self.HY)"
        # self.textureFeaturesGLCM["IMC2"] = "sum(imc2)/len(imc2)" #"self.imc2GLCM(self,)"  # produces a calculation error
        self.textureFeaturesGLCM["IDMN"] = "self.idmnGLCM(self.P_glcm, self.diffMatrix, self.Ng)"
        self.textureFeaturesGLCM["IDN"] = "self.idnGLCM(self.P_glcm, self.diffMatrix, self.Ng)"
        self.textureFeaturesGLCM["Inverse Variance"] = "self.inverseVarianceGLCM(self.P_glcm, self.diffMatrix, self.Ng)"
        self.textureFeaturesGLCM["Maximum Probability"] = "self.maximumProbabilityGLCM(self.P_glcm)"
        self.textureFeaturesGLCM["Sum Average"] = "self.sumAverageGLCM(self.pxAddy, self.kValuesSum)"
        self.textureFeaturesGLCM["Sum Entropy"] = "self.sumEntropyGLCM(self.pxAddy, self.eps)"
        self.textureFeaturesGLCM["Sum Variance"] = "self.sumVarianceGLCM(self.pxAddy, self.kValuesSum)"
        self.textureFeaturesGLCM["Variance (GLCM)"] = "self.varianceGLCM(self.P_glcm, self.ivector, self.u)"

        self.grayLevels = grayLevels
        self.parameterMatrix = parameterMatrix
        self.parameterMatrixCoordinates = parameterMatrixCoordinates
        self.parameterValues = parameterValues
        self.Ng = numGrayLevels
        self.keys = set(allKeys).intersection(list(self.textureFeaturesGLCM.keys()))
        # Callback function to stop the process if the user decided so. CalculateCoefficients can take a long time to run...
        self.checkStopProcessFunction = checkStopProcessFunction

    def CalculateCoefficients(self, printTiming=False):
        """ Calculate generic coefficients that will be reused in different markers
        IMPORTANT!! This method is VERY inefficient when the nodule is big (because
        of the function calculate_glcm at least). If these
        statistics are required it would probably need some optimizations
        :return:
        """
        # generate container for GLCM Matrices, self.P_glcm
        # make distance an optional parameter, as in: distances = numpy.arange(parameter)
        distances = numpy.array([1])
        directions = 26
        self.P_glcm = numpy.zeros((self.Ng, self.Ng, distances.size, directions))
        t1 = time.time()
        self.P_glcm = self.calculate_glcm(self.grayLevels, self.parameterMatrix, self.parameterMatrixCoordinates,
                                          distances, directions, self.Ng, self.P_glcm)
        if printTiming:
            print(("- Time to calculate glmc matrix: {0} secs".format(time.time() - t1)))
        # make each GLCM symmetric an optional parameter
        # if symmetric:
        # Pt = numpy.transpose(P, (1, 0, 2, 3))
        # P = P + Pt

        ##Calculate GLCM Coefficients
        self.ivector = numpy.arange(1, self.Ng + 1)  # shape = (self.Ng, distances.size, directions)
        self.jvector = numpy.arange(1, self.Ng + 1)  # shape = (self.Ng, distances.size, directions)
        self.eps = numpy.spacing(1)

        self.prodMatrix = numpy.multiply.outer(self.ivector, self.jvector)  # shape = (self.Ng, self.Ng)
        self.sumMatrix = numpy.add.outer(self.ivector, self.jvector)  # shape = (self.Ng, self.Ng)
        self.diffMatrix = numpy.absolute(numpy.subtract.outer(self.ivector, self.jvector))  # shape = (self.Ng, self.Ng)
        self.kValuesSum = numpy.arange(2, (self.Ng * 2) + 1)  # shape = (2*self.Ng-1)
        self.kValuesDiff = numpy.arange(0, self.Ng)  # shape = (self.Ng-1)

        # shape = (distances.size, directions)
        self.u = self.P_glcm.mean(0).mean(0)
        # marginal row probabilities #shape = (self.Ng, distances.size, directions)
        self.px = self.P_glcm.sum(1)
        # marginal column probabilities #shape = (self.Ng, distances.size, directions)
        self.py = self.P_glcm.sum(0)

        # shape = (distances.size, directions)
        self.ux = self.px.mean(0)
        # shape = (distances.size, directions)
        self.uy = self.py.mean(0)

        # shape = (distances.size, directions)
        self.sigx = self.px.std(0)
        # shape = (distances.size, directions)
        self.sigy = self.py.std(0)

        # shape = (2*self.Ng-1, distances.size, directions)
        self.pxAddy = numpy.array([numpy.sum(self.P_glcm[self.sumMatrix == k], 0) for k in self.kValuesSum])
        # shape = (self.Ng, distances.size, directions)
        self.pxSuby = numpy.array([numpy.sum(self.P_glcm[self.diffMatrix == k], 0) for k in self.kValuesDiff])

        # entropy of self.px #shape = (distances.size, directions)
        self.HX = (-1) * numpy.sum((self.px * numpy.where(self.px != 0, numpy.log2(self.px), numpy.log2(self.eps))), 0)
        # entropy of py #shape = (distances.size, directions)
        self.HY = (-1) * numpy.sum((self.py * numpy.where(self.py != 0, numpy.log2(self.py), numpy.log2(self.eps))), 0)
        # shape = (distances.size, directions)
        self.HXY = (-1) * numpy.sum(
            numpy.sum((self.P_glcm * numpy.where(self.P_glcm != 0, numpy.log2(self.P_glcm), numpy.log2(self.eps))), 0),
            0)

        self.pxy = numpy.zeros(self.P_glcm.shape)  # shape = (self.Ng, self.Ng, distances.size, directions)
        for a in range(directions):
            for g in range(distances.size):
                self.pxy[:, :, g, a] = numpy.multiply.outer(self.px[:, g, a], self.py[:, g, a])

        self.HXY1 = (-1) * numpy.sum(
            numpy.sum((self.P_glcm * numpy.where(self.pxy != 0, numpy.log2(self.pxy), numpy.log2(self.eps))), 0),
            0)  # shape = (distances.size, directions)
        self.HXY2 = (-1) * numpy.sum(
            numpy.sum((self.pxy * numpy.where(self.pxy != 0, numpy.log2(self.pxy), numpy.log2(self.eps))), 0),
            0)  # shape = (distances.size, directions)
        if printTiming:
            print(("- Time to calculate total glmc coefficients: {0} secs".format(time.time() - t1)))

    def autocorrelationGLCM(self, P_glcm, prodMatrix, meanFlag=True):
        ac = numpy.sum(numpy.sum(P_glcm * prodMatrix[:, :, None, None], 0), 0)
        if meanFlag:
            return (ac.mean())
        else:
            return ac

    def clusterProminenceGLCM(self, P_glcm, sumMatrix, ux, uy, meanFlag=True):
        # Need to validate function
        cp = numpy.sum(
            numpy.sum((P_glcm * ((sumMatrix[:, :, None, None] - ux[None, None, :, :] - uy[None, None, :, :]) ** 4)), 0),
            0)
        if meanFlag:
            return (cp.mean())
        else:
            return cp

    def clusterShadeGLCM(self, P_glcm, sumMatrix, ux, uy, meanFlag=True):
        # Need to validate function
        cs = numpy.sum(
            numpy.sum((P_glcm * ((sumMatrix[:, :, None, None] - ux[None, None, :, :] - uy[None, None, :, :]) ** 3)), 0),
            0)
        if meanFlag:
            return (cs.mean())
        else:
            return cs

    def clusterTendencyGLCM(self, P_glcm, sumMatrix, ux, uy, meanFlag=True):
        # Need to validate function
        ct = numpy.sum(
            numpy.sum((P_glcm * ((sumMatrix[:, :, None, None] - ux[None, None, :, :] - uy[None, None, :, :]) ** 2)), 0),
            0)
        if meanFlag:
            return (ct.mean())
        else:
            return ct

    def contrastGLCM(self, P_glcm, diffMatrix, meanFlag=True):
        cont = numpy.sum(numpy.sum((P_glcm * (diffMatrix[:, :, None, None] ** 2)), 0), 0)
        if meanFlag:
            return (cont.mean())
        else:
            return cont

    def correlationGLCM(self, P_glcm, prodMatrix, ux, uy, sigx, sigy, meanFlag=True):
        # Need to validate function
        uxy = ux * uy
        sigxy = sigx * sigy
        corr = numpy.sum(
            numpy.sum(((P_glcm * prodMatrix[:, :, None, None] - uxy[None, None, :, :]) / (sigxy[None, None, :, :])), 0),
            0)
        if meanFlag:
            return (corr.mean())
        else:
            return corr

    def differenceEntropyGLCM(self, pxSuby, eps, meanFlag=True):
        difent = numpy.sum((pxSuby * numpy.where(pxSuby != 0, numpy.log2(pxSuby), numpy.log2(eps))), 0)
        if meanFlag:
            return (difent.mean())
        else:
            return difent

    def dissimilarityGLCM(self, P_glcm, diffMatrix, meanFlag=True):
        dis = numpy.sum(numpy.sum((P_glcm * diffMatrix[:, :, None, None]), 0), 0)
        if meanFlag:
            return (dis.mean())
        else:
            return dis

    def energyGLCM(self, P_glcm, meanFlag=True):
        ene = numpy.sum(numpy.sum((P_glcm ** 2), 0), 0)
        if meanFlag:
            return (ene.mean())
        else:
            return ene

    def entropyGLCM(self, P_glcm, pxy, eps, meanFlag=True):
        ent = -1 * numpy.sum(numpy.sum((P_glcm * numpy.where(pxy != 0, numpy.log2(pxy), numpy.log2(eps))), 0), 0)
        if meanFlag:
            return (ent.mean())
        else:
            return ent

    def homogeneity1GLCM(self, P_glcm, diffMatrix, meanFlag=True):
        homo1 = numpy.sum(numpy.sum((P_glcm / (1 + diffMatrix[:, :, None, None])), 0), 0)
        if meanFlag:
            return (homo1.mean())
        else:
            return homo1

    def homogeneity2GLCM(self, P_glcm, diffMatrix, meanFlag=True):
        homo2 = numpy.sum(numpy.sum((P_glcm / (1 + diffMatrix[:, :, None, None] ** 2)), 0), 0)
        if meanFlag:
            return (homo2.mean())
        else:
            return homo2

    def imc1GLCM(self, HXY, HXY1, HX, HY, meanFlag=True):
        imc1 = (self.HXY - self.HXY1) / numpy.max(([self.HX, self.HY]), 0)
        if meanFlag:
            return (imc1.mean())
        else:
            return imc1

            # def imc2GLCM(self,):
            # imc2[g,a] = ( 1-numpy.e**(-2*(HXY2[g,a]-HXY[g,a])) )**(0.5) #nan value too high

            # produces Nan(square root of a negative)
            # exponent = decimal.Decimal( -2*(HXY2[g,a]-self.HXY[g,a]) )
            # imc2.append( ( decimal.Decimal(1)-decimal.Decimal(numpy.e)**(exponent) )**(decimal.Decimal(0.5)) )

            # if meanFlag:
            # return (homo2.mean())
            # else:
            # return homo2

    def idmnGLCM(self, P_glcm, diffMatrix, Ng, meanFlag=True):
        idmn = numpy.sum(numpy.sum((P_glcm / (1 + ((diffMatrix[:, :, None, None] ** 2) / (Ng ** 2)))), 0), 0)
        if meanFlag:
            return (idmn.mean())
        else:
            return idmn

    def idnGLCM(self, P_glcm, diffMatrix, Ng, meanFlag=True):
        idn = numpy.sum(numpy.sum((P_glcm / (1 + (diffMatrix[:, :, None, None] / Ng))), 0), 0)
        if meanFlag:
            return (idn.mean())
        else:
            return idn

    def inverseVarianceGLCM(self, P_glcm, diffMatrix, Ng, meanFlag=True):
        maskDiags = numpy.ones(diffMatrix.shape, dtype=bool)
        maskDiags[numpy.diag_indices(Ng)] = False
        inv = numpy.sum((P_glcm[maskDiags] / (diffMatrix[:, :, None, None] ** 2)[maskDiags]), 0)
        if meanFlag:
            return (inv.mean())
        else:
            return inv

    def maximumProbabilityGLCM(self, P_glcm, meanFlag=True):
        maxprob = P_glcm.max(0).max(0)
        if meanFlag:
            return (maxprob.mean())
        else:
            return maxprob

    def sumAverageGLCM(self, pxAddy, kValuesSum, meanFlag=True):
        sumavg = numpy.sum((kValuesSum[:, None, None] * pxAddy), 0)
        if meanFlag:
            return (sumavg.mean())
        else:
            return sumavg

    def sumEntropyGLCM(self, pxAddy, eps, meanFlag=True):
        sumentr = (-1) * numpy.sum((pxAddy * numpy.where(pxAddy != 0, numpy.log2(pxAddy), numpy.log2(eps))), 0)
        if meanFlag:
            return (sumentr.mean())
        else:
            return sumentr

    def sumVarianceGLCM(self, pxAddy, kValuesSum, meanFlag=True):
        sumvar = numpy.sum((pxAddy * ((kValuesSum[:, None, None] - kValuesSum[:, None, None] * pxAddy) ** 2)), 0)
        if meanFlag:
            return (sumvar.mean())
        else:
            return sumvar

    def varianceGLCM(self, P_glcm, ivector, u, meanFlag=True):
        vari = numpy.sum(numpy.sum((P_glcm * ((ivector[:, None] - u) ** 2)[:, None, None, :]), 0), 0)
        if meanFlag:
            return (vari.mean())
        else:
            return vari

    def calculate_glcm(self, grayLevels, matrix, matrixCoordinates, distances, directions, numGrayLevels, out):
        # VERY INEFFICIENT!!
        # 26 GLCM matrices for each image for every direction from the voxel
        # (26 for each neighboring voxel from a reference voxel centered in a 3x3 cube)
        # for GLCM matrices P(i,j;gamma, a), gamma = 1, a = 1...13

        angles_idx = 0
        distances_idx = 0
        r = 0
        c = 0
        h = 0
        rows = matrix.shape[2]
        cols = matrix.shape[1]
        height = matrix.shape[0]
        row = 0
        col = 0
        height = 0

        angles = numpy.array([(1, 0, 0),
                              (-1, 0, 0),
                              (0, 1, 0),
                              (0, -1, 0),
                              (0, 0, 1),
                              (0, 0, -1),
                              (1, 1, 0),
                              (-1, 1, 0),
                              (1, -1, 0),
                              (-1, -1, 0),
                              (1, 0, 1),
                              (-1, 0, 1),
                              (1, 0, -1),
                              (-1, 0, -1),
                              (0, 1, 1),
                              (0, -1, 1),
                              (0, 1, -1),
                              (0, -1, -1),
                              (1, 1, 1),
                              (-1, 1, 1),
                              (1, -1, 1),
                              (1, 1, -1),
                              (-1, -1, 1),
                              (-1, 1, -1),
                              (1, -1, -1),
                              (-1, -1, -1)])

        indices = list(zip(*matrixCoordinates))

        for iteration in range(len(indices)):
        #for h, c, r in indices:
            h, c, r = indices[iteration]
            for angles_idx in range(directions):
                angle = angles[angles_idx]

                for distances_idx in range(distances.size):
                    distance = distances[distances_idx]

                    i = matrix[h, c, r]
                    i_idx = numpy.nonzero(grayLevels == i)

                    row = r + angle[2]
                    col = c + angle[1]
                    height = h + angle[0]

                    # Can introduce Parameter Option for reference voxel(i) and neighbor voxel(j):
                    # Intratumor only: i and j both must be in tumor ROI
                    # Tumor+Surrounding: i must be in tumor ROI but J does not have to be
                    if row >= 0 and row < rows and col >= 0 and col < cols:
                        if tuple((height, col, row)) in indices:
                            j = matrix[height, col, row]
                            j_idx = numpy.nonzero(grayLevels == j)
                            # if i >= grayLevels.min and i <= grayLevels.max and j >= grayLevels.min and j <= grayLevels.max:
                            out[i_idx, j_idx, distances_idx, angles_idx] += 1
            # Check if the user has cancelled the process
            if iteration % 10 == 0:
                self.checkStopProcessFunction()

        return (out)

    def EvaluateFeatures(self, printTiming=False, checkStopProcessFunction=None):
        # Remove all the keys that must not be evaluated
        for key in set(self.textureFeaturesGLCM.keys()).difference(self.keys):
            self.textureFeaturesGLCM[key] = None

        if not self.keys:
            if not printTiming:
                return self.textureFeaturesGLCM
            else:
                return self.textureFeaturesGLCM, self.textureFeaturesGLCMTiming
        # normalization step:
        self.CalculateCoefficients(printTiming)

        if not printTiming:
            # Evaluate dictionary elements corresponding to user selected keys
            for key in self.keys:
                self.textureFeaturesGLCM[key] = eval(self.textureFeaturesGLCM[key])
                if checkStopProcessFunction is not None:
                    checkStopProcessFunction()
            return self.textureFeaturesGLCM
        else:
            # Evaluate dictionary elements corresponding to user selected keys
            for key in self.keys:
                t1 = time.time()
                self.textureFeaturesGLCM[key] = eval(self.textureFeaturesGLCM[key])
                self.textureFeaturesGLCMTiming[key] = time.time() - t1
                if checkStopProcessFunction is not None:
                    checkStopProcessFunction()
            return self.textureFeaturesGLCM, self.textureFeaturesGLCMTiming

max_stars_repo_name: hung-lab/SlicerCIP

max_issues_repo_name: hung-lab/SlicerCIP

max_forks_repo_name: hung-lab/SlicerCIP