<commit_before>#!/usr/bin/python
"""fit tuning curves to first level results
"""
import matplotlib as mpl
# we do this because sometimes we run this without an X-server, and this backend doesn't need one
mpl.use('svg')
import argparse
import warnings
import os
import json
import matplotlib.pyplot as plt
import pandas as pd
from scipy import optimize
import numpy as np


def log_norm_pdf(x, a, mode, sigma):
    """the pdf of the log normal distribution, with a scale factor

    we parameterize this using the mode instead of mu because we place constraints on the mode
    during optimization
    """
    mu = np.log(mode) + sigma**2
    # the normalizing term isn't necessary, but we keep it here for propriety's sake
    return a * (1/(x*sigma*np.sqrt(2*np.pi))) * np.exp(-(np.log(x)-mu)**2/(2*sigma**2))


def get_tuning_curve_xy(a, mode, sigma, x=None):
    if x is None:
        x = np.logspace(-20, 20, 20000, base=2)
    y = log_norm_pdf(x, a, mode, sigma)
    return x, y


def get_tuning_curve_xy_from_df(df, x=None):
    """given a dataframe with one associated tuning curve, return x and y of that tuning curve
    """
    params = {'x': x}
    for param, param_label in [('a', 'tuning_curve_amplitude'), ('mode', 'tuning_curve_peak'),
                               ('sigma', 'tuning_curve_sigma')]:
        if df[param_label].nunique() > 1:
            raise Exception("Only one tuning curve can be described in df %s!" % df)
        params[param] = df[param_label].unique()[0]
    return get_tuning_curve_xy(**params)


def log_norm_describe_full(a, mode, sigma):
    """calculate and return many descriptions of the log normal distribution

    returns the mu parameter, bandwidth (in octaves), low and high half max values of log normal
    curve, inf_warning, x and y.

    inf_warning is a boolean which indicates whether we calculate the variance to be infinite. this
    happens when the mode is really small as the result of an overflow and so you should probably
    examine this curve to make sure things are okay

    x and y are arrays of floats so you can plot the tuning curve over a reasonable range.
    """
    mu = np.log(mode) + sigma**2
    # we compute this because the std dev is always larger than the bandwidth, so we can use this
    # to make sure we grab the right patch of x values
    var = (np.exp(sigma**2) - 1) * (np.exp(2*mu + sigma**2))
    inf_warning = False
    if np.isinf(var):
        # if the peak of the curve would be at a *really* low or *really* high value, the variance
        # will be infinite (not really, but because of computational issues) and so we need to
        # handle it separately.
        if np.log2(mode) < 0:
            x = np.logspace(-300, 100, 100000, base=2)
        else:
            x = np.logspace(-100, 300, 100000, base=2)
        inf_warning = True
    else:
        std = np.sqrt(var)
        xmin, xmax = np.floor(np.log2(mode) - 5*std), np.ceil(np.log2(mode) + 5*std)
        x = np.logspace(xmin, xmax, 1000*(xmax - xmin), base=2)
    x, y = get_tuning_curve_xy(a, mode, sigma, x)
    half_max_idx = abs(y - (y.max() / 2.)).argsort()
    if (not (x[half_max_idx[0]] > mode and x[half_max_idx[1]] < mode) and
       not (x[half_max_idx[0]] < mode and x[half_max_idx[1]] > mode)):
        print(a, mode, sigma)
        raise Exception("Something went wrong with bandwidth calculation! halfmax x values %s and"
                        " %s must lie on either side of max %s!" % (x[half_max_idx[0]],
                                                                    x[half_max_idx[1]], mode))
    low_half_max = np.min(x[half_max_idx[:2]])
    high_half_max = np.max(x[half_max_idx[:2]])
    bandwidth = np.log2(high_half_max) - np.log2(low_half_max)
    return mu, bandwidth, low_half_max, high_half_max, inf_warning, x, y


def log_norm_describe(a, mode, sigma):
    """same as log_norm_describe_full, except we only return the mode and bandwidth
    """
    m, bw, lhm, hhm, inf, x, y = log_norm_describe_full(a, mode, sigma)
    return m, bw


def create_problems_report(fit_problems, inf_problems, save_path):
    """create html report showing problem cases
    """
    plots_save_path = os.path.join(os.path.dirname(save_path), "report_plots", "plot%03d.svg")
    if not os.path.isdir(os.path.dirname(plots_save_path)):
        os.makedirs(os.path.dirname(plots_save_path))
    with open(os.path.join(os.path.dirname(os.path.realpath(__file__)), 'report_templates.json')) as f:
            report_template_strs = json.load(f)
    report_text = report_template_strs['HEAD']
    for title, problems in zip(["fitting curves", "finding bandwidths"],
                               [fit_problems, inf_problems]):
        report_text += "<h2>Those cases that had problems %s</h2>" % title<commit_msg>two bugfixes to tuning_curves.py<commit_after>#!/usr/bin/python
"""fit tuning curves to first level results
"""
import matplotlib as mpl
# we do this because sometimes we run this without an X-server, and this backend doesn't need one
mpl.use('svg')
import argparse
import warnings
import os
import json
import matplotlib.pyplot as plt
import pandas as pd
from scipy import optimize
import numpy as np


def log_norm_pdf(x, a, mode, sigma):
    """the pdf of the log normal distribution, with a scale factor

    we parameterize this using the mode instead of mu because we place constraints on the mode
    during optimization
    """
    mu = np.log(mode) + sigma**2
    # the normalizing term isn't necessary, but we keep it here for propriety's sake
    return a * (1/(x*sigma*np.sqrt(2*np.pi))) * np.exp(-(np.log(x)-mu)**2/(2*sigma**2))


def get_tuning_curve_xy(a, mode, sigma, x=None):
    if x is None:
        x = np.logspace(-20, 20, 20000, base=2)
    y = log_norm_pdf(x, a, mode, sigma)
    return x, y


def get_tuning_curve_xy_from_df(df, x=None):
    """given a dataframe with one associated tuning curve, return x and y of that tuning curve
    """
    params = {'x': x}
    for param, param_label in [('a', 'tuning_curve_amplitude'), ('mode', 'tuning_curve_peak'),
                               ('sigma', 'tuning_curve_sigma')]:
        if df[param_label].nunique() > 1:
            raise Exception("Only one tuning curve can be described in df %s!" % df)
        params[param] = df[param_label].unique()[0]
    return get_tuning_curve_xy(**params)


def log_norm_describe_full(a, mode, sigma):
    """calculate and return many descriptions of the log normal distribution

    returns the mu parameter, bandwidth (in octaves), low and high half max values of log normal
    curve, inf_warning, x and y.

    inf_warning is a boolean which indicates whether we calculate the variance to be infinite. this
    happens when the mode is really small as the result of an overflow and so you should probably
    examine this curve to make sure things are okay

    x and y are arrays of floats so you can plot the tuning curve over a reasonable range.
    """
    mu = np.log(mode) + sigma**2
    # we compute this because the std dev is always larger than the bandwidth, so we can use this
    # to make sure we grab the right patch of x values
    var = (np.exp(sigma**2) - 1) * (np.exp(2*mu + sigma**2))
    inf_warning = False
    if np.isinf(var):
        # if the peak of the curve would be at a *really* low or *really* high value, the variance
        # will be infinite (not really, but because of computational issues) and so we need to
        # handle it separately. this really shouldn't happen anymore, since I've constrained the
        # bounds of the mode
        if np.log2(mode) < 0:
            x = np.logspace(-300, 100, 100000, base=2)
        else:
            x = np.logspace(-100, 300, 100000, base=2)
        inf_warning = True
    else:
        xmin, xmax = np.floor(np.log2(mode) - 5*sigma), np.ceil(np.log2(mode) + 5*sigma)
        x = np.logspace(xmin, xmax, 1000*(xmax - xmin), base=2)
    x, y = get_tuning_curve_xy(a, mode, sigma, x)
    half_max_idx = abs(y - (y.max() / 2.)).argsort()
    if (not (x[half_max_idx[0]] > mode and x[half_max_idx[1]] < mode) and
       not (x[half_max_idx[0]] < mode and x[half_max_idx[1]] > mode)):
        print(a, mode, sigma)
        raise Exception("Something went wrong with bandwidth calculation! halfmax x values %s and"
                        " %s must lie on either side of max %s!" % (x[half_max_idx[0]],
                                                                    x[half_max_idx[1]], mode))
    low_half_max = np.min(x[half_max_idx[:2]])
    high_half_max = np.max(x[half_max_idx[:2]])
    bandwidth = np.log2(high_half_max) - np.log2(low_half_max)
    return mu, bandwidth, low_half_max, high_half_max, inf_warning, x, y


def log_norm_describe(a, mode, sigma):
    """same as log_norm_describe_full, except we only return the mode and bandwidth
    """
    m, bw, lhm, hhm, inf, x, y = log_norm_describe_full(a, mode, sigma)
    return m, bw


def create_problems_report(fit_problems, inf_problems, save_path):
    """create html report showing problem cases
    """
    plots_save_path = os.path.join(save_path.replace('.csv', '') + "_report_plots", "plot%03d.svg")
    if not os.path.isdir(os.path.dirname(plots_save_path)):
        os.makedirs(os.path.dirname(plots_save_path))
    with open(os.path.join(os.path.dirname(os.path.realpath(__file__)), 'report_templates.json')) as f:
            report_template_strs = json.load(f)
    report_text = report_template_strs['HEAD']
    for title, problems in zip(["fitting curves", "finding bandwidths"],
                               [fit_problems, inf_problems]):
        report_text += "<h2>Those cases that had problems %s</h2>" % title
{'fraction_non_alphanumeric': 0.11091603841784571, 'fraction_numerical': 0.015491066818134875, 'mean_word_length': 3.009937888198758, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '728209', 'n_tokens_mistral': 3272, 'n_tokens_neox': 3051, 'n_words': 1059}