2018년 6월 27일 수요일

1주/부록: 천문학자 Aris Karastergiou 와의 인터뷰

[커세라 강좌 소개]자료기반 천문학(Data-Driven Astronomy)
https://www.coursera.org/learn/data-driven-astronomy

----------------------------------------------
Week 1: Thinking about data
제1주차: 자료의 개념
- Principles of computational thinking
  전자계산학에 기초한 자료처리의 원리
- Discovering pulsars in radio images
  전파망원경 관측영상에서 펄서 찾기
----------------------------------------------
----------------------------------------------
2강: 강좌의 구성 / Lesson 2: Course Overview
----------------------------------------------
3강: 펄사(Pulsars) / Lesson 3: Pulsars
----------------------------------------------
4강: 관측영상 겹쳐쌓기 / Lesson 4: Diving In: Image Stacking
----------------------------------------------
5강: 난관에 봉착하다 / Lesson 5: The Challenge: What went wrong ?
----------------------------------------------
6강: 개선책을 찾다 / Lesson 6: The Solution: Improving your method
----------------------------------------------
1주차 요약Module Summary
----------------------------------------------
1주/부록: 천문학자 Aris Karastergiou 와의 인터뷰

[펄사 사냥꾼]


네, 저의 이름은 아리스 카라스테규 입니다. 옥스포드 대학교의 물리학과에서 물리학을 가르치고 있죠. 그리고 저는 케이프타운 소재 웨스턴 케이프 대학과 남아프리카 공화국의 로드 대학교에서 시간제로 일을 하고 있어요. 저는 주로 펄사와 펄사 천문학에 관한 연구를 합니다. 거의 대부분 펄사관측 망원경의 기술개발이 주 업무죠. 그러니까 저는 펄사의 기전에 대해 관심이 많답니다

[당신의 연구에 대해 소개해 주세요.]


많은 사람들이 펄사를 연구에 활용하죠. 펄사의 회전 안정성을 연구에 활용하기도 하는데 펄사를 시계로 보는 겁니다. 저의 관심은 펄사의 근본을 이해하는 것인데, 세분하면 펄사에 관련된 근본적 물리적 원리를 이해하고자 하는 것이죠. 하지만 부분적으로 펄사를 시계처러 사용하는 여러 실험들의 개선을 돕고 있습니다. 그러니까 기본적으로 좀더 안정적인 시계 관측이 되도록 하는 겁니다.

[당신이 접하고 있는 대량의 자료로 인한 어려움은 무엇인가요]


그중 하나 예를 들면 펄사에서 방출된 전파가 은하를 통과할 때 지연이 생기는데 복사된 전파의 파장에 따라 지연의 정도가 약간씩 다릅니다. 다른말로하면 파장이 짧은 전파는 안테나에 좀더 일찍 도착합니다. 그에비해 긴 파장은 느리게 도착하죠. 이제 관측한 전 주파수 대역의 모든 신호를 합치기 원하면 관측한 모든 파장의 신호 말입니다. 지연을 보정해 줘야 하죠. 이 처리를 분산제거(de-dispersion)라고 합니다.

[대량의 자료 문제를 어떻게 풀고있죠?]


당신이 찾고자하는 신호에 대해 이와같은 처리를 하려면 이 신호가 어디에 있는지 알지 못하기 때문에 의미있는 신호가 어디에 있는지 매우 넓은 영역에 걸쳐 훑어야 하죠. 그리고 나서 간헐적인 펄사의 주기적 신호를 찾아낼 수 있길 바랄 겁니다.  하지만 새로 발견한 펄사의 주기를 알지 못한다는 문제가 있어요. 그래서 이번에는 펄사의 주기를 알고 싶어할 겁니다. 이번에는 펄사의 주기를 알기 위해 더 많은(상세한) 조건을 주고 관측을 해야 합니다. 그러면 더 많은 량의 자료가 수집 될텐데 초당 수 기가 바이트 대에 이릅니다. 이 자료를 처리하는 횟수는 관측당 수 백만 번의 계산을 해야하고 탐색영역의 모든 부분에 이 과정을 거쳐야 합니다.

당신이 풀어야할 문제를 생각해보세요. 한 펄사의 간헐적인 신호의 양상이 일관적인지 알고 싶은 겁니다. 자 이제, 앞서 얘기한 우리가 할 일은 간헐적 신호를 매달 관측 하는 일이죠. 그리고 그 자료를 쌓아놓는데, 수백 껀에 이르는 십여번 혹은 수십여번 반짝인 간헐 신호의 기록표가 됩니다. 그리고 나서 간헐 신호 기록 중 특정 부분에 대해 시간상 규칙성을 세워보려고 하겠죠. 따라서 시간적 규칙성은 관찰 기록의 중요한 부분에 해당합니다. 그리고 이 시간적 규칙성을 찾을 때 우리가 보통 가우시안 근사화(Gaussian Process Regression)라고 부르는 기법을 활용하죠. 그것은 원론적(Bayesian) 기법인데 수집한 자료의 본성을 도출해 낼때 사용하는 과학실험의 기본이라 할 수 있습니다. 어떤 원인요소도 개입하지 않고 가우시언 근사화라는 수학으로 도출한 펄사의 시간적 특성을 바탕으로 모형을 만들죠. 그러니까 이런 시간적으로 변화하는 간헐적 신호에 어떤 (이론의)법칙도 개입시키지 않은 겁니다. 이렇게 관측 자료에서 근사화한 목적은 시간간격을 두고 관측한 자료는 우주공간 잡음과 지연, 관측장비의 특성, 관측자의 기술 등의 영향으로 매우 성길 수 밖에 없습니다. 시간상 발생하는 일련의 과정을 부드럽게 모형화하기 위한 것입니다.

이런 작업을 위해 사용할 수 있는 도구(소프트웨어)들을 가지고 있는데 이미 개발되어 있는 라이브러리를 조합하여 우리의 용도에 맞게 자체 개발합니다. 이런 도구의 개발은 현재 널리 사용하는 파이썬 같은 언어와 부속 라이브러리의 덕을 많이 보죠. 파이썬은 현재 널리 사용되고 있기 때문에 사용자 활동이 활발하고 이에 따라 라이브러리들이 매우 풍부하게 개발되어 배포되고 있습니다. 최첨단 영역에서 최신의 개발이 이뤄지고 있다고 할 수 있죠. 그래서 여러분도 알다시피 우리는 이런 라이브러리를 활용합니다. 하지만 우리가 필요로하는 조건을 모두 충족해 주지 못하기에 가끔 한계에 직면하곤 하죠. 그래서 우리 자체의 도구를 개발하고 그에 대해 관련 공동체에 보고 합니다. 더 개선되고 미처 확인하지 못한 사항의 분석을 받고자 하는 겁니다. 이렇게 자료를 수집하는 단계에서 부터 과학적인 방법을 동원하는 것이 실질적인 어려움의 시작인 것입니다. 그리고 저의 개인적인 경험에 비춰보면 저의 박사과정 내내 이런 일을 했다고 기억합니다.

펄사의 신호가 저의 컴퓨터 화면에 끊임없이 반짝이는 것을 보고 있노라면 어떤 규칙적 양상을 띄고 있는지 알아차릴 수 있습니다. 척 보면 알 수 있어요. 그것은 아주 흥미로운 일입니다. 하지만 끝에가서 어떤 이는 자료와 매혹적인 신호를 찾아내고도 어떤 과학적인 결론에 도달하지 못하기도 합니다. 그러므로 여러분은 자료를 볼 때 그 이면에 여러분이 가지고 있는 질문이 과학에 근거하고 있는지 항상 자문해 봐야 합니다. 자료과학계(Data Science Community)는 이런 점을 (과학적 배경을)간과한 책임이 다소 있다는 생각이 들어요. 여러분도 알다시피 가장 기본은 과학적 사고입니다. 그것이 자료과학에도 적용되어야 한다고 생각 합니다. 그리고 나서 모형화 할 때 모든 물리법칙을 고려해 보는 겁니다.

2018년 6월 25일 월요일

1주/6강: 개선책을 찾다

[커세라 강좌 소개]자료기반 천문학(Data-Driven Astronomy)
https://www.coursera.org/learn/data-driven-astronomy

----------------------------------------------
Week 1: Thinking about data
제1주차: 자료의 개념
- Principles of computational thinking
  전자계산학에 기초한 자료처리의 원리
- Discovering pulsars in radio images
  전파망원경 관측영상에서 펄서 찾기
----------------------------------------------
----------------------------------------------
2강: 강좌의 구성 / Lesson 2: Course Overview
----------------------------------------------
3강: 펄사(Pulsars) / Lesson 3: Pulsars
----------------------------------------------
4강: 관측영상 겹쳐쌓기 / Lesson 4: Diving In: Image Stacking
----------------------------------------------
5강: 난관에 봉착하다Lesson 5: The Challenge: What went wrong ?
----------------------------------------------
6강: 개선책을 찾다
Lesson 6: The Solution: Improving your method (동영상 강좌) / 한글자막 / 영문자막

[강좌대본]

[00:05] 이제까지 우리가 다뤘던 내용은 통계적 특성을 활용하여 우리가 가진 장비의 감지능력 한계 아래에 놓인 펄사들을 드러나게 하는 방안에 대한 것이었다. 평균(mean)쌓기 알고리즘을 컴퓨터 프로그램으로 구현하였고 아주 잘 작동 하였다. 하지만 중간값(median)쌓기를 적용하므로서 예외적 신호에 대하여 통계적으로 좀더 탄탄한 알고리즘을 말들고자 한다.

[00:22] 중간값쌓기를 적용하려하자 막대한 저장장치의 용량을 필요로 한다는 문제에 직면 했다. 중간값 알고리즘을 곧이곧대로 적용하면 동시에 모든 관측 영상을 메모리에 올려 놓아야 한다. 그로인해 우리가 보유한 컴퓨터의 메모리 용량을 넘어선다. 이 문제를 해결하기 위한 몇가지 방법이 있다.


[00:38] 그중 메모리 용량을 늘릴 수 있지만 비용 들 것이고 어떤 해결책은 문제의 재구성(reframe)을 요한다. 또는 좀더 현명한 해법을 개발할 수도 있다.

[00:46] (첫째 방법은 더좋은 컴퓨터 구입) 아마도 가장 간단한 해법은 비용이 들더라도 더 성능좋은 컴퓨터를 마련하는 것이리라. 어느정도는 이 방법이 괜찮은 해법이 되겠으나 관측자료가 증가하면서 곧 같은 문제에 닥칠 것이다. 그럼 또 더좋은 컴퓨터를 마련해야 할 것이다. 심지어 우리 프로그램을 초고성능 슈퍼 컴퓨터에서 실행 시킬 수도 있다. 모르긴 해도 슈퍼 컴퓨터 관리자들은 그런 메모리 만을 많이차지하는 프로그램을 탐탁지 않아 할 것이다. 값비싼 중앙연산장치(CPU)를 이렇게 비효율적인 프로그램을 돌려야 하니까 말이다.

[01:11] (둘때 방법으로 영상크기 줄이기) 우리가 가진 문제를 세심하게 살펴보고 다른 접근방법은 없는지 찾아보기로 한다. 가로 200 세로 200 크기의 영상이 정말 다 필요한 것일까? 이는 과학자들이 엉뚱한 생각을 해보는 선택지 중의 하나이며 심지어 그들이 전혀 생각해 보지 못했던 것이다. 이번처럼 전혀 접해보지 못한 문제인 경우 유용하다. (과학자들의 사고논리에서 메모리 용량 부족 따위는 생각해보지 않았다. 구현은 그저 공학자, 프로그래머가 알아서 할일이라고 생각했다)


[01:28] 만일, 예를들어 영상을 50x50 화소로 자른다면 자료의 크기는 16분의 1의 비율로 줄어든다. 이는 총 소요 메모리 용량이 192기가 바이트에서 12기가 바이트로 준다는 뜻이다. 이 방법은 관측 영상에서 주목할 필요가 없는 영역을 제외시킬 때 유용하다.

[01:48] 세번째 문제접근 방식으로 알고리즘을 개선하는 것이다. 지금 문제는 중간값을 계산하기 위해 모든 자료를 메모리에 읽어 들이고 있다는 점이다. 중간값을 계산하는 중에 당장 필요하지 않은 자료를 전부 메모리에 읽어다 놓을 필요가 있을까? 이에대한 한 해법은 BinApprox 알고리즘으로 다음과 같이 작동한다.

[02:06] 매 영상을 읽어들인 즉시 각 화소의 값을 취해 값의 크기순을 정한 상자에 담는다. 이렇게 모든 영상이 처리되고 나면 영상에 있던 각 화소 값의 빈도수 도표(히스토그램)를 얻는다. 히스토그램에서 상자는 화소값에 대해 정렬되어 있으므로 히스토그램에서 오름차순으로 배열된 각 상자에 담개 화소의 빈도수를 더할 수 있다. 더한 수의 총 갯수가 전체 갯수의 절반을 넘을 때 더하기를 중지한다. 결국 상자에 부여된 화소값을 중간값으로 사용할 수 있다. 이에 대해 실습을 통해 좀더 자세히 알아보기로 한다.


[02:33] 실제적 숫자를 가지고 예를 들어 살펴보기 위해 정상분포(Normal Distribution)로부터 1,000개의 난수를 생성해 두었다고 하자. BinApprox 알고리즘을 적용해보면 확률변수가 0인 평균을 찾는 것으로 끝난다. 그것이 바로 우리가 예상했던 평균과 중간값이 일치하는 것이다.
---------------------------------------------------------------
[BinApprox 알고리즘의 간략한 해설]
정규분포 확률에 따라 발생된 난수의 평균을 계산을 위해 '중간값'을 찾는 가장 손수운 방법

Q: 정규분포 확률로 발생된 난수의 평균은 중간값과 같다 ?

단계 1: 정상분포 확률을 따르는 난수 1000개 생성
단계 2: 생성된 난수를 오름차순으로 정렬.
단계 3: 배열된 값에 속한 수의 갯수(같은 값을 갖는 수의 갯수) 히스토그램 작성(Bining)
단계 4: 가장 작은 수부터 빈도 갯수 더하기 시작하여 갯수의 합이 500(=1,000/2)을 넘으면 중지.
단계 5: 덧셈을 중지 했을 때  값이 중간값

---------------------------------------------------------------


[02:47] 그럼 이 알고리즘을 겹쳐쌓기에 적용하면 어떻게 될까? 결과 영상은 이와 같다. 중심의 몇개 화소가 분명하게 드러나 보인다. 잡음에서 화소 몇이 드러난 이 결과가 그리 대단해 보이지 않을 수도 있으나 한발짝 물러서서 지금 보고 있는 것이 무엇인지 생각해보자. 개별 펄사를 전혀 인지할 수 없었던 대량의 영상을가지고 있었다. 그 영상들을 줄세워 놓고 보니 전 영상의 중간에 인지 할 수 없었던 펄사가 있었고 전체 영상에 대해 중간값 계산을 한 결과 뭔가 인지할 수 있는 것이 나왔다.

[03:15] 당신이 지금보고 있는 것은 통계적으로 펄사의 모습을 보게된 것이며 너무나 희미해서 원시 자료에서는 찾아볼 수 없었던 펄사의 등장이다. 보이지 않는 우주의 일단을 탐사하기 위해 진정 단순한 기법을 적용하였고 망원경으로 볼 수 없덨던 우주를 보게되었다.

--------------------------------------------------------------------
중간값 알고리즘에 대한 해설


--------------------------------------------------------------------
참고문헌

--------------------------------------------------------------------
파이썬 예제


1주차 요약

[커세라 강좌 소개]자료기반 천문학(Data-Driven Astronomy)
https://www.coursera.org/learn/data-driven-astronomy

----------------------------------------------
Week 1: Thinking about data
제1주차: 자료의 개념
- Principles of computational thinking
  전자계산학에 기초한 자료처리의 원리
- Discovering pulsars in radio images
  전파망원경 관측영상에서 펄서 찾기
----------------------------------------------
----------------------------------------------
2강: 강좌의 구성 / Lesson 2: Course Overview
----------------------------------------------
3강: 펄사(Pulsars) / Lesson 3: Pulsars
----------------------------------------------
4강: 관측영상 겹쳐쌓기 / Lesson 4: Diving In: Image Stacking
----------------------------------------------
5강: 난관에 봉착하다Lesson 5: The Challenge: What went wrong ?
----------------------------------------------
6강: 개선책을 찾다 / Lesson 6: The Solution: Improving your method
----------------------------------------------
1주차 요약: Module Summary (동영상 강좌) / 한글자막 / 영문자막

[강좌대본]

[00:06] 이번주 강좌는 펄사에 관한 질문 부터 시작 했었다. MWA에서 관측한 영상에서 전파를 방출하는 펄사들을 얼마나 많이 찾아 낼 수 있을까? 처음 조사에서는 아무것도 발견해 낼수 없었다. 하지만 "영상쌓기"라는 재치있는 기법을 사용하여 신호대 잡음비를 개선하고 통계적 영상처리 기법으로 잠음에 감춰졌던 펄사들의 검출이 가능하게 됐다.

[00:25] 이번 작업을 시작으로 아주 단순한 평균쌓기 프로그램을 작성해 봤고 이로부터 통계적으로 매우 탄탄한 중간값 쌓기라는 기법으로 확대발전 시켰다. 자료량이 많지 않다면 중간값 쌓기를 그대로 구현하기에 문제없다. 하지만 자료의 량이 늘어나자 즉각적으로 다른 문제들이 발생했다. 이는 현대과학(측정기술의 발달로 실험관측 자료가 넘쳐나는)이 안게된 당면 과제가 되었다.

[00:47] 심지어 스스로 컴퓨터와 별 상관 없다고 여기던 이론과학자들 조차 자료를 좀더 효율적으로 분석할 수 있는 방법을 모색하지 않을 수 없게됐다. 이미 수많은 중요 기술(라이브러리화 되어있음)들이 이미 개발되어 있고 여러분들이 당면한 문제를 푸는데 매우 유용한 도구들(프로그램이나 라이브러리)이 준비되어 있다는 것을 알아두기 바란다. 하지만 이번 강좌에서는 문제를 해결해 가는 과정에 집중하기로 한다. 이것이 바로 컴퓨터를 활용 사고(Computational Thinking)라 불리는 것이다.

[01:12] "컴퓨터 활용 사고"는 문제를 인식할 수 있게하며 컴퓨터에 적합하도록 사고의 틀을 바꿀 수 있게 하며, 가능한 해법이 무엇인지 생각케하고 아울러 서로다른 해법들 사이의 장단점을 파악케한다. 어떤 해법은 간단하지만 컴퓨터 메모리의 고갈 문제를 야기할 수 있다는 점을 경험했다. 이러한 컴퓨터 활용 접근방식은 여러분이 연구를 하는동안 접하게될 자료처리 문제에 적용될 수 있다.

[01:33] 자료와 그 처리를 위한 알고리즘의 중요성을 깨닳게 됐다는 점은 다음에 무엇을 하게되던 엄청난 강점을 갖추게 되었음을 의미한다.

2018년 6월 23일 토요일

1주/5강: 난관에 봉착하다

[커세라 강좌 소개]자료기반 천문학(Data-Driven Astronomy)
https://www.coursera.org/learn/data-driven-astronomy

----------------------------------------------
Week 1: Thinking about data
제1주차: 자료의 개념
- Principles of computational thinking
  전자계산학에 기초한 자료처리의 원리
- Discovering pulsars in radio images
  전파망원경 관측영상에서 펄서 찾기
----------------------------------------------
----------------------------------------------
2강: 강좌의 구성 / Lesson 2: Course Overview
----------------------------------------------
3강: 펄사(Pulsars) / Lesson 3: Pulsars
----------------------------------------------
4강: 관측영상 겹쳐쌓기 / Lesson 4: Diving In: Image Stacking
----------------------------------------------
5강: 난관에 봉착하다
Lesson 5: The Challenge: What went wrong ?(동영상 강좌) / 한글자막 / 영문자막

[강좌대본]

지난번 강의에서 일련의 관측영상의 평균쌓기(mean stack)를 수행해 봤다. 그결과 이런 영상을 얻을 수 있었다.

[00:13] 각각의 영상에서 아무것도 관측된 것이 없어 보였지만 모두 쌓아 평균을 해보니 신호대잡음의 비율이 높아져 마침내 숨어있던 펄사들이 모습을 드러냈다.
[00:22] 이제 잡음에 대한 강인성을 높일 방법을 찾아보자. 그 방법으로 쌓기 평균대신 중간값평균을 적용하는 것이다. 말하자면 영상에서 모든 화소값 크기별 순서를 세우고 중간 값의 화소를 기준으로 새로운 영상을 생성하는 방법이다.
[00:36] 이 방법을 그대로 컴퓨터 프로그램으로 작성하기는 그리 어렵지 않았다. 그러나 학생들이 60만 매에 이르는 영상을 대상으로 이 프로그램을 실행 시켰을 때 당장 문제가 생겼다. 이 프로그램이 우리 컴퓨터의 모든 지점(병렬 분산 컴퓨터 망의 각 CPU)에서 실행 되는 것이었다. 기억장소 전부를 차지하면서 예상했던 것보다 훨씬 더 많은 실행시간이 걸렸다. 말하자면 해결책을 가늠할 수 없었다.(끝이 안보였다)
[00:56] 이런 방식의 계산 이면에 무엇이 있는지 찾아보자. 뭔가 일을 시키기 전에 최선의 결과를 얻기위해 그 이면을 살펴보자. (요즘쓰는 말로 '악마는 디테일'에 있다.) 계산의 개념을 이해하면 계산과정에서 매 단계마다 왜 그렇게 하는지 타당성을 추정할 수 있으며 컴퓨터로 구현할 때 어느 정도의 성능과 저장장치 용량을 필요할지 가늠해 볼수 있다. 정확한 숫자가 중요한 것이 아니다. 계산을 수월하게 수행할 개략적인 값을 알고 시작할 수 있다.
[01:20] 나의 해법은 이렇다. 처음에 미지의 관측영상 60만장을 가지고 있다. 각 영상은 가로 200화소 세로 200화소 크기다. 한 영상당 총 4만 화소다. 그리고 각 화소마다 8바이트 크기의 배정도 부동소수점 값으로 저장되었다. 만일 본래의 중간값평균 계산을 적용한다면 이 영상의 화소 값들을 모두 저장장치에 올려 둬야 한다. 그러니까 200곱하기 200에 8을 곱하고 다시 60만을 곱한 저장용량이 필요하다. 계산하면 무려 192기가 바이트가 필요하다. 맙소사!

[01:51] 학생들이 모르긴 해도 이정도 메모리를 사용하는 프로그램을 짰을 것이다. 우리가 가진 (슈퍼)컴퓨터 장비의 전체 용량과도 맞먹는다. 이것이 바로 계산의 이면에 숨어있는 것이다. (디테일에 숨은 악마다.) 대용량의 자료 묶음을 다룰 때 그 일을 하기 전에 어떤 문제가 있을지 예상해 봐야한다. 어떤 방법으로 이 문제를 고치거나 아니면 아주 회피할 방법은 무엇일까?

[02:09] 이것이 바로 다음에 주어진 과제에서 풀어야할 일이다. 다음 강좌에서 여러분이 작성한 프로그램의 동작 시간을 재고 계산의 효율성과 컴퓨터 자원 사용량 규모을 평가해 보기로 하자.

-------------------------------------
효율적인 컴퓨팅을 위해 알고리즘을 분석하는 일은 매우 중요하다. 계산하고자 하는 수학식을 그대로 컴퓨터 언어로 옮기는 일은 무척 쉽다. 하지만 컴퓨터의 자원과 시간을 얼마나 소요할지 평가하지 않으면 그 프로그램이 결과를 내기까지 평생이 걸릴지도 모른다. 더구나 수행시간은 곧 에너지 소모를 의미하며 컴퓨터 자원은 비용을 요구한다.

우주 탐사장비에 사용되는 컴퓨터의 경우 열악한 우주환경에 견디려면 동작 속도가 매우 느리다. 에너지 또한 매우 제한적이다. 이런 컴퓨팅 환경에서 관측 데이터를 수집하고 분석하는 소프트웨어를 가동 해야 한다. 우주 탐사선에 구현된 각종 수학식이라고 특별하지는 않다.

곧이곧대로(naively) 프로그램을 짤 일이 아니다. 수학식의 계산 절차를 분석하고 반복횟수를 최적화해야 한다. "알고리즘 분석"이 필요하다.


[참고강좌] "알고리즘 분석"

[코멘트] "프로그래밍과 코딩"
과연 초등학교부터 '코딩'을 배워야 하는가?


-------------------------------------

2018년 6월 22일 금요일

1주/4강: 관측영상 겹쳐쌓기

[커세라 강좌 소개]자료기반 천문학(Data-Driven Astronomy)
https://www.coursera.org/learn/data-driven-astronomy

----------------------------------------------
Week 1: Thinking about data
제1주차: 자료의 개념
- Principles of computational thinking
  전자계산학에 기초한 자료처리의 원리
- Discovering pulsars in radio images
  전파망원경 관측영상에서 펄서 찾기
----------------------------------------------
----------------------------------------------
2강: 강좌의 구성 / Lesson 2: Course Overview
----------------------------------------------
3강: 펄사(Pulsars)Lesson 3: Pulsars
----------------------------------------------
4강: 관측영상 겹쳐쌓기
Lesson 4: Diving In: Image Stacking(동영상 강좌) / 한글자막 / 영문자막

[강의대본]

천문학에서 잡음이 깔린 자료에서 신호를 잡아내려고 시도하는 경우가 흔하다. 다음의 관측 영상에서 얼마다 많은 펄사들을 찾을 수 있을까? 머치슨 와이드필드 배열형 전파 망원경(MWA-Murchison Widefield Array)으로 관측한 영상이다.

MWA는 서호주에 위치한 저주파수 대역의 전파 망원경이다. 주파수가 80에서 300메가 헤르츠(Mhz)사이의 신호를 관측 할 수 있다. (대개 전파망원경이 수 GHz인 것에 비하면 매우 낮은 주파수 대역이다) 이 범위의 주파수는 대부분 라디오 방송국과 비슷하다. (FM방송국,아마추어 무선 통신,위성통신,항공통신 등) 이 전파 망원경의 매우 넓은 시각특성 덕분에 넓은 영역의 탐사 관측 프로젝트에 매우 유리하다. (전파 망원경의 특성:분해능은 낮지만 넓은 관측각)


이것이 MWA의 전형적인 관측 영상이다. 천체로부터 방출되는 에너지 밀도(플럭스)의 차이가 영상에서 회색의 명도 차이로 나타나 있다. 검을 수록 방출 에너지 밀도가 강하다는 듯이며 옅은 회색은 배경잡음 이다. 영상에서 보이는 검은 점들의 대부분은 멀리 떨어진 은하들로 전파를 방출 하고 있다. 몇개의 우리은하 내부의 전파방출 천체들도 보인다. 바로 펄사나 초신성 잔재들이다.


전파 천문학에서 플럭스 밀도의 측정치 단위는 잰스키(Jy)다. 면적당(미터제곱) 주파수 마다 통과하는 에너지의 량(와트)을 플럭스라 하는데 1 잰스키는 10의 마이너스 26승이다. 그러니까 플럭스 밀도란 망원경의 감지기(수신 안테나)에 입사하는 면적당 주파수별 신호의 세기(파워 스펙트럼)를 보여주는 것이다. 이에대해 좀더 알고 싶다면 이 동영상 강좌 편에 함께 제공된 자료를 읽어보기 바란다.

[참조] "마구잡이 수학", "1.4 요약 및 추가 연습문제"의 "문제 1.12 스테판-볼츠만(Stefan-Boltzmann) 법칙"

일단 이해를 돕기 위해 우리가 측정 하려는 것이 무엇인지 알려주고 싶다. 그것은 특정 주파수에서 어떤 펄사가 내는 신호의 강도(겉보기 밝기)라는 것이다. 먼저 주목할 것은 천체에서 나오는 빛(에너지 량)이 매우 미약하다는 점이다. 휴대전화에서 방출되는 전파의 세기와 비교해 보라.


천체관측 영상은 보통 FITS라는 파일의 형식으로 저장된다. 이 형식의 영상은 DS9이나 알라딘(Aladin) 같은 온-라인 도구를 통해서 볼 수 있다.

1. SAOImage DS9
2. Aladin Sky Atlas


연구 목적에 따라 관측영상을 가상의 색으로 표시하는데 전파는 색으로 표현할 수 있는 주파수가 아니라는 점을 염두에 두자. (전파는 가시광선이 아니다.) 이렇게 색깔별로 표현되어 지도의 등고선 처럼 보이는 영상은 단지 신호강도의 차이를 색으로 보여주기 위한 것이다. 지금 보여주는 영상은 우리의 커다란 전파관측 영상의 일부분을 따온 것이다.


왼쪽에 있는 영상의 중앙에 펄사가 관측된 것이 보인다. 오른쪽 영상에는 펄사가 관측되지 않았다. 우리가 뭔가 관측됐다고 말하려면, 해당 플럭스 밀도의 세기가 주변의 잡음보다 다섯배 이상의 표준편차를 가져야 한다. 만일 우리가 MWA의 관측 주파수 대역에서 전파 방출을 관측하고 있는데 펄사가 존재한다는 위치를 향하고 있음에도 대부분 관측되지 않다가 어쩌다 한번씩 관측되기도 한다.


만일 아무것도 관측하지 못했다면 매우 여러가지 원인이 있을 것이다. 펄사는 매우 먼 천체라는 것이고, 이 천체가 내뿜는 에너지 량이 엄청 나더라도 우리가 관측 하는 주파수 영역 밖에 있을 수도 있고, 깜빡이는 중간에 사진을 찍었을 수도 있다. 우리가 펄사의 신호를 감지하지 못한 정확한 이유를 대지 못할 지라도 감지하지 못했다는 사실에서 조차 뭔가 얻어낼 것이 있지 않을까?

어쨌든 천문학자들은 비록 한순간 감지하지 못할 수도 있지만 유효한 관측의 결과를 얻기 위해 현명한 방법을 개발했다. 그중 하나가 겹쳐쌓기(stacking)인데 어쩌다 감지하지 못했던 신호의 빈도를 확률적으로 찾는 방법이다.

* 이후 "겹쳐 쌓기"의 과정이 가우시안 무작위 잡음 제거(Gaussian Random Noise Reduction)효과를 설명 하는 부분에 약간 오류가 있어 보임.

[2:51] -----------------------
Stacking works because the noise in a radio image is roughly random, with a Gaussian distribution centered on zero. When you add regions of an image that just have noise, the random numbers cancel out. But when you add regions of an image in which there are signals, the signals add together, increasing what we call the signal to noise ratio.

겹쳐쌓기는 잡음이 매우 무작위(random)로 나타나는 전파영상을 처리할 때 효과적이다. 무작위 잡음의 진폭이 될 확률은 0을 중심으로 가우시안 분포를 갖는다. 여러장의 영상을 겹쳐 더하면 잡음만 있는 영역의 무작위 잡음 신호는 사라진다. 하지만 의미있는 신호가 포함된 영역의 경우 겹쳐 더할 수록 값이 증가하는데 이를 신호대 잡음의 비율(S/N Ratio)이라고 한다.
--------------------------------

이 효과를 시각적으로 보이기 위해 1차원 신호의 예를 들어 보자.

[3:15] --------------------------

Let's take a signal that looks like a single Gaussian and then add some random noise.
단일 가우시안처럼 보이는 신호가 있다고 하자. 여기에 임의의 무작위 잡음을 더해보자.
If we add enough noise, we can no longer see the signal.
잡음을 충분히 첨가하면 더이상 신호를 알아볼 수 없게된다.

Now say we have 100 of these signals, each with random noise added.
이제 이렇게 무작위 잡음를 첨가한 신호 100개를 준비하자.

If we take the mean of these, we can see the signal to noise ratio has increased
이 신호들을 모두 겹쳐 더한 후 평균을 구하면 신호대 잡음비가 향상된 결과를 얻을 수 있다.

and we can detect the underlying signal over the whole population of 100 sources.
100개의 원시 자료에서 빈도를 취합하여 잡음에 뭍혀있던 신호를 찾아낸 것이다.

-----------------------------------------------

이제 이 기법을 활용하여 우리가 원했던 펄사 검출 문제를 해결해 보자. 그전에 집고 넘어갈 것이 있다. 탐지되지 않은 펄사는 하늘 전역에 어디든 존재한다. 따라서 겹쳐쌓기를 하기 전에 먼저 그들의 위치를 조정해 줘야 하는데 찾고자 하는 펄사를 픽셀에 맞춰 중심에 둬야 한다. 겹쳐쌓기하는 과정은 대략 다음과 같다.


겹쳐쌓기의 평균을 계산하기 위해 영상속 전체 화소의 평균을 취한다. (한 영상 전체의 평균 밝기보정) 그결과 새로운 (정규화)영상이 만들어졌다. 다음 단계의 과정도 매 영상마다 처리해 주어야 한다. 학생들에게 주어진 과제의 끝은 여기까지다. 그리고 여기까지는 모든 과정은 아주 직관적이다.

하지만 간단해 보이는 추가사항이 하나더 주어져도 일은 아주 복잡해 질 수 있다. 추가사항 이란 전체평균(mean)대신 중간값평균(median) 겹쳐쌓기를 적용하라는 것이다. 중간값평균이 통계적으로 전체평균에 비해 매우 안정적(robust)하기 때문이다. 과학적으로 이것은 어떤 의미를 가질까? 아마 고교수학 과정에서 전체평균이 중간값평균에 비해 튀는 값에 더 쉽게 영향을 받는 다는 사실을 배웠을 것이다.

만일 대상 값의 분포가 대칭이라면 전체평균이나 중간값평균은 같게 나올 것이다.

하지만 비대칭 분포를 취하는 경우 특이값에 두 평균은 상당한 차이를 보인다. 중간값평균이 좀더 중심에 위치하게 되는 것을 볼 수 있다. 다음 강좌에서 이에 대해 좀더 자세히  살펴보기로 하겠다. 어찌 이렇게 단순해 보이는 한가지 변경으로 인해 결과를 얻기까지 계산의 수월함에 큰 영향을 주게될지 알아보기로 하자.

--------------------------------------------------
[파이썬 코딩 실습]

1. 파이썬의 배열형(Array)

#-------------------------------------------------
# Write your calculate_mean function here.
def calculate_mean(a):
  return sum(a)/len(a);

# You can use this to test your function.
if __name__ == '__main__':
  # Run your `calculate_mean` function with examples:
  mean = calculate_mean([1, 2.2, 0.3, 3.4, 7.9])
  print(mean)

  mean = calculate_mean([1.2, 3.8, 2.2, 8.2, 7.1])
  print(mean)
#-------------------------------------------------


2. 'NumPy' 모듈에서 배열 (CSV, Comma-Separated-Value 파일 읽기)

#-------------------------------------------------
import numpy as np

def calc_stats(filename):
  data = np.loadtxt(filename, delimiter=',')
  print(data)

  mean = np.mean(data)
  median = np.median(data)

  return np.round(mean, 1), np.round(median, 1)

# You can use this to test your function.
if __name__ == '__main__':
  result = calc_stats('data.csv')
  print(result)

  result = calc_stats('data2.csv')
  print(result)

  result = calc_stats('data3.csv')
  print(result)
#-------------------------------------------------

3. 'NumPy'의 배열형 자료 취급법
#-------------------------------------------------
# Write your mean_datasets function here
import numpy as np

def mean_datasets(filename):
  n = len(filename);
  if (n>0):
    data = np.loadtxt(filename[0], delimiter=',')
    for i in range(1,n):
      data += np.loadtxt(filename[i], delimiter=',')
   
    #Mean
    data_mean = data/n
 
    return np.round(data_mean,1)
  else:
    print('no file')
  return

# Run your function with examples from the question:
if __name__ == '__main__':
  print(mean_datasets(['data1.csv', 'data2.csv', 'data3.csv']))
  print(mean_datasets(''))
#-------------------------------------------------

4-1. 'AstroPy' 모듈에서 FITS 형식 파일 다루기 및 영상보기 방법
#-------------------------------------------------
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 19 11:10:45 2018
@author: GoodKook
"""

# Write your load_fits function here.
from astropy.io import fits
import numpy as np

def load_fits(filename):
  hdulist = fits.open(filename)
  data = hdulist[0].data

  arg_max = np.argmax(data)
  max_pos = np.unravel_index(arg_max, data.shape)

  return max_pos

if __name__ == '__main__':
  szFITSname = 'image3.fits'

  # Run your `load_fits` function with examples:
  bright = load_fits(szFITSname)
  print(bright)

  # You can also confirm your result visually:
  #from astropy.io import fits
  import matplotlib.pyplot as plt

  hdulist = fits.open(szFITSname)
  data = hdulist[0].data

  # Plot the 2D image data
  plt.imshow(data.T, cmap=plt.cm.viridis)
  plt.colorbar()
  plt.show()
#-------------------------------------------------

4-2. 여러 FITS 형식 파일의 평균(mean) 구하기

#-------------------------------------------------
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 19 12:12:13 2018
@author: GoodKook
"""

# Write your mean_fits function here:
from astropy.io import fits
import numpy as np

def mean_fits(files):
  n = len(files)
  if n > 0:
 
    hdulist = fits.open(files[0])
    data = hdulist[0].data
    hdulist.close()
 
    for i in range(1, n):
      hdulist = fits.open(files[i])
      data += hdulist[0].data
      hdulist.close()
 
    mean = data / n
    return mean

if __name__ == '__main__':
  # Test your function with examples from the question
  data  = mean_fits(['image0.fits', 'image1.fits', 'image2.fits'])
  print(data[100, 100])
  # You can also plot the result:
  import matplotlib.pyplot as plt
  plt.imshow(data.T, cmap=plt.cm.viridis)
  plt.colorbar()
  plt.show()
#-------------------------------------------------