text:
from spectrl.main.learning import ProductMDP, Resource_Model, learn_policy, HyperParams
# from spectrl.hierarchy.construction import automaton_graph_from_spec
# from spectrl.hierarchy.reachability import HierarchicalPolicy, ConstrainedEnv
from spectrl.main.spec_compiler import ev, seq, choose, alw
from spectrl.util.io import parse_command_line_options, save_log_info
from spectrl.util.rl import print_performance, get_rollout
# from spectrl.rl.ars import HyperParams
from numpy import linalg as LA

import os
import numpy as np

from spectrl.examples.rooms_envs import GRID_PARAMS_LIST, START_ROOM, FINAL_ROOM
from spectrl.envs.rooms import RoomsEnv

# Define the resource model
# Fuel consumption proportional to distance from x-axis and the velocity
# sys_state: np.array(2)
# res_state: np.array(1)
# sys_action: np.array(2)
def resource_delta(sys_state, res_state, sys_action):
    return np.array([res_state[0] - 0.1 * abs(sys_state[0]) * LA.norm(sys_action)])

lb = [10., 20., 10., 10., 10., 9., 9., 9., 9., 9.]

fuel = [False, False, False, False, False, False, False, False, False]

# Construct Product MDP and learn policy
if __name__ == '__main__':
    flags = parse_command_line_options()
    render = flags['render']
    env_num = flags['env_num']
    folder = flags['folder']
    itno = flags['itno']
    spec_num = flags['spec_num']

    grid_params = GRID_PARAMS_LIST[env_num]

    print('\n**** Learning Policy for Spec #{} in Env #{} ****'.format(spec_num, env_num))

    # No. of steps for 3*3 per spec. 20 steps for each edge on a path.

    stepsrequired = [50, 50, 50, 50, 50, 70, 50, 70]
    # Step 1: initialize system environment
    # system = RoomsEnv(grid_params, START_ROOM[env_num],
    #                  FINAL_ROOM[env_num], max_timesteps=MAX_TIMESTEPS[env_num])

    system = RoomsEnv(grid_params, START_ROOM[env_num],
                      FINAL_ROOM[env_num], max_timesteps=stepsrequired[spec_num])

    # Step 2 (optional): construct resource model
    resource = Resource_Model(2, 2, 1, np.array([7.0]), resource_delta)

    # Step 3: List of specs.
    if env_num == 2:
        bottomright = (0, 2)
        topleft = (2, 0)
    if env_num == 3 or env_num == 4:
        bottomright = (0, 3)
        topleft = (3, 0)

    # test specs
    spec0 = ev(grid_params.in_room(FINAL_ROOM[env_num]))
    spec1 = seq(ev(grid_params.in_room(FINAL_ROOM[env_num])), ev(
        grid_params.in_room(START_ROOM[env_num])))
    spec2 = ev(grid_params.in_room(topleft))

    # Goto destination, return to initial
    spec3 = seq(ev(grid_params.in_room(topleft)), ev(grid_params.in_room(START_ROOM[env_num])))
    # Choose between top-right and bottom-left blocks (Same difficulty - learns 3/4 edges)
    spec4 = choose(ev(grid_params.in_room(bottomright)),
                   ev(grid_params.in_room(topleft)))
    # Choose between top-right and bottom-left, then go to Final state (top-right).
    # Only one path is possible (learns 5/5 edges. Should have a bad edge)
    spec5 = seq(choose(ev(grid_params.in_room(bottomright)),
                       ev(grid_params.in_room(topleft))),
                ev(grid_params.in_room(FINAL_ROOM[env_num])))
    # Add obsacle towards topleft
    spec6 = alw(grid_params.avoid_center((1, 0)), ev(grid_params.in_room(topleft)))
    # Either go to top-left or bottom-right. obstacle on the way to top-left.
    # Then, go to Final state. Only one route is possible
    spec7 = seq(choose(alw(grid_params.avoid_center((1, 0)), ev(grid_params.in_room(topleft))),
                       ev(grid_params.in_room(bottomright))),
                ev(grid_params.in_room(FINAL_ROOM[env_num])))

    specs = [spec0, spec1, spec2, spec3, spec4, spec5, spec6, spec7]

    # Step 4: construct product MDP
    env = ProductMDP(system, 2, specs[spec_num], 0.0,
                     lb[spec_num], res_model=resource) if fuel[spec_num] else ProductMDP(
                         system, 2, specs[spec_num], 0.0, lb[spec_num])

    # Step 5: Hyperparameters
    hparameterlist = [1000, 1000, 1000, 1000, 1000, 1500, 1000, 1000]
    hyperparams = HyperParams(30, 1, hparameterlist[spec_num], 30, 8, 0.05, 0.3, 0.1)

    # Step 6: Learn policy
    policy, log_info = learn_policy(env, hyperparams)

    # Save policy and log information
    logdir = os.path.join(folder, 'spec{}'.format(spec_num), 'spectrl')
    if not os.path.exists(logdir):
        os.makedirs(logdir)
    save_log_info(log_info, itno, logdir)
    # save_object('policy', policy, itno, logdir)

    # Print rollout and performance
    print_performance(env, policy)
    if render:
        get_rollout(env, policy, True)

max_stars_repo_name: luigiberducci/dirl

max_issues_repo_name: luigiberducci/dirl

max_forks_repo_name: luigiberducci/dirl


text:
# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import csv
import numpy as np
import os
import sys

from observations.util import maybe_download_and_extract

def quarrels(path):
  """Statistics of Deadly Quarrels

  *The Statistics Of Deadly Quarrels* by Lewis Fry Richardson (1960) is
  one of the earlier attempts at quantification of historical conflict
  behavior.

  The data set contains 779 dyadic deadly quarrels that cover a time
  period from 1809 to 1949. A quarrel consists of one pair of
  belligerents, and is identified by its beginning date and magnitude (log
  10 of the number of deaths). Neither actor in a quarrel is identified by
  name.

  Because Richardson took a dyad of belligerents as his unit, a given war,
  such as World War I or World War II comprises multiple observations, for
  all pairs of belligerents. For example, there are forty-four pairs of
  belligerents coded for World War I.

  For each quarrel, the nominal variables include the type of quarrel, as
  well as political, cultural, and economic similarities and
  dissimilarities between the pair of combatants.

  A data frame with 779 observations on the following 84 variables.

  `ID `
      V84: Id sequence

  `year `
      V1: Begin date of quarrel

  `international `
      V2: Nation vs nation

  `colonial `
      V3: Nation vs colony

  `revolution `
      V4: Revolution or civil war

  `nat.grp `
      V5: Nation vs gp in other nation

  `grp.grpSame `
      V6: Grp vs grp (same nation)

  `grp.grpDif `
      V7: Grp vs grp (between nations)

  `numGroups `
      V8: Num grps agst which fighting

  `months `
      V9: Num months fighting

  `pairs `
      V10: Num pairs in whole matrix

  `monthsPairs `
      V11: Num mons for all in mtrx

  `logDeaths `
      V12: Log (killed) matrix

  `deaths `
      V13: Total killed for matrix

  `exchangeGoods `
      V14: Gp sent goods to other

  `obstacleGoods `
      V15: Gp puts obstacles to goods

  `intermarriageOK `
      V16: Present intermarriages

  `intermarriageBan `
      V17: Intermarriages banned

  `simBody `
      V18: Similar body characteristics

  `difBody `
      V19: Difference in body characteristics

  `simDress `
      V20: Similarity of customs (dress)

  `difDress `
      V21: Difference of customs (dress)

  `eqWealth `
      V22: Common level of wealth

  `difWealth `
      V23: Difference in wealth

  `simMariagCust `
      V24: Similar marriage cusomst

  `difMariagCust `
      V25: Different marriage customs

  `simRelig `
      V26: Similar religeon or philosophy of life

  `difRelig `
      V27: Religeon or philisophy felt different

  `philanthropy `
      V28: General philanthropy

  `restrictMigration `
      V29: Restricted immigrations

  `sameLanguage `
      V30: Common mother tongue

  `difLanguage `
      V31: Different languages

  `simArtSci `
      V32: Similar science, arts

  `travel `
      V33: Travel

  `ignorance `
      V34: Ignorant of other/both

  `simPersLiberty `
      V35: Personal liberty similar

  `difPersLiberty `
      V36: More personal liberty

  `sameGov `
      V37: Common government

  `sameGovYrs `
      V38: Years since common govt established

  `prevConflict `
      V39: Belligerents fought previously

  `prevConflictYrs `
      V40: Years since belligerents fought

  `chronicFighting `
      V41: Chronic figthing between belligerents

  `persFriendship `
      V42: Autocrats personal friends

  `persResentment `
      V43: Leaders personal resentment

  `difLegal `
      V44: Annoyingly different legal systems

  `nonintervention `
      V45: Policy of nonintervention

  `thirdParty `
      V46: Led by 3rd group to conflict

  `supportEnemy `
      V47: Supported others enemy

  `attackAlly `
      V48: Attacked ally of other

  `rivalsLand `
      V49: Rivals territory concess

  `rivalsTrade `
      V50: Rivals trade

  `churchPower `
      V51: Church civil power

  `noExtension `
      V52: Policy not extending ter

  `territory `
      V53: Desired territory

  `habitation `
      V54: Wanted habitation

  `minerals `
      V55: Desired minerals

  `StrongHold `
      V56: Wanted strategic stronghold

  `taxation `
      V57: Taxed other

  `loot `
      V58: Wanted loot

  `objectedWar `
      V59: Objected to war

  `enjoyFight `
      V60: Enjoyed fighting

  `pride `
      V61: Elated by strong pride

  `overpopulated `
      V62: Insufficient land for population

  `fightForPay `
      V63: Fought only for pay

  `joinWinner `
      V64: Desired to join winners

  `otherDesiredWar `
      V65: Quarrel desired by other

  `propaganda3rd `
      V66: Issued of propaganda to third parties

  `protection `
      V67: Offered protection

  `sympathy `
      V68: Sympathized under control

  `debt `
      V69: Owed money to others

  `prevAllies `
      V70: Had fought as allies

  `yearsAllies `
      V71: Years since fought as allies

  `intermingled `
      V72: Had intermingled on territory

  `interbreeding `
      V73: Interbreeding between groups

  `propadanda `
      V74: Issued propaganda to other group

  `orderedObey `
      V75: Ordered other to obey

  `commerceOther `
      V76: Commercial enterprises

  `feltStronger `
      V77: Felt stronger

  `competeIntellect `
      V78: Competed succesfully intellectual occ

  `insecureGovt `
      V79: Government insecure

  `prepWar `
      V80: Preparations for war

  `RegionalError `
      V81: Regional error measure

  `CasualtyError `
      V82: Casualty error measure

  `Auxiliaries `
      V83: Auxiliaries in service of nation at war

  http://www.icpsr.umich.edu/icpsrweb/ICPSR/studies/05407

  Args:

    path: str.
      Path to directory which either stores file or otherwise file will
      be downloaded and extracted there.
      Filename is `quarrels.csv`.

  Returns:

    Tuple of np.ndarray `x_train` with 779 rows and 84 columns and
    dictionary `metadata` of column headers (feature names).
  """
  import pandas as pd
  path = os.path.expanduser(path)
  filename = 'quarrels.csv'
  if not os.path.exists(os.path.join(path, filename)):
    url = 'http://dustintran.com/data/r/HistData/Quarrels.csv'
    maybe_download_and_extract(path, url,
                               save_file_name='quarrels.csv',
                               resume=False)

  data = pd.read_csv(os.path.join(path, filename), index_col=0,
                     parse_dates=True)
  x_train = data.values
  metadata = {'columns': data.columns}
  return x_train, metadata

ext: py

reason: import numpy


text:
#!/usr/bin/python
"""
author  : Bhekimpilo Ndhlela
author  : 18998712
module  : Applied Mathematics(Numerical Methods) TW324
task    : computer assignment 04 question 3(bi) & (bii)
since   : Friday-23-03-2018
"""
def exact(steps, x=0.0, debug=False):
    fpp = lambda x, h : ((-f(x+2*h)+16*f(x+h)-30*f(x)+16*f(x-h)-f(x-2*h))/(12*h*h))
    f = lambda x : sqrt(1 - 2 * sin(x))
    exact = [abs(fpp(x, h) + 1.0) for h in steps]
    if debug is True:
        print "DEBUG MODE: [ON] QUESTION 3 bi.)"
        print "i","\t" ,"Step Size", "\t", "Exact Error"
        for i, (h, err) in enumerate(zip(steps, exact)):
            print (i+1), "\t","{:.7f}".format(h), "\t", "{:.10f}".format(err)
    return exact

def bound(steps, M=11.0, debug=False):
    machine_eps = finfo(float).eps
    bound = [(18.0 * machine_eps)/(12.0*(h**2)) + M*(h**4) for h in steps]
    if debug is True:
        print "DEBUG MODE: [ON] QUESTION 3 bii.)"
        print "i","\t" ,"Step Size", "\t", "Bound"
        for i, (h, err) in enumerate(zip(steps, bound)):
            print (i+1), "\t","{:.7f}".format(h), "\t", "{:.10f}".format(err)
    return bound

def plot_error_functions(steps, exact, bound):
    #loglog plot to display the error as function of the step size
    plt.title("Plot of the Exact Error and Bound as h(Step size) Changes")
    plt.xlabel("h (Step Size)")
    plt.ylabel("Exact vs Bound")
    plt.yscale('log')
    plt.xscale('log')
    plt.plot(steps, exact, "k-", label="Exact")
    plt.plot(steps, bound, "k--", label="Bound")
    plt.legend(bbox_to_anchor=(.65, .9))
    plt.show()

if __name__ == "__main__":
    from numpy import (abs, cos, sin, sqrt, logspace)
    from numpy import (finfo, float)
    import matplotlib.pyplot as plt
    from math import pow

    steps = logspace(-7, -1, num=100)
    exact, bound = exact(steps), bound(steps)
    plot_error_functions(steps, exact, bound)
else:
    from sys import exit
    exit("USAGE: python q3.py")

reason: from numpy


text:
#! /usr/bin/env python3
# -*- coding: utf-8 -*-

import re
import json
import argparse
import pandas as pd
import numpy as np
from datetime import datetime
import matplotlib.pyplot as plt

plt.style.use('seaborn-dark-palette')

if __name__ == '__main__':
  # parse cli and config
  parser = argparse.ArgumentParser()
  parser.add_argument('-i', '--input', type=str, required=True)
  parser.add_argument('-s', '--show', action='store_true')
  parser.add_argument('-b', '--bin', action='store_true')
  parser.add_argument('-t', '--dt', type=int, default=300)
  parser.add_argument('-r', '--range', type=str, default='')
  parser.add_argument('-d', '--dev', type=str, default='wifi')
  parser.add_argument('-f', '--filter', nargs='+', default=[])
  args = parser.parse_args()
  base = args.input[:args.input.find('_')]

  dt_fmt = '%Y%m%d-%H%M%S'
  freq = f'{args.dt}s'
  if args.range == '':
    tok = args.input[:args.input.rfind('.')].split('_')
    try:
      start = datetime.strptime(tok[-2], dt_fmt)
      stop = datetime.strptime(tok[-1], dt_fmt)
    except:
      start = datetime.strptime(tok[-3], dt_fmt)
      stop = datetime.strptime(tok[-2], dt_fmt)
  else:
    start = datetime.strptime(args.range.split('|')[0], dt_fmt)
    stop = datetime.strptime(args.range.split('|')[1], dt_fmt)
  base = f'{base}_{start.strftime(dt_fmt)}_{stop.strftime(dt_fmt)}_{args.dev}'

  df = pd.read_csv(args.input, sep=';')
  try:
    df = df.rename(columns={'mac-address':'mac_address'})
  except:
    pass
  df.date_time = pd.to_datetime(df.date_time)
  df.date_time = df.date_time.dt.tz_localize(None)
  df = df[ (df.date_time >= start) & (df.date_time < stop) ]
  df = df[ df.kind == args.dev ]

  stats = pd.DataFrame()
  t_index = pd.date_range(start=start, end=stop, freq=freq)
  for sid, dfg in df.groupby(['station_name']):
    print(f'{args.dev}) {sid}  -> {dfg.shape}')

    dfg = dfg.set_index('date_time')
    dfg.index = pd.to_datetime(dfg.index)
    dfr = dfg[['mac_address']].resample(freq).count()
    dfr.columns = [f'{sid}']

    s = pd.Series(dfg['mac_address'], index=dfg.index)
    dfu = pd.DataFrame(s.groupby(pd.Grouper(freq=freq)).value_counts())
    dfu.columns = ['repetitions_counter']
    dfu = dfu.reset_index()
    dfu = dfu.set_index('date_time')
    dfu = dfu.groupby('date_time')[['mac_address']].count()
    dfu.columns = [f'{sid}_unique']
    #print('dfu', dfu)
    #idx = [ i for i in dfu.index ]
    #print('my',idx)
    #print(dfu.index)
    #dfu = pd.DataFrame(idx, columns=[f'{sid}_unique'], index=dfu.index)#, index=[ i for idfu.index.get_level_values(0)] )
    #print(dfu)

    if len(stats) == 0:
      stats = dfr
    else:
      stats = pd.concat([stats, dfr], axis=1)
    stats = pd.concat([stats, dfu], axis=1)

  stats = stats.fillna(0)

  # plot
  w, h, d = 12, 7, 150
  fig, ax = plt.subplots(1, 1, figsize=(w, h), dpi=d)
  plt.suptitle(f'Device type {args.dev}, unique device count')

  print(stats.columns)
  print(stats)
  for cid in stats.columns:

    if len(args.filter):
      typetag = 'sel_'
      print(cid)
      check = np.array([ re.match(f'.*{f}.*', cid) != None for f in args.filter ])
      print(check)
      print(check.any())
      if not check.any():
        continue
    else:
      typetag = ''

    try:
      ntag, lbl = cid.split('-')
    except:
      ntag = lbl = cid

    #if not ntag in ['(1)', '(5)', '(6)']: continue
    print(cid)
    print(stats[cid])
    if cid.endswith('_unique'):
      ax.plot(stats.index, stats[cid], '-', label=lbl)
    else:
      continue

  ax.set_title(f'period {start} -> {stop}, resampling {freq}')
  ax.legend()
  ax.grid()
  ax.tick_params(labelrotation=45)

  if args.show:
    plt.show()
  else:
    plt.savefig(f'{base}_presence_{typetag}{freq}.png')
  plt.close()

max_stars_repo_name: physycom/slides

max_issues_repo_name: physycom/slides

max_forks_repo_name: physycom/slides


text:
# -----------------------------------------------------------------------------
# Copyright (c) 2019 Nicolas P. Rougier
# Distributed under the terms of the BSD License.
# -----------------------------------------------------------------------------
import numpy as np
import os, sys, struct

# Read the MNIST dataset (training or testing)
def read(dataset="training", path="MNIST"):
    if dataset == "training":
        filename_img = os.path.join(path, 'train-images-idx3-ubyte')
        filename_lbl = os.path.join(path, 'train-labels-idx1-ubyte')
    elif dataset == "testing":
        filename_img = os.path.join(path, 't10k-images-idx3-ubyte')
        filename_lbl = os.path.join(path, 't10k-labels-idx1-ubyte')
    else:
        raise ValueError("dataset must be 'testing' or 'training'")

    with open(filename_lbl, 'rb') as file:
        magic, count = struct.unpack(">II", file.read(8))
        labels = np.fromfile(file, dtype=np.int8)
        Y = np.zeros((len(labels), 10))
        for i in range(len(labels)): Y[i,labels[i]] = 1
        labels = Y

    with open(filename_img, 'rb') as file:
        magic, count, rows, cols = struct.unpack(">IIII", file.read(16))
        images = np.fromfile(file, dtype=np.uint8)
        images = images.reshape(count, rows, cols)
        images = (images-images.min())/(images.max()-images.min())

    return images, labels

max_stars_repo_name: rougier/VSOM

max_issues_repo_name: rougier/VSOM

max_forks_repo_name: rougier/VSOM


text:
# goals:

# prepre sac --> hdf5 file, estimating how much space it'll take
# and allowing for station selection / time frame selection
#
#
# first choose station, then choose time frame
#
# want to have tools to show uptime for the station

# which suggests having different scripts or just like

# multi_station -up -o save.png -sta TA01-TA19.txt
# to show graphically or sth -st

# multi_station -make -s input.txt -o save_folder

# multi_station -choose sac_folders -o input.txt

import matplotlib

import matplotlib.pyplot as plt

from subprocess import check_output
import argparse
from pathlib import Path
import datetime
import pandas as pd
import numpy as np
import os

import json

#get_all_files("/home/eos_data/SEISMIC_DATA_TECTONICS/RAW/ACEH/MSEED/")

def get_all_files(sac_folders, output_file):
	# look inside /tgo/SEISMIC_DATA_TECTONICS/RAW/ACEH/MSEED/
	# for .SAC files

	df_list = []

	print(sac_folders)
	for sac_folder in sac_folders:

		all_files = [str(p) for p in Path(sac_folder).rglob("*SAC")]

		df = pd.DataFrame()

		df["filepath"] = all_files

		#print(df)

		for index, row in df.iterrows():
			_sta = row.filepath.split("/")[-2] # second last entry should be the station name

			# the sac files have their channel in front and not behind so...

			# check if all 3 channels are available too?

			_file = row.filepath.split("/")[-1]

			_year = _file.split(".")[5]
			_jday = _file.split(".")[6]

			out = check_output(["saclst", "KZDATE", "KZTIME", "B", "E", "f", row.filepath])
			out = [x for x in out.decode('UTF-8').strip().split(" ") if x != ""]

			print(out)

			try:
				assert len(out) > 2
			except:
				continue

			df.at[index, 'station'] = _sta
			df.at[index, 'year'] = (_year) # it's saved as a string; so pandas probably inferred that it's an int
			df.at[index, 'jday'] = (_jday)
			df.at[index, 'start_time'] = _file.split(".")[7]

			df.at[index, "kzdate"] = out[1]
			df.at[index, "kztime"] = out[2]
			df.at[index, "B"] = out[3]
			df.at[index, "E"] = out[4]

			df.at[index, "sac_start_dt"] = datetime.datetime.strptime("{} {}".format(out[1], out[2]), "%Y/%m/%d %H:%M:%S.%f")
			df.at[index, "sac_end_dt"] = datetime.datetime.strptime("{} {}".format(out[1], out[2]), "%Y/%m/%d %H:%M:%S.%f") + datetime.timedelta(seconds = float(out[4]))

		df_list.append(df)

		# fullday doesn't give useful info because the file could start at 000000 and still be incomplete

	new_df = pd.concat(df_list)
	print(new_df)
	new_df.to_csv(output_file, index = False)

def generate_timestamps(input_csv, output_csv):
	df = pd.read_csv(input_csv)
	df["sac_start_dt"] = pd.to_datetime(df["sac_start_dt"])
	df["sac_end_dt"] = pd.to_datetime(df["sac_end_dt"])

	overlap = 0.3

	c_df = pd.DataFrame()

	for index, row in df.iterrows():

		df.at[index, "c"] = ".EHZ." in row.filepath

	df = df[df["c"] == True]

	# do for all stations by default
	gdf = df.groupby(by = "station")

	c = 0

	of = pd.DataFrame() # output dataframe

	for sta, _df in gdf:

		# there seems to be an O(n^2) slow down which is probably because the dataframe is getting so long

		timestamp_set = []
		print(sta)

		for index, row in _df.iterrows():
			new_timestamp = (row.sac_start_dt, row.sac_end_dt, row.filepath)

			#print("new timestamp", new_timestamp)
			#print(timestamp_set)

			if len(timestamp_set) == 0:
				timestamp_set.append(new_timestamp)

			else:
				overlap_set = []
				# for ~ a few months of data with station-days, an O(n^2) search is acceptable
				for test_ts in timestamp_set:
					if new_timestamp[0] < test_ts[1] and new_timestamp[1] > test_ts[0]:
						overlap_set.append(test_ts) # check for overlap, collate set of overlaps

				timestamp_set.append(new_timestamp)

		timestamp_set = sorted(timestamp_set, key = lambda x: x[0])

		for i in timestamp_set:
			print(i)

		for _ts in timestamp_set:

			print(_ts)
			# generate timestamps for each date range into the output dataframe
			# start time of 1 min, station, filepath, sac_start_time

			n_cuts = ((_ts[1] - _ts[0]).total_seconds() - (overlap * 60))/((1 - overlap)*60)

			n_cuts = int(np.floor(n_cuts))

			for i in range(n_cuts):
				_ts_start = _ts[0] + datetime.timedelta(seconds = (i * (1 - overlap) * 60))

				_dt = (1 - overlap) * 60 * i

				of.at[c, "station"] = sta
				of.at[c, "timestamp_start"] = _ts_start
				of.at[c, "time_from_file_start"] = _dt
				of.at[c, "filepath"] = _ts[2]

				c += 1

		print(of)

		of.to_csv("test.csv")

		break

def plot_all_uptime(selector_file, start_date, end_date, all_csv_path = "station/all_aceh_sac.csv"):

	# selector file: list of stations

	with open(selector_file, 'r') as f:
		station_list = f.read().split("\n")

	station_list = list(filter(lambda x: x != "", station_list))

	n_stations = len(station_list)

	_parse_char = "j"

	start_date = datetime.datetime.strptime(start_date, "%Y.%{}".format(_parse_char))
	end_date = datetime.datetime.strptime(end_date, "%Y.%{}".format(_parse_char))

	n_days = (end_date - start_date).days + 1

	#print(n_stations)
	#print(n_days)

	image = np.zeros((n_stations, n_days))
	#print(image.shape)

	df = pd.read_csv(all_csv_path)

	# for index, row in df.iterrows():
	# 	df.at[index, 'datetime'] = datetime.datetime.strptime("{}_{}".format(row.year, row.jday), "%Y_%j")

	# 	for _cha in ["EHE", "EHN", "EHZ"]:
	# 		if _cha in row.filepath:
	# 			df.at[index,'channel'] = _cha
	#

	if not 'dt' in df.columns:

		for index, row in df.iterrows():
			df.at[index, 'dt'] = datetime.datetime.strptime("{}.{}".format(row['year'], row['jday']), "%Y.%j")

	#print(df.datetime)

	for index, row in df.iterrows():

		if row.station in station_list and row['dt'] >= start_date and row['dt'] <= end_date:
			station_index = station_list.index(row.station)
			# print(row["dt"])
			# print(start_date)
			day_index = (row["dt"] - start_date).days

			image[station_index, day_index] = 1

			#print(station_index, day_index)

	print("##########################################\n")
	for i in range(n_stations):
		print(station_list[i], "{:.2f}".format(np.sum(image[i, :])/n_days))
	print("\n##########################################")

	#print()

	plt.figure(figsize=(18,12), dpi = 150)
	plt.yticks(np.arange(n_stations) + 0.5, list(station_list), fontsize = 8)
	plt.xticks(np.arange(n_days) + 0.5, np.arange(0, (n_days)), fontsize = 8)
	plt.xlabel("Days")
	plt.ylabel("Station name")
	plt.pcolormesh(image, edgecolors ='k', linewidth=2)
	plt.savefig("log/uptime/uptime_{}_{}-{}.png".format(datetime.datetime.strftime(datetime.datetime.now(), "%Y%m%d-%H%M%S"), datetime.datetime.strftime(start_date, "%Y.%j"), datetime.datetime.strftime(end_date, "%Y.%j")))

def select_files(selector_file,  all_csv_path = "station/all_aceh_sac.csv", output_file = "", patch = False):

	# start_date format = e.g. 2020_03 for month

	df = pd.read_csv(all_csv_path)

	# kinda inefficient bc 10^5 rows but it's fine bc it won't be used very often

	for index, row in df.iterrows():
		df.at[index, 'dt'] = datetime.datetime.strptime("{}.{}".format(row.year, row.jday), "%Y.%j")

		df.at[index, "_uid"] = "{}-{}-{}-{}".format(row.station, row.year, row.jday, row.start_time)

		for _cha in ["EHE", "EHN", "EHZ"]:
			if _cha in row.filepath:
				df.at[index,'channel'] = _cha

	#station_list = []

	with open(selector_file, 'r') as f:
		station_list = f.read().split("\n")

	station_list = list(filter(lambda x: x != "", station_list))

	if patch:
		_df = df[(df["station"].isin(station_list)) & (df["start_time"] != "000000") & (df["start_time"].str[-2:] != "-1")]

		# -1 at the back: almost defn a duplicate
		# anything else: likely to be a partial day file

	else:
		_df = df[df["station"].isin(station_list)]

	_df.sort_values("jday", inplace = True)

	#_df.to_csv("station/test.csv")
	# want to check if it's complete, whether it's all fulldays, if any gaps
	# which should also show me the missing files

	# if len(_df.index) < expected_files:
	# 	print("some missing, can report on the no. of missing + flag to continue")

	# 	print("expected: ", expected_files, "actual: ", len(_df.index))
	# 	#plot_all_uptime(selector_file, _startdate, _enddate)

	# elif len(_df.index) > expected_files:
	# 	print("more files than expected which is odd, have to filter so that it's only 3")

	# 	print("This shouldn't happen")

	# 	#idk what to do if this happens actually

	# else:
	# 	print("all files present, no issues :)")

	if output_file:
		_df.to_csv(output_file, index = False)

def make_station_json(station_coords, station_list, output):

	station_json = {}

	with open(station_list, "r") as f:
		station_list = f.read().split("\n")[:-1]

	for _station in station_list:
		print(_station)
		station_json[_station] = {"network": "AC", "channels":["EHZ", "EHE", "EHN"]}

	with open(station_coords, "r") as f:
		coordinates = f.read().split("\n")

		coordinates = [y.strip() for x in coordinates if len(x) > 0 for y in x.strip().split("\t") if len(y) > 0 ]

	for station in station_json:
		print(station)
		i = coordinates.index(station)

		station_json[station]["coords"] = [float(coordinates[i + 2]), 100, float(coordinates[i+1])] # i set elevation to 100 because that's the average height in sumatra and i'm not given this information + it's not that important

	with open(output, 'w') as f:
		json.dump(station_json, f)

def pbs_writer(n_nodes, output_csv, job_name, env_name, project_code, no_execute = False, yes_execute = False, use_gpu = False):

	output_pbs = os.path.join("/home/zchoong001/cy1400/cy1400-eqt/pbs", job_name +".pbs")

	with open(output_pbs, "w") as f:
		f.write("#PBS -J 0-{}\n".format(n_nodes - 1))
		if use_gpu:
			f.write("#PBS -N {}\n#PBS -P {}\n#PBS -q gpu8\n#PBS -l select={}:ncpus=1:mpiprocs=32:ngpus=1\n".format(job_name, project_code, n_nodes))
		elif no_execute and not use_gpu:
			f.write("#PBS -N {}\n#PBS -P {}\n#PBS -q q32\n#PBS -l select={}:ncpus=1:mpiprocs=32:mem=2gb -l walltime=1:00:00\n".format(job_name, project_code, n_nodes))

		elif (yes_execute or not no_execute) and not use_gpu:
			f.write("#PBS -N {}\n#PBS -P {}\n#PBS -q q32\n#PBS -l select={}:ncpus=1:mpiprocs=32:mem=16gb -l walltime=80:00:00\n".format(job_name, project_code, n_nodes))
		f.write("#PBS -e log/pbs/{0}/error.log \n#PBS -o log/pbs/{0}/output.log\n".format(job_name))
		if use_gpu:
			f.write("module load cuda/10.1\nexport HDF5_USE_FILE_LOCKING='FALSE'\n")
		f.write("module load python/3/intel/2020\nmodule load sac\ncd $PBS_O_WORKDIR\nnprocs=`cat $PBS_NODEFILE|wc -l`\ninputfile=/home/zchoong001/cy1400/cy1400-eqt/node_distributor.py\n")
		f.write("encoded_file={}\nmkdir -p log/pbs/{}\nsource activate {}\n".format(output_csv, job_name, env_name))

		if not yes_execute:
			f.write("python $inputfile -id $PBS_ARRAY_INDEX -decode $encoded_file\n")

		if not no_execute:
			f.write("/home/zchoong001/cy1400/cy1400-eqt/pbs/runtime_scripts/{0}/${{PBS_ARRAY_INDEX}}.sh >& log/pbs/{0}/$PBS_JOBID.log 2>&1\n".format(job_name))

def read_config(config_file_path):
	"""
	config_file_path should be an absolute file path for consistency
	"""

	options = {}
	with open(config_file_path, 'r') as f:
		for line in f:
			k,v = line.strip().split("=")

			options[k] = v

	print(options)
	return options

def encode_multirun(
	output_csv = "",
	station_file = "",
	n_multi = 20,
	n_nodes = None,
	job_name = "",
	start_day = "",
	end_day = "",
	model_path = "",
	hdf5_parent = "",
	detection_parent = "",
	write_hdf5 = False,
	run_eqt = False,
	plot_eqt = False,
	merge_csv = False,
	recompute_snr = False,
	filter_csv = False,
	write_headers = False,
	pbs = None,
	sac_select = None,
	make_sac_csv = False,
	no_execute = False,
	yes_execute = False,
	snr_threshold = 8,
	config = "",
	patch = False,
	json_file = "",
	use_gpu = False,
	):

	# encode everything (sac to hdf5, prediction
	# node distributor will decode according/modify a main() function in one script to control what you want to do
	#
	# also, since eqt assumes that you join all the hdf5 files together and mash them, but i'm keeping them in
	# station-named folders, so that allows me to open one node per station
	# it also means that i'll need to keep everything in an additional layer of folder
	# it could be DDMMM_experiment/TA19/TA19/TA19.csv and TA19.hdf5
	#
	# which sounds kind of stupid but it'll work fine

	# need to know how many nodes i'm using
	# plus other arguments
	# could it be a json file so i can store other args?
	#
	# what about sac writing? Will you need that many nodes as well?
	#
	#
	# i think al lthe command line arguments could be put into this file so I can control what i want at a single point

	#MAX_NODES = 20

	#n_multi = 20
	# assign one node to one station for now

	#sac_selector = "9jun_10station.csv" # idt this is needed
	# station_file = "station/random10.txt" # this means that i'll have to check somewhere if the station has any data for that period...
	# might as well solve the problem at the source i.e. for each batch that you feed in (inside multi station ) check that all stations have more than 1 file

	# and if not they are automatically dropped / with some error message
	# or just throw an error? but that's kinda disruptive i wouldn't like it
	#
	if config:
		config_options = read_config(config)

	if job_name == "":
		job_name = datetime.datetime.strftime(datetime.datetime.now(), "%Y%m%d_%H%M%S")

	if config:
		model_path = config_options["MODEL_PATH"]
		station_json = config_options["STATION_JSON"]
		if hdf5_parent == "":
			hdf5_parent = os.path.join("/scratch", (config_options["USERNAME"]), job_name)

		if detection_parent == "":
			detection_parent = os.path.join(config_options["PROJECT_ROOT"],"detections", job_name)

		if not output_csv:
			output_csv = os.path.join(config_options["PROJECT_ROOT"], "node_encode", job_name + ".csv")

		if not sac_select:
			sac_select = os.path.join(config_options["PROJECT_ROOT"], "station_time", job_name + ".csv")

		_folders = [
		os.path.join(config_options["PROJECT_ROOT"], "detections"),
		os.path.join(config_options["PROJECT_ROOT"], "node_encode"),
		os.path.join(config_options["PROJECT_ROOT"], "station_time"),
		os.path.join(config_options["PROJECT_ROOT"], "pbs", "runtime_scripts"),
		os.path.join(config_options["PROJECT_ROOT"], "log", "timing"),]

		for _f in _folders:
			if not os.path.exists(_f):
				try:
					os.makedirs(_f)
				except:
					pass

		env_name = config_options["ENV_NAME"]
		project_code = config_options["PROJECT_CODE"]

	else:
		if not model_path:
			model_path = "/home/zchoong001/cy1400/cy1400-eqt/EQTransformer/ModelsAndSampleData/EqT_model.h5"

		if not json_file:
			station_json = "/home/zchoong001/cy1400/cy1400-eqt/station/json/all_stations.json"

		else:
			station_json = json_file

		if hdf5_parent == "":
			hdf5_parent = os.path.join("/scratch/zchoong001", job_name)

		if detection_parent == "":
			detection_parent = os.path.join("/home/zchoong001/cy1400/cy1400-eqt/detections", job_name)

		if not output_csv:
			output_csv = os.path.join("/home/zchoong001/cy1400/cy1400-eqt/node_encode", job_name + ".csv")

		if not sac_select:
			sac_select = os.path.join("/home/zchoong001/cy1400/cy1400-eqt/station_time", job_name + ".csv")

		env_name = "tf2"
		project_code = "eos_shjwei"

	if make_sac_csv:
		select_files(station_file, output_file = sac_select, all_csv_path = make_sac_csv, patch = patch)

	df = pd.DataFrame(columns = ["id", "sta", "hdf5_folder", "prediction_output_folder", "merge_output_folder", "start_day", "end_day", "multi", "model_path"])

	with open(station_file, 'r') as f:
		station_list = f.read().split("\n")

	station_list = list(filter(lambda x: x != "", station_list))

	if not n_nodes:
		n_nodes = len(station_list)

	# if len(station_list) > 20:
	# 	raise ValueError("Does not support more than {} stations, please split.".format(MAX_NODES))

	for c, sta in enumerate(station_list):

		df.at[c,"id"] = c + 1 # count from 1 bc the array number starts from 1? like it's just easier

		df.at[c, "sta"] = sta

		df.at[c, "hdf5_folder"] = os.path.join(hdf5_parent,sta, sta)

		df.at[c, "prediction_output_folder"] = os.path.join(detection_parent, sta) # multiruns will have like multi_01 behind etc.
		df.at[c, "model_path"] =  model_path
		df.at[c, "merge_output_folder"] = os.path.join(detection_parent, "{}_merged".format(sta))

		df.at[c, "start_day"] = start_day
		df.at[c, "end_day"] = end_day

		df.at[c, "multi"] = n_multi
		df.at[c, "nodes"] = n_nodes
		df.at[c, "snr_threshold"] = snr_threshold

		df.at[c, "station_json"] = station_json
		if config:
			df.at[c, "project_root"] = config_options["PROJECT_ROOT"]

			df.at[c, "output_folder"] = os.path.join(config_options["PROJECT_ROOT"], "pbs", "runtime_scripts", job_name)

		else:

			df.at[c, "project_root"] = "/home/zchoong001/cy1400/cy1400-eqt/"

			df.at[c, "output_folder"] = "/home/zchoong001/cy1400/cy1400-eqt/pbs/runtime_scripts/{}".format(job_name)

		df.at[c, "job_name"] = job_name
		df.at[c, "write_hdf5"] = write_hdf5
		df.at[c, "run_eqt"] = run_eqt
		df.at[c, "plot_eqt"] = plot_eqt

		df.at[c, "merge_csv"] = merge_csv
		df.at[c, "filter_csv"] = filter_csv
		df.at[c, "recompute_snr"] = recompute_snr
		df.at[c, "write_headers"] = write_headers

		df.at[c, "sac_select"] = sac_select

	df.to_csv(output_csv, index = False)

	if pbs:
		pbs_writer(n_nodes, output_csv, job_name, env_name, project_code, no_execute = no_execute, yes_execute = yes_execute, use_gpu = use_gpu)

	# load station_list, see number of stations

	# behind hte scenes, i will split the stations into groups of 20

	# if no. of stations > n_nodes, then make a new job? or can just plan this separately

if __name__ == "__main__":

	# i'm going to make this environment so cluttered with command line arguments
	# even i won't know how to use it

	parser = argparse.ArgumentParser(description = "utils for preparing multistation hdf5 files, running eqt (future) and plotting sac files")

	parser.add_argument("--get", help = "name of actual target sac folder (month)", default = None, nargs = '+')
	parser.add_argument("-i", "--input", help = "input file")
	parser.add_argument("-o", "--output", help = "output file", default = "")

	parser.add_argument("-sf", "--selector", help = "txt file with linebreak separated station names, specifying stations of interest")

	parser.add_argument("-s", "--start_date", help = "underscore separated year with julian day e.g. 2020_085 for start date (inclusive)")
	parser.add_argument("-e", "--end_date", help = "underscore separated year with julian day e.g. 2020_108 for enddate (inclusive)")
	parser.add_argument("-m", "--month", help = "flag to use month. default is Julian day. e.g. 2020_03 to represent March", action = "store_true", default = False)
	parser.add_argument("-j", "--julian", help = "default is True, to use Julian day to specify start and end date", action = "store_true", default = True)

	parser.add_argument("-js", "--json", help = "file with the coordinates of all the stations")

	parser.add_argument("-pl", "--plot", help = "get uptime file for some start and end date", action = "store_true")

	parser.add_argument("-enc", "--encode", action="store_true", help = "flag for script to encode a csv file for multinode running on HPC array")

	parser.add_argument("-job", help = "unique string identifier for the job", default = "")

	parser.add_argument("-write_hdf5", action = "store_true", help = "flag to write from sac to hdf5 folder, default False")
	parser.add_argument("-run_eqt", action = "store_true", help = "flag to run eqt prediction and merge multiple predictions, default False")
	parser.add_argument("-merge_csv", action = "store_true", help = "flag to merge multiple csvs")
	parser.add_argument("-recompute_snr", action = "store_true", help = "flag to recompute snrs")

	parser.add_argument("-filter_csv", action = "store_true", help = "flag to filter csvs")
	parser.add_argument("-write_headers", action = "store_true", help = "flag to write sac headers")

	parser.add_argument("-plot_eqt", action = "store_true", help = "flag to plot 150s sac traces and png 3C plots, default False")

	parser.add_argument("-det_parent", help = "parent folder to keep detection output files", default = "")
	parser.add_argument("-hdf5_parent", help = "parent folder to keep hdf5 files", default = "")
	parser.add_argument("-model_path", help = "path to model", default = "")
	parser.add_argument("-n_multi", type = int, help = "no. of times to repeat prediction, default 20", default = 20)
	parser.add_argument("-n_nodes", type = int, help = "no. of HPC nodes to use", default = None)

	parser.add_argument('-p', "--patch", action = "store_true")

	parser.add_argument("-pbs", help = "flag to generate pbs script", default = False, action = "store_true")

	parser.add_argument("-msc","--make_sac_csv", help = "table of all SAC files in the target folder, filtered by the station list (-i)", default = False)

	parser.add_argument("-ss", "--sac_select", help = "string for .csv file in station time folder, generated by select_files() in multi_station.py", default = None)
	parser.add_argument("-nx", "--no_execute", action = "store_true", help = "pbs file will generate but not run", default = False)
	parser.add_argument("-yx", "--yes_execute", action = "store_true", help = "script has been generated, run script without calling node_distributor", default = False)

	parser.add_argument("-snr_threshold", type = int, help = "threshold number for S wave SNR", default = 8)
	parser.add_argument("-config", help = "path to config file from which directories are created")

	parser.add_argument("-gts", "--generate_timestamps", action = "store_true")
	parser.add_argument("-gpu", "--use_gpu", action = "store_true")

	args = parser.parse_args()

	if args.selector and not args.plot:
		select_files(args.selector, args.start_date, args.end_date, args.julian, args.month, args.input, args.output)

	elif args.generate_timestamps:
		generate_timestamps(args.input, args.output)

	elif args.get:
		get_all_files(args.get, args.output)
	elif args.json and not args.encode:
		make_station_json(args.json, args.input, args.output)

	elif args.plot:
		plot_all_uptime(args.selector, args.start_date, args.end_date, args.input)
	elif args.encode:
		encode_multirun(
			output_csv = args.output,
			station_file = args.input,
			job_name = args.job,
			start_day = args.start_date,
			end_day = args.end_date,
			model_path = args.model_path,
			hdf5_parent = args.hdf5_parent,
			detection_parent = args.det_parent,
			write_hdf5 = args.write_hdf5,
			run_eqt = args.run_eqt,
			plot_eqt = args.plot_eqt,
			pbs = args.pbs,
			n_nodes = args.n_nodes,
			n_multi = args.n_multi,
			sac_select = args.sac_select,
			make_sac_csv = args.make_sac_csv,
			merge_csv = args.merge_csv,
			recompute_snr = args.recompute_snr,
			filter_csv = args.filter_csv,
			write_headers = args.write_headers,
			no_execute = args.no_execute,
			yes_execute = args.yes_execute,
			snr_threshold = args.snr_threshold,
			config = args.config,
			patch = args.patch,
			json_file = args.json,
			use_gpu = args.use_gpu)


text:
import cv2
import math
import numpy as np
import os
import pandas as pd
import subprocess
import sys
import time
import tkinter as tk
from gaze_tracking import GazeTracking

# create the root window
root = tk.Tk()
# SCREEN_WIDTH = root.winfo_screenwidth()
# SCREEN_HEIGHT = root.winfo_screenheight()
root.withdraw()

print(f'Initializing calibration.')
start_delay = 3  # seconds
duration = 12000  # milliseconds
print(f'Calibration will start in {start_delay} seconds.')
time.sleep(start_delay)

coords = []
time_results = []
res = 0
contador = 0
temp_ini = time.time()  # tempo que começa o programa
contadorRight = 0
contadorLeft = 0
menor_x = 5000
maior_x = -5000
lista_tempo = []
lista_gaze_x = []
lista_gaze_y = []
# gaze = GazeTracking()
# webcam = cv2.VideoCapture(0)

start_time = time.time()
print('Main loop started.')
x_min = 0
x_max = 0
y_min = 0
y_max = 0

window_name = "Gaze Tracking"
cv2.namedWindow(window_name, cv2.WND_PROP_FULLSCREEN)
cv2.moveWindow(window_name, 0, 0)
cv2.setWindowProperty(window_name, cv2.WND_PROP_FULLSCREEN, cv2.WINDOW_FULLSCREEN)

def calibragem():
    print("entrou na função calibragem")
    """
    calibra bla bla
    :return:
    """
    contador = 0
    gaze = GazeTracking()
    webcam = cv2.VideoCapture(0)

    SCREEN_WIDTH = int(webcam.get(cv2.CAP_PROP_FRAME_WIDTH))
    SCREEN_HEIGHT = int(webcam.get(cv2.CAP_PROP_FRAME_HEIGHT))
    print(f'SCREEN_WIDTH = {SCREEN_WIDTH}, SCREEN_HEIGHT = {SCREEN_HEIGHT}')
    pass
    blank_image = np.zeros((SCREEN_HEIGHT,SCREEN_WIDTH,3), dtype=np.uint8)+255

    while True:
        # We get a new frame from the webcam
        _, frame = webcam.read()
        # We send this frame to GazeTracking to analyze it
        gaze.refresh(frame)

        frame = gaze.annotated_frame()

        current_time = time.time()
        elapsed_time = current_time - start_time

        ini = time.time()  # inicia tempo dentro do while
        res = ini - temp_ini  # diferença do tempo inicial e o tempo dentro do while
        time_results.append(res)  # colocar res no vetor

        text = ""
        if gaze.is_blinking():
            text = "Blinking"
        elif gaze.is_right():
            text = "Looking right"
            # contadorRight += 1
        elif gaze.is_left():
            text = "Looking left"
            # contadorLeft += 1
        elif gaze.is_center():
            text = "Looking center"

        #cv2.putText(frame, text, (90, 60), cv2.FONT_HERSHEY_DUPLEX, 1.6, (255, 0, 0), 2)
        gaze_x = gaze.horizontal_ratio()
        gaze_y = gaze.vertical_ratio()

        if res >= 0.1:

            contador = contador + 100
            if gaze_x != None and gaze_y != None:

                # +------X------+
                # |             |
                # X             X
                # |             |
                # +------X------+

                pos = ()
                if contador < 3000:
                    pos = "esquerda"
                elif 3000 <= contador < 6000:
                    pos = "cima"
                elif 6000 <= contador < 9000:
                    pos = "direita"
                elif 9000 <= contador < 12000:
                    pos = "baixo"
                else:
                    break

                cv2.imshow(window_name, cv2.imread(f'img_{pos}.png'))

                lista_gaze_x.append(gaze_x)
                lista_gaze_y.append(gaze_y)

                coords.append([pd.Timedelta(milliseconds=contador), gaze_x, gaze_y])

        #cv2.imshow("Demo", cv2.imread(path)) and elapsed_time > duration
        if cv2.waitKey(1) & 0xFF == ord('q'):
            print(f'Main loop finished after {str(int(elapsed_time))} seconds.')
            break

    # ordenar os x de forma crescente
    lista_gaze_x.append(lista_gaze_x.sort())
    lista_gaze_x.pop()

    # ordenar os y de forma crescente
    lista_gaze_y.append(lista_gaze_y.sort())
    lista_gaze_y.pop()

    # pegar 10 primeiros x
    menores_x = lista_gaze_x[:10]
    # pegar 10 ultimos x
    maiores_x = lista_gaze_x[-10:]

    # x_min = media dos 10 primeiros gaze_x
    x_min = sum(menores_x) / len(menores_x)
    print('X min ', x_min)
    # x_max = media dos 10 ultimos gaze_x
    x_max = sum(maiores_x) / len(maiores_x)
    print('X max ', x_max)

    # pegar 10 primeiros y
    menores_y = lista_gaze_y[:10]
    # pegar 10 ultimos y
    maiores_y = lista_gaze_y[-10:]

    # y_min = media dos 10 primeiros gaze_y
    y_min = sum(menores_y) / len(menores_y)
    print('Y min ', y_min)

    # y_max = media dos 10 ultimos gaze_y
    y_max = sum(maiores_y) / len(maiores_y)
    print('Y max ', y_max)

    # return [x_min, x_max, y_min, y_max]

# apenas para testes
if __name__ == '__main__':
    calibragem()


text:
import fiona
import pandas as pd
import geopandas as gpd
from tqdm import tqdm
from gerrychain import Graph
from networkx import is_connected, connected_components

# Census tracts
# See `data/src/ACS_2018_5YR_TRACT_36_NEW_YORK.gdb/TRACT_METADATA_2018.txt`
tracts_file = 'data/src/ACS_2018_5YR_TRACT_36_NEW_YORK.gdb'
keep_layers = ['X02_RACE']
geo_layer = 'ACS_2018_5YR_TRACT_36_NEW_YORK'
ej_file = 'data/src/FPTZ_NY_3/Map_data.csv'
POP_COL = 'B02001e1' # depends on what layers we load, see the metadata file above

gdb_layers = fiona.listlayers(tracts_file)
available_tract_layers = [l for l in gdb_layers if l.startswith('X')]

def load_tracts():
    tracts = gpd.read_file(tracts_file, layer=geo_layer)
    for layer in tqdm(keep_layers, desc='Merging geodatabase tract layers'):
        other = gpd.read_file(tracts_file, layer=layer)
        other['GEOID'] = other['GEOID'].str[7:]
        tracts = tracts.merge(other.drop('geometry', axis=1), on='GEOID')
    tracts['population'] = tracts[POP_COL]

    ej_df = pd.read_csv(ej_file, dtype={'Tract_ID': str})
    ej_df['Wind_Class'] = pd.factorize(ej_df['Wind_Class'], sort=True)[0] # Convert letter classes to integers
    ej_df['Solar_Class'] = pd.factorize(ej_df['Solar_Class'], sort=True)[0]
    tracts = tracts.merge(ej_df, left_on='GEOID', right_on='Tract_ID')

    return tracts

units = load_tracts()

print('Layers:', available_tract_layers)
print('Using layers:', keep_layers)
print('Tracts:', len(units))

# Load/build the node graph (node=unit)
print('Generating graph from census tracts...')
graph = Graph.from_geodataframe(units)

def find_node_by_geoid(geoid):
    for node in graph:
        if graph.nodes[node]['GEOID'] == geoid:
            return node

# Check if any disconnected components (islands), need to connect them manually
# <https://gerrychain.readthedocs.io/en/latest/user/islands.html#strategy-2-connect-the-components-manually>
print('Fixing islands (if any)...')
if graph.islands:
    components = list(connected_components(graph))
    biggest_component_size = max(len(c) for c in components)
    problem_components = [c for c in components if len(c) < biggest_component_size]
    problem_nodes = [node for component in problem_components for node in component]

    # Use this to plot out the map, to locate the problem areas (they will be yellow)
    # For NY, the problem feature is Liberty Island and Ellis Island (they're combined into a single feature)
    problem_geoids = [graph.nodes[node]['GEOID'] for node in problem_nodes]
    # is_a_problem = df['GEOID'].isin(problem_geoids)
    # df = units.to_crs({'init': 'epsg:26986'})
    # df.plot(column=is_a_problem, figsize=(10, 10))
    # plt.axis('off')
    # plt.show()

    # Manually connect them to Governor's Island
    # Look up GEOIDs here (property here is "FIPS"):
    # <https://hub.arcgis.com/datasets/esri::usa-tracts?geometry=-74.066%2C40.685%2C-74.011%2C40.696>
    to_node = find_node_by_geoid('36061000500')
    from_node = find_node_by_geoid(problem_geoids[0])
    graph.add_edge(from_node, to_node)

# Assert that all islands connected
assert is_connected(graph)

print('Saving graph...')
graph.to_json('data/gen/graph.json')

reason: from networkx


text:
#!/usr/bin/env python

# Copyright (C) 2017 Udacity Inc.
#
# This file is part of Robotic Arm: Pick and Place project for Udacity
# Robotics nano-degree program
#
# All Rights Reserved.

# Author: Harsh Pandya

# import modules
import rospy
import tf
from kuka_arm.srv import *
from trajectory_msgs.msg import JointTrajectory, JointTrajectoryPoint
from geometry_msgs.msg import Pose
from mpmath import *
from sympy import *
import numpy as np

def rot_x(q):
    R_x = Matrix([[ 1,              0,        0],
                  [ 0,         cos(q),  -sin(q)],
                  [ 0,         sin(q),  cos(q)]])
    return R_x

def rot_y(q):
    R_y = Matrix([[ cos(q),        0,  sin(q)],
                  [      0,        1,       0],
                  [-sin(q),        0, cos(q)]])
    return R_y

def rot_z(q):
    R_z = Matrix([[ cos(q),  -sin(q),       0],
                  [ sin(q),   cos(q),       0],
                  [      0,        0,       1]])
    return R_z

def rtd(r):
    return r*180/pi

def dtr(d):
    return d*pi/180

def TransMatrix(alpha, a, d, q):
    T = Matrix([[            cos(q),           -sin(q),           0,             a],
                [ sin(q)*cos(alpha), cos(q)*cos(alpha), -sin(alpha), -sin(alpha)*d],
                [ sin(q)*sin(alpha), cos(q)*sin(alpha),  cos(alpha),  cos(alpha)*d],
                [                 0,                 0,           0,             1]])
    return T

def handle_calculate_IK(req):
    rospy.loginfo("Received %s eef-poses from the plan" % len(req.poses))
    if len(req.poses) < 1:
        print "No valid poses received"
        return -1
    else:

        ### Your FK code here

        # Create symbols
        q1,q2,q3,q4,q5,q6,q7 = symbols('q1:8')
        alpha0, alpha1,alpha2,alpha3,alpha4,alpha5,alpha6 = symbols('alpha0:7')
        a0, a1,a2,a3,a4,a5,a6 = symbols('a0:7')
        d1,d2,d3,d4,d5,d6,d7 = symbols('d1:8')
        #
        #
        # Create Modified DH parameters
        DH_table = {alpha0:        0, a0:      0, d1:  0.75, q1: q1,
             alpha1: rad(-90), a1:   0.35, d2:     0, q2: q2-rad(90),
             alpha2:        0, a2:   1.25, d3:     0, q3: q3,
             alpha3: rad(-90), a3: -0.054, d4:  1.50, q4: q4,
             alpha4:  rad(90), a4:      0, d5:     0, q5: q5,
             alpha5: rad(-90), a5:      0, d6:     0, q6: q6,
             alpha6:        0, a6:      0, d7: 0.303, q7: 0
        }
        #
        #
        # Define Modified DH Transformation matrix
        T0_1 = TransMatrix(alpha0, a0, d1, q1).subs(DH_table)
        T1_2 = TransMatrix(alpha1, a1, d2, q2).subs(DH_table)
        T2_3 = TransMatrix(alpha2, a2, d3, q3).subs(DH_table)
        T3_4 = TransMatrix(alpha3, a3, d4, q4).subs(DH_table)
        T4_5 = TransMatrix(alpha4, a4, d5, q5).subs(DH_table)
        T5_6 = TransMatrix(alpha5, a5, d6, q6).subs(DH_table)
        T6_EE = TransMatrix(alpha6, a6, d7, q7).subs(DH_table)

        T0_EE = simplify(T0_1 * T1_2 * T2_3 * T3_4 * T4_5 * T5_6 * T6_EE)
        #
        #
        # Extract rotation matrices from the transformation matrices
        #
        #
            ###

        # Initialize service response
        joint_trajectory_list = []
        for x in xrange(0, len(req.poses)):
            # IK code starts here
            joint_trajectory_point = JointTrajectoryPoint()

            # Extract end-effector position and orientation from request
            # px,py,pz = end-effector position
            # roll, pitch, yaw = end-effector orientation

            px = req.poses[x].position.x
            py = req.poses[x].position.y
            pz = req.poses[x].position.z

            (roll, pitch, yaw) = tf.transformations.euler_from_quaternion(
                [req.poses[x].orientation.x, req.poses[x].orientation.y,
                    req.poses[x].orientation.z, req.poses[x].orientation.w])

            ### Your IK code here
            r, p, y = symbols('r p y')

            R_x = rot_x(r)
            R_y = rot_y(p)
            R_z = rot_z(y)

            R_EE = R_z * R_y * R_x
            R_EE, WC, theta1, theta2, theta3 = calculate_123(R_EE, px, py, pz, roll, pitch, yaw)

            R0_3 = T0_1[0:3,0:3] * T1_2[0:3,0:3] * T2_3[0:3,0:3]
            R0_3 = R0_3.evalf(subs={q1:theta1, q2:theta2, q3:theta3})

            R3_6 = R0_3.inv("LU") * R_EE

            theta6 = atan2(-R3_6[1,1], R3_6[1,0])# +0.45370228
            sq5 = -R3_6[1,1]/sin(theta6)
            cq5 = R3_6[1,2]
            theta5 = atan2(sq5, cq5)
            sq4 = R3_6[2,2]/sin(theta5)
            cq4 = -R3_6[0,2]/sin(theta5)
            theta4 = atan2(sq4, cq4)

            if x >= len(req.poses):
                theta5 = theta5_fin
                theta6 = theta6_fin

            theta1_fin = theta1
            theta2_fin = theta2
            theta3_fin = theta3
            theta4_fin = theta4
            theta5_fin = theta5
            theta6_fin = theta6
            # Compensate for rotation discrepancy between DH parameters and Gazebo
            #
            #
            # Calculate joint angles using Geometric IK method
            #
            #
            ###

            # Populate response for the IK request
            # In the next line replace theta1,theta2...,theta6 by your joint angle variables
	        joint_trajectory_point.positions = [theta1, theta2, theta3, theta4, theta5, theta6]
	        joint_trajectory_list.append(joint_trajectory_point)

        rospy.loginfo("length of Joint Trajectory List: %s" % len(joint_trajectory_list))
        return CalculateIKResponse(joint_trajectory_list)

def IK_server():
    # initialize node and declare calculate_ik service
    rospy.init_node('IK_server')
    s = rospy.Service('calculate_ik', CalculateIK, handle_calculate_IK)
    print "Ready to receive an IK request"
    rospy.spin()

if __name__ == "__main__":
    IK_server()

reason: import numpy,from sympy,from mpmath


text:
# Include dependencies
import numpy as np
from numpy import exp, pi
from numpy.linalg import norm
from . import signal as sig
from .signal import RECT

def simulate_phs(platform, points=[[0, 0, 0]], amplitudes=[1]):
    ##############################################################################
    #                                                                            #
    #  This file takes a list of target locations and amplitudes and saves the   #
    #  demodulated signal to './phase_history.npy'. It also outputs the signal   #
    #  to the function caller.                                                   #
    #                                                                            #
    ##############################################################################

    # Retrieve relevent parameters
    c = 3.0e8
    gamma = platform['chirprate']
    f_0 = platform['f_0']
    t = platform['t']
    pos = platform['pos']
    npulses = platform['npulses']
    nsamples = platform['nsamples']
    T_p = platform['T_p']

    # Simulate the phase history for each pulse, for each point
    phs = np.zeros([npulses, nsamples]) + 0j
    for i in range(npulses):
        print('simulating pulse %i' % (i + 1))

        R_0 = norm(pos[i])
        j = 0
        for p in points:
            R_t = norm(pos[i] - p)
            dr = R_t - R_0
            phase = pi * gamma * (2 * dr / c) ** 2 - \
                    2 * pi * (f_0 + gamma * t) * 2 * dr / c

            phs[i, :] += amplitudes[j] * exp(1j * phase) * RECT((t - 2 * dr / c), T_p)

            j += 1

    # np.save('./phase_history.npy', phs)

    return (phs)

def RVP_correct(phs, platform):
    ##############################################################################
    #                                                                            #
    #  Corrects Residual Video Phase using the formulation in Carrera Appendix C #
    #                                                                            #
    ##############################################################################

    # Retrieve relevent parameters
    c = 3.0e8
    gamma = platform['chirprate']
    nsamples = platform['nsamples']
    npulses = platform['npulses']
    dr = platform['delta_r']

    # Calculate frequency sample locations w.r.t. demodulated fast time
    f_t = np.linspace(-nsamples / 2, nsamples / 2, nsamples) * \
          2 * gamma / c * dr

    # Calculate correction factor
    S_c = exp(-1j * pi * f_t ** 2 / gamma)
    S_c2 = np.tile(S_c, [npulses, 1])

    # Filter original signal
    PHS = sig.ft(phs)
    phs_corr = sig.ift(PHS * S_c2)

    return (phs_corr)

def phs_to_const_ref(phs, platform, upchirp=1):
    ##############################################################################
    #                                                                            #
    #  This program converts a phase history that was demodulated using a pulse  #
    #  dependant range to scene center to a phase history that is demodulated    #
    #  using a fixed reference.  The fixed reference is defined as the minimum   #
    #  range to scene center.                                                    #
    #                                                                            #
    ##############################################################################

    # Retrieve relevent parameters
    c = 3.0e8
    f0 = platform['f_0']
    gamma = platform['chirprate']
    pos = platform['pos']
    t = platform['t']
    npulses = platform['npulses']

    # Define ranges to scene center and new reference range
    # using Carrera's notation
    R_a = norm(pos, axis=-1)
    R_a = np.array([R_a]).T
    R_s = R_a.min()
    DR = R_a - R_s

    # Derive fast time using constant reference
    t = np.tile(t, (npulses, 1))

    # Remove original reference and incorporate new reference into phs
    sgn = (-1) ** (upchirp)
    phs = phs * exp(sgn * 1j * 4 * pi * gamma / c * (f0 / gamma + t) * DR)

    return (phs)

def reMoComp(phs, platform, center=np.array([0, 0, 0])):
    ##############################################################################
    #                                                                            #
    #  This is the re-motion compensation algorithm.  It re-motion compensates   #
    #  the phase history to a new scene center.  The "center" argument is the    #
    #  3D vector (in meters) that points to the new scene center using the old   #
    #  scene center as the origin.                                               #
    #                                                                            #
    ##############################################################################

    # Retrieve relevent parameters
    k_r = platform['k_r']
    pos = platform['pos']

    R0 = np.array([norm(pos, axis=-1)]).T
    RC = np.array([norm(pos - center, axis=-1)]).T
    dr = R0 - RC
    remocomp = np.exp(-1j * k_r * dr)

    phs = phs * remocomp

    return (phs)

max_stars_repo_name: benwhittington/RITSAR

max_issues_repo_name: benwhittington/RITSAR

max_forks_repo_name: benwhittington/RITSAR

reason: import numpy,from numpy


text:
from __future__ import division, print_function, absolute_import

from math import sqrt

from numpy.testing import *

from scipy.optimize import zeros as cc

# Import testing parameters
from scipy.optimize._tstutils import methods, mstrings, functions, \
     fstrings, description

class BenchZeros(TestCase):
    def bench_run(self):
        a = .5
        b = sqrt(3)
        repeat = 2000

        print(description)

        print('TESTING SPEED\n')
        print('times in seconds for %d iterations \n'%repeat)
        for i in range(len(functions)) :
            print('function %s\n'%fstrings[i])
            func = functions[i]
            for j in range(len(methods)) :
                meth = methods[j]
                try:
                    t = measure("meth(func,a,b)",repeat)
                except:
                    print('%s : failed'%mstrings[j])
                else:
                    print('%s : %5.3f'%(mstrings[j],t))
            print('\n\n')

if __name__ == '__main__' :
    run_module_suite()

max_stars_repo_name: globusgenomics/galaxy

max_issues_repo_name: globusgenomics/galaxy

max_forks_repo_name: globusgenomics/galaxy

reason: from numpy,from scipy


text:
from collections import defaultdict
from typing import Any, Dict, List, Literal, Tuple

import numpy as np

from scripts.utils import ENTITIES, RELATIONS

class FeatureBuilder:
    def __call__(self, raw_features):
        raise NotImplementedError()

class VotingFeatures(FeatureBuilder):
    def __init__(self, voters, weighter=None):
        self._voters = voters
        self._weighter = weighter or (lambda x, y: 1)

    def __call__(self, votes_for_label):
        votes, label = votes_for_label
        n_votes = len(self._voters)
        features = np.zeros(n_votes)
        select_votes = zip(
            *[
                (i, self._weighter(submit, label))
                for i, submit in enumerate(self._voters)
                if submit in votes
            ]
        )

        try:
            voted, weights = select_votes
        except ValueError:
            voted, weights = [], []

        voted = list(voted)
        weights = list(weights)
        features[voted] = weights
        return features

class LabelFeatures(FeatureBuilder):
    def __init__(self, labels):
        self._label2index = {label: i for i, label in enumerate(labels)}

    def __call__(self, label):
        features = np.zeros(len(self._label2index))
        index = self._label2index[label]
        features[index] = 1
        return features

class WithHandler(FeatureBuilder):
    def __init__(self, builder: FeatureBuilder, handler):
        self._builder = builder
        self._handler = handler

    def __call__(self, raw_features):
        return self._builder(self._handler(raw_features))

class ConcatenatedFeatures(FeatureBuilder):
    def __init__(self, *builders_and_handlers):
        self._builders = [
            WithHandler(builder=b, handler=h) for b, h in builders_and_handlers
        ]

    def __call__(self, raw_features):
        return np.concatenate([builder(raw_features) for builder in self._builders])

class ModelHandler:
    def __call__(
        self, annotation_votes, selected_sids=None
    ) -> Dict[str, Tuple[Any, List[int], List[Any], List[str], np.ndarray]]:
        pass

class PerCategoryModel(ModelHandler):
    def __init__(
        self,
        *,
        voters,
        labels_per_category: Dict[str, list],
        model_init,
        weighter=None,
    ):
        self._builders = {
            category: self._get_builder_according_to_labels(labels, voters, weighter)
            for category, labels in labels_per_category.items()
        }
        self._models = {category: model_init() for category in labels_per_category}
        self._label2category = {
            label: category
            for category, labels in labels_per_category.items()
            for label in labels
        }
        assert sum(len(x) for x in labels_per_category.values()) == len(
            self._label2category
        )

    @classmethod
    def _get_builder_according_to_labels(cls, labels, voters, weighter):
        if len(labels) > 1:
            return ConcatenatedFeatures(
                (LabelFeatures(labels), cls._get_label),
                (VotingFeatures(voters, weighter), cls._get_votes),
            )
        else:
            return WithHandler(VotingFeatures(voters, weighter), cls._get_votes)

    @classmethod
    def _get_label(cls, item):
        label, _ = item
        return label

    @classmethod
    def _get_votes(cls, item):
        label, votes = item
        return votes, label

    def __call__(self, annotation_votes, selected_sids=None):
        per_category = defaultdict(lambda: ([], [], [], []))
        for sid, ann, votes_per_label in annotation_votes:
            if selected_sids is None or sid in selected_sids:
                # TODO: en el caso no binario estoy no hace lo esperado
                for label, votes in votes_per_label.items():
                    if label not in self._label2category:
                        print(f"Ignoring {ann} with label {label}.")
                        continue

                    category = self._label2category[label]
                    builder = self._builders[category]
                    sids, anns, labels, features = per_category[category]

                    features.append(builder((label, votes)))
                    sids.append(sid)
                    anns.append(ann)
                    labels.append(label)

        return {
            category: (self._models[category], sids, anns, labels, np.asarray(features))
            for category, (sids, anns, labels, features) in per_category.items()
        }

class AllInOneModel(PerCategoryModel):
    def __init__(self, *, voters, labels, model_init, weighter=None):
        super().__init__(
            voters=voters,
            labels_per_category={"all": labels},
            model_init=model_init,
            weighter=weighter,
        )

class PerLabelModel(PerCategoryModel):
    def __init__(self, *, voters, labels, model_init, weighter=None):
        super().__init__(
            voters=voters,
            labels_per_category={l: [l] for l in labels},
            model_init=model_init,
            weighter=weighter,
        )

def model_handler_assistant(
    *,
    voters,
    model_init,
    labels=ENTITIES + RELATIONS,
    mode: Literal["category", "all", "each"],
    weighting_table: Dict[Tuple[str, str], float] = None,
):

    weighter = (
        (lambda submit, label: weighting_table[submit, label])
        if weighting_table is not None
        else None
    )

    if mode == "category":
        labels_per_category = defaultdict(list)
        for label in labels:
            if label in ENTITIES:
                labels_per_category["ENTITIES"].append(label)
            elif label in RELATIONS:
                labels_per_category["RELATIONS"].append(label)
            else:
                raise Exception("Unknown category!")
        return lambda: PerCategoryModel(
            voters=voters,
            labels_per_category=labels_per_category,
            model_init=model_init,
            weighter=weighter,
        )
    elif mode == "all":
        return lambda: AllInOneModel(
            voters=voters, labels=labels, model_init=model_init, weighter=weighter
        )
    elif mode == "each":
        return lambda: PerLabelModel(
            voters=voters, labels=labels, model_init=model_init, weighter=weighter
        )
    else:
        raise ValueError("Unknown mode: {mode}")


text:
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",
                      size=14)    # 解决windows环境下画图汉字乱码问题

def LogisticRegression():
    data = loadtxtAndcsv_data("data2.txt", ",", np.float64)
    X = data[:, 0:-1]
    y = data[:, -1]

    plot_data(X, y)  # 作图

    X = mapFeature(X[:, 0], X[:, 1])  # 映射为多项式
    initial_theta = np.zeros((X.shape[1], 1))  # 初始化theta
    initial_lambda = 0.1  # 初始化正则化系数，一般取0.01,0.1,1.....

    # 计算一下给定初始化的theta和lambda求出的代价J
    J = costFunction(initial_theta, X, y, initial_lambda)

   # print J  #输出一下计算的值，应该为0.693147
    # result = optimize.fmin(costFunction, initial_theta, args=(X,y,initial_lambda))    #直接使用最小化的方法，效果不好
    '''调用scipy中的优化算法fmin_bfgs（拟牛顿法Broyden-Fletcher-Goldfarb-Shanno）
    - costFunction是自己实现的一个求代价的函数，
    - initial_theta表示初始化的值,
    - fprime指定costFunction的梯度
    - args是其余测参数，以元组的形式传入，最后会将最小化costFunction的theta返回
    '''
    result = optimize.fmin_bfgs(
        costFunction, initial_theta, fprime=gradient, args=(X, y, initial_lambda))
    p = predict(X, result)  # 预测
   # print u'在训练集上的准确度为%f%%'%np.mean(np.float64(p==y)*100)   # 与真实值比较，p==y返回True，转化为float

    X = data[:, 0:-1]
    y = data[:, -1]
    plotDecisionBoundary(result, X, y)  # 画决策边界

# 加载txt和csv文件
def loadtxtAndcsv_data(fileName, split, dataType):
    return np.loadtxt(fileName, delimiter=split, dtype=dataType)

# 加载npy文件

def loadnpy_data(fileName):
    return np.load(fileName)

# 显示二维图形

def plot_data(X, y):
    pos = np.where(y == 1)  # 找到y==1的坐标位置
    neg = np.where(y == 0)  # 找到y==0的坐标位置
    # 作图
    plt.figure(figsize=(15, 12))
    plt.plot(X[pos, 0], X[pos, 1], 'ro')        # red o
    plt.plot(X[neg, 0], X[neg, 1], 'bo')        # blue o
    plt.title(u"两个类别散点图", fontproperties=font)
    plt.show()

# 映射为多项式

def mapFeature(X1, X2):
    degree = 2                     # 映射的最高次方
    out = np.ones((X1.shape[0], 1))  # 映射后的结果数组（取代X）
    '''
    这里以degree=2为例，映射为1,x1,x2,x1^2,x1x2,x2^2
    '''
    for i in np.arange(1, degree+1):
        # i = 1, 2
        for j in range(i+1):
            # 矩阵直接乘相当于matlab中的点乘.*
            # 幂次的和为i
            temp = X1**(i-j)*(X2**j)
            out = np.hstack((out, temp.reshape(-1, 1)))
    return out

# 代价函数

def costFunction(initial_theta, X, y, inital_lambda):
    m = len(y)
    J = 0

    h = sigmoid(np.dot(X, initial_theta))    # 计算h(z)
    theta1 = initial_theta.copy()           # 因为正则化j=1从1开始，不包含0，所以复制一份，前theta(0)值为0
    theta1[0] = 0

    temp = np.dot(np.transpose(theta1), theta1)
    J = (-np.dot(np.transpose(y), np.log(h))-np.dot(np.transpose(1-y),
                                                    np.log(1-h))+temp*inital_lambda/2)/m   # 正则化的代价方程
    return J

# 计算梯度

def gradient(initial_theta, X, y, inital_lambda):
    m = len(y)
    grad = np.zeros((initial_theta.shape[0]))

    h = sigmoid(np.dot(X, initial_theta))  # 计算h(z)
    theta1 = initial_theta.copy()
    theta1[0] = 0

    grad = np.dot(np.transpose(X), h-y)/m+inital_lambda/m*theta1  # 正则化的梯度
    return grad

# S型函数
def sigmoid(z):
    h = np.zeros((len(z), 1))    # 初始化，与z的长度一置
    h = 1.0/(1.0+np.exp(-z))
    return h

# 画决策边界
def plotDecisionBoundary(theta, X, y):
    pos = np.where(y == 1)  # 找到y==1的坐标位置
    neg = np.where(y == 0)  # 找到y==0的坐标位置
    # 作图
    plt.figure(figsize=(15, 12))
    plt.plot(X[pos, 0], X[pos, 1], 'ro')        # red o
    plt.plot(X[neg, 0], X[neg, 1], 'bo')        # blue o
    plt.title(u"决策边界", fontproperties=font)

    #u = np.linspace(30,100,100)
    #v = np.linspace(30,100,100)

    u = np.linspace(-1, 1.5, 50)  # 根据具体的数据，这里需要调整
    v = np.linspace(-1, 1.5, 50)

    z = np.zeros((len(u), len(v)))
    for i in range(len(u)):
        for j in range(len(v)):
            z[i, j] = np.dot(mapFeature(u[i].reshape(1, -1),
                                        v[j].reshape(1, -1)), theta)    # 计算对应的值，需要map

    z = np.transpose(z)
    # 画等高线，范围在[0,0.01]，即近似为决策边界
    plt.contour(u, v, z, [0, 0.01], linewidth=2.0)
    # plt.legend()
    plt.show()

# 预测

def predict(X, theta):
    m = X.shape[0]
    p = np.zeros((m, 1))
    p = sigmoid(np.dot(X, theta))    # 预测的结果，是个概率值

    for i in range(m):
        if p[i] > 0.5:  # 概率大于0.5预测为1，否则预测为0
            p[i] = 1
        else:
            p[i] = 0
    return p

# 测试逻辑回归函数
def testLogisticRegression():
    LogisticRegression()

if __name__ == "__main__":
    testLogisticRegression()

reason: import numpy,from scipy