text:
#-----------------------------------------------------------------------------
# Copyright (c) 2014, Ryan Volz
# All rights reserved.
#
# Distributed under the terms of the BSD 3-Clause ("BSD New") license.
#
# The full license is in the LICENSE file, distributed with this software.
#-----------------------------------------------------------------------------

import numpy as np
import os
import glob
import fnmatch

from echolect.core.indexing import (find_index, slice_by_value, wrap_check_start,
                                    wrap_check_stop)

from . import raw_parsing

__all__ = ['file_times', 'find_files', 'find_files_recursive', 'map_file_blocks',
           'read_voltage',
           'voltage_reader']

#******** See raw_parsing.py for details on Jicamarca raw data format ***************

def find_files(fdir, pattern='D*.r'):
    files = glob.glob(os.path.join(fdir, pattern))
    files.sort()
    return np.asarray(files)

def find_files_recursive(fdir, pattern='D*.r'):
    files = []
    for dirpath, dirnames, filenames in os.walk(fdir):
        for filename in fnmatch.filter(filenames, pattern):
            files.append(os.path.join(dirpath, filename))
    files.sort()
    return np.asarray(files)

def file_times(files):
    file_times = []
    # find time that begins each file
    for fpath in files:
        with open(fpath, 'rb') as f:
            h = raw_parsing.read_first_header(f)
            time = raw_parsing.parse_time(h)
            file_times.append(time)

    return np.asarray(file_times)

def map_file_blocks(fpath):
    # read all of the headers, block times, and location for start of each data block
    headers = []
    block_times = []
    data_start_bytes = []
    with open(fpath, 'rb') as f:
        # get necesary info from first header
        h = raw_parsing.read_first_header(f)
        headers.append(h)
        time = raw_parsing.parse_time(h)
        block_times.append(time)
        data_start_bytes.append(f.tell())

        block_size = h['nSizeOfDataBlock']
        while True:
            # skip over the previous block of data
            f.seek(block_size, 1)

            # read the block's header
            try:
                h = raw_parsing.read_basic_header(f)
            except EOFError:
                break

            time = raw_parsing.parse_time(h)
            # check validity of header
            # assume that if time is 0, the subsequent block was
            # not written and hence EOF has been reached
            if time == 0:
                break

            headers.append(h)
            block_times.append(time)
            data_start_bytes.append(f.tell())

    headers = np.asarray(headers)
    block_times = np.asarray(block_times)
    data_start_bytes = np.asarray(data_start_bytes)

    return headers, block_times, data_start_bytes

class voltage_reader(object):
    def __init__(self, fpath):
        self.fpath = fpath
        headers, block_times, data_start_bytes = map_file_blocks(fpath)
        self.headers = headers
        self.block_times = block_times
        self.data_start_bytes = data_start_bytes

        h = headers[0]
        raw_dtype, dtype = raw_parsing.parse_dtype(h)
        self.raw_dtype = raw_dtype
        self.dtype = dtype

        block_shape = raw_parsing.parse_block_shape(h)
        self.block_shape = block_shape
        self.nprofiles_per_block = block_shape[0]
        self.nsamples_per_profile = block_shape[1]
        self.nchannels = block_shape[2]
        self.nitems_per_profile = block_shape[1]*block_shape[2]
        self.nblocks = len(block_times)
        self.nprofiles = self.nblocks*block_shape[0]
        self.profile_bytes = self.nitems_per_profile*self.raw_dtype.itemsize
        self.shape = (self.nprofiles, self.nchannels, self.nsamples_per_profile)

        self.ts = raw_parsing.parse_ts(h)
        self.ipp = raw_parsing.parse_ipp(h)
        self.r = raw_parsing.parse_range_index(h)

    def __getitem__(self, key):
        if not isinstance(key, tuple):
            pidx = key
            cidx = slice(None)
            sidx = slice(None)
        else:
            lkey = len(key)
            if lkey < 1 or lkey > 3:
                raise IndexError('Wrong number of indices')
            elif lkey == 1:
                pidx = key[0]
                cidx = slice(None)
                sidx = slice(None)
            elif lkey == 2:
                pidx = key[0]
                cidx = key[1]
                sidx = slice(None)
            else:
                pidx = key[0]
                cidx = key[1]
                sidx = key[2]

        if isinstance(pidx, int):
            return self.read_voltage(pidx, chan_idx=cidx, sample_idx=sidx)

        # pidx must be a slice object
        start, stop, step = pidx.indices(self.shape[0])
        return self.read_voltage(start, stop, step, cidx, sidx)

    def _read_from_block(self, block_num, start, stop, step, chan_idx, sample_idx):
        num = stop - start
        # so we drop dimension when num == 1
        if num == 1:
            prof_idx = 0
        else:
            prof_idx = slice(0, num, step)

        fstart = self.data_start_bytes[block_num] + self.profile_bytes*start
        with open(self.fpath, 'rb') as f:
            f.seek(fstart)
            vlt_raw = np.fromfile(f, self.raw_dtype, num*self.nitems_per_profile)
        # data arranged by channel, then sample, then profile
        try:
            vlt_raw = vlt_raw.reshape(num, self.nsamples_per_profile, self.nchannels)
        except ValueError: # we didn't get the number of samples we expected, reshape fails
            raise EOFError('End of file reached. Could not read requested data.')
        # swap axes so it matches what we want, and slice as desired
        vlt_raw = vlt_raw.swapaxes(1, 2)[prof_idx, chan_idx, sample_idx]
        vlt = np.empty(vlt_raw.shape, dtype=self.dtype)
        vlt.real = vlt_raw['real']
        vlt.imag = vlt_raw['imag']

        return vlt

    def read_from_block(self, block_num, start, stop=None, step=1, nframes=1,
                        chan_idx=slice(None), sample_idx=slice(None)):
        start = wrap_check_start(self.nprofiles_per_block, start)
        if stop is None:
            stop = start + step*nframes
        else:
            stop = wrap_check_stop(self.nprofiles_per_block, stop)
        ## change ints to lists so that we don't lose dimensions when indexing
        #if isinstance(chan_idx, int):
            #chan_idx = [chan_idx]
        #if isinstance(sample_idx, int):
            #sample_idx = [sample_idx]
        return self._read_from_block(block_num, start, stop, step, chan_idx, sample_idx)

    def _read_from_blocks(self, blocknumstart, blockstart, blocknumend, blockstop,
                          step, chan_idx, sample_idx):
        if blocknumstart == blocknumend:
            return self.read_from_block(blocknumstart, blockstart, blockstop, step,
                                        chan_idx, sample_idx)

        start = blockstart
        vlt_all = []
        for bnum in xrange(blocknumstart, blocknumend + 1):
            if bnum == blocknumend:
                vlt_all.append(self.read_from_block(bnum, start, blockstop, step,
                                                    chan_idx, sample_idx))
            else:
                vlt_all.append(self.read_from_block(bnum, start, self.nprofiles_per_block, step,
                                                    chan_idx, sample_idx))
                # set start for (possibly != 0) based on step
                # step - ((nprofiles - start) % step) == (start - nprofiles) % step
                start = (start - self.nprofiles_per_block) % step

        return np.concatenate(vlt_all, axis=0)

    def read_voltage(self, start, stop=None, step=1, nframes=1,
                     chan_idx=slice(None), sample_idx=slice(None)):
        start = wrap_check_start(self.shape[0], start)
        if stop is None:
            stop = start + step*nframes
        else:
            stop = wrap_check_stop(self.shape[0], stop)
        ## change ints to lists so that we don't lose dimensions when indexing
        #if isinstance(chan_idx, int):
            #chan_idx = [chan_idx]
        #if isinstance(sample_idx, int):
            #sample_idx = [sample_idx]

        # find blocks for start and stop
        bstart, strt = divmod(start, self.nprofiles_per_block)
        # want block of last profile to include, hence profile number end = stop - 1
        end = stop - 1
        bend, nend = divmod(end, self.nprofiles_per_block)
        stp = nend + 1

        return self._read_from_blocks(bstart, strt, bend, stp, step, chan_idx, sample_idx)

def read_voltage(fpath, key=slice(None)):
    vlt_r = voltage_reader(fpath)
    return vlt_r[key]

max_stars_repo_name: ryanvolz/echolect

max_issues_repo_name: scivision/echolect

max_forks_repo_name: scivision/echolect


text:
__copyright__ = """

    Copyright [2018] [MR.Deemarc Burakitbumrung]

    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.

"""
__license__ = "Apache 2.0"

'''
converting video into image sequence
'''

import cv2
import numpy as np
import os
import glob
# Find OpenCV version
(major_ver, minor_ver, subminor_ver) = (cv2.__version__).split('.')

# get all the video from video_in folder
videos = glob.glob('./video_in/*.avi')
# videos = glob.glob('./video_in/*.mp4')

try:
    if not os.path.exists('img'):
        os.makedirs('img')
except OSError:
    print ('Error: Creating directory of data')

for video in videos:
    filename_ext = os.path.basename(video)
    filename, file_extension  = os.path.splitext(filename_ext)
    output_folder = './img/' + filename
    try:
        if not os.path.exists(output_folder):
            os.makedirs(output_folder)
    except OSError:
        print ('Error: Creating directory of data')

    # get video handle
    cap = cv2.VideoCapture(video)
    if int(major_ver)  < 3 :
        fps = cap.get(cv2.cv.CV_CAP_PROP_FPS)
        total = int(cap.get(cv2.cv.CAP_PROP_FRAME_COUNT))
        print("Frames per second using video.get(cv2.cv.CV_CAP_PROP_FPS): {0}".format(fps))
    else :
        fps = cap.get(cv2.CAP_PROP_FPS)
        total = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))
        print("Frames per second using video.get(cv2.CAP_PROP_FPS) : {0}".format(fps))

    # #for now let's use just 1 second of video
    # second = 1
    # num_img = int(round(fps*second))
    print("fps: " + str(fps))
    print("total frame: " + str(total))
    frame_num = 0
    while(True):
        ret, frame = cap.read()
        if not ret:
            print(ret)
            break
        output_filename = output_folder +'/'+ filename + '_' + ('%04d' %frame_num)+'.jpg'
        print ('Creating...' + output_filename)
        cv2.imwrite(output_filename, frame)
        frame_num += 1

    # When everything done, release the capture
    cap.release()
    cv2.destroyAllWindows()

ext: py

reason: import numpy


text:
"""
Functions to extract synapse and nuclear features.
"""

import numpy as np
import pandas as pd
from numpy import load
from numpy import genfromtxt
import json
from synspy.analyze.util import load_segment_status_from_csv, dump_segment_info_to_csv
import tifffile as tif
from matplotlib import pyplot as plt
from matplotlib.pyplot import imshow
import matplotlib.image as mpimg
from sklearn.preprocessing import normalize
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from scipy.spatial import cKDTree

import os.path
from os import path

def nearest_pairs(v1, v2, radius):
    """
    Adopted from synspy: https://github.com/informatics-isi-edu/synspy.git
    Find nearest k-dimensional point pairs between v1 and v2 and return via output arrays.

       Inputs:
         v1: array with first pointcloud with shape (n, k)
         kdt1: must be cKDTree(v1) for correct function
         v2: array with second pointcloud with shape (m, k)
         radius: maximum euclidean distance between points in a pair

       Use greedy algorithm to assign nearest neighbors without
       duplication of any point in more than one pair.

       Outputs:
         out1: for each point in kdt1, gives index of paired point from v2 or -1
         out2: for each point in v2, gives index of paired point from v1 or -1

    """
    depth = min(max(v1.shape[0], v2.shape[0]), 100)

    out1 = np.full((v1.shape[0], 1), -1)
    out2 = np.full((v2.shape[0], 1), -1)

    kdt1 = cKDTree(v1)
    dx, pairs = kdt1.query(v2, depth, distance_upper_bound=radius)
    for d in range(depth):
        for idx2 in np.argsort(dx[:, d]):
            if dx[idx2, d] < radius:
                if out2[idx2] == -1 and out1[pairs[idx2, d]] == -1:
                    out2[idx2] = pairs[idx2, d]
                    out1[pairs[idx2, d]] = idx2
    return out2

def get_centroids(parts):
    """
    Extracts potential synapse centroids from the npz.

    Parameters:
        parts ( dict ???? ): loaded npz array
    Returns:
        centroids_tiff: centroid coordinates in the Tiff space, in pixels
        slice_origin: the top left (???) corner of the segmentation region
    """

    centroids = parts['centroids'].astype(np.int32)

    props = json.loads(parts['properties'].tostring().decode('utf8'))
    slice_origin = np.array(props['slice_origin'], dtype=np.int32)

    # convert cropped centroids back to full SPIM voxel coords
    centroids_tiff = centroids + slice_origin

    return centroids_tiff

def get_labels(parts, csv_filename):
    """
    Extracts synapse labels from the npz and csv with segmentation.

    Parameters:
        parts ( dict ???? ): loaded npz array
        csv_filename (string) : path to csv file
    Returns:
        centroids_tiff: centroid coordinates in the Tiff space, in pixels
        slice_origin: the top left (???) corner of the segmentation region
    """
    centroids = parts['centroids'].astype(np.int32)

    props = json.loads(parts['properties'].tostring().decode('utf8'))
    slice_origin = np.array(props['slice_origin'], dtype=np.int32)

    # build dense status array from sparse CSV
    statuses, _ = load_segment_status_from_csv(centroids, slice_origin, csv_filename)

    # interpret status flag values
    is_synapse = statuses == 7
    labels = is_synapse.astype(int)

    return labels

def get_centroids_and_labels(csv_filename, npz_filename):
    """
    Returns centroids with labels.

    Parameters:
        csv_filename (string): path to segmentation csv
        npz_filename (string): path to npz with candidates
    Returns:
        centroids_tiff: centroid coordinates in the Tiff space, in pixels
        labels: labels of the centroids (1=synapse/0=not)
    """
    parts = np.load(npz_filename)
    centroids_tiff = get_centroids(parts)
    labels = get_labels(parts, csv_filename)
    return centroids_tiff, labels

# get the volume around one synapse with the bounding box padding
def get_one_volume(volume, centroid, box_size):
    """
    From volume, crops one box of size (padding*2 + 1) around the center centroid.

    Parameters:
        volume (array) : 3D image to crop from
        centroid (array) : 3D coordinates, ZYX, in volume to crop around
        box_size (array) : the box size in 3D , ZYX. All 3 must be odd.
        TODO : add this check

    Returns:
        array : cropped box
    """

    centroid = np.array(centroid).astype(int)
    padding = ((np.array(box_size) - 1)//2).astype(int)

    start = centroid - padding
    end = centroid + padding + 1

    cropped_volume = volume[start[0]:end[0], start[1]:end[1], start[2]:end[2]]

    return np.array(cropped_volume)

def get_centroids_from_csv(csv_filename, img_resolution):
    """
    Get centroids in tiff space directly from csv file.
    In csv the coordinates are in the physical space, so you need to know the image
    resolution that was used when producing csv.

    Parameters:
        csv_filename (string) : path and filename of the csv.
        img_resolution : image resolution in order ZYX.

    Returns:
        centroids_tiff (np.array) : array of centroids in tiff sace in pixels, in order ZYX.
    """

    col_list = ["Z", "Y", "X"]

    centroids_phys = pd.read_csv(csv_filename, usecols=col_list)
    centroids_tiff = np.round(
        centroids_phys.to_numpy() / np.array(img_resolution)).astype(int)

    return centroids_tiff

def get_intensity_from_csv(csv_filename, img_resolution):
    """
    Get centroids in tiff space directly from csv file.
    In csv the coordinates are in the physical space, so you need to know the image
    resolution that was used when producing csv.

    Parameters:
        csv_filename (string) : path and filename of the csv.
        img_resolution : image resolution in order ZYX.

    Returns:
        centroids_tiff (np.array) : array of centroids in tiff sace in pixels, in order ZYX.
    """

    col_list = ["raw core", "raw hollow"]

    intensity_df = pd.read_csv(csv_filename, usecols=col_list)
    centroids_core_hollow = intensity_df.to_numpy()

    return centroids_core_hollow

################################ END OF GOOD ##############################################
# get_centroids_from_csv: centroid coordinates in the Tiff space directly as in csv (Z,Y,X)

# write centroid coordinates in the Tiff space to csv file
def write_centroids_to_file(csvfilename, npz_filename, save_name):
    centroids_tiff, labels = get_centroids_and_labels(csvfilename, npz_filename)
    synapses = np.array(centroids_tiff[labels == 1])

    np.savetxt(save_name, synapses, delimiter=",", header="Z,Y,X")
    return centroids_tiff, labels, synapses

def saveResultsAsTiff(filename, segmentation, predictions):
    """
    Saves results as tiff

    Parameters:
    segmentation : name of the segmentation
    filename: full filename with the extention
    image_file : file of the size for the classifiacation
    """
    csvfilename, npz_filename, image_file_green, _, _ = get_files(segmentation)
    # get file info
    im_info = tif.TiffFile(image_file_green)
    series = im_info.series[0]
    # create results as volume
    centroids_tiff, labels = get_centroids_and_labels(csvfilename, npz_filename)
    labels = labels.astype(bool)

    # all available choises
    # volume = np.zeros(series.shape,dtype = np.int16)
    # volume = make_cube(volume, centroids_tiff, 1)
    # tif.imsave(f"{filename}_all_local_max.tif", volume, photometric='minisblack')
    # print("did all available choises")

    # human
    # volume = np.zeros(series.shape,dtype = np.int16)
    # volume = make_cube(volume, centroids_tiff[labels], 1)
    # tif.imsave(f"{filename}_human_grader.tif", volume, photometric='minisblack')
    # print("did human")

    # all predictions
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[predictions], 1)
    tif.imsave(f"{filename}_all.tif", volume, photometric='minisblack')
    print("did all predictions")

    # all correct predictions
    correct = np.logical_and(predictions, labels)
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[correct], 1)
    tif.imsave(f"{filename}_correct.tif", volume, photometric='minisblack')
    print("did all correct predictions")

    # false positive
    false_pos = np.logical_and(predictions, np.logical_not(labels))
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[false_pos], 1)
    tif.imsave(f"{filename}_false_pos.tif", volume, photometric='minisblack')
    print("did all false positives")

    # false negative
    false_neg = np.logical_and(labels, np.logical_not(predictions))
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[false_neg], 1)
    tif.imsave(f"{filename}_false_neg.tif", volume, photometric='minisblack')
    print("did all false negatives")

    # TODO : make similar only for the classifier and humans separately

def applyTransform(points, transform):
    """

    """
    points = np.append(points, np.ones((points.shape[0], 1)), axis=1)
    points = np.matmul(points, transform)
    return points[:, :-1]

def transformAndSaveResultsAsTiff(filename, segmentation, predictions, transform):
    """ Saves results as tiff:
    segmentation : name of the segmentation
    filename: full filename with the extention
    image_file : file of the size for the classifiacation
    """
    csvfilename, npz_filename, image_file_green, _, _ = get_files(segmentation)
    # get file info
    im_info = tif.TiffFile(image_file_green)
    series = im_info.series[0]
    # create results as volume
    centroids_tiff, labels = get_centroids_and_labels(csvfilename, npz_filename)
    centroids_tiff = applyTransform(centroids_tiff, transform)
    labels = labels.astype(bool)

    # all available choises
    # volume = np.zeros(series.shape,dtype = np.int16)
    # volume = make_cube(volume, centroids_tiff, 1)
    # tif.imsave(f"{filename}_all_local_max.tif", volume, photometric='minisblack')
    # print("did all available choises")

    # human
    # volume = np.zeros(series.shape,dtype = np.int16)
    # volume = make_cube(volume, centroids_tiff[labels], 1)
    # tif.imsave(f"{filename}_human_grader.tif", volume, photometric='minisblack')
    # print("did human")

    # all predictions
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[predictions], 1)
    tif.imsave(f"{filename}_all.tif", volume, photometric='minisblack')
    print("did all predictions")

    # all correct predictions
    correct = np.logical_and(predictions, labels)
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[correct], 1)
    tif.imsave(f"{filename}_correct.tif", volume, photometric='minisblack')
    print("did all correct predictions")

    # false positive
    false_pos = np.logical_and(predictions, np.logical_not(labels))
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[false_pos], 1)
    tif.imsave(f"{filename}_false_pos.tif", volume, photometric='minisblack')
    print("did all false positives")

    # false negative
    false_neg = np.logical_and(labels, np.logical_not(predictions))
    volume = np.zeros(series.shape, dtype=np.int16)
    volume = make_cube(volume, centroids_tiff[false_neg], 1)
    tif.imsave(f"{filename}_false_neg.tif", volume, photometric='minisblack')
    print("did all false negatives")

    # TODO : make similar only for the classifier and humans separately

def make_cube(volume, centroids, padding):
    # radius in pixels
    for centroid in centroids:
        start = centroid - padding
        end = centroid + padding + 1
        volume[start[0]:end[0], start[1]:end[1], start[2]:end[2]] = 100
        volume[centroid[0], start[1]:end[1], start[2]:end[2]] = 200
    return volume

def get_all_volumes(image_file, centroids, padding):
    volume = tif.imread(image_file)
    dims = 2 * padding + 1
    cropped_volumes = np.zeros((centroids.shape[0], dims[0], dims[1], dims[2]))

    for i_syn in range(centroids.shape[0]):
        d0 = centroids[i_syn, 0]
        d1 = centroids[i_syn, 1]
        d2 = centroids[i_syn, 2]

        d0_left = d0 - padding[0]
        d1_left = d1 - padding[1]
        d2_left = d2 - padding[2]
        d0_right = d0 + 1 + padding[0]
        d1_right = d1 + 1 + padding[1]
        d2_right = d2 + 1 + padding[2]

        cropped_volumes[i_syn, :, :] = np.squeeze(volume[d0_left:d0_right,
                                                  d1_left:d1_right,
                                                  d2_left:d2_right])
    cropped_volumes = np.squeeze(cropped_volumes)
    # cropped_volumes = np.expand_dims(cropped_volumes, axis=3)

    return cropped_volumes

# returns mean intensity of the square box within the intensity padding
def get_mean_intensity(volumes, intensity_padding):
    num_syn, d1, d2 = volumes.shape
    intensity = np.zeros((num_syn, 1))

    for i_syn in range(num_syn):
        center1 = d1 // 2
        center2 = d2 // 2

        d1_left = center1 - intensity_padding[1]
        d2_left = center2 - intensity_padding[2]

        d1_right = center1 + 1 + intensity_padding[1]
        d2_right = center2 + 1 + intensity_padding[2]

        intensity[i_syn] = np.mean(volumes[i_syn,
                                   d1_left:d1_right,
                                   d2_left:d2_right])

    return intensity

# 1. mean image subtraction - works well if your color
# or intensity distributions is not consistent
# throughout the image (e.g. only centered objects)

# 2. per-channel normalization (subtract mean, divide by standard deviation),
# pretty standard, useful for variable sized input where you can't use 1.

# 3. per-channel mean subtraction - good for variable sized input
# where you can't use 1 and don't want to make too many assumptions about the distribution.

# 4. whitening (turn the distribution into a normal distribution,
# sometimes as easy as normalization but only if it's already normally distributed).
# Maybe others can weigh in on cases where whitening is not a good idea.

def normalize_volumes1(volumes):
    volumes = volumes - np.mean(volumes, axis=(1, 2))[:, None, None]
    volumes = np.expand_dims(volumes, axis=3)
    return volumes

def split_data_equal(data, labels):
    X_1 = data[labels == 1]
    X_0 = data[labels == 0]
    keep = np.random.choice(range(X_0.shape[0]), size=X_1.shape[0], replace=False)
    X_0 = X_0[keep, :]
    new_label = np.concatenate((np.ones((X_1.shape[0], 1)),
                                np.zeros((X_0.shape[0], 1))), axis=0)
    X_train, X_test, y_train, y_test = train_test_split(
        np.concatenate((X_1, X_0), axis=0),
        new_label, test_size=0.1)
    return X_train, X_test, y_train, y_test

def split_data_2Xmore0(data, labels):
    X_1 = data[labels == 1]
    X_0 = data[labels == 0]
    weight_0 = 2
    keep = np.random.choice(range(X_0.shape[0]), size=weight_0 * X_1.shape[0],
                            replace=False)
    X_0 = X_0[keep, :]
    new_label = np.concatenate((np.ones((X_1.shape[0], 1)),
                                np.zeros((X_0.shape[0], 1))), axis=0)
    X_train, X_test, y_train, y_test = train_test_split(
        np.concatenate((X_1, X_0), axis=0),
        new_label, test_size=0.1)
    return X_train, X_test, y_train, y_test

# get mean, max, min and median nuclear intensity for a fish
def get_nuclear_stats(segmentation):
    _, _, _, _, nucfilename = get_files(segmentation)
    if path.exists(nucfilename):
        nuclear_data = pd.read_csv(nucfilename, skiprows=[1])
        Inuc_mean = nuclear_data.loc[:, "raw core"].mean()
        Inuc_med = nuclear_data.loc[:, "raw core"].median()
        Inuc_min = nuclear_data.loc[:, "raw core"].min()
        Inuc_max = nuclear_data.loc[:, "raw core"].max()
    else:
        Inuc_mean = 'NaN'
        Inuc_med = 'NaN'
        Inuc_min = 'NaN'
        Inuc_max = 'NaN'

    return Inuc_mean, Inuc_med, Inuc_min, Inuc_max

# get mean green intensity in a volume centered at synapse and int_padding size
# returns numpy array size (NumSyn,1)
def get_green_Intensity_stats(segmentation, int_padding):
    csvfilename, npz_filename, image_file_green, _, _ = get_files(segmentation)
    centroids, _ = get_centroids_and_labels(csvfilename, npz_filename)
    volumes = get_all_volumes(image_file_green, centroids, padding)
    # mean intensity in each volume
    intensity_green = get_mean_intensity(volumes, int_padding)
    return intensity_green

# get mean red intensity in a volume centered at synapse and int_padding size
# returns numpy array size (NumSyn,1)
def get_green_Intensity_stats(segmentation, int_padding):
    csvfilename, npz_filename, _, image_file_red, _ = get_files(segmentation)
    centroids, _ = get_centroids_and_labels(csvfilename, npz_filename)
    volumes = get_all_volumes(image_file_red, centroids, padding)
    intensity_red = get_mean_intensity(volumes, int_padding)
    return intensity_red

# weight_0: how many times more 0 then 1 in the train+test
# if weight_0 = 2 it is the same as split_data_2Xmore0 function
def split_data_more0(data, labels, weight_0):
    X_1 = data[labels == 1]
    X_0 = data[labels == 0]

    keep = np.random.choice(range(X_0.shape[0]), size=weight_0 * X_1.shape[0],
                            replace=False)
    X_0 = X_0[keep, :]
    new_label = np.concatenate((np.ones((X_1.shape[0], 1)),
                                np.zeros((X_0.shape[0], 1))), axis=0)
    X_train, X_test, y_train, y_test = train_test_split(
        np.concatenate((X_1, X_0), axis=0),
        new_label, test_size=0.1)
    return X_train, X_test, y_train, y_test

def plot_slice(data, title):
    plt.figure(figsize=(10, 30))
    imshow(data, cmap='gray')
    plt.title(title)
    plt.axis('off')

# get all necessary files for the given segmentation (adds directories and extentions)
#  returns strings with full paths to the files
# segmentation is in form  ImgZfDsy20170915D3A
def get_files(segmentation):
    csvfilename = 'C:\\Users\\nadtochi\\csci599\\meExperiment\\math650\\' + segmentation + '\\Syn' + segmentation + '.synapses-only.csv'
    nucfilename = 'C:\\Users\\nadtochi\\csci599\\meExperiment\\math650\\' + segmentation + '\\Nuc' + segmentation[
                                                                                0:-1] + 'B.nuclei-only.csv'
    npz_filename = 'C:\\Users\\nadtochi\\csci599\\meExperiment\\math650\\' + segmentation + '\\Syn' + segmentation + '.npz'
    image_file_green = 'C:\\Users\\nadtochi\\csci599\\meExperiment\\math650\\' + segmentation + '\\green.tif'
    image_file_red = 'C:\\Users\\nadtochi\\csci599\\meExperiment\\math650\\' + segmentation + '\\red2green.tif'

    return csvfilename, npz_filename, image_file_green, image_file_red, nucfilename

def get_files_from_cohort(segmentation):
    csvfilename = 'D:\\TR01\\1640Cohort\\CsvFiles\\Syn' + segmentation + '.synapses-only.csv'
    npz_filename = 'D:\\TR01\\1640Cohort\\NpzFiles\\Syn' + segmentation + '.npz'
    image_file_green = 'D:\\TR01\\1640Cohort\\GreenDataFR\\' + segmentation[
                                                               :-1] + '_green.tif'
    image_file_red = 'D:\\TR01\\1640Cohort\\R2G_Images\\' + segmentation[:-1] + '_r2g.tif'
    nucfilename = 'Missing'

    return csvfilename, npz_filename, image_file_green, image_file_red, nucfilename

# returns volumes and intenseties in two channels
def get_train_data(segmentations, padding, intensity_padding):
    num_datasets = len(segmentations)
    last_dat = num_datasets - 1
    csvfilename, npz_filename, image_file_green, image_file_r2g, nucfilename = get_files(
        segmentations[last_dat])

    # get centroids and labels
    centroids, labels = get_centroids_and_labels(csvfilename, npz_filename)
    num_syn = centroids.shape[0]
    # add fish ID
    centroids_and_fish = np.concatenate(
        (centroids, np.ones((num_syn, 1)) * (last_dat + 1)), axis=1)

    # get normalized image and intensity in grenn channel
    volumes = get_all_volumes(image_file_green, centroids, padding)
    intensity_green = get_mean_intensity(volumes, intensity_padding)
    volumes = normalize_volumes1(volumes)

    # get intensity in red channel
    volumes_red = get_all_volumes(image_file_r2g, centroids, padding)
    intensity_red = get_mean_intensity(volumes_red, intensity_padding)

    # get nuclear data
    # Inuc = np.zeros((1,4))
    # Inuc[0,0],Inuc[0,1],Inuc[0,2],Inuc[0,3] = get_nuclear_stats(segmentations[last_dat])
    # Inuc = np.repeat(Inuc,num_syn, axis=0)
    # centroids_nuc = get_centroids(csv_filename,npz_filename)

    for i_fish in range(last_dat):
        csvfilename, npz_filename, image_file_green, image_file_r2g, nucfilename = get_files(
            segmentations[i_fish])

        # get centroids and labels
        centroids, i_labels = get_centroids_and_labels(csvfilename, npz_filename)
        labels = np.concatenate((labels, i_labels), axis=0)
        num_syn = centroids.shape[0]
        # add fish ID
        i_centroids_and_fish = np.concatenate(
            (centroids, np.ones((num_syn, 1)) * (i_fish + 1)), axis=1)
        centroids_and_fish = np.concatenate((centroids_and_fish, i_centroids_and_fish),
                                            axis=0)

        # get normalized image and intensity in grenn channel
        i_volumes = get_all_volumes(image_file_green, centroids, padding)
        i_intensity_green = get_mean_intensity(i_volumes, intensity_padding)
        i_volumes = normalize_volumes1(i_volumes)

        volumes = np.concatenate((volumes, i_volumes), axis=0)
        intensity_green = np.concatenate((intensity_green, i_intensity_green), axis=0)

        # get intensity in red channel
        volumes_red = get_all_volumes(image_file_r2g, centroids, padding)
        i_intensity_red = get_mean_intensity(volumes_red, intensity_padding)
        intensity_red = np.concatenate((intensity_red, i_intensity_red), axis=0)

        # get nuclear data
        # i_Inuc = np.zeros((1,4))
        # i_Inuc[0,0],i_Inuc[0,1],i_Inuc[0,2],i_Inuc[0,3] = get_nuclear_stats(segmentations[i_fish])
        # i_Inuc = np.repeat(i_Inuc,num_syn, axis=0)
        # Inuc = np.concatenate((Inuc,i_Inuc),axis = 0)

    print(f'Potential synapses:{labels.shape[0]}')

    # wrap all intensity info into a data frame
    features_labels = ["Label", "Z", "Y", "X",
                       "FishN",  # "FishID",
                       "IntGrn_avg",  # "IntGrn_med","IntGrn_min","IntGrn_max",
                       "IntRed_avg"]  # , #"IntRed_med","IntRed_min","IntRed_max",
    # "IntNuc_avg","IntNuc_med","IntNuc_min","IntNuc_max"]
    synapse_data = np.concatenate([labels[:, np.newaxis], centroids_and_fish,
                                   intensity_green, intensity_red], axis=1)  # Inuc
    synapse_dataset = pd.DataFrame(synapse_data)
    synapse_dataset.columns = features_labels
    # synapse_dataset['Label'] = synapse_dataset.Label.astype('category')
    # synapse_dataset['FishN'] = synapse_dataset.FishN.astype('category')

    return volumes, synapse_dataset, labels

def predict_and_create_csv(model, segmentation, output_path, output_filename, padding):
    csvfilename, npz_filename, image_file, image_file_red = get_files(segmentation)
    centroids_tiff, eval_labels = get_centroids_and_labels(csvfilename, npz_filename)
    print('Image file : ', image_file)
    eval_data_green = get_all_volumes(image_file, centroids_tiff, padding)
    eval_data_green = normalize_volumes1(eval_data_green)
    predictions = model.predict(eval_data_green, batch_size=128)

    # np.savetxt(output_path+segmentation+"_IntensityProb.csv",
    #           np.concatenate((intensity_green,predictions),axis = 1), delimiter=",")

    mpp = np.array([0.4, 0.26, 0.26], dtype=np.float32)  # in ZYX packed order
    parts = np.load(npz_filename)
    props = json.loads(parts['properties'].tostring().decode('utf8'))

    centroids = parts['centroids'].astype(np.int32)
    measures = parts['measures'].astype(np.float32) * np.float32(
        props['measures_divisor'])
    slice_origin = np.array(props['slice_origin'], dtype=np.int32)

    # build dense status array from sparse CSV
    statuses, saved_parms = load_segment_status_from_csv(centroids, slice_origin,
                                                         csvfilename)

    thr = 0.6
    statuses[:] = 0
    predictionsArr = (np.concatenate(predictions, axis=0) > thr)
    # statuses[predictionsArr==0] = 5
    statuses[predictionsArr == 1] = 7
    print('Sum 0.6', sum(predictionsArr))
    dump_segment_info_to_csv(centroids, measures, statuses, slice_origin,
                             output_path + "thr0p6\\" + output_filename,
                             all_segments=False)

    thr = 0.85
    statuses[:] = 0
    predictionsArr = (np.concatenate(predictions, axis=0) > thr)
    # statuses[predictionsArr==0] = 5
    statuses[predictionsArr == 1] = 7
    print('Sum 0.85', sum(predictionsArr))
    dump_segment_info_to_csv(centroids, measures, statuses, slice_origin,
                             output_path + "thr0p85\\" + output_filename,
                             all_segments=False)

    thr = 0.95
    statuses[:] = 0
    predictionsArr = (np.concatenate(predictions, axis=0) > thr)
    # statuses[predictionsArr==0] = 5
    statuses[predictionsArr == 1] = 7
    print('Sum 0.95', sum(predictionsArr))
    dump_segment_info_to_csv(centroids, measures, statuses, slice_origin,
                             output_path + "thr0p95\\" + output_filename,
                             all_segments=False)

    print('Done')

# well... plot ROC curve ;)
def plot_roc_curve(fpr, tpr, label=None):
    plt.plot(fpr, tpr, linewidth=2, label=label)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.axis([0, 1, 0, 1])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')

# writes out a csv in Karl-readable format
def create_csv(predictionsArr, segmentation, output_filename):
    csvfilename, npz_filename, _, _ = get_files(segmentation)
    # load data from *.npz
    mpp = np.array([0.4, 0.26, 0.26], dtype=np.float32)  # in ZYX packed order
    parts = np.load(npz_filename)
    props = json.loads(parts['properties'].tostring().decode('utf8'))

    centroids = parts['centroids'].astype(np.int32)
    measures = parts['measures'].astype(np.float32) * np.float32(
        props['measures_divisor'])
    slice_origin = np.array(props['slice_origin'], dtype=np.int32)

    # build dense status array from sparse CSV and override statuses
    statuses, saved_parms = load_segment_status_from_csv(centroids, slice_origin,
                                                         csvfilename)
    statuses[:] = 0
    statuses[predictionsArr == 1] = 7

    dump_segment_info_to_csv(centroids, measures, statuses, slice_origin, output_filename,
                             all_segments=False)
    print('Saved ', output_filename)

def print_acc_different_threshold(predictions):
    thr = 0.5
    predictionsArr = (np.concatenate(predictions, axis=0) > thr)
    acc = accuracy_score(predictionsArr, labels)
    print('Threshold    acc')
    print(thr, ' ', acc)

reason: import numpy,from numpy,from scipy


text:
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
from sklearn.metrics import mean_squared_error

from config import config
from scripts.utils.fit_poly import fit_poly
from scripts.utils.diffuse import diffuse

def save_to_csv(df):
    """
    Saves DataFrame to csv file.
    :param df <pandas.DataFrame> A pandas DataFrame to save
    :return: pandas.DataFrame.to_csv
    """

    filename = f'diffusion_calculation.csv'
    save_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'outputs', filename)
    return df.to_csv(save_path, index=False)

def tellme(s):
    """
    Prints a string to a plot
    :param s: <string> The title to be printed
    """

    plt.title(s, fontsize=16)
    plt.draw()

def create_dataframe_to_save(df, best_fit_and_rmse):
    """
    Saves the results as a DataFrame
    :param df: <pandas.DataFrame> The dataframe with profile distance and elevation
    :param best_fit_and_rmse: <tuple> The best fit and its rmse
    :return df_to_save: <pandas.DataFrame> The dataframe to save
    """

    df_to_save = pd.DataFrame(
        {"profile": list(set(df['profile'])),
         "age": best_fit_and_rmse[0],
         "rmse": best_fit_and_rmse[1]})

    return df_to_save

def plot_profiles_get_crest_and_base(df, use_smooth):
    """
    Plot profiles and ask user to select crest and base
    :param df: <pandas.DataFrame> The dataframe with profile distance and elevation
    :return base, crest: <tuple> The distance and height of the base and crest on the profile
    """

    plt.figure(1)
    plt.plot(df["dist"], df["elevation"])
    if use_smooth:
        plt.plot(df["dist"], df["elevation_raw"])
    plt.xlabel('Distance (m)')
    plt.ylabel('Elevation (m)')
    tellme('Select the scrap crest, then scarp base. Click to begin.')

    pts = plt.ginput(3, timeout=0)  # it will wait for three clicks
    pts = np.array(pts)
    if len(pts) == 3:
        plt.close()

    # Saves the user specified location of the crest and the base
    base = [pts[1, 0], pts[1, 1]]
    crest = [pts[0, 0], pts[0, 1]]

    return (base, crest)

def calculate_rmse(df, height):
    """
    Calculate the root mean square error between the calculated elevation and dataframe
    :param df: <pandas.DataFrame> The dataframe with profile distance and elevation
    :return height: <list> The elevation values from the diffuse calculation
    :return rmse_list: <list> The rmse list
    """

    rmse_list = []

    for n in range(0, len(df)):
        rmse_list.append(mean_squared_error(df['elevation'], height[n], squared=False))

    return rmse_list

def find_best_fit(df, diffusion_coefficient, calculation_final_time, time_steps_between_plots, use_smooth, smooth_amount,
                  plot=False):
    """
    Find the best fit between the original profile and the diffuse profiles
    :param df: <pandas.DataFrame> The dataframe with profile distance and height
    :param diffusion_coefficient: <int> The diffusion coefficient in mm per year
    :param calculation_final_time: <int> The final time in kyr
    :param time_steps_between_plots: <int> the time step between plots
    :param plot: <boolean> Whether  to plot the diffusion profiles or not
    :return best_diff, best_rmse: <tuple> The best diffusion amount and its rmse
    """

    profiles = list(set(df['profile'].to_list()))

    best_diffusion_amount = []
    best_rmse = []

    for profile in profiles:

        profile_df = df[df["profile"] == profile].reset_index()

        if use_smooth:
            # Make new column for raw data and make the smoothed data the one used for the diffusion calculations
            profile_df['elevation_raw'] = profile_df['elevation']
            profile_df['elevation'] = profile_df['elevation_raw'].rolling(smooth_amount, min_periods=1).mean()

        base_and_crest = plot_profiles_get_crest_and_base(profile_df, use_smooth)

        profile_df_simple = profile_df.copy()

        rmse_upper, upper, upper_poly_height, p_upper = fit_poly(profile_df, base_and_crest[0][0], base_and_crest[1][0],
                                                                 where="upper")
        rmse_lower, lower, lower_poly_height, p_lower = fit_poly(profile_df, base_and_crest[0][0], base_and_crest[1][0],
                                                                 where="lower")

        profile_dist_crest = profile_df['dist'].values[int(base_and_crest[1][0])]
        profile_dist_base = profile_df['dist'].values[int(base_and_crest[0][0])]
        scarp_height = upper_poly_height[profile_dist_crest] - lower_poly_height[profile_dist_base]
        scarp_height_per_distance = scarp_height / (profile_dist_base - profile_dist_crest + 1)

        new_elevation = [0] * len(profile_df)

        num = 1
        for index, row in profile_df_simple.iterrows():
            if row['dist'] < profile_dist_crest:
                new_elevation[index] = upper_poly_height[index]
            if profile_dist_crest <= row['dist'] <= profile_dist_base:
                new_elevation[index] = upper_poly_height[profile_dist_crest] - (scarp_height_per_distance * num)
                num += 1
            if row['dist'] > profile_dist_base:
                new_elevation[index] = lower_poly_height[index]

        profile_df_simple['elevation'] = new_elevation

        change_in_distance = profile_df_simple['dist'].values[1] - profile_df_simple['dist'].values[0]

        diffusion_values = diffuse(change_in_distance, new_elevation, diffusion_coefficient, calculation_final_time,
                                   time_steps_between_plots, plot=False);

        diffusion_amount = diffusion_values[2] * diffusion_coefficient

        rmse_list = calculate_rmse(profile_df, diffusion_values[0])

        best_diffusion_amount.append(diffusion_amount[rmse_list.index(min(rmse_list))])

        best_rmse.append(min(rmse_list))

    return (best_diffusion_amount, best_rmse)

def calculate_diffusion_amounts(df, use_smooth, smooth_amount, config):
    """
    Runs pipeline and saves output
    :param df: <pandas.DataFrame> The dataframe with profile distance and elevation
    """

    best_fit_and_rmse = find_best_fit(df, config.diffusion_coefficient, config.calculation_final_time,
                                      config.time_steps_between_plots, use_smooth, smooth_amount,
                                      plot=False)
    df_to_save = create_dataframe_to_save(df, best_fit_and_rmse)
    save_to_csv(df_to_save)

if __name__ == "__main__":
    df = pd.read_csv('../data/Profiles2.csv').reset_index()
    use_smooth = True
    smooth_amount = 10
    calculate_diffusion_amounts(df, use_smooth, smooth_amount, config)

max_stars_repo_name: mshodge/scarp-diffusion

max_issues_repo_name: mshodge/scarp-diffusion

max_forks_repo_name: mshodge/scarp-diffusion


text:
import numpy as np # type: ignore
import attr
# from typing import Any, Optional, Sequence, Sized, Union, Iterable, Dict, List
# import scipy # type: ignore
import scipy.optimize # type: ignore
# from numpy.linalg import norm
from collections import namedtuple
from mdh.link import mdh_params
from mdh.link import RevoluteLink
from mdh.link import JointType
from scipy.spatial.transform import Rotation as R # type: ignore
from math import atan2, pi

"""

"""

# Pose = namedtuple('Pose', "position rotation")

class UnReachable(Exception):
    pass

@attr.s
# class KinematicChain(Sized, Iterable):
class KinematicChain:
    # _links = attr.ib(type=Sequence[Any])
    _links = attr.ib()
    _size = attr.ib()

    def __len__(self):
        """Enables the length function, returns number links"""
        # return len(self._links)
        return self._size

    def __iter__(self):
        """Enables iteration of joints:  for j in chain: print(j)"""
        for l in self._links:
            yield l

    def __getitem__(self, key):
        """Return a link"""
        return self._links[key]

    def transform(self, joints):
        """Calculates the transformation and returns a 4x4 matrix"""
        if len(joints) != self.size:
            raise Exception(f"inputs don't equal number of links {len(joints)} != {self.size}")

        t = np.eye(4)
        for l, q in zip(reversed(self._links), reversed(joints)):
            m = l.transform(q)
            t = m.dot(t)
        return t

    def forward(self, joints):
        """value?
        Solve the forward kinematics and returns a 4x4 matrix

        return: (position), (euler angles)
        """
        return self.transform(joints)

    def inverse(self, pos):
            """
            Solve the inverse kinematics and returns 1xN on success or None
            on failure.

            FIXME: this only works for when the foot is straight down
            """
            joint_limits = [[],[]] # type: Sequence[List]
            for l in self:
                joint_limits[0].append(l.min)
                joint_limits[1].append(l.max)
            # # j[0][3] = -0.00001
            # # j[1][3] = 0.00001
            # joint_limits = j
            # joint_limits = np.array(j) # type: np.ndarray
            # print(joint_limits)

            # r = R.from_matrix([[1,0,0],[0,0,-1],[0,1,0]])
            # rr = R.from_euler('z', atan2(pos[1], pos[0]))
            # c = rr*r

            # make desired transform pose (position, oritentation)
            # change the -90 x-axis to something else for non-down??
            # x = pos[0]
            # y = pos[1]
            # c = R.from_euler('zyx', [-90,atan2(y,x)*180/pi,90], degrees=True)
            # ppos = np.eye(4)
            # ppos[:3,:3] = c.as_matrix()
            # ppos[:3, 3] = pos

            start = [0]*len(self._links)

            result = scipy.optimize.least_squares(
                fun=_ik_cost_function,
                x0=np.array(start),
                bounds=joint_limits,
                args=(pos, self)
            )

            if not result.success:  # pragma: no cover
                raise UnReachable(f"Can't reach: {pos}m")
            # print(">> yeah!!!")
            # print(">>", np.rad2deg(result.x))
            actual_t = self.transform(result.x)
            actual_pos = actual_t #actual_t[:3,3]
            if np.allclose(actual_pos, pos, atol=1e-3):
                # print(result)
                return result.x

            raise UnReachable(f"Can't reach: {pos}")

    @property
    def size(self):
        return self._size

    @classmethod
    def from_parameters(cls, params):
        """Builds a KinematicChain object from an input"""
        links = []
        size = 0
        for l in params:
            if not isinstance(l, dict):
                raise Exception(f"Invalid parameters: {l}")

            for key in ['alpha', 'a', 'theta', 'd']:
                if key not in l:
                    raise Exception(f"Missing parameter: {key}")

            # ll = mdh_params(l['alpha'], l['a'], l['theta'], l['d'])
            # print(f">> params: {ll}")
            # print(l['alpha'])

            # if ll.type == JointType.revolute:
            # link = RevoluteLink(a=ll.a, alpha=ll.alpha, d=ll.d, theta=ll.theta)
            # link = RevoluteLink(alpha=ll.alpha, a=ll.a, theta=ll.theta, d=ll.d)
            link = RevoluteLink(l['alpha'], l['a'], l['theta'], l['d'])
            # elif ll.type == JointType.prismatic:
            #     raise NotImplementedError(f"from_parameters: {l.type}")
            # else:
            #     raise Exception(f"Invalid parameter: {l.type}")

            # if "max_min" in l:
            #     link.min = l["max_min"][1]
            #     link.max = l["max_min"][0]
            # else:
            link.min = -pi
            link.max = pi

            links.append(link)

            # if 'fixed' in l:
            #     if l['fixed'] is not True:
            #         size += 1
            # else:
            #     size += 1
            size += 1

        ret = cls(links, size)

        return ret

def _ik_cost_function(ik_q, pose, kc):
    """Cost function for least squares algorithm.
        - ik_q: calculated 4x4 matrix from optimization
        - pose: desired 4x4 matrix with position and orientation
        - kc: KinematicChain object that is trying to find this info
    """
    actual_t = kc.transform(ik_q)
    diff = np.abs(actual_t - pose)
    return diff.ravel()

max_stars_repo_name: MultipedRobotics/dh

max_issues_repo_name: MultipedRobotics/dh

max_forks_repo_name: MultipedRobotics/dh

reason: import numpy,from numpy,import scipy,from scipy


text:
import numpy as np
from .domain import Domain

class Supersonic(Domain):
    name = 'supersonic'
    left = 0
    right = 1.43
    bottom = 0
    top = 1.43
    R1 = 1.
    R2 = 1.384

    def in_domain(self, X,Y):

        R = np.sqrt( X**2 + Y**2)
        bval = np.logical_and( np.less_equal(self.R1 , R) , np.less_equal(R , self.R2) )
        return bval.astype(int)
    def bc_id(self,bid):
        return bid
    def bc(self,x,y):
        r = np.sqrt(x**2. +  y**2.)
        r1 = -1*(np.abs(r-self.R1) < 1e-14)
        r2 = -2*(np.abs(r-self.R2) < 1e-14)

        sanity_check = np.logical_not(np.logical_xor(r1,r2))
        num_wrong = np.sum(sanity_check)
        if(num_wrong > 0):
            print("DUPLICATE BOUNDARY CONDITIONS\n")
            quit()

        return r1+r2

    # -1 reflecting, -2 is inflow/outflow
    def vertex2bc(self, x, y):
        r = np.sqrt(x**2. +  y**2.)
        r1 = -1*(np.abs(r-self.R1) < 1e-14)
        r2 = -1*(np.abs(r-self.R2) < 1e-14)
        bot= -2*(np.abs(y < 1e-14))
        left=-2*(np.abs(x < 1e-14))
        bc = r1+r2+bot+left
        decided = np.where(np.logical_xor(np.logical_xor(np.logical_xor(r1, r2), bot), left))
        out = np.zeros(x.shape)
        out[decided] = bc[decided]
        return out

    def curved_points(self,wgt,X1,Y1,X2,Y2):
        bc1 = self.bc(X1,Y1)
        bc2 = self.bc(X2,Y2)
        sanity_check = np.logical_not( np.equal(bc1,bc2) )
        num_wrong = np.sum(sanity_check)
        if(num_wrong > 0):
            print("not the same\n")
            quit()

        radii = np.where(np.equal(bc1, -1) , self.R1, self.R2)
        angle1 = np.arctan2( Y1, X1 )
        angle2 = np.arctan2( Y2, X2 )
        angles = wgt[None,:]*angle1[:,None] + (1.-wgt[None,:])*angle2[:,None]

        points = np.zeros( (angles.shape[0], angles.shape[1], 2) )
        points[:,:,0] = radii[:,None] * np.cos(angles)
        points[:,:,1] = radii[:,None] * np.sin(angles)
        return points

max_stars_repo_name: andrewgiuliani/PyGrid


text:
import torch
import torchvision
from dataset import SegmentationDataset
import utils
import argparse
import time
import os
from vision.coco_utils import get_coco_api_from_dataset
from vision.coco_eval import CocoEvaluator
import vision.transforms as T
import math
from torchvision.models.detection.faster_rcnn import FastRCNNPredictor
from torchvision.models.detection.mask_rcnn import MaskRCNNPredictor
import sys
from PIL import Image

def get_model_instance_segmentation(num_classes):
    # load an instance segmentation model pre-trained pre-trained on COCO
    # TODO, seems only work with torchvision 0.7.0 and 0.6
    # model = torchvision.models.detection.maskrcnn_resnet50_fpn(pretrained=True)
    model = torchvision.models.detection.maskrcnn_resnet50_fpn(pretrained=True, trainable_backbone_layers=5)

    # get number of input features for the classifier
    in_features = model.roi_heads.box_predictor.cls_score.in_features
    # replace the pre-trained head with a new one
    model.roi_heads.box_predictor = FastRCNNPredictor(in_features, num_classes)

    # now get the number of input features for the mask classifier
    in_features_mask = model.roi_heads.mask_predictor.conv5_mask.in_channels
    hidden_layer = 256
    # and replace the mask predictor with a new one
    model.roi_heads.mask_predictor = MaskRCNNPredictor(in_features_mask,
                                                       hidden_layer,
                                                       num_classes)

    return model

def get_transform(train):
    transforms = []
    transforms.append(T.ToTensor())
    if train:
        transforms.append(T.RandomHorizontalFlip(0.5))
    return T.Compose(transforms)

def get_data_loader(dataset_root, batch_size, is_real):
    # use our dataset and defined transformations
    dataset = SegmentationDataset(dataset_root, get_transform(train=True), is_real)
    dataset_test = SegmentationDataset(dataset_root, get_transform(train=False), is_real)

    # split the dataset in train and test set
    indices = torch.arange(0, len(dataset)).tolist()
    test_point = 5
    dataset_sub_train = torch.utils.data.Subset(dataset, indices[test_point:])
    dataset_sub_test = torch.utils.data.Subset(dataset_test, indices[:test_point])

    # define training and validation data loaders
    data_loader_train = torch.utils.data.DataLoader(
        dataset_sub_train,
        batch_size=batch_size,
        shuffle=True,
        num_workers=7,
        drop_last=True,
        collate_fn=utils.collate_fn)
    data_loader_test = torch.utils.data.DataLoader(
        dataset_sub_test,
        batch_size=batch_size,
        shuffle=False,
        num_workers=1,
        drop_last=False,
        collate_fn=utils.collate_fn)

    return data_loader_train, data_loader_test

def train_one_epoch(model, optimizer, data_loader, device, epoch, print_freq, resume=False):
    '''
    https://github.com/pytorch/vision/blob/master/references/detection/engine.py
    '''
    model.train()
    metric_logger = utils.MetricLogger(delimiter="  ")
    metric_logger.add_meter('lr', utils.SmoothedValue(window_size=1, fmt='{value:.8f}'))
    header = 'Epoch: [{}]'.format(epoch)

    lr_scheduler = None
    if epoch == 0 and not resume:
        warmup_factor = 1. / 1000
        warmup_iters = min(1000, len(data_loader) - 1)

        lr_scheduler = utils.warmup_lr_scheduler(optimizer, warmup_iters, warmup_factor)

    for images, targets in metric_logger.log_every(data_loader, print_freq, header):
        images = list(image.to(device) for image in images)
        targets = [{k: v.to(device) for k, v in t.items()} for t in targets]

        loss_dict = model(images, targets)

        losses = sum(loss for loss in loss_dict.values())

        # reduce losses over all GPUs for logging purposes
        loss_dict_reduced = utils.reduce_dict(loss_dict)
        losses_reduced = sum(loss for loss in loss_dict_reduced.values())

        loss_value = losses_reduced.item()

        if not math.isfinite(loss_value):
            print("Loss is {}, stopping training".format(loss_value))
            print(loss_dict_reduced)
            sys.exit(1)

        optimizer.zero_grad()
        losses.backward()
        optimizer.step()

        if lr_scheduler is not None:
            lr_scheduler.step()

        metric_logger.update(loss=losses_reduced, **loss_dict_reduced)
        metric_logger.update(lr=optimizer.param_groups[0]["lr"])

    return metric_logger

def _get_iou_types(model):
    model_without_ddp = model
    if isinstance(model, torch.nn.parallel.DistributedDataParallel):
        model_without_ddp = model.module
    iou_types = ["bbox"]
    if isinstance(model_without_ddp, torchvision.models.detection.MaskRCNN):
        iou_types.append("segm")
    if isinstance(model_without_ddp, torchvision.models.detection.KeypointRCNN):
        iou_types.append("keypoints")
    return iou_types

@torch.no_grad()
def evaluate(model, data_loader, device):
    n_threads = torch.get_num_threads()
    # FIXME remove this and make paste_masks_in_image run on the GPU
    torch.set_num_threads(1)
    cpu_device = torch.device("cpu")
    model.eval()
    metric_logger = utils.MetricLogger(delimiter="  ")
    header = 'Test:'

    coco = get_coco_api_from_dataset(data_loader.dataset)
    iou_types = _get_iou_types(model)
    coco_evaluator = CocoEvaluator(coco, iou_types)

    for images, targets in metric_logger.log_every(data_loader, 100, header):
        images = list(img.to(device) for img in images)

        torch.cuda.synchronize()
        model_time = time.time()
        outputs = model(images)

        outputs = [{k: v.to(cpu_device) for k, v in t.items()} for t in outputs]
        model_time = time.time() - model_time

        res = {target["image_id"].item(): output for target, output in zip(targets, outputs)}
        evaluator_time = time.time()
        coco_evaluator.update(res)
        evaluator_time = time.time() - evaluator_time
        metric_logger.update(model_time=model_time, evaluator_time=evaluator_time)

    # gather the stats from all processes
    metric_logger.synchronize_between_processes()
    print("Averaged stats:", metric_logger)
    coco_evaluator.synchronize_between_processes()

    # accumulate predictions from all images
    coco_evaluator.accumulate()
    coco_evaluator.summarize()
    torch.set_num_threads(n_threads)
    return coco_evaluator

@torch.no_grad()
def test(model, data_loader, device):
    import cv2
    import numpy as np

    model.eval()

    print(len(data_loader.dataset))
    count = 0
    for images, targets in data_loader:
        for i in range(len(images)):
            print(targets[i]['image_id'])
            image = images[i]

            # target = targets[i]
            # image = image.permute(1, 2, 0).numpy()
            # image = cv2.cvtColor(image, cv2.COLOR_RGB2BGR)
            # image *= 255
            # image = image.astype(np.uint8)
            # cv2.imwrite(str(i) + 'color.png', image)
            # masks = target['masks']
            # for mi, m in enumerate(masks):
            #     img = m.numpy()
            #     img *= 255
            #     print(np.max(img), np.min(img))
            #     img = img.astype(np.uint8)
            #     cv2.imwrite(str(i+mi) + 'mask.png', img)
            # exit(5)

            prediction = model([image.to(device)])
            print(prediction[0]['scores'])
            pred_mask = np.zeros((224, 224), dtype=np.uint8)
            for idx, mask in enumerate(prediction[0]['masks']):
                if prediction[0]['scores'][idx] > 0.9:
                    # if prediction[0]['scores'][idx] > 0.75:
                    img1 = mask[0].mul(255).byte().cpu().numpy()
                    img1[img1 > 80] = 255
                    img1[img1 <= 80] = 0
                    pred_mask[img1 > 80] = 255 - idx * 5

                    img1 = Image.fromarray(img1)
                    img1.save(str(count) + "-" + str(idx) + 'mask.png')
            for idx, mask in enumerate(targets[i]['masks']):
                img2 = Image.fromarray(mask.mul(255).byte().cpu().numpy())
                img2.save(str(count) + "-" + str(idx) + 'target.png')
            img0 = Image.fromarray(image.mul(255).permute(1, 2, 0).byte().numpy())
            img0.save(str(count) + "-" + str(idx) + 'ori.png')
            img0 = Image.fromarray(pred_mask)
            img0.save(str(count) + "-" + str(idx) + 'pred.png')
            count += 1
            # exit()

def main(args):
    data_loader, data_loader_test = get_data_loader(args.dataset_root, args.batch_size, args.is_real)

    device = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')

    model = get_model_instance_segmentation(2)
    if args.resume:
        model.load_state_dict(torch.load(os.path.join(args.dataset_root, "maskrcnn.pth")))
    model = model.to(device)

    # construct an optimizer
    params = [p for p in model.parameters() if p.requires_grad]
    optimizer = torch.optim.SGD(params, lr=1e-3, momentum=0.9, weight_decay=1e-4)

    # and a learning rate scheduler which decreases the learning rate by 10x every 1 epochs
    lr_scheduler = torch.optim.lr_scheduler.MultiStepLR(optimizer, milestones=[10, 25], gamma=0.1)

    if args.test:
        test(model, data_loader_test, device=device)
    else:
        for epoch in range(args.epochs):
            # train for one epoch, printing every 10 iterations
            train_one_epoch(model, optimizer, data_loader, device, epoch, print_freq=50, resume=args.resume)
            torch.save(model.state_dict(), os.path.join(args.dataset_root, "maskrcnn.pth"))
            # update the learning rate
            lr_scheduler.step()
            # evaluate on the test dataset
            evaluate(model, data_loader_test, device=device)

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='Train foreground')

    parser.add_argument('--dataset_root', dest='dataset_root', action='store', help='Enter the path to the dataset')
    parser.add_argument('--is_real', dest='is_real', action='store_true', default=False, help='')
    parser.add_argument('--epochs', dest='epochs', action='store', type=int, default=20, help='Enter the epoch for training')
    parser.add_argument('--batch_size', dest='batch_size', action='store', type=int,
                        default=4, help='Enter the batchsize for training and testing')
    parser.add_argument('--test', dest='test', action='store_true', default=False, help='Testing and visualizing')
    parser.add_argument('--resume', dest='resume', action='store_true', default=False, help='Enter the path to the dataset')

    args = parser.parse_args()
    if args.test:
        args.resume = True

    main(args)

max_stars_repo_name: mendelmaker/dipn

max_issues_repo_name: mendelmaker/dipn

max_forks_repo_name: mendelmaker/dipn


text:
# -*- coding: utf-8 -*-
"""
Created on 17 July 2019

@author: Jonas Karthein
@contact: jonas.karthein@cern.ch
@license: MIT license
"""

import sys, os, glob, json, csv
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.optimize

class Freq_ratio():
    # class to calculate the frequency ratio for a mass measurement (MM), optimized for ISOLTRAP PI-ICR data
    def __init__(self):
        self.degree = 3
        self.param_start = [1]*(self.degree+2)
        self.param_result = []
        self.y1df = {}
        self.y1names = []
        self.y1 = []
        self.y1err = []
        self.y2df = {}
        self.y2names = []
        self.y2 = []
        self.y2err = []
        self.x1df = {}
        self.x1helpdf = {}
        self.x1names = []
        self.x1 = []
        self.x2df = {}
        self.x2helpdf = {}
        self.x2names = []
        self.x2 = []
        self.x2folder = ''      # time stamps isotope 2
        self.isotopes = []
        self.dirname = ''
        self.run_folder = ''
        self.ratio = 0
        self.ratio_unc = 0
        self.red_chi_sq = 1
        self.covar = []
        self.correlation = 0
        self.ratio_dict = {}
        self.weighted_avg = 0
        self.total_unc = 0
        self.birge_ratio = 0
        self.merged_dict = {}
        self.lin_extrapolation = False
        self.mode = ''

        self.test = False

    def ratio_sim_fit(self, isotopes, cyc_freq_1_time_delta, cyc_freq_1, cyc_freq_1_unc, cyc_freq_2_time_delta, cyc_freq_2, cyc_freq_2_unc):
        """"
        Main function perform simultaneous polynomial fits for all for the data set possible polynomial degrees. Calculates red. chi square for each fit and returns the fit results for the one with a red. chi square clostest to 1.

        :param isotopes: list of str of isotopes
        :param cyc_freq_1_time_delta: list of time deltas in minutes to first time for isotope 1
        :param cyc_freq_1: list of cyclotron frequencies for isotope 1
        :param cyc_freq_1_unc: list of cyclotron frequency uncertainties for isotope 1
        :param cyc_freq_2_time_delta: list of time deltas in minutes to first time for isotope 2
        :param cyc_freq_2: list of cyclotron frequencies for isotope 2
        :param cyc_freq_2_unc: list of cyclotron frequency uncertainties for isotope 2

        :return param1: dict containing the fit information for isotope 1
        :return param2: dict containing the fit information for isotope 2
        :return ratio: frequency ratio deriving from best fit
        :return ratio_unc: uncertainty of frequency ratio deriving from best fit
        :return red_chi_sq: reduced chi square for best fit
        """
        self.x1 = cyc_freq_1_time_delta
        self.x2 = cyc_freq_2_time_delta
        self.y1 = cyc_freq_1
        self.y1err = cyc_freq_1_unc
        self.y2 = cyc_freq_2
        self.y2err = cyc_freq_2_unc
        self.isotopes = isotopes
        temp_result_dict = {}
        temp_result_dict2 = {}

        for deg in range(2,min(len(self.y1),len(self.y2))+2):     # calculates ratio for many polynomial degrees, as much as data allows
            try:
                self.degree = deg
                self.get_start_param()
                if len(self.y2) > 1:    # make sure there are at least two points to fit
                    print('Poly-degree: ', self.degree)
                    self.fit()
                    self.red_chi_sq_calc()
                    print('\nRed.Chi.Sq.: ', self.red_chi_sq)
                    self.covar = self.covar * self.red_chi_sq
                    self.param_result_err = np.array([np.sqrt(np.absolute(self.covar[i, i])) for i in range(len(self.param_result))])
                    self.correlation = self.covar[len(self.param_result)-2, len(self.param_result)-1]       # calculates the corellation
                    print('Corellation: ', self.correlation)
                else:       # if less than two points in data set, then a linear extrapolation is performed
                    self.lin_extrapolation = True
                    self.fit()
                    self.param_result_err = np.sqrt(np.diag(self.covar))
                # calculates the frequency ratio from the fit information
                self.ratio_calc()
                self.lin_extrapolation = False

                temp_result_dict['{}'.format(self.degree)] = {'param1': [float(i) for i in list(self.param_result[:-1])],
                                                              'param2': [float(i)/float(self.param_result[-1]) for i in list(self.param_result[:-1])],
                                                              'ratio': self.ratio,
                                                              'ratio_unc': self.ratio_unc,
                                                              'red_chi_sq': self.red_chi_sq}
                temp_result_dict2['{}'.format(self.degree)] = abs(self.red_chi_sq-1.0)      # allowing to find the red_chi_sq closest to 1 by a minimizing function
            except:     # if one fit fails, program stops
                break

        try:
            best_chi_sq = min(temp_result_dict2, key=temp_result_dict2.get)
        except:
            best_chi_sq = -1

        return(temp_result_dict[best_chi_sq]['param1'], temp_result_dict[best_chi_sq]['param2'], temp_result_dict[best_chi_sq]['ratio'], temp_result_dict[best_chi_sq]['ratio_unc'], temp_result_dict[best_chi_sq]['red_chi_sq'])

    def get_start_param(self):
        """Sets starting parameter. If fitting fails, try to adjust the starting parameter"""
        self.param_start = [1]*(self.degree+1)
        self.param_start[-2] = self.y1[0]   # linear offset of fit
        self.param_start[-3] = 0.0001           # helps p=2 fit to converge
        #     self.param_start[-3] = 0.1           # helps p=2 fit to converge
        #     print(self.y2, self.y1)
        print(self.param_start)

    def get_data(self):
        """Get frequency and timestamps"""
        # ___________________________________________
        # get frequency data for isotope 1
        self.y1df = pd.read_csv(self.y1folder)
        self.y1names = self.y1df.iloc[:, 0].tolist()
        self.y1df = self.y1df.sort_values('File')
        self.y1df = self.y1df.reset_index(drop=True)
        self.y1 = self.y1df.iloc[:, 1].tolist()
        self.y1err = self.y1df.iloc[:, 2].tolist()

        # ___________________________________________
        # get frequency data for isotope 2
        self.y2df = pd.read_csv(self.y2folder)
        self.y2df = self.y2df.sort_values('File')
        self.y2df = self.y2df.reset_index(drop=True)
        self.y2names = self.y2df.iloc[:, 0].tolist()
        self.y2 = self.y2df.iloc[:, 1].tolist()
        self.y2err = self.y2df.iloc[:, 2].tolist()
        # ___________________________________________
        # get time stamps for each file for isotope 1
        self.x1df = pd.read_csv(self.x1folder,
                                names=['names', 'date_begin', 'time_begin', 'date_end', 'time_end'],
                                parse_dates = [[1,2], [3,4]])
        self.x1df = self.x1df.sort_values('names')
        self.x1df = self.x1df.reset_index(drop=True)
        self.x1df['timedelta'] = self.x1df['date_end_time_end'] - self.x1df['date_begin_time_begin']
        self.x1names = self.x1df['names'].tolist()
        # calculate the median time stamp for each file also considering the count rate
        self.x1df['median_counts'] = self.get_timestamp_position(self.x1helpfolder, self.x1names)
        self.x1df['median_time'] = (self.x1df['date_begin_time_begin'] +
                                    self.x1df['timedelta'] * self.x1df['median_counts'])

        # ___________________________________________
        # get time stamps for each file for isotope 2
        self.x2df = pd.read_csv(self.x2folder,
                                names=['names', 'date_begin', 'time_begin', 'date_end', 'time_end'],
                                parse_dates = [[1,2], [3,4]])
        self.x2df = self.x2df.sort_values('names')
        self.x2df = self.x2df.reset_index(drop=True)

        self.x2df['timedelta'] = self.x2df['date_end_time_end'] - self.x2df['date_begin_time_begin']
        self.x2names = self.x2df['names'].tolist()
        # calculate the median time stamp for each file also considering the count rate
        self.x2df['median_counts'] = self.get_timestamp_position(self.x2helpfolder, self.x2names)
        self.x2df['median_time'] = (self.x2df['date_begin_time_begin'] +
                                    self.x2df['timedelta'] * self.x2df['median_counts'])
        # the x values for the frequency information is converted into minutes after(/before) the first
        # x1 frequency point started. This had to be changed since the fit wasn't very stable for x values
        # having the correct time converted to seconds since the values were very large (e.g. 1498777404).
        # The following statement reduces the x value to minutes in the 2-3 digit range. Important here is,
        # that x1 and x2 are normalized to the same starting time of MIN(x1[0], x2[0])
        self.x1 = [round(i.total_seconds()/60, 2) for i in (self.x1df['median_time'] - min([self.x1df['median_time'][0], self.x2df['median_time'][0]])).tolist()]
        self.x2 = [round(i.total_seconds()/60, 2) for i in (self.x2df['median_time'] - min([self.x1df['median_time'][0], self.x2df['median_time'][0]])).tolist()]

    def get_timestamp_position(self, folder, x1_or_x2_names):
        """calculates the mean position of ions to arrive during the measurement time --> if beam gate was changed during the run, the average timestamp is not in the middle of the measurement"""
        median_counts = []
        for i in x1_or_x2_names:
            tmp = folder+i+'_p1_spot_positions.csv'
            self.x1helpdf = pd.read_csv(tmp, names=['event', 'x', 'y', 'time_stamp_us'])
            # we chose the median over the mean to be less sensitive to outliers in case of a drop in count rate
            median_counts.append(float(self.x1helpdf.iloc[:, 0].median() / self.x1helpdf.iloc[:, 0].max()))
        return(median_counts)

    def fit(self):
        """Actual simultaneous polynomial fitting routine. Executed five times for stability reasons."""
        if self.lin_extrapolation == False:
            # fit ; timestampt value is reduced since resulting numbers were too large
            iter1, covar, infodict, mesg, ier = scipy.optimize.leastsq(self.simultaneous_residual, self.param_start,
                                                            args=(np.array(self.x1),
                                                                  np.array(self.x2),
                                                                  np.array(self.y1),
                                                                  np.array(self.y2),
                                                                  np.array(self.y1err),
                                                                  np.array(self.y2err)), full_output=True)#, maxfev=1000000
            iter2, covar, infodict, mesg, ier = scipy.optimize.leastsq(self.simultaneous_residual, iter1,
                                                            args=(np.array(self.x1),
                                                                  np.array(self.x2),
                                                                  np.array(self.y1),
                                                                  np.array(self.y2),
                                                                  np.array(self.y1err),
                                                                  np.array(self.y2err)), full_output=True)
            iter3, covar, infodict, mesg, ier = scipy.optimize.leastsq(self.simultaneous_residual, iter2,
                                                            args=(np.array(self.x1),
                                                                  np.array(self.x2),
                                                                  np.array(self.y1),
                                                                  np.array(self.y2),
                                                                  np.array(self.y1err),
                                                                  np.array(self.y2err)), full_output=True)
            iter4, covar, infodict, mesg, ier = scipy.optimize.leastsq(self.simultaneous_residual, iter3,
                                                            args=(np.array(self.x1),
                                                                  np.array(self.x2),
                                                                  np.array(self.y1),
                                                                  np.array(self.y2),
                                                                  np.array(self.y1err),
                                                                  np.array(self.y2err)), full_output=True)
            self.param_result, self.covar, infodict, mesg, ier = scipy.optimize.leastsq(self.simultaneous_residual, iter4,
                                                            args=(np.array(self.x1),
                                                                  np.array(self.x2),
                                                                  np.array(self.y1),
                                                                  np.array(self.y2),
                                                                  np.array(self.y1err),
                                                                  np.array(self.y2err)), full_output=True)
            if ier not in [1,2,3,4]:
                sys.exit('Fit did not converge')

            # print self.param_result
            # print self.covar

        else:       # linear extrapolation
            init_vals = [((max(self.y1)-min(self.y1))
                          /(max(self.x1)-min(self.x1))),self.y1[0]]
            fit_1, pcov = scipy.optimize.curve_fit(self.linear,
                                       self.x1, self.y1, sigma=self.y1err,
                                       absolute_sigma=True,
                                       p0=init_vals)
            fit_2, pcov = scipy.optimize.curve_fit(self.linear,
                                       self.x1, self.y1, sigma=self.y1err,
                                       absolute_sigma=True,
                                       p0=fit_1)
            fit_3, pcov = scipy.optimize.curve_fit(self.linear,
                                       self.x1, self.y1, sigma=self.y1err,
                                       absolute_sigma=True,
                                       p0=fit_2)
            fit_4, pcov = scipy.optimize.curve_fit(self.linear,
                                       self.x1, self.y1, sigma=self.y1err,
                                       absolute_sigma=True,
                                       p0=fit_3)
            self.param_result, self.covar = scipy.optimize.curve_fit(self.linear,
                                       self.x1, self.y1, sigma=self.y1err,
                                       absolute_sigma=True,
                                       p0=fit_4)
            print(self.param_result)

    def fit_function(self, x, param):
        """
        Returns a polynomial function value of degree len(param)-1 evaluated at x

        :param x: time delta value
        :param param: list of polynomial parameters

        :return: polynomial function value of degree len(param)-1 evaluated at x
        """
        return(np.polyval(list(param), x))

    def linear(self, x, m, n):
        """
        Returns a linear function value evaluated at x

        :param x: time delta value
        :param m: slope
        :param n: frequency axis intercept

        :return: linear function value evaluated at x
        """
        return(m*x+n)

    def residual(self, p, x, y, yerr):
        """
        Returns residual, which square is minimized later.

        :param p: list of polynomial parameters
        :param x: list of time delta values
        :param y: list of cyclotron frequency values
        :param yerr: list of cyclotron frequency uncertainty values

        :return: residual, which square is minimized later
        """
        return((self.fit_function(x, p) - y )/yerr)

    def simultaneous_residual(self, param_start, x1, x2, y1, y2, y1err, y2err):
        """
        Returns residual for both data sets.

        :param param_start: list of polynomial start parameters for fit
        :param x1: list of time delta values of isotope 1
        :param y1: list of cyclotron frequency values of isotope 1
        :param y1err: list of cyclotron frequency uncertainty values of isotope 1
        :param x2: list of time delta values of isotope 2
        :param y2: list of cyclotron frequency values of isotope 2
        :param y2err: list of cyclotron frequency uncertainty values of isotope 2

        :return: total residual
        """

        p_s = param_start
        p1 = p_s[:-1]  # Poly1 = R*Poly2
        p2 = [float(i)/p_s[-1] for i in p_s[:-1]]  # Poly1 = R*Poly2
        res1 = self.residual(p1, x1, y1, y1err)
        res2 = self.residual(p2, x2, y2, y2err)
        return(np.concatenate((res1, res2)))

    def final_residual_1(self, param_start, x1, y1, y1err):
        """
        Returns final residual for isotope 1.

        :param param_start: list of final polynomial parameters of fit
        :param x1: list of time delta values of isotope 1
        :param y1: list of cyclotron frequency values of isotope 1
        :param y1err: list of cyclotron frequency uncertainty values of isotope 1

        :return: final residual for isotope 1
        """
        if self.test:
            p_s = [param_start['param{}'.format(i)] for i in [len(param_start) - 1 - j for j in range(len(param_start))]]
        else:
            p_s = param_start
        p1 = p_s[:-1]  # Poly1 = R*Poly2
        res1 = self.residual(p1, x1, y1, y1err)
        return(res1)

    def final_residual_2(self, param_start, x2, y2, y2err):
        """
        Returns final residual for isotope 2.

        :param param_start: list of final polynomial parameters of fit
        :param x2: list of time delta values of isotope 2
        :param y2: list of cyclotron frequency values of isotope 2
        :param y2err: list of cyclotron frequency uncertainty values of isotope 2

        :return: final residual for isotope 2
        """
        if self.test:
            p_s = [param_start['param{}'.format(i)] for i in [len(param_start) - 1 - j for j in range(len(param_start))]]
        else:
            p_s = param_start
        p2 = [float(i)/p_s[-1] for i in p_s[:-1]]  # Poly1 = R*Poly2
        res2 = self.residual(p2, x2, y2, y2err)
        return(res2)

    def red_chi_sq_calc(self):
        """Calculates reduced chi square for fit."""
        sum_res_sq_1 = 0
        sum_res_sq_2 = 0
        for i in range(len(self.x1)):
            sum_res_sq_1 += (self.final_residual_1(self.param_result, self.x1[i], self.y1[i], self.y1err[i]))**2
        for i in range(len(self.x2)):
            sum_res_sq_2 += (self.final_residual_2(self.param_result, self.x2[i], self.y2[i], self.y2err[i]))**2
        self.red_chi_sq = (sum_res_sq_1+sum_res_sq_2)/(len(self.x1)+len(self.x2)-self.degree-1)

    def ratio_calc(self):
        """Calculates the ratio in the middle of the range of interest."""
        if self.lin_extrapolation == False:
            self.ratio = self.param_result[-1]
            self.ratio_unc = self.param_result_err[-1]
        else:
            fit_at_y2 = self.linear(self.x2[0], self.param_result[0], self.param_result[1])
            self.ratio = fit_at_y2 / self.y2[0]

            unc_lin = np.sqrt( ( self.param_result_err[0] * self.x2[0] )**2
                               + self.param_result_err[1] ** 2 )
            self.ratio_unc = np.sqrt( ( unc_lin / self.y2[0] )**2
                                      + ( -( fit_at_y2 * self.y2err[0] )
                                          / (self.y2[0])**2 )**2 )
            self.degree = 'linear'
            print('Degree:              ', self.degree)

        print('Ratio fit parameter: ', self.ratio, '+/-', self.ratio_unc)

max_stars_repo_name: jonas-ka/pi-icr-analysis

max_issues_repo_name: jonas-ka/pi-icr-analysis

max_forks_repo_name: jonas-ka/pi-icr-analysis

reason: import numpy,import scipy


text:
import gzip

import numpy as np

# Download dataset
IMAGE_SIZE = 28
NUM_CHANNELS = 1
PIXEL_DEPTH = 255
BATCH_SIZE = 128

def extract_data(filename, num_images):
    """Extract the image into a 4D tensor [image index, y, x, channels].

    Values are rescaled from [0, 255] down to [-0.5, 0.5].
    """
    print(f"Extracting {filename}")
    with gzip.open(filename) as bytestream:
        bytestream.read(16)
        buf = bytestream.read(
            IMAGE_SIZE *
            IMAGE_SIZE *
            num_images *
            NUM_CHANNELS)
        data = np.frombuffer(buf, dtype=np.uint8).astype(np.float32)
        data = (data - (PIXEL_DEPTH / 2.0)) / PIXEL_DEPTH
        data = data.reshape(num_images, IMAGE_SIZE, IMAGE_SIZE, NUM_CHANNELS)
        return data

def extract_labels(filename, num_images):
    """Extract the labels into a vector of int64 label IDs."""
    print(f"Extracting {filename}")
    with gzip.open(filename) as bytestream:
        bytestream.read(8)
        buf = bytestream.read(1 * num_images)
        labels = np.frombuffer(buf, dtype=np.uint8).astype(np.int64)
    return labels

def print_activations(t):
    print(t.op.name, f" {t.get_shape().as_list()}")

max_stars_repo_name: Lornatang/tensorflow-tutorial

max_issues_repo_name: Lornatang/tensorflow-tutorial

max_forks_repo_name: Lornatang/tensorflow-tutorial


text:
#! /usr/bin/env python

"""Wrappers for photutils' finder functions `find_peaks` and `daofind`
(technically `photutils/detection/findstars.daofind`).
I add on functionality to print output Table to a text file, similar to
IRAF's out files. The output file differs from IRAF, however, in that
it is column-based. It does include a parameter settings list in the
header comments.

Author:

    C.M. Gosmeyer, Feb. 2016

References:

About IRAF daofind:
http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?daofind

About photutils daofind:
http://photutils.readthedocs.org/en/latest/_modules/photutils/detection/findstars.html#daofind

About photutils find_peaks:
http://photutils.readthedocs.org/en/latest/api/photutils.detection.find_peaks.html#photutils.detection.find_peaks

Example for use:
http://photutils.readthedocs.org/en/latest/photutils/getting_started.html

"""
from __future__ import print_function

import multiprocessing
import os
import sys
import photutils
import numpy as np
from astropy.io import fits
from astropy.stats import mad_std
from astropy.table import Column

from photutils_plus.phot_tools import read_in_photfile
from photutils_plus.phot_tools import write_out_photfile

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

def daofind(imagename, ext=1, outname='default', threshold=3., fwhm=1.5, \
            ratio=1.0, theta=0.0, sigma_radius=1.5, sharplo=0.2, \
            sharphi=1.0, roundlo=-1.0, roundhi=1.0, sky=0.0, \
            exclude_border=True, use_bkgrd=True):
    """ Extends `photutils.daofind`, so that outputs a coordinate file
    with my formatting from `photutils_plus.phot_tools.write_out_photfile`.
    (Later, if I'm clever, use inheritance.)

    Parameters
    ----------
    imagename : string
        Name of the FITS file.
    ext : int
        The extension of the FITS file to read.
    outname : string
        Name of the output coordinate file. If "default", becomes,
        "<imagename>.coo."
    threshold : int
        Threshold for search. Default of "3."
    fwhm : float
        The full width half max. Default of "1.5."
    ratio : float
        Default of "1.0."
    theta : float
        Default of "0.0."
    sigma_radius : float
        Default of 1.5."
    sharplo : float
        Default of "0.2."
    sharphi : float
        Default of "1.0."
    roundlo : float
        Default of "-1.0."
    roundhi : float
        Default of "1.0."
    sky : float
        Default of "0.0."
    exclude_border : {True, Flase}
        When True, masks out a 10-pixel border on image from search.
    use_bkgrd : {True, False}
        Setting on will multiply the background by the threshold to
        obtain a new threshold value.

    Returns
    -------
    coo_tab : astropy.Table
        Table containing the coordinates.

    """
    # Fetch function metadata.
    current_params = locals()
    func_name = sys._getframe().f_code.co_name #daofind.__name__

    # Read in FITS file.
    hdulist = fits.open(imagename)
    data = hdulist[ext].data

    # Mask 10 pixels from border. I find the `exclude_border` in
    # photutils.daofind inadiquate.
    if exclude_border:
        data[:,0:10] = -99999.0
        data[:,-10:] = -99999.0
        data[0:10,:] = -99999.0
        data[-10:,:] = -99999.0

    # Get the background.
    if use_bkgrd:
        bkg_sigma = mad_std(data)
        threshold = threshold * bkg_sigma

    coo_tab = photutils.daofind(data=data,
                                threshold=threshold,
                                fwhm=fwhm,
                                ratio=ratio,
                                theta=theta,
                                sigma_radius=sigma_radius,
                                sharplo=sharplo,
                                sharphi=sharphi,
                                roundlo=roundlo,
                                roundhi=roundhi,
                                sky=sky,
                                exclude_border=exclude_border)

    # ??Insert background value into coo_tab??

    # Create basename for out file.
    if outname == 'default':
        baseoutname = imagename + '.coo'
    else:
        baseoutname = outname

    # Check whether file already exists. If yes, append number.
    fileoutname = baseoutname
    i=1
    while os.path.isfile(fileoutname):
        fileoutname = baseoutname
        fileoutname += ('.' + str(i))
        i += 1

    write_out_photfile(fileoutname, coo_tab, current_params, func_name)

    return coo_tab

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

def find_peaks(imagename, ext=1, outname='default', threshold=40.0,
               box_size=24.0, border_width=5.0):
#, footprint=footprint, mask=mask,
#               border_width=border_width, npeaks=npeaks, subpixel=subpixel,
#               error=error, wcs=wcs):
    """ Extends `photutils.find_peaks`.

    Parameters
    ----------
    imagename : string
        Name of the FITS file.
    ext : int
        The extension of the FITS file to read.
    outname : string
        Name of the output coordinate file. If "default", becomes,
        "<imagename>.coo."
    threshold : float
        Default of "40.0."
    box_size : float
        Default of "24.0."
    border_width : float
        Default of "5.0."

    Returns
    -------
    coo_tab : astropy.Table
        Table containing the coordinates.

    References
    ----------
    http://photutils.readthedocs.org/en/latest/api/photutils.detection.find_peaks.html#photutils.detection.find_peaks

    """
    # Fetch function metadata.
    current_params = locals()
    func_name = sys._getframe().f_code.co_name

    # Read in FITS file.
    hdulist = fits.open(imagename)
    data = hdulist[ext].data
    hdulist.close()

    coo_tab = photutils.find_peaks(data=data, \
                                   threshold=threshold, \
                                   box_size=box_size, \
                                   border_width=border_width)
                                   #footprint=footprint, \
                                   #mask=mask, \
                                   #npeaks=npeaks, \
                                   #subpixel=subpixel, \
                                   #error=error, \
                                   #wcs=wcs)

    # Create basename for out file.
    if outname == 'default':
        baseoutname = imagename + '.coo'
    else:
        baseoutname = outname

    # Check whether file already exists. If yes, append number.
    fileoutname = baseoutname
    i=1
    while os.path.isfile(fileoutname):
        fileoutname = baseoutname
        fileoutname += ('.' + str(i))
        i += 1

    write_out_photfile(fileoutname, coo_tab, current_params, func_name)

    return coo_tab

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

def multiprocess_find(func, filenames):
    """ Wrapper for multiprocessing a finder function over a list of
    FITS files.

    Parameters
    ----------
    func : find function
        Either :func:`daofind` or :func:`find_peaks`.
    filenames : list of strings
        The names of the FITS files to be multiprocssed with a
        find function.
    """

    pool = multiprocessing.Pool(processes=SETTINGS['cores'])
    pool.map(func, filenames)
    pool.close()
    pool.join()

    print("Processes complete.")

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

def test_daofind():
    """
    """
    imagename = 'n6791_11924_55106_F606W_360_NOCI_F000_c0_1.fits' #'icqu04rgq_flc.fits' #'id1l02g2q_flt.clean.fits'
    ext = 0
    daofind(imagename, ext=ext, threshold=1050., fwhm=2.5)

def test_find_peaks():
    """
    """
    imagename = 'id1l02g2q_flt.clean.fits'  #'icqu04rgq_flc.fits'
    ext = 0
    find_peaks(imagename, ext=ext, threshold=40., box_size=24.)

#-------------------------------------------------------------------------------#
# The Main, for testing.
#-------------------------------------------------------------------------------#

if __name__=='__main__':
    test_daofind()
    #data, header = read_in_photfile('id1l02g2q_flt.clean.fits.coo.1')
    #print data
    #print header

reason: import numpy,from astropy