Updated (dd3df84a) · Commits · Giorgio Calderone / CL_1ES_1927p654

src/Recipe_multi.jl

+233 −156

Original line number	Original line	Diff line number	Diff line
	@@ -7,45 +7,47 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654

	# Initialize components and guess initial values		# Initialize components and guess initial values
	println(source.log, "\nFit continuum components...")		println(source.log, "\nFit continuum components...")
	preds = Vector{Prediction}()
			multi = MultiModel()
	for id in 1:Nspec		for id in 1:Nspec
	λ = source.domain[id][:]		λ = source.domain[id][:]
	pred = Prediction(source.domain[id], :Continuum => reducer_sum([:qso_cont]),		model = Model(source.domain[id],
	:qso_cont => QSFit.qso_cont_component(TRecipe))		:qso_cont => QSFit.qso_cont_component(TRecipe),
			:Continuum => SumReducer([:qso_cont]))

	if source.options[:instr_broadening]		# TODO if source.options[:instr_broadening]
	GFit.set_instr_response!(pred, (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))		# TODO GFit.set_instr_response!(model, (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))
	end		# TODO end

	push!(preds, pred)		c = model[:qso_cont]
	c = pred[:qso_cont]
	c.x0.val = median(λ)		c.x0.val = median(λ)
	c.norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(c.x0.val)		c.norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(c.x0.val)
			push!(multi, model)
	end		end
	model = Model(preds)

	for id in 1:Nspec		for id in 1:Nspec
			model = multi[id]
	λ = source.domain[id][:]		λ = source.domain[id][:]

	# Host galaxy template		# Host galaxy template
	if source.options[:use_host_template] &&		if source.options[:use_host_template] &&
	(minimum(λ) .< 5500 .< maximum(λ))		(minimum(λ) .< 5500 .< maximum(λ))
	add!(model[id], :Continuum => reducer_sum([:qso_cont, :galaxy]),		model[:galaxy] = QSFit.hostgalaxy(source.options[:host_template])
	:galaxy => QSFit.hostgalaxy(source.options[:host_template]))		model[:Continuum] = SumReducer([:qso_cont, :galaxy])
	model[id][:galaxy].norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)		model[:galaxy].norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)
	if id != ref_id		if id != ref_id
	model[id][:galaxy].norm.fixed = true		model[:galaxy].norm.fixed = true
	@patch! model m -> m[id][:galaxy].norm = m[ref_id][:galaxy].norm		@patch! multi m -> m[id][:galaxy].norm = m[ref_id][:galaxy].norm
	end		end
	end		end

	# Balmer continuum and pseudo-continuum		# Balmer continuum and pseudo-continuum
	if source.options[:use_balmer]		if source.options[:use_balmer]
	tmp = [:qso_cont, :balmer]		tmp = [:qso_cont, :balmer]
	(:galaxy in keys(model[id])) && push!(tmp, :galaxy)		(:galaxy in keys(model)) && push!(tmp, :galaxy)
	add!(model[id], :Continuum => reducer_sum(tmp),		model[:balmer] = QSFit.balmercont(0.1, 0.5)
	:balmer => QSFit.balmercont(0.1, 0.5))		model[:Continuum] = SumReducer(tmp)
	c = model[id][:balmer]		c = model[:balmer]
	c.norm.val = 0.1		c.norm.val = 0.1
	c.norm.fixed = false		c.norm.fixed = false
	c.norm.high = 0.5		c.norm.high = 0.5
	@@ -53,38 +55,40 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	c.ratio.fixed = false		c.ratio.fixed = false
	c.ratio.low = 0.1		c.ratio.low = 0.1
	c.ratio.high = 1		c.ratio.high = 1
	@patch! model[id] m -> m[:balmer].norm *= m[:qso_cont].norm		@patch! model m -> m[:balmer].norm *= m[:qso_cont].norm
	end		end
	end		end
	bestfit = fit!(model, source.data, minimizer=mzer); show(source.log, bestfit)		bestfit = fit!(multi, source.data, minimizer=mzer); show(source.log, bestfit)

	# QSO continuum renormalization		# QSO continuum renormalization
	for id in 1:Nspec		for id in 1:Nspec
	freeze(model[id], :qso_cont)		model = multi[id]
	c = model[id][:qso_cont]		freeze(model, :qso_cont)
			c = model[:qso_cont]
	initialnorm = c.norm.val		initialnorm = c.norm.val
	if c.norm.val > 0		if c.norm.val > 0
	println(source.log, "$id: Cont. norm. (before): ", c.norm.val)		println(source.log, "$id: Cont. norm. (before): ", c.norm.val)
	while true		while true
	residuals = (model[id]() - source.data[id].val) ./ source.data[id].unc		residuals = (model() - source.data[id].val) ./ source.data[id].unc
	(count(residuals .< 0) / length(residuals) > 0.9) && break		(count(residuals .< 0) / length(residuals) > 0.9) && break
	(c.norm.val < initialnorm / 5) && break # give up		(c.norm.val < initialnorm / 5) && break # give up
	c.norm.val *= 0.99		c.norm.val *= 0.99
	evaluate!(model)		evaluate!(multi)
	end		end
	println(source.log, "$id : Cont. norm. (after) : ", c.norm.val)		println(source.log, "$id : Cont. norm. (after) : ", c.norm.val)
	else		else
	println(source.log, "$id: Skipping cont. renormalization")		println(source.log, "$id: Skipping cont. renormalization")
	end		end
	freeze(model[id], :qso_cont)		freeze(model, :qso_cont)
	(:galaxy in keys(model[id])) && freeze(model[id], :galaxy)		(:galaxy in keys(model)) && freeze(model, :galaxy)
	(:balmer in keys(model[id])) && freeze(model[id], :balmer)		(:balmer in keys(model)) && freeze(model, :balmer)
	end		end
	evaluate!(model)		evaluate!(multi)

	# Fit iron templates		# Fit iron templates
	println(source.log, "\nFit iron templates...")		println(source.log, "\nFit iron templates...")
	for id in 1:Nspec		for id in 1:Nspec
			model = multi[id]
	λ = source.domain[id][:]		λ = source.domain[id][:]

	iron_components = Vector{Symbol}()		iron_components = Vector{Symbol}()
	@@ -95,8 +99,8 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	(_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)		(_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)
	threshold = get(source.options[:min_spectral_coverage], :ironuv, source.options[:min_spectral_coverage][:default])		threshold = get(source.options[:min_spectral_coverage], :ironuv, source.options[:min_spectral_coverage][:default])
	if coverage >= threshold		if coverage >= threshold
	add!(model[id], :ironuv => comp)		model[:ironuv] = comp
	model[id][:ironuv].norm.val = 0.5		model[:ironuv].norm.val = 0.5
	push!(iron_components, :ironuv)		push!(iron_components, :ironuv)
	else		else
	println(source.log, "Ignoring ironuv component on prediction $id (threshold: $threshold)")		println(source.log, "Ignoring ironuv component on prediction $id (threshold: $threshold)")
	@@ -112,67 +116,68 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	if coverage >= threshold		if coverage >= threshold
	fwhm = 500.		fwhm = 500.
	source.options[:instr_broadening] \|\| (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))		source.options[:instr_broadening] \|\| (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
	add!(model[id],		model[:ironoptbr] = comp
	:ironoptbr => comp,		model[:ironoptna] = QSFit.ironopt_narrow(fwhm)
	:ironoptna => QSFit.ironopt_narrow(fwhm))		model[:ironoptbr].norm.val = 0.5
	model[id][:ironoptbr].norm.val = 0.5		model[:ironoptna].norm.val = 0.0
	model[id][:ironoptna].norm.val = 0.0		freeze(model, :ironoptna) # will be freed during last run
	freeze(model[id], :ironoptna) # will be freed during last run
	push!(iron_components, :ironoptbr, :ironoptna)		push!(iron_components, :ironoptbr, :ironoptna)
	else		else
	println(source.log, "Ignoring ironopt component on prediction $id (threshold: $threshold)")		println(source.log, "Ignoring ironopt component on prediction $id (threshold: $threshold)")
	end		end
	end		end
	if length(iron_components) > 0		if length(iron_components) > 0
	add!(model[id], :Iron => reducer_sum(iron_components))		model[:Iron] = SumReducer(iron_components)
	add!(model[id], :main => reducer_sum([:Continuum, :Iron]))		model[:main] = SumReducer([:Continuum, :Iron])
	evaluate!(model)		evaluate!(multi)
	bestfit = fit!(model, only_id=id, source.data, minimizer=mzer); show(source.log, bestfit)		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
	else		else
	add!(model[id], :Iron => @reducer(() -> [0.]))		model[:Iron] = @expr m -> [0.]
	add!(model[id], :main => reducer_sum([:Continuum, :Iron]))		model[:main] = SumReducer([:Continuum, :Iron])
	end		end
	(:ironuv in keys(model[id])) && freeze(model[id], :ironuv)		(:ironuv in keys(model)) && freeze(model, :ironuv)
	(:ironoptbr in keys(model[id])) && freeze(model[id], :ironoptbr)		(:ironoptbr in keys(model)) && freeze(model, :ironoptbr)
	(:ironoptna in keys(model[id])) && freeze(model[id], :ironoptna)		(:ironoptna in keys(model)) && freeze(model, :ironoptna)
	end		end
	evaluate!(model)		evaluate!(multi)

	# Add emission lines		# Add emission lines
	line_names = [collect(keys(source.line_names[id])) for id in 1:Nspec]		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]		line_groups = [unique(collect(values(source.line_names[id]))) for id in 1:Nspec]
	println(source.log, "\nFit known emission lines...")		println(source.log, "\nFit known emission lines...")
	for id in 1:Nspec		for id in 1:Nspec
			model = multi[id]
	λ = source.domain[id][:]		λ = source.domain[id][:]
	resid = source.data[id].val - model[id]() # will be used to guess line normalization		resid = source.data[id].val - model() # will be used to guess line normalization
	add!(model[id], source.line_comps[id])		for (cname, comp) in source.line_comps[id]
			model[cname] = comp
			end
	for (group, lnames) in QSFit.invert_dictionary(source.line_names[id])		for (group, lnames) in QSFit.invert_dictionary(source.line_names[id])
	add!(model[id], group => reducer_sum(lnames))		model[group] = SumReducer(lnames)
	end		end
	add!(model[id], :main => reducer_sum([:Continuum, :Iron, line_groups[id]...]))		model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]...])

	if haskey(model[id], :MgII_2798)		if haskey(model, :MgII_2798)
	model[id][:MgII_2798].voff.low = -1000		model[:MgII_2798].voff.low = -1000
	model[id][:MgII_2798].voff.high = 1000		model[:MgII_2798].voff.high = 1000
	end		end
	if haskey(model[id], :OIII_5007_bw)		if haskey(model, :OIII_5007_bw)
	model[id][:OIII_5007_bw].fwhm.val = 500		model[:OIII_5007_bw].fwhm.val = 500
	model[id][:OIII_5007_bw].fwhm.low = 1e2		model[:OIII_5007_bw].fwhm.low = 1e2
	model[id][:OIII_5007_bw].fwhm.high = 1e3		model[:OIII_5007_bw].fwhm.high = 1e3
	model[id][:OIII_5007_bw].voff.low = 0		model[:OIII_5007_bw].voff.low = 0
	model[id][:OIII_5007_bw].voff.high = 2e3		model[:OIII_5007_bw].voff.high = 2e3
	end		end
	for cname in line_names[id]		for cname in line_names[id]
	model[id][cname].norm_integrated = source.options[:norm_integrated]		model[cname].norm_integrated = source.options[:norm_integrated]
	end		end

	# Guess values		# Guess values
	evaluate!(model)		evaluate!(model)
	y = source.data[id].val - model[id]()
	for cname in line_names[id]		for cname in line_names[id]
	c = model[id][cname]		c = model[cname]
	resid_at_line = Spline1D(λ, resid, k=1, bc="nearest")(c.center.val)		resid_at_line = Spline1D(λ, resid, k=1, bc="nearest")(c.center.val)
	c.norm.val *= abs(resid_at_line) / maximum(model[id](cname))		c.norm.val *= abs(resid_at_line) / maximum(model(cname))

	# If instrumental broadening is not used and the line profile		# If instrumental broadening is not used and the line profile
	# is a Gaussian one take spectral resolution into account.		# is a Gaussian one take spectral resolution into account.
	@@ -192,54 +197,118 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	end		end

	# Patch parameters		# Patch parameters
	if haskey(model[id], :OIII_4959) &&		if haskey(model, :OIII_4959) &&
	haskey(model[id], :OIII_5007)		haskey(model, :OIII_5007)
	model[id][:OIII_4959].voff.fixed = true		model[:OIII_4959].norm.fixed = true
	@patch! model[id] m -> m[:OIII_4959].voff = m[:OIII_5007].voff		model[:OIII_4959].voff.fixed = true
	end		@patch! model m -> begin
	if haskey(model[id], :NII_6549) &&		m[:OIII_4959].norm = m[:OIII_5007].norm / 3
	haskey(model[id], :NII_6583)		m[:OIII_4959].voff = m[:OIII_5007].voff
	model[id][:NII_6549].voff.fixed = true		end
	@patch! model[id] m -> m[:NII_6549].voff = m[:NII_6583].voff		end
	end		if haskey(model, :NII_6549) &&
	if haskey(model[id], :OIII_5007_bw) &&		haskey(model, :NII_6583)
	haskey(model[id], :OIII_5007)		model[:NII_6549].voff.fixed = true
	@patch! model[id] m -> begin		@patch! model m -> m[:NII_6549].voff = m[:NII_6583].voff
			end
			if haskey(model, :OIII_5007_bw) &&
			haskey(model, :OIII_5007)
			@patch! model m -> begin
	m[:OIII_5007_bw].voff += m[:OIII_5007].voff		m[:OIII_5007_bw].voff += m[:OIII_5007].voff
	m[:OIII_5007_bw].fwhm += m[:OIII_5007].fwhm		m[:OIII_5007_bw].fwhm += m[:OIII_5007].fwhm
	end		end
	end		end
			if haskey(model, :OI_6300) &&
			haskey(model, :OI_6364)
			# model[:OI_6300].norm.fixed = true
			model[:OI_6300].voff.fixed = true
			@patch! model m -> begin
			# m[:OI_6300].norm = m[:OI_6364].norm / 3
			m[:OI_6300].voff = m[:OI_6364].voff
			end
			end
			if haskey(model, :NII_6549) &&
			haskey(model, :NII_6583)
			# model[:NII_6549].norm.fixed = true
			model[:NII_6549].voff.fixed = true
			@patch! model m -> begin
			# m[:NII_6549].norm = m[:NII_6583].norm / 3
			m[:NII_6549].voff = m[:NII_6583].voff
			end
			end
			if haskey(model, :SII_6716) &&
			haskey(model, :SII_6731)
			# model[:SII_6716].norm.fixed = true
			model[:SII_6716].voff.fixed = true
			@patch! model m -> begin
			# m[:SII_6716].norm = m[:SII_6731].norm / 1.5
			m[:SII_6716].voff = m[:SII_6731].voff
			end
			end

	if haskey(model[id], :br_Hb) &&		if haskey(model, :na_Ha) &&
	haskey(model[id], :bb_Hb)		haskey(model, :na_Hb)
	# Ensure luminosity at peak of the broad base component is		model[:na_Hb].voff.fixed = true
	# smaller than the associated broad component:		@patch! model m -> m[:na_Hb].voff = m[:na_Ha].voff
	model[id][:bb_Hb].norm.high = 1		end
	model[id][:bb_Hb].norm.val = 0.5
	@patch! model[id] m -> m[:bb_Hb].norm = m[:br_Hb].norm / m[:br_Hb].fwhm m[:bb_Hb].fwhm		# The following are required to avoid degeneracy with iron
			# template
			if haskey(model, :Hg) &&
			haskey(model, :br_Hb)
			model[:Hg].voff.fixed = true
			model[:Hg].fwhm.fixed = true
			@patch! model m -> begin
			m[:Hg].voff = m[:br_Hb].voff
			m[:Hg].fwhm = m[:br_Hb].fwhm
			end
			end
			if haskey(model, :br_Hg) &&
			haskey(model, :br_Hb)
			model[:br_Hg].voff.fixed = true
			model[:br_Hg].fwhm.fixed = true
			@patch! model m -> begin
			m[:br_Hg].voff = m[:br_Hb].voff
			m[:br_Hg].fwhm = m[:br_Hb].fwhm
			end
			end
			if haskey(model, :na_Hg) &&
			haskey(model, :na_Hb)
			model[:na_Hg].voff.fixed = true
			model[:na_Hg].fwhm.fixed = true
			@patch! model m -> begin
			m[:na_Hg].voff = m[:na_Hb].voff
			m[:na_Hg].fwhm = m[:na_Hb].fwhm
			end
	end		end

	if haskey(model[id], :br_Ha) &&
	haskey(model[id], :bb_Ha)
	# Ensure luminosity at peak of the broad base component is		# Ensure luminosity at peak of the broad base component is
	# smaller than the associated broad component:		# smaller than the associated broad component:
	model[id][:bb_Ha].norm.high = 1		if haskey(model, :br_Hb) &&
	model[id][:bb_Ha].norm.val = 0.5		haskey(model, :bb_Hb)
	@patch! model[id] m -> m[:bb_Ha].norm = m[:br_Ha].norm / m[:br_Ha].fwhm m[:bb_Ha].fwhm		model[:bb_Hb].norm.high = 1
			model[:bb_Hb].norm.val = 0.5
			@patch! model m -> m[:bb_Hb].norm = m[:br_Hb].norm / m[:br_Hb].fwhm m[:bb_Hb].fwhm
			end
			if haskey(model, :br_Ha) &&
			haskey(model, :bb_Ha)
			model[:bb_Ha].norm.high = 1
			model[:bb_Ha].norm.val = 0.5
			@patch! model m -> m[:bb_Ha].norm = m[:br_Ha].norm / m[:br_Ha].fwhm m[:bb_Ha].fwhm
	end		end

	# Avoid division by zero		# Avoid division by zero
	model[id][:br_Ha].norm.low = 1.e-10		model[:br_Ha].norm.low = 1.e-10

	model[id][:na_Hg].norm.fixed = 1		model[:na_Hg].norm.fixed = 1
	model[id][:na_Hg].fwhm.fixed = 1		model[:na_Hg].fwhm.fixed = 1
	model[id][:na_Hg].voff.fixed = 1		model[:na_Hg].voff.fixed = 1
	model[id][:bb_Hg].norm.fixed = 1		model[:bb_Hg].norm.fixed = 1
	model[id][:bb_Hg].fwhm.fixed = 1		model[:bb_Hg].fwhm.fixed = 1
	model[id][:bb_Hg].voff.fixed = 1		model[:bb_Hg].voff.fixed = 1
	model[id][:br_Hg].fwhm.fixed = 1		model[:br_Hg].fwhm.fixed = 1
	model[id][:br_Hg].voff.fixed = 1		model[:br_Hg].voff.fixed = 1
	@patch! model[id] m -> begin		@patch! model m -> begin
	m[:na_Hg].norm = m[:br_Hg].norm * (m[:na_Ha].norm / m[:br_Ha].norm)		m[:na_Hg].norm = m[:br_Hg].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
	m[:na_Hg].voff = m[:na_Ha].voff		m[:na_Hg].voff = m[:na_Ha].voff
	m[:na_Hg].fwhm = m[:na_Ha].fwhm		m[:na_Hg].fwhm = m[:na_Ha].fwhm
	@@ -252,15 +321,15 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	m[:br_Hg].fwhm = m[:br_Ha].fwhm		m[:br_Hg].fwhm = m[:br_Ha].fwhm
	end		end

	model[id][:na_Hb].norm.fixed = 1		model[:na_Hb].norm.fixed = 1
	model[id][:na_Hb].fwhm.fixed = 1		model[:na_Hb].fwhm.fixed = 1
	model[id][:na_Hb].voff.fixed = 1		model[:na_Hb].voff.fixed = 1
	model[id][:bb_Hb].norm.fixed = 1		model[:bb_Hb].norm.fixed = 1
	model[id][:bb_Hb].fwhm.fixed = 1		model[:bb_Hb].fwhm.fixed = 1
	model[id][:bb_Hb].voff.fixed = 1		model[:bb_Hb].voff.fixed = 1
	model[id][:br_Hb].fwhm.fixed = 1		model[:br_Hb].fwhm.fixed = 1
	model[id][:br_Hb].voff.fixed = 1		model[:br_Hb].voff.fixed = 1
	@patch! model[id] m -> begin		@patch! model m -> begin
	m[:na_Hb].norm = m[:br_Hb].norm * (m[:na_Ha].norm / m[:br_Ha].norm)		m[:na_Hb].norm = m[:br_Hb].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
	m[:na_Hb].voff = m[:na_Ha].voff		m[:na_Hb].voff = m[:na_Ha].voff
	m[:na_Hb].fwhm = m[:na_Ha].fwhm		m[:na_Hb].fwhm = m[:na_Ha].fwhm
	@@ -273,13 +342,13 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	m[:br_Hb].fwhm = m[:br_Ha].fwhm		m[:br_Hb].fwhm = m[:br_Ha].fwhm
	end		end

	model[id][:OIII_5007].norm.fixed = 1		model[:OIII_5007].norm.fixed = 1
	model[id][:OIII_5007].norm.val = 1.		model[:OIII_5007].norm.val = 1.

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

	for lname in line_names[id]		for lname in line_names[id]
	freeze(model[id], lname)		freeze(model, lname)
	end		end
	end		end

	@@ -287,36 +356,42 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	println(source.log, "\nFit unknown emission lines...")		println(source.log, "\nFit unknown emission lines...")
	if source.options[:n_unk] > 0		if source.options[:n_unk] > 0
	for id in 1:Nspec		for id in 1:Nspec
			model = multi[id]
	tmp = OrderedDict{Symbol, GFit.AbstractComponent}()		tmp = OrderedDict{Symbol, GFit.AbstractComponent}()
	for j in 1:source.options[:n_unk]		for j in 1:source.options[:n_unk]
	tmp[Symbol(:unk, j)] = line_component(TRecipe, QSFit.UnkLine(5e3))		tmp[Symbol(:unk, j)] = line_component(TRecipe, QSFit.UnkLine(5e3))
	tmp[Symbol(:unk, j)].norm_integrated = source.options[:norm_integrated]		tmp[Symbol(:unk, j)].norm_integrated = source.options[:norm_integrated]
	end		end
	add!(model[id], :UnkLines => reducer_sum(collect(keys(tmp))), tmp)		for (cname, comp) in tmp
	add!(model[id], :main => reducer_sum([:Continuum, :Iron, line_groups[id]..., :UnkLines]))		model[cname] = comp
	evaluate!(model)		end
			model[:UnkLines] = SumReducer(collect(keys(tmp)))
			model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]..., :UnkLines])
			evaluate!(multi)
	for j in 1:source.options[:n_unk]		for j in 1:source.options[:n_unk]
	freeze(model[id], Symbol(:unk, j))		freeze(model, Symbol(:unk, j))
	end		end
	end		end
	else		else
	# Here we need a :UnkLines reducer, even when n_unk is 0		# Here we need a :UnkLines reducer, even when n_unk is 0
	for id in 1:Nspec		for id in 1:Nspec
	add!(model[id], :UnkLines => @reducer(() -> [0.]))		model = multi[id]
	add!(model[id], :main => reducer_sum([:Continuum, :Iron, line_groups[id]...]))		model[:UnkLines] = @expr m -> [0.]
			model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]...])
	end		end
	end		end
	evaluate!(model)		evaluate!(multi)

	# Set "unknown" line center wavelength where there is a maximum in		# Set "unknown" line center wavelength where there is a maximum in
	# the fit residuals, and re-run a fit.		# the fit residuals, and re-run a fit.
	for id in 1:Nspec		for id in 1:Nspec
			model = multi[id]
	λ = source.domain[id][:]		λ = source.domain[id][:]
	λunk = Vector{Float64}()		λunk = Vector{Float64}()
	while true		while true
	(length(λunk) >= source.options[:n_unk]) && break		(length(λunk) >= source.options[:n_unk]) && break
	evaluate!(model)		evaluate!(multi)
	Δ = (source.data[id].val - model[id]()) ./ source.data[id].unc		Δ = (source.data[id].val - model()) ./ source.data[id].unc

	# Avoid considering again the same region (within 1A) TODO: within resolution		# Avoid considering again the same region (within 1A) TODO: within resolution
	for l in λunk		for l in λunk
	@@ -337,57 +412,59 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	push!(λunk, λ[iadd])		push!(λunk, λ[iadd])

	cname = Symbol(:unk, length(λunk))		cname = Symbol(:unk, length(λunk))
	model[id][cname].norm.val = 1.		model[cname].norm.val = 1.
	model[id][cname].center.val = λ[iadd]		model[cname].center.val = λ[iadd]
	model[id][cname].center.low = λ[iadd] - λ[iadd]/10. # allow to shift 10%		model[cname].center.low = λ[iadd] - λ[iadd]/10. # allow to shift 10%
	model[id][cname].center.high = λ[iadd] + λ[iadd]/10.		model[cname].center.high = λ[iadd] + λ[iadd]/10.

	thaw(model[id], cname)		thaw(model, cname)
	bestfit = fit!(model, only_id=id, source.data, minimizer=mzer); show(source.log, bestfit)		bestfit = fit!(model, source.data[id], minimizer=mzer); show(source.log, bestfit)
	freeze(model[id], cname)		freeze(model, cname)
	end		end
	end		end
	evaluate!(model)		evaluate!(multi)

	# ----------------------------------------------------------------		# ----------------------------------------------------------------
	# Last run with all parameters free		# Last run with all parameters free
	println(source.log, "\nLast run with all parameters free...")		println(source.log, "\nLast run with all parameters free...")
	for id in 1:Nspec		for id in 1:Nspec
	thaw(model[id], :qso_cont)		model = multi[id]
	(:galaxy in keys(model[id])) && thaw(model[id], :galaxy)		thaw(model, :qso_cont)
	(:balmer in keys(model[id])) && thaw(model[id], :balmer)		(:galaxy in keys(model)) && thaw(model, :galaxy)
	(:ironuv in keys(model[id])) && thaw(model[id], :ironuv)		(:balmer in keys(model)) && thaw(model, :balmer)
	(:ironoptbr in keys(model[id])) && thaw(model[id], :ironoptbr)		(:ironuv in keys(model)) && thaw(model, :ironuv)
	(:ironoptna in keys(model[id])) && thaw(model[id], :ironoptna)		(:ironoptbr in keys(model)) && thaw(model, :ironoptbr)
			(:ironoptna in keys(model)) && thaw(model, :ironoptna)

	for lname in line_names[id]		for lname in line_names[id]
	thaw(model[id], lname)		thaw(model, lname)
	end		end
	for j in 1:source.options[:n_unk]		for j in 1:source.options[:n_unk]
	cname = Symbol(:unk, j)		cname = Symbol(:unk, j)
	if model[id][cname].norm.val > 0		if model[cname].norm.val > 0
	thaw(model[id], cname)		thaw(model, cname)
	else		else
	freeze(model[id], cname)		freeze(model, cname)
	end		end
	end		end
	end		end
	bestfit = fit!(model, source.data, minimizer=mzer)		bestfit = fit!(multi, source.data, minimizer=mzer)

	# Disable "unknown" lines whose normalization uncertainty is larger		# Disable "unknown" lines whose normalization uncertainty is larger
	# than 3 times the normalization		# than 3 times the normalization
	needs_fitting = false		needs_fitting = false
	for id in 1:Nspec		for id in 1:Nspec
			model = multi[id]
	for ii in 1:source.options[:n_unk]		for ii in 1:source.options[:n_unk]
	cname = Symbol(:unk, ii)		cname = Symbol(:unk, ii)
	isfixed(model[id], cname) && continue		isfixed(model, cname) && continue
	if bestfit[id][cname].norm.val == 0.		if bestfit[id][cname].norm.val == 0.
	freeze(model[id], cname)		freeze(model, cname)
	needs_fitting = true		needs_fitting = true
	println(source.log, "Disabling $cname (norm. = 0)")		println(source.log, "Disabling $cname (norm. = 0)")
	elseif bestfit[id][cname].norm.unc / bestfit[id][cname].norm.val > 3		elseif bestfit[id][cname].norm.unc / bestfit[id][cname].norm.val > 3
	model[id][cname].norm.val = 0.		model[cname].norm.val = 0.
	freeze(model[id], cname)		freeze(model, cname)
	needs_fitting = true		needs_fitting = true
	println(source.log, "Disabling $cname (unc. / norm. > 3)")		println(source.log, "Disabling $cname (unc. / norm. > 3)")
	end		end
	@@ -395,7 +472,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	end		end
	if needs_fitting		if needs_fitting
	println(source.log, "\nRe-run fit...")		println(source.log, "\nRe-run fit...")
	bestfit = fit!(model, source.data, minimizer=mzer)		bestfit = fit!(multi, source.data, minimizer=mzer)
	end		end

	println(source.log)		println(source.log)
	@@ -405,5 +482,5 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
	println(source.log, "\nElapsed time: $elapsed s")		println(source.log, "\nElapsed time: $elapsed s")
	QSFit.close_log(source)		QSFit.close_log(source)

	return reduce(source, model, bestfit)		return QSFit.QSFitMultiResults(source, multi, bestfit)
	end		end

src/Recipe_single.jl

+40 −30

File changed.

Preview size limit exceeded, changes collapsed.