__author__ = 'Jacob Bieker'
import os, sys, random
import numpy
from multiprocessing import Pool
from astropy.io import fits

# Fit plane or line iteratively
#   if more than one cluster, each cluster in separate STSDAS table
#
# Residuals are calculated according to resalgo
#      per:  perpendicular to the fitted relation
#            delta = (y - a*x1 - b*x2 - c)/sqrt(1+a^2+b^2)
#      y  :  in y     delta = (y - a*x1 - b*x2 - c)
#      x1 :  in x1
#      x2 :  in x2
#
# Minization algorithms:  (minalgo)
#    quartile : distance between upper (0.75) quartile point
#               and lower (0.25) quartile point
#    delta100 : SUM abs(res)
#    delta60  : SUM abs(res) of the central 60 percent of distribution
#    rms100   : SUM res^2
#    rms60    : SUM res^2 of the central 60 percent of distribution
#
# Zero points : median (use with delta and quartile)
#               mean   (use with rms)
#
# Bootstrap uncertainties on coefficients if nboot>0
# Seed for random is in each step   seed-(nboot-step)
#
# WARNING: This is a rather slow algorithm.
# FP for 200 galaxies in 10 clusters takes 5min on Sparc 10
# if restart=yes and nboot=0 (no uncertainties)
#
# Version: 21.06.95
#          16.10.95  Solaris roeskva
#          14.02.96  changed resalgo=y initial values, faster conv.
#          10.11.97  use rename instead of trename in 2.11
#          21.01.98  bootstrap changed to the proper way of doing it
#          06.03.98  fixed bug in bootstrap for 2 parameter fits
#                    improved random number generation
#          27.07.98  stop if itera>3*maxiter || niterb>3*maxiter
#                    no matter the other parameter's niter
#          30.03.04  handle INDEF automatically, works under IRAF 2.12.2 Redhat9.0
# Inger Jorgensen, Gemini Observatory
# e-mail: inger@gemini.edu


def random_number(number, seed, nboot):
    rand_nums = []
    if seed <= 0 or seed <= nboot:
        seed = max(seed*seed, (seed+1)*(nboot+1))
    if seed > 0:
        seed = -seed
    random.seed(a=seed)
    for i in range(number):
        rand_num = random.randint(0, 1)
        rand_nums.append(rand_num)
    return rand_nums


def min_quartile():
    return 0


def min_delta(percentage):
    if percentage == 100:
        return 0
    elif percentage == 60:
        return 0


def min_rms(percentage):
    if percentage == 100:
        return 0
    elif percentage == 60:
        return 0


def zeropoint(type_solution):
    if type_solution.lower() == "median":
        # use with delta and quartile
        return 0
    elif type_solution.lower() == "mean":
        # use with rms
        return 0


def residuals_perpendicular():
    # TODO: Convert delta = (y - a*x1 - b*x2 - c)/sqrt(1+a^2+b^2)
    return 0


def residuals_y():
    # TODO: Convert: delta = (y - a*x1 - b*x2 - c)
    return 0


def residuals_x1():
    return 0


def residuals_x2():
    return 0


def line_solve():
    # TODO: Find the Least Squares for a line
    return 0


def plane_solve():
    # TODO: Find the least Squares for a plane
    return 0


def read_clusters(*args):
    # TODO: Read in the files containing the cluster points
    return 0


def determine_type():
    # TODO: Determine the type (i.e. Plane or Line) to solve for
    return 0


def initial_guess(largest_cluster, type_solution):
    if type_solution == "line":
        return 0
    elif type_solution == "plane":
        return 0
    # TODO: Get an inital guess from the perpendicular least squares method from the largest cluster
    return


def check_guess(cluster, type_solution, guess):
    if type_solution == "line":
        return 0
    elif type_solution == "plane":
        return 0
    # TODO: check guess with another cluster


def bootstrap_cluster(cluster, nboot):
    # TODO: bootstrap a cluster to get a different distribution to check with check_guess
    return 0


def determine_uncertainty(solutions):
    # TODO: Take a list or dict of solutions and determine the uncertainty in them
    return 0

if __name__ == "__main__":
    list_clusters = [] # prompt "List of input STSDAS tables"

    filename = input("Enter the filename containing the cluster(s): ")
    tables = input("List of input STSDAS tables: ")
    min_choice = input("Distance to minimize (delta100,delta60,rms100,rms60,quartile): " or "delta100")
    res_choice = input("Residual to minimize (per,y,x1,x2): " or "per")
    y_col = input("Column name for y: " or "lre_GR_sc")
    x1_col = input("Column name for x1: " or "lsig_re")
    x2_col = input("Column name for x2 (optional): ")
    zeropoint_choice = input("Zeropoints (median, mean): " or "median")
    galaxy_name = input("Column name for galaxy: ")
    group_name = input("Column name for group: ")
    factor_change_a = float(input("Starting factor for changes in a: ") or 0.05)
    factor_change_b = float(input("Starting factor for changes in b: ") or 0.02)
    iterations = int(input("Maximum number of iterations: ") or 0)
    restart_factor = bool(input("Restart iteration with smaller factors: ") or True)
    num_bootstrap = int(input("Number of estimates for bootstrap: ") or 0)
    rand_seed = int(input("Seed for random used in bootstrap: ") or 1)

    hdulist = fits.open(filename)
    print(hdulist.info())
    print(repr(hdulist[0].header))



@@ -4290,24 +4290,25 @@
 quartile): %22
+)
  or %22delta10
@@ -4309,17 +4309,16 @@
 elta100%22
-)
 %0A    res
@@ -4367,24 +4367,25 @@
 ,y,x1,x2): %22
+)
  or %22per%22)%0A 
@@ -4381,17 +4381,16 @@
 or %22per%22
-)
 %0A    y_c
@@ -4421,16 +4421,17 @@
 for y: %22
+)
  or %22lre
@@ -4437,17 +4437,16 @@
 e_GR_sc%22
-)
 %0A    x1_
@@ -4479,16 +4479,17 @@
 or x1: %22
+)
  or %22lsi
@@ -4493,17 +4493,16 @@
 lsig_re%22
-)
 %0A    x2_
@@ -4606,16 +4606,17 @@
 mean): %22
+)
  or %22med
@@ -4619,17 +4619,16 @@
 %22median%22
-)
 %0A    gal

Move checks around