text:
import matplotlib.pyplot as plt
import numpy as np
import emcee
import corner
import random
import math
import subprocess
from astropy.io import ascii
import pickle
from matplotlib.ticker import MaxNLocator
import sys
import idlsave
from scipy.stats.kde import gaussian_kde
import scipy.stats as stats
import matplotlib.mlab as mlab
import tables
from scipy.interpolate import interp1d
from chainconsumer import ChainConsumer
from multiprocessing import Pool
import os
import time

# -------------------------------------------------------------------------#
## load local  modules
from settle import settle
from burstrain import *
from run_model import runmodel
from get_data import get_obs
from mrprior import mr_prior
from get_data import *
from initialise import init

## Now we define the functions that emcee requires
# define likelihood as a function of theta, x, y and yerr as this is what emcee expects as the inputs

ndim, nwalkers, nsteps, run_id, theta, x, y, yerr, tref, bstart, pflux, pfluxe, tobs, numburstssim, numburstsobs, bc, ref_ind, gti_checking, fluen, restart = init()

def lnlike(theta, x, y, yerr):

    # define y = "data" parameters

    for x, i in zip(
        [x for x in range(0, len(bstart) - 1) if x != ref_ind],
        [i for i in range(0, len(bstart) - 1) if i != ref_ind],
    ):
        globals()["t%s" % i] = y[x]
    for x, i in zip(
        range(len(bstart) - 1, len(fluen) + len(bstart) - 1), range(0, len(bstart))
    ):
        globals()["Eb%s" % i] = y[x]
    for x, i in zip(
        range(len(fluen) + len(bstart) - 1, len(y)), range(0, len(bstart - 1))
    ):
        globals()["a%s" % i] = y[x]

    # define yerr as variance terms (errors) for our data parameters (listed in same order as for y)
    # *note that we have to enter three time errors for the code to work however in reality the error should be the same for all of them (st0, st2 and st3 are really dummy parameters)

    for x, i in zip(
        [x for x in range(0, len(bstart) - 1) if x != ref_ind],
        [i for i in range(0, len(bstart) - 1) if i != ref_ind],
    ):
        globals()["st%s" % i] = yerr[x]
    for x, i in zip(
        range(len(bstart) - 1, len(fluen) + len(bstart) - 1), range(0, len(bstart))
    ):
        globals()["sEb%s" % i] = yerr[x]
    for x, i in zip(
        range(len(fluen) + len(bstart) - 1, len(y)), range(0, len(bstart - 1))
    ):
        globals()["sa%s" % i] = yerr[x]

    # define theta = model parameters, which we define priors for

    X, Z, Q_b, f_a, f_E, r1, r2, r3, mass, radius = theta

    # Instead of treating s_t as a parameter, we just hardwire it here

    s_t = 10.0 / 1440.0

    # call model from IDL code defined as modeldata(base, z, x, r1, r2 ,r3)
    if gti_checking == 1:
        model, valid = runmodel(
        theta, y, tref, bstart, pflux, pfluxe, tobs, numburstssim, ref_ind, gti_checking,gti_start=st, gti_end=et
        )
    else:
        model, valid = runmodel(
            theta, y, tref, bstart, pflux, pfluxe, tobs, numburstssim, ref_ind, gti_checking
        )

    if not valid:
        return -np.inf, model

    # multiplying by scaling factors to match with the data
    model[len(bstart) - 1 : len(fluen) + len(bstart) - 1] *= r3
    model[len(fluen) + len(bstart) - 1 : len(y)] *= r2

    # To simplify final likelihood expression we define inv_sigma2 for each data parameter that describe the error.
    # The variance (eg sEb0) is underestimated by some fractional amount, f, for each set of parameters.

    sEb = yerr[len(bstart) - 1 : len(fluen) + len(bstart) - 1]
    sa = yerr[len(fluen) + len(bstart) - 1 : len(yerr)]

    inv_sigma2 = []
    for i in range(0, len(bstart) - 1):
        inv_sigma2.append(1.0 / (s_t ** 2))
    for i in range(0, len(bstart)):
        inv_sigma2.append(1.0 / ((sEb[i] * f_E) ** 2))
    for i in range(0, len(bstart) - 1):
        inv_sigma2.append(1.0 / ((sa[i] * f_a) ** 2))

    # Final likelihood expression
    cpts = (y - (model)) ** 2 * inv_sigma2 - (np.log(inv_sigma2))

    # Test if the result string is defined here. It is, so we return the selected elements of result
    # instead of the downselection in model

    base = Q_b
    z = Z
    x = X
    r1 = r1
    r2 = r2
    r3 = r3
    mass = mass
    radius = radius

    model2 = generate_burst_train(
        base,
        z,
        x,
        r1,
        r2,
        r3,
        mass,
        radius,
        bstart,
        pflux,
        pfluxe,
        tobs,
        numburstssim,
        run=run,
        double=double,
        debug=debug,
    )

    #model2 =  np.string_(model2, dtype='S1000')
    model2 = str(model2).encode('ASCII')

    # Now also return the model
    return -0.5 * np.sum(cpts), model2

# -------------------------------------------------------------------------#
# This is the expression for the prior on the model parameters.

# Define priors for theta. mr prior function is located in mrprior.py

def lnZprior(z):
    # This beta function for the metallicity prior is from Andy Casey and is an approximation of the metallicity of a mock galaxy
    # at 2.5-4.5 kpc for the location of 1808. Assuming ZCNO = 0.01 is average value.
    from scipy import stats
    import numpy as np

    beta = stats.beta
    ZCNO = 0.01

    return np.log(
        beta(10, 3).pdf((np.log10(z / ZCNO) + 3) / 3.75) / (3.75 * np.log(10) * z)
    )

def lnprior(theta):
    import numpy as np

    X, Z, Q_b, f_a, f_E, r1, r2, r3, mass, radius = theta

    if (
        0.00001 < X < 0.76
        and 0.00001 < Z < 0.056
        and 0.000001 <= Q_b < 5.0
        and 0 < f_a < 100
        and 0 < f_E < 100
        and 0.005 < r1 < 1.0
        and 0.005 < r2 < 3.0
        and 0 < r3 * 1e3 < 1000
        and 1.15 < mass < 2.5
        and 9 < radius < 17
    ):  # upper bound and lower bounds of each parameter defined here. Bounds were found by considering an estimated value for each parameter then giving reasonable limits.
        return 0.0 + lnZprior(Z) + mr_prior(mass, radius)
    return -np.inf

# -------------------------------------------------------------------------#
# Finally we combine the likelihood and prior into the overall lnprob function, called by emcee

# define lnprob, which is the full log-probability function
def lnprob(theta, x, y, yerr):
    import numpy as np

    lp = lnprior(theta)

    # Now also returns the model, to accumulate along with the likelihoods

    like, model = lnlike(theta, x, y, yerr)

    if (not np.isfinite(lp)) or (not np.isfinite(like)):
        return -np.inf, -np.inf, model

    # we return the logprobability as well as the theta parameters at this point so we can extract results later
    return lp + like, lp, model

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

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


text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
from numpy.testing import assert_allclose, assert_equal
from astropy.convolution import Box2DKernel
from ...utils.testing import requires_dependency
from .. import SkyMask

@requires_dependency('scipy')
def test_random_creation():
    exclusion = SkyMask.empty(nxpix=300, nypix=100)
    exclusion.fill_random_circles(n=6, max_rad=10)
    assert exclusion.data.shape[0] == 100

    excluded = np.where(exclusion.data == 0)
    assert excluded[0].size != 0

@requires_dependency('scipy')
def test_distance_image():
    mask = SkyMask.empty(nxpix=3, nypix=2)
    distance = mask.distance_image.data
    assert_allclose(distance, -1e10)

    mask = SkyMask.empty(nxpix=3, nypix=2, fill=1.)
    distance = mask.distance_image.data
    assert_allclose(distance, 1e10)

    data = np.array([[0., 0., 1.], [1., 1., 1.]])
    mask = SkyMask(data=data)
    distance = mask.distance_image.data
    expected = [[-1, -1, 1], [1, 1, 1.41421356]]
    assert_allclose(distance, expected)

@requires_dependency('scipy')
def test_open():
    mask = SkyMask.empty(nxpix=5, nypix=5)
    mask.data[2:3, 2:3] = 1
    structure = Box2DKernel(3).array
    mask = mask.open(structure)
    assert_equal(mask.data, 0)

@requires_dependency('scipy')
def test_close():
    mask = SkyMask.empty(nxpix=5, nypix=5)
    mask.data[1:-1, 1:-1] = 1
    mask.data[2, 2] = 0
    structure = Box2DKernel(3).array
    mask = mask.close(structure)
    assert mask.data.sum() == 9

@requires_dependency('scipy')
def test_erode():
    mask = SkyMask.empty(nxpix=5, nypix=5)
    mask.data[1:-1, 1:-1] = 1
    structure = [[1, 1, 1],
                 [1, 1, 1],
                 [1, 1, 1]]
    mask = mask.erode(structure)
    assert mask.data[2, 2] == 1
    assert mask.data.sum() == 1

@requires_dependency('scipy')
def test_dilate():
    mask = SkyMask.empty(nxpix=5, nypix=5)
    mask.data[2, 2] = 1
    structure = [[1, 1, 1],
                 [1, 1, 1],
                 [1, 1, 1]]
    mask = mask.dilate(structure)
    assert mask.data.sum() == 9

ext: py

max_stars_repo_name: grburgess/gammapy

max_issues_repo_name: grburgess/gammapy

max_forks_repo_name: grburgess/gammapy

reason: import numpy,from numpy,from astropy


text:
import sys
import numpy as np
try:
    import caffe
    caffe_imported = True
except ImportError:
    caffe_imported = None
if caffe_imported is None:
    raise ImportError("Please, install caffe.")
import os
import sys
import time
import cv2
TFEAT_PATCH_SIZE = 32
TFEAT_DESC_SIZE = 128
TFEAT_BATCH_SIZE = 256
PROTOTXT_FNAME = 'TFeatLiberty.prototxt'
WEIGHT_FNAME = 'TFeatLiberty.caffemodel'

def preprocess_patch(patch):
    out = cv2.resize(patch, (TFEAT_PATCH_SIZE, TFEAT_PATCH_SIZE)).astype(np.float32)
    return out.reshape(1, TFEAT_PATCH_SIZE, TFEAT_PATCH_SIZE)

def extract_tfeats(net,patches):
    n_desc = len(patches)
    n_batches = int(np.ceil(float(n_desc) / TFEAT_BATCH_SIZE))
    descriptors = np.zeros((n_desc,TFEAT_DESC_SIZE))
    for i in range(n_batches):
        start = i * TFEAT_BATCH_SIZE
        if i < n_batches - 1:
            end = (i + 1) * TFEAT_BATCH_SIZE
        else:
            end = n_desc
        current_batch_size = end - start;
        currect_patches = patches[start:end, :, :, :]
        net.blobs['data'].reshape(current_batch_size, 1, TFEAT_PATCH_SIZE, TFEAT_PATCH_SIZE)
        net.blobs['data'].data[...] = currect_patches
        descriptors[start:end,:] = net.forward()['tanh3'].reshape(current_batch_size, TFEAT_DESC_SIZE)
    return descriptors

if __name__ == '__main__':
    if len(sys.argv) < 2:
        print('Wrong input format. Try python extract_descriptors_from_hpatches_caffe.py ../imgs/ref.png ref_caffe.TFEAT')
        sys.exit(1)

    input_img_fname = sys.argv[1]
    output_fname = sys.argv[2]

    t = time.time()
    image = cv2.imread(input_img_fname,0) #hpatch image is patch column 65*n x 65
    h,w = image.shape
    n_patches = h/w
    print('{0} patches to describe in {1}'.format(n_patches, input_img_fname))
    patches = np.zeros((n_patches,1,TFEAT_PATCH_SIZE, TFEAT_PATCH_SIZE))
    for i in range(n_patches):
        patches[i,:,:,:] = preprocess_patch(image[i*(w): (i+1)*(w), 0:w])
    caffe.set_mode_gpu()
    net = None
    net = caffe.Net(PROTOTXT_FNAME, WEIGHT_FNAME, caffe.TEST)
    print 'Initialization and preprocessing time', time.time() - t
    t = time.time()
    out_descs = extract_tfeats(net,patches)
    print 'extraction time', time.time() - t
    np.savetxt(output_fname, out_descs, delimiter=' ', fmt='%10.5f')

max_stars_repo_name: mhaboali/tfeat

reason: import numpy


text:
# -*- coding: utf-8 -*-

import numpy as np

from .utils import nasa_ames_1001_rw as rw
from .utils import nasa_ames_1001_tools as tools

###############################################################################

class ffi_1001(object):
    """A class to work with NASA Ames files with format index 1001 in Python.

    Parameters
    ----------
    filename : str or pathlib.Path
        filepath and -name of the source.
    sep : str, optional
        General delimiter. The default is " ".
    sep_com : str, optional
        Delimiter used exclusively in comment block. The default is ";".
    sep_data : str, optional
        Delimiter used exclusively in data block. The default is "\t".
    auto_nncoml : bool, optional
        Automatically determine number of lines in comment blocks.
        The default is True.
    strip_lines : bool, optional
        Remove surrounding whitespaces from all lines before parsing.
        The default is True.
    rmv_repeated_seps : bool, optional
        Strip repeated delimiters (e.g. double space). The default is False.
    vscale_vmiss_vertical : bool, optional
        VSCALE and VMISS parameters are arranged vertically over multiple
        lines (1 entry per line) instead of in one line each.
        The default is False.
    vmiss_to_None : bool, optional
        Set True if missing values should be replaced with None.
        The default is False.
    ensure_ascii : bool, optional
        Allow only ASCII characters. The default is True.

    Returns
    -------
    FFI 1001 class instance.
    """

    # initialize class instance with some parameters:
    __INID = {
        "NLHEAD": 14,
        "ONAME": "",
        "ORG": "",
        "SNAME": "",
        "MNAME": "",
        "IVOL": 1,
        "NVOL": 1,
        "DATE": (1970, 1, 1),
        "RDATE": (1970, 1, 1),
        "DX": 0,
        "XNAME": "",
        "NV": 0,
        "VSCAL": 1.0,
        "VMISS": -9999,
        "_VNAME": "",
        "NSCOML": 0,
        "_SCOM": "",
        "NNCOML": 0,
        "_NCOM": "",
        "_X": "",
        "V": "",
        "_FFI": 1001,
        "SRC": "",
    }

    __KEYS = list(__INID.keys())

    def __init__(self, filename=None, **kwargs):
        """Constructor method"""
        for k, v in self.__INID.items():
            setattr(self, k, v)
        if filename is not None:
            self.__from_file(filename, **kwargs)

    def __repr__(self):
        """Representer"""
        s = "NASA Ames 1001\n"
        s += f"SRC: {self.SRC}\n---\n"
        s += "\n".join(f"{k}: {getattr(self, k)}" for k in self.__KEYS[3:4])
        s += "\n" + ", ".join(f"{k}: {getattr(self, k)}" for k in self.__KEYS[7:9])
        s += f"\nXNAME: {self.XNAME}"
        sv = "\n".join(v for v in self.VNAME) if self.VNAME else None
        s += f"\nNV: {self.NV}, VNAMES:\n{sv}"
        return s

    def __str__(self):
        """Stringer"""
        s = "NASA Ames 1001\n"
        s += f"SRC: {self.SRC}\n---\n"
        s += "\n".join(f"{k}: {getattr(self, k)}" for k in self.__KEYS[1:4])
        s += "\n" + ", ".join(f"{k}: {getattr(self, k)}" for k in self.__KEYS[7:9])
        return s

    @property
    def SCOM(self):
        """Special comments block"""
        return self._SCOM

    @SCOM.setter
    def SCOM(self, value):
        self._SCOM, self.NSCOML = value, len(value)
        self.NLHEAD = 14 + self.NSCOML + self.NNCOML + self.NV

    @property
    def NCOM(self):
        """Normal comments block"""
        return self._NCOM

    @NCOM.setter
    def NCOM(self, value):
        self._NCOM, self.NNCOML = value, len(value)
        self.NLHEAD = 14 + self.NSCOML + self.NNCOML + self.NV

    @property
    def VNAME(self):
        """Variables name block"""
        return self._VNAME

    @VNAME.setter
    def VNAME(self, value):
        self._VNAME, self.NV = value, len(value)
        self.NLHEAD = 14 + self.NSCOML + self.NNCOML

    @property
    def X(self):
        """Independent variable"""
        return self._X

    @X.setter
    def X(self, xarr):
        self._X = xarr
        # calculate dx as unique diffs in X,
        # unique with floats might fail, so add a round to 4 decimals:
        dx = np.unique(np.diff(np.array(xarr, dtype=float)).round(4))
        # let dx be 0 if there's more than one unique diff
        dx = dx[0] if dx.size == 1 else 0
        # use an integer if dx is close to its integer value, else float:
        dx = int(dx) if np.isclose(dx, int(dx)) else dx
        self.DX = dx

    # ------------------------------------------------------------------------------
    def __from_file(self, file, **kwargs):
        """Load NASA Ames 1001 from text file."""
        nadict = rw.na1001_cls_read(file, **kwargs)
        for k in self.__KEYS:
            setattr(self, k, nadict[k])

    # ------------------------------------------------------------------------------
    def to_file(self, file, **kwargs):
        """Write NASA Ames 1001 file from populated ffi_1001 class.

        Parameters
        ----------
        file : str or pathlib.Path
            filepath and -name of the destination.
        sep : str, optional
            General delimiter. The default is " ".
        sep_com : str, optional
            Delimiter to use in the commentary section. The default is ";".
        sep_data : str, optional
            Delimiter to separate data columns. The default is "\t".
        overwrite : int, optional
            Set to 1 to overwrite existing files. The default is 0 (no overwrite).
        verbose : bool, optional
            Verbose print output if True. The default is False.

        Returns
        -------
        int
            0 -> failed, 1 -> successful write, 2 -> successful overwrite.

        """
        io = rw.na1001_cls_write(file, self.__dict__, **kwargs)
        return io

    # ------------------------------------------------------------------------------
    def to_dict_nparray(self, **kwargs):
        """Make dictionary of numpy 1D arrays from ffi_1001 class instance.

        Parameters
        ----------
        selVnames : list of string, optional
            VNAMEs to be converted. The default is None.
        return_pddf : bool, optional
            Set to True to make the function return a Pandas dataframe.
            The default is False.
        splitVname : string, optional
            Where to split entries in VNAME. The default is ';'.
        splitIdx : int, optional
            Which part of split result to use, see splitVname. The default is 0.
        xdtype : data type, optional
            Data type for independent variable X. The default is np.float.
        vdtype : data type, optional
            Data type for dependent variable(s). The default is np.float.
        vmiss : missing value identifier, optional
            The default is np.NaN.

        Returns
        -------
        dict
            dictionary holding numpy arrays for variables from the na1001 class
            instance.

        """
        return tools.naDict_2_npndarr(self.__dict__, **kwargs)

    # ------------------------------------------------------------------------------
    def to_pddf(self, **kwargs):
        """Make a pandas DataFrame from NA 1001 data (X and V).

        Parameters
        ----------
        sep_colhdr : str, optional
            separator used in column header (last line of NCOM). Default is tab.
        idx_colhdr : int, optional
            look for column header in NCOM at index idx_colhdr. Default is -1.
        dtype : numpy array data type, optional
            data type to use for conversion to DataFrame. Default is np.float64.
        add_datetime_index: boolean, optional
            add a DateTime index to the df. The default is False.

        Returns
        -------
        pandas.DataFrame
            dataframe with a column for X and one for each parameter in V.

        """
        return tools.naDict_2_pddf(self.__dict__, **kwargs)

max_stars_repo_name: FObersteiner/pyFuppes

max_issues_repo_name: FObersteiner/pyFuppes

max_forks_repo_name: FObersteiner/pyFuppes


text:
import importlib
import solarforecast
solarforecast = importlib.reload(solarforecast)
from solarforecast import FileInf
from solarforecast import SolcastHistorical
from solarforecast import SolcastDataForecast
from solarforecast import EtapData
import os  #access files and so on
import matplotlib
import matplotlib.pyplot as plt
from solarforecast import SolarF
import numpy as np
import keras
import pandas as pd
import datetime
import time
from time import sleep
import schedule
import scipy.stats

#%%
forecasted_features = ['Ghi', 'Ghi90', 'Ghi10', 'Ebh', 'Dni', 'Dni10', 'Dni90', 'Dhi',
       'air_temp', 'Zenith', 'Azimuth', 'cloud_opacity', 'period_end',
       'Period']
historical_features = ['PeriodEnd', 'PeriodStart', 'Period', 'AirTemp', 'AlbedoDaily',
       'Azimuth', 'CloudOpacity', 'DewpointTemp', 'Dhi', 'Dni', 'Ebh', 'Ghi',
       'PrecipitableWater', 'RelativeHumidity', 'SnowDepth', 'SurfacePressure',
       'WindDirection10m', 'WindSpeed10m', 'Zenith']

whole_training_features = ['Ghi', 'Ebh', 'Dni', 'Dhi', 'AirTemp', 'CloudOpacity'] #for today which will predict tomorrow
output_feature = ['PV_power']#for the next day power generation

%reload_ext solarforecast

etap_power=EtapData(0.8)
train_times = etap_power.train_data.keys()
test_times = etap_power.test_data.keys()
#historical data from solcast
dst='data/solcast_etap_historical.csv'
hist = SolcastHistorical(dst, train_times, test_times)

#%%
def train_test_by_features(selected_features, hist, etap_power):

    train_features = []
    train_label = []
    for i in hist.train.keys():
        selected_data = hist.train[i][selected_features]
        train_features.append(selected_data.values)
        train_label.append(etap_power.train_data[i]['Avg'].values)

    train_features = np.array(train_features)
    train_label = np.array(train_label)

    test_features = []
    test_label = []
    for i in hist.test.keys():
        selected_data = hist.test[i][selected_features]
        test_features.append(selected_data.values)
        test_label.append(etap_power.test_data[i]['Avg'].values)

    test_features = np.array(test_features)
    test_label = np.array(test_label)
    #max min normalization
    powermax = np.max(train_label)
    powermin = np.min(train_label)

    feature_max = train_features.max(axis=(1,0))
    feature_min = train_features.min(axis=(1,0))

    ##normalize
    x_train = (train_features-feature_min)/(feature_max-feature_min)
    x_test = (test_features-feature_min)/(feature_max-feature_min)

    y_train = (train_label-powermin)/(powermax-powermin)
    y_test = (test_label-powermin)/(powermax-powermin)

    common_features=[feature_max,feature_min,powermax,powermin]
    return x_train, x_test, y_train, y_test,common_features

#%%
# select features for training
scenarios = {
    'whole': ['Ghi', 'Ebh', 'Dni', 'Dhi', 'AirTemp', 'CloudOpacity'],
    'radiations':['Ghi', 'Ebh', 'Dni', 'Dhi'],
    'normal':['Ghi', 'AirTemp', 'CloudOpacity'],
    'minimal':['Ghi', 'CloudOpacity'],
    'Ghi':['Ghi']
}
#%%
mses={}
for sc in scenarios:
    if sc == 'normal':
        print(sc)
        selected_features = scenarios[sc]

        feature_numbers = len(selected_features)
        resolution = 24
        x_train, x_test, y_train, y_test, common_features = train_test_by_features(selected_features, hist, etap_power)
        feature_max, feature_min, powermax, powermin = common_features
        loaded_model = keras.models.load_model('models/'+sc)
        print(loaded_model)
        predicted = loaded_model.predict(x_train)
        mse_error = loaded_model.evaluate(x_train, y_train)
        print(mse_error)
        mses[sc] = mse_error[0]

#%%
## prediciton error normal pdf
error = predicted - y_train
error = error.ravel()
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

data = error

# Fit a normal distribution to the data:
mu, std = norm.fit(data)

# Plot the histogram.
plt.hist(data,bins=30, density=True,  color='blue')

# Plot the PDF.
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mu, std)
plt.plot(x, p, 'r', linewidth=2)
title = "Error fit results: mu = %.2f,  std = %.2f" % (mu, std)
plt.title(title)
plt.legend(['normal pdf fitted for error','normalized error histogram'])
plt.savefig('figs/errorhist.png')
plt.show()
#%%
#prediciton plot
for i, k in enumerate(predicted):
    plt.plot((y_train[i] * (powermax - powermin) + powermin)*1000)
    plt.plot((predicted[i] * (powermax - powermin) + powermin)*1000)
    plt.ylabel('Power(kW)')
    plt.xlabel('Hour of the day')
    plt.legend(['ETAP real solar generation','Predicted solar generation'])
    plt.savefig('figs/realpred/' + str(i) + '.png')
    plt.show()
#%%

    #%%
from keras.models import Model, Sequential
from keras import backend as K

def create_dropout_predict_function(model, dropout):
    """
    Create a keras function to predict with dropout
    model : keras model
    dropout : fraction dropout to apply to all layers

    Returns
    predict_with_dropout : keras function for predicting with dropout
    """

    # Load the config of the original model
    conf = model.get_config()
    # Add the specified dropout to all layers
    for layer in conf['layers']:
        # Dropout layers
        if layer["class_name"] == "Dropout":
            layer["config"]["rate"] = dropout
            print('config-rate')
        # Recurrent layers with dropout
        elif "dropout" in layer["config"].keys():
            layer["config"]["dropout"] = dropout
            print(layer["config"]["dropout"],'config-droupout')
    # Create a new model with specified dropout
    if type(model) == Sequential:
        # Sequential
        newmodel = keras.models.clone_model(model)
        # newmodel.compile(optimizer='adam',
        #                     loss='mse',
        #                     metric='mse')
        # newmodel.set_weights(model.get_weights())

        model_dropout = Sequential.from_config(conf)
        # print('model1')
    else:
        # Functional
        model_dropout = Model.from_config(conf)
        # print('model2')

    model_dropout.set_weights(model.get_weights())
    # model_dropout.compile(optimizer='adam',
    #                       loss='mse')

    # final_conv_layer = get_output_layer(model, "conv5_3")
    # get_output = K.function([model.layers[0].input],
    #                         [final_conv_layer.output, model.layers[-1].output])
    # [conv_outputs, predictions] = get_output([img])
    # Create a function to predict with the dropout on
    predict_with_dropout = K.function([model_dropout.layers[0].input,K.learning_phase()],
                                      [model_dropout.layers[-1].output])

    return predict_with_dropout

#%%
def mc(x_train, iternum, drop):
    mcpred = {}
    prediction_with_dropout = create_dropout_predict_function(loaded_model, drop)
    for i in range(x_train.shape[0]):
        mcpred[i] = []

    for i in range(iternum):
        y = prediction_with_dropout([x_train, 1])
        for j in range(y[0].shape[0]):
            mcpred[j].append((y[0][j] * (powermax - powermin) + powermin)*1000)
    return mcpred
#%%
def ReLU(x):
    a=[]
    for i in x:
        a.append(max(i,0))
    return a
def mean_confidence_interval(data, zp):
    a = 1.0 * np.array(data)
    n = len(a)
    m, se = np.mean(a,axis=0), np.std(a,axis=0)
    h = se* zp/2
    return m, ReLU(m-h), m+h

def confidence_intervals(mcpred, zp):
    CItrain = {}
    for i in range(len(mcpred)):
        mn, lower, upper = mean_confidence_interval(mcpred[i],zp)
        CItrain[i] = [lower, upper ]
    return CItrain

#%%
dropout_levels = [0.1, 0.2, 0.3, 0.4]
CI_level = [0.8, 0.9, 0.95]
zp =[1.28, 1.64, 1.96]
iternum=50

results = {}
for drop in dropout_levels:
    print(drop)
    results[drop] = {}
    mcpred = mc(x_train, iternum,drop)
    for level in zp:
        print(level)
        ci = confidence_intervals(mcpred, level)
        results[drop][level] = ci

#%%
for j in zp:
    os.mkdir('figs/confidence/' + str(j))
    for i in range(16):
        temp = results[0.1][j]
        plt.plot(temp[i][0])
        plt.plot(temp[i][1])
        plt.plot((predicted[i]* (powermax - powermin) + powermin)*1000)
        plt.plot((y_train[i]* (powermax - powermin) + powermin)*1000)
        plt.legend(['lower bound', 'higher bound', 'prediction', 'real'])
        plt.ylabel('Power(kW)')
        plt.xlabel('Hour of the day')
        plt.savefig('figs/confidence/' + str(j) + '/' + str(i) + '.png')
        plt.show()

max_stars_repo_name: arminalgln/solar

max_issues_repo_name: arminalgln/solar

max_forks_repo_name: arminalgln/solar

reason: import numpy,import scipy,from scipy


text:
import numpy as np
import os.path as path
from keras.callbacks import ModelCheckpoint, ReduceLROnPlateau, EarlyStopping, TensorBoard, Callback
from keras.models import Model
from src.base.trainer import BaseTrainer
import src.utils.config as cfg
import src.utils.dirs as dirs
import src.utils.metrics as mtx
import src.data_loader.ogits as DL
import src.model.objectives as obj
import src.utils.io as io

class OGITSModelTrainer(BaseTrainer):
    """
    # Arguments
        model: Keras model, model for training
        data: BaseDataLoader, generate datasets (for train/test/validate)
        config: DotMap, settings for training 'model',
                according to 'training_mode' (i.e. pre-training/finetune)
        train_mode: String. 'pretrain_set' or 'finetune_set'.
        attrib_cls: String. Speech attribute 'manner', 'place' or 'fusion'
        verbose: Integer. 0, 1, or 2. Verbosity mode.
                    0 = silent, 1 = progress bar, 2 = one line per epoch.
        model: not trained or initially pre-trained model
    """

    def __init__(self, model, data, config, **kwargs):
        super(OGITSModelTrainer, self).__init__(model, data, config)

        self.train_mode = kwargs.get('train_mode', 'pretrain_set')
        self.fold_id = kwargs.get('fold_id', 1)
        self.attrib_cls = kwargs.get('attrib_cls', 'manner')
        self.verbose = kwargs.get('verbose', 1)
        self.overwrite = kwargs.get('overwrite', False)

        self.callbacks = []
        mpath = path.join(self.config['path'][self.train_mode], self.attrib_cls)
        self.model_file = cfg.get_model_filename(path=mpath,
                                                 train_mode=self.train_mode,
                                                 fold=self.fold_id)

    def _init_callbacks(self, cv_gen):
        valtor = ModelValidator(batch_gen=cv_gen,
                                metrics=self.config['model'][self.train_mode]['metrics'],
                                monitor=self.config['callback']['monitor'],
                                mode=self.config['callback']['mode'])

        chkp = ModelCheckpoint(
            monitor=self.config['callback']['monitor'],
            filepath=self.model_file,
            mode=self.config['callback']['mode'],
            save_best_only=self.config['callback']['chpt_save_best_only'],
            save_weights_only=self.config['callback']['chpt_save_weights_only'],
            verbose=self.verbose)

        lr_sched = ReduceLROnPlateau(monitor=self.config['callback']['monitor'],
                                     patience=self.config['callback']['lr_patience'],
                                     verbose=self.verbose,
                                     factor=self.config['callback']['lr_factor'],
                                     min_lr=self.config['callback']['lr_min'])

        e_stop = EarlyStopping(monitor=self.config['callback']['monitor'],
                               patience=self.config['callback']['estop_patience'],
                               verbose=self.verbose,
                               mode=self.config['callback']['mode'],
                               min_delta=0.001)

        log_dir = cfg.get_log_dir(self.train_mode, self.fold_id)
        print('LOG: %s' % log_dir)
        tb = TensorBoard(log_dir=log_dir,
                         write_graph=self.config['callback']['tensorboard_write_graph'])

        self.callbacks.append([valtor, chkp, lr_sched, e_stop, tb])

    @staticmethod
    def is_mfom_objective(model):
        return (model.loss in obj.MFOM_OBJECTIVES) or \
               (model.loss in obj.MFOM_OBJECTIVES.values())

    def train(self):
        if not dirs.check_file(self.model_file) or self.overwrite:
            # batch generators
            train_gen = DL.batch_handler(batch_type=self.config['model'][self.train_mode]['batch_type'],
                                         data_file=self.data.feat_file,
                                         fold_lst=self.data.meta_data.fold_list(self.fold_id, 'train'),
                                         config=self.config['model'][self.train_mode],
                                         meta_data=self.data.meta_data)

            cv_gen = DL.batch_handler(batch_type='sed_validation',
                                      data_file=self.data.feat_file,
                                      fold_lst=self.data.meta_data.fold_list(self.fold_id, 'validation'),
                                      config=self.config['model'][self.train_mode],
                                      meta_data=self.data.meta_data)
            self._init_callbacks(cv_gen)

            smp_size = train_gen.samples_number()
            print('[INFO] Epoch size: %d observations' % smp_size)

            batch_sz = self.config['model'][self.train_mode]['batch']
            nepo = self.config['model'][self.train_mode]['n_epoch']

            def mfom_batch_wrap(xy_gen):
                for x, y in xy_gen.batch():
                    yield [y, x], y

            gen = mfom_batch_wrap(train_gen) \
                if self.is_mfom_objective(self.model) \
                else train_gen.batch()

            self.model.fit_generator(gen,
                                     steps_per_epoch=smp_size // batch_sz,
                                     nb_epoch=nepo,
                                     verbose=self.verbose,
                                     workers=1,
                                     callbacks=self.callbacks)
            train_gen.stop()
            cv_gen.stop()
        else:
            print('[INFO] There is %s model: %s' % (self.train_mode.upper(), self.model_file))

    def test(self):
        test_gen = DL.batch_handler(batch_type='sed_validation',
                                    data_file=self.data.feat_file,
                                    fold_lst=self.data.meta_data.fold_list(self.fold_id, 'test'),
                                    config=self.config['model'][self.train_mode],
                                    meta_data=self.data.meta_data)

        history = ModelValidator.validate_model(model=self.model,
                                                batch_gen=test_gen,
                                                metrics=self.config['model'][self.train_mode]['metrics'],
                                                attrib_cls=self.attrib_cls)
        print(history)
        test_gen.stop()

    def evaluate(self):
        test_gen = DL.batch_handler(batch_type='sed_evaluation',
                                    data_file=self.data.feat_file,
                                    fold_lst=self.data.meta_data.fold_list(self.fold_id, 'test'),
                                    config=self.config['model'][self.train_mode],
                                    meta_data=self.data.meta_data)

        cut_model = self.model
        if OGITSModelTrainer.is_mfom_objective(self.model):
            input = self.model.get_layer(name='input').output
            preact = self.model.get_layer(name='output').output
            cut_model = Model(input=input, output=preact)

        n_class = cut_model.output_shape[-1]
        cnt = 0
        scores = {}
        for fn, X_b in test_gen.batch():
            ps = cut_model.predict_on_batch(X_b)
            ps = ps.reshape((-1, n_class))
            fid = path.basename(fn)
            fid = path.splitext(fid)[0]
            scores[fid] = ps
            print('Processed: %d' % cnt)
            cnt += 1
        print('Evaluated files: %d' % cnt)
        io.save_ark('%s.ark' % self.attrib_cls, scores)

class ModelValidator(Callback):
    def __init__(self, batch_gen, metrics, monitor, mode):
        super(ModelValidator, self).__init__()
        self.batch_gen = batch_gen
        self.metrics = metrics
        self.monitor = monitor
        self.best_epoch = 0

        if mode == 'min':
            self.monitor_op = np.less
            self.best_acc = np.Inf
        elif mode == 'max':
            self.monitor_op = np.greater
            self.best_acc = 0.
        else:
            raise AttributeError('[ERROR] ModelValidator mode %s is unknown')

    def on_train_begin(self, logs=None):
        super(ModelValidator, self).on_train_begin(logs)
        vs = ModelValidator.validate_model(self.model, self.batch_gen, self.metrics)
        for k, v in vs.items():
            logs[k] = np.float64(v)

        print(' BEFORE TRAINING: Validation loss: %.4f, Validation %s: %.4f / best %.4f' % (
            vs['val_loss'], self.monitor.upper(), vs[self.monitor], self.best_acc))
        print(logs)

        if self.monitor_op(logs[self.monitor], self.best_acc):
            self.best_acc = logs[self.monitor]
            self.best_epoch = -1

    def on_epoch_end(self, epoch, logs={}):
        super(ModelValidator, self).on_epoch_end(epoch, logs)
        vs = ModelValidator.validate_model(self.model, self.batch_gen, self.metrics)
        for k, v in vs.items():
            logs[k] = np.float64(v)

        print(' EPOCH %d. Validation loss: %.4f, Validation %s: %.4f / best %.4f' % (
            epoch, vs['val_loss'], self.monitor.upper(), vs[self.monitor], self.best_acc))
        print(logs)

        if self.monitor_op(logs[self.monitor], self.best_acc):
            self.best_acc = logs[self.monitor]
            self.best_epoch = epoch

    def on_train_end(self, logs=None):
        super(ModelValidator, self).on_train_end(logs)
        print('=' * 20 + ' Training report ' + '=' * 20)
        print('Best validation %s: epoch %s / %.4f\n' % (self.monitor.upper(), self.best_epoch, self.best_acc))

    @staticmethod
    def validate_model(model, batch_gen, metrics, attrib_cls=None):
        """
        # Arguments
            model: Keras model
            data: BaseDataLoader
            metrics: list of metrics
        # Output
            dictionary with values of metrics and loss
        """
        cut_model = model
        if OGITSModelTrainer.is_mfom_objective(model):
            input = model.get_layer(name='input').output
            preact = model.get_layer(name='output').output
            cut_model = Model(input=input, output=preact)

        n_class = cut_model.output_shape[-1]
        _, bnd, wnd, ch = batch_gen.batch_shape()
        y_true, y_pred = np.empty((0, wnd, n_class)), np.empty((0, wnd, n_class))

        loss, cnt = 0, 0
        for fn, X_b, Y_b in batch_gen.batch():
            ps = cut_model.predict_on_batch(X_b)
            y_pred = np.vstack([y_pred, ps])
            y_true = np.vstack([y_true, Y_b])
            # NOTE: it is approximated loss, caz Y is fed
            if OGITSModelTrainer.is_mfom_objective(model):
                X_b = [Y_b, X_b]
            l = model.test_on_batch(X_b, Y_b)
            loss += l
            cnt += 1
        vals = {'val_loss': loss / cnt}
        y_true = y_true.reshape((-1, n_class))
        y_pred = y_pred.reshape((-1, n_class))

        if attrib_cls == 'fusion':
            fmanner_ids = [2, 4, 9, 11, 12, 14, 15]
            y_true_m = y_true[:, fmanner_ids]
            y_pred_m = y_pred[:, fmanner_ids]

            v_max, t_max = 0, 0
            for t in [0.5]:  # np.linspace(0.2, 0.7, 50):
                p = mtx.step(y_pred_m, threshold=t)
                v = mtx.micro_f1(y_true_m, p)
                if v_max < v:
                    v_max = v
                    t_max = t
            print('F1 fusion-manner:', v_max, t_max)
            eer = np.mean(mtx.class_wise_eer(y_true_m, y_pred_m))
            print('EER fusion-manner:', eer)

            fplace_ids = [0, 1, 3, 5, 6, 7, 8, 10, 11, 13]
            y_true_p = y_true[:, fplace_ids]
            y_pred_p = y_pred[:, fplace_ids]

            v_max, t_max = 0, 0
            for t in [0.5]:  # np.linspace(0.2, 0.7, 50):
                p = mtx.step(y_pred_p, threshold=t)
                v = mtx.micro_f1(y_true_p, p)
                if v_max < v:
                    v_max = v
                    t_max = t
            print('F1 fusion-place:', v_max, t_max)
            eer = np.mean(mtx.class_wise_eer(y_true_p, y_pred_p))
            print('EER fusion-place:', eer)

        for m in metrics:
            if m == 'micro_f1':
                v_max, t_max = 0, 0
                for t in [0.5]:  # np.linspace(0.2, 0.7, 50):
                    p = mtx.step(y_pred, threshold=t)
                    v = mtx.micro_f1(y_true, p)
                    if v_max < v:
                        v_max = v
                        t_max = t
                vals[m] = v_max
                print('Optimal threshold: ', vals[m], t_max)
            elif m == 'pooled_eer':
                p = y_pred.flatten()
                y = y_true.flatten()
                vals[m] = mtx.eer(y, p)
            elif m == 'class_wise_eer':
                vals[m] = np.mean(mtx.class_wise_eer(y_true, y_pred))
            elif m == 'accuracy':
                p = np.argmax(y_pred, axis=-1)
                y = np.argmax(y_true, axis=-1)
                vals[m] = mtx.pooled_accuracy(y, p)
            else:
                raise KeyError('[ERROR] Such a metric is not implemented: %s...' % m)
        return vals


text:
import numpy as np
import pytest

from fdm import forward_fdm, backward_fdm, central_fdm, FDM
from .util import approx

def test_construction():
    with pytest.raises(ValueError):
        FDM([-1, 0, 1], 3)

@pytest.mark.parametrize("f", [forward_fdm, backward_fdm, central_fdm])
def test_correctness(f):
    approx(f(10, 1)(np.sin, 1), np.cos(1))
    approx(f(10, 2)(np.sin, 1), -np.sin(1))

    approx(f(10, 1)(np.exp, 1), np.exp(1))
    approx(f(10, 2)(np.exp, 1), np.exp(1))

    approx(f(10, 1)(np.sqrt, 1), 0.5)
    approx(f(10, 2)(np.sqrt, 1), -1 / 4)

def cosc(x):
    if x == 0:
        return 0.0
    else:
        # Does this give cancellation errors?
        return (np.cos(np.pi * x) - np.sinc(x)) / x

@pytest.mark.parametrize(
    "f,x,df,atol,atol_central",
    [
        (lambda x: 0, 0, 0, 0, 0),
        (lambda x: x, 0, 1, 5e-15, 5e-16),
        (np.exp, 0, 1, 7.5e-13, 5e-14),
        (np.sin, 0, 1, 1e-14, 5e-16),
        (np.cos, 0, 0, 5e-13, 5e-14),
        (np.exp, 1, np.exp(1), 5e-12, 5e-13),
        (np.sin, 1, np.cos(1), 5e-13, 5e-14),
        (np.cos, 1, -np.sin(1), 7.5e-13, 5e-14),
        (np.sinc, 0, 0, 5e-12, 5e-14),
        # `cosc` is hard.
        (cosc, 0, -np.pi ** 2 / 3, 5e-10, 5e-11),
    ],
)
@pytest.mark.parametrize(
    "constructor,is_central",
    [
        (central_fdm, True),
        (forward_fdm, False),
        (backward_fdm, False),
    ],
)
def test_accuracy(f, x, df, atol, atol_central, constructor, is_central):
    if is_central:
        approx(constructor(6, 1)(f, x), df, atol=atol_central, rtol=0)
        approx(constructor(7, 1)(f, x), df, atol=atol_central, rtol=0)
        if f == cosc:
            # `cosc` is hard.
            approx(constructor(14, 1)(f, x), df, atol=1e-13, rtol=0)
        else:
            approx(constructor(14, 1)(f, x), df, atol=5e-15, rtol=0)
    else:
        approx(constructor(6, 1)(f, x), df, atol=atol, rtol=0)
        approx(constructor(7, 1)(f, x), df, atol=atol, rtol=0)

def test_estimation():
    m = central_fdm(2, 1)

    assert isinstance(m.coefs, np.ndarray)
    assert isinstance(m.df_magnitude_mult, float)
    assert isinstance(m.f_error_mult, float)
    assert m.step is None
    assert m.acc is None

    m.estimate()

    assert isinstance(m.step, float)
    assert isinstance(m.acc, float)

    m(np.sin, 0, step=1e-3)

    assert isinstance(m.step, float)
    assert m.acc is None

def test_adaptation():
    def f(x):
        return 1 / x

    def df(x):
        return -1 / x ** 2

    err1 = np.abs(forward_fdm(3, 1, adapt=0)(f, 1e-3) - df(1e-3))
    err2 = np.abs(forward_fdm(3, 1, adapt=2)(f, 1e-3) - df(1e-3))

    # Check that adaptation helped.
    assert err2 <= 0.01 * err1

    # Check that adaptation gets it right.
    assert err2 <= 1e-4

def test_order_monotonicity():
    err_ref = 1e-4

    for i in range(3, 10):
        method = central_fdm(i, 2, adapt=1)

        # Check order of method.
        assert method.order == i

        err = np.abs(method(np.sin, 1) + np.sin(1))

        # Check that it did better than the previous estimator.
        assert err <= err_ref

        err_ref = err

def test_default_step():
    approx(
        central_fdm(2, 1).estimate().step, np.sqrt(2 * np.finfo(np.float64).eps / 10)
    )

def test_zero_bound_fixed():
    m = central_fdm(2, 1)
    f = np.sin
    x = 0

    assert m.bound_estimator(f, x) == 0
    assert m.bound_estimator(f, x, magnitude=True)[0] > 0
    assert m.bound_estimator(f, x, magnitude=True)[1] > 0

    approx(m(f, x), 1)

def test_zero_bound_zero_error_not_fixed():
    m = central_fdm(2, 1)
    f = lambda _: 0
    x = 0

    assert m.bound_estimator(f, x) == 0
    assert m.bound_estimator(f, x, magnitude=True)[0] == 0
    assert m.bound_estimator(f, x, magnitude=True)[1] == 0

    approx(m(f, x), 0)

def test_step_limiting():
    m = central_fdm(2, 1)
    f = lambda x: np.exp(1e-3 * x)
    x = 0

    step_max = m.estimate().step * 1000
    assert m.estimate(f, x).step == step_max

def test_factor():
    assert (
        central_fdm(2, 1, factor=4).estimate().step
        == 2 * central_fdm(2, 1, factor=1).estimate().step
    )

def test_condition():
    assert (
        central_fdm(2, 1, condition=1).estimate().step
        == 2 * central_fdm(2, 1, condition=4).estimate().step
    )

class _TrackEvals:
    def __init__(self, f):
        self.evals_x = []
        self.f = f

    def __call__(self, x):
        self.evals_x.append(x)
        return self.f(x)

    def true_range(self, x0):
        return np.max(np.abs([x - x0 for x in self.evals_x]))

@pytest.mark.parametrize("factor", [10, 100, 1000])
def test_range_max(factor):
    f = _TrackEvals(np.sin)

    # Determine the range.
    central_fdm(8, 1, adapt=2)(f, 1)
    true_range = f.true_range(1)

    # Limit the range.
    f.evals_x.clear()
    central_fdm(8, 1, adapt=2)(f, 1, max_range=true_range / factor)
    approx(f.true_range(1), true_range / factor)

max_stars_repo_name: wesselb/fdm

max_issues_repo_name: wesselb/fdm

max_forks_repo_name: wesselb/fdm


text:
"""
@author: karthikrao

Partition a graph into two blocks using the Kernighan-Lin algorithm.

"""

import numpy as np
import networkx as nx
from util import *

def get_init_partition(nodes):
    """

    Parameters
    ----------
    nodes : arr_like.
        Vertices of the netlist's graph.

    Returns
    -------
    list, list
        Nodes in partition 1 and 2.

    """
    return list(nodes[:len(nodes)//2]), list(nodes[len(nodes)//2:])

def kernighan_lin_bisection(g, max_iter=10):
    """

    Parameters
    ----------
    g : Networkx undirected graph
    max_iter : int, optional
        Maximum iterations the algorithm has to search for
        optimal partitions before giving up. The default is 10.

    Returns
    -------
    av : list
        Nodes in partition 1.
    bv : list
        Nodes in partition 2.

    """
    nodes = list(g.nodes(data=False))
    a, b = get_init_partition(nodes)
    for epoch in range(max_iter):
        av, bv, gv = [], [], []
        g1 = g.copy()
        for i in range(len(nodes)//2):
            th = -np.inf
            for j in range(len(a)):
                t1 = a[j]
                for k in range(len(b)):
                    t2 = b[k]
                    t3 = compute_gain(g1, t1, t2, a, b)
                    if t3>th:
                        th = t3
                        rnode = []
                        rnode.extend((t1, t2))
            av.append(rnode[0])
            bv.append(rnode[1])
            gv.append(th)
            a.remove(rnode[0])
            b.remove(rnode[1])
            g1.remove_nodes_from(rnode)
        max_sum = np.argmax(np.array(gv))
        av, bv = swap_from_list(av, bv, max_sum)
        a, b = av, bv
    return av, bv

max_stars_repo_name: karthik-r-rao/vlsi-physical-design

max_issues_repo_name: karthik-r-rao/vlsi-physical-design

max_forks_repo_name: karthik-r-rao/vlsi-physical-design

reason: import numpy,import networkx


text:
import os
import numpy as np
from dotmap import DotMap
from collections import OrderedDict

import torch
import torchvision
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.data import DataLoader

from src.datasets import datasets
from src.models.resnet import resnet18
from src.objectives.memory import MemoryBank
from src.models.logreg import LogisticRegression
from src.objectives.instdisc import InstDisc, NCE, Ball, Ring, SimCLR
from src.objectives.moco import MoCo
from src.utils import utils

import pytorch_lightning as pl
torch.autograd.set_detect_anomaly(True)

def create_dataloader(dataset, config, shuffle=True):
    loader = DataLoader(
        dataset,
        batch_size=config.optim_params.batch_size,
        shuffle=shuffle,
        pin_memory=True,
        drop_last=True,
        num_workers=config.data_loader_workers,
    )
    return loader

class PretrainSystem(pl.LightningModule):

    def __init__(self, config):
        super().__init__()
        self.config = config
        self.train_dataset, self.val_dataset = datasets.get_datasets(
            config.data_params.dataset)
        train_labels = self.train_dataset.dataset.targets
        self.train_ordered_labels = np.array(train_labels)
        self.model = self.create_encoder()
        self.memory_bank = MemoryBank(len(self.train_dataset),
                                      self.config.model_params.out_dim)

    def create_encoder(self):
        model = resnet18(low_dim=self.config.model_params.out_dim)
        if self.config.model_params.projection_head:
            mlp_dim = model.fc.weight.size(1)
            model.fc = nn.Sequential(
                nn.Linear(mlp_dim, mlp_dim),
                nn.ReLU(),
                model.fc,
            )
        return model

    def configure_optimizers(self):
        optim = torch.optim.SGD(self.model.parameters(),
                                lr=self.config.optim_params.learning_rate,
                                momentum=self.config.optim_params.momentum,
                                weight_decay=self.config.optim_params.weight_decay)
        return [optim], []

    def forward(self, img):
        return self.model(img)

    def get_losses_for_batch(self, batch, train=True):
        indices, img, _, = batch
        outputs = self.forward(img)
        loss_fn = NCE(indices, outputs, self.memory_bank,
                      k=self.config.loss_params.k,
                      t=self.config.loss_params.t,
                      m=self.config.loss_params.m)
        loss = loss_fn.get_loss()

        if train:
            with torch.no_grad():
                new_data_memory = loss_fn.updated_new_data_memory()
                self.memory_bank.update(indices, new_data_memory)

        return loss

    def get_nearest_neighbor_label(self, batch):
        """
        NOTE: ONLY TO BE USED FOR VALIDATION.

        For each image in validation, find the nearest image in the
        training dataset using the memory bank. Assume its label as
        the predicted label.
        """
        _, img, label = batch
        outputs = self.forward(img)

        all_dps = self.memory_bank.get_all_dot_products(outputs)
        _, neighbor_idxs = torch.topk(all_dps, k=1, sorted=False, dim=1)
        neighbor_idxs = neighbor_idxs.squeeze(1)
        neighbor_idxs = neighbor_idxs.cpu().numpy()

        neighbor_labels = self.train_ordered_labels[neighbor_idxs]
        neighbor_labels = torch.from_numpy(neighbor_labels).long()

        num_correct = torch.sum(neighbor_labels.cpu() == label.cpu()).item()
        return num_correct, img.size(0)

    def training_step(self, batch, batch_idx):
        loss = self.get_losses_for_batch(batch, train=True)
        metrics = {'loss': loss}
        return {'loss': loss, 'log': metrics}

    def validation_step(self, batch, batch_idx):
        loss = self.get_losses_for_batch(batch, train=False)
        num_correct, batch_size = self.get_nearest_neighbor_label(batch)
        output = OrderedDict({'val_loss': loss,
                              'val_num_correct': num_correct,
                              'val_num_total': batch_size})
        return output

    def validation_epoch_end(self, outputs):
        metrics = {}
        for key in outputs[0].keys():
            metrics[key] = torch.tensor([elem[key] for elem in outputs]).float().mean()
        num_correct = sum([out['val_num_correct'] for out in outputs])
        num_total = sum([out['val_num_total'] for out in outputs])
        val_acc = num_correct / float(num_total)
        metrics['val_acc'] = val_acc
        return {'val_loss': metrics['val_loss'], 'log': metrics, 'val_acc': val_acc}

    @pl.data_loader
    def train_dataloader(self):
        return create_dataloader(self.train_dataset, self.config)

    @pl.data_loader
    def val_dataloader(self):
        return create_dataloader(self.val_dataset, self.config, shuffle=False)

class PretrainTwoViewsSystem(PretrainSystem):

    def __init__(self, config):
        super().__init__(config)

        self.loss_name = self.config.loss_params.name
        if self.loss_name == 'MoCo':
            self.model_k = self.create_encoder()

            for param_q, param_k in zip(self.model.parameters(), self.model_k.parameters()):
                param_k.data.copy_(param_q.data)  # initialize
                param_k.requires_grad = False     # do not update

            # create queue (k x out_dim)
            moco_queue = torch.randn(
                self.config.loss_params.k,
                self.config.model_params.out_dim,
            )
            self.register_buffer("moco_queue", moco_queue)
            self.moco_queue = utils.l2_normalize(moco_queue, dim=1)
            self.register_buffer("moco_queue_ptr", torch.zeros(1, dtype=torch.long))

    def get_losses_for_batch(self, batch, train=True):
        indices, img1, img2, _, = batch
        outputs1 = self.forward(img1)

        if self.loss_name == 'SimCLR':
            outputs2 = self.forward(img2)
            loss_fn = SimCLR(outputs1, outputs2,
                             t=self.config.loss_params.t)
            loss = loss_fn.get_loss()
        elif self.loss_name == 'MoCo':
            with torch.no_grad():
                self._momentum_update_key_encoder()
                if self.use_ddp or self.use_ddp2:
                    img2, idx_unshuffle = self._batch_shuffle_ddp(img2)
                outputs2 = self.model_k(img2)
                if self.use_ddp or self.use_ddp2:
                    outputs_k = self._batch_unshuffle_ddp(outputs_k, idx_unshuffle)

            loss_fn = MoCo(outputs1, outputs2,
                           self.moco_queue.clone().detach(),
                           t=self.config.loss_params.t)
            loss = loss_fn.get_loss()

            if train:
                outputs_k = utils.l2_normalize(outputs2, dim=1)
                self._dequeue_and_enqueue(outputs_k)
        else:
            raise Exception(f'Loss {self.loss_name} not supported.')

        if train:
            with torch.no_grad():
                new_data_memory = utils.l2_normalize(outputs1)
                self.memory_bank.update(indices, new_data_memory)

        return loss

    def get_nearest_neighbor_label(self, batch):
        """
        NOTE: ONLY TO BE USED FOR VALIDATION.

        For each image in validation, find the nearest image in the
        training dataset using the memory bank. Assume its label as
        the predicted label.
        """
        _, img, _, label = batch
        outputs = self.forward(img)

        all_dps = self.memory_bank.get_all_dot_products(outputs)
        _, neighbor_idxs = torch.topk(all_dps, k=1, sorted=False, dim=1)
        neighbor_idxs = neighbor_idxs.squeeze(1)
        neighbor_idxs = neighbor_idxs.cpu().numpy()

        neighbor_labels = self.train_ordered_labels[neighbor_idxs]
        neighbor_labels = torch.from_numpy(neighbor_labels).long()

        num_correct = torch.sum(neighbor_labels.cpu() == label.cpu()).item()
        return num_correct, img.size(0)

    @torch.no_grad()
    def _momentum_update_key_encoder(self):
        m = self.config.loss_params.m
        for param_q, param_k in zip(self.model.parameters(), self.model_k.parameters()):
            param_k.data = param_k.data * m + param_q.data * (1. - m)

    @torch.no_grad()
    def _dequeue_and_enqueue(self, keys):
        # gather keys before updating queue
        if self.use_ddp or self.use_ddp2:
            keys = utils.concat_all_gather(keys)

        batch_size = keys.size(0)

        k = self.config.loss_params.k
        ptr = int(self.moco_queue_ptr)
        assert k % batch_size == 0  # why?

        # replace keys at ptr
        self.moco_queue[ptr:ptr+batch_size] = keys
        # move config by full batch size even if current batch is smaller
        ptr = (ptr + batch_size) % k

        self.moco_queue_ptr[0] = ptr

    @torch.no_grad()
    def _batch_shuffle_ddp(self, x):  # pragma: no-cover
        """
        Batch shuffle, for making use of BatchNorm.
        *** Only support DistributedDataParallel (DDP) model. ***
        """
        # gather from all gpus
        batch_size_this = x.shape[0]
        x_gather = concat_all_gather(x)
        batch_size_all = x_gather.shape[0]

        num_gpus = batch_size_all // batch_size_this

        # random shuffle index
        idx_shuffle = torch.randperm(batch_size_all).cuda()

        # broadcast to all gpus
        torch.distributed.broadcast(idx_shuffle, src=0)

        # index for restoring
        idx_unshuffle = torch.argsort(idx_shuffle)

        # shuffled index for this gpu
        gpu_idx = torch.distributed.get_rank()
        idx_this = idx_shuffle.view(num_gpus, -1)[gpu_idx]

        return x_gather[idx_this], idx_unshuffle

    @torch.no_grad()
    def _batch_unshuffle_ddp(self, x, idx_unshuffle):  # pragma: no-cover
        """
        Undo batch shuffle.
        *** Only support DistributedDataParallel (DDP) model. ***
        """
        # gather from all gpus
        batch_size_this = x.shape[0]
        x_gather = concat_all_gather(x)
        batch_size_all = x_gather.shape[0]

        num_gpus = batch_size_all // batch_size_this

        # restored index for this gpu
        gpu_idx = torch.distributed.get_rank()
        idx_this = idx_unshuffle.view(num_gpus, -1)[gpu_idx]

        return x_gather[idx_this]

class TransferSystem(pl.LightningModule):

    def __init__(self, config):
        super().__init__()
        self.config = config
        self.device = f'cuda:{config.gpu_device}' if config.cuda else 'cpu'
        self.train_dataset, self.val_dataset = datasets.get_datasets(
            config.data_params.dataset)
        self.encoder = self.load_pretrained_model()
        utils.frozen_params(self.encoder)
        self.model = self.create_model()

    def load_pretrained_model(self):
        base_dir = self.config.pretrain_model.exp_dir
        checkpoint_name = self.config.pretrain_model.checkpoint_name

        config_path = os.path.join(base_dir, 'config.json')
        config_json = utils.load_json(config_path)
        config = DotMap(config_json)
        # overwrite GPU to load on same as current agent
        config.gpu_device = self.config.gpu_device

        SystemClass = globals()[config.system]
        system = SystemClass(config)
        checkpoint_file = os.path.join(base_dir, 'checkpoints', checkpoint_name)
        checkpoint = torch.load(checkpoint_file)
        system.load_state_dict(checkpoint['state_dict'])

        encoder = system.model.eval()
        return encoder

    def create_model(self):
        dataset_name = self.config.data_params.dataset
        NUM_CLASS_DICT = {'cifar10': 10, 'imagenet': 1000}
        model = LogisticRegression(512*7*7, NUM_CLASS_DICT[dataset_name])
        return model

    def forward(self, img):
        batch_size = img.size(0)
        embs = self.encoder(img, layer=5)
        embs = embs.view(batch_size, -1)
        return self.model(embs)

    def get_losses_for_batch(self, batch):
        _, img, label = batch
        logits = self.forward(img)
        return F.cross_entropy(logits, label)

    def get_accuracies_for_batch(self, batch):
        _, img, label = batch
        logits = self.forward(img)
        preds = torch.argmax(F.log_softmax(logits, dim=1), dim=1)
        preds = preds.long().cpu()
        num_correct = torch.sum(preds == label.long().cpu()).item()
        num_total = img.size(0)
        return num_correct, num_total

    def training_step(self, batch, batch_idx):
        loss = self.get_losses_for_batch(batch)
        with torch.no_grad():
            num_correct, num_total = self.get_accuracies_for_batch(batch)
            metrics = {
                'train_loss': loss,
                'train_num_correct': num_correct,
                'train_num_total': num_total,
                'train_acc': num_correct / float(num_total),
            }
        return {'loss': loss, 'log': metrics}

    def validation_step(self, batch, batch_idx):
        loss = self.get_losses_for_batch(batch)
        num_correct, num_total = self.get_accuracies_for_batch(batch)
        return OrderedDict({
            'val_loss': loss,
            'val_num_correct': num_correct,
            'val_num_total': num_total,
            'val_acc': num_correct / float(num_total)
        })

    def validation_epoch_end(self, outputs):
        metrics = {}
        for key in outputs[0].keys():
            metrics[key] = torch.tensor([elem[key] for elem in outputs]).float().mean()
        num_correct = sum([out['val_num_correct'] for out in outputs])
        num_total = sum([out['val_num_total'] for out in outputs])
        val_acc = num_correct / float(num_total)
        metrics['val_acc'] = val_acc
        return {'val_loss': metrics['val_loss'], 'log': metrics, 'val_acc': val_acc}

    def configure_optimizers(self):
        optim = torch.optim.SGD(
            self.model.parameters(),
            lr=self.config.optim_params.learning_rate,
            momentum=self.config.optim_params.momentum,
            weight_decay=self.config.optim_params.weight_decay,
        )
        return [optim], []

    @pl.data_loader
    def train_dataloader(self):
        return create_dataloader(self.train_dataset, self.config)

    @pl.data_loader
    def val_dataloader(self):
        return create_dataloader(self.val_dataset, self.config, shuffle=False)


text:
## =========================================================================
## @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
## =========================================================================

import os, sys
sys.path.append( os.path.join( os.getcwd( ), '../../lib/python3' ) )

import numpy
import PUJ.Model.Linear, PUJ.Regression.MSE
import PUJ.Optimizer.GradientDescent
import PUJ.Debug.Polynomial

## -------------------------------------------------------------------------
learning_rate = 1e-4
maximum_iterations = 1000
epsilon = 1e-8
debug_step = 10
init_theta = 0

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

## Read data
D = numpy.loadtxt( open( sys.argv[ 1 ], 'rb' ), delimiter = ',' )
numpy.random.shuffle( D )

## -- Separate X and y
x0 = D[ : , : -1 ]
y0 = numpy.matrix( D[ : , -1 ] ).T

# Prepare regression
cost = PUJ.Regression.MSE( x0, y0 )

# Analitical solution
tA = cost.AnalyticSolve( )

# Prepare debug
debug = PUJ.Debug.Polynomial( x0[ : , 0 ], y0 )

# Iterative solution
tI, nI = PUJ.Optimizer.GradientDescent(
  cost,
  learning_rate = learning_rate,
  init_theta = init_theta,
  maximum_iterations = maximum_iterations,
  epsilon = epsilon,
  debug_step = debug_step,
  debug_function = debug
  )

print( '=================================================================' )
print( '* Analitical solution  :', tA )
print( '* Iterative solution   :', tI )
print( '* Difference           :', ( ( tI - tA ) @ ( tI - tA ).T )[ 0, 0 ] ** 0.5 )
print( '* Number of iterations :', nI )
print( '=================================================================' )

debug.KeepFigures( )

## eof - $RCSfile$

max_stars_repo_name: GioLopez/PUJ-Buscaminas

max_issues_repo_name: GioLopez/PUJ-Buscaminas

max_forks_repo_name: GioLopez/PUJ-Buscaminas


text:
import pdb

import numpy as np
from scipy.stats import norm

from bgtools.utils.dates import parse_date
from bgtools.math.solvers import Secant

from pybg.settings import get_eval_date
from pybg.enums import (
    BusinessDayConventions, Calendars, DayCounters, TimeUnits
    )

def nprime(x):
     return np.exp(-x*x/2.0)/np.sqrt(2.*np.pi)

class OptionType:
    CALL = 1
    PUT = -1

class EquityOption(object):
    """
    Equity Option: European exercise
    strike:             strike price
    expiration:         date of expiration
    dividend:           dividend rate
    putcall:            1 = call, -1 = put (see OptionType class)
    settlementDays:     CBOE standard = 1
    calendar:           default UnitedStates(NYSE)

    """
    def __init__(self,
                 strike,
                 expiration,
                 dividend=None,
                 putcall=OptionType.CALL,
                 settlementDays = 1,
                 calendar=Calendars.UnitedStates(Calendars.NYSE)
                 ):

        self._strike = float(strike)
        self._expiration = parse_date(expiration)
        self._dividend = dividend if dividend else 0.0
        self._putcall = putcall
        self._calendar = calendar
        self._settlementDays = settlementDays
        self._settlement = None

    @property
    def settlementDate(self):
        dt = getattr(self, "_settlement", None)
        if not dt:
            dt = Calendars.advance(get_eval_date(),
                                   self.settlementDays,
                                   TimeUnits.Days,
                                   self._calendar,
                                   BusinessDayConventions.Following)
        self._settlement = dt

        return dt

    @property
    def strike(self):
        return self._strike

    @property
    def expiration(self):
        return self._expiration

    @property
    def dividend(self):
        return self._dividend

    @property
    def settlementDays(self):
        return self._settlementDays

    @settlementDays.setter
    def settlementDays(self, ndays):
        self._settlementDays = ndays

    def calc(self, vol, stock_px, short_rate, greeks=True):
        """
        vol, volatility (sigma)
        stock_px, stock price
        short_rate, short rate in continous compounding (=log(1+discount))

        """
        # rename variables to shorten formulas for readability
        v = vol
        S = stock_px
        r = short_rate

        d = self._putcall
        X = self.strike
        q = self.dividend

        cdf = norm.cdf
        exp = np.exp

        T = DayCounters.year_fraction(self.settlementDate,
                                      self.expiration,
                                      DayCounters.ActualActual())

        pv_q = exp(-q * T)

        d1 = (np.log(S/X) + (r - q + v * v/2.) * T) / (v*np.sqrt(T))
        d2 = d1 - v * np.sqrt(T)

        value = d * (S * pv_q * cdf(d * d1) - X * exp(-r * T) * cdf(d * d2))

        if greeks:
            N = nprime(d1)
            sqrtT = np.sqrt(T)
            pvr = exp(-r * T)

            gamma = pv_q * N / (S * v * sqrtT)
            vega = S * sqrtT * N * pv_q

            if d > 0:
                delta = pv_q * cdf(d1)
                theta = S*pv_q*(q*cdf(d1) - N*v/(2.*sqrtT)) - r*X*pvr*cdf(d2)
                rho = X * T * exp(-r*T) * cdf(d2)
            else:
                delta = pv_q * (cdf(d1) - 1)
                theta = -S*pv_q*(q*cdf(-d1) + N*v/(2.*sqrtT)) - r*X*pvr*cdf(-d2)
                rho = -X * T * np.exp(-r*T) * cdf(-d2)

            self.greeks = {
                'value': value,
                'T': T,
                'delta': delta,
                'gamma': gamma,
                'theta': theta,
                'vega': vega,
                'rho':rho
                }

        else:
            self.greeks = {}

        return value

    def implVol(self, option_px, stock_px, short_rate, greeks=True):

        self.greeks = {}
        obj_px = lambda v: self.calc(v, stock_px, short_rate, False)

        iv = Secant(.05, .20, obj_px, option_px)

        if greeks:
            px = self.calc(stock_px, iv, short_rate, True)

        return iv

max_stars_repo_name: bondgeek/pybg

max_stars_repo_licenses: Unlicense MIT

max_issues_repo_name: bondgeek/pybg

max_issues_repo_licenses: Unlicense MIT

max_forks_repo_name: bondgeek/pybg

max_forks_repo_licenses: Unlicense MIT

reason: import numpy,from scipy


text:
from __future__ import print_function

import argparse
import cv2
import glob
import numpy as np
import os
import re
import sys
from fpdf import FPDF
from multiprocessing import cpu_count
from multiprocessing.pool import Pool
from os.path import join, isfile
from subprocess import check_call, check_output

import algorithm
import binarize
import dewarp
from crop import crop
from geometry import Crop
from lib import debug_imwrite
import lib

extension = '.png'
def process_image(original, dpi=None):
    original_rot90 = original

    for i in range(args.rotate // 90):
        original_rot90 = np.rot90(original_rot90)

    # original_rot90 = cv2.resize(original_rot90, (0, 0), None, 1.5, 1.5)
    im_h, im_w = original_rot90.shape[:2]
    # image height should be about 10 inches. round to 100
    if not dpi:
        dpi = int(round(im_h / 1100.0) * 100)
        print('detected dpi:', dpi)

    split = im_w > im_h # two pages

    cropped_images = []
    if args.dewarp:
        lib.debug_prefix.append('dewarp')
        dewarped_images = dewarp.kim2014(original_rot90)
        for im in dewarped_images:
            bw = binarize.binarize(im, algorithm=binarize.sauvola, resize=1.0)
            lib.debug_prefix.append('crop')
            _, [lines] = crop(im, bw, split=False)
            lib.debug_prefix.pop()
            c = Crop.from_lines(lines)
            if c.nonempty():
                cropped_images.append(Crop.from_whitespace(bw).apply(im))
        lib.debug_prefix.pop()
    else:
        bw = binarize.binarize(original_rot90, algorithm=binarize.adaptive_otsu, resize=1.0)
        debug_imwrite('thresholded.png', bw)
        AH, line_sets = crop(original_rot90, bw, split=split)

        for lines in line_sets:
            c = Crop.from_lines(lines)
            if c.nonempty():
                lib.debug = False
                bw_cropped = c.apply(bw)
                orig_cropped = c.apply(original_rot90)
                angle = algorithm.skew_angle(bw_cropped, original_rot90, AH, lines)
                if not np.isfinite(angle): angle = 0.
                rotated = algorithm.safe_rotate(orig_cropped, angle)

                rotated_bw = binarize.binarize(rotated, algorithm=binarize.adaptive_otsu)
                _, [new_lines] = crop(rotated, rotated_bw, split=False)

                # dewarped = algorithm.fine_dewarp(rotated, new_lines)
                # _, [new_lines] = crop(rotated, rotated_bw, split=False)
                new_crop = Crop.union_all([line.crop() for line in new_lines])

                if new_crop.nonempty():
                    # cropped = new_crop.apply(dewarped)
                    cropped = new_crop.apply(rotated)
                    cropped_images.append(cropped)

    out_images = []
    lib.debug_prefix.append('binarize')
    for i, cropped in enumerate(cropped_images):
        lib.debug_prefix.append('page{}'.format(i))
        if lib.is_bw(original_rot90):
            out_images.append(binarize.otsu(cropped))
        else:
            out_images.append(
                binarize.ng2014_fallback(binarize.grayscale(cropped))
            )
        lib.debug_prefix.pop()
    lib.debug_prefix.pop()

    return dpi, out_images

def process_file(file_args):
    (inpath, outdir, dpi) = file_args
    outfiles = glob.glob('{}/{}_*{}'.format(outdir, inpath[:-4], extension))
    if outfiles:
        print('skipping', inpath)
        return outfiles
    else:
        print('processing', inpath)

    original = lib.imread(inpath)
    dpi, out_images = process_image(original, dpi=dpi)
    for idx, outimg in enumerate(out_images):
        outfile = '{}/{}_{}{}'.format(outdir, inpath[:-4], idx, extension)
        print('    writing', outfile)
        cv2.imwrite(outfile, outimg)
        outfiles.append(outfile)

    return outfiles

def pdfimages(pdf_filename):
    assert pdf_filename.endswith('.pdf')
    dirpath = pdf_filename[:-4]
    if not os.path.isdir(dirpath):
        os.makedirs(dirpath)
    if not os.listdir(dirpath):
        check_call(['pdfimages', '-png', pdf_filename, join(dirpath, 'page')])
    return dirpath

def unzip(zip_filename):
    assert zip_filename.endswith('.zip')
    directory = os.path.dirname(zip_filename)
    filenames = check_output(['zipinfo', '-1', zip_filename]).splitlines()
    paths = [os.path.join(directory, f) for f in filenames]
    if not any([os.path.isfile(f) for f in paths]):
        check_call(['unzip', zip_filename])
    return [f for f in paths if os.path.isfile(f)]

def sorted_numeric(strings):
    return sorted(strings, key=lambda f: list(map(int, re.findall('[0-9]+', f))))

def accumulate_paths(target, accum):
    for path in target:
        assert os.path.exists(path)
        if os.path.isfile(path):
            if path.endswith('.pdf'):
                accumulate_paths([pdfimages(path)], accum)
            elif path.endswith('.zip'):
                accumulate_paths(unzip(path), accum)
            elif re.match(r'.*\.(png|jpg|tif|dng)', path):
                accum.append(path)
        else:
            assert os.path.isdir(path)
            files = [os.path.join(path, base) for base in sorted_numeric(os.listdir(path))]
            accumulate_paths(files, accum)

def run(args):
    if args.single_file:
        lib.debug = True
        im = lib.imread(args.single_file)
        _, out_images = process_image(im, dpi=args.dpi)
        for idx, outimg in enumerate(out_images):
            cv2.imwrite('out{}.png'.format(idx), outimg)
        return

    if args.concurrent:
        pool = Pool(cpu_count())
        map_fn = pool.map
    else:
        map_fn = map

    files = []
    accumulate_paths(args.indirs, files)
    files = sorted_numeric(list(set(files)))
    print('Files:', files)

    for p in files:
        d = os.path.dirname(p)
        if not os.path.isdir(join(args.outdir, d)):
            os.makedirs(join(args.outdir, d))

    outfiles = map_fn(process_file, list(zip(files,
                        [args.outdir] * len(files),
                        [args.dpi] * len(files))))

    outfiles = sum(outfiles, [])
    outfiles.sort(key=lambda f: list(map(int, re.findall('[0-9]+', f))))

    # outtif = join(args.outdir, 'out.tif')
    outpdfpath = join(args.outdir, 'out.pdf')
    if not isfile(outpdfpath):
        print('making pdf:', outpdfpath)
        pdf = FPDF(unit='in', format='Letter')
        pdf.set_margins(0, 0, 0)
        for outfile in outfiles:
            print('.', end='')
            sys.stdout.flush()
            pdf.add_page()
            im = cv2.imread(outfile)
            inches_w = float(im.shape[1]) / args.dpi
            inches_h = float(im.shape[0]) / args.dpi
            x = (8.5 - inches_w) / 2
            y = (11.0 - inches_h) / 2
            pdf.image(outfile, w=inches_w, x=x, y=y)
        print()

        pdf.output(name=outpdfpath)

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Batch-process for PDF')
    parser.add_argument('outdir', nargs='?', help="Output directory")
    parser.add_argument('indirs', nargs='+', help="Input directory")
    parser.add_argument('-f', '--file', dest='single_file', action='store',
                        help="Run on single file instead")
    parser.add_argument('-c', '--concurrent', action='store_true',
                        help="Run w/ threads.")
    parser.add_argument('-d', '--dpi', action='store', type=int,
                        help="Force a particular DPI")
    parser.add_argument('--dewarp', action='store_true', help="Dewarp pages.")
    parser.add_argument('--rotate', action='store', type=int, choices=[0, 90, 180, 270],
                        default=0, help="Rotate CCW by 90, 180, or 270 degrees.")

    global args
    args = parser.parse_args()
    run(args)

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
#
# Reference : https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3864968/
#
# Written by SK Park , NIRAL, UNC
# 2021-04-18

import dmri.preprocessing as prep

import yaml,traceback
import BASELINE_Average.computations as computations
import numpy as np
import time
from pathlib import Path

logger=prep.logger.write

class BASELINE_Average(prep.modules.DTIPrepModule):
    def __init__(self,config_dir,*args,**kwargs):
        super().__init__(config_dir)

    def generateDefaultProtocol(self,image_obj):
        super().generateDefaultProtocol(image_obj)
        ## todos
        return self.protocol

    def process(self,*args,**kwargs): ## variables : self.source_image, self.image (output) , self.result_history , self.result (output) , self.protocol, self.template
        super().process()
        inputParams=self.getPreviousResult()['output']
        opts=args[0]
        self.baseline_threshold=opts['baseline_threshold']

        output=None
        output_image_path=Path(self.output_dir)
        if self.image.image_type.lower()=='nrrd':
            output_image_path=str(output_image_path.joinpath('output.nrrd'))
        elif self.image.image_type.lower()=='nifti':
            output_image_path=str(output_image_path.joinpath('output.nii.gz'))

        output_filename=Path(self.computation_dir).joinpath('computations.yml')
        if output_filename.exists() and not self.options['recompute']:
            ## pass recomputation
            logger("Computing ommited",prep.Color.INFO)
            pass
        else: ##computed parameters doesn't exist or recompute is true
            ## compute or recompute
            logger("Computing ... ",prep.Color.PROCESS)
            #self.image.deleteGradientsByOriginalIndex([49, 65, 97, 129, 145])#([0, 17, 49, 65, 97, 129, 145]) For test
            new_image, excluded_original_indexes=computations.baseline_average(self.image, opt=None ,
                                                      averageInterpolationMethod=self.protocol['averageInterpolationMethod'],
                                                      averageMethod=self.protocol['averageMethod'],
                                                      b0Threshold=self.baseline_threshold,
                                                      stopThreshold=self.protocol['stopThreshold'],
                                                      maxIterations=self.protocol['maxIterations'])

        if new_image is not None:
            self.image=new_image
            ### if image is changed, next module should load the file. So set image_object to None and write the file instead
        self.image.setSpaceDirection(target_space=self.getSourceImageInformation()['space'])
        self.writeImage(output_image_path,dest_type=self.image.image_type)

        ### output preparation
        self.result['output']['excluded_gradients_original_indexes']=excluded_original_indexes
        self.result['output']['success']=True
        #raise Exception("User Exception for development ...")
        return self.result

max_stars_repo_name: NIRALUser/DTIPlayground

max_issues_repo_name: NIRALUser/DTIPlayground

max_forks_repo_name: NIRALUser/DTIPlayground


text:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize

def Att(f, ft, a):
    return a/((1+(f/ft)**2)**0.5)

def derivata(f, ft, a):
    return -a*(f/ft)*((1+(f/ft)**2)**(-1.5))

def Bode(f, ft):
    return 20*(np.log10(ft/f))

def derivataBode(f, ft):
    return np.abs(-20/(np.log(10)*f))

Vout, dVout, f, df = np.genfromtxt('data.txt', skip_header = 1, unpack = True)
Vi = 6.08
dVi = 0.19
A = Vout/Vi
dA = (1/Vi)*dVout + (Vout/(Vi**2))*dVi

# Fit normale
dA_eff = dA
popt = (174,1) # parametri iniziali
for i in range(5):
    popt, pcov = scipy.optimize.curve_fit(Att, f, A, popt, dA_eff, absolute_sigma = False)
    chi_2 = np.sum(((A - Att(f,*popt))/dA_eff)**2)
    print(chi_2)
    dA_eff = np.sqrt(((derivata(f,*popt))*df)**2 + dA**2)

ndof = len(Vout)-2
print('\nil chi quadro è', chi_2)
print('il chi quadro ridotto è', chi_2/ndof)
print('ft (frequenza di taglio) =',popt[0], '+-', pcov[0][0]**0.5)
print('a =',popt[1], '+-', pcov[1][1]**0.5)
print("Cov normalizzata =", pcov[1][0]/(pcov[0][0]*pcov[1][1])**0.5, "\n")

# Fit Bode
mask = f>400
ABode = 20*np.log10(A[mask])
fBode = f[mask]
dABode = 20*(dA[mask]/A[mask])/np.log(10)
dA_effB = dABode
dfBode = df[mask]
for i in range(5):
    poptBode, pcovBode = scipy.optimize.curve_fit(Bode, fBode, ABode, (1), dA_effB, absolute_sigma = False)
    chi_2 = np.sum(((ABode - Bode(fBode,*poptBode))/dA_effB)**2)
    print(chi_2)
    dA_effB = np.sqrt(((derivataBode(fBode, *poptBode))*dfBode)**2 + dABode**2)

print("ft =",poptBode, "+-", np.diag(pcovBode)**0.5, "\n")

# Grafici
# Plot funzione
x = np.logspace(1, 5, 4000)
plt.figure()
plt.subplot(211)
plt.xscale('log')
plt.yscale('log')
plt.errorbar(f, A, yerr=dA, xerr=df, fmt='.', label="Misure")
plt.plot(x, Att(x, *popt), label="Fit")
plt.xlabel("f [Hz]")
plt.ylabel("A")
plt.legend()

# Plot residui normalizzati
plt.subplot(212)
plt.xscale('log')
plt.errorbar(f, (A -Att(f, *popt))/dA_eff, fmt='.', label="Res. norm.")
plt.plot(x, x*0)
plt.xlabel("f [Hz]")
plt.legend()
plt.show()

# Bode plot
xB = np.linspace(0, 8, 2000)
plt.figure()
plt.subplot(211)
plt.xscale('log')
plt.errorbar(f, 20*np.log10(A), yerr=20*(dA/A)/np.log(10), xerr=df, fmt='.', label="Misure")
plt.plot(x, Bode(x, *poptBode), label="Fit")
plt.xlabel("f [Hz]")
plt.ylabel("A")
plt.legend()

# Bode residui normalizzati
plt.subplot(212)
plt.xscale('log')
plt.errorbar(fBode, (ABode - Bode(fBode, *poptBode))/dA_effB, fmt='.', label="Res. norm.")
plt.plot(x, x*0)
plt.xlabel("f [Hz]")
plt.legend()
plt.show()

max_stars_repo_path: Lab 2/7 filtri RC/filtriRC.py

max_issues_repo_path: Lab 2/7 filtri RC/filtriRC.py

max_forks_repo_path: Lab 2/7 filtri RC/filtriRC.py

reason: import numpy,import scipy