text:
from tcga_encoder.utils.helpers import *
from tcga_encoder.data.data import *
from tcga_encoder.definitions.tcga import *
#from tcga_encoder.definitions.nn import *
from tcga_encoder.definitions.locations import *
from tcga_encoder.analyses.dna_functions import *

#from tcga_encoder.algorithms import *
import seaborn as sns
from sklearn.manifold import TSNE, locally_linear_embedding
from scipy import stats, special

def view_results( results_location, store, gene, n_permutations, source, method, disease_string, title_str = "", max_nbr = 100, zoom = True ):
  mean_aucs = store["/%s/%s/%s/%s/labels_0/xval_aucs"%(disease_string,gene,source, method)]
  mean_auc = store["/%s/%s/%s/%s/labels_0/xval_aucs"%(disease_string,gene,source, method)].mean()
  var_auc  = store["/%s/%s/%s/%s/labels_0/xval_aucs"%(disease_string,gene,source, method)].var()

  barcodes = store["/%s/%s/%s/%s/labels_0/xval_predictions"%(disease_string,gene,source, method)].index.values
  diseases = np.array( [s.split("_")[0] for s in barcodes])
  u_diseases = np.unique( diseases )

  disease_aucs = store[ "/%s/%s/%s/%s/labels_0/xval_disease_aucs"%(disease_string,gene,source, method)]
  mean_disease_aucs = disease_aucs.mean(1)
  var_disease_aucs = disease_aucs.var(1)

  std_auc  = np.sqrt( var_auc )
  ordered_mean_aucs = store["/%s/%s/%s/%s/labels_0/xval_aucs_elementwise"%(disease_string,gene,source, method)].mean(1).sort_values(ascending=False)
  ordered_source_genes = ordered_mean_aucs.index.values
  ordered_var_aucs = store["/%s/%s/%s/%s/labels_0/xval_aucs_elementwise"%(disease_string,gene,source, method)].loc[ordered_source_genes].var(1)
  order_std_aucs = np.sqrt(ordered_var_aucs)
  D = len(ordered_mean_aucs.values)

  element_aucs = store["/%s/%s/%s/%s/labels_0/xval_aucs_elementwise"%(disease_string,gene,source, method)]
  element_aucs=element_aucs.T
  element_aucs["ALL"] = mean_aucs
  element_aucs=element_aucs.T
  orientation = "horizontal"
  if zoom is True:
    marker = 'o'
  else:
    marker = '.'

  nD = np.minimum( D, max_nbr )

  if orientation == "vertical":
    f1=pp.figure( figsize=(6,16))
  else:
    f1=pp.figure( figsize=(16,6))

  ax11 = f1.add_subplot(111)
  #
  # disease_aucs = []
  # for disease in u_diseases:
  #   disease_aucs.append( store[ "/%s/%s/%s/%s/labels_0/diseases/%s/xval_aucs_elementwise"%(disease_string,gene,source, method, disease)].mean(1) )
  #
  # disease_aucs = pd.concat(disease_aucs,axis=1)
  #pdb.set_trace()
  if orientation == "vertical":
    # for d_idx in range(len(u_diseases)):
   #    disease = u_diseases[d_idx]
   #    results_store[ "/%s/%s/%s/labels_%d/diseases/%s/xval_aucs_elementwise"%(dna_gene,source, method,label_permutation_idx,disease)] = pd.DataFrame( elementwise_aucs_by_disease[d_idx,:,:], index = source_data.columns, columns=xval_columns )
   #    results_store[ "/%s/%s/%s/labels_%d/diseases/%s/xval_aucs"%(dna_gene,source, method,label_permutation_idx,disease)] = pd.DataFrame( np.array(aucs_by_disease).T, index = u_diseases, columns=xval_columns )
   #
    ax11.vlines( mean_disease_aucs.values, 0, nD-1, color='g' )

    for disease in u_diseases:
      aucs =store[ "/%s/%s/%s/%s/labels_0/diseases/%s/xval_aucs_elementwise"%(disease_string,gene,source, method, disease)].mean(1)
      ax11.plot( aucs.loc[ordered_source_genes].values[:nD], nD-np.arange(nD)-1, '.-', mec = 'k', label = "%s"%(disease) )

    #pdb.set_trace()
    ax11.plot( ordered_mean_aucs.values[:nD], nD-np.arange(nD)-1, 'b'+marker+"-", mec = 'k', label = "True" )

    ax11.fill_betweenx( nD-np.arange(nD), \
                        ordered_mean_aucs.values[:nD] + 2*order_std_aucs.values[:nD], \
                        ordered_mean_aucs.values[:nD] - 2*order_std_aucs.values[:nD], facecolor='blue', edgecolor = 'k', alpha=0.5 )
    ax11.plot( ordered_mean_aucs.values[:nD], nD-np.arange(nD)-1, 'b'+marker+"-", mec = 'k', label = "True" )

    ax11.fill_betweenx( nD-np.arange(nD), \
                        mean_auc*np.ones(nD) -2*std_auc, \
                        mean_auc*np.ones(nD) +2*std_auc, facecolor='blue',edgecolor='k', alpha=0.5 )

    ax11.vlines( mean_auc, 0, nD-1, color='b' )
    if zoom is True:
      ax11.set_yticks( nD-1-np.arange(nD) )
      ax11.set_yticklabels( ordered_source_genes[:nD], rotation='horizontal', fontsize=8 )

  else:
    #ax11.fill_between( 2+np.arange(nD), \
    #                    ordered_mean_aucs.values[:nD] + 2*order_std_aucs.values[:nD], \
    #                    ordered_mean_aucs.values[:nD] - 2*order_std_aucs.values[:nD], facecolor='blue', edgecolor = 'k', alpha=0.5 )
    #ax11.plot( np.arange(nD)+2, ordered_mean_aucs.values[:nD], 'b'+marker+"-", mec = 'k', label = "True" )
    ax11.plot( np.arange(nD)+2, ordered_mean_aucs.values[:nD], 'b-', mec = 'k', label = "True" )

    # ax11.fill_between( 1+np.arange(nD), \
    #                     mean_auc*np.ones(nD) -2*std_auc, \
    #                     mean_auc*np.ones(nD) +2*std_auc, facecolor='blue',edgecolor='k', alpha=0.5 )
    #
    # ax11.hlines( mean_auc, 1, nD, color='b' )
    if zoom is True:
      ax11.set_xticks( 2+np.arange(nD) )
      ax11.set_xticklabels( ordered_source_genes[:nD], rotation='vertical', fontsize=8 )

  #
  #pdb.set_trace()
  #ax11.plot( np.ones( len(mean_aucs.values)), mean_aucs.values, 'o', ms=10, color='orange', mec='k', alpha=0.75)
  #ax11.plot( [1], [mean_auc], 'd', color='orchid',mec='orchid' ,ms=30, mew=2, lw=2, alpha=0.75 )
  permutations = []
  combined_permutations = []
  for permutation_idx in range(n_permutations):
    mean_auc_p = store["/%s/%s/%s/%s/labels_%d/xval_aucs"%(disease_string,gene,source, method,permutation_idx+1)].mean()
    combined_permutations.append( mean_auc_p)
  combined_permutations = pd.Series( np.array(combined_permutations), index = np.arange(n_permutations) )

  #permutations.append(combined_permutations )
  for permutation_idx in range(n_permutations):
    mean_auc_p = store["/%s/%s/%s/%s/labels_%d/xval_aucs"%(disease_string,gene,source, method,permutation_idx+1)].mean()
    var_auc_p  = store["/%s/%s/%s/%s/labels_%d/xval_aucs"%(disease_string,gene,source, method, permutation_idx+1)].var()
    std_auc_p  = np.sqrt( var_auc_p )

    mean_aucs = store["/%s/%s/%s/%s/labels_%d/xval_aucs_elementwise"%(disease_string,gene,source, method,permutation_idx+1)].loc[ordered_source_genes].mean(1)

    #permutations.append( store["/%s/%s/%s/%s/labels_%d/xval_aucs_elementwise"%(disease_string,gene,source, method,permutation_idx+1)].loc[ordered_source_genes] )
    permutations.append( mean_aucs )

    # if orientation == "vertical":
    #   ax11.vlines( mean_auc_p, 0, nD-1, color='r' )
    #   #ax11.plot( mean_aucs[:nD], nD-1-np.arange(nD),  'o', color='orange', mec='k', alpha=0.5)
    # else:
    #   ax11.hlines( mean_auc_p, 0, nD-1, color='r' )
    #   #ax11.plot( nD-1-np.arange(nD), mean_aucs[:nD], 'o', color='orange', mec='k', alpha=0.5)
  #

  permutations = pd.concat( permutations,axis=1 )
  permutations = permutations.T
  permutations["ALL"] = combined_permutations
  new_order = ["ALL"]
  new_order.extend(ordered_source_genes[:nD] )
  permutations = permutations.T.loc[new_order]
  element_aucs=element_aucs.loc[new_order]
  print permutations
  #pdb.set_trace()
  correct_labels = store["/%s/%s/%s/%s/labels_%d/xval_aucs_elementwise"%(disease_string,gene,source, method,0)].loc[ordered_source_genes]
  if orientation == "vertical":
    color = dict(boxes='DarkRed', whiskers='DarkOrange',  medians='Red', caps='Black')
    color2 = dict(boxes='DarkBlue', whiskers='DarkBlue',  medians='DarkBlue', caps='Cyan')
    permutations.T.boxplot(ax=ax11,color=color)
    element_aucs.T.boxplot(ax=ax11,color=color2)
  else:
    color = dict(boxes='LightCoral', whiskers='DarkRed',  medians='DarkRed', caps='LightCoral')
    color2 = dict(boxes='SkyBlue', whiskers='DarkBlue',  medians='DarkBlue', caps='SkyBlue')
    permutations.T.plot.box(ax=ax11,color=color,patch_artist=True)
    element_aucs.T.plot.box(ax=ax11,color=color2,patch_artist=True, widths=0.25)
    if zoom is True:
      ax11.set_xticks( 1+np.arange(len(new_order)) )
      ax11.set_xticklabels( new_order, rotation='vertical', fontsize=8 )

  #pdb.set_trace()
  t_tests = []
  for this_gene in ordered_source_genes[:nD]:
    k = n_permutations

    p_value = ( np.sum( correct_labels.loc[this_gene].values.mean() < permutations.loc[this_gene].values ) + 1.0 )/ (k+1.0)
    #t_tests.append( [gene,stats.ttest_ind( permutations.loc[gene].values, correct_labels.loc[gene], equal_var=False )] )
    t_tests.append(p_value)
  #pdb.set_trace()
  pp.grid('on')
  pp.title( "%s using %s of %s with %s mean AUC = %0.3f"%(gene,disease_string, source, method, mean_auc))
  pp.subplots_adjust(bottom=0.2)
  figname1 = os.path.join( HOME_DIR, os.path.dirname(results_location) ) + "/aucs_%s_%s_%s_%s_%s.png"%(gene,source, method,disease_string,title_str)
  f1.savefig(  figname1, dpi=300 )

def main( data_location, results_location, alpha=0.02 ):
  data_path    = os.path.join( HOME_DIR ,data_location ) #, "data.h5" )
  results_path = os.path.join( HOME_DIR, results_location )

  data_filename = os.path.join( data_path, "data.h5")
  fill_filename = os.path.join( results_path, "full_vae_fill.h5" )

  dna_dir = os.path.join( results_path, "dna_prediction_by_tissue" )
  check_and_mkdir(dna_dir)

  aucs_dir = os.path.join( dna_dir, "aucs_viz" )
  check_and_mkdir(aucs_dir)
  print "HOME_DIR: ", HOME_DIR
  print "data_filename: ", data_filename
  print "fill_filename: ", fill_filename

  #dna_genes = ["TP53","APC","PIK3CA"]
  dna_genes = ["APC","TP53","KRAS","PIK3CA","FBXW7","SMAD4","NRAS","ARID1A","ATM","CTNNB1"]
  print "LOADING stores"
  data_store = pd.HDFStore( data_filename, "r" )
  fill_store = pd.HDFStore( fill_filename, "r" )

  Z_train = fill_store["/Z/TRAIN/Z/mu"]
  Z_val = fill_store["/Z/VAL/Z/mu"]

  Z = np.vstack( (Z_train.values, Z_val.values) )
  n_z = Z.shape[1]
  #pdb.set_trace()
  z_names = ["z_%d"%z_idx for z_idx in range(Z.shape[1])]
  Z = pd.DataFrame( Z, index = np.hstack( (Z_train.index.values, Z_val.index.values)), columns = z_names )

  barcodes = np.union1d( Z_train.index.values, Z_val.index.values )

  dna_observed_bcs = data_store["/CLINICAL/observed"][ data_store["/CLINICAL/observed"]["DNA"]==1 ].index.values

  barcodes = np.intersect1d( dna_observed_bcs, barcodes )

  Z=Z.loc[barcodes]
  Z_values = Z.values

  dna = data_store["/DNA/channel/0"].loc[barcodes] #[dna_genes]
  tissues = data_store["/CLINICAL/TISSUE"].loc[barcodes]

  tissue_names = tissues.columns
  tissue_idx = np.argmax( tissues.values, 1 )

  n = len(Z)
  n_tissues = len(tissue_names)
  n_trials = 10
  trial_names = ["r_%d"%(trial_idx) for trial_idx in range(n_trials)]

  aucs_true  = {} #np.ones( (n_tissues,n_z), dtype=float)
  aucs_true_not  = {}
  aucs_random  = {} #np.ones( (n_tissues,n_trials), dtype=float)
  se_auc_true = {}
  se_auc_random=  {}
  se_auc_true_not={}
  true_y = np.ones(n, dtype=int)

  auc_pvalues = {}
  for dna_gene in dna_genes:
    aucs_true[dna_gene] = 0.5*np.ones((n_tissues,n_z))
    aucs_true_not[dna_gene] = 0.5*np.ones((n_tissues,n_z))
    aucs_random[dna_gene] = 0.5*np.ones((n_tissues,n_z,n_trials))
    se_auc_true[dna_gene] = np.ones((n_tissues,n_z))
    se_auc_true_not[dna_gene] = np.ones((n_tissues,n_z))
    se_auc_random[dna_gene] = np.ones(n_tissues)

    auc_pvalues[dna_gene] = np.ones((n_tissues,n_z))

  aucs_true       = OrderedDict()
  aucs_true_not   = OrderedDict()
  se_auc_true     = OrderedDict()
  se_auc_true_not = OrderedDict()
  se_auc_random   = OrderedDict()

  all_auc_series = []
  all_p_series = []
  for t_idx in range(n_tissues):
    tissue_name = tissue_names[t_idx]
    if tissue_name == "gbm" or tissue_name == "meso":
      print "skipping ",tissue_name
      continue
    print "working %s"%(tissue_name)

    t_ids_cohort = tissue_idx == t_idx

    bcs_cohort = barcodes[pp.find(t_ids_cohort)]

    z_cohort   = Z.loc[bcs_cohort]
    dna_cohort = dna.loc[bcs_cohort]
    dna_genes = dna_cohort.sum(0).sort_values(ascending=False)[:100].index.values
    #pdb.set_trace()
    for dna_gene in dna_genes:

      dna_values = dna_cohort[dna_gene]
      n_1 = dna_values.sum()
      n_0 = len(dna_values)-n_1

      if n_1 == 0:
        continue
      true_y = dna_values.values

      print "     %s"%(dna_gene)

      aucs        = np.zeros( n_z, dtype=float )
      p_values        = np.zeros( n_z, dtype=float )
      #aucs_not    = np.zeros( n_z, dtype=float )
      se_aucs     = np.zeros( n_z, dtype=float )
      #se_aucs_not = np.zeros( n_z, dtype=float )
      se_random   = np.zeros( n_z, dtype=float )
      for z_idx in range(n_z):
        z_values = z_cohort.values[:,z_idx]
        #pdb.set_trace()
        if true_y.sum() == 0 or true_y.sum() == len(true_y):
          pdb.set_trace()

        aucs[z_idx]        = roc_auc_score( true_y, z_values )
        aucs[z_idx] = max(aucs[z_idx],1.0-aucs[z_idx])
        se_aucs[z_idx]     = auc_standard_error( aucs[z_idx], n_1, n_0 )
        se_random[t_idx]   = auc_standard_error( 0.5, n_1, n_0 )
        p_values[z_idx] = auc_p_value( aucs[z_idx], 0.5, se_aucs[z_idx], se_random[t_idx] )

      aucs_series = pd.Series( aucs, index = z_names, name=dna_gene)
      p_series = pd.Series( p_values, index = z_names, name=dna_gene)

      all_auc_series.append( [tissue_name,aucs_series] )
      all_p_series.append( [tissue_name,p_series] )

  z_names = all_auc_series[0][1].index.values
  multi_index = []
  for tissue, series  in all_auc_series:
    #tissue = x[0]; series = x[1]
    multi_index.append( (tissue, series.name) )
  #multi_index = np.array(multi_index).T
  index = pd.MultiIndex.from_tuples(multi_index, names=['tissue', 'gene'])

  auc_df = pd.DataFrame( index = index, columns = z_names )
  p_df = pd.DataFrame( index = index, columns = z_names )

  for tissue, series  in all_auc_series:
    gene = series.name
    auc_df.loc[tissue,gene] = series
  for tissue, series  in all_p_series:
    gene = series.name
    p_df.loc[tissue,gene] = series

  auc_df.to_csv( dna_dir + "/auc_values.csv")
  p_df.to_csv( dna_dir + "/p_values.csv")

  f=pp.figure()
  ax1 = f.add_subplot(121)
  ax2 = f.add_subplot(122)
  ax1.hist( p_df.values.flatten(), bins = np.linspace(0,0.5,21), lw=2, histtype="step", normed=True, range=(0,0.5) )
  ax2.hist( auc_df.values.flatten(), bins = np.linspace(0.5,1,21), lw=2, histtype="step", normed=True, range=(0.5,1) )
  pp.savefig( dna_dir + "/p_values_and_auc.png", fmt = "png", dpi=300)
  f=pp.figure()
  ax1 = f.add_subplot(121)
  ax2 = f.add_subplot(122)
  ax1.hist( p_df.values.flatten(), bins = np.linspace(0,0.5,61), lw=2, histtype="step", normed=True, range=(0,0.5) )
  ax2.hist( auc_df.values.flatten(), bins = np.linspace(0.5,1,61), lw=2, histtype="step", normed=True, range=(0.5,1) )
  pp.savefig( dna_dir + "/p_values_and_auc2.png", fmt = "png", dpi=300)

  f=pp.figure()
  for z in auc_df.columns:
    pp.semilogy( auc_df[z].values.flatten(), p_df[z].values.flatten(), '.' )
  pp.xlabel('AUC');pp.ylabel('p-value')
  pp.savefig( dna_dir + "/p_values_and_auc_scatter.png", fmt = "png", dpi=300)
  pp.close('all')
  return auc_df, p_df
  #pdb.set_trace()
  # sorted_pan = pd.concat(sorted_pan,axis=1)
  # sorted_pan.columns = dna_genes
  # sorted_pan.to_csv( dna_dir + "/pan_logpvals_z_for_dna.csv" )
  #
  # sorted_pan_sig = pd.concat(sorted_pan_sig,axis=1)
  # sorted_pan_sig.columns = dna_genes
  # sorted_pan_sig.to_csv( dna_dir + "/pan_sig_z_for_dna.csv" )
  # size1 = max( int( n_z*size_per_unit ), 12 )
  # size2 = max( int( len(dna_genes)*size_per_unit ), 12 )
  # f = sns.clustermap( sorted_pan.T, figsize=(size1,size2) )
  #
  # pp.setp(f.ax_heatmap.yaxis.get_majorticklabels(), rotation=0)
  # pp.setp(f.ax_heatmap.xaxis.get_majorticklabels(), rotation=90)
  # pp.setp(f.ax_heatmap.yaxis.get_majorticklabels(), fontsize=12)
  # pp.setp(f.ax_heatmap.xaxis.get_majorticklabels(), fontsize=12)
  # pp.savefig( dna_dir + "/z_for_dna_clustermap_logpval.png", format="png")
  #
  # f = sns.clustermap( sorted_pan_sig.T, figsize=(size1,size2) )
  #
  # pp.setp(f.ax_heatmap.yaxis.get_majorticklabels(), rotation=0)
  # pp.setp(f.ax_heatmap.xaxis.get_majorticklabels(), rotation=90)
  # pp.setp(f.ax_heatmap.yaxis.get_majorticklabels(), fontsize=12)
  # pp.setp(f.ax_heatmap.xaxis.get_majorticklabels(), fontsize=12)
  # pp.savefig( dna_dir + "/z_for_dna_clustermap_sig.png", format="png")

  #pp.close('all')

if __name__ == "__main__":

  data_location = sys.argv[1]
  results_location = sys.argv[2]

  a, p = main( data_location, results_location )
  # z_names = a[0][1].index.values
  # multi_index = []
  # for tissue, series  in a:
  #   #tissue = x[0]; series = x[1]
  #   multi_index.append( (tissue, series.name) )
  # #multi_index = np.array(multi_index).T
  # index = pd.MultiIndex.from_tuples(multi_index, names=['tissue', 'gene'])
  #
  # auc_df = pd.DataFrame( index = index, columns = z_names )
  # p_df = pd.DataFrame( index = index, columns = z_names )
  #
  # for tissue, series  in a:
  #   gene = series.name
  #   auc_df.loc[tissue,gene] = series
  # for tissue, series  in p:
  #   gene = series.name
  #   p_df.loc[tissue,gene] = series


text:
import numpy as np
import pytest

from cblearn.datasets import fetch_food_similarity

@pytest.mark.remote_data
def test_fetch_food(tmp_path):
    data_home = tmp_path / 'cblearn_datasets'
    bunch = fetch_food_similarity(data_home=data_home, shuffle=False)

    assert bunch.data.shape == (190376, 3)
    assert bunch.image_names.shape == (100, )
    assert bunch.image_names[bunch.data[0, 0]] == 'images/214649bfd7ea489b8daf588e6fed45aa.jpg'

    triplets = fetch_food_similarity(data_home=data_home, shuffle=False, return_triplets=True)
    np.testing.assert_equal(bunch.data, triplets)

ext: py

max_stars_repo_name: JFHoelscher/cblearn

max_issues_repo_name: JFHoelscher/cblearn

max_forks_repo_name: JFHoelscher/cblearn

reason: import numpy


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

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
from tensor2tensor.layers import discretization
import tensorflow as tf

class DiscretizationTest(tf.test.TestCase):
  """Tests for discretization layers."""

  def setUp(self):
    tf.set_random_seed(1234)
    np.random.seed(123)

  def testBitToIntZeros(self):
    x_bit = tf.zeros(shape=[1, 10], dtype=tf.float32)
    x_int = tf.zeros(shape=[1], dtype=tf.int32)
    diff = discretization.bit_to_int(x_bit, num_bits=10) - x_int
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      d = sess.run(diff)
      self.assertEqual(d, 0)

  def testBitToIntOnes(self):
    x_bit = tf.ones(shape=[1, 3], dtype=tf.float32)
    x_int = 7 * tf.ones(shape=[1], dtype=tf.int32)
    diff = discretization.bit_to_int(x_bit, num_bits=3) - x_int
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      d = sess.run(diff)
      self.assertEqual(d, 0)

  def testIntToBitZeros(self):
    x_bit = tf.zeros(shape=[1, 10], dtype=tf.float32)
    x_int = tf.zeros(shape=[1], dtype=tf.int32)
    diff = discretization.int_to_bit(x_int, num_bits=10) - x_bit
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      d = sess.run(diff)
      self.assertTrue(np.all(d == 0))

  def testIntToBitOnes(self):
    x_bit = tf.ones(shape=[1, 3], dtype=tf.float32)
    x_int = 7 * tf.ones(shape=[1], dtype=tf.int32)
    diff = discretization.int_to_bit(x_int, num_bits=3) - x_bit
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      d = sess.run(diff)
      self.assertTrue(np.all(d == 0))

  def testProjectHidden(self):
    hidden_size = 60
    block_dim = 20
    num_blocks = 3
    x = tf.zeros(shape=[1, hidden_size], dtype=tf.float32)
    projection_tensors = tf.random_normal(
        shape=[num_blocks, hidden_size, block_dim], dtype=tf.float32)
    x_projected = discretization.project_hidden(x, projection_tensors,
                                                hidden_size, num_blocks)
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      x_projected_eval = sess.run(x_projected)
      self.assertEqual(np.shape(x_projected_eval), (1, num_blocks, block_dim))
      self.assertTrue(np.all(x_projected_eval == 0))

  def testSliceHiddenZeros(self):
    hidden_size = 60
    block_dim = 20
    num_blocks = 3
    x = tf.zeros(shape=[1, hidden_size], dtype=tf.float32)
    x_sliced = discretization.slice_hidden(x, hidden_size, num_blocks)
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      x_sliced_eval = sess.run(x_sliced)
      self.assertEqual(np.shape(x_sliced_eval), (1, num_blocks, block_dim))
      self.assertTrue(np.all(x_sliced_eval == 0))

  def testSliceHiddenOnes(self):
    hidden_size = 60
    block_dim = 20
    num_blocks = 3
    x = tf.ones(shape=[1, hidden_size], dtype=tf.float32)
    x_sliced = discretization.slice_hidden(x, hidden_size, num_blocks)
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      x_sliced_eval = sess.run(x_sliced)
      self.assertEqual(np.shape(x_sliced_eval), (1, num_blocks, block_dim))
      self.assertTrue(np.all(x_sliced_eval == 1))

  def testNearestNeighbors(self):
    x = tf.constant([[0, 0.9, 0], [0.8, 0., 0.]], dtype=tf.float32)
    x = tf.expand_dims(x, axis=0)
    means = tf.constant(
        [[1, 0, 0], [0, 1, 0], [0, 0, 1], [9, 9, 9]], dtype=tf.float32)
    means = tf.stack([means, means], axis=0)
    x_means_hot = discretization.nearest_neighbor(x, means, block_v_size=4)
    x_means_hot_test = np.array([[0, 1, 0, 0], [1, 0, 0, 0]])
    x_means_hot_test = np.expand_dims(x_means_hot_test, axis=0)
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      x_means_hot_eval = sess.run(x_means_hot)
      self.assertEqual(np.shape(x_means_hot_eval), (1, 2, 4))
      self.assertTrue(np.all(x_means_hot_eval == x_means_hot_test))

  def testGetVQBottleneck(self):
    bottleneck_bits = 2
    bottleneck_size = 2**bottleneck_bits
    hidden_size = 3
    means, _, ema_count = discretization.get_vq_bottleneck(bottleneck_size,
                                                           hidden_size)
    assign_op = means.assign(tf.zeros(shape=[bottleneck_size, hidden_size]))
    means_new, _, _ = discretization.get_vq_bottleneck(bottleneck_size,
                                                       hidden_size)
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      sess.run(assign_op)
      self.assertTrue(np.all(sess.run(means_new) == 0))
      self.assertTrue(np.all(sess.run(ema_count) == 0))

  def testVQNearestNeighbors(self):
    x = tf.constant([[0, 0.9, 0], [0.8, 0., 0.]], dtype=tf.float32)
    means = tf.constant(
        [[1, 0, 0], [0, 1, 0], [0, 0, 1], [9, 9, 9]], dtype=tf.float32)
    x_means_hot, _ = discretization.vq_nearest_neighbor(x, means)
    x_means_hot_test = np.array([[0, 1, 0, 0], [1, 0, 0, 0]])
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      x_means_hot_eval = sess.run(x_means_hot)
      self.assertEqual(np.shape(x_means_hot_eval), (2, 4))
      self.assertTrue(np.all(x_means_hot_eval == x_means_hot_test))

  def testVQDiscreteBottleneck(self):
    x = tf.constant([[0, 0.9, 0], [0.8, 0., 0.]], dtype=tf.float32)
    x_means_hot, _ = discretization.vq_discrete_bottleneck(x, bottleneck_bits=2)
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      x_means_hot_eval = sess.run(x_means_hot)
      self.assertEqual(np.shape(x_means_hot_eval), (2, 4))

  def testVQDiscreteUnbottlenck(self):
    x = tf.constant([[1, 0, 0, 0], [0, 0, 1, 0]], dtype=tf.int32)
    x_means = discretization.vq_discrete_unbottleneck(x, hidden_size=3)
    with self.test_session() as sess:
      tf.global_variables_initializer().run()
      x_means_eval = sess.run(x_means)
      self.assertEqual(np.shape(x_means_eval), (2, 3))

if __name__ == '__main__':
  tf.test.main()


text:
import os, sys, glob, scipy
import numpy as np
import pandas as pd
from scipy.stats import spearmanr

# Arguments
screen = str(sys.argv[1])
cond = str(sys.argv[2])
nparcel = int(sys.argv[3])
parcel = int(sys.argv[4])
metric_model = str(sys.argv[5])
metric_brain = str(sys.argv[6])
permutation_method = str(sys.argv[7])
n_permute = int(sys.argv[8])

# Select participants
base_dir = 'yourprojectfolderhere'
oldSubs = np.array(pd.read_csv(os.path.join(base_dir,'Data/subjectsIncluded_batch1.csv'),header=None,index_col=None))[0]
newSubs = np.array(pd.read_csv(os.path.join(base_dir,'Data/subjectsIncluded_batch2.csv'),header=None,index_col=None))[0]
fullSample = np.hstack([oldSubs,newSubs])
subNums = fullSample
params = pd.read_csv(os.path.join(base_dir,'Results/2.Behavior-ClusterModel/ParticipantClustering.csv'),index_col=0)
subIndices = np.where(pd.DataFrame(fullSample).isin(subNums))[0]

# Compute model distance
modelDist = scipy.spatial.distance.pdist(params.loc[params['sub'].isin(subNums),['theta','phi']],metric=metric_model)

# Data
brainDist = pd.read_csv(os.path.join(base_dir,'Results/3.fMRI-ISRSA/BrainDist/BrainDist_%s_%s%s_roi%i.csv' % (
				metric_brain,screen,cond,nparcel)))
brainDist = brainDist.iloc[:,parcel]
brainDist_mat = pd.DataFrame(scipy.spatial.distance.squareform(brainDist)).iloc[subIndices,subIndices]
brainDist_vec = brainDist_mat.values[np.triu_indices(len(subIndices),k=1)]

# Operations
corr = spearmanr(modelDist,brainDist_vec)[0]

if n_permute is None:
	pval = spearmanr(modelDist,brainDist_vec)[1]
elif isinstance(n_permute, int):
    perm = []
    if permutation_method == 'vector':
        for p in range(n_permute):
            perm.append(spearmanr(np.random.permutation(modelDist),brainDist_vec)[0])
    if corr>=0:
        perm_p = np.mean(perm>=corr)
    else:
        perm_p = np.mean(perm<=corr)
    pval = perm_p

out = pd.DataFrame([[corr,pval]],columns=['r','p'])
pathCur = os.path.join(base_dir,
        'Results/3.fMRI-ISRSA/IS-RSA/IS-RSA_nparcel-%i_perm-%s_%s%s/'%(
        nparcel,permutation_method,screen,cond))
# if ~os.path.exists(pathCur):
#     os.mkdir(pathCur)
out.to_csv(os.path.join(pathCur,'parcel%03d.csv'%(parcel)))

reason: import numpy,from scipy


text:
#!/usr/bin/python

import sys
import time
import threading
import numpy
import string
import copy
from scipy.optimize import curve_fit
from math import sqrt,exp,log,pi,acos,atan,cos,asin

def g_CIMP(x):

	x=x[0,:]
	g=numpy.zeros(len(x))
	#g=2.691*(1-0.2288964/x)*1/((1+0.16*x)*(1+1.35*numpy.exp(-x/0.2)))

	if (max(x)<1):
		g=-18.743261164767357+101.6507221241339*x**(0.33333333)-104.59646433814892*numpy.sqrt(x)+33.73393945878933*x-10.325598001906716*x**(1.5)
	elif (min(x)>1):
		g=1.1976693536243692*(1-0.2660904859953754/x)*1/((1+0.04920104690300144*x)*(1-0.5874697493344921*numpy.exp(-x*0.09913039025775051)))
	else:
		g=2.691*(1-0.2288964/x)*1/((1+0.16*x)*(1+1.35*numpy.exp(-x/0.2)))

	return g

def gauss_function(x, a, x0, sigma):
    return a*numpy.exp(-(x-x0)**2/(2*sigma**2))

def calc_sigma(pos,profile):

	max_val=max(profile)
	min_val=min(profile)
	profile=(profile-min_val)/max_val
	aux=profile*pos
	ycm=numpy.sum(aux)/numpy.sum(profile)
	aux=profile*(pos-ycm)**2
	yvar=sqrt(numpy.sum(aux)/numpy.sum(profile))

	return (ycm,yvar)

def Calc_Growth(twiss,param):

	#Define parameters
	brel=sqrt(1-1/param['gamma']**2)
	ex=param['exi']
	ey=param['eyi']
	ss=param['ssi']
	sp=param['spi']

	#Define twiss arrays

	s=numpy.zeros(len(twiss))
	betax=numpy.zeros(len(twiss))
	alphax=numpy.zeros(len(twiss))
	betay=numpy.zeros(len(twiss))
	alphay=numpy.zeros(len(twiss))
	Dx=numpy.zeros(len(twiss))
	Dpx=numpy.zeros(len(twiss))
	Dy=numpy.zeros(len(twiss))
	Dpy=numpy.zeros(len(twiss))

	s=twiss[:,0]
	betax=twiss[:,2]
	alphax=twiss[:,3]
	betay=twiss[:,6]
	alphay=twiss[:,7]
	Dx=twiss[:,4]
	Dpx=twiss[:,5]
	Dy=twiss[:,8]
	Dpy=twiss[:,9]

	#Calculate the parameters
	A=param['cluz']*param['Np']*param['r0']**2/(64*numpy.pi**2*brel**3*param['gamma']**4*ex*ey*ss*sp)
	logCIMP=numpy.log(param['gamma']**2*ex*numpy.sqrt(betay*ey)/(param['r0']*betax))
	Hx=1/betax*[Dx**2+(betax*Dpx+alphax*Dx)**2]
	Hy=1/betay*[Dy**2+(betay*Dpy+alphay*Dy)**2]
	SigH=numpy.sqrt(1/sp**2+Hx/ex+Hy/ey)**(-1)
	aCIMP=SigH/param['gamma']*numpy.sqrt(betax/ex)
	bCIMP=SigH/param['gamma']*numpy.sqrt(betay/ey)

	#Calculate Function g
	g_ab=g_CIMP(aCIMP/bCIMP)
	g_ba=g_CIMP(bCIMP/aCIMP)

	#Saves values for the ration a/b and b/a
	#f=open('RatioAB.txt','w')
	#for j in range(len(aCIMP[0,:])):
	#	f.write(str(aCIMP[0,j]/bCIMP[0,j])+'\t\t'+str(bCIMP[0,j]/aCIMP[0,j])+'\n')
	#f.close()

	#Calculate Growth Rates
	fp=A*logCIMP*(SigH**2/sp**2)*(g_ba/aCIMP+g_ab/bCIMP)
	fx=A*logCIMP*(-aCIMP*g_ba+Hx*SigH**2/ex*(g_ba/aCIMP+g_ab/bCIMP))
	fy=A*logCIMP*(-bCIMP*g_ab+Hy*SigH**2/ey*(g_ba/aCIMP+g_ab/bCIMP))

	#Integrate along the s coordinate
	invTp=2*pi**(3.0/2.0)*numpy.trapz(fp,s)
	invTx=2*pi**(3.0/2.0)*numpy.trapz(fx,s)
	invTy=2*pi**(3.0/2.0)*numpy.trapz(fy,s)

	#Calculate growth
	Tp=invTp
	Tx=invTx
	Ty=invTy

	return (Tx,Ty,Tp)

# Fucntion that iterates emittances for the case with no harmonic system (simple calculation of the bunch length)
def Iterate_emittances(twiss,param):

	#Define differences
	i=1
	time=0
	diff1=1
	diff2=1
	diff3=1
	diff4=1
	difftot=diff1+diff2+diff3+diff4

	#Calculate U0
	U0=param['Cgamma']/(2*pi)*(param['En']/1e+9)**4*param['I2']*1e+9
	#print U0

	#Calculate damping partition numbers
	Jx=1-param['I4']/param['I2']
	Jy=1
	Jp=2+param['I4']/param['I2']
	#print Jx,Jy,Jp

	# Caluclate damping times
	taux=(2*param['En']*param['C'])/(Jx*U0*param['cluz'])
	tauy=(2*param['En']*param['C'])/(Jy*U0*param['cluz'])
	taup=(2*param['En']*param['C'])/(Jp*U0*param['cluz'])
	#print taux,tauy,taup

	#Define step for iteration
	tt=taux/5

	# Synchrotron tune
	Qs0=sqrt(param['ap']*param['hh']*sqrt(param['Vrf']**2-U0**2)/(2*pi*param['En']))

	#Cretaes an array that's a subgroup of param
	inter={}
	inter['exi']=param['ex0']
	inter['eyi']=(param['k_dw']+param['k_beta'])*param['ex0']
	inter['ssi']=param['ss0']
	inter['spi']=param['sp0']
	inter['gamma']=param['gamma']
	inter['r0']=param['r0']
	inter['Np']=param['Np']
	inter['cluz']=param['cluz']

	while (difftot>10**(-7)):
		(Tx,Ty,Tp)=Calc_Growth(twiss,inter)
		Tx=float(Tx)/param['C']
		Ty=float(Ty)/param['C']
		Tp=float(Tp)/param['C']
		#print Tx,Ty,Tp

		exx=(-param['ex0']+exp(2*tt*(Tx-1/taux))*(param['ex0']+inter['exi']*(-1+Tx*taux)))/(-1+Tx*taux)

		eyy=(-(param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+exp(2*tt*(Ty-1/tauy))*((param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+inter['eyi']*(-1+Ty*tauy)))/(-1+Ty*tauy)

		spp=(-param['sp0']+exp(tt*(Tp-1/taup))*(param['sp0']+inter['spi']*(-1+Tp*taup)))/(-1+Tp*taup)
		# Accelerating cavity system only

		sss=inter['spi']*param['C']*sqrt(param['ap']*param['En']/(2*pi*param['hh']*(param['Vrf']**2-U0**2)**0.5));
		#print exx,eyy,spp,sss

		diff1=abs(exx-inter['exi'])/inter['exi']
		diff2=abs(eyy-inter['eyi'])/inter['eyi']
		diff3=abs(spp-inter['spi'])/inter['spi']
		diff4=abs(sss-inter['ssi'])/inter['ssi']
		difftot=diff1+diff2+diff3+diff4
		#print difftot

		inter['exi']=exx;

		inter['eyi']=eyy;

		inter['spi']=spp;

		inter['ssi']=sss;

		time=i*tt;
		i=i+1

	return (exx,eyy,spp,sss)

# Function that iterates emittances using the results from tracking to calculate bunch length
def Iterate_emittances3HC(twiss,param,phimain,Vmain,phiharm,Vharm):

	#Define differences
	i=1
	time=0
	diff1=1
	diff2=1
	diff3=1
	diff4=1
	difftot=diff1+diff2+diff3+diff4

	#Calculate U0
	U0=param['Cgamma']/(2*pi)*(param['En']/1e+9)**4*param['I2']*1e+9

	#Calculate synchronous phase
	Phi_sync_nat=asin(U0/param['Vrf'])

	#Calculate damping partition numbers
	Jx=1-param['I4']/param['I2']
	Jy=1
	Jp=2+param['I4']/param['I2']
	#print Jx,Jy,Jp

	# Caluclate damping times
	taux=(2*param['En']*param['C'])/(Jx*U0*param['cluz'])
	tauy=(2*param['En']*param['C'])/(Jy*U0*param['cluz'])
	taup=(2*param['En']*param['C'])/(Jp*U0*param['cluz'])
	#print taux,tauy,taup

	#Define step for iteration
	tt=taux/5

	#RF frequency
	w_rf =2*pi*(param['hh']*param['cluz']/param['C']-param['Detune0'])	#Generator Frequency

	#Creates arrays for 3HC calculation
	posz=numpy.zeros(5000)
	perfil=numpy.zeros(5000)
	pot=numpy.zeros(5000)

	#Define longitudinal scale array
	posz=numpy.arange(0,5000.)/10-250 # in milimiters

	#Cretaes an array that's a subgroup of param
	inter={}
	inter['exi']=param['ex0']
	inter['eyi']=(param['k_dw']+param['k_beta'])*param['ex0']
	inter['spi']=param['sp0']
	inter['gamma']=param['gamma']
	inter['r0']=param['r0']
	inter['Np']=param['Np']
	inter['cluz']=param['cluz']

	pot=1/(param['En']*param['C'])*param['cluz']/w_rf*(Vmain*1e3*(cos(Phi_sync_nat-phimain)-numpy.cos(posz/1000*w_rf/param['cluz']+Phi_sync_nat-phimain))+Vharm*1e3/param['mharm']*(cos(param['mharm']*pi-phiharm)-numpy.cos(param['mharm']*posz/1000*w_rf/param['cluz']+param['mharm']*pi-phiharm)))-1/(param['En']*param['C'])*U0*posz/1000
	perfil=numpy.exp(-pot/(param['ap']*param['sp0']**2))
	(pos0,sigma_mm)=calc_sigma(posz,perfil)
	inter['ssi']=sigma_mm/1000

	while (difftot>10**(-7)):
		(Tx,Ty,Tp)=Calc_Growth(twiss,inter)
		Tx=float(Tx)/param['C']
		Ty=float(Ty)/param['C']
		Tp=float(Tp)/param['C']
		#print Tx,Ty,Tp

		exx=(-param['ex0']+exp(2*tt*(Tx-1/taux))*(param['ex0']+inter['exi']*(-1+Tx*taux)))/(-1+Tx*taux)

		eyy=(-(param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+exp(2*tt*(Ty-1/tauy))*((param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+inter['eyi']*(-1+Ty*tauy)))/(-1+Ty*tauy)

		spp=(-param['sp0']+exp(tt*(Tp-1/taup))*(param['sp0']+inter['spi']*(-1+Tp*taup)))/(-1+Tp*taup)

		#Calculate bunch length according to the RF potential (Main RF + 3HC)
		pot=1/(param['En']*param['C'])*param['cluz']/w_rf*(Vmain*1e3*(cos(Phi_sync_nat-phimain)-numpy.cos(posz/1000*w_rf/param['cluz']+Phi_sync_nat-phimain))+Vharm*1e3/param['mharm']*(cos(param['mharm']*pi-phiharm)-numpy.cos(param['mharm']*posz/1000*w_rf/param['cluz']+param['mharm']*pi-phiharm)))-1/(param['En']*param['C'])*U0*posz/1000
		perfil=numpy.exp(-pot/(param['ap']*spp**2))
		(pos0,sigma_mm)=calc_sigma(posz,perfil)
		sss=sigma_mm/1000

		#print exx,eyy,spp,sss

		diff1=abs(exx-inter['exi'])/inter['exi']
		diff2=abs(eyy-inter['eyi'])/inter['eyi']
		diff3=abs(spp-inter['spi'])/inter['spi']
		diff4=abs(sss-inter['ssi'])/inter['ssi']
		difftot=diff1+diff2+diff3+diff4
		#print difftot

		inter['exi']=exx;

		inter['eyi']=eyy;

		inter['spi']=spp;

		inter['ssi']=sss;

		time=i*tt;
		i=i+1

	return (exx,eyy,spp,sss)

# Function that iterates emittances for the case with no harmonic system (simple calculation of the bunch length) but
# takes into account the longitudinal growth rate due to microwave instability
def Iterate_emittancesMW(twiss,param,sigS,Curr,GT):

	#Define differences
	i=1
	time=0
	diff1=1
	diff2=1
	diff3=1
	diff4=1
	difftot=diff1+diff2+diff3+diff4

	#Calculate U0
	U0=param['Cgamma']/(2*pi)*(param['En']/1e+9)**4*param['I2']*1e+9
	#print U0

	#Calculate damping partition numbers
	Jx=1-param['I4']/param['I2']
	Jy=1
	Jp=2+param['I4']/param['I2']
	#print Jx,Jy,Jp

	# Caluclate damping times
	taux=(2*param['En']*param['C'])/(Jx*U0*param['cluz'])
	tauy=(2*param['En']*param['C'])/(Jy*U0*param['cluz'])
	taup=(2*param['En']*param['C'])/(Jp*U0*param['cluz'])
	#print taux,tauy,taup

	#Define step for iteration
	tt=taux/5

	# Synchrotron tune
	Qs0=sqrt(param['ap']*param['hh']*sqrt(param['Vrf']**2-U0**2)/(2*pi*param['En']))

	# Define the interpolation function for Microwave Instability
	microwave=interp2d(sigS,Curr,GT,kind='linear')

	#Cretaes an array that's a subgroup of param
	inter={}
	inter['exi']=param['ex0']
	inter['eyi']=(param['k_dw']+param['k_beta'])*param['ex0']
	inter['ssi']=param['ss0']
	inter['spi']=param['sp0']
	inter['gamma']=param['gamma']
	inter['r0']=param['r0']
	inter['Np']=param['Np']
	inter['cluz']=param['cluz']
	sss=param['ss0']

	while (difftot>10**(-7)):
		#Add the Microwave growth rate to the longitudinal plane
		DTp=microwave(sss,param['Np'])
		#print DTp

		(Tx,Ty,Tp)=Calc_Growth(twiss,inter)
		Tx=float(Tx)/param['C']
		Ty=float(Ty)/param['C']
		Tp=float(Tp)/param['C']+DTp

		exx=(-param['ex0']+exp(2*tt*(Tx-1/taux))*(param['ex0']+inter['exi']*(-1+Tx*taux)))/(-1+Tx*taux)

		#eyy=(-param['ey0']+exp(2*tt*(Ty-1/tauy))*(param['ey0']+inter['eyi']*(-1+Ty*tauy)))/(-1+Ty*tauy)
		eyy=(-(param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+exp(2*tt*(Ty-1/tauy))*((param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+inter['eyi']*(-1+Ty*tauy)))/(-1+Ty*tauy)

		spp=(-param['sp0']+exp(tt*(Tp-1/taup))*(param['sp0']+inter['spi']*(-1+Tp*taup)))/(-1+Tp*taup)
		# Accelerating cavity system only

		sss=inter['spi']*param['C']*sqrt(param['ap']*param['En']/(2*pi*param['hh']*(param['Vrf']**2-U0**2)**0.5));

		diff1=abs(exx-inter['exi'])/inter['exi']
		diff2=abs(eyy-inter['eyi'])/inter['eyi']
		diff3=abs(spp-inter['spi'])/inter['spi']
		diff4=abs(sss-inter['ssi'])/inter['ssi']
		difftot=diff1+diff2+diff3+diff4
		#print difftot

		inter['exi']=exx;

		inter['eyi']=eyy;

		inter['spi']=spp;

		inter['ssi']=sss;

		time=i*tt;
		i=i+1

	return (exx,eyy,spp,sss)

# Function that iterates emittances using the results from tracking to calculate bunch length
# and also takes into account the longitudinal growth rate due to microwave instability
def Iterate_emittances3HC_MW(twiss,param,phimain,Vmain,phiharm,Vharm,sigS,Curr,GT):

	#Definde differences
	i=1
	time=0
	diff1=1
	diff2=1
	diff3=1
	diff4=1
	difftot=diff1+diff2+diff3+diff4

	#Calculate U0
	U0=param['Cgamma']/(2*pi)*(param['En']/1e+9)**4*param['I2']*1e+9

	#Calculate synchronous phase
	Phi_sync_nat=asin(U0/param['Vrf'])

	#Calculate damping partition numbers
	Jx=1-param['I4']/param['I2']
	Jy=1
	Jp=2+param['I4']/param['I2']
	#print Jx,Jy,Jp

	# Caluclate damping times
	taux=(2*param['En']*param['C'])/(Jx*U0*param['cluz'])
	tauy=(2*param['En']*param['C'])/(Jy*U0*param['cluz'])
	taup=(2*param['En']*param['C'])/(Jp*U0*param['cluz'])
	#print taux,tauy,taup

	#Define step for iteration
	tt=taux/5

	# Synchrotron tune
	Qs0=sqrt(param['ap']*param['hh']*sqrt(param['Vrf']**2-U0**2)/(2*pi*param['En']))

	# Define the interpolation function for Microwave Instability
	microwave=interp2d(sigS,Curr,GT,kind='linear')

	#RF frequency
	w_rf =2*pi*(param['hh']*param['cluz']/param['C']-param['Detune0'])	#Generator Frequency

	#Creates arrays for 3HC calculation
	posz=numpy.zeros(5000)
	perfil=numpy.zeros(5000)
	pot=numpy.zeros(5000)

	#Define longitudinal scale array
	posz=numpy.arange(0,5000.)/10-250 # in milimiters

	#Cretaes an array that's a subgroup of param
	inter={}
	inter['exi']=param['ex0']
	inter['eyi']=(param['k_dw']+param['k_beta'])*param['ex0']
	inter['spi']=param['sp0']
	inter['gamma']=param['gamma']
	inter['r0']=param['r0']
	inter['Np']=param['Np']
	inter['cluz']=param['cluz']

	pot=1/(param['En']*param['C'])*param['cluz']/w_rf*(Vmain*1e3*(cos(Phi_sync_nat-phimain)-numpy.cos(posz/1000*w_rf/param['cluz']+Phi_sync_nat-phimain))+Vharm*1e3/param['mharm']*(cos(param['mharm']*pi-phiharm)-numpy.cos(param['mharm']*posz/1000*w_rf/param['cluz']+param['mharm']*pi-phiharm)))-1/(param['En']*param['C'])*U0*posz/1000
	perfil=numpy.exp(-pot/(param['ap']*param['sp0']**2))
	(pos0,sigma_mm)=calc_sigma(posz,perfil)
	inter['ssi']=sigma_mm/1000

	while (difftot>10**(-7)):
		#Add the Microwave growth rate to the longitudinal plane
		DTp=microwave(sss,param['Np'])
		#print DTp

		(Tx,Ty,Tp)=Calc_Growth(twiss,inter)
		Tx=float(Tx)/param['C']
		Ty=float(Ty)/param['C']
		Tp=float(Tp)/param['C']+DTp

		exx=(-param['ex0']+exp(2*tt*(Tx-1/taux))*(param['ex0']+inter['exi']*(-1+Tx*taux)))/(-1+Tx*taux)

		eyy=(-(param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+exp(2*tt*(Ty-1/tauy))*((param['k_dw']*param['ex0']+param['k_beta']*exx*(1-tauy/Ty))+inter['eyi']*(-1+Ty*tauy)))/(-1+Ty*tauy)
		spp=(-param['sp0']+exp(tt*(Tp-1/taup))*(param['sp0']+inter['spi']*(-1+Tp*taup)))/(-1+Tp*taup)

		#Calculate bunch length according to the RF potential (Main RF + 3HC)
		pot=1/(param['En']*param['C'])*param['cluz']/w_rf*(Vmain*1e3*(cos(Phi_sync_nat-phimain)-numpy.cos(posz/1000*w_rf/param['cluz']+Phi_sync_nat-phimain))+Vharm*1e3/param['mharm']*(cos(param['mharm']*pi-phiharm)-numpy.cos(param['mharm']*posz/1000*w_rf/param['cluz']+param['mharm']*pi-phiharm)))-1/(param['En']*param['C'])*U0*posz/1000
		perfil=numpy.exp(-pot/(param['ap']*spp**2))
		(pos0,sigma_mm)=calc_sigma(posz,perfil)
		sss=sigma_mm/1000

		#print exx,eyy,spp,sss

		diff1=abs(exx-inter['exi'])/inter['exi']
		diff2=abs(eyy-inter['eyi'])/inter['eyi']
		diff3=abs(spp-inter['spi'])/inter['spi']
		diff4=abs(sss-inter['ssi'])/inter['ssi']
		difftot=diff1+diff2+diff3+diff4
		#print difftot

		inter['exi']=exx;

		inter['eyi']=eyy;

		inter['spi']=spp;

		inter['ssi']=sss;

		time=i*tt;
		i=i+1

	return (exx,eyy,spp,sss)

max_stars_repo_name: carneirofc/lnls

max_issues_repo_name: carneirofc/lnls

max_forks_repo_name: carneirofc/lnls


text:
import cv2
import numpy as np
import os
import glob

# TODO: 1. check existence of the dataset file
#    or 2. download dataset automatically
def resize2square(img, size, inter_img=cv2.INTER_CUBIC,
                  inter_mask=cv2.INTER_LINEAR):
    """
    Detect if a given argument is the source image or its mask by its channel dims
    and resize it while preserving its aspect ratio
    """
    h, w = img.shape[:2]
    c = None if len(img.shape) < 3 else img.shape[2]
    if h == w:
        if c is None:
            return cv2.resize(img, (size, size), inter_mask)
        else:
            return cv2.resize(img, (size, size), inter_img)
    if h > w:
        dif = h
    else:
        dif = w
    x_pos = int((dif - w)/2.)
    y_pos = int((dif - h)/2.)
    if c is None:
        mask = np.zeros((dif, dif), dtype=img.dtype)
        mask[y_pos:y_pos+h, x_pos:x_pos+w] = img[:h, :w]
        return cv2.resize(mask, (size, size), inter_mask)
    else:
        mask = np.zeros((dif, dif, c), dtype=img.dtype)
        mask[y_pos:y_pos+h, x_pos:x_pos+w, :] = img[:h, :w, :]
        return cv2.resize(mask, (size, size), inter_img)

def check_dir(dir_path):
    """
    Check if a directory exists and create it if it doesn't exist
    """
    check_folder = os.path.isdir(dir_path)

    # If folder doesn't exist, then create it.
    if not check_folder:
        os.makedirs(dir_path)
        print(f"created folder:  '{dir_path}'")
    else:
        print(f"'{dir_path}' folder already exists.")

def rgb_hist_equalize(img):
    """
    Perform histogram equalization on a RGB image
    """
    ycrcb = cv2.cvtColor(img, cv2.COLOR_RGB2YCR_CB)
    channels = cv2.split(ycrcb)

    cv2.equalizeHist(channels[0], channels[0])

    cv2.merge(channels, ycrcb)
    return cv2.cvtColor(ycrcb, cv2.COLOR_YCR_CB2RGB)

def create_dataset_folder(input_path, output_path, image_size=224):
    """
    Create the smallleaves dataset by extracting and saving the image portions that contain
    large contours.
    """
    check_dir(output_path)

    # get masks and their corresponding img paths
    mask_paths = glob.glob(input_path + '/*seg.png')
    img_paths = list(map(lambda st: st.replace("seg", "img"), mask_paths))

    img_count = {
        'seg': 0,
        'og': 0,
        'error': 0,
        'black': 0
    }

    for i, (mask_path, img_path) in enumerate(zip(mask_paths, img_paths)):
        # read image and convert to RGB
        og_img = cv2.imread(img_path)
        og_img = cv2.cvtColor(og_img, cv2.COLOR_BGR2RGB)

        og_img = rgb_hist_equalize(og_img)
        og_img = cv2.cvtColor(og_img, cv2.COLOR_BGR2GRAY)

        # read mask and convert to grayscale
        mask = cv2.imread(mask_path)
        gray = cv2.cvtColor(mask, cv2.COLOR_BGR2GRAY)

        # threshold and otsu
        thresh = cv2.threshold(
            gray, 1, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1]

        # extract four contours with largest areas
        cnts = cv2.findContours(
            thresh.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
        cnts = cnts[0] if len(cnts) == 2 else cnts[1]
        cnts = sorted(cnts, key=cv2.contourArea, reverse=True)[:5]

        border = 10

        # loop over the contours
        for j, c in enumerate(cnts):
            # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
            (x, y, w, h) = cv2.boundingRect(c)

            # checking for errors because some contours are too small
            try:
                roi = thresh[y - border:y + h + border,
                             x - border:x + w + border].copy()
                og_roi = og_img[y - border:y + h + border,
                                x - border:x + w + border].copy()

                # resize to square with padded black pixels
                if roi.shape == (0, 0):
                    raise ValueError

                # prevent saving all black images
                if np.sum(roi == 255) < 10:
                    img_count['black'] += 1
                    continue

                roi = resize2square(roi, image_size)
                og_roi = resize2square(og_roi, image_size)

                og_filename = "smleaf" + \
                    str(i+1).zfill(4) + "_" + str(j+1).zfill(2) + "_img.jpg"
                mask_filename = "smleaf" + \
                    str(i+1).zfill(4) + "_" + str(j+1).zfill(2) + "_seg.jpg"

                cv2.imwrite(os.path.join(output_path, og_filename), og_roi)
                img_count['og'] += 1
                cv2.imwrite(os.path.join(output_path, mask_filename), roi)
                img_count['seg'] += 1

                # ax1.imshow(og_roi)
                # ax2.imshow(roi, cmap='gray')
            except ValueError:
                img_count['error'] += 1
                print("[ValueError] img probably too small")
                continue
        print(f"saved {img_count['seg']} images")
    print(img_count)

if __name__ == "__main__":
    img_size = 224

    # generate train folder
    folder_path = './DenseLeaves/train'
    save_path = '../dataset/train'
    create_dataset_folder(folder_path, save_path, img_size)

    # generate train folder
    folder_path = './DenseLeaves/test'
    save_path = '../dataset/test'
    create_dataset_folder(folder_path, save_path, img_size)

max_stars_repo_name: gvil-research/leaf-segmentation-unet

max_issues_repo_name: gvil-research/leaf-segmentation-unet

max_forks_repo_name: gvil-research/leaf-segmentation-unet


text:
from PyQt5.QtCore import QObject, Qt
from PyQt5.QtGui import QPainterPath, QFont, QPen
import numpy as np

from urh import constants
from urh.SceneManager import SceneManager
from urh.signalprocessing.Signal import Signal
from urh.ui.ZoomableScene import ZoomableScene
from urh.cythonext import path_creator, util, signalFunctions

class SignalSceneManager(SceneManager):
    def __init__(self, signal: Signal, parent):
        super().__init__(parent)
        self.signal = signal
        self.text_item = self.scene.addText("Loading...", QFont("Helvetica"))
        self.scene_type = 0  # 0 = Analog Signal, 1 = QuadDemodView

    def show_scene_section(self, x1: float, x2: float, subpath_ranges=None, colors=None):
        self.plot_data = self.signal.real_plot_data if self.scene_type == 0 else self.signal.qad
        super().show_scene_section(x1, x2, subpath_ranges=subpath_ranges, colors=colors)

    def init_scene(self):
        stored_minimum, stored_maximum = self.minimum, self.maximum

        if self.scene_type == 0:
            # Ensure Real plot have same y Axis
            self.plot_data = self.signal.real_plot_data
        else:
            noise_val = signalFunctions.get_noise_for_mod_type(self.scene_type - 1)
            # Bypass Min/Max calculation
            if noise_val == 0:
                # ASK
                self.minimum, self.maximum = 0, self.padding * np.max(self.signal.qad)
            else:
                self.minimum, self.maximum = 0, self.padding * noise_val

            self.plot_data = self.signal.qad

        super().init_scene(apply_padding=self.scene_type == 0)
        self.minimum, self.maximum = stored_minimum, stored_maximum

        self.line_item.setLine(0, 0, 0, 0)  # Hide Axis

        if self.scene_type == 0:
            self.scene.draw_noise_area(self.signal.noise_min_plot, self.signal.noise_max_plot - self.signal.noise_min_plot)
        else:
            self.scene.draw_sep_area(-self.signal.qad_center)

max_stars_repo_name: dspmandavid/urh

max_issues_repo_name: dspmandavid/urh

max_forks_repo_name: dspmandavid/urh


text:
import numpy as np
import matplotlib.pyplot as plt
csv = """1896,4.47083333333333
1900,4.46472925981123
1904,5.22208333333333
1908,4.1546786744085
1912,3.90331674958541
1920,3.5695126705653
1924,3.8245447722874
1928,3.62483706600308
1932,3.59284275388079
1936,3.53880791562981
1948,3.6701030927835
1952,3.39029110874116
1956,3.43642611683849
1960,3.2058300746534
1964,3.13275664573212
1968,3.32819844373346
1972,3.13583757949204
1976,3.07895880238575
1980,3.10581822490816
1984,3.06552909112454
1988,3.09357348817
1992,3.16111703598373
1996,3.14255243512264
2000,3.08527866650867
2004,3.1026582928467
2008,2.99877552632618
2012,3.03392977050993"""

if sys.version_info[0] >= 3:
    import io # Python3
    olympics = np.genfromtxt(io.BytesIO(csv.encode()), delimiter=",")
else:
    from StringIO import StringIO  # Python2
    olympics = np.genfromtxt(StringIO(csv), delimiter=',') #Python 2

#print(olympics)
x = olympics[:, 0:1]
y = olympics[:, 1:2]
plt.plot(x, y, 'rx')
plt.show()

num_data = x.shape[0]
num_pred_data = 100 # how many points to use for plotting predictions
x_pred = linspace(1890, 2016, num_pred_data)[:, None] # input locations for predictions
order = 4 # The polynomial order to use.
print ('Num of training samples: ',num_data)
print('Num of testing samples: ',num_pred_data)

Phi = np.zeros((num_data, order+1))
Phi_pred = np.zeros((num_pred_data, order+1))
for i in range(0, order+1):
    Phi[:, i:i+1] = x**i
    Phi_pred[:, i:i+1] = x_pred**i

w = np.linalg.solve(np.dot(Phi.T, Phi), np.dot(Phi.T, y)) # back to you

f = Phi*w
f_pred = Phi_pred*w

SSR = sum((y-f)**2)
print(SSR)

#

max_stars_repo_name: inverthermit/machine-learning-playboard

max_issues_repo_name: inverthermit/machine-learning-playboard

max_forks_repo_name: inverthermit/machine-learning-playboard


text:
import numpy as np

from autumn.tools.project import Project, ParameterSet, TimeSeriesSet, build_rel_path, get_all_available_scenario_paths, use_tuned_proposal_sds
from autumn.tools.calibration import Calibration
from autumn.tools.project.params import read_yaml_file
from autumn.tools.calibration.priors import UniformPrior, TruncNormalPrior
from autumn.tools.calibration.targets import NormalTarget, PoissonTarget, TruncNormalTarget
from autumn.models.covid_19 import base_params, build_model
from autumn.settings import Region, Models
import os

CLUSTERS = [Region.to_filename(r) for r in Region.VICTORIA_SUBREGIONS]

# Just calibrate to June, July, August and September for now (but run for some lead in time at the start)
TARGETS_START_TIME = 153  # 1st June
TARGETS_END_TIME = 305  # 31st October
TARGETS_RANGE = (TARGETS_START_TIME, TARGETS_END_TIME)

# Load and configure model parameters
default_path = build_rel_path("params/default.yml")
mle_path = build_rel_path("params/mle-params.yml")
scenario_dir_path = build_rel_path("params/")
scenario_paths = get_all_available_scenario_paths(scenario_dir_path)
baseline_params = base_params.update(default_path).update(mle_path, calibration_format=True)
scenario_params = [baseline_params.update(p) for p in scenario_paths]
param_set = ParameterSet(baseline=baseline_params, scenarios=scenario_params)

# Add calibration targets and priors
ts_set = TimeSeriesSet.from_file(build_rel_path("targets.secret.json"))

# For all the cluster targets, a universal calibrated parameter called "target_output_ratio" is used to scale the
# dispersion parameter of the targets' normal likelihoods.
cluster_targets = []
for cluster in CLUSTERS:
    notifs_ts = ts_set.get(f"notificationsXcluster_{cluster}").moving_average(4)
    target = NormalTarget(notifs_ts)
    cluster_targets.append(target)

# Request calibration targets
targets = [
    PoissonTarget(ts_set.get("notifications").truncate_times(*TARGETS_RANGE).round_values()),
    PoissonTarget(ts_set.get("infection_deaths").truncate_times(*TARGETS_RANGE).moving_average(7)),
    PoissonTarget(ts_set.get("hospital_admissions").truncate_times(*TARGETS_RANGE)),
    PoissonTarget(ts_set.get("icu_admissions").truncate_times(*TARGETS_RANGE)),
    *cluster_targets,
]

# Add multiplier for most services, except use South Metro for South East Metro, use North Metro for West Metro
cluster_priors = []
regions_for_multipliers = Region.VICTORIA_METRO
regions_for_multipliers.append(Region.BARWON_SOUTH_WEST)

for region in regions_for_multipliers:
    region_name = region.replace("-", "_")
    name = f"victorian_clusters.contact_rate_multiplier_{region_name}"
    # Shouldn't be too peaked with these values.
    prior = TruncNormalPrior(name, mean=1.0, stdev=0.5, trunc_range=[0.5, np.inf], jumping_stdev=0.15)
    cluster_priors.append(prior)

priors = [
    # Global COVID priors, but with jumping sds adjusted
    TruncNormalPrior(
        "sojourn.compartment_periods_calculated.exposed.total_period",
        mean=5.5,
        stdev=0.97,
        trunc_range=[1.0, np.inf],
        jumping_stdev=0.5,
    ),
    TruncNormalPrior(
        "sojourn.compartment_periods_calculated.active.total_period",
        mean=6.5,
        stdev=0.77,
        trunc_range=[4.0, np.inf],
        jumping_stdev=0.4,
    ),
    # Cluster specific priors.
    *cluster_priors,
    # Victorian regional priors.
    TruncNormalPrior(
        f"victorian_clusters.contact_rate_multiplier_regional",
        mean=1.0,
        stdev=0.5,
        trunc_range=[0.5, np.inf],
        jumping_stdev=0.15,
    ),
    UniformPrior("contact_rate", [0.04, 0.07], jumping_stdev=0.008),
    UniformPrior("victorian_clusters.intercluster_mixing", [0.005, 0.05], jumping_stdev=0.01),
    UniformPrior("infectious_seed", [20., 70.], jumping_stdev=4.),
    UniformPrior("clinical_stratification.non_sympt_infect_multiplier", [0.2, 0.8], jumping_stdev=0.05),
    UniformPrior("infection_fatality.top_bracket_overwrite", [0.05, 0.3], jumping_stdev=0.04),
    UniformPrior("clinical_stratification.props.hospital.multiplier", [0.5, 5.0], jumping_stdev=0.4),
    UniformPrior("testing_to_detection.assumed_cdr_parameter", [0.05, 0.3], jumping_stdev=0.04),
    TruncNormalPrior(
        "sojourn.compartment_periods.icu_early",
        mean=12.7,
        stdev=4.0,
        trunc_range=[5.0, np.inf],
        jumping_stdev=4.
    ),
    UniformPrior(
        "victorian_clusters.metro.mobility.microdistancing.behaviour_adjuster.parameters.effect",
        [0.0, 0.6],
        jumping_stdev=0.075,
    ),
    UniformPrior(
        "victorian_clusters.metro.mobility.microdistancing.face_coverings_adjuster.parameters.effect",
        [0.0, 0.6],
        jumping_stdev=0.04,
    ),
    UniformPrior(
        "victorian_clusters.metro.mobility.microdistancing.home_reduction.parameters.effect",
        [0.0, 0.4],
        jumping_stdev=0.04,
    ),
    UniformPrior("target_output_ratio", [0.2, 0.7], jumping_stdev=0.04),
    UniformPrior("contact_tracing.assumed_trace_prop", [0.2, 0.5], jumping_stdev=0.04),
]

# Load proposal sds from yml file
use_tuned_proposal_sds(priors, build_rel_path("proposal_sds.yml"))

calibration = Calibration(
    priors,
    targets,
    metropolis_init="current_params",
    metropolis_init_rel_step_size=0.05,
    fixed_proposal_steps=500,
    jumping_stdev_adjustment=0.8,
)

# FIXME: Replace with flexible Python plot request API.
import json

plot_spec_filepath = build_rel_path("targets.secret.json")
with open(plot_spec_filepath) as f:
    plot_spec = json.load(f)

project = Project(
    Region.VICTORIA_2020, Models.COVID_19, build_model, param_set, calibration, plots=plot_spec
)

# Write parameter table to tex file
main_table_params_list = [
    "clinical_stratification.icu_prop",
    "sojourn.compartment_periods_calculated.exposed.total_period",
    "contact_rate"
]
# project.write_params_to_tex(main_table_params_list, project_path=build_rel_path(''))

# from autumn.tools.calibration.proposal_tuning import perform_all_params_proposal_tuning
# perform_all_params_proposal_tuning(project, calibration, priors, n_points=100, relative_likelihood_reduction=0.2)

max_stars_repo_name: monash-emu/AuTuMN

max_issues_repo_name: monash-emu/AuTuMN

max_forks_repo_name: monash-emu/AuTuMN


text:
"""
Unit tests for the spline interpolator component.
"""
from copy import deepcopy
import unittest

import numpy as np

from openmdao.components.interp_util.interp import InterpND, SPLINE_METHODS, TABLE_METHODS
from openmdao.components.interp_util.interp_semi import InterpNDSemi, INTERP_METHODS
from openmdao.components.interp_util.outofbounds_error import OutOfBoundsError
from openmdao.utils.assert_utils import assert_near_equal, assert_equal_arrays

def rel_error(actual, computed):
    return np.linalg.norm(actual - computed) / np.linalg.norm(actual)

scipy_gte_019 = True
try:
    from scipy.interpolate._bsplines import make_interp_spline
except ImportError:
    scipy_gte_019 = False

class InterpNDStandaloneFeatureTestcase(unittest.TestCase):

    def test_interp_spline_akima(self):

        xcp = np.array([1.0, 2.0, 4.0, 6.0, 10.0, 12.0])
        ycp = np.array([5.0, 12.0, 14.0, 16.0, 21.0, 29.0])
        n = 50
        x = np.linspace(1.0, 12.0, n)

        interp = InterpND(method='akima', points=xcp, x_interp=x, delta_x=0.1)
        y = interp.evaluate_spline(ycp)

        assert_near_equal(y,
                          np.array([ 5.       ,  7.20902005,  9.21276849, 10.81097162, 11.80335574,
                                     12.1278001 , 12.35869145, 12.58588536, 12.81022332, 13.03254681,
                                     13.25369732, 13.47451633, 13.69584534, 13.91852582, 14.14281484,
                                     14.36710105, 14.59128625, 14.81544619, 15.03965664, 15.26399335,
                                     15.48853209, 15.7133486 , 15.93851866, 16.16573502, 16.39927111,
                                     16.63928669, 16.8857123 , 17.1384785 , 17.39751585, 17.66275489,
                                     17.93412619, 18.21156029, 18.49498776, 18.78433915, 19.07954501,
                                     19.38053589, 19.68724235, 19.99959495, 20.31752423, 20.64096076,
                                     20.96983509, 21.37579297, 21.94811407, 22.66809748, 23.51629844,
                                     24.47327219, 25.51957398, 26.63575905, 27.80238264, 29.        ]),
                          tolerance=1e-6)

    def test_interp_spline_akima_derivs(self):

        xcp = np.array([1.0, 2.0, 4.0, 6.0, 10.0, 12.0])
        ycp = np.array([5.0, 12.0, 14.0, 16.0, 21.0, 29.0])
        n = 5
        x = np.linspace(1.0, 12.0, n)

        interp = InterpND(method='akima', points=xcp, x_interp=x, delta_x=0.1)
        y, dy_dycp = interp.evaluate_spline(ycp, compute_derivative=True)

        assert_near_equal(dy_dycp,
                          np.array([[ 1.00000000e+00,  0.00000000e+00,  0.00000000e+00,
                                      0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
                                    [-1.86761492e-06,  3.31278014e-02,  1.05874907e+00,
                                     -9.18750000e-02,  0.00000000e+00,  0.00000000e+00],
                                    [ 0.00000000e+00,  0.00000000e+00, -2.10964627e-01,
                                      1.19119941e+00,  2.02602810e-02, -4.95062934e-04],
                                    [ 0.00000000e+00,  0.00000000e+00, -2.64126732e-01,
                                      5.82784977e-01,  6.83151998e-01, -1.81024253e-03],
                                    [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
                                      0.00000000e+00,  0.00000000e+00,  1.00000000e+00]]),
                          tolerance=1e-6)

    def test_interp_spline_bsplines(self):

        xcp = np.array([1.0, 2.0, 4.0, 6.0, 10.0, 12.0])
        ycp = np.array([5.0, 12.0, 14.0, 16.0, 21.0, 29.0])
        n = 50
        x = np.linspace(1.0, 12.0, n)

        interp = InterpND(method='bsplines', num_cp=6, x_interp=x)

        y = interp.evaluate_spline(ycp)

        assert_near_equal(y,
                          np.array([ 9.21614583,  9.90911525, 10.52244151, 11.06231159, 11.53491244,
                                     11.94643105, 12.30305438, 12.61096939, 12.87636305, 13.10542234,
                                     13.30433422, 13.47928566, 13.63646363, 13.7820551 , 13.92203064,
                                     14.05954727, 14.19579437, 14.33192094, 14.46907599, 14.60840854,
                                     14.75106758, 14.89820214, 15.05096121, 15.2104938 , 15.37794893,
                                     15.5544756 , 15.74122282, 15.9393396 , 16.14997495, 16.37427787,
                                     16.61339737, 16.86848247, 17.14102103, 17.43486416, 17.75486932,
                                     18.10589772, 18.49281052, 18.92046894, 19.39373414, 19.91746734,
                                     20.4965297 , 21.13578243, 21.8400867 , 22.61430372, 23.46329467,
                                     24.39192074, 25.40504312, 26.507523  , 27.70422156, 29.        ]),
                          tolerance=1e-6)

    def test_table_interp(self):

        # create input param training data, of sizes 25, 5, and 10 points resp.
        p1 = np.linspace(0, 100, 25)
        p2 = np.linspace(-10, 10, 5)
        p3 = np.linspace(0, 1, 10)

        # can use meshgrid to create a 3D array of test data
        P1, P2, P3 = np.meshgrid(p1, p2, p3, indexing='ij')
        f = np.sqrt(P1) + P2 * P3

        interp = InterpND(method='lagrange3', points=(p1, p2, p3), values=f)

        x = np.array([55.12, -2.14, 0.323])
        f, df_dx = interp.interpolate(x, compute_derivative=True)

        actual = np.array([6.73306794])
        deriv_actual = np.array([[ 0.06734927, 0.323 , -2.14]])

        assert_near_equal(f, actual, tolerance=1e-7)
        assert_near_equal(df_dx, deriv_actual, tolerance=1e-7)

    def test_cs_across_interp(self):
        # The standalone interpolator is used inside of components, so the imaginary part must
        # be carried through all operations to the outputs.
        xcp = np.array([1.0, 2.0, 4.0, 6.0, 10.0, 12.0])
        ycp = np.array([5.0, 12.0, 14.0, 16.0, 21.0, 29.0])
        n = 50
        x = np.linspace(1.0, 12.0, n)

        ycp = np.array([[5.0 + 1j, 12.0, 14.0, 16.0, 21.0, 29.0],
                        [5.0, 12.0 + 1j, 14.0, 16.0, 21.0, 29.0]])

        for method in SPLINE_METHODS:

            # complex step not supported on scipy methods
            if method.startswith('scipy'):
                continue

            interp = InterpND(method=method, points=xcp, x_interp=x)
            y, dy = interp.evaluate_spline(ycp, compute_derivative=True)

            self.assertTrue(y.dtype == complex)

            if method in ['akima']:
                # Derivs depend on values only for akima.
                self.assertTrue(dy.dtype == complex)

        p1 = np.linspace(0, 100, 25)
        p2 = np.linspace(-10, 10, 5)
        p3 = np.linspace(0, 1, 10)

        # can use meshgrid to create a 3D array of test data
        P1, P2, P3 = np.meshgrid(p1, p2, p3, indexing='ij')
        f = np.sqrt(P1) + P2 * P3

        x = np.array([[55.12 + 1j, -2.14, 0.323],
                      [55.12, -2.14 + 1j, 0.323],
                      [55.12, -2.14, 0.323 + 1j]])

        for method in TABLE_METHODS:

            # complex step not supported on scipy methods
            if method.startswith('scipy'):
                continue

            interp = InterpND(method=method, points=(p1, p2, p3), values=f)
            y, dy = interp.interpolate(x, compute_derivative=True)

            self.assertTrue(y.dtype == complex)
            self.assertTrue(dy.dtype == complex)

class TestInterpNDSemiPython(unittest.TestCase):
    """Tests for the standalone semi structured interp."""

    def setUp(self):
        self.interp_configs = {
            "slinear": 2,
            "lagrange2": 3,
            "lagrange3": 4,
            "akima": 4,
        }
        self.interp_methods = self.interp_configs.keys()

        self.tol = {
            "slinear": 5e-2,
            "lagrange2": 5e-2,
            "lagrange3": 1e-4,
            "akima": 1e-3,
        }

    def _get_sample_2d(self):
        # test problem with enough points for smooth spline fits
        def f(u, v):
            return u * np.cos(u * v) + v * np.sin(u * v)

        def df(u, v):
            return (-u * v * np.sin(u * v) + v**2 * np.cos(u * v) +
                    np.cos(u * v),
                    -u**2 * np.sin(u * v) + u * v * np.cos(u * v) +
                    np.sin(u * v))

        # uniformly spaced axis
        u = np.linspace(0, 3, 50)
        # randomly spaced axis
        np.random.seed(7590)
        v = np.random.uniform(0, 3, 50)
        v.sort()

        points = [u, v]
        values = f(*np.meshgrid(*points, indexing='ij'))
        return points, values, f, df

    def _get_sample_4d_large(self):
        def f(x, y, z, w):
            return x**2 + y**2 + z**2 + w**2
        X = np.linspace(-10, 10, 6)
        Y = np.linspace(-10, 10, 7)
        np.random.seed(0)
        Z = np.random.uniform(-10, 10, 6)
        Z.sort()
        W = np.linspace(-10, 10, 8)
        points = [X, Y, Z, W]
        values = f(*np.meshgrid(*points, indexing='ij'))
        return points, values

    def test_2d(self):
        # test interpolated values
        points, values, func, df = self._get_sample_2d()
        np.random.seed(1)
        X, Y = np.meshgrid(*points, indexing='ij')
        X = X.ravel()
        Y = Y.ravel()
        grid = np.array([X, Y]).T

        test_pt = np.random.uniform(0, 3, 2)
        actual = func(*test_pt)
        for method in self.interp_methods:
            interp = InterpNDSemi(grid, values.ravel(), method=method)
            computed = interp.interpolate(test_pt)
            r_err = rel_error(actual, computed)
            assert r_err < self.tol[method]

    def test_minimum_required_gridsize(self):
        for method in self.interp_methods:

            k = self.interp_configs[method] - 1
            x = np.linspace(0, 1, k)
            y = np.linspace(0, 1, k)

            points = [x, y]
            X, Y = np.meshgrid(*points, indexing='ij')
            X = X.ravel()
            Y = Y.ravel()
            values = X + Y
            grid = np.array([X, Y]).T
            with self.assertRaises(ValueError) as cm:
                interp = InterpNDSemi(grid, values, method=method)

            msg = 'There are {} points in a data dimension, but method'.format(k)
            self.assertTrue(str(cm.exception).startswith(msg))

    def test_NaN_exception(self):
        np.random.seed(1234)
        x = np.linspace(0, 2, 5)
        y = np.linspace(0, 1, 7)
        points = [x, y]
        X, Y = np.meshgrid(*points, indexing='ij')
        X = X.ravel()
        Y = Y.ravel()
        grid = np.array([X, Y]).T
        values = np.random.rand(5, 7).ravel()
        interp = InterpNDSemi(grid, values, method='slinear', extrapolate=False)

        with self.assertRaises(OutOfBoundsError) as cm:
            interp.interpolate(np.array([1, np.nan]))

        err = cm.exception

        self.assertEqual(str(err), 'One of the requested xi contains a NaN')
        self.assertEqual(err.idx, 1)
        self.assertTrue(np.isnan(err.value))

    def test_error_messages(self):
        points, values = self._get_sample_4d_large()
        X, Y, Z, W = np.meshgrid(*points, indexing='ij')
        X = X.ravel()
        Y = Y.ravel()
        Z = Z.ravel()
        W = W.ravel()
        grid = np.array([X, Y, Z, W]).T
        values = values.ravel()

        with self.assertRaises(ValueError) as cm:
            interp = InterpNDSemi(grid, values, method='junk')

        msg = ('Interpolation method "junk" is not defined. Valid methods are')
        self.assertTrue(cm.exception.args[0].startswith(msg))

        with self.assertRaises(ValueError) as cm:
            interp = InterpNDSemi(grid, values, method=points)

        msg = ("Argument 'method' should be a string.")
        self.assertTrue(cm.exception.args[0].startswith(msg))

        with self.assertRaises(ValueError) as cm:
            interp = InterpNDSemi(grid, values[:-1], method='slinear')

        msg = ('There are 2016 point arrays, but 2015 values.')
        self.assertEqual(cm.exception.args[0], msg)

        badgrid = deepcopy(grid)
        badgrid[0][0] = -6.0
        with self.assertRaises(ValueError) as cm:
            interp = InterpNDSemi(badgrid, values, method='slinear')

        msg = ('The points in dimension 0 must be strictly ascending.')
        self.assertEqual(cm.exception.args[0], msg)

        with self.assertRaises(KeyError) as cm:
            interp = InterpNDSemi(grid, values, method='slinear', bad_arg=1)

        msg = ("Option 'bad_arg' cannot be set because it has not been declared.")
        self.assertTrue(cm.exception.args[0].endswith(msg))

class TestInterpNDPython(unittest.TestCase):
    """Tests for the non-scipy interpolation algorithms."""

    def setUp(self):
        self.interp_configs = {
            "slinear": 2,
            "cubic": 3,
            "lagrange2": 3,
            "lagrange3": 4,
            "akima": 4,
            "scipy_slinear": 1,
            "scipy_cubic": 3,
            "scipy_quintic": 5,
        }
        self.spline_configs = {
            "slinear": 2,
            "cubic": 3,
            "lagrange2": 3,
            "lagrange3": 4,
            "akima": 4,
            "bsplines": 4,
            "scipy_slinear": 1,
            "scipy_cubic": 3,
            "scipy_quintic": 5,
        }
        self.interp_methods = self.interp_configs.keys()
        self.spline_methods = self.spline_configs.keys()

        self.tol = {
            "slinear": 5e-2,
            "lagrange2": 5e-2,
            "lagrange3": 1e-4,
            "cubic": 1e-4,
            "akima": 1e-3,
            "bsplines": 1e-1,
            "scipy_slinear": 5e-2,
            "scipy_cubic": 1e-4,
            "scipy_quintic": 1e-6,
        }

    def _get_sample_1d(self):
        # test problem with enough points for smooth spline fits
        def f(u):
            return 2.0 * np.cos(u)

        def df(u):
            return -2.0 * np.sin(u)

        # uniformly spaced axis
        u = np.linspace(0, 3, 50)

        points = [u]
        values = f(u)
        return points, values, f, df

    def _get_sample_2d(self):
        # test problem with enough points for smooth spline fits
        def f(u, v):
            return u * np.cos(u * v) + v * np.sin(u * v)

        def df(u, v):
            return (-u * v * np.sin(u * v) + v**2 * np.cos(u * v) +
                    np.cos(u * v),
                    -u**2 * np.sin(u * v) + u * v * np.cos(u * v) +
                    np.sin(u * v))

        # uniformly spaced axis
        u = np.linspace(0, 3, 50)
        # randomly spaced axis
        np.random.seed(7590)
        v = np.random.uniform(0, 3, 50)
        v.sort()

        points = [u, v]
        values = f(*np.meshgrid(*points, indexing='ij'))
        return points, values, f, df

    def _get_sample_4d_large(self):
        def f(x, y, z, w):
            return x**2 + y**2 + z**2 + w**2
        X = np.linspace(-10, 10, 6)
        Y = np.linspace(-10, 10, 7)
        np.random.seed(0)
        Z = np.random.uniform(-10, 10, 6)
        Z.sort()
        W = np.linspace(-10, 10, 8)
        points = [X, Y, Z, W]
        values = f(*np.meshgrid(*points, indexing='ij'))
        return points, values

    def test_minimum_required_gridsize(self):
        for method in self.interp_methods:

            # Scipy does order reduction as needed.
            if method.startswith('scipy'):
                continue

            k = self.interp_configs[method] - 1
            x = np.linspace(0, 1, k)
            y = np.linspace(0, 1, k)
            points = [x, y]
            X, Y = np.meshgrid(*points, indexing='ij')
            values = X + Y
            #self.assertRaises(ValueError, InterpND, points, values, method)
            with self.assertRaises(ValueError) as cm:
                interp = InterpND(method=method, points=points, values=values)

            msg = 'There are {} points in a data dimension, but method'.format(k)
            self.assertTrue(str(cm.exception).startswith(msg))

    def test_spline_single_dim(self):
        # test interpolated values
        points, values, func, df = self._get_sample_1d()
        test_pt = np.array([[0.76], [.33]])
        actual = func(test_pt).flatten()
        for method in self.interp_methods:
            interp = InterpND(method=method, points=points, values=values)
            computed = interp.interpolate(test_pt)
            r_err = rel_error(actual, computed)
            assert r_err < self.tol[method]

    def test_spline_xi1d(self):
        # test interpolated values
        points, values, func, df = self._get_sample_2d()
        np.random.seed(1)
        test_pt = np.random.uniform(0, 3, 2)
        actual = func(*test_pt)
        for method in self.interp_methods:
            interp = InterpND(method=method, points=points, values=values)
            computed = interp.interpolate(test_pt)
            r_err = rel_error(actual, computed)
            assert r_err < self.tol[method]

    def test_spline_out_of_bounds_extrap(self):
        points, values, func, df = self. _get_sample_2d()
        np.random.seed(5)
        test_pt = np.random.uniform(3, 3.1, 2)
        actual = func(*test_pt)
        gradient = np.array(df(*test_pt))
        tol = 1e-1
        for method in self.interp_methods:
            k = self.interp_configs[method]
            if method == 'slinear':
                tol = 2
            interp = InterpND(method=method, points=points, values=values,
                              extrapolate=True)
            computed, computed_grad = interp.interpolate(test_pt, compute_derivative=True)
            computed_grad = interp.gradient(test_pt)
            r_err = rel_error(actual, computed)
            assert r_err < tol

            r_err = rel_error(gradient, computed_grad)
            # extrapolated gradients are even trickier, but usable still
            assert r_err < 2 * tol

    def test_spline_xi3d(self):
        points, values, func, df = self. _get_sample_2d()
        np.random.seed(1)
        test_pt = np.random.uniform(0, 3, 6).reshape(3, 2)
        actual = func(*test_pt.T)
        for method in self.interp_methods:
            interp = InterpND(method=method, points=points, values=values)
            computed = interp.interpolate(test_pt)
            r_err = rel_error(actual, computed)
            #print(method, computed, actual, r_err)
            assert r_err < self.tol[method]

    def test_spline_xi3d_akima_delta_x(self):
        points, values, func, df = self. _get_sample_2d()
        np.random.seed(1)
        test_pt = np.random.uniform(0, 3, 6).reshape(3, 2)
        actual = func(*test_pt.T)
        interp = InterpND(method='akima', points=points, values=values, delta_x=0.01)
        computed = interp.interpolate(test_pt)
        r_err = rel_error(actual, computed)
        #print('akima', computed, actual, r_err)
        assert r_err < self.tol['akima']

    def test_spline_deriv_xi1d(self):
        # tests gradient values
        points, values, func, df = self._get_sample_2d()
        np.random.seed(1234)
        test_pt = np.random.uniform(0, 3, 2)
        actual = np.array(df(*test_pt))
        for method in self.interp_methods:
            interp = InterpND(method=method,points=points, values=values)
            computed = interp.gradient(test_pt)
            r_err = rel_error(actual, computed)
            assert r_err < 2.5 * self.tol[method]

            # test that gradients have been cached
            assert_equal_arrays(interp._xi.flatten(), test_pt.flatten())
            assert_equal_arrays(interp._d_dx.flatten(), computed.flatten())

    def test_gradients_returned_by_xi(self):
        # verifies that gradients with respect to xi are returned if cached
        points, values, func, df = self._get_sample_2d()
        np.random.seed(4321)
        for method in self.interp_methods:
            interp = InterpND(method=method, points=points, values=values)
            x = np.array([0.9, 0.1])
            interp._xi = x
            dy = np.array([0.997901, 0.08915])
            interp._d_dx = dy
            assert_near_equal(interp.gradient(x), dy, tolerance=1e-7)

    def test_akima_interpolating_spline(self):
        n_cp = 80
        n_point = 160

        t = np.linspace(0, 3.0*np.pi, n_cp)
        tt = np.linspace(0, 3.0*np.pi, n_point)

        x = np.sin(t)

        # Now, test newer interface for order_reducing spline.

        interp = InterpND(method='akima', points=t, x_interp=tt)
        computed = interp.evaluate_spline(x)

        x_expected = np.sin(tt)
        delta = computed.flatten() - x_expected

        # Here we test that we don't have crazy interpolation error.
        self.assertLess(max(delta), .15)
        # And that it gets middle points a little better.
        self.assertLess(max(delta[15:-15]), .06)

    def test_bsplines_interpolating_spline(self):
        n_cp = 80
        n_point = 160

        t = np.linspace(0, 3.0*np.pi, n_cp)
        tt = np.linspace(0, 3.0*np.pi, n_point)

        x = np.sin(t)

        # Now, test newer interface for order_reducing spline.

        interp = InterpND(method='bsplines', points=t, x_interp=tt)
        computed = interp.evaluate_spline(x)

        x_expected = np.sin(tt)
        delta = computed.flatten() - x_expected

        # Here we test that we don't have crazy interpolation error.
        self.assertLess(max(delta), .15)
        # And that it gets middle points a little better.
        self.assertLess(max(delta[15:-15]), .06)

    def test_interpolating_spline_derivs(self):
        n_cp = 12
        n_point = 21

        t = np.linspace(0, 3.0*np.pi, n_cp)
        tt = np.linspace(0, 3.0*np.pi, n_point)

        x = np.sin(t)

        for method in self.spline_methods:
            interp = InterpND(method=method, points=t, x_interp=tt)
            computed, deriv = interp.evaluate_spline(x, compute_derivative=True)

            x_expected = np.sin(tt)
            delta = computed.flatten() - x_expected

            # Here we test that we don't have crazy interpolation error.
            self.assertLess(max(delta), .25)

    def test_interp_akima_derivs(self):

        xcp = np.array([1.0, 2.0, 4.0, 6.0, 10.0, 12.0])
        ycp = np.array([5.0, 12.0, 14.0, 16.0, 21.0, 29.0])
        n = 23
        x = np.linspace(1.0, 12.0, n)

        interp = InterpND(method='akima', points=xcp, x_interp=x, delta_x=0.1)
        y, dy_dycp = interp.evaluate_spline(ycp, compute_derivative=True)

        deriv = np.array([[ 1.00000000e+00,  0.00000000e+00,  0.00000000e+00,
                            0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
                          [ 3.12508538e-01,  7.81237193e-01, -9.37457312e-02,
                            0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
                          [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00,
                            0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
                          [-1.92097534e-05,  7.05028815e-01,  3.39990395e-01,
                           -4.50000000e-02,  0.00000000e+00,  0.00000000e+00],
                          [-1.70753364e-05,  3.80025613e-01,  7.39991462e-01,
                           -1.20000000e-01,  0.00000000e+00,  0.00000000e+00],
                          [-6.40325114e-06,  1.15009605e-01,  1.01999680e+00,
                           -1.35000000e-01,  0.00000000e+00,  0.00000000e+00],
                          [ 0.00000000e+00,  1.11022302e-16,  1.00000000e+00,
                            0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
                          [ 0.00000000e+00, -5.62500000e-03,  7.60412946e-01,
                            2.45667949e-01, -5.30632175e-04,  7.47369260e-05],
                          [ 0.00000000e+00, -5.00000000e-03,  5.07767857e-01,
                            4.98447864e-01, -1.41501913e-03,  1.99298469e-04],
                          [ 0.00000000e+00, -1.87500000e-03,  2.51238839e-01,
                            7.52003846e-01, -1.59189652e-03,  2.24210778e-04],
                          [ 0.00000000e+00,  0.00000000e+00,  1.11022302e-16,
                            1.00000000e+00,  0.00000000e+00,  0.00000000e+00],
                          [ 0.00000000e+00,  0.00000000e+00, -2.10964627e-01,
                            1.19119941e+00,  2.02602810e-02, -4.95062934e-04],
                          [ 0.00000000e+00,  0.00000000e+00, -3.55003720e-01,
                            1.28195585e+00,  7.41473347e-02, -1.09946322e-03],
                          [ 0.00000000e+00,  0.00000000e+00, -4.35442243e-01,
                            1.27731946e+00,  1.59810592e-01, -1.68780891e-03],
                          [ 0.00000000e+00,  0.00000000e+00, -4.55605159e-01,
                            1.18234038e+00,  2.75399483e-01, -2.13470805e-03],
                          [ 0.00000000e+00,  0.00000000e+00, -4.18817429e-01,
                            1.00206876e+00,  4.19063438e-01, -2.31476868e-03],
                          [ 0.00000000e+00,  0.00000000e+00, -3.28404018e-01,
                            7.41554727e-01,  5.88951889e-01, -2.10259885e-03],
                          [ 0.00000000e+00,  0.00000000e+00, -1.87689887e-01,
                            4.05848427e-01,  7.83214266e-01, -1.37280661e-03],
                          [ 0.00000000e+00,  0.00000000e+00,  2.22044605e-16,
                            0.00000000e+00,  1.00000000e+00,  4.33680869e-19],
                          [ 0.00000000e+00,  0.00000000e+00,  1.18164062e-01,
                           -2.58789062e-01,  1.05371094e+00,  8.69140625e-02],
                          [ 0.00000000e+00,  0.00000000e+00,  1.05034722e-01,
                           -2.50868056e-01,  8.32465278e-01,  3.13368056e-01],
                          [ 0.00000000e+00,  0.00000000e+00,  3.93880208e-02,
                           -1.17513021e-01,  4.44986979e-01,  6.33138021e-01],
                          [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
                            0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])

        assert_near_equal(deriv, dy_dycp, tolerance=1e-6)

    def test_scipy_auto_reduce_spline_order(self):
        # if a spline method is used and spline_dim_error=False and a dimension
        # does not have enough points, the spline order for that dimension
        # should be automatically reduced
        np.random.seed(314)

        # x dimension is too small for cubic, should fall back to linear
        x = [0, 1]
        y = np.linspace(-10, 4, 10)
        z = np.linspace(1000, 2000, 20)

        points = [x, y, z]
        values = np.random.randn(2, 10, 20)

        interp = InterpND(method='scipy_cubic', points=points, values=values)

        # first dimension (x) should be reduced to k=1 (linear)
        self.assertEqual(interp.table._ki[0], 1)

        # should operate as normal
        x = np.array([0.5, 0, 1001])
        result = interp.interpolate(x)
        assert_near_equal(result, -0.046325695741704434, tolerance=1e-5)

        interp = InterpND(method='scipy_slinear', points=points, values=values)
        value1 = interp.interpolate(x)

        # cycle through different methods that require order reduction
        # in the first dimension
        interp = InterpND(method='scipy_quintic', points=points, values=values)
        value2, grad2 = interp.interpolate(x, compute_derivative=True)
        interp = InterpND(method='scipy_cubic', points=points, values=values)
        value3, grad3 = interp.interpolate(x, compute_derivative=True)

        # values from different methods should be different
        self.assertTrue(value1[0] != value2[0])
        self.assertTrue(value2[0] != value3[0])

    def test_NaN_exception(self):
        np.random.seed(1234)
        x = np.linspace(0, 2, 5)
        y = np.linspace(0, 1, 7)
        values = np.random.rand(5, 7)
        interp = InterpND(method='slinear', points=(x, y), values=values)

        with self.assertRaises(OutOfBoundsError) as cm:
            interp.interpolate(np.array([1, np.nan]))

        err = cm.exception

        self.assertEqual(str(err), 'One of the requested xi contains a NaN')
        self.assertEqual(err.idx, 1)
        self.assertTrue(np.isnan(err.value))
        self.assertEqual(err.lower, 0)
        self.assertEqual(err.upper, 1)

    def test_interp_1Dflat(self):

        x = np.array([0.,1.,2.,3.,4.])
        y = x**2
        f = InterpND(points=x, values=y)

        computed = f.interpolate(np.array([2.5, 3.5]))
        assert_equal_arrays(computed, np.array([6.5, 12.5]))

    def test_interp_1Dflat_list_points(self):

        x = np.array([0.,1.,2.,3.,4.])
        y = x**2
        f = InterpND(points=[x], values=y)

        computed = f.interpolate(np.array([2.5, 3.5]))
        assert_equal_arrays(computed, np.array([6.5, 12.5]))

    def test_interp_1Dflat_list_x(self):

        x = np.array([0.,1.,2.,3.,4.])
        y = x**2
        f = InterpND(points=x, values=y)

        computed = f.interpolate([np.array([2.5]), np.array([3.5])])
        assert_equal_arrays(computed, np.array([6.5, 12.5]))

    def test_error_messages(self):
        points, values = self._get_sample_4d_large()

        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='junk', points=points, values=values.tolist())

        msg = ('Interpolation method "junk" is not defined. Valid methods are')
        self.assertTrue(str(cm.exception).startswith(msg))

        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method=points, points=points, values=values.tolist())

        msg = ("Argument 'method' should be a string.")
        self.assertTrue(str(cm.exception).startswith(msg))

        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='slinear', points=points, values=values.tolist()[1])

        msg = ('There are 4 point arrays, but values has 3 dimensions')
        self.assertEqual(str(cm.exception), msg)

        badpoints = deepcopy(points)
        badpoints[0][0] = 55.0
        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='slinear', points=badpoints, values=values.tolist())

        msg = ('The points in dimension 0 must be strictly ascending')
        self.assertEqual(str(cm.exception), msg)

        badpoints[0] = np.vstack((np.arange(6), np.arange(6)))
        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='slinear', points=badpoints, values=values.tolist())

        msg = ('The points in dimension 0 must be 1-dimensional')
        self.assertEqual(str(cm.exception), msg)

        badpoints[0] = (np.arange(4))
        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='slinear', points=badpoints, values=values.tolist())

        msg = ('There are 4 points and 6 values in dimension 0')
        self.assertEqual(str(cm.exception), msg)

        badvalues = np.array(values, dtype=complex)
        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='slinear', points=badpoints, values=badvalues.tolist())

        msg = ("Interpolation method 'slinear' does not support complex values.")
        self.assertEqual(str(cm.exception), msg)

        interp = InterpND(method='slinear', points=points, values=values.tolist())
        x = [0.5, 0, 0.5, 0.9]

        with self.assertRaises(KeyError) as cm:
            interp = InterpND(method='slinear', points=points, values=values.tolist(), bad_arg=1)

        msg = ("\"InterpLinear: Option 'bad_arg' cannot be set because it has not been declared.")
        self.assertTrue(str(cm.exception).startswith(msg))

        # Bspline not supported for tables.
        points, values, func, df = self. _get_sample_2d()
        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='bsplines', points=points, values=values)

        msg = "Method 'bsplines' is not supported for table interpolation."
        self.assertTrue(str(cm.exception).startswith(msg))

        with self.assertRaises(ValueError) as cm:
            interp = InterpND(method='bsplines', points=points)

        msg = "Either 'values' or 'x_interp' must be specified."
        self.assertTrue(str(cm.exception).startswith(msg))

    def test_derivative_hysteresis_bug(self):
        alt = np.array([-1000, 0, 1000], dtype=float)
        rho = np.array([0.00244752, 0.00237717, 0.00230839])

        rho_interp = InterpND(method='slinear', points=alt, values=rho, extrapolate=True)

        x = 0.0
        _, dval1 = rho_interp.interpolate([x], compute_derivative=True)

        x = 0.5
        _, dval2 = rho_interp.interpolate([x], compute_derivative=True)

        x = 0.0
        _, dval3 = rho_interp.interpolate([x], compute_derivative=True)

        x = -0.5
        _, dval4 = rho_interp.interpolate([x], compute_derivative=True)

        x = 0.0
        _, dval5 = rho_interp.interpolate([x], compute_derivative=True)

        assert_near_equal(dval3 - dval1, np.array([[0.0]]))
        assert_near_equal(dval5 - dval1, np.array([[0.0]]))

if __name__ == '__main__':
    unittest.main()

max_stars_repo_name: kaushikponnapalli/OpenMDAO

max_issues_repo_name: kaushikponnapalli/OpenMDAO

max_forks_repo_name: kaushikponnapalli/OpenMDAO