text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# File            : live_visualisation.py
# Author          : Duy Anh Pham <duyanh.y4n.pham@gmail.com>
# Date            : 10.04.2020
# Last Modified By: Duy Anh Pham <duyanh.y4n.pham@gmail.com>

from djitellopy.realtime_plot.RealtimePlotter import *
import redis
import numpy as np
import traceback
import matplotlib

# define data to get from db
# sensorMeshList = ['baro', 'h', 'tof', 'runtime']
# row = len(sensorMeshList)
data_len = 300
plot_update_interval = 0.005
datasource = redis.StrictRedis(host='localhost', port=6379, db=0)

plt.figure()

baro_axes = plt.subplot(3, 1, 1)
plt.title('tello_edu sensors')
baro_data_list = ['baro', 'runtime']
baro_ylim = [-47, -57]
baro_option = DataplotOption.TIMESTAMP_CUSTOM
baro_dataplot = DataPlot(2, data_len, option=baro_option)
baro_plot = RealtimePlotter(baro_dataplot)
baro_plot.config_plots(baro_axes, y_labels=baro_data_list, ylim=baro_ylim)
baro_plot.axes.set_xlabel('time in ms')
baro_plot.axes.set_ylabel('barometer in cmHg')

tof_axes = plt.subplot(3, 1, 2)
tof_data_list = ['tof', 'runtime']
tof_ylim = [-10, 500]
tof_option = DataplotOption.TIMESTAMP_CUSTOM
tof_dataplot = DataPlot(2, data_len, option=tof_option)
tof_plot = RealtimePlotter(tof_dataplot)
tof_plot.config_plots(tof_axes, y_labels=tof_data_list, ylim=tof_ylim)
tof_plot.axes.set_xlabel('time in ms')
tof_plot.axes.set_ylabel('vertical distance in cm')

h_axes = plt.subplot(3, 1, 3)
h_ylim = [-50, 300]
h_data_list = ['h', 'runtime']
h_option = DataplotOption.TIMESTAMP_CUSTOM
h_dataplot = DataPlot(2, data_len, option=h_option)
h_plot = RealtimePlotter(h_dataplot)
h_plot.config_plots(h_axes, y_labels=h_data_list, ylim=h_ylim)
h_plot.axes.set_xlabel('time in ms')
tof_plot.axes.set_ylabel('height in cm')

if __name__ == "__main__":
    while True:
        # get new data from database and plot
        # baro
        baro_plot.dataplot.clear_data_regs()
        new_data = []
        for sensor in baro_data_list:
            new_sensor_data = datasource.lrange(sensor, 0, data_len)
            # reverse, bc first element is the newest (not the oldest like deque)
            new_sensor_data.reverse()
            new_data.append(new_sensor_data)
        try:
            baro_y = np.array(new_data[:-1], dtype=np.float)
            baro_x = np.array(new_data[-1], dtype=np.int64)
            baro_plot.dataplot.append(
                y=baro_y, x=baro_x, single=False)
            baro_plot.plot_data()
        except Exception as e:
            print(e)
        # tof
        tof_plot.dataplot.clear_data_regs()
        new_data = []
        for sensor in tof_data_list:
            new_sensor_data = datasource.lrange(sensor, 0, data_len)
            # reverse, bc first element is the newest (not the oldest like deque)
            new_sensor_data.reverse()
            new_data.append(new_sensor_data)
        try:
            tof_y = np.array(new_data[:-1], dtype=np.float)
            tof_x = np.array(new_data[-1], dtype=np.int64)
            tof_plot.dataplot.append(
                y=tof_y, x=tof_x, single=False)
            tof_plot.plot_data()
        except Exception as e:
            print(e)
        # height
        h_plot.dataplot.clear_data_regs()
        new_data = []
        for sensor in h_data_list:
            new_sensor_data = datasource.lrange(sensor, 0, data_len)
            # reverse, bc first element is the newest (not the oldest like deque)
            new_sensor_data.reverse()
            new_data.append(new_sensor_data)
        try:
            h_y = np.array(new_data[:-1], dtype=np.float)
            h_x = np.array(new_data[-1], dtype=np.int64)
            h_plot.dataplot.append(
                y=h_y, x=h_x, single=False)
            h_plot.plot_data()
        except Exception as e:
            print(e)

        plt.pause(plot_update_interval)
    input("Exit(press any key)?")


text:
"""Non-Linear Independent Components Estimation with Tensorflow

This script is implementation of Non-Linear Independent components estimation algorithm with Tensorflow

Command- python nice.py

This file contains the following class & functions:
	* ScalingLayer - This class rescales the output of coupling layers.
    * CouplingLayer - This class implement additive coupling layer.
    * NICE - This class train the model using NICE algorithm and generate new images.
"""

# Load Libraries
import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = "2"

# Load Libraries
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
from tensorflow.keras.datasets import mnist
import numpy as np
import matplotlib.pyplot as plt

class ScalingLayer(layers.Layer):
	""" This method define the scaling layer with normal distribution.

	Parameters
	----------
	input_dim : int
	- Number of Pixels - width * height
	"""
	def __init__(self,input_dim):
		super(ScalingLayer,self).__init__()
		self.scaling_layer = tf.Variable(tf.random.normal([input_dim]),trainable=True)

	""" This method rescale the final out of coupling layers.

	Parameters
	----------
	x : Tensor Array
	- Tensors to be rescale.
	"""
	def call(self,x):
		return x*tf.transpose(tf.math.exp(self.scaling_layer))

	""" This method is inverse implementation of the call method.

	Parameters
	----------
	z : Tensor Array
	    - Tensors to be rescale.
	"""
	def inverse(self,z):
	    return z*tf.transpose(tf.math.exp(-self.scaling_layer))

class CouplingLayer(layers.Layer):

	""" This method define the neural network to learn function F

    Parameters
    ----------
    input_dim : int
        - Number of Pixels - width * height
    hidden_dim : int
        - Number of Units per hidden layer
    partition : String
        - Odd or even indexing scheme
    num_layers : int
        - Number of hidden layer

    """

	def __init__(self,input_dim, hidden_dim, partition, num_hidden_layers):

		super(CouplingLayer, self).__init__()

		assert input_dim % 2 == 0

		self.partition = partition # Posssible values - odd or even

		# Define Neural Network

		inputs = keras.Input(shape = (input_dim))

		x = inputs
		for layer in range(num_hidden_layers):
			x = layers.Dense(hidden_dim)(x)
			x = layers.LeakyReLU()(x)

		outputs = layers.Dense(input_dim, activation='linear')(x)
		self.m = keras.Model( inputs = inputs, outputs = outputs)

		# Define lambda functions to get odd or even indexed data.
		if self.partition == 'even':
			self._first = lambda xs: xs[:,0::2]
			self._second = lambda xs: xs[:,1::2]
		else:
			self._first = lambda xs: xs[:,1::2]
			self._second = lambda xs: xs[:,0::2]

	""" This method implement additive coupling layer operations.
	y1 = x1
	y2 = x2 + m(x1)

    Parameters
    ----------
    inputs : Tensor Array
        Input Tensors.

    Returns
        -------
        Tensor Array
            Latent representation.
    """

	def call(self,inputs):

		# Split input into two parts x1 and x2
		x1 = self._first(inputs)
		x2 = self._second(inputs)

		# Inference latent representation
		y1 = x1
		y2 = x2 + self.m(x1)

		# Merge both y1 an y2 using interleave method. e.g y1 = [1,3,5,7,9] and y2 = [2,4,6,8,10] then y = [1,2,3,4,5,6,7,8,9,10]
		if self.partition == 'even':
		      y = tf.reshape(tf.concat([y1[...,tf.newaxis], y2[...,tf.newaxis]], axis=-1), [tf.shape(y1)[0],-1])
		else:
		      y = tf.reshape(tf.concat([y2[...,tf.newaxis], y1[...,tf.newaxis]], axis=-1), [tf.shape(y2)[0],-1])
		return y

	"""This function will generate new image using latent variable

    Parameters
    ----------
    latent_variable : Tensor Array
       . latent variable from logistic distribution.

    Returns
    -------
    List
        Pixels value.
    """
	def inverse(self,latent_variable):

		y1 = self._first(latent_variable)
		y2 = self._second(latent_variable)
		x1,x2 = y1, y2 - self.m(y1)

		if self.partition == 'even':
		      x = tf.reshape(tf.concat([x1[...,tf.newaxis], x2[...,tf.newaxis]], axis=-1), [tf.shape(x1)[0],-1])
		else:
		      x = tf.reshape(tf.concat([x2[...,tf.newaxis], x1[...,tf.newaxis]], axis=-1), [tf.shape(x2)[0],-1])

		return x

class NICE(keras.Model):

	""" This function define trainable model with multiple coupling layers and a scaling layer.

    Parameters
    ----------
    input_dim : int
        - Number of Pixels - width * height
    hidden_dim : int
        - Number of Units per hidden layer
    num_hidden_layers : int
        - Number of hidden layer

    """

	def __init__(self,input_dim,hidden_dim,num_hidden_layers):

		super(NICE,self).__init__()
		self.input_dim = input_dim
		# Coupling layer will have half dimension of the input data.
		half_dim = int(self.input_dim / 2)

		# Define 4 coupling layers.
		self.coupling_layer1 = CouplingLayer(input_dim=half_dim, hidden_dim = hidden_dim, partition='odd', num_hidden_layers=num_hidden_layers)
		self.coupling_layer2 = CouplingLayer(input_dim=half_dim, hidden_dim = hidden_dim, partition='even', num_hidden_layers=num_hidden_layers)
		self.coupling_layer3 = CouplingLayer(input_dim=half_dim, hidden_dim = hidden_dim, partition='odd', num_hidden_layers=num_hidden_layers)
		self.coupling_layer4 = CouplingLayer(input_dim=half_dim, hidden_dim = hidden_dim, partition='even', num_hidden_layers=num_hidden_layers)

		# Define scaling layer which rescaling the output for more weight variations.
		self.scaling_layer = ScalingLayer(self.input_dim)

	"""This function calculates the log likelihood.

    Parameters
    ----------
    z : Tensor Array
        Latent space representation.

    Returns
        -------
        List
            Log likelihood.
    """

	def log_likelihood(self, z):
		log_likelihood = tf.reduce_sum(-(tf.math.softplus(z) + tf.math.softplus(-z)), axis=1) + tf.reduce_sum(self.scaling_layer.scaling_layer)
		return log_likelihood

	""" This function passes input through coupling layer. Scale the output using scaling layer and return latent space z and loglikelihood

    Parameters
    ----------
    x : Tensor Array
        Input data from training samples.

    Returns
        -------
        Tensor Array
            Latent space representation.
        List
        	Log-likelihood for each image in the batch
    """

	def call(self, x):

		z = self.coupling_layer1(x)
		z = self.coupling_layer2(z)
		z = self.coupling_layer3(z)
		z = self.coupling_layer4(z)
		z = self.scaling_layer(z)

		log_likelihood = self.log_likelihood(z)
		return z,log_likelihood

	""" Generate sample data using latent space. This is exact inverse of the call method.

    Parameters
    ----------
    z : Tensor Array
        Latent Spaec.

    Returns
    -------
    Tensor Array
        Newly generated data
    """

	def inverse(self,z):
		x  = z
		x  = self.scaling_layer.inverse(x)
		x  = self.coupling_layer4.inverse(x)
		x  = self.coupling_layer3.inverse(x)
		x  = self.coupling_layer2.inverse(x)
		x  = self.coupling_layer1.inverse(x)

		return x

	""" Sample out new data from learnt probability distribution

    Parameters
    ----------
    num_samples : int
        - Number of sample images to generate.

    Returns
    -------
    List
    	- List of generated sample data.

    """

	def sample(self, num_samples):
		z = tf.random.uniform([num_samples, self.input_dim])
		z = tf.math.log(z) - tf.math.log(1.-z)
		return self.inverse(z)

""" This function show the generated images """
def show_samples():
    # Generate 100 sample digits
    ys = model.sample(21)
    plt.figure(figsize=(5, 5))
    plt.rcParams.update({
        "grid.color": "black",
        "figure.facecolor": "black",
        "figure.edgecolor": "black",
        "savefig.facecolor": "black",
        "savefig.edgecolor": "black"})
    new_images = tf.reshape(ys,(-1,28,28))
    for i in range(1,20):
      plt.subplot(2,10,i)
      plt.axis('off')
      plt.imshow(new_images[i],cmap="gray")
    plt.show()

# Basic Configuration
input_data_size = 784 # for mnist dataset. 28*28
num_epochs = 1
batch_size = 256
num_hidden_layer = 5
num_hidden_units = 1000

# Prepare training dataset
(x_train, y_train), (x_test, y_test) = mnist.load_data()
x_train = (x_train.reshape(-1,input_data_size) + tf.random.uniform([input_data_size])) / 256.0
dataloader = tf.data.Dataset.from_tensor_slices((x_train, y_train)).batch(batch_size)

# Custom Training Loop for the NICE Model
model = NICE(input_dim = input_data_size, hidden_dim = num_hidden_units, num_hidden_layers = num_hidden_layer)

opt = keras.optimizers.Adam(learning_rate=0.001,clipnorm=1.0)

# Iterate through training epoch
for i in range(num_epochs):

	# Iterate through dataloader
	for batch_id,(x,_) in enumerate(dataloader):

		# Forward Pass
		with tf.GradientTape() as tape:
			z, likelihood = model(x)
			log_likelihood_loss = tf.reduce_mean(likelihood)
			# Convert the maximum likelihood problem into negative of minimum likelihood problem
			loss = -log_likelihood_loss
		# Backward Pass
		gradient = tape.gradient(loss,model.trainable_weights)
		opt.apply_gradients(zip(gradient,model.trainable_weights))

	# Log the maxmimum likelihood after each training loop
	print('Epoch {} completed. Log Likelihood:{}'.format(i,log_likelihood_loss))
	# Generate 20 new images after 10 epochs.
	if i % 10 == 0 and i!=0:
		show_samples()

model.summary()
# Save model
model.save("nicegenerative")

ext: py

reason: import numpy


text:
"""
port octave code
"""
import numpy as np

class model():

    # state vector
    s1 = np.array([[0], [1]])
    s0 = np.array([[1], [0]])

    def __init__(self, n, num_layers, theta = None, theta_import = False):
        """
        n is the n x n matirx
        num layers is the number of cascading H Rz Cz gates
        theta should be a np array
        """
        self._n = n
        self._num_layers = num_layers
        # self.theta = theta #.npz file for theta parameters
        self._H = self.H()
        self._CZ = self.CZ()

        self.theta_list = self.random_theta(self._n, self._num_layers)

    @property
    def n(self):
        return self._n

    @n.setter
    def n(self, new_n):
        if new_n > 14:
            print("n cannot be greater than 14, n is set to 14")
            self._n = new_n
        elif new_n < 3:
            print("n cannot be smaller than 3, n is set to 3")
            self._n = 3

        else:
            self._n = new_n

    def get_num_layers(self):
        return self.num_layers

    def set_num_layers(self, new_num_layers):
        self.num_layers = new_num_layers

    def random_theta(self, n, num_layers):
        output = []
        for layer in range(num_layers):
            output.append(np.random.rand(n,n))

        return output

    def CZ(self):
        """
        n wire CZ, CZ acts on n wires only to modify
        last ancilla qubit
        matrix will be in the form
        1 0 0 0 ... 0
        0 1 0 0 ... 0
        0 0 1 0 ... 0
        .     .
        .       .
        .         .
                   1 0
        0 ...      0-1

        """
        I = np.identity(2**self._n)
        I[2**self._n - 1, 2**self._n -1] -= 2
        return I

    def H(self):
        """
        Returns a H gate
        """
        output = None
        H = (1/np.sqrt(2))*(np.matrix([[1,1],[1,-1]]))

        for x in range(self._n - 1):
            output = np.kron(H, H)

        return output

    def Rz(self, theta):
        """
        Returns R(theta) gate, or rotation matrix
        """
        return np.matrix([[exp(-1j * theta /2.0), 0],[0, exp(1j * theta/ 2.0)]])
        #TODO Kron the Rz theta gate with theta (n x n) input
        # Find a more efficient method

    def run_model(self, num_epoch, save_theta = True, location = "./theta.npz"):
        #for i in xrange(num_epoch):
        #    #print loss every 1000 epoch

        # RZ kron, if there are 4 lines, each layer should have 4 theta values.
        # I need to find some way of having some list of theta value and
        # converting to some Rz matrix and then kron those matrices

        graph = CZ * Rz1 * H * CZ
    def predict(theta):
        pass

if __name__ == "__main__":
    A = model(2,3)

    print(A._H)
    print(A.theta_list)


text:
'''
Created on Jul 26, 2021

@author: mhindle
'''
from pedigree.pedigree_dag import PedigreeDAG

import unittest
import matplotlib.pyplot as plt
from networkx.drawing.nx_agraph import graphviz_layout

class TestPedigreeDAG(unittest.TestCase):

    def testme(self):
        pedigree = PedigreeDAG.from_file("/gensys/mhindle/quick_sim/poc3/allPed-sex.txt")
        subgraph = pedigree.subgraph(pedigree.get_decendents(3788146172147))

        #pos=graphviz_layout(subgraph, prog='dot')
        #nx.draw(subgraph, pos, with_labels=False, arrows=False)
        #plt.savefig("/home/mhindle/path.pdf")

        sire, dam = pedigree.get_parents(3788146191054)
        print("sire %s dam %s " % (sire, dam))
        self.assertEqual(sire, 3746107111635)
        self.assertEqual(dam, 3746101701052)

        sire, dam = pedigree.get_parents(3743113491101)
        print("sire %s dam %s " % (sire, dam))
        self.assertEqual(sire, None)

        self.assertEqual(pedigree.get_relationship(3779143521305,3809179231803), "daughter")
        self.assertEqual(pedigree.get_relationship(3809179231803,3779143521305), "dam")
        self.assertEqual(pedigree.get_relationship(3779143521305,3848109802129), "granddaughter")
        self.assertEqual(pedigree.get_relationship(3779143521305,3779143501120), "husband")

        decendA = pedigree.get_decendents(3788146172147)
        left_commontree = pedigree.get_left_common_tree(3788146172147,[3785146191219])

        self.assertEqual(sum([not x in left_commontree for x in decendA]),0)
        self.assertEqual(len(left_commontree) > len(decendA),0)

        print(pedigree.get_all_siblings(3809179231803))
        print(pedigree.get_full_siblings(3890148892054))
        self.assertTrue(len(pedigree.get_full_siblings(3890148892054)) == 1)
        #self.assertEqual()
        #print(pedigree.get_decendents(3788146211035))

if __name__ == "__main__":
    #import sys;sys.argv = ['', 'Test.testName']
    unittest.main()

max_stars_repo_name: andreaskranis/genofix

max_issues_repo_name: andreaskranis/genofix

max_forks_repo_name: andreaskranis/genofix

reason: from networkx


text:
# Import packages.
import pandas as pd
import numpy as np
from .functions.timeseries_stats import rolling_mean, rolling_std
from bokeh.plotting import figure
from bokeh.models import (HoverTool, ColumnDataSource, Panel)
from bokeh.models.widgets import (Slider, CheckboxButtonGroup,
                                  Select, TextInput, Div)
from bokeh.layouts import WidgetBox, gridplot
from bokeh.palettes import Category20_16

# Defining colour palette
ts_colors = Category20_16

# Timeseries Tab.
def timeseries_tab(ts, hist_ts):
    # Set up plot
    def make_lineplot(source):
        ts_list = source.column_names
        ts_list.remove('index')
        ts_list.remove('time')

        ttp = [("Time", "$x"), ("Value", "$y")]

        plot = figure(plot_height=400, plot_width=600, tooltips=ttp,
                      title="Plot Name", x_axis_label="Axis label",
                      y_axis_label="Axis label",
                      tools="hover, pan, zoom_in, zoom_out, reset, save")
        for i, name in enumerate(ts_list):
            plot.line('time', name, source=source, line_width=3,
                      line_color=ts_colors[i], legend_label=name)

        plot.legend.location = "top_right"
        plot.legend.click_policy = "hide"
        plot.title.text_font_size = '14pt'
        plot.xaxis.axis_label_text_font_size = '12pt'
        plot.xaxis.axis_label_text_font_style = 'bold'
        plot.yaxis.axis_label_text_font_size = '12pt'
        plot.yaxis.axis_label_text_font_style = 'bold'

        return plot

    # Set up callbacks
    def update_title(attrname, old, new):
        plot.title.text = text.value

    def update_x_ax_label(attrname, old, new):
        plot.xaxis.axis_label = text2.value

    def update_y_ax_label(attrname, old, new):
        plot.yaxis.axis_label = text3.value

    def update_window(attrname, old, new):
        if window_select.value == 0:

            raw_source = ColumnDataSource(data=ts)
            source.data.update(raw_source.data)
        else:
            if rolling_method.value == "Mean":
                # Apply moving average to pandas dataframe, export to CDS type.
                ts_rmean = rolling_mean(ts, window=window_select.value)
                new_source = ColumnDataSource(data=ts_rmean)
            else:
                # Apply moving std to pandas dataframe, export to CDS type.
                ts_rstd = rolling_std(ts, window=window_select.value)
                new_source = ColumnDataSource(data=ts_rstd)

            source.data.update(new_source.data)

    # Set up widgets
    text = TextInput(title="Title of Plot", value='Plot Name')
    text.on_change('value', update_title)
    text2 = TextInput(title="x-axis Label", value='Axis Label')
    text2.on_change('value', update_x_ax_label)
    text3 = TextInput(title="y-axis Label", value='Axis Label')
    text3.on_change('value', update_y_ax_label)

    rolling_method = Select(value='Mean', title='Statistic',
                            options=['Mean', 'Std. Dev.'])
    rolling_method.on_change('value', update_window)

    window_select = Slider(start=0, end=50, step=1, value=10,
                           title='Window Size')
    window_select.on_change('value', update_window)

    div1 = Div(text="""Plot Attributes:
 Modify axis labels.
""")
    div2 = Div(text="""Rolling Window Functions:
 Specify the size and statistic for a rolling window function.
""")
    div3 = Div(text="""Data Distribution:
 Examine the distribution of time series data points.
""")

    source = ColumnDataSource(data=ts)
    plot = make_lineplot(source)

    # Set up layouts and add to document
    # Put controls in a single element.
    lineplot_controls = WidgetBox(text, div1, text2, text3, div2, window_select,
                                  rolling_method)

# Make plot with histogram
    # Find maximum range of ts.

    hist_ts_range = hist_ts.apply(lambda x: x.max() - x.min())
    range_max = hist_ts_range.max()

    # Function to make a dataset for histogram based on a list of
    # timeseries, using a default histogram bin width of 5
    def make_dataset(timeseries_list, bin_width=5):
        # Dataframe to hold information
        df_to_plot = pd.DataFrame(columns=['proportion', 'left', 'right',
                                           'f_proportion', 'f_interval',
                                           'name', 'color'])

        # Iterate through all the TS.
        for i, timeseries_name in enumerate(timeseries_list):
            # Subset to the carrier
            subset = hist_ts[timeseries_name]

            # Create a histogram with default 5 minute bins
            arr_hist, edges = np.histogram(subset,
                                           bins=int(range_max / bin_width))

            # Divide the counts by the total to get a proportion
            arr_df = pd.DataFrame({'proportion': arr_hist / np.sum(arr_hist),
                                   'left': edges[:-1], 'right': edges[1:]})

            # Format the proportion
            arr_df['f_proportion'] = ['%0.5f' % proportion for
                                      proportion in arr_df['proportion']]

            # Format the interval
            arr_df['f_interval'] = ['%d to %d units' % (left, right) for left,
                                    right in zip(arr_df['left'],
                                    arr_df['right'])]

            # Assign the carrier for labels
            arr_df['name'] = timeseries_name

            # Color each carrier differently
            arr_df['color'] = Category20_16[i]

            # Add to the overall dataframe
            df_to_plot = df_to_plot.append(arr_df)

        # Overall dataframe
        df_to_plot = df_to_plot.sort_values(['name', 'left'])

        return ColumnDataSource(df_to_plot)

    def style(hist):
        # Title
        hist.title.align = 'left'
        hist.title.text_font_size = '14pt'

        # Axis titles
        hist.xaxis.axis_label_text_font_size = '12pt'
        hist.xaxis.axis_label_text_font_style = 'bold'
        hist.yaxis.axis_label_text_font_size = '12pt'
        hist.yaxis.axis_label_text_font_style = 'bold'

        # Tick labels
        hist.xaxis.major_label_text_font_size = '12pt'
        hist.yaxis.major_label_text_font_size = '12pt'

        return hist

    def make_plot(hist_src):
        # Blank plot with correct labels
        hist = figure(plot_width=700, plot_height=700,
                      title='Histogram of Times Series Data Points',
                      x_axis_label='Bin width', y_axis_label='Proportion')

        # Quad glyphs to create a histogram
        hist.quad(source=hist_src, bottom=0, top='proportion', left='left',
                  right='right',
                  color='color',
                  fill_alpha=0.7, hover_fill_color='color', legend='name',
                  hover_fill_alpha=1.0, line_color='black')

        # Hover tool with vline mode
        hover = HoverTool(tooltips=[('Times Series',
                                     '@name'),
                                    ('x', '@f_interval'),
                                    ('Proportion',
                                    '@f_proportion')], mode='vline')

        hist.add_tools(hover)

        # Styling.
        hist = style(hist)

        return hist

    def update(attr, old, new):
        hist_ts_to_plot = [hist_ts_selection.labels[i]
                           for i in hist_ts_selection.active]
        new_hist_src = make_dataset(hist_ts_to_plot,
                                    bin_width=binwidth_select.value)
        hist_src.data.update(new_hist_src.data)

    # Widgets for interactivity.
    # Carriers and colors
    available_hist_ts = hist_ts.columns.tolist()
    available_hist_ts.sort()

    hist_ts_selection = CheckboxButtonGroup(labels=available_hist_ts,
                                            active=[0, 1])
    hist_ts_selection.on_change('active', update)

    binwidth_select = Slider(start=1, end=10,
                             step=1, value=5,
                             title='Bin Width')
    binwidth_select.on_change('value', update)

    # Call to set initial dataset.
    print(hist_ts_selection.labels, hist_ts_selection.active)
    # Initial ts and data source
    initial_hist_ts = available_hist_ts

    hist_src = make_dataset(initial_hist_ts, bin_width=binwidth_select.value)

    hist = make_plot(hist_src)

    # Put controls in a single element.
    histogram_controls = WidgetBox(div3, hist_ts_selection, binwidth_select)

    # Create a row layout
    grid = gridplot([[lineplot_controls, plot], [histogram_controls, hist]],
                    plot_width=500, plot_height=500)

    # Make a tab with the layout
    tab = Panel(child=grid, title='Time Series Plots')

    return tab

max_stars_repo_name: goodteamname/spino

max_issues_repo_name: goodteamname/spino

max_forks_repo_name: goodteamname/spino


text:
import pathlib
pathlib.Path.ls = lambda x: list(x.iterdir())

import os
from datetime import datetime

import numpy as np
from torch.nn.utils import clip_grad_norm_
from torch.optim.lr_scheduler import StepLR
from torch.utils.data import DataLoader
from torch.utils.tensorboard import SummaryWriter
from tqdm import tqdm

from evaluate import evaluate
from onsets_and_frames import *

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

leave_one_out=None
batch_size = 8
sequence_length = 327680
model_complexity = 48

if torch.cuda.is_available() and torch.cuda.get_device_properties(torch.cuda.current_device()).total_memory < 10e9:
    batch_size //= 2
    sequence_length //= 2
    print(f'Reducing batch size to {batch_size} and sequence_length to {sequence_length} to save memory')

train_groups, validation_groups = ['train'], ['validation', 'test']

if leave_one_out is not None:
    all_years = {'2004', '2006', '2008', '2009', '2011', '2013', '2014', '2015', '2017'}
    train_groups = list(all_years - {str(leave_one_out)})
    validation_groups = [str(leave_one_out)]

dataset = MAESTRO(groups=validation_groups, sequence_length=sequence_length)
dataset.purplearrow_process(groups=validation_groups, id=2)

max_stars_repo_name: purplearrow/onsets-and-frames

max_issues_repo_name: purplearrow/onsets-and-frames

max_forks_repo_name: purplearrow/onsets-and-frames


text:
import torch
import torch.nn as nn
import torch.nn.functional as F
from collections import OrderedDict
import numpy as np
import os
from os.path import expanduser
import collections
from torch.autograd import Variable
from resnext_features import resnext101_32x4d_features
from resnext_features import resnext101_64x4d_features

class Conv(nn.Module):
    def __init__(self, n_in, n_out, stride=1):
        super(Conv, self).__init__()
        self.conv = nn.Conv2d(n_in, n_out, kernel_size=3, stride=stride, padding=1, bias=False)
        self.bn = nn.BatchNorm2d(n_out)
        self.relu = nn.ReLU(inplace=True)

    def forward(self, x):
        out = self.conv(x)
        out = self.bn(out)
        out = self.relu(out)
        return out

class Bottleneck(nn.Module):
    def __init__(self, n_in, n_out, stride=1, expansion=4):
        super(Bottleneck, self).__init__()
        self.conv1 = nn.Conv2d(n_in, n_out, kernel_size=1, bias=False)
        self.bn1 = nn.BatchNorm2d(n_out)
        self.conv2 = nn.Conv2d(n_out, n_out, kernel_size=3, stride=stride, padding=1, bias=False)
        self.bn2 = nn.BatchNorm2d(n_out)
        self.conv3 = nn.Conv2d(n_out, n_out * expansion, kernel_size=1, bias=False)
        self.bn3 = nn.BatchNorm2d(n_out * expansion)

        self.downsample = None
        if stride != 1 or n_in != n_out * expansion:
            self.downsample = nn.Sequential(
                nn.Conv2d(n_in, n_out * expansion, kernel_size=1, stride=stride, bias=False),
                nn.BatchNorm2d(n_out * expansion))

        self.relu = nn.ReLU(inplace=True)

    def forward(self, x):
        residual = x

        out = self.conv1(x)
        out = self.bn1(out)
        out = self.relu(out)
        out = self.conv2(out)
        out = self.bn2(out)
        out = self.relu(out)
        out = self.conv3(out)
        out = self.bn3(out)

        if self.downsample is not None:
            residual = self.downsample(x)

        out += residual
        out = self.relu(out)

        return out

class Null(nn.Module):
    def __init__(self):
        super(Null, self).__init__()

    def forward(self, x):
        return x

class Denoise(nn.Module):
    def __init__(self, h_in, w_in, block, fwd_in, fwd_out, num_fwd, back_out, num_back):
        super(Denoise, self).__init__()

        h, w = [], []
        for i in range(len(num_fwd)):
            h.append(h_in)
            w.append(w_in)
            h_in = int(np.ceil(float(h_in) / 2))
            w_in = int(np.ceil(float(w_in) / 2))

        if block is Bottleneck:
            expansion = 4
        else:
            expansion = 1

        fwd = []
        n_in = fwd_in
        for i in range(len(num_fwd)):
            group = []
            for j in range(num_fwd[i]):
                if j == 0:
                    if i == 0:
                        stride = 1
                    else:
                        stride = 2
                    group.append(block(n_in, fwd_out[i], stride=stride))
                else:
                    group.append(block(fwd_out[i] * expansion, fwd_out[i]))
            n_in = fwd_out[i] * expansion
            fwd.append(nn.Sequential(*group))
        self.fwd = nn.ModuleList(fwd)

        upsample = []
        back = []
        n_in = (fwd_out[-2] + fwd_out[-1]) * expansion
        for i in range(len(num_back) - 1, -1, -1):
            upsample.insert(0, nn.Upsample(size=(h[i], w[i]), mode='bilinear'))
            group = []
            for j in range(num_back[i]):
                if j == 0:
                    group.append(block(n_in, back_out[i]))
                else:
                    group.append(block(back_out[i] * expansion, back_out[i]))
            if i != 0:
                n_in = (back_out[i] + fwd_out[i - 1]) * expansion
            back.insert(0, nn.Sequential(*group))
        self.upsample = nn.ModuleList(upsample)
        self.back = nn.ModuleList(back)

        self.final = nn.Conv2d(back_out[0] * expansion, fwd_in, kernel_size=1, bias=False)

    def forward(self, x):
        out = x
        outputs = []
        for i in range(len(self.fwd)):
            out = self.fwd[i](out)
            if i != len(self.fwd) - 1:
                outputs.append(out)

        for i in range(len(self.back) - 1, -1, -1):
            out = self.upsample[i](out)
            out = torch.cat((out, outputs[i]), 1)
            out = self.back[i](out)
        out = self.final(out)
        out += x
        return out

class ResNeXt101_64x4d(nn.Module):
    def __init__(self, denoise, shrink, noise, num_classes=1000):
        super(ResNeXt101_64x4d, self).__init__()
        self.features = resnext101_64x4d_features
        self.fc = nn.Linear(2048, num_classes)

        self.denoise = denoise
        self.shrink = shrink
        self.noise = noise

    def forward(self, x, defense=False):
        outputs = []
        if defense:
            x = self.shrink * self.denoise(x)
            r = torch.autograd.Variable(
                self.noise * (torch.rand(x.size()[0], x.size()[1], x.size()[2], x.size()[3]).cuda() - 0.5),
                requires_grad=False)
            x = x + r
        outputs.append(x)
        x = self.features(x)
        outputs.append(x)
        size = x.size()
        x = x.view(size[0], size[1], -1)
        x = x.mean(2)
        x = self.fc(x)
        outputs.append(x)
        return outputs

class DenoiseLoss(nn.Module):
    def __init__(self, n, hard_mining=0, norm=False):
        super(DenoiseLoss, self).__init__()
        self.n = n
        assert (hard_mining >= 0 and hard_mining <= 1)
        self.hard_mining = hard_mining
        self.norm = norm

    def forward(self, x, y):
        loss = torch.pow(torch.abs(x - y), self.n) / self.n
        if self.hard_mining > 0:
            loss = loss.view(-1)
            k = int(loss.size(0) * self.hard_mining)
            loss, idcs = torch.topk(loss, k)
            y = y.view(-1)[idcs]

        loss = loss.mean()
        if self.norm:
            norm = torch.pow(torch.abs(y), self.n)
            norm = norm.data.mean()
            loss = loss / norm
        return loss

class Loss(nn.Module):
    def __init__(self, n, hard_mining=0, norm=False):
        super(Loss, self).__init__()
        self.loss = DenoiseLoss(n, hard_mining, norm)

    def forward(self, x, y):
        z = []
        for i in range(len(x)):
            z.append(self.loss(x[i], y[i]))
        return z

class Net(nn.Module):
    def __init__(self, input_size, block, fwd_out, num_fwd, back_out, num_back, n, shrink, noise, hard_mining=0,
                 loss_norm=False):
        super(Net, self).__init__()
        denoise = Denoise(input_size[0], input_size[1], block, 3, fwd_out, num_fwd, back_out, num_back)
        self.net = ResNeXt101_64x4d(denoise, shrink, noise)
        self.loss = Loss(n, hard_mining, loss_norm)

    def forward(self, orig_x, adv_x, requires_control=True, train=True):
        ratio = 1
        orig_outputs = self.net(ratio * orig_x)

        if requires_control:
            control_outputs = self.net(ratio * adv_x)
            control_loss = self.loss(control_outputs, orig_outputs)

        if train:
            adv_x.volatile = False
            for i in range(len(orig_outputs)):
                orig_outputs[i].volatile = False
        adv_outputs = self.net(ratio * adv_x, defense=True)
        loss = self.loss(adv_outputs, orig_outputs)

        if not requires_control:
            return orig_outputs[-1], adv_outputs[-1], loss
        else:
            return orig_outputs[-1], adv_outputs[-1], loss, control_outputs[-1], control_loss

def get_net(input_size, block, fwd_out, num_fwd, back_out, num_back, n, shrink=1, noise=0, hard_mining=0,
            loss_norm=False):
    net = Net(input_size, block, fwd_out, num_fwd, back_out, num_back, n, shrink, noise, hard_mining, loss_norm)

    return net


text:
import numpy as np
import csv
from matplotlib import pyplot as plt

def evolvingProfs(state):
	# Compiling csv input into individual variables
	dataArr = []
	with open('output_runModel/weekly_results.csv', 'r') as weeklyFile:
		csvRead = csv.reader(weeklyFile)
		for row in csvRead:
			dataArr.append(row)
	dataArr =  [ele for ele in dataArr if ele != []]
	olrs, bdry_tempDiff, surfT, netEn = dataArr[0], dataArr[1], dataArr[2], dataArr[3]

	# OLR
	fig = plt.figure()
	olr = [float(i) for i in olrs]
	plt.plot(olr)
	plt.xlabel("Weeks")
	plt.ylabel("Outgoing Longwave Radiation (W m^-2)")
	plt.savefig("graphs/weeklyOLR.png")

	# Plotting bdry_tempDiff
	fig = plt.figure()
	bdry_tempDiff = [float(i) for i in bdry_tempDiff]
	plt.plot(bdry_tempDiff)
	plt.xlabel("Weeks")
	plt.ylabel("Surface-Air Temperature Gradient (K)")
	plt.savefig("graphs/weeklybdry_tempDiff.png")

	# Plotting surfT
	fig = plt.figure()
	surfT = [float(i) for i in surfT]
	plt.plot(surfT)
	plt.xlabel("Weeks")
	plt.ylabel("Surface Temperature (K)")
	plt.savefig("graphs/weekly_surfT.png")

	#Plotting netEn
	fig = plt.figure()
	netEn = [float(i) for i in netEn]
	plt.plot(netEn)
	plt.xlabel("Weeks")
	plt.ylabel("Net Energy Level in Atmospheric Column (W m^-2)")
	plt.savefig("graphs/weekly_netEn.png")

	return 0.

max_stars_repo_name: Mihir-DG/Modelling-Planetary-Climate

max_issues_repo_name: Mihir-DG/Modelling-Planetary-Climate

max_forks_repo_name: Mihir-DG/Modelling-Planetary-Climate


text:
# This script creates an instrument for the MLBozone project

# Importing required modules

import pandas as pd
import numpy as np
import os

# Define username + filepath

username = ''
filepath = 'C:/Users/' + username + '/Documents/Data/MLBozone/raw_data2/'

# Reading in the data

data = pd.read_csv(filepath[:-10] + 'maybe_final_data_set.csv')
files = os.listdir(filepath)
df = pd.DataFrame()

for f in files:

    tmp = pd.read_csv(filepath + f, header = None)
    df = pd.concat([df, tmp], axis = 0)

# Sort by date

df = df.sort_values([0])

# Remove games before 2009

df = df[df[0] > 20090000].reset_index(drop = True)

# Convert FLO to MIA so same team is identified

df[[3,6]] = df[[3,6]].replace(['FLO'], ['MIA'])

# Create a list of teams

teams = list(df[3].unique())

# Main loop creating the instrument

iv10 = []
iv20 = []
iv50 = []
iv100 = []
dates10 = []
dates20 = []
dates50 = []
dates100 = []
ts10 = []
ts20 = []
ts50 = []
ts100 = []

for t in teams:

    print(t)

    tmpa = df[df[3] == t]
    tmph = df[df[6] == t]
    tmp = pd.concat([tmpa, tmph], axis = 0)
    tmp = tmp.sort_values([0]).reset_index(drop = True)
    a = [int(tmp[3][x] == t) for x in range(len(tmp))]
    h = [int(tmp[6][x] == t) for x in range(len(tmp))]
    aw = [1 if tmp[9][i] > tmp[10][i] else 0 for i in range(len(tmp))]
    hw = [1 if tmp[10][i] > tmp[9][i] else 0 for i in range(len(tmp))]
    w = [a[i]*aw[i] + h[i]*hw[i] for i in range(len(a))]

    iv10x = [np.mean(w[x-10:x]) for x in range(10,len(w))]
    iv20x = [np.mean(w[x-20:x]) for x in range(20,len(w))]
    iv50x = [np.mean(w[x-50:x]) for x in range(50,len(w))]
    iv100x = [np.mean(w[x-100:x]) for x in range(100,len(w))]

    dates10x = [tmp[0][x] for x in range(10,len(tmp))]
    dates20x = [tmp[0][x] for x in range(20,len(tmp))]
    dates50x = [tmp[0][x] for x in range(50,len(tmp))]
    dates100x = [tmp[0][x] for x in range(100,len(tmp))]

    ts10x = [t]*len(dates10x)
    ts20x = [t]*len(dates20x)
    ts50x = [t]*len(dates50x)
    ts100x = [t]*len(dates100x)

    iv10 = iv10 + iv10x
    iv20 = iv20 + iv20x
    iv50 = iv50 + iv50x
    iv100 = iv100 + iv100x
    dates10 = dates10 + dates10x
    dates20 = dates20 + dates20x
    dates50 = dates50 + dates50x
    dates100 = dates100 + dates100x
    ts10 = ts10 + ts10x
    ts20 = ts20 + ts20x
    ts50 = ts50 + ts50x
    ts100 = ts100 + ts100x

# Create some dataframes

iv10 = pd.Series(iv10, name = 'IV10')
iv20 = pd.Series(iv20, name = 'IV20')
iv50 = pd.Series(iv50, name = 'IV50')
iv100 = pd.Series(iv100, name = 'IV100')
dates10 = pd.Series(dates10, name = 'Date')
dates20 = pd.Series(dates20, name = 'Date')
dates50 = pd.Series(dates50, name = 'Date')
dates100 = pd.Series(dates100, name = 'Date')
ts10 = pd.Series(ts10, name = 'Team')
ts20 = pd.Series(ts20, name = 'Team')
ts50 = pd.Series(ts50, name = 'Team')
ts100 = pd.Series(ts100, name = 'Team')

newdf10 = pd.concat([iv10,dates10,ts10], axis = 1)
newdf20 = pd.concat([iv20,dates20,ts20], axis = 1)
newdf50 = pd.concat([iv50,dates50,ts50], axis = 1)
newdf100 = pd.concat([iv100,dates100,ts100], axis = 1)

x10 = [str(dates10[x]) + str(ts10[x]) for x in range(len(dates10))]
x20 = [str(dates20[x]) + str(ts20[x]) for x in range(len(dates20))]
x50 = [str(dates50[x]) + str(ts50[x]) for x in range(len(dates50))]
x100 = [str(dates100[x]) + str(ts100[x]) for x in range(len(dates100))]

x10 = pd.Series(x10, name = 'x')
x20 = pd.Series(x20, name = 'x')
x50 = pd.Series(x50, name = 'x')
x100 = pd.Series(x100, name = 'x')

newdf10 = pd.concat([newdf10, x10], axis = 1).reset_index(drop = True)
newdf20 = pd.concat([newdf20, x20], axis = 1).reset_index(drop = True)
newdf50 = pd.concat([newdf50, x50], axis = 1).reset_index(drop = True)
newdf100 = pd.concat([newdf100, x100], axis = 1).reset_index(drop = True)

# Create final instruments data

i10 = []
i20 = []
i50 = []
i100 = []
obs = [str(data.Date[i]) + str(data.City[i]) for i in range(len(data))]
prev_t = ''

for o in obs:

    print(o)
    t = o[-3:]

    if t == prev_t:

        tmp10 = newdf10[newdf10.x == o].reset_index(drop = True)
        tmp20 = newdf20[newdf20.x == o].reset_index(drop = True)
        tmp50 = newdf50[newdf50.x == o].reset_index(drop = True)
        tmp100 = newdf100[newdf100.x == o].reset_index(drop = True)

        try:

            i10.append(tmp10.IV10[0])
            i20.append(tmp20.IV20[0])
            i50.append(tmp50.IV50[0])
            i100.append(tmp100.IV100[0])

        except:

            i10.append(None)
            i20.append(None)
            i50.append(None)
            i100.append(None)

    else:

        tmp10 = newdf10[newdf10.Team == t]
        tmp10 = tmp10[tmp10.Date < 20100000].reset_index(drop = True)
        i10.append(tmp10.IV10[len(tmp10)-1])

        tmp20 = newdf20[newdf20.Team == t]
        tmp20 = tmp20[tmp20.Date < 20100000].reset_index(drop = True)
        i20.append(tmp20.IV20[len(tmp20)-1])

        tmp50 = newdf50[newdf50.Team == t]
        tmp50 = tmp50[tmp50.Date < 20100000].reset_index(drop = True)
        i50.append(tmp50.IV50[len(tmp50)-1])

        tmp100 = newdf100[newdf100.Team == t]
        tmp100 = tmp100[tmp100.Date < 20100000].reset_index(drop = True)
        i100.append(tmp100.IV100[len(tmp100)-1])

# Removing Nones

for i in range(len(i10)):

    if i10[i] == None:

        i10[i] = i10[i-1]

for i in range(len(i20)):

    if i20[i] == None:

        i20[i] = i20[i-1]

for i in range(len(i50)):

    if i50[i] == None:

        i50[i] = i50[i-1]

for i in range(len(i100)):

    if i100[i] == None:

        i100[i] = i100[i-1]

# Add instruments to dataframe

i10 = pd.Series(i10, name = 'IV10')
i20 = pd.Series(i20, name = 'IV20')
i50 = pd.Series(i50, name = 'IV50')
i100 = pd.Series(i100, name = 'IV100')
data = pd.concat([data, i10, i20, i50, i100], axis = 1)

# Write data to file

data.to_csv(filepath[:-10] + 'mlbozone_data.csv', index = False)

# This works as an instrument because better team performance can
# bring people to the stadium who otherwise wouldn't have been out

max_stars_repo_name: cat-astrophic/MLBozone

max_issues_repo_name: cat-astrophic/MLBozone

max_forks_repo_name: cat-astrophic/MLBozone


text:
import numpy as np
import scipy.sparse as sp

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

class SNOPT_work(object):
    def __init__(self,lencw=500,leniw=500,lenrw=500):
        self.setup(lencw,leniw,lenrw)

    def setup(self,lencw=500,leniw=500,lenrw=500):
        self.lencw   = lencw
        self.leniw   = leniw
        self.lenrw   = lenrw
        self.cw      = np.empty((lencw,8),dtype='|S1')
        self.iw      = np.zeros(leniw,'i',order='F')
        self.rw      = np.zeros(lenrw,float,order='F')

    def work_resize(self,mincw,miniw,minrw):
        if mincw > self.lencw:
            tmp = np.copy(self.cw)
            self.cw = np.empty((mincw,8),dtype='|S1')
            self.cw[0:self.lencw] = tmp[0:self.lencw]
            del tmp

            self.lencw = mincw
            self.iw[6] = mincw

        if miniw > self.leniw:
            tmp = self.iw
            self.iw = np.zeros(miniw,'i',order='F')
            self.iw[0:self.leniw] = tmp[0:self.leniw]
            del tmp

            self.leniw = miniw
            self.iw[4] = miniw

        if minrw > self.lenrw:
            tmp = self.rw
            self.rw = np.zeros(minrw,float,order='F')
            self.rw[0:self.lenrw] = tmp[0:self.lenrw]
            del tmp

            self.lenrw = minrw
            self.iw[2] = minrw

max_stars_repo_name: jainajinkya/snopt-python

max_issues_repo_name: jainajinkya/snopt-python

max_forks_repo_name: jainajinkya/snopt-python

reason: import numpy,import scipy


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,py:light
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 1.0.0
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %matplotlib inline
# %autosave 0
# %load_ext autoreload
# %autoreload 2

import ast
import logging
import numpy as np
import pandas as pd
import pkg_resources
import seaborn as sns
import matplotlib.pyplot as plt
from natsort import natsorted
from crispy.QCPlot import QCplot
from crispy.DataImporter import Sample
from crispy.CrispyPlot import CrispyPlot
from scipy.stats import spearmanr, pearsonr, ttest_ind
from crispy.CRISPRData import CRISPRDataSet, Library, ReadCounts
from minlib.Utils import project_score_sample_map, density_interpolate
from sklearn.metrics import (
    recall_score,
    precision_score,
    average_precision_score,
    precision_recall_curve,
)

LOG = logging.getLogger("Crispy")
DPATH = pkg_resources.resource_filename("crispy", "data/")
RPATH = pkg_resources.resource_filename("notebooks", "minlib/reports/")

# Project Score samples acquired with Kosuke_Yusa v1.1 library
#

avana_cmap = pd.read_csv(f"{DPATH}/crispr_manifests/Avana_DepMap20Q1_sample_map.csv.gz", index_col=0)

avana = CRISPRDataSet("Avana_DepMap19Q2")

avana_counts = avana.counts.remove_low_counts(avana.plasmids)

avana_fc = avana_counts.norm_rpm().foldchange(avana.plasmids)

# KY v1.1 library
#

av_lib = Library.load_library("MasterLib_v1.csv.gz").query("Library == 'Avana'")
av_lib = av_lib[av_lib.index.isin(avana_counts.index)]

# Minimal library (top 2)
#

NGUIDES, REMOVE_DISCORDANT = 2, False
ml_lib_name = (
    f"MinimalLib_top{NGUIDES}{'_disconcordant' if REMOVE_DISCORDANT else ''}_Avana.csv.gz"
)
ml_lib = Library.load_library(ml_lib_name).query("Library == 'Avana'")
ml_lib = ml_lib[ml_lib.index.isin(avana_counts.index)]
ml_lib = ml_lib.loc[[i for i in ml_lib.index if not i.startswith("CTRL0")]]

# Genes overlap
#

genes = set(av_lib["Approved_Symbol"]).intersection(ml_lib["Approved_Symbol"])
LOG.info(f"Genes={len(genes)}")

# Comparisons
#

FDR_THRES = 0.01

libraries = dict(
    All=dict(name="All", lib=av_lib), Minimal=dict(name="Minimal", lib=ml_lib)
)

# Minimal counts
#

for ltype in libraries:
    lib = libraries[ltype]["lib"].copy()

    # sgRNA fold-changes
    fc = avana_counts.copy().loc[lib.index].norm_rpm().foldchange(avana.plasmids)

    # Gene fold-changes
    fc_gene = fc.groupby(lib["Approved_Symbol"]).mean()

    # Gene average fold-changes
    fc_gene_avg = fc_gene.groupby(avana_cmap["DepMap_ID"], axis=1).mean()

    # Gene average scaled fold-changes
    fc_gene_avg_scl = ReadCounts(fc_gene_avg).scale()

    # FDR threshold
    fpr = pd.DataFrame(
        {
            s: QCplot.aroc_threshold(fc_gene_avg[s], fpr_thres=FDR_THRES)
            for s in fc_gene_avg
        },
        index=["auc", "thres"],
    ).T
    fpr["auc"] = [
        QCplot.aroc_threshold(fc_gene_avg[s], fpr_thres=0.2)[0] for s in fpr.index
    ]

    # Gene binarised
    fc_gene_avg_bin = (fc_gene_avg[fpr.index] < fpr["thres"]).astype(int)

    # Store
    libraries[ltype]["fc"] = fc
    libraries[ltype]["fc_gene"] = fc_gene
    libraries[ltype]["fc_gene_avg"] = fc_gene_avg
    libraries[ltype]["fc_gene_avg_scl"] = fc_gene_avg_scl
    libraries[ltype]["fc_gene_avg_bin"] = fc_gene_avg_bin
    libraries[ltype]["fpr"] = fpr

# Recovered dependencies per cell line
#

deprecovered = pd.DataFrame(
    [
        dict(
            samples=s,
            recall=recall_score(
                libraries["All"]["fc_gene_avg_bin"].loc[genes, s],
                libraries["Minimal"]["fc_gene_avg_bin"].loc[genes, s],
            ),
            ap=average_precision_score(
                libraries["All"]["fc_gene_avg_bin"].loc[genes, s],
                -libraries["Minimal"]["fc_gene_avg"].loc[genes, s],
            ),
            precision=precision_score(
                libraries["All"]["fc_gene_avg_bin"].loc[genes, s],
                libraries["Minimal"]["fc_gene_avg_bin"].loc[genes, s],
            ),
            spearman=spearmanr(
                libraries["All"]["fc_gene_avg"].loc[genes, s],
                libraries["Minimal"]["fc_gene_avg"].loc[genes, s],
            )[0],
        )
        for s in libraries["All"]["fc_gene_avg_bin"]
    ]
)
deprecovered = deprecovered.dropna(subset=["ap"]).sort_values("ap", ascending=False)

# Disconcordance between number of dependent cell lines found per gene
#

gene_deps = pd.concat(
    [libraries[l]["fc_gene_avg_bin"].sum(1).rename(l) for l in libraries],
    axis=1,
    sort=False,
).dropna()
gene_deps["difference"] = gene_deps.eval("All - Minimal")
gene_deps = gene_deps.loc[gene_deps.sort_values("difference").index]
gene_deps["match"] = (gene_deps["difference"].abs() > 100).replace(
    {True: "Disconcordant", False: "Concordant"}
)
gene_deps["n_guides"] = (
    av_lib["Approved_Symbol"].dropna().value_counts().loc[gene_deps.index].values
)
gene_deps.to_excel(
    f"{RPATH}/{ml_lib_name.split('.')[0]}_benchmark_gene_concordance_avana.xlsx"
)

# Gene fold-change correlation
#

plt.figure(figsize=(2.5, 1.5), dpi=600)
sns.distplot(
    deprecovered["spearman"], hist=False, kde_kws={"cut": 0, "shade": True}, color="k"
)
plt.xlabel("Spearman's R")
plt.title(
    f"Samples gene fold-change correlation\n(mean={deprecovered['spearman'].mean():.2f})"
)
plt.grid(True, ls=":", lw=0.1, alpha=1.0, zorder=0, axis="x")
plt.savefig(
    f"{RPATH}/{ml_lib_name.split('.')[0]}_benchmark_dependencies_fc_corr_histogram_avana.pdf",
    bbox_inches="tight",
    transparent=True,
)
plt.close("all")

# Effect sizes of recap dependencies
#

plot_df = pd.concat(
    [
        libraries["All"]["fc_gene_avg_bin"].loc[genes].unstack().rename("All_bin"),
        libraries["Minimal"]["fc_gene_avg_bin"]
        .loc[genes]
        .unstack()
        .rename("Minimal_bin"),
        libraries["All"]["fc_gene_avg_scl"].loc[genes].unstack().rename("All_fc"),
    ],
    axis=1,
    sort=False,
).reset_index()
plot_df["Recap"] = [
    a if a != m else sum([a, m]) for a, m in plot_df[["All_bin", "Minimal_bin"]].values
]
plot_df["Recap"] = plot_df["Recap"].replace({2: "Yes", 1: "No", 0: "NA"})

ttest_ind(
    plot_df.query("Recap == 'Yes'")["All_fc"],
    plot_df.query("Recap == 'No'")["All_fc"],
    equal_var=False,
)

order = ["Yes", "No"]
pal = {"Yes": CrispyPlot.PAL_DBGD[1], "No": CrispyPlot.PAL_DBGD[0]}
samples = set(plot_df["model_id"])

sns.boxplot(
    "Recap",
    "All_fc",
    data=plot_df,
    order=order,
    palette=pal,
    saturation=1,
    showcaps=False,
    boxprops=dict(linewidth=0.3),
    whiskerprops=dict(linewidth=0.3),
    flierprops=dict(
        marker="o",
        markerfacecolor="black",
        markersize=1.0,
        linestyle="none",
        markeredgecolor="none",
        alpha=0.6,
    ),
    notch=True,
)
plt.grid(True, ls=":", lw=0.1, alpha=1.0, zorder=0, axis="y")
plt.xlabel("Dependency recapitulated")
plt.ylabel("All sgRNAs fold-change (scaled)")
plt.gcf().set_size_inches(0.75, 2.0)
plt.savefig(
    f"{RPATH}/{ml_lib_name.split('.')[0]}_benchmark_dependencies_effect_sizes_boxplots_avana.pdf",
    bbox_inches="tight",
    transparent=True,
)
plt.close("all")

# AROC and FC threshold scatters
#

for metric in ["auc", "thres"]:
    plot_df = pd.concat(
        [libraries[l]["fpr"][metric].rename(l) for l in libraries], axis=1
    )
    plot_df["density"] = density_interpolate(plot_df["All"], plot_df["Minimal"])
    plot_df = plot_df.sort_values("density")

    x_min, x_max = plot_df.iloc[:, :2].min().min(), plot_df.iloc[:, :2].max().max()

    plt.scatter(
        plot_df["All"],
        plot_df["Minimal"],
        c=plot_df["density"],
        marker="o",
        edgecolor="",
        cmap="Spectral_r",
        s=3,
        alpha=0.7,
    )

    lims = [x_max, x_min]
    plt.plot(lims, lims, "k-", lw=0.3, zorder=0)

    plt.grid(True, ls=":", lw=0.1, alpha=1.0, zorder=0)

    plt.xlabel("All sgRNAs")
    plt.ylabel("Minimal Library - 2 sgRNAs")
    plt.title(f"AROC essential genes" if metric == "auc" else f"Fold-change threshold")

    plt.gcf().set_size_inches(2, 2)
    plt.savefig(
        f"{RPATH}/{ml_lib_name.split('.')[0]}_benchmark_scatter_{metric}_avana.pdf",
        bbox_inches="tight",
        transparent=True,
    )
    plt.close("all")

# Number of dependencies scatter
#

gene_deps["density"] = density_interpolate(gene_deps["All"], gene_deps["Minimal"])
gene_deps["n_guides"] = (
    av_lib["Approved_Symbol"].dropna().value_counts().loc[gene_deps.index].values
)
gene_deps = gene_deps.sort_values("density")

x_min, x_max = gene_deps.iloc[:, :2].min().min(), gene_deps.iloc[:, :2].max().max()

ax = plt.gca()

for t, df in gene_deps.groupby("match"):
    ax.scatter(
        df["All"],
        df["Minimal"],
        c=df["density"],
        marker="o",
        edgecolor="#e34a33" if t == "Disconcordant" else "",
        cmap="Spectral_r",
        s=3,
        alpha=0.7,
    )

ax.grid(True, ls=":", lw=0.1, alpha=1.0, zorder=0)

lims = [x_max, x_min]
ax.plot(lims, lims, "k-", lw=0.3, zorder=0)

cor, pval = pearsonr(plot_df["All"], plot_df["Minimal"])
annot_text = f"Spearman's R={cor:.2g}"

ax.text(0.05, 0.95, annot_text, fontsize=4, transform=ax.transAxes, ha="left")

ax.set_xlabel("All sgRNAs")
ax.set_ylabel("Minimal Library - 2 sgRNAs")
ax.set_title("Number of dependencies per gene")

plt.gcf().set_size_inches(2, 2)
plt.savefig(
    f"{RPATH}/{ml_lib_name.split('.')[0]}_benchmark_dependencies_recap_sctter_avana.pdf",
    bbox_inches="tight",
    transparent=True,
)
plt.close("all")

# Recover plot
#

plot_df = deprecovered.reset_index(drop=True)

_, ax1 = plt.subplots(figsize=(3, 2))

ax1.scatter(plot_df.index, plot_df["ap"], c=QCplot.PAL_DBGD[0], s=5, label=t)

for t in [0.8, 0.9]:
    s_t = plot_df[plot_df["ap"] < t].index[0]
    ax1.axvline(s_t, lw=0.1, c="k", zorder=0)

    s_t_text = f"{((s_t / plot_df.shape[0]) * 100):.1f}%"
    ax1.text(s_t, plt.ylim()[0], s_t_text, fontsize=6, c="k")

plt.grid(True, ls=":", lw=0.1, alpha=1.0, zorder=0, axis="y")

plt.xlabel("Number of cell lines")
plt.ylabel(f"Average precission score")

plt.savefig(
    f"{RPATH}/{ml_lib_name.split('.')[0]}_benchmark_recall_plot_avana.pdf",
    bbox_inches="tight",
    transparent=True,
)
plt.close("all")

max_stars_repo_name: QuantitativeBiology/crispy

reason: import numpy,from scipy


text:
# Copyright (c) OpenMMLab. All rights reserved.
import importlib
import os.path as osp
import tempfile

import mmcv
import pytest

from mmdeploy.apis import torch2torchscript
from mmdeploy.utils import IR, Backend
from mmdeploy.utils.test import get_random_name

ts_file = tempfile.NamedTemporaryFile(suffix='.pt').name
input_name = get_random_name()
output_name = get_random_name()

def get_deploy_cfg(input_name, output_name):
    return mmcv.Config(
        dict(
            ir_config=dict(
                type=IR.TORCHSCRIPT.value,
                input_names=[input_name],
                output_names=[output_name],
                input_shape=None),
            codebase_config=dict(type='mmedit', task='SuperResolution'),
            backend_config=dict(type=Backend.TORCHSCRIPT.value)))

def get_model_cfg():
    return mmcv.Config(
        dict(
            model=dict(
                pretrained=None,
                type='BasicRestorer',
                generator=dict(
                    type='RRDBNet',
                    in_channels=3,
                    out_channels=3,
                    mid_channels=64,
                    num_blocks=23,
                    growth_channels=32),
                pixel_loss=dict(
                    type='L1Loss', loss_weight=1.0, reduction='mean')),
            test_cfg=dict(metrics='PSNR'),
            test_pipeline=[
                dict(
                    type='LoadImageFromFile',
                    io_backend='disk',
                    key='lq',
                    flag='unchanged'),
                dict(
                    type='LoadImageFromFile',
                    io_backend='disk',
                    key='gt',
                    flag='unchanged'),
                dict(type='RescaleToZeroOne', keys=['lq', 'gt']),
                dict(
                    type='Normalize',
                    keys=['lq', 'gt'],
                    mean=[0, 0, 0],
                    std=[1, 1, 1],
                    to_rgb=True),
                dict(
                    type='Collect',
                    keys=['lq', 'gt'],
                    meta_keys=['lq_path', 'lq_path']),
                dict(type='ImageToTensor', keys=['lq', 'gt'])
            ]))

@pytest.mark.parametrize('input_name', [input_name])
@pytest.mark.parametrize('output_name', [output_name])
@pytest.mark.skipif(
    not importlib.util.find_spec('mmedit'), reason='requires mmedit')
def test_torch2torchscript(input_name, output_name):
    import numpy as np
    deploy_cfg = get_deploy_cfg(input_name, output_name)
    torch2torchscript(
        np.random.rand(8, 8, 3),
        '',
        ts_file,
        deploy_cfg,
        model_cfg=get_model_cfg(),
        device='cpu')

    print(ts_file)
    assert osp.exists(ts_file)

max_stars_repo_name: grimoire/mmdeploy

max_issues_repo_name: grimoire/mmdeploy

max_forks_repo_name: grimoire/mmdeploy


text:
import cv2
import scipy.io as sio
import numpy as np
from dehaze_patchMap_dehaze import dehaze_patchMap

hazy_img_dir = 'image/'
patchMap_dir = 'patchMap/'
save_dir = 'result/'
imgname= 'girls.tif'

image = cv2.imread(hazy_img_dir + imgname)
if image.shape[0] != 480 or image.shape[1] != 640:
    print('resize image tp 640*480')
    image = cv2.resize(image,(640,480))

print(imgname)
patchMapname = imgname[:-4] + '.mat'

patchMap = sio.loadmat(patchMap_dir + patchMapname)
patchMap = np.array(patchMap['patchMap'])

recover_result, tx = dehaze_patchMap(image, 0.95, patchMap)

savename_result = save_dir + 'py_recover_' + imgname

cv2.imwrite(savename_result, recover_result)


text:
import numpy as np
import parameters as param

from art.estimators.classification import TensorFlowV2Classifier
from art.defences.trainer import AdversarialTrainer
from art.attacks.evasion import FastGradientMethod
from art.attacks.evasion import CarliniLInfMethod
from art.attacks.evasion import ProjectedGradientDescent
from art.attacks.evasion import BasicIterativeMethod
from art.attacks.evasion import SaliencyMapMethod
from art.attacks.evasion import AutoProjectedGradientDescent
from art.attacks.evasion import DeepFool, NewtonFool
from art.attacks.evasion import SquareAttack, SpatialTransformation
from art.attacks.evasion import ShadowAttack, Wasserstein

import tensorflow as tf
tf.compat.v1.logging.set_verbosity(tf.compat.v1.logging.ERROR)
from tensorflow.python.client import device_lib

DATA_DIR = param.DATA_DIR
DATASET_NAMES = param.DATASET_NAMES

def load_data(dataset_name):
    assert dataset_name in DATASET_NAMES
    x_train = np.load(DATA_DIR + dataset_name + '/benign/x_train.npy')
    y_train = np.load(DATA_DIR + dataset_name + '/benign/y_train.npy')
    x_test = np.load(DATA_DIR + dataset_name + '/benign/x_test.npy')
    y_test = np.load(DATA_DIR + dataset_name + '/benign/y_test.npy')
    return x_train, y_train, x_test, y_test

def call_function_by_attack_name(attack_name):
    if attack_name not in param.ATTACK_NAMES:
        print('Unsupported attack: {}'.format(attack_name))
        sys.exit(1)
    return {
        param.APGD: AutoProjectedGradientDescent,
        param.BIM: BasicIterativeMethod,
        param.CW: CarliniLInfMethod,
        param.DF: DeepFool,
        param.FGSM: FastGradientMethod,
        param.JSMA: SaliencyMapMethod,
        param.NF: NewtonFool,
        param.PGD: ProjectedGradientDescent,
        param.SA: SquareAttack,
        param.SHA: ShadowAttack,
        param.ST: SpatialTransformation,
        param.WA: Wasserstein
    }[attack_name]

def get_available_gpus():
    local_device_protos = device_lib.list_local_devices()
    return [x.name for x in local_device_protos if x.device_type == 'GPU']


text:
"""DEFINES THE INVERSEDYNAMICS SOLVER, A Solver for solving the joint based model of a dog."""

from scipy import optimize, signal

from data.data_loader import C3DData, load_force_plate_data, ForcePlateData, SMALData, get_delay_between, DataSources, \
	path_join
from vis.utils import *
from vis import visualisations
from dynamics.footfall_detector import FootfallDetector
from tqdm import tqdm

# pure constants (no optimisation needed)
g = 9.81
freq_forceplate = 100  # Hz
foot_joint_labels = ["front left", "front right", "rear left", "rear right"]
foot_joint_indices = [0, 9, 23, 20]  # for set 2 3r3

class Model:
	"""ID Model, with all parameters derived/optimised"""

	def __init__(self):
		# CONSTANTS
		self.paws = {}
		self.bone_density = 1950  # Estimate - needs refining! From paper: Development of a neuromusculoskeletal computer model in a chondrodystrophic dog.
		self.muscle_density = 1060  # From above

		# params to optimise
		self.bone_length_definitions = {
			"normal": lambda l: dict(inner_radius=0.01, outer_radius=0.05, displacement=0),
			"body": lambda l: dict(inner_radius=l / 20, outer_radius=l / 7, displacement=l / 4 - l / 20), }

		# Paw parameters. All scaled to be in standard form - exponent in separate dict.
		self.paw_params_normalised = {
			"L0_front": 6.9,  # 6.9 # in .1mm
			"L0_rear": 6.9,  # in .1mm
			"k_front": 3.42 * .18,  # in kN/m
			"k_rear": 2.0 * .21,  # in kN/m
			"c_front": 20,
			"c_rear": 20,
			"k_rear_prop": 0.85,  # k = k_rear * m **.85
			"frame_delay": 0  # Used for analysis of paw treadmill forces. Not used for normal ID solver
		}
		self.paw_exponents = {
			"L0_front": -4,
			"L0_rear": -4,
			"k_front": 3,
			"k_rear": 3,
			"c_front": 0,
			"c_rear": 0,
			"k_rear_prop": 0,
			"frame_delay": 0
		}

		self.calc_paw_params()

		self.freq_par_data = 200

		# weightings used in dynamics calculations
		self.equation_weighting = {
			"Inertial": 2,
			"Rotational": 2,
			"Leg spring": 0.5,
			"Paw spring": 1,
		}

	def calc_paw_params(self):
		"""Calculates paw parameters (separate function for optimisation purposes)"""
		for param, val in self.paw_params_normalised.items():
			self.paws[param] = val * 10 ** (self.paw_exponents[param])

	def edit_paw_param(self, param, val):
		"""Edit paw parameter (separate for optimisation purposes)"""
		self.paw_params_normalised[param] = val
		self.calc_paw_params()

model = Model()

def time_deriv(X, dt):
	"""Finds the time derivative of a given series of data.
	Always treats the first dimension as time - works for any number of dimensions (n_frames, M, N, O, ...).
	For all except first and last val, calcs difference over 2 timesteps"""

	diff = np.zeros_like(X)

	diff[0] = X[1] - X[0]
	diff[1:-1] = (X[2:] - X[:-2]) / 2
	diff[-1] = X[-1] - X[-2]

	return diff * 1 / dt

def nth_time_deriv(X, dt, n=2):
	"""Recursively get the nth time derivative"""

	if n == 1:
		return time_deriv(X, dt)
	else:
		return time_deriv(nth_time_deriv(X, dt, n=n - 1), dt)

def get_principal_axes(vector=Vector(1, 0, 0), cardinal=np.identity(3)):
	"""Given a vector, devise a basis of principle axis with any two perpendicular vectors (for application of an
	axisymmetric object - cylinder) """

	i, j, k = cardinal
	K = vector.unit()
	# Now find any two perp vectors to K
	if not K.is_parallel(i):
		I = K.cross(i).unit()
		J = K.cross(I).unit()
	else:
		I = K.cross(j).unit()
		J = K.cross(I).unit()

	return np.array([I, J, K])

def I_cylinder(density, length, radius):
	mass = density * np.pi * (radius ** 2) * length
	Ixx, Izz = (length ** 2) / 12 + (radius ** 2) / 4, radius ** 2 / 2
	return mass * np.diag([Ixx, Ixx, Izz])

class DoubleCylinder:
	"""An object comprised of a cylinder of given length between two end points, of radius inner_radius and density bone_density,
	and an outer cylinder that does NOT share the same central axis, of radius outer_radius, displaced by a distance <displacement> normally from the centerline.

	Cylinder is defined with the centerline vertical (z direction), and the displacement always in the normal closest to the z direction downwards.

	For InverseDynamics calculations, this object will have a start and end index, which correspond to the joint indices in which the end point data is held.
	"""

	def __init__(self, start, end, length, inner_radius, outer_radius, displacement, freq=50.0, name=""):

		self.name = name
		self.freq = freq  # Frequency, in Hz

		self.start = start
		self.end = end

		self.length = length
		self.displacement = displacement

		if outer_radius is None: outer_radius = inner_radius

		self.inner_mass = model.bone_density * np.pi * inner_radius ** 2 * self.length
		self.outer_mass = model.muscle_density * np.pi * self.length * (outer_radius ** 2 - inner_radius ** 2)
		self.mass = self.inner_mass + self.outer_mass

		I_bone = I_cylinder(model.bone_density, length, inner_radius)
		I_muscle = I_cylinder(model.muscle_density, length, outer_radius) - I_cylinder(model.muscle_density, length,
																					   inner_radius)

		# By parallel axis theorem, add component of I due to outer radius being displaced from the centerline axis
		I_axis_displacement = np.zeros((3, 3))
		I_axis_displacement[0, 0] = self.outer_mass * displacement ** 2

		self.I = I_bone + I_muscle + I_axis_displacement  # Inertia tensor in a reference frame in which the bone is lengthwise facing upwards

	def get_kinematics(self, data):
		"""Given a numpy array of time, data, of shape (n_frames, 2, 3),
		giving the position data of both ends of the cylinder over time, compute the kinematics of the cylinder"""

		X = self.X = np.array(data)  # positions
		V = self.V = time_deriv(X, 1 / self.freq)  # velocities
		A = self.A = time_deriv(V, 1 / self.freq)  # accelerations

		self.XG = np.mean(X, axis=1)  # average over X
		self.VG = np.mean(V, axis=1)  # average over V
		self.AG = np.mean(A, axis=1)  # average over A

		# Rotational
		R = self.R = [Vector(*x[1]) - Vector(*x[0]) for x in X]  # Vector from bone start to end in each frame
		local_axes = [get_principal_axes(r) for r in R]  # Get principal axes for each frame

		# theta_g = (n_frame, 3) of angular rotation about i, j, k for each frame
		# angular rotation about each axis is defined as 0 for the next vector in the cycle
		# i.e. angular rotation about i = 0 for a vector parallel to j
		zero_angles = [[0, 1, 0], [0, 0, 1], [1, 0, 0]]  # definition of 'zero angle' vector for i, j, k
		theta_g = []

		# Compute theta_g in local axes first, where K is the unit vector
		for n_frame in range(len(X) - 1):
			local_ax = local_axes[n_frame]
			# representation as a a single rotation theta about an axis e (https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula)
			a = R[n_frame]  # rotation from one frame...
			b = R[n_frame + 1]  # ...to the next
			if np.array_equal(a, b):
				theta_g += [[0, 0, 0]]  # If no rotation, return 0
			else:
				axis = np.cross(a, b) / (np.linalg.norm(np.cross(a, b)))  # unit vector of omega
				with np.errstate(invalid='raise'):
					try:
						alignment = np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))
						alignment = np.clip(alignment, a_min=-1,
											a_max=1)  # clip between -1 and 1 to deal with rounding errors
						angle = np.arccos(alignment)  # magnitude of theta
					except:
						print((np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))))
						raise ValueError("INVALID ANGLE", a, b)

				theta_g += [axis * angle]

		theta_g = np.array(theta_g)
		self.theta_g = signal.savgol_filter(theta_g, window_length=19, polyorder=2, axis=0)
		self.omega_g = time_deriv(self.theta_g, dt=1 / self.freq)
		self.alpha_g = time_deriv(self.omega_g, dt=1 / self.freq)  # angular acceleration

		self.I_fixed = [la.T @ self.I @ la for la in local_axes]  # compute I in fixed reference frame at each frame

	def get_dynamics(self):
		"""Compute dynamics (F_net, Torque_net) at each frame"""
		self.F_net = [self.mass * a_g for a_g in self.AG]
		self.tau_net = [I_f @ alpha for (I_f, alpha) in zip(self.I_fixed, self.alpha_g)]

class Body(DoubleCylinder):
	"""A unique case of double cylinder, where the (multiple) joints connect at the cylindrical surface at either end.
	These joint attachments are defined by an angle from the i direction normal to the centerline at initialisation.

	Dynamics for the body then must be calculated using a separate set of equations. Define the body such that all
	joints bones go into it rather than out of it (i.e. all input forces are positive on the body) """

	def __init__(self, start_joints, end_joints, all_joint_positions, **cylinder_kwaargs):
		"""From the indices given by start_joints and end_joints, identify a cylinder shape that best fits these
		points on either side, and create that as the cylinder. """

		self.start_joints = start_joints
		self.end_joints = end_joints

		start_pos = Vector(*np.mean(all_joint_positions[40, start_joints], axis=0))
		end_pos = Vector(*np.mean(all_joint_positions[40, end_joints], axis=0))

		length = start_pos > end_pos

		super().__init__(start=None, end=None, length=length, **model.bone_length_definitions["body"](length),
						 **cylinder_kwaargs)

	def get_centre_of_gravity(self, start: 'Vector', end: 'Vector'):
		"""Calculates the centre of gravity based on the displacement from the centerline."""
		centreline_g = 0.5 * (start + end)

		# to find a normal that is closest to z, find N possible equispaced normals, and see which one has the greatest .k product
		normal = (start - end).find_normal()
		N = 20  # number of normals to consider
		all_normals = [normal.rotate_about((start - end).unit(), angle=(n * 2 * np.pi / N)) for n in range(N)]

		idx = np.argmax([v.dot(Vector(0, 0, -1)) for v in all_normals])
		chosen_normal = all_normals[idx]  # choose most downwards normal

		return centreline_g + self.displacement * chosen_normal

	def get_kinematics(self, data):
		"""For body, data is of shape (n_frames, 2, 2, 3), where it is split by rear and front.
		So average across rear and front to get behaviour of centerline, and then run normal get_kinematics"""
		super().get_kinematics(np.mean(data, axis=2))

def weighted_bound_least_squares(A, b, weights=None, bounds=None, **kwargs):
	"""Completes a least squares solve of the equation A x = b, to solve N unknowns from M equations
	where A is an M x N matrix, x is an N x 1 vector, and b is an M x 1 vector.
	Applies weightings to each row to favour certain datapoints. weights is an M x 1 vector.

	Applies bounds where bounds is an M x 2 array. each tuple in the array gives the LB and UB for the given equation"""

	if weights is None: weights = np.ones_like(b)  # If no weight given, equal weight to all
	w = np.array(weights)

	weighted_A, weighted_b = np.array(A) * w[:, np.newaxis], np.array(b) * w  # Apply weights to A, b

	try:
		solve = optimize.lsq_linear(weighted_A, weighted_b, bounds=list(zip(*bounds)), tol=1e-2)
		return solve["x"]

	except np.linalg.LinAlgError as e:
		out = f"SVD did not converge in Lin Least Sq. Printing params: {A}, {b}"
		raise ArithmeticError(out)

class InverseDynamicsSolver:
	"""Through scipy optimisation, Skeleton finds a set of force data that corresponds to the correct kinematic data.
	Takes a skeleton, and the relevant bones and joints, and solves the set of forces that correspond to correct kinematics."""

	def __init__(self, joint_data, target_bones, body_joints, no_torque_joints=None, no_reaction_joints=None,
		foot_joints=None, leg_spring_joints=None, model=Model(),
		freq=50.0, name="output", is_mocap=True):

		for var in [foot_joints, leg_spring_joints, no_reaction_joints, no_torque_joints]:
			if var is None:
				var = []

		self.name = name
		self.freq = freq
		self.n_frames, self.n_joints, _ = joint_data.shape

		self.model = model
		self.is_mocap = is_mocap

		# Preprocess joint data - basic smoothing
		if is_mocap:
			window_length = self.freq // 2
		else:
			window_length = 0.75 * self.freq

		if window_length % 2 == 0: window_length -= 1

		self.T = self.n_frames / self.freq

		self.smooth = lambda X, p=5: signal.savgol_filter(X, window_length=int(window_length), polyorder=p, axis=0)

		p = 5 if self.is_mocap else 2

		self.unsmoothed_data = joint_data  # save unsmoothed data for other uses
		self.joint_pos = self.smooth(joint_data, p=p)
		self.joint_vel = time_deriv(self.joint_pos, 1 / freq)
		self.joint_accel = time_deriv(self.smooth(self.joint_vel), 1 / freq)

		self.foot_joints = foot_joints
		self.body_joints = body_joints
		self.get_foot_joint_from_index = {}  # Identify which foot from the index
		for fj in self.foot_joints:
			for bone, (j1, j2) in target_bones.items():
				if fj in [j1, j2]:
					self.get_foot_joint_from_index[fj] = bone

		self.no_torque_joints = no_torque_joints
		self.no_reaction_joints = no_reaction_joints

		self.target_bones_dict = target_bones  # for use in plotting
		self.target_bones = []

		self.total_mass = 0

		for bone, (joint1, joint2) in target_bones.items():
			# Calculate length using the initial positions of jointA and B.
			# Smoothing functions can cause the issues for the first few frames, so take avg of later frames

			frames = [50, 51, 52, 53, 54, 55, 56]
			n_averaging = len(frames)
			length = 0

			for frame in frames:
				posA = Vector(*self.joint_pos[frame, joint1])
				posB = Vector(*self.joint_pos[frame, joint2])

				if posA.length() == 0 or posB.length() == 0:
					n_averaging -= 1
				else:
					length += posA > posB

			length = length / n_averaging  # avg of all the frames data taken from

			if length == 0:
				print(f"Warning: Error in calculating length of '{bone}'")
				length = 0.01

			b = DoubleCylinder(start=joint1, end=joint2, length=length, name=bone, freq=freq,
							   **self.model.bone_length_definitions["normal"](length))

			self.target_bones.append(b)  # add bone to list
			self.total_mass += b.mass

		self.body = Body(*body_joints, self.joint_pos, freq=freq, name="body")

		self.body.get_kinematics(
			np.stack([self.joint_pos[:, body_joints[0]], self.joint_pos[:, body_joints[1]]], axis=1))
		self.body.get_dynamics()

		self.total_mass += self.body.mass

		# Paw parameters
		m = self.total_mass
		paw_d = self.model.paws
		self.L0_paws = {"front": paw_d["L0_front"] * m, "rear": paw_d["L0_rear"] * m}
		self.k_paws = {"front": paw_d["k_front"] * m, "rear": paw_d["k_rear"] * m ** paw_d["k_rear_prop"]}
		self.c_paws = {"front": paw_d["c_front"] * m, "rear": paw_d["c_rear"] * m}

		# if self.model.equation_weighting['Paw spring'] > 0:
		self.set_paw_equilibrium()

		self.get_dynamics()
		self.leg_spring_joints = leg_spring_joints

		self.calc_leg_lengths()

		self.equation_weighting = model.equation_weighting

	def get_dynamics(self):
		"""Gets dynamics of centre of mass of each bone & body"""
		for bone in self.target_bones:
			bone.get_kinematics(self.joint_pos[:, [bone.start, bone.end]])
			bone.get_dynamics()

		body = self.body
		body.get_kinematics(
			np.stack([self.joint_pos[:, body.start_joints], self.joint_pos[:, body.end_joints]], axis=1))
		body.get_dynamics()

	def calculate_forces(self, n_frame, report_equations=True):
		"""
		Sets up a system of linear equations governing the motion of the skeleton at a given frame.

		These equations are:

		- FREE JOINTS:
		The torques at free joints are zero. Free joints are joints only connected to one bone, on the end of the body eg the feet

		- INERTIA:
		On each bone, the sum of the two joint forces is equal to the mass * acceleration of the bone

		- ROTATION:
		On each bone, the net torque about the bone is equal to the I * alpha_g of the bone

		- BODY:
		The body is set up as a slightly different type of bone, in which it has several joints connected at either end, and its position is dictated by all of those joints.
		See the code for it below, it has its own set of inertial and rotational equations.

		This is set up as a least squares problem Ax = b, where A is a matrix of coefficients to multiply the unknowns by,
		x is the unknowns (in the form [F_1_x, F_1_y, F_1_z, F_2_x, ... T_1, T, ...]
		b is the result of the equations.

		A weighting is also applied to each row to weight the least squares problem (eg to priorities free joint equations)

		The problem also has bounds applied to it. For now, these bounds are simply that foot joint vertical reaction forces are non negative.

		Improvements:
		- Replace the current spinal system with a large non axisymmetric cylinder to represent the body
		- Add a sphere to represent the head

		"""

		# Consult report for explanation of system

		A = []
		b = []
		weights = []  # Collect weightings for each equation as they are added to the system
		equation_weighting = self.equation_weighting

		# Reasonable bounds for each force, and for each torque. Current limits set at 10 * weight for mass, 10 * mass at one metre for torque
		max_force = 3 * self.total_mass * g
		max_torque = 3 * self.total_mass * g
		# bounds can be adjusted further for specific joints (eg no downards reaction at the feet)
		bounds = [(-max_force, max_force)] * (3 * self.n_joints) + [(-max_torque, max_torque)] * (self.n_joints)

		def A_row(vals={}):
			"""Returns a row of 0s length 4 * self.n_joints, with other vectors in any indices in vals.
			vals is a dict of index:vector"""
			row = [0.0] * 4 * self.n_joints
			for index, val in vals.items():
				row[index] = val
			return row

		def add_blank_row():
			A.append(A_row({}))
			b.append(0)
			weights.append(0)

		def add_n_blank_rows(n=1):
			for i in range(n): add_blank_row()

		null, unit, g_vec = Vector(0, 0, 0), Vector(1, 1, 1), Vector(0, 0, -g)

		n_joints = self.n_joints

		def get_index(joint, dimension=0, is_force=True):
			"""Get correct index of D"""
			return (3 * n_joints * int(not is_force)) + ([1, 3][is_force] * joint) + dimension

		# dimension = 0 for x, 1 for y, 2 for z

		# First, add the equations to show that the torques in each of the foot joints are zero
		for no_torque_joint in self.no_torque_joints:
			# Set up the equation 1 * tau_{foot_joint} = 0
			# BOUNDARY CONDITIONS ARE FIXED, RATHER THAN AN ADDITIONAL EQUATION. SO INCORPORATE THEM INTO BOUNDS
			bounds[get_index(no_torque_joint, is_force=False)] = (0, 1e-10)

		for no_reaction_joint in self.no_reaction_joints:  # BC : no reactions
			for dim in [0, 1, 2]:
				bounds[get_index(no_reaction_joint, dimension=dim, is_force=True)] = (0, 1e-10)

		for foot_joint in self.foot_joints:
			## If the feet are a certain amount off the ground for that foot, also assign the reaction forces to be zero
			bone_name = self.get_foot_joint_from_index[foot_joint]

			end = bone_name.split(" ")[1]  # get 'front' or 'rear'
			L0 = self.L0_paws[end]  # get stiffness from 'front' or 'rear' in bone name
			# L0 = self.paw_equilibrium_values[foot_joint]

			k_paw = self.k_paws[end]
			c_paw = self.c_paws[end]

			paw_disp = self.paw_disps[foot_joint][n_frame]

			paw_off_ground = self.joint_pos[n_frame, foot_joint, 2] >= L0  # BC: no reaction in foot off ground
			paw_off_ground = paw_disp == 0

			if paw_off_ground:  # BC: no reaction in foot off ground
				for dim in [0, 1, 2]:
					bounds[get_index(foot_joint, dimension=dim, is_force=True)] = (0, 1e-10)

				add_n_blank_rows(4)  # for consistency of number of eqns

			else:  # If paw near ground, add force due to spring

				height = self.unsmoothed_data[n_frame, foot_joint, 2]

				eps = L0 - height  # min((L0 - height), L0/2)
				eps_dot = self.joint_vel[n_frame, foot_joint, 2]

				F_damp = 0  # c_paw * eps_dot

				if self.model.equation_weighting['Paw spring'] > 0:
					## PAW SPRING MODEL
					eps = paw_disp
					F_spring = k_paw * eps + c_paw * eps_dot

					if foot_joint != 20:
						A.append(A_row({get_index(foot_joint, dimension=2, is_force=True): 1}))
						b.append(F_spring + F_damp)
						weights.append(equation_weighting["Paw spring"])

				if self.model.equation_weighting['Leg spring'] > 0:
					## LEG SPRING MODEL
					K = 3000 if end == "front" else 2000
					for dim in [0, 1, 2]:
						# component = self.leg_vecs[foot_joint][n_frame][dim]
						F_spring = self.leg_disps[foot_joint][n_frame] * K  # * component
						A.append(A_row({get_index(foot_joint, dimension=dim, is_force=True): 1}))
						b.append(F_spring + F_damp)
						weights.append(equation_weighting["Leg spring"])

				# Set bounds for foot joints to only have positive vertical reactions
				bounds[get_index(foot_joint, dimension=2, is_force=True)] = (0, max_force)
				bounds[get_index(foot_joint, dimension=1, is_force=True)] = (0, 1e-10)  # set Fy=0

		for bone in self.target_bones:
			j_1, j_2 = bone.start, bone.end
			x_1, x_2 = bone.X[n_frame]

			# F_1 + F_2 + F_grav = F_net
			F_net = bone.F_net[n_frame]
			for dim in [0, 1, 2]:
				A.append(A_row({get_index(j_1, dim): 1, get_index(j_2, dim): - 1}))
				b.append((F_net - bone.mass * g_vec)[dim])
				weights.append(equation_weighting["Inertial"])

			tau_net = bone.tau_net[n_frame]
			x_g = bone.XG[n_frame]
			r_1, r_2 = (x_1 - x_g), (x_2 - x_g)

			# direction of each T is perpendicular to the bones that the joint is on
			adjacent_1_bone = [b for b in self.target_bones if b.end == j_1 and b != bone]
			if len(adjacent_1_bone) == 1:  # if there is an adjacent bone
				adj_bone = adjacent_1_bone[0]
				T_1_dir = Vector(*r_1).cross((adj_bone.X[n_frame, 1] - adj_bone.XG[n_frame])).unit()
			if len(adjacent_1_bone) == 0 or np.isnan(T_1_dir).any():  # if no adjacent, or if above calc causes error
				T_1_dir = (0, 1, 0)  # Improve later, for now say all torques about y axis

			adjacent_2_bone = [b for b in self.target_bones if b.start == j_2 and b != bone]
			if len(adjacent_2_bone) == 1:  # if there is an adjacent bone
				adj_bone = adjacent_2_bone[0]
				T_2_dir = Vector(*r_2).cross((adj_bone.X[n_frame, 0] - adj_bone.XG[n_frame])).unit()
			if len(adjacent_2_bone) == 0 or np.isnan(T_2_dir).any():  # if no adjacent, or if above calc causes error
				T_2_dir = (0, 1, 0)  # Improve later, for now say all torques about y axis

			for dim in [0, 1, 2]:
				# This loop essentially writes out the following equations into A and b for each dimension (x,y,z):
				# r1 x F1 + r2 x F2 + T1 + T2 = T_net

				# The cross product of r = (x,y,z) and F = (Fx, Fy, Fz) yields (Fz*y - Fy*z, ...)
				# Take the x component, x -> Fz*y - Fy*z
				# Notice that Fy is negative, and Fz is positive. This is always true, that, for the forces, one lower dimension than the current is positive, and one higher is negative (cyclical relations)
				# use this below

				# Get dim above and below, wrapping round for below x and above z
				dim_below = (dim - 1) % 3
				dim_above = (dim + 1) % 3

				coeff_dict = {
					get_index(j_1, dim): 0,
					# eg no effect of F_x in the x directional torque (not relevant statement, only here for readability)
					get_index(j_1, dim_above): - r_1[dim_below],  # eg multiply - z by Fy in the x direction
					get_index(j_1, dim_below): r_1[dim_above],  # eg multiply y by Fz in the x direction

					# Reversed polarity for joint 2 as the desired force is - F2
					get_index(j_2, dim_above): r_2[dim_below],
					get_index(j_2, dim_below): - r_2[dim_above],

					# Add the torques on each joint
					get_index(j_1, is_force=False): T_1_dir[dim],
					get_index(j_2, is_force=False): -T_2_dir[dim]

				}

				A.append(A_row(coeff_dict))
				b.append(tau_net[dim])
				weights.append(equation_weighting["Rotational"])

		### SOLVE FORCES ON BODY. Note body defined so all joint forces/torques on it are positive
		body = self.body
		F_net = body.F_net[n_frame]

		#  BODY INERTIAL FORCES
		for dim in [0, 1, 2]:
			A.append(A_row({get_index(j, dim): 1 for j in self.body.start_joints + self.body.end_joints}))
			b.append((F_net - body.mass * g_vec)[dim])
			weights.append(equation_weighting["Inertial"])

		# BODY ROTATIONAL FORCES - same as for bones
		x_g = body.XG[n_frame]
		tau_net = body.tau_net[n_frame]

		# Improve above later, for now say all torques about y axis
		T_dir = (0, 1, 0)

		for dim in [0, 1, 2]:
			coeff_dict = {}
			for joint in body.start_joints + body.end_joints:
				x_j = self.joint_pos[n_frame, joint]
				r_j = (x_j - x_g)  # position vector to centre

				# Get dim above and below, wrapping round for below x and above z
				dim_below, dim_above = (dim - 1) % 3, (dim + 1) % 3

				coeff_dict[get_index(joint, dim_above)] = -r_j[dim_below]  # eg multiply - z by Fy in the x direction
				coeff_dict[get_index(joint, dim_below)] = r_j[dim_above]  # eg multiply y by Fz in the x direction

				coeff_dict[get_index(joint, is_force=False)] = T_dir[dim]  # Add pure torque of pin

			A.append(A_row(coeff_dict))
			b.append(tau_net[dim])
			weights.append(equation_weighting["Rotational"])

		# print each line of the equations defined by A, b, with the final result
		# Only print variables with both non-zero values, and non-zero coefficients
		if report_equations:
			print(f"----Frame {n_frame}----")
			params = []

			for joint in range(self.n_joints):
				for dim in "xyz":
					params.append(F"F_{joint}_{dim}")  # Add forces by joint

			for joint in range(self.n_joints):
				params.append(F"T_{joint}")  # Add torques by joint

			for n, (coeffs, result) in enumerate(zip(A, b)):
				s = []
				for j, (coeff, param) in enumerate(zip(coeffs, params)):
					if coeff != 0:
						s.append(f"{round(coeff, 3)} * {param}")

				# b_actual = np.dot(A[n], D)
				# pct_error = abs(100 * (b_actual - result) / b_actual)
				if n <= 7:
					print(f"{' + '.join(s)} = {round(result, 3)}")  # ({round(b_actual, 3)}) [{round(pct_error, 2)}%]")

		return A, b, weights, bounds

	def solve_forces(self, report_equations=False, end_frames_disregarded=5, prefix="",
					 save=True):
		"""Solves the forces at each frame for the system, collects them and saves them to .npy files.

		Note: Currently, due to smoothing, the first 5 and last 5 frames are disregarded"""

		self.get_dynamics()
		n_joints = self.n_joints

		if report_equations:
			print("Solving system...")
			print(f"Total mass {round(self.total_mass, 2)} kg.")

		# If dir doesn't exist, make it
		dir = path_join(DataSources.dynamics_data, self.name)
		if self.name not in os.listdir(DataSources.dynamics_data):
			os.mkdir(dir)

		forces, torques = [], []

		f_shape, t_shape = (self.n_joints, 3), (self.n_joints,)
		# Add zeros either end due to not being able to calculate for the first or last 2 frames
		for i in range(end_frames_disregarded):
			forces.append(np.zeros(f_shape))
			torques.append(np.zeros(t_shape))

		calc_forces = []
		calc_torques = []

		progress = tqdm(total=self.n_frames - 2 * end_frames_disregarded)
		for n_frame in range(end_frames_disregarded, self.n_frames - end_frames_disregarded):
			A, b, weights, bounds = self.calculate_forces(n_frame, report_equations=report_equations)

			D = weighted_bound_least_squares(A, b, weights, bounds, rcond=None)

			f, tau = D[:(3 * n_joints)], D[(3 * n_joints):]

			f, tau = f.reshape((n_joints, 3)), tau.reshape((n_joints))

			calc_forces.append(f)
			calc_torques.append(tau)

			progress.update()

		forces[end_frames_disregarded: - end_frames_disregarded] = calc_forces
		torques += calc_torques

		for i in range(end_frames_disregarded):
			forces.append(np.zeros(f_shape))
			torques.append(np.zeros(t_shape))

		if save:
			np.save(path_join(dir, prefix + "forces.npy"), forces)
			np.save(path_join(dir, prefix + "torques.npy"), torques)

		return np.array(forces), np.array(torques)

	def get_com_position(self):
		"""Calculates the position of the centre of mass of the whole system at each timestep"""
		return sum(b.XG * b.mass for b in self.target_bones + [self.body]) / self.total_mass

	def return_equations(self, end_frames_disregarded=5):
		"""For each frame, return the equation vector b"""
		self.get_dynamics()
		bs = []

		for n_frame in range(end_frames_disregarded, self.n_frames - end_frames_disregarded):
			A, b, weights, bounds = self.calculate_forces(n_frame, report_equations=False)
			bs.append(b)

		return np.array(bs)

	def set_paw_equilibrium(self):
		"""Get paw equilibrium from mocap data by finding the drop of the paw.
		This method will work for the current dataset, but is likely not robust, so can be replaced with
		a better method of finding the paw equilibrium at a later date"""

		if self.is_mocap:
			paw_z_heights = self.unsmoothed_data[:, self.foot_joints, 2]

		else:
			paw_z_heights = self.unsmoothed_data[:, self.foot_joints, 2]

		self.paw_disps = {}  # paw joint: displacement over time, for paw spring model

		min_contacts_detected = 3  # minimum requirement to use peak detection mode

		plot = True
		if plot:
			fig, axes = plt.subplots(nrows=2, ncols=2)

		footfall_detector = FootfallDetector(train=False, load=True, name=["smal", "mocap"][self.is_mocap])
		for n, paw in enumerate(self.foot_joints):
			contact_ends_failed = False
			disp = np.zeros((self.n_frames))  # will give eps - the displacement of the paw from equilibrium
			# for when the paw is in contact with the ground

			Z = paw_z_heights[:, n]

			on_ground = footfall_detector.process_clip(Z)
			on_ground_idxs = np.where(on_ground > 0)[0]

			if plot:
				axes[n // 2, n % 2].plot(Z.mean() * (on_ground), color="red", alpha=0.3)

			min_footfall_width = 3  # 3 frames long minimum to count as a footfall
			footfalls = consecutive(on_ground_idxs)
			trigger_height = np.percentile(np.array([Z[ff].max() for ff in footfalls]), 25)  # mean trigger height
			for footfall in footfalls:
				if len(footfall) > min_footfall_width:
					# disp[footfall] = Z[footfall].max() - Z[footfall] # old
					disp[footfall] = np.clip(trigger_height - Z[footfall], a_min=0, a_max=None)

			self.paw_disps[paw] = disp

			if plot:
				ax = axes[n // 2, n % 2]
				ax.plot(Z)

				Z_on_ground = Z.copy()
				Z_on_ground[disp == 0] = np.nan
				ax.plot(Z_on_ground, color="green")
				ax.plot(disp)
				Z_smoothed = self.joint_pos[:, paw, 2]

				ax.set_title(n)

		if plot:
			plt.show(block=False)
			plt.draw()
			plt.pause(1e-8)

	def view_ground_displacements(self, deriv=0):
		"""Plot and show a graph of vertical displacement against frames for each paw - identifying L0 for each paw"""

		fig, axes = plt.subplots(nrows=4)
		for n, j in enumerate(self.foot_joints):
			label = foot_joint_labels[n]
			ax = axes[n]
			if deriv == 0:
				X = self.joint_pos[:, j, 2]
				X_unsmoothed = self.unsmoothed_data[:, j, 2]
				ax.plot(X)
				ax.plot(X_unsmoothed, alpha=.6)
				# ax.axhline(self.paw_equilibrium_values[j], ls = "--")
				ax.axhline(self.L0_paws[label.split(" ")[0]])
			elif deriv == 1:
				ax.plot(self.joint_vel[:, j, 2])

			ax.set_title(label)

		plt.show()

	def view_com_displacements(self, deriv=0):
		"""Plot and show graph of X, Y, and Z motion of CoM of dog.
		If deriv > 0, plot that derivative of the displacement"""

		fig, ax = plt.subplots()
		com_data = self.get_com_position()
		if deriv > 0:
			com_data = nth_time_deriv(com_data, 1 / self.freq, n=deriv)

		for i in [0, 1, 2]:
			ax.plot(com_data[:, i], label="xyz"[i])

		ax.legend()
		plt.show()

	def calc_leg_lengths(self):
		"""Uses the compliant-legged walking model estimation to work out the average length of legs.
		Assume legs are undeformed while off ground. Work out avg distance from leg to COM"""

		self.leg_disps = {}  # length of leg over time for each paw
		self.leg_vecs = {}  # normalised vector of leg spring direction for each paw

		plot = True
		if plot: fig, axes = plt.subplots(nrows=2, ncols=2, sharex="all", sharey="row")

		for n, paw in enumerate(self.foot_joints):
			is_front = n < 2  # Assumes order of f left, f right, r left, r right

			tol = 1e-3
			on_ground = self.paw_disps[paw] > tol
			off_ground = self.paw_disps[paw] <= tol

			# centre_of_rot = self.body.XG[:]#self.body.X[:, int(is_front)]
			# centre_of_rot = self.unsmoothed_data[:, self.body_joints[is_front][n%2]]
			if self.is_mocap:
				centre_of_rot = self.unsmoothed_data[:, self.leg_spring_joints[n]]
				paw_pos = self.unsmoothed_data[:, paw]

			else:
				centre_of_rot = self.unsmoothed_data[:, self.leg_spring_joints[n]]
				paw_pos = self.unsmoothed_data[:, paw]

			X, Z = np.swapaxes(centre_of_rot[:, [0, 2]], 0, 1)  # get X, Z position of CoM
			X_PAW, Z_PAW = np.swapaxes(paw_pos[:, [0, 2]], 0, 1)  # get X, Z position of CoM

			THETA = np.arctan((X_PAW - X) / (Z - Z_PAW))  # angle between spring and vertical

			L = ((X - X_PAW) ** 2 + (Z - Z_PAW) ** 2) ** .5
			L0 = (L).max()

			z_disp = (L - L0) * np.cos(THETA)
			x_disp = (L - L0) * np.sin(THETA)

			# get z displacement by footfall
			disp = np.zeros(self.n_frames)

			# if self.is_mocap:
			for ff in consecutive(np.where(on_ground)[0]):
				if len(ff) < 3: continue  # min width of footfall required
				disp[ff] = z_disp[ff].max() - z_disp[ff]

			# else:
			#     disp = -z_disp

			self.leg_disps[paw] = disp

			if plot:
				ax = axes[n // 2, n % 2]

				# ax.plot(L)
				ax.plot(L - L0)
				ax.plot(disp, color="green")

		if plot:
			plt.tight_layout()
			# plt.show()
			plt.show(block=False)
			plt.draw()
			plt.pause(1e-8)

def norm_kin_data(kin_data, targ_markers=None):
	"""Normalise kinematic data.
	If targ_markers given, normalise so these markers are at desired height"""

	norm_height = 0.4  # 0.635 # fixed to Ally height for now

	# scale so minimum is at (0,0,0)
	for dim in [0, 1, 2]:
		kin_data[:, :, dim] -= kin_data[:, :, dim].min()

	if targ_markers is None:
		kin_data = norm_height * kin_data / np.max(kin_data[:, :, 2])

	elif targ_markers is not None:
		height_target = kin_data[:, targ_markers, 2].mean()
		kin_data = norm_height * kin_data / height_target

	return kin_data

def get_dyn_data(dynamic_src, clip_length, mass, is_mocap=True, target_freq=100):
	"""Loads and returns kinematic data"""

	force_plate_data, force_plate_tdelay = load_force_plate_data(dynamic_src, is_mocap)
	raw_dyn_data = force_plate_data
	raw_dyn_data *= 1 / (mass * 9.81)

	# resample if requested
	if target_freq != freq_forceplate:
		target_frames = int(len(raw_dyn_data) * target_freq / freq_forceplate)
		dyn_data = signal.resample(raw_dyn_data, target_frames)

		# this resampling causes a jumpiness for the periods of zero value. Fix that here:
		tol = 1e-4
		for paw in range(dyn_data.shape[1]):
			# get indices where should be 0
			antifootfalls = consecutive(np.where(raw_dyn_data[:, paw] < tol)[0])
			min_width = 10  # in frames

			for aff in antifootfalls:
				if len(aff) < min_width: continue
				start, end = aff[0] * target_freq / freq_forceplate, aff[-1] * target_freq / freq_forceplate
				# ^ start and end indices, in remapped frame
				dyn_data[int(start):int(end), paw] = 0  # set to 0

		freq = target_freq

	else:
		freq = freq_forceplate
		dyn_data = raw_dyn_data

	frame_delay = int(freq * force_plate_tdelay)
	n_frames_forceplate = int(clip_length * freq)  # number of frames for forceplate to be same time length as mocap

	if frame_delay == 0:
		return dyn_data[:n_frames_forceplate]

	if frame_delay > 0:  # crop forceplate data
		return dyn_data[frame_delay: frame_delay + n_frames_forceplate]  # crop forceplate data to match mocap/SMAL data

	else:  # fdelay <0, pad forceplate data
		return np.pad(dyn_data, ((int(-frame_delay), 0), (0, 0)))[:n_frames_forceplate]

kin_src_to_solver_name = lambda s: s.replace("/", " ").replace(" ", "_").replace(".c3d", "")

def load_solver(kin_src, clip_length, mocap=True, resample_freq=100):
	if mocap:
		joint_data = C3DData(ax=None, src=kin_src, interpolate=True, crop=clip_length,
							 fix_rotations="3 kph" in kin_src)  # only fix rotations for 3 kph for now

	else:
		joint_data = SMALData(kin_src, freq=30, norm=True, crop=clip_length, smooth=True)

	joint_data.resample_at(resample_freq)  ### TRY RESAMPLING DATA TO 100 Hz
	target_bones, body_joints, no_torque_joints, leg_spring_joints = joint_data.generate_skeleton_mapping()

	# Normalise data based on z data, so that the dog is roughly 0.5m high. Also smooth data
	kin_data = np.array(joint_data.all_data)
	kin_data = norm_kin_data(kin_data, targ_markers=leg_spring_joints)

	solver_kwargs = dict(target_bones=target_bones,
						 body_joints=body_joints, no_torque_joints=no_torque_joints,
						 foot_joints=no_torque_joints, leg_spring_joints=leg_spring_joints,
						 freq=joint_data.freq,
						 name=kin_src_to_solver_name(kin_src))

	solver = InverseDynamicsSolver(joint_data=kin_data, **solver_kwargs, is_mocap=mocap)
	print(f"Solver loaded. Mass = {solver.total_mass:.1f} kg.")
	return solver

max_stars_repo_name: OllieBoyne/dog-dynamics

max_issues_repo_name: OllieBoyne/dog-dynamics

max_forks_repo_name: OllieBoyne/dog-dynamics

reason: from scipy


text:
import sympy as sy

x = sy.symbols('x')
y = sy.symbols('y')
rho = sy.sqrt(x**2 + y**2)
phi = sy.atan2(y,x)

f = rho**(2./3)*sy.sin((2./3)*phi)
b = sy.diff(f,x)
d = sy.diff(f,y)
curlB = sy.diff(d, x) - sy.diff(b, y)

Error = 0
for j in range(21):
    for i in range(41):
        if abs(curlB.subs(x,float(i)/20-1).subs(y,float(j)/20).evalf()) > 1e-10:
            print "curl(b) != 0 for x = ", float(i)/100, ' and y = ', float(j)/100
            Error = 1

if Error != 1:
    print "Curl(b) in upper region = 0 for all (x,y)"
else:
    print "Curl(b) in upper region != 0 for all (x,y)"

curlBlower = curlB.xreplace({sy.atan2(y, x):(sy.atan2(y, x)+2*sy.pi)})
Error = 0

for j in range(21):
    for i in range(41):
        if abs(curlBlower.subs(x,float(i)/20-1).subs(y,float(j)/20-1).evalf()) > 1e-10:
            print "curl(b) != 0 for x = ", float(i)/100, ' and y = ', float(j)/100
            Error = 1

if Error != 1:
    print "Curl(b) in lower region = 0 for all (x,y)"
else:
    print "Curl(b) in lower region != 0 for all (x,y)"

max_stars_repo_name: wathen/PhD

max_issues_repo_name: wathen/PhD

max_forks_repo_name: wathen/PhD

reason: import sympy