text:
import os
import sys

import numpy as np

from GenNet_utils.hase.config import MAPPER_CHUNK_SIZE, basedir, CONVERTER_SPLIT_SIZE, PYTHON_PATH

os.environ['HASEDIR'] = basedir
if PYTHON_PATH is not None:
    for i in PYTHON_PATH: sys.path.insert(0, i)
from GenNet_utils.hase.hdgwas.tools import Timer, Checker, study_indexes, Mapper, HaseAnalyser, merge_genotype, \
    check_np, check_converter
from GenNet_utils.hase.hdgwas.converter import GenotypePLINK, GenotypeMINIMAC, GenotypeVCF
from GenNet_utils.hase.hdgwas.data import Reader, MetaParData, MetaPhenotype
from GenNet_utils.hase.hdgwas.fake import Encoder
from GenNet_utils.hase.hdgwas.hdregression import HASE, A_inverse, B4
import argparse
import gc
from GenNet_utils.hase.hdgwas.pard import partial_derivatives
from GenNet_utils.hase.hdgwas.regression import haseregression
import pandas as pd
import time
from GenNet_utils.hase.hdgwas.protocol import Protocol

__version__ = '1.1.0'

HEAD = "*********************************************************************\n"
HEAD += "* HASE: Framework for efficient high-dimensional association analyses \n"
HEAD += "* Version {V}\n".format(V=__version__)
HEAD += "* (C) 2015-2017 Gennady Roshchupkin and Hieab Adams\n"
HEAD += "* Erasmus MC, Rotterdam /  Department of Medical Informatics, Radiology and Epidemiology \n"
HEAD += "* GNU General Public License v3\n"
HEAD += "*********************************************************************\n"

if __name__ == '__main__':

    start = time.time()

    parser = argparse.ArgumentParser(description='Script to use HASE in command line')

    # MAIN SETTINGS
    parser.add_argument("-thr", type=float,
                        help="predefined threshold for t-value, if not defined save everything")

    parser.add_argument("-o", "--out", type=str, required=True, help="path to save result folder")

    parser.add_argument("-mode", required=True, type=str,
                        choices=['regression', 'converting', 'single-meta', 'meta-stage', 'encoding'],
                        help='Main argument, specify type of analysis'
                             '*****************'
                             'converting - before start to use any kind of analysis you need to convert data, it could take some time, '
                             'you need to do it ones and then everything will be super fast'
                             '*****************'
                             'single-meta - if you are the part of meta-analysis and will send data to main center, use this mode'
                             '*****************'
                             'encoding - if you are doing HD association, before send data to main center create fake data'
                             '*****************'
                             'meta-stage - if you are the main center of meta-analysis and already got precomputed data from all sites, use it'
                             '*****************'
                             'regression - if you want run everything yourself and do not use precomputed data use it'

                             'If you want to run single site full association analysis run (1) single-meta mode and then (2) meta-stage'

                             '*****************'
                        )

    parser.add_argument("-g", "--genotype", nargs='+', type=str, help="path/paths to genotype data folder")
    parser.add_argument("-ph", "--phenotype", nargs='+', type=str, help="path to phenotype data folder")
    parser.add_argument("-cov", "--covariates", type=str, help="path to covariates data folder")

    parser.add_argument("-derivatives", nargs='+', type=str, help="path to derivatives data folder")

    parser.add_argument("-protocol", type=str, help="path to study protocol file")

    parser.add_argument('-study_name', type=str, required=True, nargs='+',
                        help=' Name for saved genotype data, without ext')

    parser.add_argument('-mapper', type=str, help='Mapper data folder')

    parser.add_argument('-ref_name', type=str, default='1000Gp1v3_ref', help='Reference panel name')

    # ADDITIONAL SETTINGS
    parser.add_argument("-snp_id_inc", type=str, help="path to file with SNPs id to include to analysis")
    parser.add_argument("-snp_id_exc", type=str, help="path to file with SNPs id to exclude from analysis")
    parser.add_argument("-ph_id_inc", type=str, help="path to file with phenotype id to exclude from analysis")
    parser.add_argument("-ph_id_exc", type=str, help="path to file with phenotype id to exclude from analysis")

    # parser.add_argument("-ind_id_inc", type=str, help="path to file with individuals id to include to analysis") #TODO (low)
    # parser.add_argument("-ind_id_exc", type=str, help="path to file with individuals id to exclude from analysis")#TODO (low)
    # parser.add_argument("-cov_name_inc", type=str, help="path to file with covariates names to include to analysis")#TODO (low)
    # parser.add_argument("-cov_name_exc", type=str, help="path to file with covariates names to exclude from analysis")#TODO (low)

    parser.add_argument('-intercept', type=str, default='y', choices=['y', 'n'],
                        help='include intercept to regression, default yes')

    parser.add_argument('-maf', type=float, default=0.0, help='MAF for genetics data')

    parser.add_argument('-encoded', nargs='+', type=int, help='Value per every study, 1 - if encoded, 0 - if not')
    ###

    # FLAGS
    parser.add_argument('-hdf5', action='store_true', default=True, help='flag for genotype data format')
    parser.add_argument('-id', action='store_true', default=False,
                        help='Flag to convert minimac data to genotype per subject files first (default False)')
    parser.add_argument('-pd_full', action='store_true', default=False, help='For not HD association study')
    parser.add_argument('-effect_intercept', action='store_true', default=False,
                        help='Flag for add study effect to PD regression model')
    parser.add_argument('-permute_ph', action='store_true', default=False, help='Flag for phenotype permutation')
    parser.add_argument('-vcf', action='store_true', default=False, help='Flag for VCF data to convert')
    parser.add_argument('-np', action='store_true', default=True, help='Check BLAS/LAPACK/MKL')

    # TODO (low) save genotype after MAF
    ###

    ###CLUSTER SETTING
    parser.add_argument('-cluster', type=str, default='n', choices=['y', 'n'],
                        help=' Is it parallel cluster job, default no')
    parser.add_argument('-node', nargs='+', type=int, help='number of nodes / this node number, example: 10 2 ')
    ###

    # ADVANCED SETTINGS
    parser.add_argument('-interaction', type=str,
                        help='path to file with data for genotype interaction test')  # TODO (low)
    parser.add_argument('-mapper_chunk', type=int, help='Change mapper chunk size from config file')
    ###
    args = parser.parse_args()
    if not args.thr:
        print('WARNING!!! You did not set threshold for t-value, all results will be saved')
    if args.mapper_chunk:
        MAPPER_CHUNK_SIZE = args.mapper_chunk
    ARG_CHECKER = Checker()
    print(args)
    os.environ['HASEOUT'] = args.out

    if args.cluster == 'y':
        if args.node is not None:
            if args.node[1] > args.node[0]:
                raise ValueError('Node # {} > {} total number of nodes'.format(args.node[1], args.node[0]))

    if not os.path.isdir(args.out):
        print("Creating output folder {}".format(args.out))
        os.mkdir(args.out)

    if args.np:
        check_np()

    ################################### CONVERTING ##############################
    if args.mode == 'converting':

        # ARG_CHECKER.check(args,mode='converting')

        R = Reader('genotype')
        R.start(args.genotype[0], vcf=args.vcf)

        with Timer() as t:
            if R.format == 'PLINK':
                G = GenotypePLINK(args.study_name[0], reader=R)
                G.split_size = CONVERTER_SPLIT_SIZE
                G.plink2hdf5(out=args.out)

            elif R.format == 'MINIMAC':
                G = GenotypeMINIMAC(args.study_name[0], reader=R)
                if args.cluster == 'y':
                    G.cluster = True
                G.split_size = CONVERTER_SPLIT_SIZE
                G.MACH2hdf5(args.out, id=args.id)

            elif R.format == 'VCF':
                G = GenotypeVCF(args.study_name[0], reader=R)
                if args.cluster == 'y':
                    G.cluster = True
                G.split_size = CONVERTER_SPLIT_SIZE
                G.VCF2hdf5(args.out)
            else:
                raise ValueError('Genotype data should be in PLINK/MINIMAC/VCF format and alone in folder')

        check_converter(args.out, args.study_name[0])
        print(('Time to convert all data: {} sec'.format(t.secs)))

    ################################### ENCODING ##############################

    elif args.mode == 'encoding':

        # ARG_CHECKER.check(args,mode='encoding')
        mapper = Mapper()
        mapper.genotype_names = args.study_name
        mapper.chunk_size = MAPPER_CHUNK_SIZE
        mapper.reference_name = args.ref_name
        mapper.load_flip(args.mapper)
        mapper.load(args.mapper)

        phen = Reader('phenotype')
        phen.start(args.phenotype[0])

        gen = Reader('genotype')
        gen.start(args.genotype[0], hdf5=args.hdf5, study_name=args.study_name[0], ID=False)

        e = Encoder(args.out)
        e.study_name = args.study_name[0]

        row_index, ids = study_indexes(phenotype=phen.folder._data, genotype=gen.folder._data)
        with Timer() as t:

            e.matrix(len(ids), save=True)
            N_snps_read = 0
            while True:
                with Timer() as t_gen:
                    genotype = gen.get_next()
                    if isinstance(genotype, type(None)):
                        break

                    flip = mapper.flip[args.study_name[0]][N_snps_read:N_snps_read + genotype.shape[0]]
                    N_snps_read += genotype.shape[0]
                    flip_index = (flip == -1)
                    genotype = np.apply_along_axis(lambda x: flip * (x - 2 * flip_index), 0, genotype)
                    genotype = genotype[:, row_index[0]]
                    encode_genotype = e.encode(genotype, data_type='genotype')
                    e.save_hdf5(encode_genotype, os.path.join(args.out, 'encode_genotype'), info=gen.folder,
                                index=row_index[0])
                    encode_genotype = None
                    gc.collect()

                print(('Time to create fake genotype is {}sec'.format(t_gen.secs)))

            while True:
                with Timer() as t_phen:
                    phenotype = phen.get_next(index=row_index[1])
                    if isinstance(phenotype, type(None)):
                        break
                    encode_phenotype = e.encode(phenotype, data_type='phenotype')
                    if phen.folder.format == '.npy':
                        e.save_npy(encode_phenotype, save_path=os.path.join(args.out, 'encode_phenotype'),
                                   info=phen.folder, index=row_index[1])
                    if phen.folder.format in ['.csv', '.txt']:
                        e.save_csv(encode_phenotype, save_path=os.path.join(args.out, 'encode_phenotype'),
                                   info=phen.folder, index=row_index[1])
                    encode_phenotype = None
                    gc.collect()
            if phen.folder.format == '.npy':
                np.save(os.path.join(os.path.join(args.out, 'encode_phenotype', 'info_dic.npy')), e.phen_info_dic)
            print(('Time to create fake phenotype is {}sec'.format(t_phen.secs)))

        print(('Time to encode all data: {} sec'.format(t.secs)))

    ################################### SINGLE META STAGE ##############################

    elif args.mode == 'single-meta':

        # ARG_CHECKER.check(args,mode='single-meta')
        mapper = Mapper()
        mapper.genotype_names = args.study_name
        mapper.chunk_size = MAPPER_CHUNK_SIZE
        mapper.reference_name = args.ref_name
        mapper.load_flip(args.mapper)
        mapper.load(args.mapper)
        mapper.cluster = args.cluster
        mapper.node = args.node

        phen = Reader('phenotype')
        phen.start(args.phenotype[0])

        cov = Reader('covariates')
        cov.start(args.covariates)

        if cov.folder.n_files > 1:
            raise ValueError('In covariates folder should be only one file!')

        gen = Reader('genotype')
        gen.start(args.genotype[0], hdf5=args.hdf5, study_name=args.study_name[0], ID=False)

        with Timer() as t:
            partial_derivatives(save_path=args.out, COV=cov, PHEN=phen, GEN=gen,
                                MAP=mapper, MAF=args.maf, R2=None, B4_flag=args.pd_full,
                                study_name=args.study_name[0], intercept=args.intercept)
        print(('Time to compute partial derivatives : {} sec'.format(t.secs)))

    ################################### MULTI META STAGE ##############################

    elif args.mode == 'meta-stage':

        # ARG_CHECKER.check(args,mode='meta-stage')

        ##### Init data readers #####
        if args.derivatives is None:
            raise ValueError('For meta-stage analysis partial derivatives data are required!')
        mapper = Mapper()
        mapper.chunk_size = MAPPER_CHUNK_SIZE
        mapper.genotype_names = args.study_name
        mapper.reference_name = args.ref_name
        if args.snp_id_inc is not None:
            mapper.include = pd.DataFrame.from_csv(args.snp_id_inc, index_col=None)
            print('Include:')
            print(mapper.include.head())
            if 'ID' not in mapper.include.columns and (
                    'CHR' not in mapper.include.columns or 'bp' not in mapper.include.columns):
                raise ValueError('{} table does not have ID or CHR,bp columns'.format(args.snp_id_inc))
        if args.snp_id_exc is not None:
            mapper.exclude = pd.DataFrame.from_csv(args.snp_id_exc, index_col=None)
            print('Exclude:')
            print(mapper.exclude.head())
            if 'ID' not in mapper.exclude.columns and (
                    'CHR' not in mapper.exclude.columns or 'bp' not in mapper.exclude.columns):
                raise ValueError('{} table does not have ID or CHR,bp columns'.format(args.snp_id_exc))
        mapper.load(args.mapper)
        mapper.load_flip(args.mapper, encode=args.encoded)
        mapper.cluster = args.cluster
        mapper.node = args.node

        Analyser = HaseAnalyser()

        pard = []

        with Timer() as t:
            for i, j in enumerate(args.derivatives):
                pard.append(Reader('partial'))
                pard[i].start(j, study_name=args.study_name[i])
                pard[i].folder.load()

        print("time to set PD is {}s".format(t.secs))

        PD = [False if isinstance(i.folder._data.b4, type(None)) else True for i in pard]

        if np.sum(PD) != len(pard) and np.sum(PD) != 0:
            raise ValueError('All study should have b4 data for partial derivatives!')

        if args.protocol is not None:
            protocol = Protocol(args.protocol)
        else:
            protocol = None

        meta_pard = MetaParData(pard, args.study_name, protocol=protocol)

        if np.sum(PD) == 0:
            phen = []

            with Timer() as t:
                for i, j in enumerate(args.phenotype):
                    phen.append(Reader('phenotype'))
                    phen[i].start(j)
            print("Time to set pheno {} s".format(t.secs))
            meta_phen = MetaPhenotype(phen, include=args.ph_id_inc, exclude=args.ph_id_exc)

            N_studies = len(args.genotype)

            gen = []
            with Timer() as t:
                for i, j in enumerate(args.genotype):
                    gen.append(Reader('genotype'))
                    gen[i].start(j, hdf5=args.hdf5, study_name=args.study_name[i], ID=False)
            print("Time to set gen {}s".format(t.secs))

            row_index, ids = study_indexes(phenotype=tuple(i.folder._data for i in phen),
                                           genotype=tuple(i.folder._data for i in gen),
                                           covariates=tuple(i.folder._data.metadata for i in pard))
            if row_index[2].shape[0] != np.sum([i.folder._data.metadata['id'].shape[0] for i in pard]):
                raise ValueError(
                    'Partial Derivatives covariates have different number of subjects {} than genotype and phenotype {}'.format(
                        row_index[2].shape[0],
                        np.sum([i.folder._data.metadata['id'].shape[0] for i in pard])))
        while True:
            if mapper.cluster == 'n':
                SNPs_index, keys = mapper.get()
            else:
                ch = mapper.chunk_pop()
                if ch is None:
                    SNPs_index = None
                    break
                SNPs_index, keys = mapper.get(chunk_number=ch)

            if isinstance(SNPs_index, type(None)):
                break

            Analyser.rsid = keys
            if np.sum(PD) == 0:
                genotype = np.array([])
                with Timer() as t_g:
                    genotype = merge_genotype(gen, SNPs_index, mapper)
                    genotype = genotype[:, row_index[0]]
                print("Time to get G {}s".format(t_g.secs))
            # TODO (low) add interaction

            a_test = np.array([])
            b_cov = np.array([])
            C = np.array([])
            a_cov = np.array([])
            b4 = np.array([])

            if args.protocol is not None:
                if protocol.enable:
                    regression_model = protocol.regression_model()
            else:
                regression_model = None

            with Timer() as t_pd:
                if np.sum(PD) == 0:
                    a_test, b_cov, C, a_cov = meta_pard.get(SNPs_index=SNPs_index, regression_model=regression_model,
                                                            random_effect_intercept=args.effect_intercept)
                else:
                    a_test, b_cov, C, a_cov, b4 = meta_pard.get(SNPs_index=SNPs_index, B4=True,
                                                                regression_model=regression_model,
                                                                random_effect_intercept=args.effect_intercept)

            print("Time to get PD {}s".format(t_pd.secs))

            MAF = meta_pard.maf_pard(SNPs_index=SNPs_index)

            if args.maf != 0:
                filter = (MAF > args.maf) & (MAF < 1 - args.maf) & (MAF != 0.5)
                a_test = a_test[filter, :]
                Analyser.MAF = MAF[filter]
                if np.sum(PD) == 0:
                    genotype = genotype[filter, :]
                else:
                    b4 = b4[filter, :]
                Analyser.rsid = Analyser.rsid[filter]
                if a_test.shape[0] == 0:
                    print('NO SNPs > MAF')
                    continue
            else:
                Analyser.MAF = MAF

            a_inv = A_inverse(a_cov, a_test)
            N_con = a_inv.shape[1] - 1
            print('There are {} subjects in study.'.format(meta_pard.get_n_id()))
            DF = (meta_pard.get_n_id() - a_inv.shape[1])

            if np.sum(PD) == 0:

                while True:
                    phenotype = np.array([])
                    with Timer() as t_ph:
                        phenotype, phen_names = meta_phen.get()
                    print("Time to get PH {}s".format(t_ph.secs))

                    if isinstance(phenotype, type(None)):
                        meta_phen.processed = 0
                        print('All phenotypes processed!')
                        break
                    print(("Merged phenotype shape {}".format(phenotype.shape)))
                    # TODO (middle) select phen from protocol
                    phenotype = phenotype[row_index[1], :]
                    print(("Selected phenotype shape {}".format(phenotype.shape)))
                    keys = list(meta_pard.phen_mapper.dic.keys())
                    phen_ind_dic = {k: i for i, k in enumerate(keys)}
                    phen_ind = np.array([phen_ind_dic.get(i, -1) for i in phen_names])
                    if np.sum(phen_ind == -1) == len(phen_ind):
                        print('There is no common ids in phenotype files and PD data!')
                        break
                    else:
                        print('There are {} common ids in phenotype files and PD data!'.format(np.sum(phen_ind != -1)))
                    C_test = C[phen_ind]
                    b_cov_test = b_cov[:, phen_ind]

                    b4 = B4(phenotype, genotype)
                    print(("B4 shape is {}".format(b4.shape)))
                    t_stat, SE = HASE(b4, a_inv, b_cov_test, C_test, N_con, DF)

                    if mapper.cluster == 'y':
                        Analyser.cluster = True
                        Analyser.chunk = ch
                        Analyser.node = mapper.node[1]
                    Analyser.t_stat = t_stat
                    Analyser.SE = SE
                    Analyser.threshold = args.thr
                    Analyser.out = args.out
                    Analyser.save_result(phen_names[(phen_ind != -1)])

                    t_stat = None
                    Analyser.t_stat = None
                    gc.collect()

            else:
                t_stat, SE = HASE(b4, a_inv, b_cov, C, N_con, DF)

                Analyser.t_stat = t_stat
                Analyser.SE = SE
                Analyser.threshold = args.thr
                Analyser.out = args.out
                Analyser.save_result(np.array(list(meta_pard.phen_mapper.dic.keys())))

                t_stat = None
                Analyser.t_stat = None
                gc.collect()

    ################################### TO DO EVERYTHING IN ONE GO ##############################

    elif args.mode == 'regression':

        # ARG_CHECKER.check(args,mode='regression')

        print('START regression mode...')
        if args.mapper is None:
            os.environ['MAPPER'] = "False"  # need to set it here, before genotype Reader start

        phen = Reader('phenotype')
        phen.start(args.phenotype[0])
        phen.permutation = args.permute_ph

        cov = Reader('covariates')
        cov.start(args.covariates)

        interaction = None
        if args.interaction is not None:
            interaction = Reader('interaction')
            interaction.start(args.interaction)

        if (cov.folder.n_files > 1 and cov.folder.format != '.npy') or (
                cov.folder.n_files > 2 and cov.folder.format == '.npy'):  # TODO (middle) test
            raise ValueError('In covariates folder should be only one file!')

        gen = []
        for i, j in enumerate(args.genotype):
            gen.append(Reader('genotype'))
            gen[i].start(j, hdf5=args.hdf5, study_name=args.study_name[i], ID=False)

        if args.mapper is not None:
            mapper = Mapper()
            mapper.chunk_size = MAPPER_CHUNK_SIZE
            mapper.genotype_names = args.study_name
            mapper.reference_name = args.ref_name
            if args.snp_id_inc is not None:
                mapper.include = pd.DataFrame.from_csv(args.snp_id_inc, index_col=None)
                print('Include:')
                print(mapper.include.head())
                if 'ID' not in mapper.include.columns and (
                        'CHR' not in mapper.include.columns or 'bp' not in mapper.include.columns):
                    raise ValueError('{} table does not have ID or CHR,bp columns'.format(args.snp_id_inc))
            if args.snp_id_exc is not None:
                mapper.exclude = pd.DataFrame.from_csv(args.snp_id_exc, index_col=None)
                print('Exclude:')
                print(mapper.exclude.head())
                if 'ID' not in mapper.exclude.columns and (
                        'CHR' not in mapper.exclude.columns or 'bp' not in mapper.exclude.columns):
                    raise ValueError('{} table does not have ID or CHR,bp columns'.format(args.snp_id_exc))
            mapper.load(args.mapper)
            mapper.load_flip(args.mapper)
            mapper.cluster = args.cluster
            mapper.node = args.node
        else:
            if len(args.genotype) == 1:
                mapper = Mapper()
                mapper.chunk_size = MAPPER_CHUNK_SIZE
                mapper.genotype_names = args.study_name
                mapper.reference_name = args.ref_name
                mapper.cluster = args.cluster
                mapper.node = args.node
                mapper.n_study = 1
                mapper.n_keys = gen[0].folder._data.names.shape[0]
                mapper.keys = np.array(gen[0].folder._data.names.tolist())
                mapper.values = np.array(list(range(mapper.n_keys))).reshape(-1, 1)
                mapper.flip[args.study_name[0]] = np.array([1] * mapper.n_keys)
                if args.snp_id_exc is not None or args.snp_id_inc is not None:
                    raise ValueError('You can not exclude or include variants to analysis without mapper!')
            else:
                raise ValueError('You can not run regression analysis with several genotype data without mapper!')
        # mapper=None

        Analyser = HaseAnalyser()
        Analyser.threshold = args.thr
        Analyser.out = args.out
        haseregression(phen, gen, cov, mapper, Analyser, args.maf, intercept=args.intercept, interaction=interaction)

end = time.time()

print(('experiment finished in {} s'.format((end - start))))

max_stars_repo_name: reneevdw/GenNet

max_issues_repo_name: reneevdw/GenNet

max_forks_repo_name: reneevdw/GenNet


text:
import scipy.stats as ss
import numpy as np
from itertools import combinations_with_replacement
import pandas as pd
from scipy.sparse.csgraph import reverse_cuthill_mckee
from scipy.sparse import csr_matrix

def cramers_corrected_stat(confusion_matrix):
    """
    Calculate Cramers V statistic for categorial-categorial association.
    Uses correction from Bergsma and Wicher,
    Journal of the Korean Statistical Society 42 (2013): 323-328

    See:
    https://stackoverflow.com/questions/20892799/
    using-pandas-calculate-cram%C3%A9rs-coefficient-matrix

    Parameters
    --------------------------
    confusion_matrix : numpy confusion matrix

    Returns
    ---------------------------
    float : Cramer's V statistic with bias correction
    """
    chi2 = ss.chi2_contingency(confusion_matrix)[0]
    n = confusion_matrix.sum()
    phi2 = chi2/n
    r,k = confusion_matrix.shape
    phi2corr = max(0, phi2 - ((k-1)*(r-1))/(n-1))
    rcorr = r - ((r-1)**2)/(n-1)
    kcorr = k - ((k-1)**2)/(n-1)
    return(np.sqrt(phi2corr / min( (kcorr-1), (rcorr-1))))

def cramers_corrected_matrix(df, reorder_cuthill_mckee = True):
    """
    Calculate Cramers V statistic with bias correction for all
    combinations of columns in pandas DataFrame df

    Parameters
    --------------------------
    df : pandas DataFrame - all columns must be categorical

    reorder_cuthill_mckee : boolean - whether to reorder to the columns
        based on the reverse Cuthill McKee algorithm applied to the
        matrix of Cramers V statistics

    Returns
    ---------------------------
    Z : numpy array with Cramers V statistics
    """
    cols = df.columns.tolist()
    Z = np.zeros((len(cols),len(cols)))

    for i,j in combinations_with_replacement(range(len(cols)),2):
        z = cramers_corrected_stat(
                pd.crosstab(df[cols[i]],df[cols[j]]).values
            )
        Z[i,j] = Z[j,i] = z

    if reorder_cuthill_mckee is True:
        perm = reverse_cuthill_mckee(
            csr_matrix(Z),
            symmetric_mode = True
        )
        cols = [cols[i] for i in perm]
        gx,gy = np.meshgrid(perm,perm)
        Z = (csr_matrix(Z)[gx,gy]).toarray()

    Z = pd.DataFrame(
            Z,
            index = cols,
            columns = cols
        )

    return(Z)

def cramers_v(df,x,y):
    if isinstance(x,str): x = np.array([x])
    elif isinstance(x,list): x = np.array(x)

    @np.vectorize
    def _cramers_v(x):
        z = cramers_corrected_stat(
                pd.crosstab(df[x],df[y]).values)
        return(z)

    return(_cramers_v(x))

ext: py

max_stars_repo_name: Strabes/helpers

max_issues_repo_name: Strabes/helpers

max_forks_repo_name: Strabes/strappy

reason: import numpy,import scipy,from scipy


text:
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from scipy.optimize import linprog
from scipy.optimize import minimize

x = np.sort(np.random.rand(10))
y = 1 + 2*x + np.random.rand(x.size) - 0.5
 # 添加异常值
x = np.append(x, [0.1,0.5,0.9])
k = -5
y = np.append(y, [k, k ,k])
plt.scatter(x,y,marker='o',color='black')

n = len(x)
Xtrain = np.expand_dims(x, axis=1)
ytrain = np.expand_dims(y, axis=1)

# 使用高斯分布作为似然函数
lr = LinearRegression()
lr.fit(Xtrain,ytrain)
Xtest = np.arange(0,1,0.1)
gaussYhat = lr.predict(np.expand_dims(Xtest, axis=1))
plt.plot(Xtest, gaussYhat, color='r',label="least squares")
plt.scatter(Xtest, gaussYhat, marker='o', linewidths=2, color='', edgecolors='r')

#plt.show()

# 使用拉普拉斯分布作为似然函数
f = np.append(np.zeros(2), np.ones(2*n))
Aeq = np.hstack((np.ones((n,1)), Xtrain, np.eye(n), -1*np.eye(n)))
beq = ytrain
l = [None]*(2) + [0]*(2*n)
u = [None]*(2+2*n)
lu = zip(l,u)
res = linprog(f, A_eq=Aeq, b_eq=beq, bounds=tuple(lu))
laplace_w = res.x[0:2]
laplaceYhat = laplace_w[1]*Xtest + laplace_w[0]
plt.plot(Xtest, laplaceYhat, color='b',label="laplace")
plt.scatter(Xtest, laplaceYhat, marker='s', linewidths=2, color='', edgecolors='b')

# 使用huber损失作为目标函数,这一块内容目前不够正确，待定
def HuberLoss(w):
    #r = np.dot(np.hstack((np.ones((n,1)),Xtrain)), w) - ytrain
    r = w[0] + w[1]*Xtrain-ytrain
    temp1 = r[np.abs(r)<= 1]
    if len(temp1) == 0: temp1=0
    temp2 = r[np.abs(r) > 1]
    if len(temp2) == 0: temp2=0
    temp3 = np.abs(r) > 1
    if len(temp3) == 0: temp3=0
    f= (1/2)*np.sum(temp1**2) + 1*np.sum(np.abs(temp2)) - 1/2*np.sum(temp3)*(1**2)
    return f
w_int = (lr.intercept_[0], lr.coef_[0][0])
loss = HuberLoss(w_int)
res = minimize(lambda w:HuberLoss(w), w_int)
huber_w = res.x
huberYhat = huber_w[0] + huber_w[1]*Xtest
plt.plot(Xtest, huberYhat, color='green', label='huberloss')
plt.scatter(Xtest, huberYhat, marker='s', color='', edgecolors='green', linewidths=2)

plt.legend()
plt.xlim([min(Xtest),max(Xtest)])
plt.show()

max_stars_repo_name: xishansnow/MLAPP

max_issues_repo_name: xishansnow/MLAPP

max_forks_repo_name: xishansnow/MLAPP

reason: import numpy,from scipy


text:
'''
Created on 4 Sep 2015

@author: maxz
'''
import unittest
import numpy as np
import GPy

class BGPLVMTest(unittest.TestCase):

    def setUp(self):
        np.random.seed(12345)
        X, W = np.random.normal(0,1,(100,6)), np.random.normal(0,1,(6,13))
        Y = X.dot(W) + np.random.normal(0, .1, (X.shape[0], W.shape[1]))
        self.inan = np.random.binomial(1, .1, Y.shape).astype(bool)
        self.X, self.W, self.Y = X,W,Y
        self.Q = 3
        self.m_full = GPy.models.BayesianGPLVM(Y, self.Q)

    def test_lik_comparisons_m1_s0(self):
        # Test if the different implementations give the exact same likelihood as the full model.
        # All of the following settings should give the same likelihood and gradients as the full model:
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=True, stochastic=False)
        m[:] = self.m_full[:]
        np.testing.assert_almost_equal(m.log_likelihood(), self.m_full.log_likelihood(), 7)
        np.testing.assert_allclose(m.gradient, self.m_full.gradient)
        assert(m.checkgrad())

    def test_predict_missing_data(self):
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=True, stochastic=True, batchsize=self.Y.shape[1])
        m[:] = self.m_full[:]
        np.testing.assert_almost_equal(m.log_likelihood(), self.m_full.log_likelihood(), 7)
        np.testing.assert_allclose(m.gradient, self.m_full.gradient)

        self.assertRaises(NotImplementedError, m.predict, m.X, full_cov=True)

        mu1, var1 = m.predict(m.X, full_cov=False)
        mu2, var2 = self.m_full.predict(self.m_full.X, full_cov=False)
        np.testing.assert_allclose(mu1, mu2)
        np.testing.assert_allclose(var1, var2)

        mu1, var1 = m.predict(m.X.mean, full_cov=True)
        mu2, var2 = self.m_full.predict(self.m_full.X.mean, full_cov=True)
        np.testing.assert_allclose(mu1, mu2)
        np.testing.assert_allclose(var1[:,:,0], var2)

        mu1, var1 = m.predict(m.X.mean, full_cov=False)
        mu2, var2 = self.m_full.predict(self.m_full.X.mean, full_cov=False)
        np.testing.assert_allclose(mu1, mu2)
        np.testing.assert_allclose(var1[:,[0]], var2)

    def test_lik_comparisons_m0_s0(self):
        # Test if the different implementations give the exact same likelihood as the full model.
        # All of the following settings should give the same likelihood and gradients as the full model:
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=False, stochastic=False)
        m[:] = self.m_full[:]
        np.testing.assert_almost_equal(m.log_likelihood(), self.m_full.log_likelihood(), 7)
        np.testing.assert_allclose(m.gradient, self.m_full.gradient)
        assert(m.checkgrad())

    def test_lik_comparisons_m1_s1(self):
        # Test if the different implementations give the exact same likelihood as the full model.
        # All of the following settings should give the same likelihood and gradients as the full model:
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=True, stochastic=True, batchsize=self.Y.shape[1])
        m[:] = self.m_full[:]
        np.testing.assert_almost_equal(m.log_likelihood(), self.m_full.log_likelihood(), 7)
        np.testing.assert_allclose(m.gradient, self.m_full.gradient)
        assert(m.checkgrad())

    def test_lik_comparisons_m0_s1(self):
        # Test if the different implementations give the exact same likelihood as the full model.
        # All of the following settings should give the same likelihood and gradients as the full model:
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=False, stochastic=True, batchsize=self.Y.shape[1])
        m[:] = self.m_full[:]
        np.testing.assert_almost_equal(m.log_likelihood(), self.m_full.log_likelihood(), 7)
        np.testing.assert_allclose(m.gradient, self.m_full.gradient)
        assert(m.checkgrad())

    def test_gradients_missingdata(self):
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=True, stochastic=False, batchsize=self.Y.shape[1])
        assert(m.checkgrad())

    def test_gradients_missingdata_stochastics(self):
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=True, stochastic=True, batchsize=1)
        assert(m.checkgrad())
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=True, stochastic=True, batchsize=4)
        assert(m.checkgrad())

    def test_gradients_stochastics(self):
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=False, stochastic=True, batchsize=1)
        assert(m.checkgrad())
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=False, stochastic=True, batchsize=4)
        assert(m.checkgrad())

    def test_predict(self):
        # Test if the different implementations give the exact same likelihood as the full model.
        # All of the following settings should give the same likelihood and gradients as the full model:
        m = GPy.models.bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch(self.Y, self.Q, missing_data=True, stochastic=True, batchsize=self.Y.shape[1])
        m[:] = self.m_full[:]
        np.testing.assert_almost_equal(m.log_likelihood(), self.m_full.log_likelihood(), 7)
        np.testing.assert_allclose(m.gradient, self.m_full.gradient)
        assert(m.checkgrad())

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

max_stars_repo_name: beckdaniel/GPy

max_issues_repo_name: beckdaniel/GPy

max_forks_repo_name: beckdaniel/GPy

reason: import numpy


text:
"""

File input/output for localization data in SMLM files.

File specifications are provided at https://github.com/imodpasteur/smlm-file-format/blob/master/specification.md.

Code is adapted from https://github.com/imodpasteur/smlm-file-format/blob/master/implementations/Python/smlm_file.py.
(MIT license)
"""

import json
import logging
import time
import zipfile

import numpy as np
import pandas as pd
from google.protobuf import json_format

import locan.constants
from locan.data import metadata_pb2
from locan.data.locdata import LocData
from locan.locan_io.locdata import manifest_pb2
from locan.locan_io.locdata.utilities import (
    convert_property_names,
    convert_property_types,
)
from locan.utils.format import _time_string

__all__ = [
    "manifest_format_from_locdata",
    "manifest_file_info_from_locdata",
    "manifest_from_locdata",
    "save_SMLM",
    "load_SMLM_manifest",
    "load_SMLM_header",
    "load_SMLM_file",
]

logger = logging.getLogger(__name__)

dtype2length = {
    "int8": 1,
    "uint8": 1,
    "int16": 2,
    "uint16": 2,
    "int32": 4,
    "uint32": 4,
    "int64": 8,
    "uint64": 8,
    "float32": 4,
    "float64": 8,
}

def manifest_format_from_locdata(locdata):
    """
    Prepare a manifest["format"] protobuf from locdata.
    The manifest holds metadata for smlm files.

    Parameters
    ----------
    locdata : LocData
        The LocData object.

    Returns
    -------
    locan.locan_io.locdata.manifest_pb2.Format
        The manifest format
    """
    format = manifest_pb2.Format()

    # required by
    format.name = "smlm-table(binary)-0"
    format.type = manifest_pb2.Type.TABLE
    format.mode = manifest_pb2.Mode.BINARY
    format.columns = len(locdata.data.columns)
    format.headers.extend(locdata.data.columns)
    format.dtype.extend(
        [manifest_pb2.Dtype.Value(str(dt).upper()) for dt in locdata.data.dtypes.values]
    )
    format.shape.extend([1] * len(locdata.data.columns))
    format.units.extend([""] * len(locdata.data.columns))  # todo: add correct units

    # optional by specifications
    format.description = "localization table generated from locan.LocData"

    return format

def manifest_file_info_from_locdata(locdata):
    """
    Prepare a manifest["file_info"] protobuf from locdata.
    The manifest holds metadata for smlm files.

    Parameters
    ----------
    locdata : LocData
        The LocData object.

    Returns
    -------
    locan.locan_io.locdata.manifest_pb2.FileInfo
        The manifest file information
    """
    file_info = manifest_pb2.FileInfo()

    # required by specifications
    file_info.name = "table-0.bin"
    file_info.type = manifest_pb2.TABLE
    file_info.format = "smlm-table(binary)-0"
    file_info.channel = "default"
    file_info.rows = len(locdata)
    for column in locdata.data.columns:
        file_info.offset[column] = 0

    # optional by specifications
    for column in locdata.data.columns:
        file_info.min[column] = locdata.data[column].min()
    for column in locdata.data.columns:
        file_info.max[column] = locdata.data[column].max()
    # file_info.exposure = 20

    return file_info

def manifest_from_locdata(locdata, return_json_string=False):
    """
    Prepare a manifest protobuf from locdata.
    The manifest holds metadata for smlm files.

    Parameters
    ----------
    locdata : LocData
        The LocData object.
    return_json_string : bool
        Flag for returning json string

    Returns
    -------
    locan.locan_io.locdata.manifest_pb2.Manifest, str
        The manifest
    """
    manifest = manifest_pb2.Manifest()
    format = manifest_format_from_locdata(locdata)
    file_info = manifest_file_info_from_locdata(locdata)

    # required by specifications
    manifest.format_version = "0.2"
    manifest.formats["smlm-table(binary)-0"].CopyFrom(format)
    manifest.files.append(file_info)

    # recommended by specifications
    manifest.name = locdata.meta.identifier
    manifest.description = locdata.meta.comment
    manifest.tags.append("SMLM")
    manifest.thumbnail = ""
    manifest.sample = ""
    manifest.labeling = ""
    manifest.date = _time_string(time.time())

    # optional by specifications
    manifest.author = ""
    manifest.citation = ""
    manifest.email = ""
    manifest.locdata_meta.CopyFrom(locdata.meta)

    if return_json_string:
        json_string = json_format.MessageToJson(
            manifest,
            preserving_proto_field_name=True,
            including_default_value_fields=False,
        )
        json_string = _change_upper_to_lower_keys(json_string)
        return json_string
    else:
        return manifest

def _change_upper_to_lower_keys(json_string: str) -> str:
    """ Switch selected key words in json_string from upper to lower letters."""
    json_string = json_string.replace("BINARY", "binary")
    json_string = json_string.replace("TEXT", "text")
    json_string = json_string.replace("TABLE", "table")
    json_string = json_string.replace("IMAGE", "image")
    json_string = json_string.replace("INT", "int")
    json_string = json_string.replace("UINT", "uint")
    json_string = json_string.replace("FLOAT", "float")
    return json_string

def save_SMLM(locdata, path, manifest=None):
    """
    Save LocData attributes in a SMLM single-molecule
    localization (zip) file with manifest.json (version 0.2).

    In the smlm file format we store metadata as a human-
    readable manifest. The data is stored as byte string.

    Note
    ----
    Only selected LocData attributes are saved.
    These are: 'data', 'columns', 'meta'.

    Parameters
    ----------
    locdata : LocData
        The LocData object to be saved.
    path : str, bytes, os.PathLike, file-like
        File path including file name to save to.
    manifest : locan.locan_io.locdata.manifest_pb2.Manifest
        Protobuf with manifest to use instead of an autogenerated manifest.
    """
    if manifest is None:
        manifest_json = manifest_from_locdata(locdata, return_json_string=True)
    elif isinstance(manifest, str):
        manifest_json = manifest
    elif isinstance(manifest, manifest_pb2.Manifest):
        manifest_json = json_format.MessageToJson(
            manifest,
            preserving_proto_field_name=True,
            including_default_value_fields=False,
        )
        manifest_json = _change_upper_to_lower_keys(manifest_json)
    else:
        raise TypeError("Type of manifest is not correct.")

    byte_string = b"".join(
        [
            locdata.data[column][i].tobytes()
            for i in range(len(locdata))
            for column in locdata.data.columns
        ]
    )

    # save zip file
    with zipfile.ZipFile(
        path, "w", compression=zipfile.ZIP_DEFLATED, allowZip64=True
    ) as zf:
        # write manifest
        zf.writestr("manifest.json", manifest_json)
        # write files
        zf.writestr("table-0.bin", byte_string)

def load_SMLM_manifest(path):
    """
    Read manifest.json (version 0.2) from a SMLM single-molecule localization (zip) file.

    Parameters
    ----------
    path : str, bytes, os.PathLike, file-like
        File path for a SMLM file to load.

    Returns
    -------
    dictionary
        manifest in json format
    """
    zf = zipfile.ZipFile(path, "r")
    file_names = zf.namelist()
    if "manifest.json" not in file_names:
        raise Exception("invalid file: no manifest.json found in the smlm file.")
    manifest = json.loads(zf.read("manifest.json"))
    assert manifest["format_version"] == "0.2"
    return manifest

def load_SMLM_header(path):
    """
    Read header (manifest) from a SMLM single-molecule localization file and identify column names.

    Parameters
    ----------
    path : str, bytes, os.PathLike, file-like
        File path for a SMLM file to load.

    Returns
    -------
    list(str)
        A list of dataset property keys as derived from the SMLM identifiers.
    """
    zf = zipfile.ZipFile(path, "r")
    file_names = zf.namelist()
    if "manifest.json" not in file_names:
        raise Exception("invalid file: no manifest.json found in the smlm file.")
    manifest = json.loads(zf.read("manifest.json"))
    assert manifest["format_version"] == "0.2"

    locdata_columns_list = []
    for file_info in manifest["files"]:
        if file_info["type"] != "table":
            logger.info("ignore file with type: %s", file_info["type"])
        else:
            name = file_info["name"]
            logger.info(f"loading table {name}....")
            format_key = file_info["format"]
            file_format = manifest["formats"][format_key]
            if name not in file_names:
                logger.error("ERROR: Did not find %s in zip file", file_info["name"])
            else:
                headers = file_format["headers"]

                column_keys = convert_property_names(
                    properties=headers, property_mapping=locan.constants.SMLM_KEYS
                )
                locdata_columns_list.append(column_keys)

    if len(locdata_columns_list) == 1:
        return locdata_columns_list[0]
    else:
        return locdata_columns_list

def load_SMLM_file(path, nrows=None, convert=True):
    """
    Load data from a SMLM single-molecule localization file.

    Parameters
    ----------
    path : str, bytes, os.PathLike, file-like
        File path for a SMLM file to load.
    nrows : int, None
        The number of localizations to load from file. None means that all available rows are loaded.
    convert : bool
        If True convert types by applying type specifications in locan.constants.PROPERTY_KEYS.

    Returns
    -------
    LocData, list(LocData)
        A new instance of LocData with all localizations. Returns a list of LocData if multiple tables are found.
    """
    zf = zipfile.ZipFile(path, "r")
    file_names = zf.namelist()
    if "manifest.json" not in file_names:
        raise Exception("invalid file: no manifest.json found in the smlm file.")
    manifest = json.loads(zf.read("manifest.json"))
    assert manifest["format_version"] == "0.2"

    locdatas = []
    for file_info in manifest["files"]:
        if file_info["type"] != "table":
            logger.info("ignore file with type: %s", file_info["type"])
        else:
            name = file_info["name"]
            logger.debug(f"start loading {name} ...")
            format_key = file_info["format"]
            file_format = manifest["formats"][format_key]
            if file_format["mode"] != "binary":
                raise Exception(f"format mode {file_format['mode']} not supported.")
            else:
                try:
                    table_file = zf.read(file_info["name"])
                except KeyError:
                    logger.error(
                        "ERROR: Did not find %s in zip file", file_info["name"]
                    )
                    continue
                else:
                    logger.debug(f"loading {len(table_file)} bytes")
                    headers = file_format["headers"]
                    dtype = file_format["dtype"]
                    shape = file_format["shape"]
                    cols = len(headers)
                    rows = nrows if nrows is not None else int(file_info["rows"])

                    logger.debug("columns: %s", headers)
                    logger.debug("rows: %s, columns: %s", rows, cols)
                    assert len(headers) == len(dtype) == len(shape)

                    rowLen = sum(
                        dtype2length[dtype[i]] for i, header in enumerate(headers)
                    )
                    tableDict = {}
                    byteOffset = 0
                    for i, header in enumerate(headers):
                        tableDict[header] = np.ndarray(
                            (rows,),
                            buffer=table_file,
                            dtype=dtype[i],
                            offset=byteOffset,
                            order="C",
                            strides=(rowLen,),
                        )
                        byteOffset += dtype2length[dtype[i]]
                    logger.debug(f"finished loading {name}")

                    dataframe = pd.DataFrame.from_dict(tableDict)

                    column_keys = {}
                    for header in headers:
                        if header in locan.constants.SMLM_KEYS:
                            column_keys[header] = locan.constants.SMLM_KEYS[header]
                        elif header in locan.constants.RAPIDSTORM_KEYS:
                            column_keys[header] = locan.constants.RAPIDSTORM_KEYS[
                                header
                            ]
                        else:
                            logger.warning(
                                f"Column {header} is not a Locan property standard."
                            )
                            column_keys[header] = header
                    dataframe.rename(columns=column_keys, inplace=True)

                    if convert:
                        dataframe = convert_property_types(
                            dataframe, types=locan.constants.PROPERTY_KEYS
                        )

                    locdata = LocData.from_dataframe(dataframe=dataframe)

                    # todo: check on taking correct metadata
                    locdata.meta.source = metadata_pb2.EXPERIMENT
                    locdata.meta.state = metadata_pb2.RAW
                    locdata.meta.file.type = metadata_pb2.SMLM
                    locdata.meta.file.path = str(path)

                    for property_ in sorted(
                        list(
                            set(dataframe.columns).intersection(
                                {"position_x", "position_y", "position_z"}
                            )
                        )
                    ):
                        locdata.meta.localization_properties.add(
                            name=property_, unit="nm", type="float"
                        )

                    del locdata.meta.history[:]
                    locdata.meta.history.add(
                        name="load_SMLM_file", parameter=f"path={path}, nrows={nrows}"
                    )

                    locdatas.append(locdata)

        if len(locdatas) == 1:
            return locdatas[0]
        else:
            return locdatas

max_stars_repo_name: super-resolution/Locan

max_issues_repo_name: super-resolution/Locan

max_forks_repo_name: super-resolution/Locan


text:
import pygame
import numpy
from os import system
from world import *
from robot import*
import time

clock = pygame.time.Clock()
np.random.seed(100)
h=600
w=600
gameDisplay=pygame.display.set_mode((h,w))
pygame.display.set_caption("Title")
gameDisplay.blit(pygame.transform.flip(gameDisplay, False, True), dest=(0, 0))

crashed = False; done = False

sim = World(h,w)

robo = RobotNonHolonomic_differential_drive(sim.spawn, sim.goal)
# robo = RobotNonHolonomic_tricycle_drive(sim.spawn, sim.goal)
# robo = RobotHolonomic(sim.spawn, sim.goal)
robot_size = robo.robot_size + 3

obstacles = [[0,0,100,100],
             [100,150,10,100],
             [250,250,10,300],
             [100, 245, 150,10],
             [350, 0, 10, 300],
             [350, 295, 170, 10],
             [250, 515, 300, 10],
             [450, 325, 10, 140],
             [350, 395, 200, 10],
             [220, 50, 10, 150],
             [145, 125, 150, 10]]

sim.add_obstacles(obstacles)
sim.grow_obstacles(robot_size)

num_steps = 100000
step = 0
while not (crashed):
    print(step)
    step = step+1
    # time.sleep(1)

    # crashing the simulation manually
    for i in pygame.event.get():
      if i.type == pygame.KEYDOWN:
        if i.unicode == "q":
          crashed = True

    if step%5 == 0:
      goal = sim.goal
    else:
      goalx = np.random.uniform(0,h)
      goaly = np.random.uniform(0,w)
      goal = (goalx, goaly)

    done = robo.RRT_step(goal, sim)

    gameDisplay.fill((255,255,255))
    sim.print_obstacles(gameDisplay, robot_size)
    robo.print_paths(gameDisplay, type='center', done=done)
    # robo.print_paths(gameDisplay, type='wheels', done=done)
    pygame.display.update()
    clock.tick(60)

    # if step>num_steps:
    #   done = True


text:
import IMLearn.learners.regressors.linear_regression
from IMLearn.learners.regressors import PolynomialFitting
from IMLearn.utils import split_train_test

import numpy as np
import pandas as pd
import plotly.express as px
import plotly.io as pio
import plotly.graph_objects as go
pio.templates.default = "simple_white"

def load_data(filename: str) -> pd.DataFrame:
    """
    Load city daily temperature dataset and preprocess data.
    Parameters
    ----------
    filename: str
        Path to house prices dataset

    Returns
    -------
    Design matrix and response vector (Temp)
    """
    df = pd.read_csv(filename, parse_dates=["Date"])
    df["DayOfYear"] = df["Date"].dt.dayofyear
    df.dropna(inplace=True)
    return df

if __name__ == '__main__':
    np.random.seed(0)
    # Question 1 - Load and preprocessing of city temperature dataset
    df = load_data("C:\\Users\\DvirAssulin\\Desktop\\university\\IML.HUJI\\datasets\\City_Temperature.csv")

    # Question 2 - Exploring data for specific country
    df_Israel = df[df["Country"] == "Israel"]
    df_Israel = df_Israel[df_Israel["Temp"] > -5]
    fig1 = go.Figure(go.Scatter(x=df_Israel["DayOfYear"], y=df_Israel["Temp"], mode='markers', marker=dict(color=df_Israel["Year"]), showlegend=False))
    fig1.update_layout(title_text="$\\text{Temperature as Function of Day of Year}$", title_x = 0.5)
    fig1.update_xaxes(title_text="$\\text{Day of Year}$")
    fig1.update_yaxes(title_text="$Temp$")
    fig1.show()
    monthly_temp_Israel = df_Israel.groupby("Month").agg("std")
    fig2 = go.Figure(go.Bar(x=monthly_temp_Israel.index, y=monthly_temp_Israel["Temp"], marker=dict(color="blue"), showlegend=False))
    fig2.update_layout(title_text="$\\text{Standard Deviation of Temperature by Month}$", title_x = 0.5)
    fig2.update_xaxes(title_text="$\\text{Month}$")
    fig2.update_yaxes(title_text="$Stdev$")
    fig2.show()

    # Question 3 - Exploring differences between countries
    monthly_temp_std = df.groupby("Month").agg("std")
    monthly_temp_avg = df.groupby("Month").agg("mean")
    month_country_temp_std = df.groupby(["Country", "Month"]).agg("std")
    month_country_temp_avg = df.groupby(["Country", "Month"]).agg("mean")
    month_country_temp_avg.reset_index(inplace=True)
    fig3 = px.line(month_country_temp_avg,"Month", "Temp", color="Country", error_y=month_country_temp_std["Temp"])
    fig3.add_trace(go.Scatter(x=monthly_temp_avg.index,y=monthly_temp_avg["Temp"],name="All Countries"))
    fig3.update_layout(title_text="$\\text{Average Temperature of All Countries by Month}$", title_x = 0.5)
    fig3.update_xaxes(title_text="$\\text{Month}$")
    fig3.update_yaxes(title_text="$Mean$")
    fig3.show()

    # Question 4 - Fitting model for different values of `k`
    y = df_Israel["Temp"]
    X = df_Israel["DayOfYear"]
    train_x = X.sample(frac=0.75)
    train_y = y.loc[train_x.index]
    X.drop(train_x.index, inplace=True)
    y.drop(train_y.index, inplace=True)
    test_x = X
    test_y = y
    loss_pred = []
    for k in range(1, 11):
        poly_model = PolynomialFitting(k)
        poly_model.fit(train_x.to_numpy(), train_y.to_numpy())
        loss = poly_model.loss(test_x.to_numpy(), test_y.to_numpy())
        loss = loss.round(2)
        loss_pred.append(loss)
        print(f"The loss for degree of {k} is {loss}\n")
    fig4 = go.Figure(go.Bar(x=[k for k in range(1, 11)], y=loss_pred, marker=dict(color="green"), showlegend=False))
    fig4.update_layout(title_text="$\\text{Loss Per Degree of K}$", title_x = 0.5)
    fig4.update_xaxes(title_text="$\\text{Degree}$")
    fig4.update_yaxes(title_text="$Loss$")
    fig4.show()

    # Question 5 - Evaluating fitted model on different countries
    poly_model = PolynomialFitting(5)
    poly_model.fit(df_Israel["DayOfYear"], df_Israel["Temp"])
    countries = ["The Netherlands", "South Africa", "Jordan"]
    country_loss = []
    for country in countries:
        sub_df = df[df["Country"] == country]
        X = sub_df["DayOfYear"]
        y = sub_df["Temp"]
        loss = poly_model.loss(X, y)
        loss.round(2)
        country_loss.append(loss)

    fig5 = go.Figure(go.Bar(x=countries, y=country_loss, marker=dict(color="red"), showlegend=False))
    fig5.update_layout(title_text="$\\text{Loss Per Country}$", title_x=0.5)
    fig5.update_xaxes(title_text="$\\text{Country}$")
    fig5.update_yaxes(title_text="$Loss$")
    fig5.show()


text:
import imageio
import numpy as np
from glob import glob
import os
from datetime import datetime
import tqdm
import random
import pickle
import modules
from skimage.util.shape import view_as_blocks
training_folder="./flower/"
now = datetime.now()

def read_images():
    x_truth=[]
    y_truth=[]
    file_order=[]
    i=0
    print("Reading Images\n")
    for image in tqdm.tqdm(glob(training_folder+"*.jpg")):
        basename=os.path.basename(image)
        filename=os.path.splitext(basename)[0]
        x_truth.append(imageio.imread(training_folder+filename+"_bw.jpeg")/256)
        y_truth.append(imageio.imread(image)/256)

        file_order.append(filename)
        i=i+1
        if i>3001:
            break
    return x_truth, y_truth, file_order

def generate_validate_set(size, prob=0.1):
    validation_list=random.sample(range(size),int(size*prob*2))
    training_list = [x for x in list(range(len(validation_list))) if x not in validation_list]
    test_list=validation_list[:len(validation_list)//2]
    validation_list = validation_list[len(validation_list) // 2:]

    return training_list,validation_list,test_list

def write_current_configuration(validation_list,test_list,file_order):
    pickle.dump(validation_list,open("validate_list_"+str(now)+".p", "wb"))
    pickle.dump(test_list, open("test_list_" + str(now) + ".p", "wb"))
    pickle.dump(file_order, open("file_order" + str(now) + ".p", "wb"))

def normalize_image(x_image,y_image):
    x_image=x_image/256
    y_image=y_image/256
    return x_image,y_image

def prepare_data(x_truth,y_truth,training_list,validation_list,test_list):
    x_data=[]
    y_data=[]
    x_validation_data=[]
    y_validation_data = []
    x_test_data=[]
    y_test_data = []
    for i in range(len(x_truth)):
        if i in training_list:
            x_data.append(np.array(x_truth[i]).flatten())
            y_data.append(np.array(y_truth[i]))
        else:
            if i in validation_list:
                x_validation_data.append(np.array(x_truth[i]).flatten())
                y_validation_data.append(np.array(y_truth[i]))
            else:
                x_test_data.append(np.array(x_truth[i]).flatten())
                y_test_data.append(np.array(y_truth[i]))
    return x_data,y_data,x_validation_data,y_validation_data,x_test_data,y_test_data

def forward_pro_2(x_data, networks):
    memory = {}
    i = 0
    for layer in networks:
        if i == 0:
            lo = layer.forward(x_data)
            memory[str(i) + "_prev"] = x_data
            # print("forward_x_data_shape")
            # print(x_data.shape)
        else:
            lo = layer.forward(prev)
            memory[str(i) + "_prev"] = prev
            ##print("forward_prev_shape")
            # print(prev.shape)
        # print(lo.shape)
        memory[str(i) + "_middle"] = lo
        prev = modules.relu(lo)
        memory[str(i) + "_after"] = prev
        # print(prev.shape)
        # print(lo.shape)

        i = i + 1
    return prev, memory, i

def back_pro_2(dLoss, memory, networks, i):
    i = i - 1
    for layers in reversed(networks):
        after = memory[str(i) + "_middle"]
        # print(str(i)+"_middle")
        # print(str(i)+"_prev")
        # print("after")
        # print(after)
        d_prev = modules.relu_backward(dLoss, after)
        # print(d_prev)
        before = memory[str(i) + "_prev"]
        # print(before)
        # print(len(before))
        # print(len(d_prev))
        dLoss = layers.backward(d_prev)
        # print(d_prev)
        # print(layers.weights)
        i = i - 1
    return networks

def forward_pro(x_data,networks):
    memory={}
    i=0
    for layer in networks:
        if i==0:
            lo=layer.forward(x_data)
            memory[str(i)+"_prev"]=x_data
            #print("forward_x_data_shape")
            #print(x_data.shape)
        else:
            lo=layer.forward(prev)
            memory[str(i) + "_prev"] = prev
            ##print("forward_prev_shape")
            #print(prev.shape)
        #print("forward")
        #print(lo.shape)
        memory[str(i) + "_middle"]=lo
        prev=modules.sigmoid(lo)
        memory[str(i) + "_after"] = prev
        #print(prev.shape)
        #print(lo.shape)

        i=i+1
    return prev, memory,i

def back_pro(dLoss,memory,networks,i):
    i=i-1
    for layers in reversed(networks):
        after=memory[str(i)+"_middle"]
        #print(str(i)+"_middle")
        #print(str(i)+"_prev")
        # print("after")
        # print(after)
        d_prev=modules.sigmoid_backward(dLoss,after)
        #print(d_prev)
        before=memory[str(i)+"_prev"]
        #print(before)
        #print(len(before))
        #print(len(d_prev))
        dLoss=layers.backward(d_prev)
        #print(d_prev)
        #print(layers.weights)
        i=i-1
    return networks

# def begin_training_patch(x_data, y_data, x_validation_data, y_validation_data, x_test_data, y_test_data):
#     networks_r = []
#     networks_g = []
#     networks_b = []
#     lr = 0.1
#     basesize = 64 #8*8 patch
#     # print(basesize)
#     networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     # print(1)
#     networks_g.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     networks_g.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     networks_g.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     # print(2)
#     networks_b.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     networks_b.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     networks_b.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
#     # print(3)
#     #minibatch_size = 32
#     # Forward Propogation
#     for x in tqdm.tqdm(range(100)):
#         sum_r = 0
#         sum_g = 0
#         sum_b = 0
#         for p in range(0, 30):
#             patch = view_as_blocks(x_data[p], block_shape=(8, 8))
#             p
#             #x_to_train = x_data[p]  # copy 3 times for R,G,B channel
#             # print(x_to_train.shape)
#             # print(x_data[p].shape)
#             #x_to_train = x_to_train.reshape((basesize, 1))
#
#         batch=np.concatenate(batch,axis=-1)
#         batch_y = np.stack(batch_y,axis=-1)
#         # print(batch.shape)
#         # print(batch_y.shape)
#         # print(batch)
#         # print("data_prepared")
#         y_pred_r, memory_r, i_r = forward_pro(batch, networks_r)
#         y_pred_g, memory_g, i_g = forward_pro(batch, networks_g)
#         y_pred_b, memory_b, i_b = forward_pro(batch, networks_b)
#         y_r_pred=np.zeros((basesize,basesize))
#         y_g_pred=np.zeros((basesize,basesize))
#         y_b_pred=np.zeros((basesize,basesize))
#         dMSE_r_s=np.zeros((basesize,minibatch_size))
#         dMSE_g_s = np.zeros((basesize, minibatch_size))
#         dMSE_b_s = np.zeros((basesize, minibatch_size))
#         for p in range(0,minibatch_size):
#             y_r_pred += y_pred_r[:,p]
#             y_g_pred += y_pred_g[:, p]
#             y_b_pred += y_pred_b[:, p]
#             #a picture
#             MSE_r = np.square(
#                 np.subtract(batch_y[:, :, 0,p].reshape((basesize, 1)), y_pred_r)).mean()  # accuracy
#             MSE_g = np.square(np.subtract(batch_y[:, :, 1,p].reshape((basesize, 1)), y_pred_g)).mean()
#             MSE_b = np.square(np.subtract(batch_y[:, :, 2,p].reshape((basesize, 1)), y_pred_b)).mean()
#             dMSE_r = np.subtract(y_pred_r[p], batch_y[:, :, 0,p].reshape((basesize, 1)))
#             dMSE_g = np.subtract(y_pred_g[p], batch_y[:, :, 1,p].reshape((basesize, 1)))
#             dMSE_b = np.subtract(y_pred_b[p], batch_y[:, :, 2,p].reshape((basesize, 1)))
#             dMSE_r_s += dMSE_r
#             dMSE_g_s += dMSE_g
#             dMSE_b_s += dMSE_b
#             # print(dMSE_r.shape)
#             # print(dMSE_b.shape)
#             # print(dMSE_g.shape)
#             sum_r+=MSE_r
#             sum_g+=MSE_g
#             sum_b+=MSE_b
#
#
#         dMSE_r_s=dMSE_r_s/minibatch_size
#         dMSE_g_s = dMSE_g_s / minibatch_size
#         dMSE_b_s = dMSE_b_s / minibatch_size
#
#             # y_pred_image=y_pred.reshape(64,64,3)
#             # print(y_pred)
#         y_data_to_compare = y_data[p]
#
#             # print("dMSE_r")
#             # print(dMSE_r)
#             # print("dMSE_g")
#             # print(dMSE_g)
#             # print("dMSE_b")
#             # print(dMSE_b)
#             # y_back=dMSE.reshape(basesize*3) #linearize the stuff
#
#             # Start backward propogation
#         networks_r = back_pro(dMSE_r_s, memory_r, networks_r, i_r)
#         networks_g = back_pro(dMSE_g_s, memory_g, networks_g, i_g)
#         networks_b = back_pro(dMSE_b_s, memory_b, networks_b, i_b)
#             # print(MSE_r)
#             # print(MSE_g)
#             # print(MSE_b)
#         # sum_r += MSE_r
#         #         # sum_g += MSE_g
#         #         # sum_b += MSE_b
#             # print(y_pred)
#             # print(y_data_to_compare)
#         pickle.dump(networks_r, open("networks_r_dump_" + str(x) + ".p", "wb"))
#         pickle.dump(networks_g, open("networks_g_dump_" + str(x) + ".p", "wb"))
#         pickle.dump(networks_b, open("networks_b_dump_" + str(x) + ".p", "wb"))
#         print(sum_r / 30)
#         print(sum_g / 30)
#         print(sum_b / 30)
#

def begin_training_minibatch(x_data, y_data, x_validation_data, y_validation_data, x_test_data, y_test_data):
    networks_r = []
    networks_g = []
    networks_b = []
    lr = 0.1
    basesize = x_data[0].size
    # print(basesize)
    networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    # print(1)
    networks_g.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_g.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_g.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_g.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    # print(2)
    networks_b.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_b.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_b.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_b.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    # print(3)
    minibatch_size = 32
    # Forward Propogation
    for x in tqdm.tqdm(range(100)):
        sum_r = 0
        sum_g = 0
        sum_b = 0
        batch=[]
        batch_y=[]
        for p in range(0, minibatch_size):
            batch.append(x_data[p].reshape((basesize,1)))
            batch_y.append(y_data[p])
            #x_to_train = x_data[p]  # copy 3 times for R,G,B channel
            # print(x_to_train.shape)
            # print(x_data[p].shape)
            #x_to_train = x_to_train.reshape((basesize, 1))

        batch=np.concatenate(batch,axis=-1)
        batch_y = np.stack(batch_y,axis=-1)
        # print(batch.shape)
        # print(batch_y.shape)
        # print(batch)
        # print("data_prepared")
        y_pred_r, memory_r, i_r = forward_pro(batch, networks_r)
        y_pred_g, memory_g, i_g = forward_pro(batch, networks_g)
        y_pred_b, memory_b, i_b = forward_pro(batch, networks_b)
        y_r_pred=np.zeros((basesize,basesize))
        y_g_pred=np.zeros((basesize,basesize))
        y_b_pred=np.zeros((basesize,basesize))
        dMSE_r_s=np.zeros((basesize,minibatch_size))
        dMSE_g_s = np.zeros((basesize, minibatch_size))
        dMSE_b_s = np.zeros((basesize, minibatch_size))
        for p in range(0,minibatch_size):
            y_r_pred += y_pred_r[:,p]
            y_g_pred += y_pred_g[:, p]
            y_b_pred += y_pred_b[:, p]
            #a picture
            MSE_r = np.square(
                np.subtract(batch_y[:, :, 0,p].reshape((basesize, 1)), y_pred_r)).mean()  # accuracy
            MSE_g = np.square(np.subtract(batch_y[:, :, 1,p].reshape((basesize, 1)), y_pred_g)).mean()
            MSE_b = np.square(np.subtract(batch_y[:, :, 2,p].reshape((basesize, 1)), y_pred_b)).mean()
            dMSE_r = np.subtract(y_pred_r[p], batch_y[:, :, 0,p].reshape((basesize, 1)))
            dMSE_g = np.subtract(y_pred_g[p], batch_y[:, :, 1,p].reshape((basesize, 1)))
            dMSE_b = np.subtract(y_pred_b[p], batch_y[:, :, 2,p].reshape((basesize, 1)))
            dMSE_r_s += dMSE_r
            dMSE_g_s += dMSE_g
            dMSE_b_s += dMSE_b
            # print(dMSE_r.shape)
            # print(dMSE_b.shape)
            # print(dMSE_g.shape)
            sum_r+=MSE_r
            sum_g+=MSE_g
            sum_b+=MSE_b

        dMSE_r_s=dMSE_r_s/minibatch_size
        dMSE_g_s = dMSE_g_s / minibatch_size
        dMSE_b_s = dMSE_b_s / minibatch_size

            # y_pred_image=y_pred.reshape(64,64,3)
            # print(y_pred)
        y_data_to_compare = y_data[p]

            # print("dMSE_r")
            # print(dMSE_r)
            # print("dMSE_g")
            # print(dMSE_g)
            # print("dMSE_b")
            # print(dMSE_b)
            # y_back=dMSE.reshape(basesize*3) #linearize the stuff

            # Start backward propogation
        networks_r = back_pro(dMSE_r_s, memory_r, networks_r, i_r)
        networks_g = back_pro(dMSE_g_s, memory_g, networks_g, i_g)
        networks_b = back_pro(dMSE_b_s, memory_b, networks_b, i_b)
            # print(MSE_r)
            # print(MSE_g)
            # print(MSE_b)
        # sum_r += MSE_r
        #         # sum_g += MSE_g
        #         # sum_b += MSE_b
            # print(y_pred)
            # print(y_data_to_compare)
        pickle.dump(networks_r, open("networks_r_dump_" + str(x) + ".p", "wb"))
        pickle.dump(networks_g, open("networks_g_dump_" + str(x) + ".p", "wb"))
        pickle.dump(networks_b, open("networks_b_dump_" + str(x) + ".p", "wb"))
        print(sum_r / 30)
        print(sum_g / 30)
        print(sum_b / 30)

def begin_training(x_data,y_data,x_validation_data,y_validation_data,x_test_data,y_test_data):
    networks_r=[]
    networks_g = []
    networks_b = []
    lr=0.05
    basesize = x_data[0].size
    #print(basesize)
    networks_r.append(modules.Dense_layer(basesize,basesize,learning_rate=lr))
    networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    networks_r.append(modules.Dense_layer(basesize, basesize, learning_rate=lr))
    #print(1)
    networks_g.append(modules.Dense_layer(basesize,basesize,learning_rate=lr))
    networks_g.append(modules.Dense_layer(basesize,basesize,learning_rate=lr))
    networks_g.append(modules.Dense_layer(basesize,basesize,learning_rate=lr))
    #print(2)
    networks_b.append(modules.Dense_layer(basesize,basesize,learning_rate=lr))
    networks_b.append(modules.Dense_layer(basesize,basesize,learning_rate=lr))
    networks_b.append(modules.Dense_layer(basesize,basesize,learning_rate=lr))
    #print(3)
    minibatch_size=32
    #Forward Propogation
    for x in tqdm.tqdm(range(100)):
        sum_r=0
        sum_g = 0
        sum_b = 0
        for p in range(0, 30):

            x_to_train=x_data[p] #copy 3 times for R,G,B channel
            #print(x_to_train.shape)
            #print(x_data[p].shape)
            x_to_train=x_to_train.reshape((basesize,1))
            y_pred_r,memory_r,i_r =forward_pro(x_to_train,networks_r)
            y_pred_g, memory_g, i_g = forward_pro(x_to_train, networks_g)
            y_pred_b, memory_b, i_b = forward_pro(x_to_train, networks_b)
            #y_pred_image=y_pred.reshape(64,64,3)
            #print(y_pred)
            y_data_to_compare=y_data[p]
            MSE_r = np.square(np.subtract(y_data_to_compare[:,:,0].reshape((basesize,1)), y_pred_r)).mean() #accuracy
            MSE_g = np.square(np.subtract(y_data_to_compare[:, :, 1].reshape((basesize,1)), y_pred_g)).mean()
            MSE_b = np.square(np.subtract(y_data_to_compare[:, :, 2].reshape((basesize,1)), y_pred_b)).mean()
            dMSE_r=np.subtract(y_pred_r,y_data_to_compare[:,:,0].reshape((basesize,1)))
            dMSE_g =np.subtract(y_pred_g, y_data_to_compare[:, :, 1].reshape((basesize,1)))
            dMSE_b =np.subtract(y_pred_b, y_data_to_compare[:, :, 2].reshape((basesize,1)))
            # print("dMSE_r")
            # print(dMSE_r)
            # print("dMSE_g")
            # print(dMSE_g)
            # print("dMSE_b")
            # print(dMSE_b)
            #y_back=dMSE.reshape(basesize*3) #linearize the stuff

            #Start backward propogation
            networks_r=back_pro(dMSE_r,memory_r,networks_r,i_r)
            networks_g = back_pro(dMSE_g, memory_g, networks_g, i_g)
            networks_b = back_pro(dMSE_b, memory_b, networks_b, i_b)
            # print(MSE_r)
            # print(MSE_g)
            # print(MSE_b)
            sum_r+=MSE_r
            sum_g += MSE_g
            sum_b += MSE_b
            #print(y_pred)
            #print(y_data_to_compare)
        pickle.dump(networks_r, open("networks_r_dump_" + str(x) + ".p", "wb"))
        pickle.dump(networks_g, open("networks_g_dump_" + str(x) + ".p", "wb"))
        pickle.dump(networks_b, open("networks_b_dump_" + str(x) + ".p", "wb"))
        output_loss_r = forward_pro()

    #try 1-> Linearize and use 4 dense layer each with 3 densely connected layers

if(os.path.exists("x_truth.p")):
    x_truth=pickle.load(open("x_truth.p","rb"))
    y_truth=pickle.load(open("y_truth.p","rb"))
else:
    x_truth,y_truth,file_order=read_images()
    pickle.dump(x_truth,open("x_truth.p","wb"))
    pickle.dump(y_truth,open("y_truth.p","wb"))
    pickle.dump(file_order,open("file_order.p","wb"))
##print(len(x_truth))
#print(len(y_truth))
training_list,validation_list,test_list=generate_validate_set(len(x_truth))
x_data,y_data,x_validation_data,y_validation_data,x_test_data,y_test_data=prepare_data(x_truth,y_truth,training_list,validation_list,test_list)
#write_current_configuration(training_list,validation_list,test_list)
#begin_training(x_data,y_data,x_validation_data,y_validation_data,x_test_data,y_test_data)
begin_training_minibatch(x_data,y_data,x_validation_data,y_validation_data,x_test_data,y_test_data)


text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst

# TEST_UNICODE_LITERALS

import sys

import pytest
import numpy as np

from ... import table
from ...table import Row
from ...extern.six.moves import zip
from .conftest import MaskedTable

def test_masked_row_with_object_col():
    """
    Numpy < 1.8 has a bug in masked array that prevents access a row if there is
    a column with object type.
    """
    t = table.Table([[1]], dtype=['O'], masked=True)
    t['col0'].mask = False
    assert t[0]['col0'] == 1
    t['col0'].mask = True
    assert t[0]['col0'] is np.ma.masked

@pytest.mark.usefixtures('table_types')
class TestRow():
    def _setup(self, table_types):
        self._table_type = table_types.Table
        self._column_type = table_types.Column

    @property
    def t(self):
        # py.test wants to run this method once before table_types is run
        # to set Table and Column.  In this case just return None, which would
        # cause any downstream test to fail if this happened in any other context.
        if self._column_type is None:
            return None
        if not hasattr(self, '_t'):
            a = self._column_type(name='a', data=[1, 2, 3], dtype='i8')
            b = self._column_type(name='b', data=[4, 5, 6], dtype='i8')
            self._t = self._table_type([a, b])
        return self._t

    def test_subclass(self, table_types):
        """Row is subclass of ndarray and Row"""
        self._setup(table_types)
        c = Row(self.t, 2)
        assert isinstance(c, Row)

    def test_values(self, table_types):
        """Row accurately reflects table values and attributes"""
        self._setup(table_types)
        table = self.t
        row = table[1]
        assert row['a'] == 2
        assert row['b'] == 5
        assert row[0] == 2
        assert row[1] == 5
        assert row.meta is table.meta
        assert row.colnames == table.colnames
        assert row.columns is table.columns
        with pytest.raises(IndexError):
            row[2]
        if sys.byteorder == 'little':
            assert str(row.dtype) == "[('a', '<i8'), ('b', '<i8')]"
        else:
            assert str(row.dtype) == "[('a', '>i8'), ('b', '>i8')]"

    def test_ref(self, table_types):
        """Row is a reference into original table data"""
        self._setup(table_types)
        table = self.t
        row = table[1]
        row['a'] = 10
        if table_types.Table is not MaskedTable:
            assert table['a'][1] == 10

    def test_left_equal(self, table_types):
        """Compare a table row to the corresponding structured array row"""
        self._setup(table_types)
        np_t = self.t.as_array()
        if table_types.Table is MaskedTable:
            with pytest.raises(ValueError):
                self.t[0] == np_t[0]
        else:
            for row, np_row in zip(self.t, np_t):
                assert np.all(row == np_row)

    def test_left_not_equal(self, table_types):
        """Compare a table row to the corresponding structured array row"""
        self._setup(table_types)
        np_t = self.t.as_array()
        np_t['a'] = [0, 0, 0]
        if table_types.Table is MaskedTable:
            with pytest.raises(ValueError):
                self.t[0] == np_t[0]
        else:
            for row, np_row in zip(self.t, np_t):
                assert np.all(row != np_row)

    def test_right_equal(self, table_types):
        """Test right equal"""
        self._setup(table_types)
        np_t = self.t.as_array()
        if table_types.Table is MaskedTable:
            with pytest.raises(ValueError):
                self.t[0] == np_t[0]
        else:
            for row, np_row in zip(self.t, np_t):
                assert np.all(np_row == row)

    def test_convert_numpy_array(self, table_types):
        self._setup(table_types)
        d = self.t[1]

        np_data = np.array(d)
        if table_types.Table is not MaskedTable:
            assert np.all(np_data == d.as_void())
        assert np_data is not d.as_void()
        assert d.colnames == list(np_data.dtype.names)

        np_data = np.array(d, copy=False)
        if table_types.Table is not MaskedTable:
            assert np.all(np_data == d.as_void())
        assert np_data is not d.as_void()
        assert d.colnames == list(np_data.dtype.names)

        with pytest.raises(ValueError):
            np_data = np.array(d, dtype=[(str('c'), 'i8'), (str('d'), 'i8')])

    def test_format_row(self, table_types):
        """Test formatting row"""
        self._setup(table_types)
        table = self.t
        row = table[0]
        assert repr(row).splitlines() == ['<{0} {1}{2}>'
                                          .format(row.__class__.__name__,
                                                  'index=0',
                                                  ' masked=True' if table.masked else ''),
                                          '  a     b  ',
                                          'int64 int64',
                                          '----- -----',
                                          '    1     4']
        assert str(row).splitlines() == [' a   b ',
                                         '--- ---',
                                         '  1   4']

        assert row._repr_html_().splitlines() == ['<{0} {1}{2}>'
                                                  .format(row.__class__.__name__,
                                                          'index=0',
                                                          ' masked=True' if table.masked else ''),
                                                  '
'.format(id(table)),
                                                  '

a	b

',
                                                  '

int64	int64

',
                                                  '
1	4
',
                                                  '
']

    def test_as_void(self, table_types):
        """Test the as_void() method"""
        self._setup(table_types)
        table = self.t
        row = table[0]

        # If masked then with no masks, issue numpy/numpy#483 should come
        # into play.  Make sure as_void() code is working.
        row_void = row.as_void()
        if table.masked:
            assert isinstance(row_void, np.ma.mvoid)
        else:
            assert isinstance(row_void, np.void)
        assert row_void['a'] == 1
        assert row_void['b'] == 4

        # Confirm row is a view of table but row_void is not.
        table['a'][0] = -100
        assert row['a'] == -100
        assert row_void['a'] == 1

        # Make sure it works for a table that has masked elements
        if table.masked:
            table['a'].mask = True

            # row_void is not a view, need to re-make
            assert row_void['a'] == 1
            row_void = row.as_void()  # but row is a view
            assert row['a'] is np.ma.masked

    def test_row_and_as_void_with_objects(self, table_types):
        """Test the deprecated data property and as_void() method"""
        t = table_types.Table([[{'a': 1}, {'b': 2}]], names=('a',))
        assert t[0][0] == {'a': 1}
        assert t[0]['a'] == {'a': 1}
        assert t[0].as_void()[0] == {'a': 1}
        assert t[0].as_void()['a'] == {'a': 1}

    def test_bounds_checking(self, table_types):
        """Row gives index error upon creation for out-of-bounds index"""
        self._setup(table_types)
        for ibad in (-5, -4, 3, 4):
            with pytest.raises(IndexError):
                self.t[ibad]


text:
import networkx as nx
import numpy as np
import math
from tqdm import tqdm
import numba
from numba.experimental import jitclass
from numba import jit

steadyspec = [
    ('adj_matrix',numba.float64[:,:]),
    ('graph_size',numba.int32),
    ('background_field',numba.float64[:]),
    ('fixed_point_iter',numba.int32),
    ('fp_tol_fac',numba.float64),
]

@jitclass(steadyspec)
class steady_state(object):
    """
    A class to calculate steady state of magnetisation.

    ...

    Attributes
    ----------
    adj_matrix : numpy.array
        Adjacency matrix of the graph
    fixed_point_iter : float, optional
        Max number of iterations used in self-consistency equations (default 1e5)
    fp_tol_fac : float, optional
        Tolerance factor in stoppping condition for consistency equations (default 1e-6)

    """

    def __init__(self,adj_matrix,fixed_point_iter=10000,fp_tol_fac=1e-6):

        self.adj_matrix = adj_matrix
        self.graph_size = self.adj_matrix.shape[0]
        self.fixed_point_iter=fixed_point_iter
        self.fp_tol_fac=fp_tol_fac

    def single_mag(self,i,m,beta,field):

        """
        Calculates magnetisation for a single node. Subfunction of magnetisation function.

        Parameters
        ------------
        i : int
            Index of the node in question.
        m : numpy.array
            magnetisation array for all nodes.
        beta: float
            Interaction strength
        field: numpy.array
            Array of agent's control field for each node

        """

        gamma=1.0
        spin_field = self.adj_matrix[i].dot(m)
        term = math.tanh(beta*(spin_field+field[i]))
        return (1.0-gamma)*m[i] + gamma*term

    def magnetisation(self,mag,beta,field):

        """
        Calculates magnetisation for the whole system

        Parameters
        ------------
        m : numpy.array
            magnetisation array for all nodes.
        beta: float
            Interaction strength
        field: numpy.array
            Array of agent's control field for each node

        """

        m_old = mag
        m_new = np.zeros(len(m_old))
        for i in range(self.graph_size):
            m_new[i]=self.single_mag(i,m_old,beta,field)
        return m_new

    def aitken_method(self,mag0,beta,field):

        """
        Solves self-consistency equation by following Aitken method* for accelerating convergence.
        * Numerical Analysis Richard L.Burden 9th Edition, p. 105

        Parameters
        ------------
        m0 : numpy.array
            Initial guess of magnetisation for all nodes.
        beta: float
            Interaction strength
        field: numpy.array
            Array of agent's control field for each node

        """

        mag1=self.magnetisation(mag0,beta,field)
        for i in range(self.fixed_point_iter):
            mag2=self.magnetisation(mag1,beta,field)
            if ((mag0+mag2-2*mag1)!=0).all():
                mag_d = mag0 - (mag1-mag0)**2/(mag0+mag2-2*mag1)
            else:
                mag_d = mag1

            if abs(np.sum(mag0)-np.sum(mag_d))<self.fp_tol_fac:
                break
            mag0=mag1
            mag1=mag2
            if i+1==self.fixed_point_iter:
                mag_d = mag1

        return mag_d

@jit(nopython=True)
def isclose(a,b):
    return abs(a-b) <= max(1e-9 * max(abs(a), abs(b)), 1e-5)

@jit(nopython=True)
def susc_grad(beta,mag,adj_matrix):

    """
    Calculates mean field susceptibility.

    Parameters
    ------------
    beta: float
        Interaction strength.
    mag : numpy.array
        Magnetisation array for all nodes.
    adj_matrix : numpy.array
        Adjacency matrix of the network.

    """
    D=np.identity(mag.shape[0])*np.array([(1-i**2) for i in mag])
    inv = np.linalg.inv(np.identity(mag.shape[0])-beta*D.dot(adj_matrix))
    susc_matrix = beta*inv.dot(D)
    gradient = np.sum(susc_matrix,axis=1).flatten()

    return gradient

def mag_grad(beta,mag,adj_matrix):

    """
    Calculates gradient of the magnetisation with respect to change in the external control field. Nominally mean field susceptibility.

    Parameters
    ------------
    beta : float
        Interaction strength.
    mag : numpy.array
        Magnetisation array for all nodes.
    adj_matrix : numpy.array
        Adjacency matrix of the network.

    """

    if np.all([isclose(i,j) for i,j in zip(mag,np.ones(mag.shape[0]))]):
        return np.zeros(len(mag))
    else:
        return susc_grad(beta,mag,adj_matrix)

@jit(nopython=True)
def projection_simplex_sort(v, z):

    """
    Bounds control field to agent's magnetic field budget.
    ...

    Parameters
    ----------
    v : numpy.array
        Control field allocation of the agent.
    z : float
        Magnetic field budget.

    """

    n_features = v.shape[0]
    v = np.abs(v)
    u = np.sort(v)[::-1]
    cssv = np.cumsum(u) - z
    ind = np.arange(n_features) + 1
    cond = u - cssv / ind > 0
    rho = ind[cond][-1]
    theta = cssv[cond][-1] / float(rho)
    w = np.maximum(v - theta, 0)
    return w

@jit(nopython=True)
def lr_1(x,iim_iter=5000):
    return np.exp(-x/(0.1*iim_iter))

@jit(nopython=True)
def lr_2(x,iim_iter=5000):
    return np.exp(-x/(0.1*iim_iter))

@jit(nopython=True)
def adam(grad,it,typ,ms,vs,iim_iter,beta1=0.9,beta2=0.999,eps=0.1):

    """
        Adam optimiser.

        Parameters
        ------------
        grad : numpy.array
            Gradient array for all nodes.
        it : int
            Iteration index.
        typ : string
            'pos' if calculating for positive agent; 'neg' if calculating for negative agent.
        ms : numpy.array
            Momentum array with linear gradients
        vs : numpy.array
            Momentum array with squared gradients
        iim_iter : int
            Max number of iteration in optimisation algorithm.
        beta1 : float
            Default 0.9
        beta2 : float
            Default 0.999
        eps : float
            Default 0.1

    """

    if typ=='pos':
        lr=lr_1(it)
    elif typ=='neg':
        lr=lr_2(it)
    ms_new = beta1 * ms + (1.0 - beta1) * grad
    vs_new = beta2 * vs + (1.0 - beta2) * grad**2
    mhat = ms_new / (1.0 - beta1**(it+1))
    vhat = vs_new / (1.0 - beta2**(it+1))
    ms = ms_new
    vs = vs_new
    change = lr* mhat/(np.sqrt(vhat) + eps)
    return change,ms,vs

class mf_ising_system():

    """
    A class to calculate mean field Ising Maximisation Influence problem for a two agents.

    ...

    Attributes
    ----------
    graph : networkx graph
        undirected graph
    background_field : numpy.array
        Background field applied to the system
    iim_iter: float, optional
        Number of gradient ascent iterations (default 1000)
    iim_tol_fac: float, optional
        Tolerance factor in stopping condition for gradient ascent algorithm (default 1e-5)
    """

    def __init__(self,graph,background_field,iim_iter=10000,iim_tol_fac=1e-3):

        self.graph=graph
        self.adj_matrix = nx.to_numpy_matrix(graph).astype(np.float64)
        self.graph_size = self.adj_matrix.shape[0]
        self.background_field = background_field.astype(np.float64)
        self.iim_iter=iim_iter
        self.iim_tol_fac = iim_tol_fac

    def positive_agent(self,mag_i,it,pos_budget,beta,neg_change):

        """
        Single move by the positive agent.

        ...

        Parameters
        ----------
        mag_i : numpy.array
            Magnetisation array for all nodes.
        it : int
            Iteration of the algorithm.
        pos_budget : float
            Magnetic field budget for the positive agent.
        beta: float
            Interaction strength
        neg_change : numpy.array
            Change in allocation of the negative agent.

        """

        mag_i_grad = mag_grad(beta,mag_i,self.adj_matrix)
        control_field = self.control_field_history_pos[-1]
        change,ms,vs = adam(mag_i_grad,it,'pos',self.ms_pos,self.vs_pos,self.iim_iter)
        self.ms_pos=ms
        self.vs_pos=vs
        control_field_update = control_field  + (change+neg_change)
        control_field_new = projection_simplex_sort(control_field_update.T,pos_budget)
        self.control_field_history_pos.append(control_field_new)
        return control_field_new,mag_i_grad

    def negative_agent(self,mag_i,it,neg_budget,beta):

        """
        Single move by the negative agent.

        ...

        Parameters
        ----------
        mag_i : numpy.array
            Magnetisation array for all nodes.
        it : int
            Iteration of the algorithm.
        neg_budget : float
            Magnetic field budget for the negative agent.
        beta: float
            Interaction strength

        """
        mag_i = -1.0*mag_i
        mag_i_grad = -mag_grad(beta,mag_i,self.adj_matrix)
        control_field = self.control_field_history_neg[-1]
        change,ms,vs=adam(mag_i_grad,it,'neg',self.ms_neg,self.vs_neg,self.iim_iter)
        self.ms_neg=ms
        self.vs_neg=vs
        control_field_update = control_field - change
        control_field_new = projection_simplex_sort(control_field_update.T,neg_budget)
        self.control_field_history_neg.append(control_field_new)

        return control_field_new,mag_i_grad,change

    def second_partial_dffs(self,state,mag_ii,tot_field,beta,a=1e-5):

        """
        Calculates 2nd gradients of magnetisation with respect to each agents' control field.
        Calculated using central difference formula.

        ...

        Parameters
        ----------

        state : class instance
            steady_state class instance
        mag_ii : numpy.array
            Magnetisation array for all nodes.
        tot_field : numpy.array
            Total net magnetic field experienced by each node.
        beta: float
            Interaction strength
        a : float, optional
            Used in central difference formula. Specifies magnitude of change of control field.

        """

        update = a*np.ones(self.graph_size)
        upper_change=tot_field+update
        mag_plus= -state.aitken_method(mag_ii,beta,upper_change)
        grad_plus = -mag_grad(beta,mag_plus,self.adj_matrix)

        lower_change = tot_field-update
        mag_minus= -state.aitken_method(mag_ii,beta,lower_change)
        grad_minus = -mag_grad(beta,mag_minus,self.adj_matrix)
        second_total_grad =  (grad_plus - grad_minus)/(2*update) # central difference formula
        curv_player_neg = - second_total_grad # minus because product rule : H_pos = H_pos - H_neg
        curv_player_pos = curv_player_neg
        return np.array([curv_player_pos,curv_player_neg])

    def init_lists(self):
        """
        Initialises lists for storing variables.

        """

        self.control_field_history_pos =[]
        self.control_field_history_neg = []
        self.mag_history = []
        self.pos_gradient_history=np.zeros((self.iim_iter,self.graph_size))
        self.neg_gradient_history=np.zeros((self.iim_iter,self.graph_size))

        self.ms_pos = np.zeros(self.graph_size,dtype=np.float64)
        self.vs_pos = np.zeros(self.graph_size,dtype=np.float64)
        self.ms_neg = np.zeros(self.graph_size,dtype=np.float64)
        self.vs_neg = np.zeros(self.graph_size,dtype=np.float64)

    def MF_IIM(self,pos_budget,neg_budget,beta,init_alloc='random',progress=True):

        """
        Calculates competitive MF-IIM by following stochastic gradient ascent optimisation with
        Adam optimiser.

        Parameters
        ------------
        pos_budget : float
            Maximum magnetic field budget to be spent by the positive agent.
        neg_budget : float
            Maximum magnetic field budget to be spent by the negative agent.
        beta : float
            Interaction strength
        init_alloc : string or numpy.array, optional
            Either 'uniform' which corresponds to uniform spread of financial budget equaly among nodes.
            'random' corresponds to random initialisations. Alternatively, provide custom numpy.array
            allocation of your own choice. Default 'random'.
        progress : boolean
            If True shows progress bar; False otherwise.

        Outputs
        -----------
        control_field_pos : numpy.array
            Positive agent's control field allocation that results in the equilibrium.
        control_field_neg : numpy.array
            Negative agent's control field allocation that results in the equilibrium.
        final_mag : numpy.array
            Final magnetisation of the system.

        """

        if isinstance(init_alloc,(np.ndarray, np.generic)):
            control_field_pos = init_alloc[0,:]
            control_field_neg = init_alloc[1,:]
        elif isinstance(init_alloc,str):
            if  init_alloc=='aligned':
                control_field_pos =( pos_budget /self.graph_size)*np.ones(self.graph_size)
                control_field_neg = ( neg_budget /self.graph_size)*np.ones(self.graph_size)
            elif init_alloc=='random':
                control_field_pos  = np.random.dirichlet(np.ones(self.graph_size))*pos_budget
                control_field_neg  = np.random.dirichlet(np.ones(self.graph_size))*neg_budget

        init_mag = np.array([np.random.choice([-1,1]) for i in range(self.graph_size)]).astype(np.float64)

        self.init_lists()

        self.control_field_history_pos.append(control_field_pos)
        self.control_field_history_neg.append(control_field_neg)

        state = steady_state(self.adj_matrix)
        self.state=state

        tot_field = control_field_pos-control_field_neg

        mag_i = state.aitken_method(init_mag,beta,tot_field)
        self.mag_history.append(mag_i)

        for it in tqdm(range(self.iim_iter)) if progress else range(self.iim_iter):

            gradients=[]

            if neg_budget!=0:
                control_neg,neg_gradient,neg_change = self.negative_agent(mag_i,it,neg_budget,beta)
                tot_field-=control_neg
                self.neg_gradient_history[it]=neg_gradient
                gradients.append(neg_gradient)

            if pos_budget!=0:
                control_pos,pos_gradient = self.positive_agent(mag_i,it,pos_budget,beta,neg_change)
                tot_field +=control_pos
                self.pos_gradient_history[it]=pos_gradient
                gradients.append(pos_gradient)

            mag_ii= state.aitken_method(mag_i,beta,tot_field)

            self.mag_history.append(mag_ii)

            if np.all([(abs(gradient)<self.iim_tol_fac).all() for gradient in gradients]):
                second_dffs=self.second_partial_dffs(state,mag_ii,tot_field,beta)
                if (second_dffs[0]<0).all() and (second_dffs[1]<0).all():
                    break
            mag_i=mag_ii
            tot_field=0.0
        if it==self.iim_iter-1:
            print('Failed to converge after {} iterations'.format(self.iim_iter))
            final_mag = mag_ii

        elif it < self.iim_iter-1:
            final_mag = mag_ii

        self.control_field_history_pos = np.array(self.control_field_history_pos)
        self.control_field_history_neg = np.array(self.control_field_history_neg)
        self.mag_history = np.array(self.mag_history)

        return self.control_field_history_pos[-1],self.control_field_history_neg[-1],final_mag

max_stars_repo_name: kikatuso/Master-Thesis

max_issues_repo_name: kikatuso/Master-Thesis

max_forks_repo_name: kikatuso/Master-Thesis

reason: import numpy,import numba,from numba,import networkx