Skip to content

Instantly share code, notes, and snippets.

@ellielinc
Last active November 29, 2018 23:01
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ellielinc/fce4bb2e14d3653ede5422b50666f441 to your computer and use it in GitHub Desktop.
Save ellielinc/fce4bb2e14d3653ede5422b50666f441 to your computer and use it in GitHub Desktop.
Uses multiple gridsearches to minimize the chi square of the best-fit diameter and limb darkening coefficient of a star with the linear visibililty function and generates errorbars for the best-fit parameters by calculating the probability density of those parameters.
function errorbar_LDlin_prob(data::OIdata,
dataname="", diam_min=1.0, diam_max=9.0, LD_min=0.01, LD_max=0.9, iter=200)
chisq=param->norm((abs2.(visibility_ldlin(param,data.v2_baseline))-data.v2)./data.v2_err).^2;
###taken from bic.jl in numericalmethods
function logsumexp(X)
# Compute log(sum(exp(X))) while avoiding numerical over/underflow.
return log(sum(exp.(X-maximum(X))))+maximum(X)
end
#initial grid search
LD = linspace(LD_min,LD_max,iter);
diam = linspace(diam_min, diam_max, iter);
chi2val1 = zeros(length(diam),length(LD));
for i=1:length(diam)
for j=1:length(LD)
chi2val1[i,j]=chisq([diam[i],LD[j]]);
end
end
#chi2 map
indx = find(chi2val1.>(minimum(chi2val1)+1000));
chi2val_filtered = copy(chi2val1);
chi2val_filtered[indx] = 0;
minchisq1 = find(chi2val_filtered.>0);
chi2_1 = minimum(chi2val_filtered[minchisq1]);
fig1 = imshow(rotl90(chi2val_filtered./data.nv2),ColorMap("magma"), extent=[diam[1], diam[end], 0, 1])
colorbar(); title("Initial Chi2 Map $dataname");
println("Best Chi2 Value of initial $dataname gridsearch = $chi2_1")
#refine params
best1=ind2sub(chi2val1,indmin(chi2val1));
diam_min=diam[best1[1]]-1;
diam_plus=diam[best1[1]]+1;
LD_min=LD[best1[2]]-0.1
LD_plus=LD[best1[2]]+0.1
#2nd gridsearch
LD = linspace(LD_min, LD_plus, iter);
diam = linspace(diam_min, diam_plus, iter);
chi2val2=zeros(length(diam),length(LD));
for i=1:length(diam)
for j=1:length(LD)
chi2val2[i,j]=chisq([diam[i],LD[j]]);
end
end
#refined chi2 map
indx = find(chi2val2.>(minimum(chi2val2)+500));
chi2val_filtered = copy(chi2val2);
chi2val_filtered[indx] = 0;
minchisq2 = find(chi2val_filtered.>0)
chi2_2 = minimum(chi2val_filtered[minchisq2])
figure()
fig2 = imshow(rotl90(chi2val_filtered./data.nv2),ColorMap("magma"), extent=[diam[1], diam[end], 0, 1])
colorbar(); title("Refined Chi2 Map $dataname");
println("Best Chi2 Value of refined $dataname gridsearch = $chi2_2")
#best diam param
best = ind2sub(chi2val2,indmin(chi2val2));
d_slice = diam[find(mapslices(logsumexp, -chi2val_filtered./2, dims=[2]).>0)]
diam_best = diam[best[1]]
diam_best_min = diam_best-minimum(d_slice)
diam_best_max = maximum(d_slice)-diam_best
LD_slice = LD[find(mapslices(logsumexp, -chi2val_filtered./2, dims=[1]).>0)]
LD_best = LD[best[2]]
LD_best_min = LD_best-minimum(LD_slice)
LD_best_max = maximum(LD_slice)-LD_best
bestchi2 = minimum(chi2val2)./data.nv2
println("Best Chi2: $bestchi2")
println("Best Diameter: $diam_best, +$diam_best_max, -$diam_best_min")
println("Best LDC: $LD_best, +$LD_best_max, -$LD_best_min")
return bestchi2, diam_best, diam_best_max, diam_best_min, LD_best, LD_best_max, LD_best_min
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment