ULABDoppler/Inverse.py at main · cliulinnaeus/ULABDoppler · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
from Simulator import *
from Forward import *
from load_model import *
import numpy as np
from scipy.optimize import minimize
from sklearn.metrics import mean_squared_error

import autograd.numpy as np
from autograd import grad
from scipy.optimize import fsolve


#TODO:
# add gaussian noise
# define accuracy measure
# make doppler shift more visible
# implement maximum entropy method

def R_transpose(R, D):

    return R.T @ D

def pseudo_inverse(R, D):

    #print(np.linalg.pinv(R).shape)
    #print(D.shape)

    return np.linalg.pinv(R).T @ D

def rmse(I,computed_I):

    return mean_squared_error(I, computed_I, squared = False)

def mse(I, computed_I):

    return mean_squared_error(I, computed_I, squared = True)

def maximum_entropy(X0, D, l = 1000):

    computed_I = X0

    # phi_list = list(range(0, 10))
    # R_guess_lst = []

    # for i in phi_list:

    #     I = s_R.rotate(np.pi * 2 / len(phi_list))
    #     stellar_disk_R = s_R.get_stellar_disk(I)
    #     max_wavelength = 5
    #     R_guess = get_R(s_R, 400, max_wavelength=max_wavelength)

    #     R_guess_lst.append(R_guess)


    # R = np.hstack(tuple(R_guess_lst))

    chi2 = mse(D, R.T @ computed_I)

    S = 0
    p = computed_I/np.sum(computed_I)
    S = -np.sum(p*np.log(p))

    Q = S - l*chi2

    return -Q


if __name__ == '__main__':

    s_R = Star(np.pi/4.2, 4.5, 3.4e6, 0.5e1, 500, guess = True)
    s_D = Star(np.pi/4.2, 4.5, 3.4e6, 0.5e1, 500)
    I_file_path = './I/I_vector.csv'
    R_file_path = './R/R_matrix.csv'
    D_file_path = './D/flux_vs_wavelength_data.csv'

    R = load_R(R_file_path)
    print(R.shape)
    D, wavelengths = load_D(D_file_path)
    I = load_I(I_file_path)

    D = D + np.random.normal(0, 1, D.shape)


    computed_stellar_disk = pseudo_inverse(R,D)

    #true_stellar_disk = s.get_stellar_disk()
    s_R.plot_on_sphere(computed_stellar_disk, savefig = True, parent_directory='./stellar_disk_computed/')
    s_D.plot_on_sphere(I, savefig = True, parent_directory = './stellar_disk_ground_truths/')

    #print(computed_stellar_disk[len(computed_stellar_disk)-2])
    #print(rmse(true_stellar_disk, computed_stellar_disk))

    s_D.plot_on_sphere(computed_stellar_disk)

    '''optimize'''
    s_mem_guess = Star(np.pi/4.2, 4.5, 3.4e6, 0.5e1, 500, guess = True)
    I0 = s_mem_guess.I
    result = minimize(maximum_entropy, x0 = I0, args = (D))
    optimized_I = result.x
    s_D.plot_on_sphere(optimized_I)


    #print(s_mem_guess.I)
    #print(optimized_I)


#MAKE D AND R HAVE DIFFERENT BUT CLOSE TEST STARS