text:
\documentclass{beamer}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loading packages that alter the style
\usepackage[farmed, numberded]{graphicx}
\usepackage[]{color}
\usepackage{alltt}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{dirtytalk}
\usepackage{ragged2e}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{physics}
\usepackage{listings}
\usepackage{epigraph}
\usepackage{blindtext} %This package generates automatic text
\usepackage[framed, numbered]{matlab-prettifier}
\setcounter{secnumdepth}{3}
\setcounter{tocdepth}{3}
\setlength{\parskip}{\smallskipamount}
\setlength{\parindent}{0pt}
\newtheorem{teo}{Teorema}[section]
\newtheorem{coro}{Corolario}[subsection]
\newtheorem{prop}{Proposición}[subsection]
\newtheorem{dfn}{Definición}
\newtheorem{examp}{Ejemplo}
\usetheme{Boadilla}
\usepackage{lipsum}
\usepackage{hyperref}
\hypersetup{hidelinks,linkcolor = black}
\usepackage{listings}

\raggedbottomhttps%://www.overleaf.com/project/61a28f32cde4ed5c68942eea

\newcommand\Fontvi{\fontsize{6}{7.5}\selectfont}

%%%%% Portada %%%%%%%
\title[Advanced MicroEconometrics]
{Estimating the Return to Schooling: Progress on Some Persistent Econometric Problems}

\author[Luis Martínez]
{Luis Gerardo Martínez Valdés}
\institute{ITAM}
\date{Spring 2022}
\logo{\includegraphics[height=.5cm]{itamlogo.png}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\frame{\titlepage}

\begin{frame}{Table of Contents}
\Fontvi
\tableofcontents
\end{frame}

%%% Sintaxis de las diapositivas %%%%%%
\section{Objetive}
\begin{frame}{Objective}
 \begin{itemize}
     \item To measure the \textbf{causal effect of education} on \textcolor{red}{labor market earnings} by using institutional features of the \textcolor{blue}{supply side} of educations system as an exogeneous determinants of schooling outcomes.

 \begin{block}{Key Finding}
Author made a review of studies that have used compulsory schooling laws, differences in the accessibility of schools, and similar features as \textbf{IVs} for completed education, reveals that the resulting estimates of the return to schooling are typically as big or bigger than the corresponding \textbf{OLS'} estimates.
\end{block}

    \item We will give a few explanations or interpretations (if you will) of why this happens.
 \end{itemize}
\end{frame}

\section{Introduction}
\subsection{Background}
\begin{frame}{Background}

   \begin{itemize}
       \item Resurgence of interest in the study of causal links between education and labor market success.
       \item Qualified vs Under-Qualified Workers.
       \item Revival in the determinants of economic growth.
       \item New focus on the role of human capital.
       \item Concerns about the relative costs and benefits of higher education for those who were not previously receiving it.
   \end{itemize}

\end{frame}

\subsection{Introduction}
\begin{frame}{Introduction}
Institutional features of the education system and supply-side variables.
\begin{itemize}
    \item How are they used?
    \item Why?
    \item What is their impact today?
\end{itemize}

\begin{block}{Aim}
Card aims to present a present a survey and partial synthesis of the recent literature that has used \say{supply-side} features of the education system to help identify the causal effect of education. By doing this, he will try to reconcile various findings in the literature, and also provide a useful framework for generating new hypotheses and insights about the connection between education and earnings.
\end{block}

\end{frame}

\subsection{Structure}
\begin{frame}{Structure}
\begin{itemize}
    \item Presentation of simple theoretical model of endogenous schooling.
    \item Use this model to to motivate an extended discussion of various econometric issues.
    \item Present a selective review of the recent literature on estimating the economic returns to education, drawing on studies of the U.S. and other developed economies, as well as a handful of studies of developing economies.
\end{itemize}

\end{frame}

\section{A Model of Endogenous Schooling}
\subsection{Framework}
\begin{frame}{Framework}
    Card builds on Becker's (1967) model of endogenous schooling, since most conceptual issues underlying the interpretation of recent studies of the return to education can be illustrated in the framework presented by the latter.

    \begin{itemize}
        \item Individuals face a market opportunity locus that gives the level of earnings associated with alternative schooling choices.
        \item Thus, they reach an optimal schooling decision by balancing the benefits of higher schooling against the costs.
    \end{itemize}

    \begin{block}{Strong Assumption}
        Traditionally, it is assumed that individuals seek to maximize the discounted present value of earnings, net of schooling costs

    \end{block}

\end{frame}
\subsection{Model}
\begin{frame}{Model}
    \begin{itemize}
        \item Assume that individuals have an infinite planning horizon that starts at the minimum school-leaving age ($t=0$).
        \item They accumulate a flow of utility in period $t$ that depends on:
        \begin{itemize}
            \item consumption $c(t)$ at period $t$
            \item Wether they are in school (and working part time) or out of school and working full time.
        \end{itemize}
        \item \textbf{Utility while in school is}:
                        \begin{equation*}
                                u(c(t))-\phi(t)
                        \end{equation*}
        \item \textbf{Utility out of school is}:
                        \begin{equation*}
                                u(c(t))
                        \end{equation*}
    \end{itemize}

\end{frame}

\begin{frame}{}
    \begin{block}{Note:}
        Note that $u(\cdot)$ is an increasing concave function and $\phi(t)$ is a convex function that represents \textcolor{blue}{relative disutility of school vs. work for the t-th year of schooling}.
    \end{block}
 Finally,
 \begin{itemize}
     \item Individuals discount future utility flows at a subjective discount rate $\rho$, and make a once-for-all decision on when to leave school.
 \end{itemize}

 Therefore, one can define \textbf{lifecycle utility, conditional on schooling $S$} and a given consumption profile as
 \begin{equation*}
     V(S,c(t))=\int_{0}^{S}\Big[ u(c(t)-\phi(t))\Big]e^{-\rho t} dt + \int_{S}^{\infty}u(c(t))e^{-\rho t} dt
 \end{equation*}

\end{frame}

\begin{frame}{Intertemporal Budget Constraint}
\begin{itemize}
    \item Let $y(S,t)$ denote real earnings at age $t$ of an individual who has completed $S$ years of post-compulsory schooling (with $0\leq S \leq t$).
    \item Assume that individuals who are in school at time $t$ work part time and earn $p(t)$ and pay tuition costs of $T(t)$.
    \item Also, assume that individuals can borrow or lend freely at a \textbf{fixed} interest rate R.
\end{itemize}

Thus, we can define the \textbf{intertemporal budget constraint}:
\begin{equation*}
    \int_{0}^{\infty}c(t)e^{-Rt}dt= \int_{0}^{S}\Big[ p(t)- T(t)\Big]e^{-R t} dt + \int_{S}^{\infty}y(S,t)e^{-R t} dt
\end{equation*}

\end{frame}

\begin{frame}{Optimal Schooling choice and Optimal Consumption Path}

Therefore, an individual’s optimal schooling choice and optimal consumption path maximize
\begin{align*}
    \Omega(S,c(t),\lambda) &  =  V(S,c(t)) \\
                           &  - \lambda \Bigg[ \int_{0}^{\infty}c(t)e^{-Rt}dt~- \int_{0}^{S}\Big[ p(t)- T(t)\Big]e^{-R t} dt \\
                           &  - \int_{S}^{\infty}y(S,t)e^{-R t} dt \Bigg]\\
\end{align*}

We can see that the FOC with respect to $S$ yields:
\begin{equation*}
    \Omega_{S}(S,c(t),\lambda)=\lambda e^{-RS}[MB(S)-MC(S)]
\end{equation*}
\end{frame}

\begin{frame}{}
    Where,
    \begin{itemize}
        \item The Marginal Benefit of the $S$-th unit of schooling is defined by
        \begin{equation*}
            MB(S) = \int_{0}^{\infty} \frac{\partial y(S,S+\tau)}{\partial Se^{-R \tau}} d \tau
        \end{equation*}
        \item And, the Marginal Cost of $S$-th unit of schooling is given by
        \begin{equation*}
            MC(S) = y(S,S) - p(S) + T(S) +\frac{1}{\lambda e^{-(\rho - R)S}}\phi(S)
        \end{equation*}
    \end{itemize}
\begin{block}{What happens if $MC(S)$ rises faster than $MB(S)$ ?}
We would need a necessary and sufficient condition for the optimal schooling choice $S$, i.e. $MC(S)=MB(S)$
\end{block}
\end{frame}

\begin{frame}{}
Card solves this problem by assuming that log earnings are additively separable in education and years of post-schooling experience. Therefore, one can write the earnings function as $y(S,t)= f(S)h(t-S)$ (recall that $0 \leq S \leq t$), then
\begin{equation*}
    MB(S)= f'(S) \int_{0}^{\infty} h(\tau) e^{R \tau}d\tau \equiv f'(S) H(R)
\end{equation*}

 where $H(R)$ is a decreasing function of the interest rate. Therefore, 2 conclusions arise:
 \begin{enumerate}
     \item If earnings are fixed after completion of schooling: $H(R)=\frac{1}{R}$
     \item If earnings follow a \textit{concave} lifecycle profile and taking g as a constant growth rate(equivalent to the lifecycle profile): $H(R)=\frac{1}{(R-g)}$

 \end{enumerate}

 Therefore, under separability(check Minkowski's Theorem), the marginal costs and benefits of additional schooling are equated when
 \begin{equation*}
     \frac{f'(S)}{f(S)}= \frac{1}{H(R)}\Bigg[ 1+ \frac{(T(S)-p(S))}{f(S)} + \frac{1}{\lambda e^{-(\rho - R)S}}\frac{\phi(S)}{f(S)}\Bigg]
 \end{equation*}
\end{frame}

 \begin{frame}{General Case: $U(c(t)) = log[c(t)]$}
     The FOC's for an optimal consumption profile, together with the lifecycle budget constraint, imply that
     \begin{equation*}
         \frac{1}{\lambda}= \rho \Bigg\{ e^{-RS}f(S) H(R) + \int_{0}^{S} \Big[ p(t) -T(t)\Big] e^{-R t} dt \equiv \rho W(S) \Bigg\}
     \end{equation*}

where $W(S)$ is the value of lifecycle wealth associated with the schooling choice S. If $T(S) \approx p(S)$ then
\begin{equation}
    \frac{f'(S)}{f(S)}= R - g +\rho e^{-\rho S}\phi(S) \equiv d(S)
\end{equation}

\begin{block}{Observations}
 The interest rate must be adjusted to reflect lifecycle earnings growth and the marginal cost has to account for the relative disutility of attending the Sth year of schooling.

\end{block}
\end{frame}

\begin{frame}{Heterogeinity in the Optimal Schooling Choice}
    Individual heterogeinity can arise from 2 sources:
    \begin{enumerate}
        \item differences in the economic benefits of schooling,
        \item differences in the marginal costs of schooling.
    \end{enumerate}

Then we would have
% Equation 2
\begin{equation}
    \frac{f'(S)}{f(S)}= b_{i} - k_{1}S,~~~ k_{1}>0
\end{equation}
%Equation 3
\begin{equation}
    d(S) = r_{i} + k_{2}S, ~~~k_{2}>0
\end{equation}

where $b_{i}\sim F(\bar{b},\sigma_{b}^{2})$, $r_{i}\sim F(\bar{r},\sigma_{r}^{2})$. Equating (2) with (3) and solving for $S$ will yield,
%Equation 4
\begin{equation}
    S_{i} =(b_{i}-r_{i})/k,~~~k=(k_{1}+k_{2})>0
\end{equation}

\end{frame}

\begin{frame}{}
At the equilibrium level of schooling described by equation (4) individual i’s marginal return to schooling is
\begin{equation*}
    \beta_{i} \equiv b_{i} - k_{1}S = b_{i}(1-\frac{k_{1}}{k}) + r_{i}\frac{k_{1}}{k}
\end{equation*}

The average marginal return to education is therefore,
\begin{equation*}
    \bar{\beta} = \mathbb{E}[\beta_{i}]= \mathbb{E}[b_{i}-k_{1}S]= \bar{b} + k_{1}\bar{S}
\end{equation*}

\begin{block}{Interpretation of equation (4)}
 Card interprets (4) as partial equilibrium description of the relative education choices of a cohort, given the institutional environment and economic conditions that prevailed during their late teens and early twenties.
\end{block}
\end{frame}

\section{Econometric Issues Raised By Endogenous Schooling}
\subsection{OLS Estimates of the Return to Schooling}
\begin{frame}{Econometric Issues Raised By Endogenous Schooling: OLS Estimates of the Return to Schooling}
Note that equation (2) implies a model for log earnings of the form
%%%% EQUATION 5
    \begin{equation}
        log~y_{i}= a_{0} + \bar{b}S_{i}-\frac{1}{2}k_{1}S_{i}^{2} + a_{i} + (b_{i}-\bar{b})S_{i},
    \end{equation}

where $a_{i} \equiv \alpha_{i}-a_{0}$ has mean $0$. Together with some algebraic manipulation andconsidering the linear projectos of the underlying random variables $a_{i},~~~ b_{i}$ and $r_{i}$, one can rewrite expression (5) as

 \begin{equation}
        log~y_{i}= K + (\bar{b}+\lambda_{0}-\psi_{0}\bar{S})S_{i} + (\psi_{0} - \frac{1}{2}k_{1})S_{i}^{2} + u_{i} + v_{i}S_{i},
    \end{equation}

\end{frame}

\begin{frame}{}
\begin{block}{}

Equations (4) and (6) together describe a two-equation system for schooling and earnings in terms of the underlying random variables $a_{i}, b_{i}$, and $r_{i}$.
Ignoring other covariates (or assuming these have already been conditioned out) it is straightforward to derive the implications of this model for conventional OLS estimates of the return to schooling.
\end{block}
 Moreover, under the assumption that the \texbf{third central moment of the joint distribution of $b_{i}$ and $r_{i}$ we have that $\mathbb{E}[(S_{i}-\bar{S})^{3}]=0$}. Then, the probability limit of the OLS regression coefficient $b_{OLS}$ from a regression of equation (6) is

 \begin{equation}
     plim~~b_{OLS}= \bar{\beta} + \lambda_{0} + \psi_{0}\bar{S}
 \end{equation}

If there is no heterogeinity in MB(S) and that the log earnings are linear in schooling (i.e. $k_{1}=0$). Then (7) implies
\begin{equation*}
    plim~~b_{OLS} - \bar{\beta}= \lambda_{0}
\end{equation*}
\end{frame}

\begin{frame}{}
\begin{itemize}
    \item Generally, people with a higher return to education have an incentive to acquire more schooling, a cross-sectional regression of earnings on schooling is likely to yield an upward-biased estimate of the average marginal return to schooling.

    \item With a larger bias the more important are the comparative advantage incentives that lead individuals with higher returns to schooling to acquire mor schooling.

    \item Problems arising from the measurment of schooling (Griliches 1977) Downward bias in the OLS estimate of the effect of schooling on earnings vs. Upward bias of the effect of ability on earnings.
\end{itemize}

\end{frame}

\subsection{IV Estimates of the Return to Schooling}
\begin{frame}{IV Estimates of the Return to Schooling}
\begin{block}{Trends}
 Recently, much attention has focused on supply-side sources of variation in schooling, attributable to such features as the minimum school-leaving age, tuition costs, or the geographic proximity of schools.
\end{block}

Suppose that the marginal cost component $r_{i}$ is \textit{linearly independent} realted to a \textbf{set of observable variables} $Z_{i}$: $r_{i}= Z_{i}\pi_{1} + \eta_{i}$, where $\eta_{i}$ takes into account other unobserved taste and cost factors; also, $\eta_{i}$ is uncorrelated with $Z_{i}$.

The Optimal Schooling Choice is

\begin{equation} \tag{4'}
S_{i} = \frac{(b_{i}-r_{i})}{k}  =   \frac{\bar{b}}{k}-\frac{Z_{i}\pi_{1}}{k} + \frac{(b_{i}-\bar{b}-\eta_{i})}{k}  =   Z_{i}\pi + \xi_{i}

\end{equation}

\end{frame}

\begin{frame}{}
\begin{itemize}
    \item If $a_{i}$ is the only individual-specific component of ability, then equations (4') and (5) constitute a standard simultaneous equations system.
    \item It is sufficient to suppose  $\mathbb{E}[a_{i}Z_{i}]=0$, to ensure that an IV estimator based on $Z_{i}$ will yield a consistent estimate of the average return to schooling $\bar{b}$.
    \item So, we have that $Z_{i}$ is independent from individual abilities and the reduced form of schooling residual $\xi_{i}$
\end{itemize}

Then,
\begin{align*}
   \mathbb{E}[log~~y_{i}~~|~~Z_{i}] & =  \mathbb{E}\big[a_{0} + \bar{b}S_{i}-\frac{1}{2}k_{1}S_{i}^{2} + a_{i} + (b_{i}-\bar{b})S_{i}~ | ~Z_{i}\big]\\
                                    & = a_{0}+ \bar{b}Z_{i} - \frac{1}{2}k_{1}(Z_{i}\pi)^{2} - \frac{1}{2}k_{1}\mathbb{E}[\xi_{i}^{2}~|~Z_{i}]\\
                                    & + \mathbb{E}[(b_{i}-\bar{b})\xi_{i}~|~Z_{i}]
\end{align*}

 \begin{itemize}
     \item Thus, the average marginal return to schooling can be consistently estimated by IV. Why?
 \end{itemize}

\end{frame}

\begin{frame}{Treatments}
\begin{block}{Note:}
Unfortunately, the assumptions proposed previously and weaker ones proposed by Woolridge (1997), are likely to be violated when $Z_{i}$ is a variable representing exposure to institutional structures on the supply-side of the education system.
\end{block}

Assume that the joint distribution of abilities and tastes ($a_{i}.b_{i},r_{i}$) is the same for individuals who attended the reformed schools (i.e. $Z_{i}$=1) and those who did not (i.e. $Z_{i}$=0); however, for reformed schools the optimal choice of schooling is given by
\begin{equation}\tag{4''}
    S_{i}= \frac{(b_{i}-\theta r_{i})}{k};~~\theta<1
\end{equation}

\end{frame}

\begin{frame}{}
    Take $r_{i}= \bar{r}+\eta_{i}$; then,
    \begin{itemize}
        \item Unreformed Schools: $S_{i}^{u}= \pi_{0} + \xi_{i0}$,
        \item Reformed Schools: $S_{i}^{R}= \pi_{0} + \pi_{1} + \eta_{i1}$,
        \item Reduced form schooling equation is therefore
              \begin{equation*}
                  S_{i}= \pi_{0} + Z_{i}\pi_{1} + \xi_{i}~~; ~~~\xi_{i} = (1-Z_{i})\xi_{i0} + Z_{i}\xi_{i1}
              \end{equation*}
    \end{itemize}

    \begin{block}{Note:}
     Since the schooling reform lowers the effect of cost differences in the optimal schooling decision, $Var[\xi_{i}~|~Z_{1}] \leq~Var[\xi_{i}~|~Z_{0}]=0 $.

    \end{block}
\end{frame}

\begin{frame}{}
    Some evidence that changes in the institutional structure of the education system affect the mapping between ability and schooling outcomes:
    \begin{figure}
        \centering
        \includegraphics[scale= .5]{Table1.png}
    \end{figure}

    \begin{block}{Conclusion}
    Assuming the presence of nearby college is uncorrelated with ability, \textbf{college proximity is a potential IV for schooling}.

    \end{block}
\end{frame}

\subsection{Alternatives to IV}
\begin{frame}{Alternatives to IV}
    \begin{itemize}
        \item Control Function Approach(Garen, 1984).
        \begin{block}{Basic Idea}
        Make some assumptions about the nature of the covariances between the unobserved ability components $a_{i}$ and $b_{i}$ and the observable variables $S_{i}$ and $Z_{i}$, and include additional terms in the earnings model that capture these relationships.
        \end{block}

        \item Assume mean-independence between unobserved ability and taste components and Z; also, he conditional expectations of the unobserved ability components $a_{i}$ and $b_{i}$ are linear in the schooling residual:
        \begin{equation}\tag{8a}
           \mathbb{E}[(b_{i}-\bar{b})~|~S_{i},Z_{i}]= \psi_{0}\xi_{i},
        \end{equation}
        \begin{equation}\tag{8b}
           \mathbb{E}[a_{i}~|~S_{i},Z_{i}]= \lambda_{0}\xi_{i}.
        \end{equation}
    \end{itemize}
\end{frame}

\begin{frame}{}
  Together with eqaution (5), this equations imply that
  \begin{equation*}
      \mathbb{E}[log~y_{i}~|~Z_{i}]=a_{0}+\bar{b}S_{i}-\frac{1}{2}k_{1}S_{i}^{2} + \lambda_{0}\xi_{i} + \psi_{0}\xi_{i}S_{i}
  \end{equation*}

 This equation can be estimated by a two-step procedure in which the estimated residual \hat{\xi_{i}} from the reduced form schooling equation is substituted for $\xi_{i}$, and is substituted for $\xi_{i}$, $\hat{\xi_{i}}S_{i}$ is substituted for $\xi_{i}S_{i}$.

 \begin{itemize}
     \item Assumption (8a) IS PROBLEMATIC! Why?
     \begin{itemize}
         \item $Cov[b_{i},\xi_{1}~|~Z_{i}]$ and $Var[\xi_{i}~|~Z_{i}]$ potentially vary with $Z_{i}$
     \end{itemize}
     \item A simple extension of the control function approach may be appropiate if $Z_{i}$ is an indicator variable. Remember the college proximity idea?
     \item Taking this into account, and with some algebraic manipulation of equations (8a), (8b) then
     \begin{align*}
      \mathbb{E}[log~y_{i}~|~Z_{i}] & =a_{0}+\bar{b}S_{i}-\frac{1}{2}k_{1}S_{i}^{2} + \lambda_{00}\xi_{i} + (\lambda_{01}-\lambda_{00})Z_{i}\xi_{i} + \psi_{00}\xi_{i}S_{i} \\      & + (\psi_{01}-\psi_{00})Z_{i}S_{i}\xi_{i}
  \end{align*}
 \end{itemize}
\end{frame}

\begin{frame}{A more radical approach}
Maximum likelihood estimation of a structural model of earnings and schooling, based on a complete specification of the unobservable components in the earnings function and the utility function.

\begin{block}{Edge}
An advantage of this approach is that the earnings function can be made quite general; for example, by allowing the returns to different years of schooling to vary in a flexible manner with individual ability.
\end{block}

\end{frame}

\subsection{What Does IV Estimate?}
\begin{frame}{What Does a Conventional IV Estimate?}
Suppose that:
\begin{itemize}
    \item A given individual would have schooling level $S_{i}^{c}$ and earning $y_{i}^{c}$ if he or she attended regular school system.
    \item If he or she attended a reformed school, then the individual would have a schooling outcome of $S_{i}^{c} + \Delta S_{i}$
\end{itemize}

Let $\beta \equiv$ individual $i$'s marginal return to schooling. Then the effect of schooling reform on earnings for individual $i$ is
\begin{equation*}
    \Delta log~y_{i}= \beta_{i}\cdot \Delta S_{i}
\end{equation*}

Implying,
\begin{align*}
   (\ast)~~ plim~b_{IV}~= \frac{Cov[log~y_{i},Z_{i}]}{Cov[S_{i},Z_{i}]} & = \frac{\mathbb{E}[log~y_i~|~Z_{i}=1]-\mathbb{E}[log~y_{i}~|~Z_{i}=0]}{\mathbb{E}[S_{i}~|~Z_{i}=1]-\mathbb{E}[S_{i}~|~Z_{i}=0]}\\
    & = \frac{\mathbb{E}[\beta_{i}\cdot \Delta S_{i}]}{\mathbb{E}[\Delta S_{i}]}
\end{align*}
\end{frame}

\begin{frame}{}
\begin{block}{Note: }
If marginal return to schooling and treatment of schooling ($\Delta S_{i}$), then the IV estimator is a consistent estimate of the avg. marginal return to eduction $\bar{\beta}=\mathbb{E}[\beta_{i}]$
\end{block}

\begin{itemize}
    \item An IV procedure based on a school reform that leads to bigger changes in the education choices of people with relatively high marginal returns to education will tend to produce an over-estimate of the average marginal return to education.
    \item Important to reconsider that the effects of a supply-side change that causes a porportional reduction in MC(S). This induced change is represented by
    \begin{equation*}
        \Delta S_{i}= r_{i}\frac{1-\theta}{k}= \bar{r}\frac{(1-\theta)}{k} + \eta_{i}\frac{(1-\theta)}{k}>0
    \end{equation*}
    Thus, the monoticity assumption required for LATE is satisfied.
\end{itemize}

\end{frame}

\begin{frame}{}
Then, individual $i$'s marginal return to schooling (in the absence intervention) is
\begin{equation*}
    \beta_{i}= \bar{\beta} + (b_{i}-\bar{b})(\frac{1-k_{1}}{k}) + \eta_{i}\frac{k_{1}}{k}
\end{equation*}
    Susbtituting in equation ($\ast$),
    \begin{equation*}
        plim ~b_{IV}= \bar{\beta} + \Bigg[ \sigma_{\eta}^{2} \frac{k_{1}}{k} + \frac{\sigma_{b \eta}}{\bar{r}}(\frac{1-k_{1}}{k})\Bigg]
    \end{equation*}

Two features are worth noticing:
\begin{enumerate}
    \item The probability limit of the IV estimator is unaffected by classical measurement error in schooling.
    \item The validity of a particular IV estimator depends crucially on the assumption that the instruments are uncorrelated with other latent characteristics of individuals that may affect their earnings.
\end{enumerate}
\end{frame}

\section{IV Estimates of the Return to Schooling}
\begin{frame}{IV Estimates of the Return to Schooling}
\begin{figure}
    \centering
    \includegraphics[scale= .42]{Table2.1.png}
\end{figure}
\end{frame}

\begin{frame}{}
\begin{figure}
    \centering
    \includegraphics[scale= .4]{Table2.2.png}
\end{figure}
\end{frame}

\begin{frame}{}
\begin{figure}
    \centering
    \includegraphics[scale= .4]{Table2.3.png}
\end{figure}
\end{frame}

\section{Conclusions}
\begin{frame}{Conclusions}
   \begin{itemize}
       \item We have attempted to measure the causal effect of education on labor market earnings by using institutional features on the supply side of the education system as exogenous determinants of schooling outcomes.
       \item Card believes it is helpful to place the returns to education literature in a standard \say{supply and demand} framework, which leads to a somewhat richer econometric model for schooling and earnings than is usually adopted in the applied literature.
       \begin{block}{NOTE: }
       Different individuals finish their schooling at a point where the marginal return to the last unit of education may be either above or below the average marginal return in the population as a whole.

       \end{block}
   \end{itemize}
\end{frame}

\begin{frame}{Conclusions}
   \begin{itemize}
       \item IV estimation will typically recover a weighted avg. of returns to education for people whose choices were affected by the instrument, rather than the marginal return to education in the population.
       \item IV estimates of the return to schooling that are at least as big and sometimes substantially bigger than the corresponding OLS estimates.
       \item In many cases the IV estimates are relatively imprecise, and none of the empirical strategies is based on true randomization.
       \begin{itemize}
           \item Thus, no individual study is likely to be decisive in the debate over the magnitude of ability biases in OLS estimates of the return to schooling.
       \end{itemize}
       \item Seems that Griliches' (1977) assessment was spot on.
   \end{itemize}

\end{frame}
\section{Bibliography}
\begin{frame}{Bibliography}
    \begin{itemize}
        \item Card, David. \say{Estimating the return to schooling: Progress on some persistent econometric problems}. \textit{Econometrica} 69, no. 5 (2001): 1127-1160.
    \end{itemize}
\end{frame}

\end{document}

max_stars_repo_name: LMtzV/AdvancedMicroeconometrics

max_issues_repo_name: LMtzV/AdvancedMicroeconometrics

max_forks_repo_name: LMtzV/AdvancedMicroeconometrics


text:
\documentclass{article}
\usepackage[utf8]{inputenc}
\usepackage{graphicx}
\usepackage{amsmath}
\graphicspath{ {./Images/} }
\usepackage{float}

\title{CTA200 Computing Project}
\author{Caleb Lammers}
\date{May 19 2020}

\begin{document}

\maketitle

\section*{Introduction}

The light emitted from a galaxy provides information about its physical properties, including both
its star formation history and stellar mass. Light emitted from galaxies in different wavelengths
contain different information about the galaxy's composition. For instance, it is well known that
younger stars emit shorter wavelength light; that is, bluer light than old stars. Thus, if a galaxy
emits much more shorter wavelength light, this galaxy likely has many younger stars and thus was forming
stars not so long ago. Because different wavelengths help provide different insights, it is useful
to represent a galaxy's emitted light as a graph of luminosity/flux density vs. wavelength, called
the galaxy's Spectral Energy Distribution (SED).

\begin{figure}[h]
  \centering
\includegraphics[scale=0.4]{SED Example}
\caption{Example SED of flux density vs. wavelength}
\end{figure}

By combining knowledge on stellar evolution, calibrated stellar libraries, initial mass functions
(IMF), and dust in the interstellar medium, it is possible to predict the SED of galaxies, based on
their physical properties. This is known as stellar population synthesis (SPS); a model that
combines these ingredients is called a stellar population synthesis model. In this project,
the FSPS (Conroy and Gunn 2009) SPS model was used, particularly the python package python-FSPS.
Being able to compute a galaxy's SED based on its physical properties is a step in the right direction,
however, as astronomers, we would like to be able to determine a galaxy's physical properties from
its observable SED. This project will consist of two parts. In the first, the FSPS model was set up
and functions were created so that SEDs and other important graphs could be made. In the second,
the goal was to work backwards and determine the physical properties of two galaxies, given their SEDs, which was done using
brute force and Markov chain Monte Carlo (MCMC) methods.

\section*{Description of Script}
\subsection*{Part 1}

A StellarPopulation object was initialized in accordance with the settings provided in question 4.
Now, by using the built in get spectrum method, SEDs of (bolometric) luminosity in $L_\odot/Hz$
vs. wavelength in $\AA$ can be made.

\begin{figure}[h]
  \centering
\includegraphics[scale=0.3]{SED Default Units}
\caption{SED of luminosity vs. wavelength}
\end{figure}

SEDs are more commonly presented as flux density (in $\mu Jy$) vs wavelength (in $\AA$), so, next a
function was written to convert the SED from units of luminosity into units of spectral flux,
according to the equation below.

\begin{figure}[H]
\[\textrm{f}_v = \frac{\textrm{L}_{(1 + z)v}}{4\pi D_L^2} \]
\caption{Flux density from luminosity}
\end{figure}

Because $\textrm{f}_v$ depends on $(1 + z)v$ rather than $v$, instead of converting each $\textrm{f}_v$
into $\textrm{L}_v$, it is easier to convert $\textrm{f}_v$ to $\textrm{L}_{(1 + z)v}$
and scale each $\lambda$ by $(1 + z)$. Now, SEDs of this galaxy can be created in units of spectral flux by specifying a redshift. Note
that at all redshifts we expect the SED to have similar behaviour as above because we have only
rescaled the x and y axes.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.3]{SED redshift 1}
\caption{SED of flux density vs. wavelength at redshift 1}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.3]{SED redshift 2}
\caption{SED of flux density vs. wavelength at redshift 2}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.3]{SED redshift 3}
\caption{SED of flux density vs. wavelength at redshift 3}
\end{figure}

Next, the goal was to fix $\lambda = 8000 \AA$ and make a plot of flux density vs. redshift. Note
that at each z the generated SED is not guaranteed to have a flux density value at $\lambda = 8000 \AA$. To
get around this, the wavelength closest to $8000 \AA$ was linearly interpolated to find a reasonable
value for flux density, giving the following graph. As you would expect, flux density decreases as redshift increases.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.3]{flux density vs redshift}
\caption{Flux density vs. redshift at $\lambda = 8000 \AA$}
\end{figure}

Finally, a function to generate the
SED in flux density ($\mu Jy$) vs. wavelength ($\AA$) for a galaxy
of any age at any redshift was created. Using this function, the following graphs were made.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED age 0.01 redshift 1}
\caption{SED of flux density vs. wavelength for galaxy of age 0.01 Gyr at redshift 1}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED age 0.1 redshift 1}
\caption{SED of flux density vs. wavelength for galaxy of age 0.1 Gyr at redshift 1}
\label{fig:9}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED age 1 redshift 1}
\caption{SED of flux density vs. wavelength for galaxy of age 1 Gyr at redshift 1}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED age 3 redshift 1}
\caption{SED of flux density vs. wavelength for galaxy of age 3 Gyr at redshift 1}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED age 10 redshift 1}
\caption{SED of flux density vs. wavelength for galaxy of age 10 Gyr at redshift 1}
\end{figure}

\subsection*{Part 2}

In order to fit a straight line to the given data, it was assumed that the error is Gaussian
so that the use of weighted linear least-square fitting is justified (see answer to 2.1 for more motivation). Using linear
algebra, a function to determine the slope and y-intercept (along with uncertainties) was made.
Additionally, a function to calculate $\chi ^2$ was created. Then, the line may be plotted and
its fit may be quantitatively measured.

$$b = -7.390^{+2.158}_{-2.158}$$
$$m = 5.277^{+0.451}_{-0.451}$$

\begin{figure}[H]
  \centering
\includegraphics[scale=0.3]{data and line}
\caption{Optimal line to fit the given data. $\chi ^2 = 17,847$}
\end{figure}

Similarily, in order to estimate the age and redshift of the galaxy, we must have a way to quantify
the goodness of fit of a theoretical SED to a noisy observed SED. By the analogy to fitting a
straight line, a similar likelihood function was used:

\begin{figure}[H]
\[Y ^2 = - \sum_{i=1}^{N = 1} \frac{(y_i - f(x_i))^2}{\sigma_{i}^2} \]
\caption{Likelihood function for the fit between theoretical and observational SEDs}
\label{fig:14}
\end{figure}

Where now $y_i$ (Fig.~\ref{fig:14}) represents an observed flux density and $f(x_i)$ is the theoretical flux density.
Note that the negative sign is required because a higher likelihood must imply a closer fit. With
this definition in place, we can now compare how a galaxy with age 0.1 Gyr at z = 1.0 compares with
galaxy 1 and galaxy 2's SEDs. This galaxy had a likelihood of -3121221711061.225 with galaxy 1 and a
likelihood of -364732.5275292772 with galaxy 2.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED galaxy 1}
\caption{given SED of galaxy 1}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED galaxy 2}
\caption{given SED of galaxy 2}
\end{figure}

%Include figure
Clearly, a galaxy with age 0.1 Gyr (Fig.~\ref{fig:9}) and z = 1.0 has a more similar SED to galaxy 2, which makes sense
when comparing the SED of galaxy 2 to the SED of this galaxy (above, include figure). By using the
likelihood function, we may also try to brute force our way to an estimate of the proper age and z
of these galaxies. That is, we may loop through a number of combinations of age and z to find which combination
has the highest likelihood. Using this implemented function, we get that galaxy 1 has age of around
6.4 Gyr and redshift of 0.8 and galaxy 2 has age of around 3.4 Gyr and redshift 0.2.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED closest galaxy 1}
\caption{SED of closest galaxy to galaxy 1 found by brute force}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.25]{SED closest galaxy 2}
\caption{SED of closest galaxy to galaxy 2 found by brute force}
\end{figure}

As you can see, the graphs are somewhat similar but have their discrepancies. Before continuing, one
should note that this is not a great way of approximating the age and redshifts of these galaxies.
In fact, we have no guarantees that we are even close to correct, we may have simply found a local
max of the likelihood function. To do this in a more intelligent way, Markov chain Monte Carlo
(MCMC) methods were used. First, MCMC was applied to the example of fitting a line. Of course,
applying MCMC to the the $\chi ^2$ model is pointless because we have an analytic solution. Instead,
we now adopt a more complex model with the assumption that the errors were somewhat underestimated.

\begin{figure}[H]
\[\ln({p}) = -\frac{1}{2} \sum_{i=1}^{N = 1} \frac{(y_i - m(x_i) - b)^2}{s_i^2} + \ln({2\pi s_i^2})\]
\[ s_i^2 = \sigma_{i}^2 + f^2 (m(x_i) + b)^2\]
\caption{Likelihood function for new model}
\end{figure}

Notice that this model has a so-called nuisance parameter $f^2$. Because it is now not clear that
there is an analytical solution or even how f, m, and b should
relate, MCMC is useful. To use MCMC, we define a uniform prior probability function, and use python
emcee's features. With 32 walkers centered around the values found from weighted linear least-square fitting and
10000 steps we get the following Markov chain.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.3]{Markov Chain}
\caption{Markov chain with emcee}
\end{figure}

According to the autocorrection time, it takes about 40 steps for the Markov chain to burn in.
Removing the first 40 steps, we get the following corner plots.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.3]{Corner Graphs}
\caption{Markov chain with emcee}
\end{figure}

These plots show that b and m are closely related while b and f/m and f/b are not. The fact that f
takes on low values shows that the the errors are likely not underestimated by very much, if at all. Nonetheless, MCMC gives the following values for b, m, and f and graph.

$$m = 5.261^{+0.444}_{-0.455}$$

$$b = -7.391^{+2.172}_{-2.108}$$

$$f = 0.052^{+0.056}_{-0.037}$$

\begin{figure}[H]
  \centering
\includegraphics[scale=0.2]{Data and line 2}
\caption{Markov chain with emcee}
\end{figure}

The slope, y-int and uncertainties of this line are similar to the those from the weighted linear
least-square fitting line. A similar method was used to find optimal age and redshift for galaxy
1 and 2. Again, a uniform prior was used and the same likelihood function from above was used. When run with 2500 steps and 16 walkers centered around the estimated ages and redshifts,
the walkers did not move around the space - they got stuck around the estimated values. So, instead,
the walkers were started in random locations distributed around the estimated values. Here are the Markov
chains and corner plots.

\begin{figure}[H]
  \centering
\includegraphics[scale=0.5]{Markov chain Galaxy 2}
\caption{Markov chain of galaxy 1 data with emcee}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.5]{Corner Graph Galaxy 1}
\caption{Corner graph of galaxy 1 data with emcee}
\end{figure}

$$ln{\textrm(age)} = 1.608^{+0.105}_{-0.002}$$

$$z = 0.853^{+0.016}_{-0.024}$$

Note that the error is likely somewhat higher than calculated here (which is due to how the distribution looks).

\begin{figure}[H]
  \centering
\includegraphics[scale=0.5]{Markov chain Galaxy 2}
\caption{Markov chain of galaxy 2 data with emcee}
\end{figure}

\begin{figure}[H]
  \centering
\includegraphics[scale=0.5]{Corner Graph Galaxy 2}
\caption{Corner graph of galaxy 2 data with emcee}
\end{figure}

$$ln{\textrm(age)} = 1.843^{+0.559}_{-0.001}$$

$$z = 0.162^{+0.000}_{-0.029}$$

Due to the difficulty of calculation, the Markov chain and corner plots were not as continuous but it is
interesting how they still find a similar answer to the brute force method

The next step would be to add stellar mass and dust attenuation to the free variables in MCMC.
Though it would be easy to extend the code to do this, the MCMC calculations above already
took several hours, so, this would undoubtedly take unreasonably long in the current setup. To make
this feasable, one could use the Scinet cluster or try to further optimize the code.

\subsection*{Question Answers}

\paragraph{1.0.} As briefly mentioned in the introduction, stellar population synethsis (SPS) refers to
combining information on stellar evolution, calibrated stellar libraries, initial mass functions
(IMF), and dust in the interstellar medium to construct the SED of a galaxy. Many different SPS models exist,
and they  have their own assumptions, tools, and calibrations. By using an SPS model, specifying a
galaxy's physical properties allows for reconstructiion of its SED. Of course, this is only half the battle
because we would like to go from the SED to a galaxy's physical properties. There have been many approaches to
solving this problem, including Monte Carlo methods, machine learning, and numerical methods. No matter
the tool used, the role of the SPS model is the same. The SPS model will generate SEDs given specified physical
parameters which will be compared to the SED of the desired galaxy (there are many ways to "compare," of course).
Then, adjustments will be made to the physical parameters in hopes of finding a better fit between the generated
and desired SED. In this way, SPS models provide the foundation for this method of determining a galaxy's
physical properties from its emitted light.

\paragraph{1.3.} The StellarPopulation can take in a variety of inputs about the galaxy, including its dust
content, initial mass function, metallicity, the type of star formation history and more.

\paragraph{1.4.} As discussed above, the units of the output spectrum are solar luminosity per hertz
($L_\odot/Hz$). The StellarPopulation object is using MESA Isochrones and Stellar Tracks (MIST) and
the MILES spectral library. At 100 Myr, the stellar mass of the galaxy is about 0.92 $M_\odot$.

\paragraph{1.7.} As galaxies cannot be older than the universe, when a galaxy is observed at some redshift
$z$, we are effectively looking back in time $t$ at its light. Then, $t + t_g < t_u$ (there are equations
to calculate $t$ from $z$) where $t_u$ is about 13.8 billion years. So, the larger a galaxy's redshift,
the smaller its possible maximum age is.

\paragraph{2.1} There are many possible techniques to fit a straight line to a set of sampled
points. As Hogg, Lang, and Bovy point out, the standard method when there is negligible error in the x
data and known Gaussian error in the y data is weighted linear least-square fitting. This type of fitting minimizes $\chi ^2$ which can be intuitively understood as a measure of
the line from the data. It is mathematically defined as:

$$\chi ^2 = \sum_{i=1}^{N = 1} \frac{(y_i - f(x_i))^2}{\sigma_{yi}^2}$$

Part of the motivation for using this measure of goodness of fit is that it has an easily calculable
minimum with linear algebra. With this method, we get a reasonable looking line of best fit. It is
difficult to say whether a line is the best curve to represent the data without any physical
knowledge about the data, however, a line appears to fit the data well with a relatively low $\chi
^2$.

\end{document}

ext: tex

max_stars_repo_path: Final Project/writeup.tex

max_issues_repo_path: Final Project/writeup.tex

max_forks_repo_path: Final Project/writeup.tex


text:
\documentclass[12pt]{article}
\usepackage{a4wide,alltt,xspace,times}
\usepackage{skaff}
\usepackage{optionman}
\newcommand{\Size}[1]{|#1|}
\newcommand{\GPC}{\texttt{gpcluster}\xspace}

\author{Stefan Kurtz\thanks{\SKaffiliation}}
\title{\textbf{Finding Clusters of Close Genepairs}\\[2mm]
       \textbf{Manual}}
\begin{document}
\maketitle
The program \GPC reads files specifying gene locations and pairs of
duplicated genes. It clusters two duplicated genes, if the gene locations
are close together.

\section{Program Options}
The program is called as follows:
\par
\noindent\GPC [\emph{options}]
\par
The options of \GPC are as follows:

\begin{Justshowoptions}

\Option{genefile}{$\Showoptionarg{genefile}$}{
Specify the file containing gene specifications of the following form:
\begin{alltt}
geneid loc frompos topos str
\end{alltt}
\texttt{geneid} is the unique identifier of the gene,
\texttt{frompos} and \texttt{topos} are the coordinates
(start an end position) of the gene, while \texttt{loc} is a
location identifier (i.e.\ a string without white spaces).
\texttt{str} is \texttt{f} (for forward) or \texttt{r} (for reverse)
denoting the direction. The coordinates are relative to the location specified.
The columns in the given $\Showoptionarg{genefile}$ are separated by
white spaces. The file may contain comment lines
starting with the symbol \texttt{\symbol{35}} in the first column.
This option is mandatory.
}

\Option{genepairfile}{$\Showoptionarg{genepairfile}$}{
Specify the file containing lines with pairs of duplicated genes. Each
line is of the following form:
\begin{footnotesize}
\begin{alltt}
genespec$_{1}$ genespec$_{2}$ Evalue
\end{alltt}
\end{footnotesize}
where \texttt{genespec}$_{1}$ and \texttt{genespec}$_{2}$
both specify a pair of duplicated genes (each in five columns) called
genepair in the following.
%(geneid$_{1}$ < geneid$_{2}$)
%There are 46 overlaps between consecutive genes.
An assessment of their
similarity is given by an E-value in column 10. The positions and
locations given for two genes must be consistent with the information
in the corresponding gene specification file.
\Showoptionarg{genepairfile} may contain comment lines
starting with the symbol \texttt{\symbol{35}} in the first column.
This option is mandatory.
}

\Option{maxevalue}{$\Showoptionarg{maxevalue}$}{
Specify a maximum E-value threshold $\Showoptionarg{maxevalue}$
to discard all genepairs with E-value larger than the threshold.
}

\Option{gapsize}{$\Showoptionarg{gapsize}$}{
Specify the main clustering-parameter $\Showoptionarg{gapsize}$. We
perform single linkage clustering. Two genepairs
\((\mathit{g1},\mathit{g2})\) and
\((\mathit{g3},\mathit{g4})\) are put into one cluster,
if and only if the following is true:
\begin{itemize}
\item
\(\emph{g1}.\mathit{loc}=\emph{g3}.\mathit{loc}\)
\item
\(\emph{g2}.\mathit{loc}=\emph{g4}.\mathit{loc}\)
\item
If \(\emph{g1}.\mathit{frompos}\leq
     \emph{g3}.\mathit{frompos}\) then
\(\emph{g1}.\mathit{topos}+gapsize\geq
  \emph{g3}.\mathit{frompos}\).
\item
If \(\emph{g3}.\mathit{frompos}\leq
     \emph{g1}.\mathit{frompos}\) then
\(\emph{g3}.\mathit{topos}+gapsize\geq
  \emph{g1}.\mathit{frompos}\).
\item
If \(\emph{g2}.\mathit{frompos}\leq
     \emph{g4}.\mathit{frompos}\) then
\(\emph{g2}.\mathit{topos}+gapsize\geq
  \emph{g4}.\mathit{frompos}\).
\item
If \(\emph{g4}.\mathit{frompos}\leq
     \emph{g2}.\mathit{frompos}\) then
\(\emph{g4}.\mathit{topos}+gapsize\geq
  \emph{g2}.\mathit{frompos}\).
\end{itemize}

If this option is not specified, then the default gap size is 0.
}

\Option{loose}{$\Showoptionarg{gapsize}$}{
The standard clustering strategy links a new gene pair to an existing
cluster if there is one gene pair in the existing cluster within the
prescribed gap size of the new gene pair.  This is tight clustering.
If this option is used, then tight clustering is followed by
loose clustering according to the given \Showoptionarg{gapsize}.
Loose clustering iterates the same clustering procedure
with new units being the current clusters with their collective start and
end positions (replacing the original gene pairs as units). The iteration
continues until no further links can be made. The output will be in the
same format as when this option is not used.
}

\Option{outfile}{$\Showoptionarg{filename}$}{
Specify a filename prefix for the output file.
If the appropriate option is chosen, then the singleton genepairs are output
in file \emph{filename}\texttt{.}\texttt{single}.
If the option \Showoption{multioutfile}
is specified, then the clusters are output in files named
\emph{filename}\texttt{.}\emph{clsize}\texttt{.}\emph{clnum}, where
\emph{clsize} is the size of cluster number \emph{clnum}.
If the option \Showoption{multioutfile}
is not specified, then the clusters are output in a single file
\emph{filename}\texttt{.}\texttt{all}.
}

\Option{separator}{$c$}{
Specify a character to seperate columns in the output. The default is a blank.
To, for example, use a tabulator as separator, specify
\Showoption{separator}~\texttt{\symbol{39}\symbol{92}t\symbol{39}}.
}

\Option{mincls}{$\Showoptionarg{mincls}$}{
Specify the minimal size of the clusters to be output. If this option is not
used, then the default value for \Showoptionarg{mincls} is 2.
If $\Showoptionarg{mincls}=1$, then also singleton genepairs are output.
This option requires the use of the option \Showoption{outfile}.
}

\Option{maxcls}{$\Showoptionarg{maxcls}$}{
Specify the maximum size of the clusters to be reported.
If this option is not used, then the default value for \Showoptionarg{maxcls}
is \(\infty\), i.e., there is no upper bound on the size of clusters to be
reported. The value of $\Showoptionarg{maxcls}$ must be larger or equal to
$\Showoptionarg{mincls}$. This option requires the use of the option
\Showoption{outfile}.
}

\Option{multioutfile}{}{
Output each cluster in an extra file. If this option is not used,
then all clusters are output in a single file.
This option requires the use of the option \Showoption{outfile}.
}

\Option{help}{}{
Show a summary of all options and terminate.
}

\end{Justshowoptions}

\noindent
Each of the options can be specified only once.

\section{Output Format}
Each cluster of close genepairs is output as follows:
\begin{small}
\begin{alltt}
# cluster 1312 of 9 genepairs
At3g58360 3 21470537 21471898 r At3g58400 3 21479699 21481095 r 1e-86
At3g58360 3 21470537 21471898 r At3g58380 3 21474997 21476227 r 1e-65
At3g58360 3 21470537 21471898 r At3g58410 3 21481903 21483261 r 4e-77
At3g58350 3 21468650 21470261 r At3g58410 3 21481903 21483261 r 1e-62
At3g58350 3 21468650 21470261 r At3g58400 3 21479699 21481095 r 4e-85
At3g58350 3 21468650 21470261 r At3g58360 3 21470537 21471898 r 3e-94
At3g58340 3 21466103 21467433 r At3g58350 3 21468650 21470261 r 1e-66
At3g58340 3 21466103 21467433 r At3g58360 3 21470537 21471898 r 6e-73
At3g58340 3 21466103 21467433 r At3g58410 3 21481903 21483261 r 1e-83
# matching ranges: (3,21466103,21471898) vs (3,21468650,21483261) with overlap
\end{alltt}
\end{small}
The first line shows the cluster number and the number of genepairs in the
cluster. The genepairs in clusters are shown in the same format
as the input. Columns are tab-separated.
The last line shows a minimal region containing all genepairs. The
attribute \texttt{with overlap} appears if for the matching ranges
$(l_{1},p_{1},q_{1})$ and $(l_{2},p_{2},q_{2})$ we either have
\(p_{1}\leq p_{2}\leq q_{2}\) or
\(q_{1}\leq p_{2}\leq q_{2}\).

\section{Example}
The file \texttt{ATpep.in} contains 25,545 lines with gene specfications.
Here are the first five
lines:
\begin{alltt}
At1g01010 1 3760   5630  f
At1g01020 1 7729   8666  r
At1g01030 1 11864  12940 r
At1g01040 1 23519  31079 f
At1g01050 1 31382  32670 r
\end{alltt}
The file \texttt{ATpairs.in} contains
71,846 lines with genepairs. The first five lines are as follows:
\begin{alltt}
At1g01010 1 3760 5630 f At1g02230 1 433031 436775 r 1e-41
At1g01010 1 3760 5630 f At4g01550 4 673025 675225 r 7e-49
At1g01020 1 7729 8666 r At4g01510 4 641683 643190 f 2e-34
At1g01030 1 11864 12940 r At1g13260 1 4542386 4543420 f 8e-40
At1g01030 1 11864 12940 r At1g25560 1 8981886 8982971 r 1e-36
\end{alltt}

The program call
\begin{alltt}
\GPC.x -genefile ATpep.in -genepairfile ATpairs.in
            -gapsize 10000 -maxevalue 1e-60
            -outfile Cluster/cl -mincls 1
\end{alltt}
gives to the following output:
\begin{small}
\begin{alltt}
# read gene specifications in file "ATpep.in"
# sort 25545 gene specifications
# read genepairs with E-value <= 1e-60 in "ATpairs.in"
# sorting 16664 genepairs:
# tight clustering of 16664 genepairs: gapsize is 10000
# output clusters of size >= 1
# output singlets
# cluster statistics:
# 1772 clusters
# 4912 elements out of 16664 (29.48%) are in clusters
# 11752 elements out of 16664 (70.52%) are singlets
# 1213 clusters of size 2
# 295 clusters of size 3
# 122 clusters of size 4
# 48 clusters of size 5
# 37 clusters of size 6
# 11 clusters of size 7
# 6 clusters of size 8
# 9 clusters of size 9
# 8 clusters of size 10
# 6 clusters of size 11
# 3 clusters of size 12
# 2 clusters of size 13
# 2 clusters of size 14
# 1 cluster of size 15
# 3 clusters of size 16
# 1 cluster of size 18
# 1 cluster of size 20
# 1 cluster of size 22
# 1 cluster of size 28
# 2 clusters of size 29
\end{alltt}
\end{small}
It shows the steps of the computation and the distribution of the
cluster sizes. Additionally, files
\texttt{Cluster/cl.all} and
\texttt{Cluster/cl.single} are
produced containing the selected clusters and all singleton genepairs.

\section{Things to be done}
\begin{enumerate}
\item
The direction of the genes is not considered in the clustering. This should be
removed. Ask Volker how to incorporate this into the clustering.
\item
Volker, is the clustering criteria as you had it in mind?
\item
Define other clustering criteria, like overlap. This may be specified by
a negative gap size.
\end{enumerate}
\end{document}

max_stars_repo_name: gordon/vstree

max_issues_repo_name: gordon/vstree

max_forks_repo_name: gordon/vstree


text:
\documentclass[11pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[english]{babel}
\usepackage[english]{isodate}
\usepackage[paper=a4paper]{geometry}
\newgeometry{top=3.5cm,bottom=2.5cm,right=2.5cm,left=2.5cm}
\usepackage{graphicx}
\usepackage{comment}
\usepackage{fancyhdr}
\usepackage{framed}
\usepackage{lastpage}
\usepackage[hidelinks]{hyperref}
\usepackage{tabularx}
\usepackage[table]{xcolor}
\usepackage{enumitem}
\usepackage{mdwlist}
\usepackage{placeins}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{xcolor}
\usepackage{listings}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINIZIONI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\titolo}  {Assignment 3}
\newcommand{\nome}    {Matrix Library}
\newcommand{\versione}{2.0}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SETUP DOCUMENTO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\input{advancedAl_setup.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SEZIONI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section {Introduction}
The goal of this document is to describe strategies adopted to solve this third assignment. In particular, the request of this one is to extend the matrix library, developed until now, so that the operations are performed in parallel when it is possible. The execution of multi operations in parallel has the advantage of improving the performance, but there also some particular issues that will be discussed during the following document.

\section{Problem Analysis}
In some cases matrix operations can be performed parallelly. For this reason providing an architecture that takes advantage from this could bring to optimal results in terms of performance. Here we present some cases in which operations can be performed in parallel:
\begin{itemize}
	\item $(A+B) + (C+D)$
	\item $(A+B) * (C+D)$
	\item $A*B*C*D$
\end{itemize}
On the first case we consider the sum between two matrices produced by other sum expressions, the idea of parallelism, in this case, consists in performing the sum $(A+B)$ and $(C+D)$ at the same time since they are independent and it does not change the final result. When the two parallel sums are finished, the main thread has to combine the two results and execute the final sum. The second case is very similar to the first one, but with the difference that at the end the main thread has to perform the matrix product of the two results instead of a sum.
The third case is a little bit difficult, since now there is not special order given by the usage of parentheses so we rely on the product optimization realized on the previous assignment. The idea is that $A*B*C*D$ can be optimized by the product optimizer such that the final expression could be for instance $(A*B)*(C*D)$, and in this case we reduce the problem to the second case. \\

The second requirement for this assignment is to apply another improvement to the matrix product, such that each matrix $A$ is composed of several submatrices $A_{ij}$ and the product $C= A*B$ can be expressed like:
$$C_{ij} = \sum_k A_{ik}*B_{kj}$$
Where $*$ represents the product operator between submatrices. The advantage of this technique is that each component $C_{ij}$ can be computed independently, taking advantage of possible parallelism algorithms. Note that $k$ is an arbitrary constant that define the size of each block of the matrix.
So at the end of this assignment matrix product operator should be able to apply multiple product in parallel and each pair of expression can be performed in parallel using the previous formula.
\newpage
\section{Implementation}
In this section we are going to discuss about code implementation, class responsibilities, motivations and limits behind different design choices.

\subsection{Design}
During the design phase of the library we thought about several solutions, but some of them did not satisfy all the requirements of the project or they had strong drawbacks.\\
For instance our first idea was to encapsulate a future object inside the \verb|matrix| class. In this way the internal structure of a matrix is substituted with the resulting matrix only during the evaluation. The problem of this solution was to define correctly the evaluation time. We decided to make transparent to the user as much as possible technical knowledge of the library, in order to provide a readable and clear syntax and semantics.

\begin{lstlisting}[language=C++]
matrix<int> out = A + B
matrix<int> ev = (A+B).eval()
\end{lstlisting}
On the first case the expression is by evaluated hiding technical aspects of library implementation. Instead the second one requires the knowledge of certain behaviors inside the library and it is more verbose. But this choice brought to some problems for the previous solution, since it was necessary to define a strategy that converts a \verb|matrix| object to another \verb|matrix| of the same type; but if the original matrix has a future it must be evaluated first. Our idea was to implement this strategy inside the move constructor but doing that we encountered the problem of \textbf{copy/move elision}. From \textit{C++11} it was introduced an optimization to C++ compiler that omits copy and move constructors when they are unecessary. For instance, when a temporary object with the same type should be created, it is constructed directly in the location where it would otherwise be copied or moved. In this case if the result of $A+B$ is a matrix object it will be created directly inside the final object, without calling move or copy constructor. At the end it was necessary to apply eval operation to the final matrix. There is the possibility to remove this optimization by setting a specific flag, but the resulting performance was not so good. For this reason we decided to adopt another strategy, based on the separation of complete matrix and future one, discussed on the following subsections.

\subsubsection{Parallel sum}
In order to improve the performance of sum operator, we decided to adopt a strategy that can take advantage of parallel expressions.
Considering for instance the following expression:
$$(A+B) + (C+D)$$
We can see that $A+B$ and $C+D$ can be performed parallelly, and when results are ready, the sum of the two resulting matrices is computed. By doing this, a problem arose: how can we distinguish the previous expression from this one: $A+B$?
At compile time we have no idea of the total structure of the expression, since they are evaluated considering possible parentheses and operators from left to right. In this case the first term that the compiler can see is $(A+B)$ which have no difference with respect to $A+B$, so it will execute the classic sum operation without performing any parallel operations. So it is necessary to introduce a temporary object that accumulates or keeps track of the expression seen before and perform operation when it is necessary.
For this reason we have introduced the concept of \textbf{future\_matrix}.
\newpage
\begin{lstlisting}[language= C++]
template<typename T>
class future_matrix<T, Plain> : public matrix_ref<T, Plain> {
	...
	matrix_ref<T, Plain> eval() {
	    return future.get();
    };

	std::shared_future<matrix<T>> future;
};
\end{lstlisting}
A future\_matrix is like a "promise", it will return the result of an expression when the \verb|eval| method is called. Its role is to encapsulate a future operation and when a caller asks for the evaluation the future matrix waits until the operation is not completely performed and later on returns the final resulting matrix. It is used by operator function that, instead of returning the final result will return a "matrix that will give the final result".  Just to provide an example of usage we want to introduce the sum operator between dynamic matrices:

\begin{lstlisting}[language= C++]
operator + (const matrix_ref<T,LType>& left, const matrix_ref<U,RType>& right) {
	...
	/*Create future */
	typedef matrix<typename op_traits<T, U>::sum_type> matrix_type;
	std::shared_future<matrix_type> ft =
			std::async(std::launch::async,
			sum<T, LType, U, RType, matrix_type>,
			left,right)
	.share();
	/* return future dynamic matrix */
	return future_matrix<typename op_traits<T, U>::sum_type, Plain>(std::move(ft), left.get_height(), left.get_width());
}
\end{lstlisting}
As we can see the sum operator creates the future operation \verb|ft| that will be used by the future matrix in order to synchronize eval operation with function \verb|sum|. The function \verb|sum| performs the sum between two matrices but it is assigned to an asynchronous thread that will deal the operation. When an entity asks for the evaluation of the operation to the future matrix, it can happen that for some reasons (like scheduling or complexity of the operation) the asynchronous thread has not completed the operation. In this case is essential for the future matrix to synchronize its operations with the future object returned by the asynchronous thread and wait until the operation is not completed.\\
Very similar is the behavior of this sum operator, but now it is applied to a standard dynamic matrix and a future one:
\begin{lstlisting}[language=C++]
operator+(const matrix_ref<T, LType> &left, future_matrix<U, Plain> &&right) {
	auto sum=[](const matrix_ref<T,LType>&left, future_matrix<U,Plain> &&right){
		matrix_ref<T,Plain> right_mat = right.eval();
		return sum<T,LType,U,Plain,matrix_type>(left, right_mat);
	};
	std::shared_future<matrix_type> ft =
	std::async(std::launch::async, sum, left, std::move(right)).share();

	return future_matrix<typename op_traits<T, U>::sum_type, Plain>(std::move(ft), left.get_height(), left.get_width());
}
\end{lstlisting}
We can see that more or less the two sum operators have the same structure, but in this case, since the right matrix is a \verb|future_matrix|, it is necessary to request its evaluation for performing the operation. The following example can explain better what is happening on this function:
$$A + B + C$$
Considering the above expression, in a first moment the evaluation of the sum between $A$ and $B$ is performed. During the evaluation, a brand new \verb|future_matrix| which returns the result of $A+B$ is created. We define as $t_1$ the thread responsible for the sum between $A$ and $B$. The expression $future\_matrix(A,B) + C$ is immediately evaluated, now the new thread $t_2$ must synchronize its operation with the thread $t_1$, since before applying the final sum, $t_2$ requires the result of the previous sum. This synchronization is performed using the eval of the \verb|future_matrix|, that will block the execution of $t_2$ until $t_1$ has completed the operation.\\

Methods discussed until now are defined and used for dynamic matrices but we know that the library allows the usage of static matrices, which provide information about their dimension also at compile time. For this reason there are two version of \verb|future_matrix|:
\begin{lstlisting}[language=C++]
class future_matrix<T, Plain>: public matrix_ref<T, Plain>{ ... }

class future_matrix<T, Sized<h, w>>: public matrix_ref<T, Sized<h, w>>{ ... }
\end{lstlisting}
The first one is used for dynamic operations and the second one for static operations. This allows also to maintain, when it is possible, static information inside the \verb|future_matrix|. In any case the two types of \verb|future_matrix| are converted to dynamic or static complete matrices when they are assigned to an object of \verb|matrix| type. When the assignment is performed, the \verb|future_matrix| is evaluated and after that, data is moved from the result matrix to the destination object.
\begin{lstlisting}[language=C++]
matrix(future_matrix<T, Plain> &&X) {
	matrix<T> result = X.eval();
	height = result.get_height();
	width = result.get_width();
	data = std::move(result.data);
}
\end{lstlisting}

A good question could be: Why don't use a \verb|matrix_wrap|?
It can be reasonable to use a \verb|matrix_wrap| and defining a unique \verb|future_matrix| that uses it. But our choice was to avoid the usage of \verb|matrix_wrap| since we know that it introduces some overhead that can affect drastically the performance of sum operator.  \\

We decided to define \verb|future_matrix| as a subclass of \verb|matrix_ref|. The reason is strictly related to the possibility of the library to be extended to include other more complex operations like:
$$(A+B)^T$$
For the moment the result of $A+B$ must be stored or evaluated explicitly and only after this, the transpose operator can be used. So, this expression can be computed with this code: $(A+B).eval().transpose()$, but it is necessary to give knowledge about the library to the user. It can be interesting to extend primitive operations/transformations of matrices hiding this technical knowledge to users. For instance, the \verb|transpose| method of a \verb|future_matrix| has to evaluate first the future operation, after that, the transpose operation can be applied to the result.\\

The general behavior can be summarized with the following schema:
\image{img/parallelSum.png}{Parallel sum schema}{0.6}

\subsubsection{Parallel product}\label{parallel_prod}
The second operation that can be performed using parallel execution is the multiplication between matrices. As we anticipated before there are essentially three cases that can be summarized with this two examples:
\begin{itemize}
	\item $(A+B) * (C+D)$
	\item $A*B*C*D$
\end{itemize}
The first case is an extension of the parallel sum and it is solved simply defining this function:
\begin{lstlisting}[language=C++]
operator*(future_matrix<T, LType> &&left, future_matrix<T, RType> &&right) {
	...
	return left.eval() * right.eval();
}
\end{lstlisting}
So it is necessary to apply the standard product to the two evaluated matrices. The product is not performed until the expression on the left and on the right are completed.\\
The general behavior of parallel sum is so summarized in the following image:
\image{img/parallelProduct1.png}{Parallel sum}{0.5}

For the second case we know from the previous assignment that this sequence of products can be optimized considering the size of the matrices. Optimization reduces the number of iteration necessary to compute the result. Another characteristic of product operation is that it is performed using \verb|matrix_wrap| and for this reason it is not possible to use the structures defined before with \verb|future_matrix|, since it is not possible to get directly the right type of matrix wrapped in a \verb|matrix_wrap|. For this reason we decided to take advantage of the structure defined from the \verb|matrix_wrap| and define a new implementation \verb|pimpl| that follows the same idea of a \verb|future_matrix|. In order to do that we have defined a new virtual method \textit{eval} inside the \verb|matrix_wrap_impl| interface. The role of \textit{eval} method is to return an implementation to a physical and evaluated matrix. All the standard implementations return the current implementation, instead the new future implementation uses the evaluation of the future operation and returns a pointer to the new implementation.

\begin{lstlisting}[language=C++]
class concrete_matrix_wrap_impl<T, Future> : public matrix_wrap_impl<T> {
	...
	std::unique_ptr<matrix_wrap_impl<T>> eval() const override {
		matrix<T> mat = future.get();
		return std::make_unique<concrete_matrix_wrap_impl<T, Plain>>(future.get());
	}
private:
	std::shared_future<matrix<T>> future;
};
\end{lstlisting}
The new \textit{pimpl} allows now the \verb|matrix_wrap| to wrap not physical/complete matrices but also future matrices, produced by product operation.

\begin{lstlisting}[language=C++]
matrix<T> do_multiply(matrix_wrap<T> lhs, matrix_wrap<T> rhs) {
	assert(lhs.get_width() == rhs.get_height());

	lhs.eval();
	rhs.eval();
	...
}
\end{lstlisting}
The function \verb|do_multiply| perform first the evaluation of the two matrices given in input. This allows to decide correctly the right implementation. If at least one of the two wraps a future, its pimpl will be substituted with the result, meaning that \verb|do_multiply| is synchronized and has to wait for the evaluation. For example, imagine that the \textit{lhs} matrix wraps a future and the \textit{rhs} matrix is complete. Assume also that the thread responsible for the operation of \textit{lhs} is identified by $t_1$. During the evaluation of \textit{lhs}, the current thread has to wait that $t_1$ completes the operation. When the result is ready the new \textit{pimpl} is created and now \textit{lhs} is a complete matrix. Evaluation of \textit{rhs} returns the \textit{pimpl} to itself. This structure allows also to maintain the same code developed until now and include only the parallelism with the usage of threads. The thread development is added to the method \verb|resolve_one|.
\begin{lstlisting}[language=C++]
void resolve_one() {
	typename std::list<matrix_wrap<T>>::iterator lhs = find_max();
	typename std::list<matrix_wrap<T>>::iterator rhs = lhs;
	++rhs;

	auto run_multiply = [this](matrix_wrap<T> lhs, matrix_wrap<T> rhs) {
		return do_multiply(lhs, rhs);
	};

	std::shared_future<matrix<T>> ft =
		async(std::launch::async,run_multiply, std::move(*lhs),std::move(*rhs));

	matrix_wrap<T> mat(std::move(ft),lhs->get_height(),rhs->get_width());
	typename std::list<matrix_wrap<T>>::iterator result=matrices.emplace(lhs,mat);

	matrices.erase(lhs);
	matrices.erase(rhs);
}
\end{lstlisting}
We can see that the future is wrapped in a \verb|matrix_wrap| that at the end is pushed inside the expression list. Later on, when it will be picked-up in order to perform a new operation it is first evaluated inside the \verb|do_multiply| method and the result is provided and repeat the procedure. Only the last step does not push inside the optimization list and returns directly the final result.

\image{img/parallelProduct2.png}{Parallel product}{0.65}

%%%%%

\subsubsection{Block multiplication}\label{block_structure}
The second task of the assignment requires the implementation of a parallel
matrix multiplication. This can be obtained by dividing the matrix in blocks
of size $K$ (where $K$ is an arbitrary constant); every pair of block-matrices is
multiplied in a dedicated thread. The provided implementation extends the one
described in Section \ref{parallel_prod}.\\ Once the block size $K$ has been
chosen, the number of row-blocks and column-blocks that compose the resulting
matrix can be calculated by taking the maximum size resulting from the division
by $K$:
\[
    \#\textit{rows}=\max\left(\left\lceil \frac{H_r}{K}\right\rceil, \left\lceil \frac{H_l}{K}\right\rceil \right),~\#\textit{cols}=\max\left(\left\lceil \frac{W_r}{K}\right\rceil, \left\lceil \frac{W_l}{K}\right\rceil \right)
\]
where:
\begin{itemize}
    \item $H_l$ is the left hand matrix height.
    \item $W_l$ is the left hand matrix width.
    \item $H_r$ is the right hand matrix height.
    \item $H_r$ is the right hand matrix width.
\end{itemize}
This strategy is valid only if matrices can be divided in equally sized partitions,
this is not always possible; for this reason they are padded with zeros.  In fact
every cell of a block matrix is copied in a brand new $K\times K$ matrix that is
filled with zeros.\\ Every block of the resulting matrix $C=A\times B$ can be
computed with the following formula:
\[
    C_{ij}=\sum_k A_{ik}*B_{kj}
\]
for instance, given $K=2$:

\image{img/matmul.png}{Multithread product}{0.9}

In order to perform a parallel block-multiplication, threads are required. There
are several approaches that can be used:
\begin{itemize}
    \item Simple threads
    \item Promise/Future semantics
    \item Async function
    \item Packaged task
\end{itemize}
The full usage of threads is not required, so the \textit{Promise/Future} semantics
can be used in order to add an abstraction layer to thread management. An improvement
to this approach is the \verb!std::async! function. It reduces further more the
capabilities of threads, but simplifies the \textit{Promise/Future} semantics, since
the creation of a \textit{Promise} is not required. In this case the function to run
is passed as functor or as a lambda-abstraction to \verb!std::async! this will result
in a more clear code.
\begin{lstlisting}[language= C++]
std::async(std::launch::async, foo, x, y)
\end{lstlisting}
Given an arbitrary size $K$, when a multiplication between two matrices is performed,
every block $A_{ik}$ and $B_{kj}$ is copied to a plain \verb!matrix<T>! object. The
standard matrix product is performed between \verb!matrix_wrap<T>! objects (plain
matrices are wrapped in a \verb!matrix_wrap<T>! list) that have an high polymorphic
overhead because they need to pass through the \textit{pimpl}.\\ A repeated and
alternated access to two possibly huge matrices leads to a continuous access to the
memory; the performance is very low if the access need to pass through a stack of
classes. The process can be sped up by copying the matrices to two plain
\verb!matrix<T>! objects:
\newpage
\begin{lstlisting}[language= C++]
void multiply(const matrix_wrap<T> &lhs, const matrix_wrap<T> &rhs, matrix<T> &result, unsigned start_row,
              unsigned start_column) {

    ...

    matrix<T> block(K, K);
    matrix<T> Aij(K, K);
    matrix<T> Bji(K, K);
    for (unsigned j = 0; j < col_blocks; ++j) {
        matrix<T> Aij =
            lhs.get_block({start_row,start_row + K,curr_col,curr_col + K});
        matrix<T> Bji =
            rhs.get_block({curr_col,curr_col + K,start_column,start_column + K});
        curr_col = curr_col + K;
        matrix<T> prod = simple_prod(Aij , Bji);
        block = block + prod;
    }

    unsigned max_height = std::min(start_row + K, result.get_height());
    unsigned max_width = std::min(start_column + K, result.get_width());
    for (unsigned i = start_row; i < max_height; ++i) {
        for (unsigned j = start_column; j < max_width; ++j) {
            result(i, j) = block(i - start_row, j - start_column);
        }
    }
}
\end{lstlisting}
By doing this, the overhead of the access operation due to polymorphism is reduced,
in fact the two matrices are retrieved only once and later on, only small portions of
them are accessed. A possible strategy is to use the \verb!window! function, but it
cannot be used on \verb!matrix_wrap<T>! objects because of cyclic type expansion
that leads to an infinite loop. A simpler solution is to a perform a deep copy of
the matrix portion that is required. This approach is more inefficient than an
access by reference, but is better than an access by reference in a
\verb!matrix_ref<T>! object.
\newpage
\begin{lstlisting}[language= C++]
matrix<T> get_block(window_spec spec) const override {
    matrix<T> block(K, K);
    for (unsigned i = spec.row_start; i < spec.row_end; ++i) {
            for (unsigned j = spec.col_start; j < spec.col_end; ++j)
                    block(i - spec.row_start, j - spec.col_start) = mat(i, j);
    }
    return block;
}
\end{lstlisting}
When only the \textit{pimpl} is visible, the \verb!get_block! operation is
delegated to it in order to avoid an access through \verb!matrix_wrap<T>!:
\begin{lstlisting}[language=C++]
matrix<T> get_block(window_spec spec) const {
    unsigned max_height = std::min(spec.row_start + K, height);
    unsigned max_width = std::min(spec.col_start + K, width);

    return pimpl->get_block({spec.row_start,max_height,spec.col_start,max_width});
}
\end{lstlisting}
In order to reduce the overhead more and more, the product between two blocks does
not rely on the multiplication infrastructure built for the library usage
(the lazy-like one). In that implementation, matrices are wrapped in
\verb!matrix_wrap<T>! objects and are inserted in a list, to avoid this, a simple
matrix product is used. After this computation, the block is copied to the output
matrix. Every product between blocks is performed in a dedicated thread.
\begin{lstlisting}[language=C++]
matrix<T> do_multiply(matrix_wrap<T> lhs, matrix_wrap<T> rhs) {
    ...
    matrix<T> result(lhs.get_height(), rhs.get_width());

    std::vector<std::future<bool>> workers;
    ...

    for (unsigned i = 0; i < row_blocks; ++i) {
        for (unsigned j = 0; j < col_blocks; ++j) {
            workers.emplace_back(std::async(std::launch::async, product, lhs, rhs, std::ref(result), i * K, j * K));
        }
    }

    for (auto &worker: workers)
        worker.get();

    return result;
};
\end{lstlisting}
Generally, when more threads access concurrently to the same object, it is necessary
to introduce some synchronisation mechanisms like \textit{mutexes} and
\textit{monitors}. In this case, every thread accesses to the original matrices
only for reading, whereas even if the output matrix is accessed for writing by several
threads at the same time, every thread write to a portion of the matrix that is
needed only by itself, so the behaviour is closer to parallelism than to concurrency.
We can consider the access to the matrix as non-atomic.
\par\medskip\noindent
Even if $K$ is arbitrary, its value can affect the overall performance of the
program. For instance, let $K=5$, it can be considered as a small block size if
$n,m\approx 30$, in this case matrices are split in a lot of partitions, several
threads run in parallel; if the matrix is smaller than 5 then only a small set of
cells is wasted by padded zeros. In the case of matrices that are greater than a
multiple of five, another block has to be added, but the number of padded cells is
at most $5\times \max(n,m)$.\\ On the other hand, choosing an high $K$ implies a
lower number threads involved in the computation, but the load of every thread will
be heavier. If the matrix is small, several cells are wasted by padded
zeros. In the case of matrices that are greater than a multiple of $K$, a great amount
of memory will likely be wasted by padding.\\ In order to get the best performance
from the library, a profiling should be performed trying to tune the block size $K$.
The worst cases is so described:
\[
    n,m \ll K
\]
\[
    n,m\in I^{+}(K,\varepsilon),\qquad K,\varepsilon\in\mathbb{N}
\]
where $I^{+}(K,\varepsilon)$ is defined as a neighborhood of $K$:
\[
    I^{+}(K,\varepsilon)=\{k\in\mathbb{N}\quad \textit{s.t.}\quad K<x<K+\varepsilon\}
\]
The worst case is given when $n = m = K+1$, the reason is that we create 4 matrices in which only one is dense, and the other 3 have only one column or 1 row with non-zero values. The following picture can explain better this case:
\image{img/InefficientBlock.png}{Worst case for K}{0.7}

\par\medskip\noindent
So $K$ should be tuned in order to use those matrices as little as possible. A possible
improvement of the library in terms of performance could consist on an algorithm that
every time that a product is performed, chooses the block size depending on the size of
the two matrices, selecting the value that allow to reduce as much as possible the
number of padded zeros.

%%%%%%

\newpage
\section{Assumptions and Conclusion}
Our extension essentially satisfies the requests for this assignment and introduce a structure that can be easily extended introducing new possible combinations. In fact parallel sum was not necessary for this assumption but the idea used for product was essentially the same so we decide to take advantage of this similarity including also parallelism for sum operator. More complex expression can be executed and performed as much as possible with multiple thread, for we don't provide all the expression, here we show few examples:

$$E * (F * E)\qquad (A+B)^T$$
On the first case, on the left, is is necessary to define a new operator that get in input a matrix on the left and an optimizer on the right. This is no provided from the previous version because this case is given by the usage of parenthesis inside expressions that introduce different priorities. The second case, on the right, as we have seen is a simple transformation a \verb|future_matrix|, for the moment is not available but it is necessary only to override the original methods performing first the \verb|eval| method and later return the matrix transformed. These operations are not required for this assignment but the entire structure allows to do that assuming that the matrix will be completed in the future. The constraints of this structure are: re-definition of the same structure for the \verb|future_matrix|, to distinguish the static or dynamic matrices, given by decision of avoiding possible wraps; large set of possible combinations for operations, each case has be implemented separatly also if they have more or less the same code and behavior. \\
In order to give a complete library could be interesting to introduce also this operations and possibly other ones. And as we have seen another interesting improvement could be to develop an dynamic algorithm that adjust the value of $K$ for the block parallel product between matrices. We have seen that the values of $K$ have an important role in term of performance, for that reason a possible heuristic or strategy could be reasonable to adapt.
\end{document}


text:
\section{Modules}
\subsection{Tensor product}
\begin{Definition}
    Let $R$ be a ring and $M, N$ be $R$-modules.
    A \emph{tensor product of $M$ and $N$} is a pair $(T, \phi)$, where
    \begin{itemize}
        \item $T$ is a $R$-module and
        \item $\phi: M \oplus N \to T$ is $R$-bilinear,
    \end{itemize}
    and the following universal property is fulfilled:

    If $L$ is a $R$-module and $\chi: M \oplus N \to T$ is $R$-bilinear, there is a unique $R$-homomorphism
    $g: T \to L$ such that $\chi = g \circ \phi$.
\end{Definition}

\begin{Lemma}
    Let $R$ be a ring and $M, N$ be $R$-modules.
    \begin{enumerate}
        \item The tensor product of $M$ and $N$ is unique up to isomorphism.
        \item The tensor product of $M$ and $N$ exists.
    \end{enumerate}
\end{Lemma}

max_stars_repo_name: lugggas/math

max_issues_repo_name: lugggas/math

max_forks_repo_name: lugggas/math


text:
%!TEX root = ../thesis.tex
% ******************************* Thesis Appendix A ****************************
\chapter{Detailed breakdown of systematic uncertainties}
\label{app:uncert}

\begin{longtable}{ l | c | c | c }%[htbp!]
%\center
\caption{Extracted cross sections with detailed
  list of uncertainties. Besides the contribution to the total
  uncertainty in \%, the fitted value of the nuisance parameter (pull), as well as the ratio of the estimated
  uncertainty over the uncertainty from a 1 $\sigma$ variation, called
  constr/$\sigma$, are shown.
  \label{tab:lh_res_eightfull}}\\
\hline
Name & Pull & Constr / $\sigma$ & Contribution [\%] \\
\hline
\endfirsthead
\hline
Name & Pull & Constr / $\sigma$ & Contribution [\%] \\
\hline
\endhead
\input{Results/Figures/FitPlots/xsecFit_tab13TeV.tex}
\end{longtable}

max_stars_repo_name: TillArndt/PhD-Thesis

max_issues_repo_name: TillArndt/PhD-Thesis

max_forks_repo_name: TillArndt/PhD-Thesis


text:
\chapter{The Story of Guṇanidhi}

Sūta said:

1. O great sages, after hearing these words of Brahmā, Nārada once again bowed
to him and asked humbly.

Nārada said:

2-3. When did Śiva favourably disposed to His devotees go to Kailāsa? Where did
he have the intimate acquaintance with Kubera\footnote{Kubera is the son of
“viśravas” by “iḍaviḍā”. He is the chief of the “yakṣas” and a friend of Rudra.
He is mythologised as having three legs and eight teeth.} of great and noble
soul? What did Śiva of auspicious form do there? Please narrate all these things
to me. I am deeply interested in it.

Brahmā said:

4. O Nārada, listen. I shall tell you the story of the moon-crested lord, how he
went to Kailāsa and how he contracted friendship of Kubera.

5. In the city of Kāmpilya\footnote{The country known to \tl{Vājasaneyi Saṃhitā}
(xxiii, 18) and \tl{Śatapatha Brāhmaṇa} (xiii. 2.8.3) can be identified with
the city of Kāmpīla in the Furrukhabad district, Uttara Pradeśa. It was
the Southern Capital of Pāñcāladeśa in ancient India. Dr. Awasthi however,
places it in the Ānartadeśa, a region of the Western India.} there was
a sacrificer named Yajñadatta. Born of Somayāji family he was an adept in
the performance of sacrifice.

6. He knew Vedas and Vedāṅgas. He was a great scholar of Vedānta \etc He was
honoured by the king. He was a liberal-minded donor and as such his fame had
spread far and wide.

7-8. He assiduously maintained the sacrificial fire and was devoted to the study
of the Vedas. His son (Guṇanidhi) was of a very handsome complexion and shone
like the moon’s disc. After the investiture with the sacred thread he learned
all the eight lores\footnote{The eight sciences included (1) the triple Veda
(trayī), (2) logic and metaphysics (ānvīkṣikī), (3) the science of Government
(daṇḍanīti), (4) practical arts such as agriculture, commerce, medicine \etc
(vārtā), (5) ancient historical and mythological tradition, (6) science of
rituals, (7) logic and (8) dharma or law.} over and over again. Yet, unknown to
his father he indulged in gambling.

9. Ever and anon he took plenty of sums from his mother and gave them over to
other gamblers with whom he contracted great intimacy.

10. He eschewed all brahminical ways and conduct of life. He was averse to
the performance of \tl{sandhyā} prayers and ceremonial ablutions. He began to
speak ill of the Vedas, sacred texts, \tl{devas} and \tl{brāhmaṇas}.

11. He did not follow the conventions and injunctions of the Smṛti code, He
indulged in singing and playing. Actors, heretics \etc were his beloved friends.

12-15. Although his mother wanted him to meet his father now and then, he never
went near his father. Engaged in extra-domestic activities Yajñadatta used to
ask his wife “Dear good woman, what is our son Guṇanidhi doing? He is not at
home.” Then the woman used to say, “He has gone out just now. So long he had
been taking his bath and worshipping the deities. He has finished his Vedic
studies and has just gone out in the company of two or three friends for
the purpose of learning somewhere”. The poor woman in view of the fact that she
had only one son deceived her husband thus.

16. The simple husband did not know anything about the nefarious activities of
his son or his bad conduct. All sacred rites ending with \tl{keśakarma}\footnote{
The religious ceremony \tl{keśānta} in which the hair were cut off was performed
upon \tl{brāhmaṇas} at 16 years of age, \tl{kṣattriyas} at 22 and \tl{vaiśyas}
at 24. Cf. \tl{Manu. U.} 65, \tl{Yājñavalkya} 1.36.} were performed in
the sixteenth year of the son.

17. Thereafter Yajñadatta performed the marriage rite of the son in accordance
with the rules prescribed in the \tl{Gṛhya Sūtras}.

18. O Nārada, every day the woman with her heart melting with motherly affection
used to make her son sit up and gently upbraid him.

19. “Dear son, your father is surely a great man, but he is of rash temperament.
If he comes to know of your activities he will beat you and will not spare me
too.

20. I conceal your nefarious activities from your father every day. Due to his
good conduct and his affluent circumstances he is honoured by all the people.

21. Dear child, a good learning and association with men of saintly character
constitute a great asset for \tl{brāhmaṇas}. How is it that you do not gladly
take interest in such things?

22. Your ancestors and grandfathers had all earned the reputation of being good
Vedic scholars, well learned in Śāstras, and performers of sacrifices,
especially \tl{somayāgas}.

23. Shun the company of the wicked people, associate with good men, turn your
attention to good learning and strictly adhere to brāhmaṇical conventions.

24. Emulate your father in form, fame and traditional activity. Why don’t you
feel ashamed? Cast off your wickedness.

25. You are nineteen now. This girl is sixteen years old. She is a good girl.
Take her. Protect her. Above all be devoted to your father.

26. You shall respect your father-in-law also, in view of his good qualities and
conduct. How is it that you do not feel ashamed of wickedness?

27. Dear son, your maternal uncles too are matchless in learning, conduct and
pedigree and other things. You are not afraid even of them. Your paternal and
maternal lineages are equally pure.

28. See the \tl{brāhmaṇas} boys of your neighbourhood. Even in our house see
the disciples of your father. How humbly do they behave?

29. Dear son, if the king hears of your evil propensities, he will cease to
respect your father and may even suspend the regular maintenance allowance.

30. Till now people used to call your activities the foolish blunders of an
ignorant boy. Hereafter they may take away the traditional title of \tl{dīkṣita}.

31. People will curse and cavil at your father and me saying such evil words as
“The son has adopted the wickedness of the mother.”

32. Your father has never been a sinner. He strictly follows the path of
the Vedas and Smṛtis. Lord Śiva is my witness for the purity of my mind that is
riveted to his feet.

33. I have not seen the face of any wicked man after my menstrual bath. Powerful
indeed is Fate whence a boy like you is born of my womb!”

34. Although constantly advised thus by his mother, the wicked boy did not
abandon his evil ways. For, an idiot indulging in vice is beyond redemption.

35. Who is he that is not broken up by the evil influences of hunting, wine,
slander, untruth, theft, gambling and prostitutes?

36. The wicked fellow (Guṇanidhi) used to lay his hands on whatever he could see
in the house, a cloth, a base metal \etc and take it to the gambling den, there
to lose the same to his brother gamblers.

37. Once he stole a very valuable ring of his father set with precious stones
and gave it to one the of gamblers.

38. It chanced that one day the \tl{dīkṣita} saw it in the hand of the gambler.
He asked the fellow—“Where did you get this ring from?”

39-40. First the gambler did not say anything. When repeatedly asked he said—“O
\tl{brāhmaṇa}, you are unnecessarily accusing me of theft. It was your son who
gave it to me. On the previous day I had won his mother’s upper garment.

41. Do not think that I alone was the winner of this ring. He has lost many
costly things to other gamblers as well.

42. He has thus given gems, metals, silk garments, vessels, golden vases, and
different sorts of copper and bell metal pots.

43. Everyday he is being bound stark naked by the gamblers. In the whole world
you cannot see such a useless poor gambler as he (your son).

44. How is it that till now, O \tl{brāhmaṇa}, you have not realised that your
son is a ring leader of base gamblers, very clever in misdemeanour and unfair
means?

45. On hearing these words, the poor \tl{dīkṣita}’s head bent down with shame.
He covered his face and head with a cloth and quietly slipped back into his
house.

46. Yajñadatta, the sacrificer, well versed in Vedic rites spoke thus to his
wife who was a very chaste lady.

Yajñadatta said:

47-48. O mistress! where is that gambling rogue of a son, Guṇanidhi? Or let it
be. Why should I ask for him? Where is that auspicious ring which you took off
at the time of applying unguents on my body? Bring it quickly and give it to me.

49-51. The mistress was frightened at these words. While she was engaged in
arranging for bath and midday sacred rites she replied—“O lord, I am busy
arranging the various articles of offerings for worship. O lord, fond of guests,
the guests may be unnecessarily detained. While I was busy cooking the pudding
I kept the ring somewhere in some vessel just now. What a pity! I have forgotten
it. I do not know where it has been kept.

Dīkṣita said:

52-53. O truthful lady who has given birth to a base boy, whenever I asked
“Where has the son gone?” you used to say, “Dear lord, just now he has gone out
after finishing his lesson of the Vedas, in the company of two or three friends
for revision of the lesson”.

54. Where is your silk saree red like madder which I had presented to you and
which used to hang down here in the house always? Tell me the truth. Do not be
afraid.

55. That gem-set golden vase which I had given you is also missing. That tripod
with a velvet cushion which I had given you is nowhere to be seen.

56. Where is that bell metal pot made in the South? Where is that copper pot
made in Bengal? Where is that ivory casket intended for curios and trinkets?

57. Where is that wonderfully fine statuette of a lady lighting a lamp, shining
like the moon, and brought from the hilly province?

58. Why should I unnecessarily speak much? O lady of a noble family, it is
futile to be angry with you. I shall take food—only after I marry again!

59. I am childless now since that wicked fellow has defiled the whole family.
Get up and fetch me some water. Let me offer libations to him with gingelly
seeds\footnote{It is customary among the orthodox Hindu families in India to
offer libations of water mixed with gingelly seeds to the manes on particular
days.}.

60. Better to be issueless than have a wicked son who defiles the entire family.
It is the traditional policy to abandon one to save the family.

61. The \tl{brāhmaṇa} took his bath, performed his daily rites and married
the daughter of a Vedic scholar the same day.

max_stars_repo_name: elvendrim/shivamahapurana

max_issues_repo_name: elvendrim/shivamahapurana

max_forks_repo_name: elvendrim/shivamahapurana


text:
\section{Problem Statement}

The hereby \textbf{Report 5} will state an essay for literature knowledge that might support our research. We evaluate several research work. This report will focus on the topic towards (i) \textbf{Gamification \& Crowdsourcing} application on \textit{Health Informatics}; (ii) \textbf{Gathering Annotated Data}; and (iii) \textbf{Interactive Learning}. All topics are of chief importance to our research work, therefore should be analysed.

max_stars_repo_name: mida-project/reading-reports

max_issues_repo_name: mida-project/reading-reports

max_forks_repo_name: mida-project/reading-reports