Updated (9abef609) · Commits · Giorgio Calderone / CL_1ES_1927p654

run.jl

+19 −18

Original line number	Diff line number	Diff line
		@@ -10,7 +10,8 @@ using CL_1ES_1927p654
		using Dates
		using TextParse

		mkdir("output")
		path = "output-" * readlines(`git branch --show-current`)[1]
		mkdir("$(path)")
		epoch_filenames = Vector{String}()
		for (root, dirs, files) in walkdir("AT2018zf")
		for file in files
		@@ -48,15 +49,15 @@ spec = read_spec("B03", resolution=200.)
		add_spec!(source, spec);
		res = fit(source);
		viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
		filename="output/results_boller.html")
		filename="$(path)/results_boller.html")
		println(res.bestfit[:OIII_5007].norm.val )
		# 0.00029085925218136
		# 0.00029150559312419105
		println(res.bestfit[:OIII_5007].norm.val / res.bestfit[:galaxy].norm.val )
		# 59.406476855719376
		# 63.87625309927029



		if !isfile("output/scale_fwhm.dat")
		if !isfile("$(path)/scale_fwhm.dat")
		all_scale = fill(1., length(epoch_filenames))
		all_fwhm = fill(NaN,length(epoch_filenames))
		for id in 1:length(epoch_filenames)
		@@ -71,7 +72,7 @@ if !isfile("output/scale_fwhm.dat")
		add_spec!(source, spec);
		res = fit(source);
		viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
		filename="output/results_$(id).html")
		filename="$(path)/results_$(id).html")

		# Update all_scale and store FWHM
		@info res.bestfit[:OIII_5007].norm.val
		@@ -81,15 +82,15 @@ if !isfile("output/scale_fwhm.dat")
		end
		end
		all_scale[.!isfinite.(all_fwhm)] .= NaN
		serialize("output/scale_fwhm.dat", (all_scale, all_fwhm))
		serialize("$(path)/scale_fwhm.dat", (all_scale, all_fwhm))
		else
		(all_scale, all_fwhm) = deserialize("output/scale_fwhm.dat")
		(all_scale, all_fwhm) = deserialize("$(path)/scale_fwhm.dat")
		end

		@gp :prenorm "set key bottom right" "set grid" xlab="Epoch" ylab="Value" :-
		@gp :- :prenorm all_scale./1e-18 "w lp t '[OIII] norm.'" :-
		@gp :- :prenorm all_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)"
		save(:prenorm, term="png size 800,600", output="output/oiii_norm_fwhm.png")
		save(:prenorm, term="png size 800,600", output="$(path)/oiii_norm_fwhm.png")

		w = 10 .^(2:0.01:log10(400));
		r = sqrt.(400^2 .- w.^2) ./ 2.355;
		@@ -113,9 +114,9 @@ for id in 1:length(epoch_filenames)
		@gp :- :allepochs x./(1+Z) y "u 1:(\$2+$id):($id) w l notit lc pal" :-
		end
		@gp :- :allepochs yr=[1, 29]
		save(:allepochs, term="png size 800,600", output="output/all_epochs.png")
		save(:allepochs, term="png size 800,600", output="$(path)/all_epochs.png")
		@gp :- :allepochs xr=[4900, 5100]
		save(:allepochs, term="png size 800,2000", output="output/all_epochs_zoom.png")
		save(:allepochs, term="png size 800,2000", output="$(path)/all_epochs_zoom.png")



		@@ -136,9 +137,9 @@ resolution = Dict(
		job = :all
		chosen_epochs = dict_chosen_epochs[job]
		Nloop = 6
		(all_scale, all_fwhm) = deserialize("output/scale_fwhm.dat")
		(all_scale, all_fwhm) = deserialize("$(path)/scale_fwhm.dat")
		for loop in 1:Nloop
		file = "output/results_$(job)_$(loop).dat"
		file = "$(path)/results_$(job)_$(loop).dat"
		if !isfile(file)
		source = QSO{q1927p654}("1ES 1927+654", Z, ebv=EBV);
		source.options[:min_spectral_coverage][:OIII_5007] = 0.5
		@@ -150,9 +151,9 @@ for loop in 1:Nloop
		end
		res = multi_fit(source);
		viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
		filename="output/results_$(job)_$(loop).html")
		filename="$(path)/results_$(job)_$(loop).html")
		viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
		filename="output/results_$(job)_$(loop)_rebin4.html", rebin=4)
		filename="$(path)/results_$(job)_$(loop)_rebin4.html", rebin=4)

		# Find best normalization for [OIII]
		multi2 = deepcopy(res.multi);
		@@ -183,7 +184,7 @@ end

		OIII_norm_evol = fill(NaN, length(chosen_epochs))
		for loop in 1:Nloop
		file = "output/results_$(job)_$(loop).dat"
		file = "$(path)/results_$(job)_$(loop).dat"
		(res, all_scale, OIII_norm) = deserialize(file);
		OIII_norm_evol = hcat(OIII_norm_evol, OIII_norm)
		println(res.bestfit.cost / res.bestfit.dof)
		@@ -309,7 +310,7 @@ for i in 1:nrow(tab)
		ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
		@gp :- domain(res.model[i])[:] res.model[i](:qso_cont) ss
		end
		save(:prenorm, term="png size 800,600", output="output/evolution.png")
		save(:prenorm, term="png size 800,600", output="$(path)/evolution.png")


		@gp "set grid"
		@@ -321,4 +322,4 @@ for i in 1:length(chosen_epochs)
		ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
		@gp :- x y .- Spline1D(x, y0)(5007.) ss
		end
		save(:prenorm, term="png size 800,600", output="output/evolution_oiii_norm.png")
		save(:prenorm, term="png size 800,600", output="$(path)/evolution_oiii_norm.png")

src/Recipe_multi.jl

+24 −24

Original line number	Diff line number	Diff line
		@@ -6,7 +6,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		mzer.config.ftol = mzer.config.gtol = mzer.config.xtol = 1.e-6

		# Initialize components and guess initial values
		println(source.log, "\nFit continuum components...")
		println(logio(source), "\nFit continuum components...")

		multi = MultiModel()
		for id in 1:Nspec
		@@ -60,7 +60,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		@patch! model m -> m[:balmer].norm *= m[:qso_cont].norm
		end
		end
		bestfit = fit!(multi, source.data, minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(multi, source.data, minimizer=mzer); show(logio(source), bestfit)

		# QSO continuum renormalization
		for id in 1:Nspec
		@@ -69,7 +69,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		c = model[:qso_cont]
		initialnorm = c.norm.val
		if c.norm.val > 0
		println(source.log, "$id: Cont. norm. (before): ", c.norm.val)
		println(logio(source), "$id: Cont. norm. (before): ", c.norm.val)
		while true
		residuals = (model() - source.data[id].val) ./ source.data[id].unc
		(count(residuals .< 0) / length(residuals) > 0.9) && break
		@@ -77,9 +77,9 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		c.norm.val *= 0.99
		evaluate!(model)
		end
		println(source.log, "$id : Cont. norm. (after) : ", c.norm.val)
		println(logio(source), "$id : Cont. norm. (after) : ", c.norm.val)
		else
		println(source.log, "$id: Skipping cont. renormalization")
		println(logio(source), "$id: Skipping cont. renormalization")
		end
		freeze(model, :qso_cont)
		(:galaxy in keys(model)) && freeze(model, :galaxy)
		@@ -88,7 +88,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		evaluate!(multi)

		# Fit iron templates
		println(source.log, "\nFit iron templates...")
		println(logio(source), "\nFit iron templates...")
		for id in 1:Nspec
		model = multi[id]
		λ = source.domain[id][:]
		@@ -105,7 +105,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		model[:ironuv].norm.val = 0.5
		push!(iron_components, :ironuv)
		else
		println(source.log, "Ignoring ironuv component on prediction $id (threshold: $threshold)")
		println(logio(source), "Ignoring ironuv component on prediction $id (threshold: $threshold)")
		end
		end

		@@ -125,14 +125,14 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		freeze(model, :ironoptna) # will be freed during last run
		push!(iron_components, :ironoptbr, :ironoptna)
		else
		println(source.log, "Ignoring ironopt component on prediction $id (threshold: $threshold)")
		println(logio(source), "Ignoring ironopt component on prediction $id (threshold: $threshold)")
		end
		end
		if length(iron_components) > 0
		model[:Iron] = SumReducer(iron_components)
		model[:main] = SumReducer([:Continuum, :Iron])
		evaluate!(model)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
		else
		model[:Iron] = @expr m -> [0.]
		model[:main] = SumReducer([:Continuum, :Iron])
		@@ -146,7 +146,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		# Add emission lines
		line_names = [collect(keys(source.line_names[id])) for id in 1:Nspec]
		line_groups = [unique(collect(values(source.line_names[id]))) for id in 1:Nspec]
		println(source.log, "\nFit known emission lines...")
		println(logio(source), "\nFit known emission lines...")
		for id in 1:Nspec
		model = multi[id]
		λ = source.domain[id][:]
		@@ -193,7 +193,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		if isa(c, QSFit.SpecLineGauss)
		c.spec_res_kms = source.spectra[id].resolution
		else
		println(source.log, "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
		println(logio(source), "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
		end
		end
		end
		@@ -286,7 +286,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		model[:OIII_5007].norm.fixed = 1
		model[:OIII_5007].norm.val = 1.

		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)

		for lname in line_names[id]
		freeze(model, lname)
		@@ -294,7 +294,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		end

		# Add unknown lines
		println(source.log, "\nFit unknown emission lines...")
		println(logio(source), "\nFit unknown emission lines...")
		if source.options[:n_unk] > 0
		for id in 1:Nspec
		model = multi[id]
		@@ -359,7 +359,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		model[cname].center.high = λ[iadd] + λ[iadd]/10.

		thaw(model, cname)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
		freeze(model, cname)
		end
		end
		@@ -367,7 +367,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654

		# ----------------------------------------------------------------
		# Last run with all parameters free
		println(source.log, "\nLast run with all parameters free...")
		println(logio(source), "\nLast run with all parameters free...")
		for id in 1:Nspec
		model = multi[id]
		thaw(model, :qso_cont)
		@@ -402,26 +402,26 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
		if bestfit[id][cname].norm.val == 0.
		freeze(model, cname)
		needs_fitting = true
		println(source.log, "Disabling $cname (norm. = 0)")
		println(logio(source), "Disabling $cname (norm. = 0)")
		elseif bestfit[id][cname].norm.unc / bestfit[id][cname].norm.val > 3
		model[cname].norm.val = 0.
		freeze(model, cname)
		needs_fitting = true
		println(source.log, "Disabling $cname (unc. / norm. > 3)")
		println(logio(source), "Disabling $cname (unc. / norm. > 3)")
		end
		end
		end
		if needs_fitting
		println(source.log, "\nRe-run fit...")
		println(logio(source), "\nRe-run fit...")
		bestfit = fit!(multi, source.data, minimizer=mzer)
		end

		println(source.log)
		show(source.log, bestfit)
		println(logio(source))
		show(logio(source), bestfit)

		out = QSFit.QSFitMultiResults(source, multi, bestfit)
		elapsed = time() - elapsed
		println(source.log, "\nElapsed time: $elapsed s")
		QSFit.close_log(source)

		return QSFit.QSFitMultiResults(source, multi, bestfit)
		println(logio(source), "\nElapsed time: $elapsed s")
		QSFit.close_logio(source)
		return out
		end

src/Recipe_single.jl

+24 −24

Original line number	Diff line number	Diff line
		@@ -4,7 +4,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		mzer.config.ftol = mzer.config.gtol = mzer.config.xtol = 1.e-6

		# Initialize components and guess initial values
		println(source.log, "\nFit continuum components...")
		println(logio(source), "\nFit continuum components...")
		λ = source.domain[id][:]
		model = Model(source.domain[id],
		:qso_cont => QSFit.qso_cont_component(TRecipe),
		@@ -43,13 +43,13 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		@patch! model m -> m[:balmer].norm *= m[:qso_cont].norm
		end

		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)

		# QSO continuum renormalization
		freeze(model, :qso_cont)
		c = model[:qso_cont]
		if c.norm.val > 0
		println(source.log, "Cont. norm. (before): ", c.norm.val)
		println(logio(source), "Cont. norm. (before): ", c.norm.val)
		initialnorm = c.norm.val
		while true
		residuals = (model() - source.data[id].val) ./ source.data[id].unc
		@@ -58,9 +58,9 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		c.norm.val *= 0.99
		evaluate!(model)
		end
		println(source.log, "Cont. norm. (after) : ", c.norm.val)
		println(logio(source), "Cont. norm. (after) : ", c.norm.val)
		else
		println(source.log, "Skipping cont. renormalization")
		println(logio(source), "Skipping cont. renormalization")
		end
		freeze(model, :qso_cont)
		(:galaxy in keys(model)) && freeze(model, :galaxy)
		@@ -68,7 +68,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		evaluate!(model)

		# Fit iron templates
		println(source.log, "\nFit iron templates...")
		println(logio(source), "\nFit iron templates...")
		iron_components = Vector{Symbol}()
		if source.options[:use_ironuv]
		fwhm = 3000.
		@@ -81,7 +81,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		model[:ironuv].norm.val = 0.5
		push!(iron_components, :ironuv)
		else
		println(source.log, "Ignoring ironuv component (threshold: $threshold)")
		println(logio(source), "Ignoring ironuv component (threshold: $threshold)")
		end
		end

		@@ -101,14 +101,14 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		freeze(model, :ironoptna) # will be freed during last run
		push!(iron_components, :ironoptbr, :ironoptna)
		else
		println(source.log, "Ignoring ironopt component (threshold: $threshold)")
		println(logio(source), "Ignoring ironopt component (threshold: $threshold)")
		end
		end
		if length(iron_components) > 0
		model[:Iron] = SumReducer(iron_components)
		model[:main] = SumReducer([:Continuum, :Iron])
		evaluate!(model)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
		else
		model[:Iron] = @expr m -> [0.]
		model[:main] = SumReducer([:Continuum, :Iron])
		@@ -121,7 +121,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		# Add emission lines
		line_names = collect(keys(source.line_names[id]))
		line_groups = unique(collect(values(source.line_names[id])))
		println(source.log, "\nFit known emission lines...")
		println(logio(source), "\nFit known emission lines...")
		resid = source.data[id].val - model() # will be used to guess line normalization
		for (cname, comp) in source.line_comps[id]
		model[cname] = comp
		@@ -165,7 +165,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		if isa(c, QSFit.SpecLineGauss)
		c.spec_res_kms = source.spectra[id].resolution
		else
		println(source.log, "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
		println(logio(source), "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
		end
		end
		end
		@@ -255,13 +255,13 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		m[:br_Hb].fwhm = m[:br_Ha].fwhm
		end

		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
		for lname in line_names
		freeze(model, lname)
		end

		# Add unknown lines
		println(source.log, "\nFit unknown emission lines...")
		println(logio(source), "\nFit unknown emission lines...")
		if source.options[:n_unk] > 0
		tmp = OrderedDict{Symbol, GFit.AbstractComponent}()
		for j in 1:source.options[:n_unk]
		@@ -312,7 +312,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		model[cname].center.high = λ[iadd] + λ[iadd]/10.

		thaw(model, cname)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
		bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
		freeze(model, cname)
		end
		end
		@@ -320,7 +320,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654

		# ----------------------------------------------------------------
		# Last run with all parameters free
		println(source.log, "\nLast run with all parameters free...")
		println(logio(source), "\nLast run with all parameters free...")
		thaw(model, :qso_cont)
		(:galaxy in keys(model)) && thaw(model, :galaxy)
		(:balmer in keys(model)) && thaw(model, :balmer)
		@@ -350,25 +350,25 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
		if bestfit[cname].norm.val == 0.
		freeze(model, cname)
		needs_fitting = true
		println(source.log, "Disabling $cname (norm. = 0)")
		println(logio(source), "Disabling $cname (norm. = 0)")
		elseif bestfit[cname].norm.unc / bestfit[cname].norm.val > 3
		model[cname].norm.val = 0.
		freeze(model, cname)
		needs_fitting = true
		println(source.log, "Disabling $cname (unc. / norm. > 3)")
		println(logio(source), "Disabling $cname (unc. / norm. > 3)")
		end
		end
		if needs_fitting
		println(source.log, "\nRe-run fit...")
		println(logio(source), "\nRe-run fit...")
		bestfit = fit!(model, source.data[id], minimizer=mzer)
		end

		println(source.log)
		show(source.log, bestfit)
		println(logio(source))
		show(logio(source), bestfit)

		out = QSFit.QSFitResults(source, model, bestfit)
		elapsed = time() - elapsed
		println(source.log, "\nElapsed time: $elapsed s")
		QSFit.close_log(source)

		return QSFit.QSFitResults(source, model, bestfit)
		println(logio(source), "\nElapsed time: $elapsed s")
		QSFit.close_logio(source)
		return out
		end