text:
"""
Tests the various implementations of Restricted Boltzmann Machines.
"""

import numpy as np
import bernoulli_lib as blib
import vector_lib as vlib
from bernoulli_rbm import StandardBernoulliRBM

if "__main__" == __name__:
    # Test case for BRBM: outputs = inputs
    VEC_SIZE = 5
    X = blib.binary_matrix(range(2 ** VEC_SIZE), VEC_SIZE)

    # Hand-build a model designed to replicate the inputs
    a = -0.5 * np.ones(VEC_SIZE)
    b = -0.5 * np.ones(VEC_SIZE)
    W = np.diag(np.ones(VEC_SIZE))
    gen_model = StandardBernoulliRBM(params=(a, W, b), n_iter=25, n_report=1)

    # Test that the generative model reconstructs the inputs
    q = gen_model.predict(X)
    Y = blib.binary_decision(q)
    assert np.all(Y == X)
    print("Hand-crafted model:")
    gen_model.fit(X)

    # Learn random models of the inputs
    print("Independent training runs:")
    for i in range(10):
        lrn_model = StandardBernoulliRBM(VEC_SIZE)
        lrn_model.fit(X)
        report = lrn_model._report_card(**lrn_model.score(X))
        print("Trained random model %d - %s" % (i + 1, report))

ext: py

reason: import numpy


text:
'''Classes for interfacing with metadata on GEO
'''
import os
import warnings
from rpy2.rinterface import RRuntimeError

with warnings.catch_warnings():
	warnings.simplefilter("ignore")
	try:
		import geo_query as gq
	except (ImportError, RRuntimeError) as e:
		pass

import pandas as pd
import numpy as np

SCRIPT_DIR = os.path.dirname(os.path.realpath(__file__))

class GSE(object):
	"""docstring for GSE"""
	coll = 'geo'
	id_field = 'geo_accession'
	def __init__(self, geo_id):
		self.id = geo_id
		self.meta = {}
		self.GSMs = _GSMList(self.id)

	def retrieve(self):
		'''Retrieve metadata of GSE and GSMs from GEO'''
		with warnings.catch_warnings():
			warnings.simplefilter("ignore")
			gse_obj = gq.download_gse(self.id)
			gse_meta = gq.retrieve_gse_meta(gse_obj)
			gse_meta = gq.convert_r_list(gse_meta)
			# make sure the id_field exist
			gse_meta[self.id_field] = self.id
			self.meta = gse_meta
			self.GSMs.retrieve(gse_obj, gse_meta)

	def construct_sample_meta_df(self):
		'''Construct a dataframe for metadata based on GSMs'''
		meta_df = self.GSMs.construct_sample_meta_df()
		meta_df = meta_df.loc[self.meta['sample_id']]
		return meta_df

	def save(self, db):
		'''Save to database'''
		insert_result = db[self.coll].insert_one(self.meta)
		self.GSMs.save(db)
		return insert_result.inserted_id

	@classmethod
	def load(cls, geo_id, db):
		'''Load from the database.'''
		doc = db[cls.coll].find_one({cls.id_field: geo_id}, {'_id':False})
		obj = cls(geo_id)
		obj.meta = doc
		obj.GSMs = _GSMList.load(obj, db)
		return obj

def unpack_gse_meta(d):
	parsed_d = {}
	for key, value in d.items():
		if type(value) == list:
			for val in value:
				if val.count(': ') == 1:
					k, v = val.split(': ')
					parsed_d[k] = v
		else:
			parsed_d[key] = value
	return parsed_d

class _GSMList(object):
	"""docstring for _GSMList, this class will only be used internally for GSE"""
	coll = 'gsm'
	id_field = 'geo_accession'
	def __init__(self, gse_id):
		self.gse_id = gse_id
		self.GSMs = []

	def retrieve(self, gse_obj, gse_meta):
		'''Retrieve GSM metas from gse_obj'''
		gsm_metas = gq.retrieve_gsms_meta_from_gse(gse_obj)
		gsm_metas = gq.convert_gsms_meta_r_list(gsm_metas)
		gsm_metas = gq.drop_duplicated_meta_in_gsm(gse_meta, gsm_metas)

		for gsm_id, gsm_meta in gsm_metas.items():
			gsm_meta[self.id_field] = gsm_id
			self.GSMs.append(gsm_meta)

	def __len__(self):
		return len(self.GSMs)

	def construct_sample_meta_df(self):
		''''''
		# unpack GSM meta fields that are lists
		meta_df = pd.DataFrame.from_records(map(unpack_gse_meta, self.GSMs))\
			.set_index(self.id_field)
		# remove columns with only 1 unique value
		meta_df = meta_df.loc[:, meta_df.nunique() > 1]
		# replace column names with '.'
		meta_df.columns = meta_df.columns.map(lambda x:x.replace('.', '_'))
		# attempt to convert columns to numerical type
		meta_df.replace('NA', np.nan, inplace=True)
		for col in meta_df.columns:
			meta_df[col] = pd.to_numeric(meta_df[col], errors='ignore')
		return meta_df

	def save(self, db):
		'''Save to database'''
		insert_result = db[self.coll].insert_many(self.GSMs)
		return insert_result.inserted_ids

	@classmethod
	def load(cls, gse, db):
		'''Load from the database. Need a GSE instance as param.'''
		cur = db[cls.coll].find({cls.id_field: {'$in': gse.meta['sample_id']}}, {'_id': False})
		obj = cls(gse.id)
		obj.GSMs = [doc for doc in cur]

		return obj

max_stars_repo_name: MaayanLab/scAVI

max_issues_repo_name: MaayanLab/scAVI

max_forks_repo_name: MaayanLab/scAVI


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

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

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

    noise_ratio: float
        Ratio of labels to invert

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

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

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

    # Question 1: Train- and test errors of AdaBoost in noiseless case
    model = AdaBoost(DecisionStump, n_learners).fit(train_X, train_y)
    train_errs, test_errs = [], []
    learners = np.arange(1, n_learners)
    for iteration in learners:
        train_errs.append(model.partial_loss(train_X, train_y, iteration))
        test_errs.append(model.partial_loss(test_X, test_y, iteration))

    make_subplots() \
        .add_traces([go.Scatter(x=learners, y=test_errs, mode='lines',
                                marker=dict(color='red', line_width=1), name='Test error'),
                     go.Scatter(x=learners, y=train_errs, mode='lines',
                                marker=dict(color='green', line_width=1), name='Train error')]) \
        .update_layout(xaxis_title="Number of weak learners", yaxis_title="Error rate",
                       title="Error rate as an output to number of weak learners") \
        .show()

    # Question 2: Plotting decision surfaces
    T = [5, 50, 100, 250]
    limits = np.array([np.r_[train_X, test_X].min(axis=0), np.r_[train_X, test_X].max(axis=0)]).T + np.array([-.1, .1])
    fig = make_subplots(rows=2, cols=2, subplot_titles=[f"Number of weak learners - {n}" for n in T],
                        horizontal_spacing=0.05, vertical_spacing=.2)

    for ind, t in enumerate(T):
        fig.add_traces([decision_surface(lambda x: model.partial_predict(x, t), limits[0], limits[1],
                                         showscale=False, colorscale=["red", "blue"]),
                        go.Scatter(x=train_X[:, 0], y=train_X[:, 1], mode="markers", showlegend=False,
                                   marker=dict(color=train_y, colorscale=["red", "blue"],
                                               line=dict(color="black", width=1)))],
                       rows=(ind // 2) + 1, cols=(ind % 2) + 1) \
            .update_layout(title='decision surface as an output to number of weak learners',
                           margin=dict(t=100)).update_xaxes(visible=False).update_yaxes(visible=False)
    fig.show()

    # Question 3: Decision surface of best performing ensemble
    test_error = []
    for T in learners:
        test_error.append(model.partial_loss(test_X, test_y, T))
    best_learner_idx = np.argmin(test_errs) + 1

    make_subplots() \
        .add_traces(
        [decision_surface(lambda x: model.partial_predict(x, best_learner_idx), limits[0], limits[1],
                          showscale=False),
         go.Scatter(x=test_X[:, 0], y=test_X[:, 1], mode="markers", showlegend=False,
                    marker=dict(color=test_y,
                                colorscale=[custom[0], custom[-1]],
                                line=dict(color="black", width=1)))]) \
        .update_layout(title=f"Decision Surface as an output to number of weak learners = {best_learner_idx}, accuracy is "
                             f"{1 - test_error[best_learner_idx]}",
                       font=dict(family="Arial", size=20),
                       margin=dict(t=50)) \
        .show()

    # Question 4: Decision surface with weighted samples
    D = model.D_ / np.max(model.D_) * 5
    make_subplots().add_traces([decision_surface(model.predict, limits[0], limits[1], showscale=False),
         go.Scatter(x=train_X[:, 0], y=train_X[:, 1], mode="markers", showlegend=False,
                    marker=dict(colorscale=[custom[0], custom[-1]], color=train_y,
                                size=D, line=dict(color="black", width=1)))]) \
        .update_layout(margin=dict(t=100)) \
        .show()

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


text:
import os
import glob
import time
import numpy as np
from PIL import Image
from pathlib import Path
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
from skimage.color import rgb2lab, lab2rgb

import torch
from torch import nn, optim
from torchvision import transforms
from torchvision.utils import make_grid
from torch.utils.data import Dataset, DataLoader

from dataset import *
#from dataset import make_dataloaders
from loss import *
from  models import *
from utils import *
from pre_train import *

#device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

def train_model(model, train_dl, epochs, display_every=200):
    data = next(iter(val_dl)) # getting a batch for visualizing the model output after fixed intrvals
    for e in range(epochs):
        loss_meter_dict = create_loss_meters() # function returing a dictionary of objects to
        i = 0                                  # log the losses of the complete network
        for data in train_dl:
            model.setup_input(data)
            model.optimize()
            update_losses(model, loss_meter_dict, count=data['L'].size(0)) # function updating the log objects
            i += 1
            if i % display_every == 0:
                print(f"\nEpoch {e+1}/{epochs}")
                print(f"Iteration {i}/{len(train_dl)}")
                log_results(loss_meter_dict) # function to print out the losses
                #visualize(model, data, save=False) # function displaying the model's outputs

def run_experiment(resnet, learningRate,weight_decay):
    train_dl = make_dataloaders(batch_size=8, n_workers=1,paths=train_paths, split='train')
    val_dl = make_dataloaders(batch_size=8, n_workers=1, paths=val_paths, split='val')

    net_G = build_res_unet(resnet, n_input=1, n_output=2, size=256)
    opt = optim.Adam(net_G.parameters(), lr=learningRate,weight_decay=weight_decay)
    criterion = nn.L1Loss()
    pretrain_generator(net_G, train_dl,val_dl, opt, criterion, 20, resnet, learningRate,weight_decay)

#net_G = build_res_unet(n_input=1, n_output=2, size=256)
#opt = optim.Adam(net_G.parameters(), lr=1e-4)
#criterion = nn.L1Loss()
#pretrain_generator(net_G, train_dl,val_dl, opt, criterion,40)
#torch.save(net_G.state_dict(), "res18-unet.pt")

#net_G = build_res_unet(n_input=1, n_output=2, size=256)
#net_G.load_state_dict(torch.load("res18-unet.pt", map_location=device))
#model = MainModel(net_G=net_G)
#train_model(model, train_dl, 1)


text:
import numpy as np
from scipy import special

sqrtpi = np.sqrt(np.pi)
sqrt2 = np.sqrt(2)
ln10 = np.log(10)
log2onsqrtpi = np.log(2.0/sqrtpi)
log2sqrt2 = np.log(2*sqrt2)
ln2 = np.log(2)

def erf(x,deriv = False):

    y = special.erf(x)
    if not deriv: return y

    return y, lambda dy: (2.0/sqrtpi)*dy*np.exp(-x**2)

def erf_and_logderiv(x, deriv = False):

    y = special.erf(x)
    logderiv = log2onsqrtpi - x**2

    if not deriv: return y, logderiv

    def back(dy, dlogderiv):
        return dy*np.exp(logderiv) - 2*dlogderiv*x

    return y, logderiv, back

def cs_erfinv(p):
    """Complex-step friendly version of scipy.special.erfinv.

    For complex p, it returns the complex number:

        x + j*dx

    where x = erfinv(real(p)) and dx = dx/dp * imag(p)
    """
    if not np.iscomplexobj(p):
        return special.erfinv(p)
    rp = np.real(p)
    x = special.erfinv(rp)
    dp = (sqrtpi/2.0)*np.exp(x**2)*np.imag(p)
    return x + 1.0j*dp

def erfinv(p,deriv = False):

    if not deriv: return cs_erfinv(p)

    x = cs_erfinv(p)
    back = lambda dy: (sqrtpi/2.0)*dy*np.exp(x**2)
    return x, back

def Phi(x,deriv = False):
    """Normal cumulative distribution."""
    arg = x/sqrt2
    if not deriv: return (1 + erf(arg) ) / 2.0

    e, back1 = erf(arg,deriv=True)
    y = (1+e)/2
    def back(dy): return back1(dy/2.0)/sqrt2
    return y, back

def Phi_and_logderiv(x,deriv = False):

    arg = x/sqrt2
    #if deriv: e, logderiv1, back1 = erf_and_logderiv(arg, deriv = deriv)
    e, logderiv1 = erf_and_logderiv(arg)
    y = (1+e)/2.0
    logderiv = logderiv1 - log2sqrt2
    if not deriv: return y, logderiv

    def back(dy, dlogderiv):
        return dy*np.exp(logderiv) -sqrt2*dlogderiv*arg
        #return back1(dy/2.0,dlogderiv)/sqrt2

    return y, logderiv, back

def Phi_inv(p,deriv = False):

    if not deriv: return sqrt2*erfinv(2.0*p - 1.0)

    arg, back1 = erfinv(2.0*p - 1.0, deriv=True)
    x = sqrt2*arg
    def back(dx): return 2.0*back1(sqrt2*dx)
    return x, back

probit = Phi_inv

def dprime2EER(dprime,deriv = False):

    if not deriv: return Phi(-0.5*dprime)

    EER, back1 = Phi(-0.5*dprime,deriv=True)
    def back(dEER): return -0.5*back1(dEER)
    return EER, back

def logdprime2EER(logdprime,deriv = False):

    arg = -0.5*np.exp(logdprime)
    if not deriv: return Phi(arg)

    EER, back1 = Phi(arg,deriv=True)
    def back(dEER): return arg*back1(dEER)
    return EER, back

def logdprime2EER_and_logabsderiv(logdprime,deriv = False):

    arg = -0.5*np.exp(logdprime)
    if not deriv: EER, logderiv1 = Phi_and_logderiv(arg)
    else: EER, logderiv1, back1 = Phi_and_logderiv(arg, deriv = True)
    logderiv = logderiv1 + logdprime - ln2
    if not deriv: return EER, logderiv

    def back(dEER, dlogderiv):
        return dlogderiv + arg*back1(dEER,dlogderiv)

    return EER, logderiv, back

def EER2logdprime(EER,deriv = False):

    if not deriv: arg = Phi_inv(EER)
    else: arg, back1 = Phi_inv(EER,deriv=True)
    dprime = -2.0*arg
    logdprime = np.log(dprime)
    if not deriv: return logdprime
    def back(dlogprime): return back1(-2.0*dlogprime/dprime)
    return logdprime, back

def EER2dprime(EER,deriv = False):

    if not deriv: arg = Phi_inv(EER)
    else: arg, back1 = Phi_inv(EER,deriv=True)
    dprime = -2.0*arg
    if not deriv: return dprime
    def back(ddprime): return back1(-2.0*ddprime)
    return dprime, back

def reparam_pdf(x_pdf, y2x, y, cs = True):
    """
    Scalar PDF reparametrization.

    If y = f(x), where x ~ x_pdf
    then this function returns y_pdf(y).

    params:
        x_pdf: the probability density function for x
        y2x: the inverse function of f, equipped to compute derivatives
        y: a scalar or vector at which to evaluate y_pdf.

    """
    if np.iscomplexobj(y) or not cs:
        x, back = y2x(y, deriv=True)
        return x_pdf(x) * np.abs(back(1.0))
    else:
        cx = y2x(y + 1e-20j)
        x, dxdy = np.real(cx) , 1e20*np.imag(cx)
        return x_pdf(x) * np.abs(dxdy)

if __name__ == "__main__":
    print("Running test script for module erf_tools\n")

    from numpy.random import randn, rand
    from lib.deriv.adtools import cstest

    print("\nTesting erf:")
    x = randn(5)
    delta = cstest(erf,x)
    print(delta)

    print("\nTesting Phi:")
    x = randn(5)
    delta = cstest(Phi,x)
    print(delta)

    print("\nTesting dprime2EER:")
    x = randn(5)
    delta = cstest(dprime2EER,x)
    print(delta)

    print("\nTesting EER2dprime:")
    EER = rand(5)/2
    delta = cstest(EER2dprime,EER)
    print(delta)
    EER2 = dprime2EER(EER2dprime(EER))
    print(abs(EER-EER2).max())

    print("\nTesting logdprime2EER:")
    x = randn(5)
    delta = cstest(logdprime2EER,x)
    print(delta)

    print("\nTesting erfinv:")
    p = rand(5)
    delta = cstest(erfinv,p)
    print(delta)

    print("\nTesting Phi_inv:")
    p = rand(5)
    delta = cstest(Phi_inv,p)
    print(delta)

    print("\nTesting EER2logdprime:")
    p = rand(5)/2
    delta = cstest(EER2logdprime,p)
    print(delta)

    print("\nTesting erf_and_logderiv:")
    x = randn(5)
    delta = cstest(erf_and_logderiv,x)
    print(delta)

    print("\nTesting Phi_and_logderiv:")
    x = randn(5)
    delta = cstest(Phi_and_logderiv,x)
    print(delta)

    print("\nTesting logdprime2EER_and_logabsderiv:")
    x = randn(5)
    delta = cstest(logdprime2EER_and_logabsderiv,x)
    print(delta)
    EER1, logabsderiv = logdprime2EER_and_logabsderiv(x)
    EER2, back = logdprime2EER(x,deriv=True)
    print(abs(EER1-EER2).max())
    print(abs(logabsderiv-np.log(-back(1))).max())

    if False:

        import matplotlib.pyplot as plt
        fig, ax = plt.subplots(1, 1)
        from scipy import stats

    #    y = np.linspace(-5,EER2logdprime(0.000001),200)
    #    pdf_x = stats.beta(1,1,scale=0.5).pdf      #uniform between 0 and 0.5
    #    y2x = logdprime2EER
    #    pdf_y = reparam_pdf(pdf_x,y2x,y)
    #    ax.plot(y,pdf_y,label="log(d')")

    #    y = np.linspace(np.exp(-5),np.exp(EER2logdprime(0.000001)),200)
    #    pdf_x = stats.beta(1,4,scale=0.5).pdf      #uniform between 0 and 0.5
    #    y2x = dprime2EER
    #    pdf_y = reparam_pdf(pdf_x,y2x,y)
    #    ax.plot(y,pdf_y,label="d'")

        y = np.linspace(0,0.49,200)
        pdf_x = stats.invgauss(mu = 1, scale = 1).pdf
        y2x = EER2dprime
        pdf_y = reparam_pdf(pdf_x,y2x,y)
        ax.plot(y,pdf_y,label="EER")

        ax.legend(loc='best', frameon=False)
        plt.xlabel("EER")
        plt.grid()
        plt.show()

reason: import numpy,from numpy,from scipy


text:
"""
Test cases for aggregate_project

E.g. to run from Jupyter enter

!python test_cases.py

"""

import sys

sys.path.insert(0,'/s/telos/python/aggregate_project/')

from aggregate import Aggregate, Severity, Portfolio, Underwriter
import numpy as np
import unittest

# import warnings
# warnings.filterwarnings("ignore")

def kword_expander(**kwargs):
    pass

class TestAggregateModule(unittest.TestCase):

    def test_Aggregate(self):
        sig = .75
        d1b1 = Aggregate('single', exp_en=10, sev_name='lognorm', sev_a=[sig, sig / 2, sig / 5],
                         sev_scale=100 * np.exp(-sig ** 2 / 2), sev_wt=1, freq_name='gamma', freq_a=0.81)
        xs = np.linspace(0, 1023, 1024)
        d1b1.update(xs, verbose=False)
        self.assertTrue(d1b1.statistics_total_df.T.loc['el', :].sum() == 5146.408347357843)

        sig = 0.75
        d1b2 = Aggregate('sublines', exp_en=[10, 10, 10, 10], exp_attachment=[0, 0, 100, 50],
                         exp_limit=[np.inf, 100, np.inf, 100], sev_name='lognorm', sev_a=sig,
                         sev_scale=100 * np.exp(-sig ** 2 / 2),
                         sev_wt=1, freq_name='gamma', freq_a=0.81)
        self.assertTrue(d1b2.statistics_total_df.T.loc['el', :].sum() == 6152.9037739696396)

        d2b1 = Aggregate('mixed', exp_en=50, exp_attachment=[0, 50, 100, 150, 200],
                         exp_limit=[250, 250, np.inf, 100, 5000],
                         sev_name='lognorm', sev_mean=100, sev_cv=np.linspace(0.05, 2, 5), sev_wt=0.2 * np.ones(5),
                         freq_name='gamma', freq_a=1.2)
        self.assertTrue(d2b1.statistics_total_df.T.loc['el', :].sum() == 42977.102931408575)

        fixed = Aggregate('chistogram', exp_en=1, sev_name='dhistogram', sev_xs=[0, 1, 2, 3, 4], sev_ps=.2,
                          freq_name='fixed')
        self.assertTrue(fixed.statistics_df['sev_cv'].sum() +
                        fixed.statistics_total_df['sev_cv'].sum()== 2.1213203435596428)

    def test_Severity(self):
        fixed = Severity('dhistogram', sev_xs=[0, 1, 2, 3, 4], sev_ps=[.2, .3, .4, .05, .05])
        self.assertTrue(fixed.moms() == (1.45, 3.1500000000000004, 8.0500000000000007))

        fixed = Severity('chistogram', sev_xs=[0, 1, 2, 3, 4], sev_ps=[.1, .2, .3, 0, .4])
        self.assertTrue(fixed.moms() == (2.9000000000000004, 10.449999999999999, 41.825000000000003))

        fixed = Severity('fixed', sev_xs=2)
        self.assertTrue(fixed.moms() == (2, 4, 8))

    def test_big_example(self):
        # bigger agg example
        sa = Aggregate('my test', exp_en=np.ones(50), exp_premium=0, exp_lr=.75,
                       exp_attachment=20, exp_limit=np.linspace(10, 500, 50), sev_name='lognorm', sev_mean=100,
                       sev_cv=0.8, freq_name='poisson', freq_a=0)

        bs = 1  # 2 and 13
        log2 = 13
        N = 1 << log2
        MAXL = N * bs
        xs = np.linspace(0, MAXL, N, endpoint=False)
        audit = sa.update(xs, 1, None, 'exact', sev_calc='discrete', discretization_calc='survival', verbose=True)
        # average of square errors is small:
        self.assertTrue(
            np.sum(np.abs(audit.iloc[0:-3, :]['rel sev err']) ** 2) ** .5 / (len(audit) - 3) < 1e-6)
        self.assertTrue(
            np.sum(np.abs(audit.iloc[0:-3, :]['abs sev err']) ** 2) ** .5 / (len(audit) - 3) < 1e-4)

    def test_excel_read(self):
        port = Portfolio.from_Excel('example', 'Examples.xlsx', 'Generic')
        port.update(11, 50)
        # overall mean error is less than 0.001
        self.assertTrue( port.audit_df['MeanErr'].abs().sum() < 0.001)

class TestUnderwriter(unittest.TestCase):
    def setUp(self):
        self.uw = Underwriter()

    def test_severity(self):
        cc = self.uw('liaba')
        self.assertTrue(np.allclose(cc.stats(), (50, 2500)))

    def test_portfolio(self):
        p = self.uw('THREE~LINE~EXAMPLE')
        p.update(14, 1)
        self.assertTrue(np.all(p.audit_df.loc['total', ['MeanErr', 'CVErr']].abs() < 2e-4))

        self.assertTrue(p.audit_df.MeanErr.abs().sum() < 1.5e-4)
        self.assertTrue(p.audit_df.CVErr.abs().sum() < 5e-4)
        a, p, test, params, dd, table, stacked = p.uat()
        self.assertTrue(a['lr err'].abs().sum() < 1e-8)
        self.assertTrue(np.all(test.filter(regex='err[_s]', axis=1).abs().sum() < 1e-8))

    def test_parser(self):
        portfolio_program = """port test~portfolio
    agg big~mixture  50 claims              [50, 100, 150, 200] xs 0  sev lognorm 12 cv [1,2,3,4] wts [0.25 .25 .25 .25]  poisson
    agg A1a          500 premium at 0.5                               sev gamma 12 cv .30   poisson
    agg A1b          500 premium at 0.5 lr                            sev gamma 12 cv .30   poisson
    agg A2           50  claims             30 xs 10                  sev gamma 12 cv .30   poisson
    agg A3           50  claims                                       sev gamma 12 cv .30   poisson
    agg hcmp         1e-8 * agg.CMP
            """
        p = self.uw.write(portfolio_program, update=True, log2=13, remove_fuzz=True)
        self.assertTrue(np.all(p.audit_df.iloc[:-1, :].CVErr.abs() < 0.005))
        self.assertTrue(np.all(p.audit_df.iloc[:-1, :].MeanErr.abs() < 0.001))

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

max_stars_repo_name: mynl/aggregate

max_issues_repo_name: mynl/aggregate

max_forks_repo_name: mynl/aggregate


text:
#!/usr/bin/env python3
"""
Display some information like pie charts...etc... in order to analyse the
results of the evaluation part of the algorithm.
"""

import tensorflow as tf
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import yaml
import numpy as np
import preprocessing as pp

# Constants
# ==================================================
matplotlib.rcParams['font.size'] = 5.0

# Definitions
# ==================================================

def pie_chart_support_distribution(classification_report, title, folder):
    """
    Plot a pie chart which describes the distribution of each class.
    :param classification_report: Sliced classification report : classes,
    toPlot, support. toPlot must be a tuple (precision, recall, f1-score)
    """

    classes, toPlot, support = slice_classification_report(
            classification_report)

    # Don't take into account the last column which is the total number
    # of each class
    labels = classes[0:len(classes)-1]
    sizes = support[0:len(classes)-1]

    fig1, ax1 = plt.subplots()
    patches, texts, _ = ax1.pie(sizes, labels=labels, autopct='%1.1f%%',
                                startangle=90)
    # Equal aspect ratio ensures that pie is drawn as a circle.
    ax1.axis('equal')
    ax1.set_title(title)
    ax1.legend(patches, labels, loc="best")

    plt.savefig(folder+"/"+title.replace(" ", "_")+".png", format="png",
                dpi=1000)

def bar_chart_classification_report(classification_report, title, folder):
    """
    Plot a bar graph which sums up the classification report of the scikit
    learn tool.
    :param classification_report: Sliced classification report : classes,
    toPlot, support. toPlot must be a tuple (precision, recall, f1-score)
    """
    classes, toPlot, support = slice_classification_report(
            classification_report)

    N = 3
    bar_width = 0.05

    ind = np.arange(N)
    fig, ax = plt.subplots()

    # Enumerate over each class except the last one which represent the average
    # and total
    bars = []
    for i in range(len(classes)):
        bar_i = ax.bar(ind + i * bar_width, toPlot[i], bar_width)
        bars.append(bar_i)

    # Add some text for labels, title and axes ticks
    ax.set_ylabel("Percent")
    ax.set_title(title)
    ax.set_xticks(ind + bar_width / len(classes))
    ax.set_xticklabels(("Precision", "Recall", "F1-score"))

    ax.legend(bars, classes, loc="best")

    plt.savefig(folder+"/"+title.replace(" ", "_")+".png", format="png",
                dpi=1000)

def slice_classification_report(classification_report):
    """
    Plot scikit-learn classification report.
    Extension based on https://stackoverflow.com/a/31689645/395857
    """
    lines = classification_report.split('\n')

    classes = []
    plotMat = []
    support = []
    class_names = []

    for line in lines[2: (len(lines) - 2)]:
        t = line.strip().split()

        if len(t) < 2:
            continue

        classes.append(t[0])
        v = [float(x) for x in t[1: len(t) - 1]]
        support.append(int(t[-1]))
        class_names.append(t[0])
        plotMat.append(v)

    # Save the average precision/recall/F1-score and total support
    t = lines[len(lines) - 2].strip().split()
    classes.append(t[0] + t[1] + t[2])
    v = [float(x) for x in t[3: len(t) - 1]]
    support.append(int(t[-1]))
    class_names.append(t[0] + t[1] + t[2])
    plotMat.append(v)

    print("\n")
    print("plotMat: {0}".format(plotMat))
    print("support: {0}".format(support))

    return classes, plotMat, support

def display_stat(filepath):
    """
    Statistics from the SemEval 2016 competition, Task 5, Subtask 1 dataset.
    :param filepath: Path of the dataset SemEval. The path must leads to a
    folder containing both the training and testing sets.
    :type filepath: string
    :return: Pandas.dataframe with the following columns : review_id,
    sentence_id, text, feature, polarity
    """

    training_set = pp.parse_XML(filepath+"/train.xml")
    testing_set = pp.parse_XML(filepath+"/test/test_gold.xml")

    # Some opinions concerns various food, drinks...etc... but the opinion
    # is the same while the target differ. So deleting duplicates as the scope
    # of this study does not imply target (OPE in SemEval)
    training_set = training_set.drop_duplicates()
    testing_set = testing_set.drop_duplicates()

    # Count # of opinions for each sentence
    count_opinions_train = training_set['sentence_id'].value_counts()
    count_opinions_train = count_opinions_train.value_counts()
    count_opinions_test = testing_set['sentence_id'].value_counts()
    count_opinions_test = count_opinions_test.value_counts()

    # Display pie charts
    count_dict_train = count_opinions_train.to_dict()
    labels = list(count_dict_train.keys())
    sizes = list(count_dict_train.values())

    fig, ax = plt.subplots()
    ax.pie(sizes, labels=labels, autopct='%1.1f%%', startangle=90)
    ax.axis('equal')
    ax.set_title('[TRAIN] Percentage of opinion occurences in a sentence')

    count_dict_test = count_opinions_test.to_dict()
    labels = list(count_dict_test.keys())
    sizes = list(count_dict_test.values())

    fig1, ax1 = plt.subplots()
    ax1.pie(sizes, labels=labels, autopct='%1.1f%%', startangle=90)
    ax1.axis('equal')
    ax1.set_title('[TEST] Percentage of opinion occurences in a sentence')

    plt.show()

def display_pie(data, folder):

    # Feature distribution
    # --------------------

    distribution_train = data.groupby('feature').size()
    labels = distribution_train.index.tolist()
    sizes = distribution_train.tolist()

    fig, ax = plt.subplots()
    patches, texts, autotexts = ax.pie(sizes, labels=labels, autopct='%1.1f%%',
                               startangle=90)
    # Equal aspect ratio ensures that pie is drawn as a circle.
    ax.axis('equal')

    this_folder = folder+"/feature_distribution"
    plt.savefig(this_folder+".png", format="png",
                dpi=1000)

    # Polarity distribution
    # ---------------------

    distribution_train = data.groupby('polarity').size()
    labels = distribution_train.index.tolist()
    sizes = distribution_train.tolist()

    fig, ax = plt.subplots()
    patches, texts, autotexts = ax.pie(sizes, labels=labels, autopct='%1.1f%%',
                               startangle=90)
    # Equal aspect ratio ensures that pie is drawn as a circle.
    ax.axis('equal')

    this_folder = folder+"/polarity_distribution"
    plt.savefig(this_folder+".png", format="png",
                dpi=1000)

    # Combination distribution
    # ------------------------

    distribution_train = data.groupby(['feature', 'polarity']).size()
    labels = distribution_train.index.tolist()
    sizes = distribution_train.tolist()

    fig, ax = plt.subplots()
    patches, texts, autotexts = ax.pie(sizes, labels=labels, autopct='%1.1f%%',
                               startangle=90)
    # Equal aspect ratio ensures that pie is drawn as a circle.
    ax.axis('equal')

    this_folder = folder+"/comb_distribution"
    plt.savefig(this_folder+".png", format="png",
                dpi=1000)

def display_distrib_data():

    folder = "Figures/Data_distribution"

    # Aspect solution
    # ===============

    # Restaurant domain
    filepath = "../data/SemEval/Subtask1/restaurant"

    training_set = pp.parse_XML(filepath+"/train.xml", True)
    testing_set = pp.parse_XML(filepath+"/test/test_gold.xml", True)

    training_set = training_set.drop_duplicates()
    testing_set = testing_set.drop_duplicates()

    display_pie(training_set, folder+"/Train/Aspect/restaurant")
    display_pie(testing_set, folder+"/Test/Aspect/restaurant")

    # Laptop domain
    filepath = "../data/SemEval/Subtask1/laptop"

    training_set = pp.parse_XML(filepath+"/train.xml", True)
    testing_set = pp.parse_XML(filepath+"/test/test_gold.xml", True)

    training_set = training_set.drop_duplicates()
    testing_set = testing_set.drop_duplicates()

    display_pie(training_set, folder+"/Train/Aspect/laptop")
    display_pie(testing_set, folder+"/Test/Aspect/laptop")

    # Entity solution
    # ===============

    # Restaurant domain
    filepath = "../data/SemEval/Subtask1/restaurant"

    training_set = pp.parse_XML(filepath+"/train.xml")
    testing_set = pp.parse_XML(filepath+"/test/test_gold.xml")

    training_set = training_set.drop_duplicates()
    testing_set = testing_set.drop_duplicates()

    display_pie(training_set, folder+"/Train/Entity/restaurant")
    display_pie(testing_set, folder+"/Test/Entity/restaurant")

    # Laptop domain
    filepath = "../data/SemEval/Subtask1/laptop"

    training_set = pp.parse_XML(filepath+"/train.xml")
    testing_set = pp.parse_XML(filepath+"/test/test_gold.xml")

    training_set = training_set.drop_duplicates()
    testing_set = testing_set.drop_duplicates()

    display_pie(training_set, folder+"/Train/Entity/laptop")
    display_pie(testing_set, folder+"/Test/Entity/laptop")

def slot3_accuracy():

    # Filepath for restaurant domain predictions
    filepath = "runs/aspects/word2vec_200-epochs/1501773565/predictions.csv"

    # Load the predictions into a DataFrame
    predictions = pd.DataFrame.from_csv(filepath)
    predictions = predictions[['text', 'feature', 'polarity', 'pred_polarity']]

    # For each aspect, determine the accuracy
    # ---------------------------------------
    predictions = predictions.groupby('feature')

    acc = {}
    for key, item in predictions:
        predictions_one_aspect = predictions.get_group(key)
        total_number = len(predictions_one_aspect)
        df = pd.DataFrame(columns=['correct_polarities'])
        df['correct_polarities'] = (predictions_one_aspect.polarity == predictions_one_aspect.pred_polarity)
        correct_polarities = df['correct_polarities'].sum()
        acc[key] = float(correct_polarities)/total_number

    # Get the mean accuracy
    mean_acc = float(sum(acc.values())) / len(acc)
    acc["mean"] = mean_acc

    print("Slot 3 - Accuracy measure :")
    for key, value in acc.items():
        print(key + " : " + str(value))

if __name__ == '__main__':

    with open("config.yml", 'r') as ymlfile:
        cfg = yaml.load(ymlfile)

    # Parameters
    # ==================================================

    # Data Parameters

    # Eval Parameters
    tf.flags.DEFINE_boolean("display_stat", False,
                            "Display statistics of SemEval dataset")

    tf.flags.DEFINE_boolean("slot3", False,
                            "Display accuracy following the accuracy measure "+
                            "of SemEval competition")

    FLAGS = tf.flags.FLAGS
    FLAGS._parse_flags()

    dataset_filepath_REST = "../data/SemEval/Subtask1/restaurant"
    dataset_filepath_LAPT = "../data/SemEval/Subtask1/laptop"

    if FLAGS.display_stat:
        display_stat(dataset_filepath_REST)
        display_stat(dataset_filepath_LAPT)

    if FLAGS.slot3:
        slot3_accuracy()


text:
# Author: Pierre Ablin <pierre.ablin@ens.fr>
#
# License: MIT
import torch
import numpy as np
from .kernels import (imq_kernel, gaussian_stein_kernel_single,
                      linear_stein_kernel)

from scipy.optimize import fmin_l_bfgs_b
from time import time

def ksdd_gradient(x0, score, step, kernel='gaussian', max_iter=1000, bw=1,
                  store=False, verbose=False, clamp=None, beta=0.2):
    '''Kernel Stein Discrepancy descent with gradient descent

    Perform Kernel Stein Discrepancy descent with gradient descent.
    Since it uses gradient descent, a step size must be specified.

    Parameters
    ----------
    x0 : torch.tensor, size n_samples x n_features
        initial positions

    score : callable
        function that computes the score

    step : float
        step size

    max_iter : int
        max numer of iters

    bw : float
        bandwidth of the stein kernel

    stores : None or list of ints
        whether to stores the iterates at the indexes in the list

    verbose: bool
        wether to print the current loss

    clamp:
        if not None, should be a tuple (a, b). The points x are then
        constrained to stay in [a, b]

    Returns
    -------
    x: torch.tensor, size n_samples x n_features
        The final positions

    loss_list : list of floats
        List of the loss values during iterations

    References
    ----------
    A.Korba, P-C. Aubin-Frankowski, S.Majewski, P.Ablin.
    Kernel Stein Discrepancy Descent
    International Conference on Machine Learning, 2021.
    '''
    x = x0.clone().detach()
    n_samples, p = x.shape
    x.requires_grad = True
    if store:
        storage = []
        timer = []
        t0 = time()
    loss_list = []
    n = None
    for i in range(max_iter + 1):
        if store:
            timer.append(time() - t0)
            storage.append(x.clone())
        scores_x = score(x)
        if kernel == 'gaussian':
            K = gaussian_stein_kernel_single(x, scores_x, bw)
        else:
            K = imq_kernel(x, x, scores_x, scores_x, g=bw, beta=beta)
        loss = K.sum() / n_samples ** 2
        loss.backward()
        loss_list.append(loss.item())
        if verbose and i % 100 == 0:
            print(i, loss.item())
        with torch.no_grad():
            x[:] -= step * x.grad
            if n is not None:
                x[:] -= n
            if clamp is not None:
                x = x.clamp(clamp[0], clamp[1])
            x.grad.data.zero_()
        x.requires_grad = True
    x.requires_grad = False
    if store:
        return x, storage, timer
    else:
        return x

def ksdd_lbfgs(x0, score, kernel='gaussian', bw=1.,
               max_iter=10000, tol=1e-12, beta=.5,
               store=False, verbose=False):
    '''Kernel Stein Discrepancy descent with L-BFGS

    Perform Kernel Stein Discrepancy descent with L-BFGS.
    L-BFGS is a fast and robust algorithm, that has no
    critical hyper-parameter.

    Parameters
    ----------
    x0 : torch.tensor, size n_samples x n_features
        initial positions

    score : callable
        function that computes the score

    kernl : 'gaussian' or 'imq'
        which kernel to choose

    max_iter : int
        max numer of iters

    bw : float
        bandwidth of the stein kernel

    tol : float
        stopping criterion for L-BFGS

    store : bool
        whether to stores the iterates

    verbose: bool
        wether to print the current loss

    Returns
    -------
    x: torch.tensor, size n_samples x n_features
        The final positions

    References
    ----------
    A.Korba, P-C. Aubin-Frankowski, S.Majewski, P.Ablin.
    Kernel Stein Discrepancy Descent
    International Conference on Machine Learning, 2021.
    '''
    x = x0.clone().detach().numpy()
    n_samples, p = x.shape
    if store:
        class callback_store():
            def __init__(self):
                self.t0 = time()
                self.mem = []
                self.timer = []

            def __call__(self, x):
                self.mem.append(np.copy(x))
                self.timer.append(time() - self.t0)

            def get_output(self):
                storage = [torch.tensor(x.reshape(n_samples, p),
                                        dtype=torch.float32)
                           for x in self.mem]
                return storage, self.timer
        callback = callback_store()
    else:
        callback = None

    def loss_and_grad(x_numpy):
        x_numpy = x_numpy.reshape(n_samples, p)
        x = torch.tensor(x_numpy, dtype=torch.float32)
        x.requires_grad = True
        scores_x = score(x)
        if kernel == 'gaussian':
            stein_kernel = gaussian_stein_kernel_single(x, scores_x, bw)
        elif kernel == 'imq':
            stein_kernel = imq_kernel(x, x, scores_x, scores_x, bw, beta=beta)
        else:
            stein_kernel = linear_stein_kernel(x, x, scores_x, scores_x)
        loss = stein_kernel.sum()
        loss.backward()
        grad = x.grad
        return loss.item(), np.float64(grad.numpy().ravel())

    t0 = time()
    x, f, d = fmin_l_bfgs_b(loss_and_grad, x.ravel(), maxiter=max_iter,
                            factr=tol, epsilon=1e-12, pgtol=1e-10,
                            callback=callback)
    if verbose:
        print('Took %.2f sec, %d iterations, loss = %.2e' %
              (time() - t0, d['nit'], f))
    output = torch.tensor(x.reshape(n_samples, p), dtype=torch.float32)
    if store:
        storage, timer = callback.get_output()
        return output, storage, timer
    return output

max_stars_repo_name: pierreablin/ksddescent

max_issues_repo_name: pierreablin/ksddescent

max_forks_repo_name: pierreablin/ksddescent

reason: import numpy,from scipy


text:
# This code was developed by Man Yi Yim (manyi.yim@gmail.com) under Python 2.7.13.
# Note that the figure order, label and the content may have changed during revision.

import numpy as np
import pylab
import matplotlib as mpl
import pickle
import itertools
import os
figpath = './figure/'
if not os.path.exists(figpath):
    os.makedirs(figpath)
datapath = './data/'
if not os.path.exists(datapath):
    os.makedirs(datapath)
exec(open('gridplacefunc.py').read())
exec(open('mlfunc.py').read())
exec(open('mathfunc.py').read())

font_size = 14
mpl.rcParams['axes.titlesize'] = font_size
mpl.rcParams['xtick.labelsize'] = font_size-2
mpl.rcParams['ytick.labelsize'] = font_size-2
mpl.rcParams['axes.labelsize'] = font_size-2
mpl.rcParams['legend.fontsize'] = font_size-7
mpl.rcParams['font.size'] = font_size-1
new_rc_params = {'text.usetex': False,"svg.fonttype": 'none'}
mpl.rcParams.update(new_rc_params)
fs = 16

### Figures
color = ['#1e90ff','#ff8c00','#3cb371']
color4 = ['#DC143C','#3cb371','#9932CC','#FFD700','#1e90ff','#ff8c00','#3cb371','m']

def fig1_2():
    rng = np.random.RandomState(1)
    fig = pylab.figure(figsize=[7,5])
    fig.text(0.02,0.95,'A',fontsize=fs)
    fig.text(0.33,0.95,'B',fontsize=fs)
    fig.text(0.6,0.95,'C',fontsize=fs)
    fig.text(0.02,0.48,'D',fontsize=fs)
    fig.text(0.33,0.48,'E',fontsize=fs)
    fig.text(0.6,0.48,'F',fontsize=fs)
    # A
    ax = pylab.subplot(231)
    #ax.text(1,5,'pattern set $x_j, j=1,\dots,S$')
    temp = np.arange(1.75,2.25,0.05)
    #for j in [0,1,3]:
        #pylab.plot(j+temp,4+0.1*rng.randn(temp.size),'m.',markersize=3)
    pylab.plot([2,3,5],3*[3],'o',c='0.8',ms=20)
    ax.set_title('random input')
    ax.text(2,3,'1      2   $\cdots$      $N$')
    ax.text(2,2,'w x',fontweight='bold')
    pylab.plot([4],[2],'o',c='0.5',mfc='w',ms=22)
    #ax.text(3,2,'$\sum$')
    ax.text(2,1,'$f(x_j)$ = sgn$(w\cdot x_j)$')
    ax.text(3,0.5,'or')
    #ax.text(1,0,r'Learning rule: $\Delta w = \alpha(y_j-f(x_j))x_j$')
    pylab.xlim(1,7)
    pylab.ylim(0,5)
    ax.axis('off')
    # B
    ax = pylab.subplot(232)
    ax.text(0,4,'realizable',fontsize=14)
    ax.text(0,3,'unrealizable',fontsize=14)
    ax.text(0,1,'$+$',color='g',fontsize=15)
    ax.text(0,2,'$-$',color='r',fontsize=15)
    ax.text(1,1,'Field',color='g',fontsize=10)
    ax.text(1,2,'~Field',color='r',fontsize=10)
    ax.set_xlim(-1,4)
    ax.set_ylim(0,5)
    ax.axis('off')
    # C
    ax = pylab.subplot(233)
    N = 5
    R = 3*N
    pcover = []
    for x in range(1,R+1):
        pcover.append(cover(N,x)/float(2**x))
    pylab.plot(np.arange(1,R+1,1)/float(N),pcover,'0.8',label='$N$='+str(N))
    N = 50
    R = 3*N
    pcover = []
    for x in range(1,R+1):
        pcover.append(cover(N,x)/float(2**x))
    pylab.plot(np.arange(1,R+1,1)/float(N),pcover,'0.6',label='$N$='+str(N))
    N = 500
    R = 3*N
    pcover = []
    xrange = range(N/10,R+1,N/10)
    xrange.extend(range(2*N-100,2*N+100))
    xrange.sort()
    for x in xrange:
        if x <= 1024:
            pcover.append(float(cover(N,x)/2)/2**(x-1))
        else:
            pcover.append(1./(2**x/cover(N,x)))
    pylab.plot(np.array(xrange)/float(N),pcover,'k',label='$N$='+str(N))
    #pylab.plot([0,float(R)/N],[0,0],'k--',lw=1)
    pylab.plot([0,float(R)/N],[0.5,0.5],'k--',lw=1)
    #pylab.plot([0,float(R)/N],[1,1],'k--',lw=1)
    pylab.plot([1]*2,[0,1],'k--',lw=1)
    pylab.plot([2]*2,[0,1],'k--',lw=1)
    pylab.xticks(np.array([1,N,2*N])/float(N),('1','$N$','2$N$'))
    pylab.yticks(np.arange(0,1.1,0.5))
    pylab.ylim(0,1.05)
    pylab.xlabel('number of patterns $P$')
    pylab.ylabel('realizable fraction')
    pylab.legend(bbox_to_anchor=(0.67,0.6),frameon=False)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    # D
    ax = pylab.subplot(234)
    R = 6
    x = np.arange(-0.5,R-0.5+0.01,0.01)
    temp = grid(x,2,0.5)+grid(x,3,1./3)-6
    pylab.plot(x,grid(x,2,0.5)-0.5,c=color[0])
    pylab.plot(x,grid(x,3,1./3)-2.5,c=color[1])
    pylab.plot(x,temp,'r')
    pylab.plot(x[111:189],temp[111:189],'g')
    pylab.plot([x[0],x[-1]],-4.5*np.ones(2),'k--')
    #pylab.plot(x,grid(x,2,0)-7.5,c=color[0])
    #pylab.plot(x,grid(x,3,0)-9.5,c=color[1])
    #pylab.plot(x,grid(x,3,2./3)-11.5,c=color[1])
    #pylab.plot(x,grid(x,2,0)+0.5+grid(x,3,0)+0.5*grid(x,3,2./3)-14,'k')
    pylab.ylim(-7,1)
    pylab.xlim(-2,6)
    ax.set_title('place cell output')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['left'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    pylab.xticks([])
    pylab.yticks([])
    ax.axis('off')
    ax = pylab.subplot(235)
    temp = np.arange(1.75-0.5/12,2.25-0.5/12,0.01/12.)
    pylab.plot(4+0.3*grid(temp-1.75,2./12,0),temp,c=color[0])
    pylab.plot(4+0.3*grid(temp-1.75,2./12,0.5),temp+1,c=color[0])
    pylab.plot(4+0.3*grid(temp-1.75,3./12,2./3),temp+3,c=color[1])
    pylab.ylim(1,7)
    pylab.xlim(0,5)
    ax.set_title('grid cell input')
    ax.text(1,4,'place cell')
    ax.text(1,2.5,'perceptron')
    ax.axis('off')
    # F
    ax = pylab.subplot(236)
    for j in range(5):
        pylab.plot([0,R],np.ones(2)-j*1.5,c=color[int(j>1)],lw=1)
        pylab.plot([0,R],np.zeros(2)-j*1.5,c=color[int(j>1)],lw=1)
        for k in range(R+1):
            pylab.plot([k]*2,np.array([0,1])-j*1.5,c=color[int(j>1)],lw=1)
    mat = [['$x_1$','$x_2$','$x_3$','$x_4$','$x_5$','$x_6$'],['1','0','1','0','1','0'],['0','1','0','1','0','1'],['1','0','0','1','0','0'],['0','1','0','0','1','0'],['0','0','1','0','0','1']]
    for j in range(6):
        for k in range(R):
            ax.text(k+0.32,1.7-j*1.5,mat[j][k],fontsize=12)
    ax.text(0,1,'x',fontweight='bold',fontsize=12)
    ax.text(0,3,'$x_P$',fontsize=12)
    pylab.ylim(-7,2)
    pylab.xlim(-0.1,R+2)
    ax.axis('off')
    pylab.subplots_adjust(left=0.05,top=0.95,right=0.95,bottom=0.1,hspace=0.5,wspace=0.3)
    pylab.savefig(figpath+'f12.svg')

def fig6():
    pylab.figure(figsize=[7,7])
    ax = pylab.subplot(222)
    seed = 5
    rng = np.random.RandomState(seed)
    wp1 = 0
    wp2 = 1
    bin = 1
    l = [31,43,59]
    Ng = 15
    Np = 100
    sig = 0.16
    pid = 0
    R = l[0]
    for j in range(len(l)-1):
        R *= l[j+1]
    x = np.arange(R)
    w = rng.lognormal(wp1,wp2,size=[Np,Ng])
    p1d = rng.rand(Np,Ng)
    v = np.zeros(((Np,Ng,R)))
    a = np.zeros((Np,R))
    print '1D'
    for iN in range(Np):
        for iM in range(Ng):
            v[iN,iM,:] = grid(x,l[iM/(Ng/len(l))],p1d[iN,iM],phsh=0.,gridmodel='gau',sig=sig)
        a[iN,:] = np.dot(w[iN,:],v[iN,:,:])
    # 1D
    nth = 2
    af,tharr = detect_field(a,nth,fmerge=10,fwidthmin=0)
    fid1d = []
    nf = np.zeros(Np)
    for iN in range(Np):
        fid1d.append(list(pylab.find(af[iN,:]==1)))
        nf[iN] = np.sum(af[iN])
    print nf
    ifiall = []
    for iN in range(Np):   # choose one neuron for the IFI histogram
        ifiall.extend(np.diff(fid1d[iN]))
    y = np.histogram(ifiall,np.arange(0.5,np.max(ifiall)+1,bin))
    if 1:
        pylab.figure()
        for iNg in range(Ng):
            pylab.plot(v[0,iNg,:]+iNg+1,'k')
        pylab.plot(a[0,:]/np.max(a[0,:])-1,'r')
        for k in np.where(af[0]==1)[0]:
           pylab.plot([k]*2,[-1,16],'b--')
        pylab.xlim(0,1500)
        pylab.yticks([])
        pylab.xlabel('location (up to 78647 cm)')
        pylab.ylabel('15 GCs in black; sum in red, rescaled')
        pylab.title('blue lines denote place fields')
        pylab.subplots_adjust(left=0.05,right=0.95,bottom=0.15)
        pylab.savefig(figpath+'fig6A.svg')
    count = np.zeros(Ng+1)
    for iNp in range(Np):
        for k in np.where(af[iNp]==1)[0]:
            idx = 0
            for iNg in range(Ng):
                gpeaks = pylab.find(v[iNp,iNg,:]>0.95)
                idx += np.any(k==gpeaks)
            count[idx] += 1
    pylab.bar(range(Ng+1),count/np.sum(count),color='k')
    pylab.yticks(np.arange(0.,0.4,0.1))
    pylab.ylabel('fraction')
    pylab.xlabel('number of coincidence peaks')
    pylab.subplots_adjust(left=0.05,right=0.95,hspace=0.4)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    # IFI
    ax = pylab.subplot(223)
    for j in range(len(l)):
        for k in range(1,10):
            pylab.plot([k*l[j]]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    pylab.hist(ifiall,np.arange(0.5,np.max(ifiall)+1,bin),color='k')
    acf1 = np.correlate(v[pid,0,:],v[pid,0,:],'same')
    acf2 = np.correlate(v[pid,Ng/3,:],v[pid,Ng/3,:],'same')
    acf3 = np.correlate(v[pid,2*Ng/3,:],v[pid,2*Ng/3,:],'same')
    acf1 /= np.max(acf1)
    acf2 /= np.max(acf2)
    acf3 /= np.max(acf3)
    pylab.plot(range(-R/2,R/2),np.max(ifiall)*(acf1+acf2+acf3)/3.,'k',label='ACF')
    pylab.xlim(0,300)
    pylab.xlabel('IFI')
    pylab.title('1D periodic')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    # slices through 2D
    ax = pylab.subplot(224)
    l = [31,43]
    Np = 100
    Ng = Ng*2/3
    px = rng.rand(Np,Ng)
    py = rng.rand(Np,Ng)
    w = rng.lognormal(wp1,wp2,size=[Np,Ng])
    ori = [0.39,0.46]
    R = 20000
    x = np.arange(R)
    v = np.zeros(((Np,Ng,R)))
    a = np.zeros((Np,R))
    for iN in range(Np):
        for iM in range(Ng):
            v[iN,iM,:] = grid1d_orient(x,l[iM/(Ng/len(l))],ori[iM/(Ng/len(l))],px[iN,iM],py[iN,iM],sig=sig)
        a[iN,:] = np.dot(w[iN,:],v[iN,:,:])
    print '2D'
    af,tharr = detect_field(a,nth,fmerge=10,fwidthmin=0)
    print np.sum(af[0])
    fid = []
    for iN in range(Np):
        fid.append(list(pylab.find(af[iN,:]==1)))
    ifiall = []
    for iN in range(Np):
        ifiall.extend(np.diff(fid[iN]))
    y = np.histogram(ifiall,np.arange(0.5,np.max(ifiall)+1,bin))
    ymax = max(y[0])
    j = 0
    pylab.plot([l[j]*np.sqrt(3)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    #pylab.plot([l[j]*np.sqrt(3)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(7)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]+2*l[j]*np.sqrt(3)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    #pylab.plot([2*l[j]+l[j]*np.sqrt(3)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(3)+2*l[j]*np.sqrt(7)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    j = 1
    #pylab.plot([l[j]]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(3)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    pylab.plot([2*l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    #pylab.plot([l[j]*np.sqrt(7)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    #pylab.plot([l[j]+2*l[j]*np.sqrt(3)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    #pylab.plot([l[j]+2*l[j]*np.sqrt(7)]*2,[0,ymax+1],'--',c=color[j],lw=0.5)
    pylab.hist(ifiall,np.arange(0.5,np.max(ifiall)+1,bin),color='k')
    acf1 = np.correlate(v[pid,0,:],v[pid,0,:],'same')
    acf2 = np.correlate(v[pid,Ng/2,:],v[pid,Ng/2,:],'same')
    acf1 /= np.max(acf1)
    acf2 /= np.max(acf2)
    pylab.plot(range(-R/2,R/2),np.max(ifiall)*(acf1+acf2)/2.,'k',label='ACF')
    pylab.xlim(0,300)
    #pylab.yticks([0,5,10])
    pylab.xlabel('IFI')
    pylab.title('1D slices through 2D grid')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    pylab.subplots_adjust(top=0.9,wspace=0.3,hspace=0.4)
    pylab.savefig(figpath+'f6_Ng_'+str(Ng/2)+'_seed'+str(seed)+'.svg')

def fig6b():
    Ng = 10  # dinstinct for each neuron
    Np = 5
    N2plot = 5
    l = [31,43] #[53,87]
    R = 10000
    sig = 0.16
    bin = 1
    nth = 2
    #ori = [0.2153,0.3162] # ori2
    ori = [0.39,0.46] # ori1
    rng = np.random.RandomState(52)
    x = np.arange(R)
    px = rng.rand(Np,Ng)
    py = rng.rand(Np,Ng)
    w = rng.lognormal(0,0.25,size=[Np,Ng])
    v = np.zeros(((Np,Ng,R)))
    a = np.zeros((Np,R))
    acf = np.zeros((Np,2*R/bin-1))
    for iN in range(Np):
        for iM in range(Ng):
            v[iN,iM,:] = grid1d_orient(x,l[iM/(Ng/len(l))],ori[iM/(Ng/len(l))],px[iN,iM],py[iN,iM],sig=sig)
            #v[iN,iM,:] = grid(x,l[iN],float(iM)/l[iN],phsh=0.,gridmodel='del'')    #
            a[iN,:] += w[iN,iM]*v[iN,iM,:]
    af,tharr = detect_field(a,nth,fmerge=10,fwidthmin=0)
    fid = []
    for iN in range(Np):
        fid.append(list(pylab.find(af[iN,:]==1)))
    # Activity
    fig = pylab.figure(figsize=[12,7])
    for iN in range(N2plot):
        ax = pylab.subplot(N2plot,1,iN+1)
        pylab.plot(x,a[iN,:])
        pylab.plot(fid[iN],[tharr[iN]]*len(fid[iN]),'ko')
        pylab.plot([0,R],[tharr[iN]]*2,'k',lw=1)
        if iN != Np-1:
            pylab.xticks([])
        if iN == 0:
            pylab.title('Activity on track: $\lambda$='+str(l)+';#grid='+str(Ng)+';R='+str(R)+';#STD='+str(nth))
    pylab.savefig(figpath+'slices_act_l'+str(l)+'g'+str(Ng)+'R'+str(R)+'nstd'+str(nth)+'.png')
    # IFI
    fig = pylab.figure(figsize=[12,7])
    ifiall = []
    for iN in range(Np):
        ifiall.extend(np.diff(fid[iN]))
    for iN in range(N2plot):
        ax = pylab.subplot(2,3,iN+1)
        ifi = np.diff(fid[iN])
        pylab.hist(ifi,np.arange(0.5,np.max(ifiall)+1,bin))
        if iN == 1:
            pylab.title('IFI: 5 neurons & all; bin='+str(bin)+';$\lambda$='+str(l)+';#grid='+str(Ng)+';R='+str(R)+';#STD='+str(nth))
        pylab.xlim(0,200)
    ax = pylab.subplot(236)
    y = np.histogram(ifiall,np.arange(0.5,np.max(ifiall)+1,bin))

    j = 0
    pylab.plot([l[j]*np.sqrt(3)]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(7)]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(3)+2*l[j]*np.sqrt(7)]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    j = 1
    pylab.plot([l[j]*np.sqrt(3)]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    pylab.plot([l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    pylab.plot([2*l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,max(y[0])+1],'--',c=color[j],lw=0.5)
    pylab.hist(ifiall,np.arange(0.5,np.max(ifiall)+1,bin))
    pylab.xlim(0,200)
    pylab.savefig(figpath+'slices_ifi_l'+str(l)+'g'+str(Ng)+'R'+str(R)+'b'+str(bin)+'nstd'+str(nth)+'.png')
    # ACF
    x = np.arange(0,R,bin)
    fig = pylab.figure(figsize=[12,7])
    for iN in range(Np):
        afbin = []
        for j in range(R/bin):
            afbin.append(np.sum([af[iN,bin*j:bin*(j+1)]]))
        acf[iN,:] = np.correlate(afbin,afbin,'full')
        if iN < N2plot:
            ax = pylab.subplot(2,3,iN+1)
            pylab.plot(x[1:],acf[iN,R/bin:])
            if iN == 1:
                pylab.title('ACF: 5 neurons & all; bin='+str(bin)+';$\lambda$='+str(l)+';#grid='+str(Ng)+';R='+str(R)+';#STD='+str(nth))
            pylab.xlim(0,200)
    ax = pylab.subplot(236)
    acfmax = np.max(np.sum(acf[:,R/bin:],0))
    for j in range(len(l)):
        for k in range(1,10):
            pylab.plot([k*l[j]]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
            pylab.plot([k*l[j]*np.sqrt(3)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
            pylab.plot([k*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        # 2
        pylab.plot([l[j]+l[j]*np.sqrt(3)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        # 3
        pylab.plot([2*l[j]+l[j]*np.sqrt(3)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([2*l[j]+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([2*l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+2*l[j]*np.sqrt(3)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+2*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]*np.sqrt(3)+2*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        # 4
        pylab.plot([3*l[j]+l[j]*np.sqrt(3)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([3*l[j]+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([3*l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+3*l[j]*np.sqrt(3)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+3*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]*np.sqrt(3)+3*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([2*l[j]+2*l[j]*np.sqrt(3)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([2*l[j]+2*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([2*l[j]*np.sqrt(3)+2*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([2*l[j]+l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+2*l[j]*np.sqrt(3)+l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
        pylab.plot([l[j]+l[j]*np.sqrt(3)+2*l[j]*np.sqrt(7)]*2,[0,acfmax+1],'--',c=color[j],lw=0.5)
    pylab.plot(x[1:],np.sum(acf[:,R/bin:],0))
    pylab.xlim(0,200)
    pylab.savefig(figpath+'slices_acf_l'+str(l)+'g'+str(Ng)+'R'+str(R)+'b'+str(bin)+'nstd'+str(nth)+'.png')
    # PSD
    bin = 1
    acf = np.zeros((Np,2*R/bin-1))
    x = np.arange(2*R/bin)/(2.*R/bin)
    for iN in range(Np):
        afbin = []
        for j in range(R/bin):
            afbin.append(np.sum([af[iN,bin*j:bin*(j+1)]]))
        acf[iN,:] = np.correlate(afbin,afbin,'full')
    fig = pylab.figure(figsize=[12,7])
    for iN in range(N2plot):
        ax = pylab.subplot(2,3,iN+1)
        psd = np.abs(np.fft.fft(np.append([0],acf[iN,:])))
        psdmax = np.max(psd[1:])
        for j in range(len(l)):
            pylab.plot([2.*np.sin(ori[j])/(np.sqrt(3)*l[j])]*2,[0,psdmax+10],'--',c=color[j],lw=1)
            pylab.plot([1./l[j]*(np.cos(ori[j])-np.sin(ori[j])/np.sqrt(3))]*2,[0,psdmax+10],'--',c=color[j],lw=1)
            pylab.plot([1./l[j]*(np.cos(ori[j])+np.sin(ori[j])/np.sqrt(3))]*2,[0,psdmax+10],'--',c=color[j],lw=1)
            pylab.plot([2*2.*np.sin(ori[j])/(np.sqrt(3)*l[j])]*2,[0,psdmax+10],'--',c=color[j],lw=1)
            pylab.plot([2*1./l[j]*(np.cos(ori[j])-np.sin(ori[j])/np.sqrt(3))]*2,[0,psdmax+10],'--',c=color[j],lw=1)
            pylab.plot([2*1./l[j]*(np.cos(ori[j])+np.sin(ori[j])/np.sqrt(3))]*2,[0,psdmax+10],'--',c=color[j],lw=1)
        pylab.plot(x[1:],psd[1:])
        if iN == 1:
            pylab.title('PSD: 5 neurons & all; bin='+str(bin)+';$\lambda$='+str(l)+';#grid='+str(Ng)+';R='+str(R)+';#STD='+str(nth))
        pylab.xlim(0,0.05)
    ax = pylab.subplot(236)
    psd = np.abs(np.fft.fft(np.append([0],np.sum(acf,0))))
    psdmax = np.max(psd[1:])
    for j in range(len(l)):
        for k in range(1,10):
            pylab.plot([k*2.*np.sin(ori[j])/(np.sqrt(3)*l[j])]*2,[0,psdmax+10],'--',c=color[j],lw=1)
            pylab.plot([k*1./l[j]*(np.cos(ori[j])-np.sin(ori[j])/np.sqrt(3))]*2,[0,psdmax+10],'--',c=color[j],lw=1)
            pylab.plot([k*1./l[j]*(np.cos(ori[j])+np.sin(ori[j])/np.sqrt(3))]*2,[0,psdmax+10],'--',c=color[j],lw=1)
    pylab.plot(x[1:],psd[1:])
    pylab.xlim(0,0.05)
    pylab.savefig(figpath+'slices_psd_l'+str(l)+'g'+str(Ng)+'R'+str(R)+'nstd'+str(nth)+'.png')

def fig3():
    fig = pylab.figure(figsize=[7,7])
    fig.text(0.05,0.96,'A',fontsize=fs)
    fig.text(0.05,0.65,'B',fontsize=fs)
    fig.text(0.05,0.35,'C',fontsize=fs)
    # A
    ax = pylab.subplot(711)
    pylab.plot([-8,-8],[1,3],c=color[0],lw=1)
    pylab.plot([-8,-8],[1,3],'ko',ms=15,mfc='w')
    pylab.plot([-5.5,-3.5,-4.5,-5.5],[2.5,2.5,1.5,2.5],c=color[1],lw=1)
    pylab.plot([-5.5,-3.5,-4.5,-5.5],[2.5,2.5,1.5,2.5],'ko',ms=15,mfc='w')
    pylab.plot([-1,-1],[1,3],c=color[0],lw=1)
    pylab.plot([1,1],[1,3],c=color[0],lw=1)
    pylab.plot([0,0],[0,2],c=color[0],lw=1)
    pylab.plot([-1,1,0,-1],[1,1,0,1],c=color[1],lw=1)
    pylab.plot([-1,1,0,-1],[3,3,2,3],c=color[1],lw=1)
    pylab.plot([-1,1,0,-1],[1,1,0,1],'ko',ms=15,mfc='w')
    pylab.plot([-1,1,0,-1],[3,3,2,3],'ko',ms=15,mfc='w')
    pylab.text(-8.5,3.5,'(1,0)')
    pylab.text(-8.5,0,'(0,1)')
    pylab.text(-6,3,'(1,0,0)')
    pylab.text(-4,3,'(0,1,0)')
    pylab.text(-5,0,'(0,0,1)')
    pylab.text(-2.5,3.5,'(1,0,1,0,0)')
    pylab.text(0.5,3.5,'(1,0,0,1,0)')
    pylab.text(0.5,1.5,'(1,0,0,0,1)')
    pylab.text(-2.5,0,'(0,1,1,0,0)')
    pylab.text(0.5,0,'(0,1,0,1,0)')
    pylab.text(0,-1,'(0,1,0,0,1)')
    pylab.text(-7,2,'x',fontsize=20)
    pylab.text(-2.5,2,'=',fontsize=20)
    pylab.xlim(-9,3.5)
    pylab.ylim(-0.5,4)
    ax.axis('off')
    # B
    ax = pylab.subplot(745)
    pylab.plot([-1,-1],[1,3],c=color[0],lw=1)
    pylab.plot([1,1],[1,3],c=color[0],lw=1)
    pylab.plot([0,0],[0,2],c=color[0],lw=1)
    pylab.plot([-1,1,0,-1],[1,1,0,1],c=color[1],lw=1)
    pylab.plot([-1,1,0,-1],[3,3,2,3],c=color[1],lw=1)
    pylab.plot([-1,1,0,-1],[1,1,0,1],'ko',ms=10,mfc='w')
    pylab.plot([-1,1,0,-1],[3,3,2,3],'ko',ms=10,mfc='w')
    #pylab.plot([-1],[3],'go',ms=15,alpha=0.5)
    pylab.xlim(-2,2)
    pylab.ylim(-1,4.5)
    pylab.text(-1.45,-1,'$------$',color='r',fontsize=11,fontweight='bold')
    pylab.text(-1.9,-1,'$+$',color='g',fontsize=11,fontweight='bold')
    ax.axis('off')
    ax = pylab.subplot(713)
    pylab.text(0,5,'Remove 1st\n location',fontsize=10)
    pylab.text(0,2,'Remove 2nd\n location \n (adjacent to 1st)',fontsize=10)
    pylab.xlim(-1,6)
    pylab.ylim(-1,6)
    ax.axis('off')
    ax = pylab.subplot(223)
    R = 6
    x = np.arange(-0.5,R-0.5+0.01,0.01)
    temp = grid(x,2,0)+grid(x,2,0.5)+grid(x,3,0)+grid(x,3,2./3)-11
    pylab.plot(x,grid(x,2,0)-0.5,c=color[0])
    pylab.plot(x,grid(x,2,0.5)-2.5,c=color[0])
    pylab.plot(x,grid(x,3,0)-4.5,c=color[1])
    pylab.plot(x,grid(x,3,2./3)-6.5,c=color[1])
    pylab.plot(x,temp,'r')
    fw = 59
    pylab.plot(x[:fw],temp[:fw],'g')
    pylab.plot(x[-fw:],temp[-fw:],'g')
    pylab.plot(x[300-fw:300+fw],temp[300-fw:300+fw],'g')
    pylab.plot([x[0],x[-1]],-8.7*np.ones(2),'k--')
    #pylab.plot(x,grid(x,2,0)-7.5,c=color[0])
    #pylab.plot(x,grid(x,3,0)-9.5,c=color[1])
    #pylab.plot(x,grid(x,3,2./3)-11.5,c=color[1])
    #pylab.plot(x,grid(x,2,0)+0.5+grid(x,3,0)+0.5*grid(x,3,2./3)-14,'k')
    pylab.ylim(-10,1)
    pylab.xlim(-2,6)
    #ax.set_title('place cell output')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['left'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    pylab.xticks([])
    pylab.yticks([])
    ax.axis('off')
    pylab.subplots_adjust(left=0.1,top=0.97,right=0.95,bottom=0.1,wspace=0.2)
    pylab.savefig(figpath+'f3.svg')

def fig5():
    mpl.rcParams['legend.fontsize'] = font_size-5
    mpl.rcParams['axes.titlesize'] = font_size-3
    rng = np.random.RandomState(3)
    fig = pylab.figure(figsize=[7,7])
    fig.text(0.02,0.95,'A',fontsize=fs)
    fig.text(0.5,0.95,'B',fontsize=fs)
    fig.text(0.02,0.64,'C',fontsize=fs)
    fig.text(0.5,0.64,'D',fontsize=fs)
    fig.text(0.02,0.32,'E',fontsize=fs)
    fig.text(0.5,0.32,'F',fontsize=fs)
    # A
    l = [3,4]
    N = len(l)
    Sc = int(np.round(testrange(l)[0]))
    R = l[0]
    for j in range(N-1):
        R = lcm(R,l[j+1])
    fname = 'pcorr'+str(l[0])+str(l[1])
    if os.path.isfile(datapath+fname+'.txt'):
        with open(datapath+fname+'.txt','rb') as f:
            pall,p2,p3,p4,p5,p6 = pickle.load(f)
    else:
        pall,p2,p3,p4,p5,p6 = frac_vs_S(l,R,return6=1)
        with open(datapath+fname+'.txt','wb') as f:
            pickle.dump((pall,p2,p3,p4,p5,p6),f)
    ax = pylab.subplot(321)
    pylab.plot(range(1,R+1),np.ones(R),'o-',ms=5,label='1')
    pylab.plot(range(2,R+1),p2,'o-',ms=5,label='2')
    pylab.plot(range(3,R+1),p3,'o-',ms=5,label='3')
    pylab.plot(range(4,R+1),p4,'o-',ms=5,label='4')
    pylab.plot(range(5,R+1),p5,'o-',ms=5,label='5')
    pylab.plot(range(6,R+1),p6,'o-',ms=5,label='6')
    pylab.plot([R]*2,[0,1],'k--',lw=1)
    pylab.text(11,0.85,'$L$')
    pylab.ylim(0,1.05)
    pylab.xticks([1,Sc,R],('1',str(Sc),str(R)))
    pylab.yticks([0,1])
    pylab.legend(loc=3,frameon=False)
    pylab.ylabel('realizable fraction')
    pylab.xlabel('$l$')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.set_title('$\lambda=\{3,4\}$')
    # B
    ax = pylab.subplot(322)
    mpl.rcParams['legend.fontsize'] = font_size-7
    pylab.plot(range(1,R+1),pall,'ko-',ms=5,label='grid')
    pylab.plot([Sc]*2,[0,1],'k--',lw=1)
    pylab.plot([R]*2,[0,1],'k--',lw=1)
    pcover = []
    for x in range(1,R+1):
        pcover.append(cover(Sc,x)/float(2**x))
    pylab.plot(np.arange(1,R+1,1),pcover,'o-',c='c',ms=5,label='random (rank)')
    pcover = []
    for x in range(1,R+1):
        pcover.append(cover(np.sum(l),x)/float(2**x))
    pylab.plot(np.arange(1,R+1,1),pcover,'o-',c='#8A2BE2',ms=5,label='random (dim)')
    for k in range(5):
        pwcon = random_weightcon(np.sum(l),R,1,k+1)
        if k == 0:
            pylab.plot(np.arange(1,R+1,1),pwcon,'o-',c='#B0C4DE',ms=5,label='random,\nconstrained')
        else:
            pylab.plot(np.arange(1,R+1,1),pwcon,'o-',c='#B0C4DE',ms=5)
    pylab.legend(loc=6,frameon=False)
    pylab.plot(range(1,R+1),pall,'ko-',ms=5,label='grid')
    pylab.ylim(0,1.05)
    pylab.xticks([1,Sc,R],('1',str(Sc),str(R)))
    pylab.yticks([0,1])
    pylab.text(6.2,0.05,'$l^*$')
    pylab.text(11,0.05,'$L$')
    pylab.xlabel('$l$')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    # C
    mpl.rcParams['legend.fontsize'] = font_size-5
    ax = pylab.subplot(323)
    lsum = []
    ctrackq0 = []
    ctrackq2 = []
    ctrackq4 = []
    # [lmin,lmax)
    lmin = 3
    lmax = 20
    for j in range(100):
        M = rng.randint(2,7)
        l = lmin + rng.rand(M)*(lmax-lmin)
        temp = testrange(l)
        lsum.append(np.sum(l))
        ctrackq0.append(temp[0])
        ctrackq2.append(temp[2])
        ctrackq4.append(temp[4])
        #ctrackq8.append(temp[-1])
        pylab.plot(np.array(temp)/np.sum(l),'-',lw=1)
    pylab.xlabel('resolution $q$')
    pylab.ylabel('$R^q_{re}/\Sigma$')
    pylab.ylim(0,1.05)
    pylab.yticks([0,1])
    pylab.xticks(range(0,9,2))
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    # F
    ax = pylab.subplot(324)
    pylab.plot(lsum,ctrackq0,'.',lw=0,ms=5,label='$q$=0')
    pylab.plot(lsum,ctrackq2,'.',lw=0,ms=5,label='$q$=2')
    pylab.plot(lsum,ctrackq4,'.',lw=0,ms=5,label='$q$=4')
    #pylab.plot(lsum,ctrackq8,'.',lw=0,ms=5,label='$S_{real}^{q=8}$')
    pylab.legend(loc=2,frameon=False)
    pylab.plot([1,120],[1,120],'k--',lw=1)
    pylab.xlim(0,110.5)
    pylab.ylim(0,110.5)
    pylab.xticks(range(0,101,50))
    pylab.yticks(range(0,101,50))
    pylab.xlabel('period sum $\Sigma$')
    pylab.ylabel('$R^q_{re}$')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    pylab.subplots_adjust(left=0.15,top=0.95,right=0.95,bottom=0.1,wspace=0.3,hspace=0.4)
    pylab.savefig(figpath+'f5.svg')

def fig8():
    is_qp = 0   # use quadratic programming with the options of constrained weights
    if is_qp:
        print 'Using quadratic programming'
    else:
        print 'Using sklearn SVM'
    mpl.rcParams['legend.fontsize'] = font_size-6
    rng = np.random.RandomState(4)
    import seaborn as sns
    gridmodel = 'del'
    mth = np.arange(0,1.0001,0.01)
    sym = ['^','+','x']
    num = 10
    # A
    fig = pylab.figure(figsize=[7,8.5*0.8])
    fig.text(0.02,0.65,'A',fontsize=fs)
    fig.text(0.48,0.65,'B',fontsize=fs)
    fig.text(0.02,0.3,'C',fontsize=fs)
    fig.text(0.48,0.3,'D',fontsize=fs)
    ax = fig.add_subplot(323)
    l = [31,43]
    N = len(l)
    R = l[0]
    for j in range(N-1):
        R = lcm(R,l[j+1])
    K = 6
    u = act_mat_grid_binary(l)
    u /= 2
    mode = 's1000'
    if os.path.exists(datapath+'f8_s1000.txt'):
        with open(datapath+'f8_s1000.txt','rb') as f:
            margin,rmargin,smargin,numKarr,rnumKarr,snumKarr = pickle.load(f)
    else:
        margin,rmargin,smargin,numKarr,rnumKarr,snumKarr = margin_gridvsrandom(l,K=K,num=num,mode=mode)
    margin = margin[:K]
    rmargin = rmargin[:K]
    smargin = smargin[:K]
    numKarr = numKarr[:K]
    rnumKarr = rnumKarr[:K]
    snumKarr = snumKarr[:K]
    # for violin plot: random
    kmat1 = []
    for k in range(K):
        kmat1.extend([k+1]*np.sum(rnumKarr[k]))
    mmat1 = [item for sublist in rmargin for subsublist in sublist for item in subsublist]
    nmat1 = ['random']*np.sum(np.sum(rnumKarr))
    # for violin plot: shuffle
    kmat2 = []
    for k in range(K):
        kmat2.extend([k+1]*np.sum(snumKarr[k]))
    mmat2 = [item for sublist in smargin for subsublist in sublist for item in subsublist]
    nmat2 = ['shuffled']*np.sum(np.sum(snumKarr))
    kmat1 = np.array(kmat1)
    kmat2 = np.array(kmat2)
    #sns.violinplot(np.append(kmat1,kmat2),np.append(mmat1,mmat2),np.append(nmat1,nmat2),inner=None,linewidth=.4,scale='width',width=0.5,bw=.2,gridsize=100)
    sns.violinplot(kmat1,mmat1,inner=None,linewidth=.4,scale='width',width=0.5,bw=.2,gridsize=100,color='m')
    sns.violinplot(kmat2,mmat2,inner=None,linewidth=.4,scale='width',width=0.5,bw=.2,gridsize=100,color='b')
    #sns.violinplot(kmat1,mmat1,nmat1,inner=None,linewidth=.4,bw=.2)
    #sns.violinplot(kmat2,mmat2,nmat2,inner=None,linewidth=.4,bw=.2,color='#ff7f0e')
    for k in range(1,K+1):
        for mu in margin[k-1]:
            pylab.plot(k-1+np.array([-0.3,0.3]),2*[mu],'k')
    ax.set_yticks(np.arange(0,0.5,0.2))
    ax.set_xlim(-0.5,K-0.5)
    ax.set_ylim(0,0.4)
    ax.set_xlabel('number of fields ($K$)')
    ax.set_ylabel('margin')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    # B
    ax = fig.add_subplot(324)
    m_inset = []
    for k in range(1,K+1):
        count,temp = np.histogram(rng.rand(int(mode[1:])*num),np.array([0.]+list(np.cumsum(numKarr[k-1]))+[nCr(R,k)])/nCr(R,k)) # randomly draw field arrangements based on fractions of margins
        for j in range(len(count)-1):
            m_inset.extend([margin[k-1][j]]*count[j])
        print np.array([0.]+list(np.cumsum(numKarr[k-1]))+[nCr(R,k)]),count,m_inset
    count0L = []
    for m in mth:
        count0L.append(np.sum(np.array(m_inset)>=m))
    countr = []
    for m in mth:
        countr.append(np.sum(mmat1>=m))
    counts = []
    for m in mth:
        counts.append(np.sum(mmat2>=m))
    temp = float(K)*int(mode[1:])*num
    ax.plot(mth,np.array(count0L)/temp,'k')
    ax.plot(mth,np.array(countr)/temp,color='m',lw=1.5) #1f77b4
    ax.plot(mth,np.array(counts)/temp,color='b',lw=1.5)    #ff7f0e
    ax.plot([0,1],2*[1],'k--',lw=1)
    ax.set_ylim(0,1)
    ax.set_xlim(0,0.4)
    ax.set_xlabel('margin $\kappa$')
    ax.set_ylabel('CDF')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    # CD
    color = ['b','g','r','c','m']
    inp = [0,100]  # first entry is zero
    numr = len(inp)
    u = act_mat_grid_binary(l)
    if 1:
        for id in range(2):
            margin_spatial = []
            margin_new = []
            for n in range(num):
                for r in range(numr): # from no addition to numr additions
                    if n == 0:
                        margin_spatial.append([])
                        if r > 0:
                            margin_new.append([])
                    if n > 0 and r == 0:
                        continue
                    if inp[r] > 0:
                        if id == 0:
                            # uniform
                            randinp = 0.054*0.2*rng.rand(inp[r],u.shape[1])    #2*2/74 *0.2
                            # gaussian
                            #randinp = 0.1*rng.randn(inp[r],u.shape[1])
                        elif id == 1:
                            randinp = np.zeros((inp[r],u.shape[1]))
                            for j in range(inp[r]):
                                #randinp[j,rng.choice(range(u.shape[1]),10,replace=False)] = 1
                                randinp[j,rng.choice(range(u.shape[1]),7,replace=False)] = 1    # 1333*2/74.*0.2
                        v = np.append(u,randinp,axis=0)
                    else:
                        v = np.copy(u)
                    # Normalization
                    for jj in range(u.shape[1]):
                        # L1
                        v[:,jj] = v[:,jj]/np.sum(v[:,jj])
                        # L2
                        #v[:,jj] = v[:,jj]/pylab.norm(v[:,jj])
                    for k in range(1,K+1):
                        if n == 0:
                            margin_spatial[r].append([])
                            if r > 0:
                                margin_new[r-1].append([])
                        partition = partitions(k)
                        part = []
                        for p in partition:
                            if np.all(np.array(p)<=np.min(l)):
                                part.append(list(p))
                                # Young diagram
                                mat = np.zeros((l[0],l[1]))
                                for j in range(len(p)):
                                    mat[:p[j],j] = 1
                                #pylab.figure()
                                #pylab.imshow(mat,aspect='auto')
                                i1 = np.tile(range(l[0]),l[1])
                                i2 = np.tile(range(l[1]),l[0])
                                Y = mat[i1,i2]
                                print Y
                                m,w,b = svm_margin(v.T,Y)
                                margin_spatial[r][k-1].append(m)
                                dec = np.sign(np.dot(w.T,v)+b)
                                dec[dec<0] = 0
                                print k,p,r,abs(np.sum(np.abs(Y-dec))),m
                        is_fit = 0
                        while r>0 and k>1 and not is_fit:     # unrealizable becomes realizable
                            print n,r,k
                            Y = np.zeros(u.shape[1])
                            Y[0:k] = 1
                            rng.shuffle(Y)
                            m,w,b = svm_margin(u.T,Y)
                            dec = np.sign(np.dot(w.T,u)+b)
                            dec[dec<0] = 0
                            if abs(np.sum(np.abs(Y-dec)))>1e-10: # not realizable
                                print 'found not realizable with u'
                                m,w,b = svm_margin(v.T,Y)
                                dec = np.sign(np.dot(w.T,v)+b)
                                dec[dec<0] = 0
                                if abs(np.sum(np.abs(Y-dec)))<1e-10:    # realizable
                                    print 'found realizable with v'
                                    margin_new[r-1][k-1].append(m)
                                    is_fit = 1
            with open(datapath+'fig8CD_'+str(id)+'.txt','wb') as f:
                pickle.dump((margin_spatial,margin_new),f)
    color = ['b','#0f9b8e','#0cfc73']
    for id in range(2):
        ax = pylab.subplot(3,2,5+id)
        with open(datapath+'fig8CD_'+str(id)+'.txt','r') as f:
            margin_spatial,margin_new = pickle.load(f)
        # violin
        if 1:
            mmat1 = [item for sublist in margin_spatial[numr-1][:K] for item in sublist]
            mmat2 = [item for sublist in margin_new[numr-2][:K] for item in sublist]
            kmat1 = []
            kmat2 = []
            for k in range(1,K+1):
                kmat1.extend([k]*len(margin_spatial[numr-1][k-1]))
                kmat2.extend([k]*len(margin_new[numr-2][k-1]))
            nmat1 = ['exiting']*len(kmat1)
            nmat2 = ['new']*len(kmat2)
            sns.violinplot(kmat1,mmat1,inner=None,linewidth=.4,scale='width',width=0.5,bw=.2,gridsize=100,color=color4[0])
            sns.violinplot(np.append([1],kmat2),np.append([-1],mmat2),inner=None,linewidth=.4,scale='width',width=0.5,bw=.2,gridsize=100,color='#ff7f0e')
        for k in range(1,K+1):
            for r in [0,numr-1]:#range(numr):
                if r > 0:
                    # mean
                    if 0:
                        pn = len(margin_spatial[r][k-1])/num
                        for j in range(pn):
                            pylab.plot([k],[np.mean(margin_spatial[r][k-1][j::pn])],'.',ms=10,fillstyle='none',c=color[r-1],label=str(inp[r]))
                        pylab.plot([k],[np.mean(margin_new[r-1][k-1])],'^',ms=6,fillstyle='none',c=color[r-1])
                    # all trials
                    if 0:
                        pylab.plot([k]*len(margin_spatial[r][k-1]),margin_spatial[r][k-1],'x',c=color[r-1],label=str(inp[r])+' random input')
                else:
                    for mu in margin_spatial[0][k-1]:
                        pylab.plot(k-1+np.array([-0.3,0.3]),2*[mu],'k',label='grid')
            #if k == 1 and id == 0:
                #pylab.legend(loc=1,frameon=False)
        pylab.xlim(-0.5,K-0.5)
        pylab.ylim(0,0.45)
        pylab.yticks([0,0.2,0.4])
        #pylab.xticks(range(1,K+1))
        if id == 0:
            pylab.ylabel('margin')
        pylab.xlabel('# fields ($K$)')
        ax.spines['right'].set_visible(False)
        ax.spines['top'].set_visible(False)
    fig.subplots_adjust(left=0.1,top=0.95,right=0.95,bottom=0.1,wspace=0.4,hspace=0.5)
    fig.savefig(figpath+'f8'+'qp'*is_qp+'.svg')

def fig7(seed):
    rng = np.random.RandomState(seed)
    wp1 = 0
    wp2 = 1
    bin = 1
    #l = [31,43,59]
    l = [31,43]
    ori = [0.39,0.46]
    Ng = 600
    Np = 10#0
    sig = 0.16
    R = l[0]
    for j in range(len(l)-1):
        R *= l[j+1]
    R = 2000
    nK = 4
    nadd = 4
    nex = 10
    x = np.arange(R)
    w = rng.lognormal(wp1,wp2,size=[Np,Ng])
    for iN in range(Np):
        if np.sum(w[iN,:]) > 0:
            w[iN,:] /= sum(w[iN,:])
    p1d = rng.rand(Np,Ng)
    p1dy = rng.rand(Np,Ng)
    v = np.zeros(((Np,Ng,R)))
    a = np.zeros((Np,R))
    a1 = np.zeros((Np,R))
    fieldshift = []
    nfieldchange = [] # number of fields to begin with
    nfieldstay = []
    for j in range(nK):
        nfieldchange.append([])
        nfieldstay.append([])
        for jj in range(nadd):
            nfieldchange[j].append([])
            nfieldstay[j].append([])
    for iN in range(Np):
        for iM in range(Ng):
            #v[iN,iM,:] = grid(x,l[iM/(Ng/len(l))],p1d[iN,iM],phsh=0.,gridmodel='gau',sig=sig)
            v[iN,iM,:] = grid1d_orient(x,l[iM/(Ng/len(l))],ori[iM/(Ng/len(l))],xph=p1d[iN,iM],yph=p1dy[iN,iM],sig=0.16,full=1,mode='exp')
        a[iN,:] = np.dot(w[iN,:],v[iN,:,:])
        print 'Neuron ',iN
        idx = fieldloc(a[iN])
        pk = a[iN,idx]
        sid = pylab.argsort(pk)
        sid = sid[-1::-1]
        pkarr = pk[sid]
        idall = idx[sid]
        pool = range(R)
        for K in range(1,nK+1):   # number of existing fields
            print 'K = ',K
            idarr = idall[:K]
            th = np.sum(pkarr[K-1:K+1])/2  # set threshold by searching for the bumps
            pool = [elem for elem in pool if (elem < idarr[K-1]-10 or elem > idarr[K-1]+10)]
            for k in range(1,nadd+1):   # added fields
                for j in range(nex):  # total number of examples at difficult target locations
                    print 'Example '+str(j)
                    pylab.figure(figsize=[6,7])
                    ax = pylab.subplot(511) # original
                    #pylab.plot(a[iN,:])
                    pylab.plot(np.dot(w[iN,:],v[iN,:,:]))
                    pylab.plot([0,R],[th]*2,'k',lw=1)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(a[iN,:]),np.max(a[iN,:])],'r--',lw=1)
                    ax = pylab.subplot(512) # original thresholded
                    ath = a[iN,:]-th
                    ath[ath<0] = 0
                    pylab.plot(ath)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(ath),np.max(ath)],'r--',lw=1)
                    # weight perturbation
                    ax = pylab.subplot(513)
                    dw = rng.rand(Ng)
                    w1 = w[iN,:] + dw*2/1000. # 0.0001
                    w1 /= np.sum(w1)
                    a1 = np.dot(w1,v[iN,:,:])
                    ath = a1-th
                    ath[ath<0] = 0
                    pylab.plot(ath)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(ath),np.max(ath)],'r--',lw=1)
                    ax = pylab.subplot(514)
                    dw = rng.rand(Ng)
                    w1 = w[iN,:] + dw*2/1000.  # 0.001
                    w1 /= np.sum(w1)
                    a1 = np.dot(w1,v[iN,:,:])
                    ath = a1-th
                    ath[ath<0] = 0
                    pylab.plot(ath)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(ath),np.max(ath)],'r--',lw=1)
                    ax = pylab.subplot(515)
                    dw = rng.rand(Ng)
                    w1 = w[iN,:] + dw*2/1000.  # 0.001
                    w1 /= np.sum(w1)
                    a1 = np.dot(w1,v[iN,:,:])
                    ath = a1-th
                    ath[ath<0] = 0
                    pylab.plot(ath)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(ath),np.max(ath)],'r--',lw=1)
                    ax = pylab.subplot(513) # added activity
                    loc = rng.choice(pool,k,replace=False)
                    dw = np.sum(v[iN,:,loc],0) #np.sum(v[iN,:,loc-25:loc+26],1)  # size of the induced plateau
                    alpha = 1e-5
                    w1 = np.copy(w[iN,:])
                    while np.min(np.dot(v[iN,:,loc],w1)) < th: # when field is not there
                        if alpha > 0.1:
                            print 'alpha = ',alpha,' --> reset in ',iN,K,k,j
                            loc = rng.choice(pool,k,replace=False)
                            dw = np.sum(v[iN,:,loc],0) #np.sum(v[iN,:,loc-25:loc+26],1)  # size of the induced plateau
                            alpha = 1e-5
                        alpha += 1e-5
                        w1 = w[iN,:] + alpha*dw
                        if np.sum(w1) > 0:
                            w1 /= np.sum(w1)
                    a_dw = np.dot(alpha*dw/np.sum(w1),v[iN,:,:])
                    pylab.plot(a_dw)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(a_dw),np.max(a_dw)],'r--',lw=1)
                    for id in loc:
                        pylab.plot([id]*2,[np.min(a_dw),np.max(a_dw)],'g--',lw=1)
                    #print alpha,np.dot(v[iN,:,loc],w1),th
                    ax = pylab.subplot(514) # after learning
                    a1 = np.dot(w1,v[iN,:,:])   # a1 is activity after learning
                    pylab.plot(a1)
                    pylab.plot([0,R],[th]*2,'k',lw=1)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(a[iN,:]),np.max(a[iN,:])],'r--',lw=1)
                    for id in loc:
                        pylab.plot([id]*2,[np.min(a1),np.max(a1)],'g--',lw=1)
                    ax = pylab.subplot(515) # after learning thresholded
                    ath1 = a1-th
                    ath1[ath1<0] = 0
                    pylab.plot(ath1)
                    for id in idarr:
                        pylab.plot([id]*2,[np.min(ath1),np.max(ath1)],'r--',lw=1)
                    for id in loc:
                        pylab.plot([id]*2,[np.min(ath1),np.max(ath1)],'g--',lw=1)
                    loc1 = fieldloc(ath1)
                    for id in idarr:
                        fieldshift.extend(loc1-id)
                    nfieldchange[K-1][k-1].append(loc1.size-k-np.sum(ath1[idarr]>0)+np.sum(ath1[idarr]==0))
                    temp = 0
                    for id in idarr:
                        if np.any(id==loc1):
                            temp += 1
                    nfieldstay[K-1][k-1].append(temp)
                pylab.savefig(figpath+'seed{}pc{}K{}A{}eg{}.svg'.format(seed,iN,K,k,j))
                pylab.close('all')
                with open(datapath+'fieldchange_seed'+str(seed)+'.txt','wb') as f:
                    pickle.dump((fieldshift,nfieldchange,nfieldstay),f)

def readfig7():
    nK = 4
    nadd = 4
    R = 2000
    x = np.arange(R)
    fieldshift = []
    nfieldchange = []
    nfieldstay = []
    l = [31,43]
    ori = [0.39,0.46]
    v1 = grid1d_orient(x,l[0],ori[0],xph=0,yph=0,sig=0.16,full=1,mode='exp')
    v2 = grid1d_orient(x,l[1],ori[1],xph=0,yph=0,sig=0.16,full=1,mode='exp')
    for j in range(nK):
        nfieldchange.append([])
        nfieldstay.append([])
        for jj in range(nadd):
            nfieldchange[j].append([])
            nfieldstay[j].append([])
    for seed in range(1,11):
        with open(datapath+'fieldchange_seed'+str(seed)+'.txt','rb') as f:
            a,b,c = pickle.load(f)
        fieldshift.extend(a)
        for j in range(nK):
            for jj in range(nadd):
                nfieldchange[j][jj].extend(b[j][jj])
                nfieldstay[j][jj].extend(c[j][jj])
    pylab.figure(figsize=[9,3])
    ax = pylab.subplot(131)
    fieldshifthist,dump = pylab.histogram(fieldshift,np.arange(-R+0.5,R))
    pylab.plot(range(-R+1,R),fieldshifthist/(100.*10*4*(1+2+3+4)),label='field %')
    meanfieldchange = np.zeros((nK,nadd))
    pylab.xlim(-70,70)
    pylab.xticks(range(-70,71,35))
    pylab.xlabel('centroid shift')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax = pylab.subplot(132)
    fieldshifthist,dump = pylab.histogram(fieldshift,np.arange(-R+0.5,R))
    pylab.plot(range(-R+1,R),fieldshifthist/float(np.max(fieldshifthist)),label='shift')
    acf1 = np.correlate(v1,v1,'same')
    acf2 = np.correlate(v2,v2,'same')
    acf1 /= np.max(acf1)
    acf2 /= np.max(acf2)
    pylab.plot(range(-R/2,R/2),(acf1+acf2)/2.,'k',label='ACF')
    pylab.legend(loc=2,frameon=False)
    pylab.xlim(-300,300)
    pylab.xticks(range(-300,301,150))
    pylab.xlabel('centroid shift / lag')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax = pylab.subplot(133)
    meanfieldchange = np.zeros((nK,nadd))
    for j in range(nK):
        for jj in range(nadd):
            meanfieldchange[j,jj] = np.mean(nfieldchange[j][jj])
    nadd_vec = []
    nfc_vec = []
    nK_vec = []
    for j in range(nK):
        for jj in range(nadd):
            nadd_vec.extend([jj+1]*len(nfieldchange[j][jj]))
            nfc_vec.extend(nfieldchange[j][jj])
            nK_vec.extend(['$K$='+str(j+1)]*len(nfieldchange[j][jj]))
    import seaborn as sns
    sns.violinplot(nadd_vec,nfc_vec,nK_vec,inner=None,bw=.2,gridsize=100)
    pylab.legend(loc=2,frameon=False)
    pylab.xticks(range(nadd),range(1,nadd+1))
    pylab.ylabel('num of changes')
    pylab.xlabel('num of added fields')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    pylab.subplots_adjust(left=0.1,right=0.9,wspace=0.4,hspace=0.4,bottom=0.2)
    pylab.savefig(figpath+'fig7.svg')

#fig1_2()
#fig3()
#fig5()
#fig5()
#fig6()
#fig6b()
fig8()
#for j in range(1,11):
#    fig7(j)
#readfig7()

max_stars_repo_name: myyim/placecellperceptron

max_issues_repo_name: myyim/placecellperceptron

max_forks_repo_name: myyim/placecellperceptron


text:
# Copyright © VASP Software GmbH,
# Licensed under the Apache License 2.0 (http://www.apache.org/licenses/LICENSE-2.0)
from py4vasp.data._base import DataBase, RefinementDescriptor
import py4vasp._util.documentation as _documentation
import py4vasp.exceptions as exception
import py4vasp._util.convert as _convert
import functools
from fractions import Fraction
import numpy as np

_kpoints_opt_source = """
source : str, optional
    If you used a KPOINTS_OPT file to use a second k-point mesh, you can provide
    a keyword argument `source="kpoints_opt"` to use the k-points defined in that
    file instead of the one specified in KPOINTS.
""".strip()

class Kpoint(DataBase):
    """The **k** points used in the Vasp calculation.

    This class provides utility functionality to extract information about the
    **k** points used by Vasp. As such it is mostly used as a helper class for
    other postprocessing classes to extract the required information, e.g., to
    generate a band structure.

    Parameters
    ----------
    raw_kpoint : RawKpoint
        Dataclass containing the raw **k**-points data used in the calculation.
    """

    to_dict = RefinementDescriptor("_to_dict")
    read = RefinementDescriptor("_to_dict")
    line_length = RefinementDescriptor("_line_length")
    number_lines = RefinementDescriptor("_number_lines")
    distances = RefinementDescriptor("_distances")
    mode = RefinementDescriptor("_mode")
    labels = RefinementDescriptor("_labels")
    __str__ = RefinementDescriptor("_to_string")

    def _to_string(self):
        text = f"""k-points
{len(self._raw_data.coordinates)}
reciprocal"""
        for kpoint, weight in zip(self._raw_data.coordinates, self._raw_data.weights):
            text += "\n" + f"{kpoint[0]} {kpoint[1]} {kpoint[2]}  {weight}"
        return text

    @_documentation.add(
        f"""Read the **k** points data into a dictionary.

Parameters
----------
{_kpoints_opt_source}

Returns
-------
dict
    Contains the coordinates of the **k** points (in crystal units) as
    well as their weights used for integrations. Moreover, some data
    specified in the input file of Vasp are transferred such as the mode
    used to generate the **k** points, the line length (if line mode was
    used), and any labels set for specific points."""
    )
    def _to_dict(self):
        return {
            "mode": self._mode(),
            "line_length": self._line_length(),
            "coordinates": self._raw_data.coordinates[:],
            "weights": self._raw_data.weights[:],
            "labels": self._labels(),
        }

    @_documentation.add(
        f"""Get the number of points per line in the Brillouin zone.

Parameters
----------
{_kpoints_opt_source}

Returns
-------
int
    The number of points used to sample a single line."""
    )
    def _line_length(self):
        if self._mode() == "line":
            return self._raw_data.number
        return len(self._raw_data.coordinates)

    @_documentation.add(
        f"""Get the number of lines in the Brillouin zone.

Parameters
----------
{_kpoints_opt_source}

Returns
-------
int
    The number of lines the band structure contains. For regular meshes this is set to 1."""
    )
    def _number_lines(self):
        return len(self._raw_data.coordinates) // self._line_length()

    @_documentation.add(
        f"""Convert the coordinates of the **k** points into a one dimensional array

For every line in the Brillouin zone, the distance between each **k** point
and the start of the line is calculated. Then the distances of different
lines are concatenated into a single list. This routine is mostly useful
to plot data along high-symmetry lines like band structures.

Parameters
----------
{_kpoints_opt_source}

Returns
-------
np.ndarray
    A reduction of the **k** points onto a one-dimensional array based
    on the distance between the points."""
    )
    def _distances(self):
        cell = self._raw_data.cell.lattice_vectors[-1]
        cartesian_kpoints = np.linalg.solve(cell, self._raw_data.coordinates[:].T).T
        kpoint_lines = np.split(cartesian_kpoints, self._number_lines())
        kpoint_norms = [_line_distances(line) for line in kpoint_lines]
        concatenate_distances = lambda current, addition: (
            np.concatenate((current, addition + current[-1]))
        )
        return functools.reduce(concatenate_distances, kpoint_norms)

    @_documentation.add(
        f"""Get the **k**-point generation mode specified in the Vasp input file

Parameters
----------
{_kpoints_opt_source}

Returns
-------
str
    A string representing which mode was used to setup the k-points."""
    )
    def _mode(self):
        mode = _convert.text_to_string(self._raw_data.mode).strip() or "# empty string"
        first_char = mode[0].lower()
        if first_char == "a":
            return "automatic"
        elif first_char == "e":
            return "explicit"
        elif first_char == "g":
            return "gamma"
        elif first_char == "l":
            return "line"
        elif first_char == "m":
            return "monkhorst"
        else:
            raise exception.RefinementError(
                f"Could not understand the mode '{mode}' when refining the raw kpoints data."
            )

    @_documentation.add(
        f"""Get any labels given in the input file for specific **k** points.

Parameters
----------
{_kpoints_opt_source}

Returns
-------
list[str]
    A list of all the k-points explicitly named in the file or the coordinates of the
    band edges if no name was provided."""
    )
    def _labels(self):
        if self._raw_data.label_indices is not None:
            return self._labels_from_file()
        elif self._mode() == "line":
            return self._labels_at_band_edges()
        else:
            return None

    def _labels_from_file(self):
        labels = [""] * len(self._raw_data.coordinates)
        for label, index in zip(self._raw_data.labels, self._raw_indices()):
            labels[index] = _convert.text_to_string(label.strip())
        return labels

    def _raw_indices(self):
        indices = np.array(self._raw_data.label_indices)
        if self._mode() == "line":
            line_length = self._line_length()
            return line_length * (indices // 2) - (indices + 1) % 2
        else:
            return indices - 1  # convert from Fortran to Python indices

    def _labels_at_band_edges(self):
        line_length = self._line_length()
        band_edge = lambda index: not (0 < index % line_length < line_length - 1)
        return [
            _kpoint_label(kpoint) if band_edge(index) else ""
            for index, kpoint in enumerate(self._raw_data.coordinates)
        ]

def _line_distances(coordinates):
    distances = np.zeros(len(coordinates))
    norms = np.linalg.norm(coordinates[1:] - coordinates[:-1], axis=1)
    distances[1:] = np.cumsum(norms)
    return distances

def _kpoint_label(kpoint):
    fractions = [_to_latex(coordinate) for coordinate in kpoint]
    return f"$[{fractions[0]} {fractions[1]} {fractions[2]}]$"

def _to_latex(float):
    fraction = Fraction.from_float(float).limit_denominator()
    if fraction.denominator == 1:
        return str(fraction.numerator)
    else:
        return f"\\frac{{{fraction.numerator}}}{{{fraction.denominator}}}"