text:
\documentclass[12pt,a4paper]{article}
% AUTHOR: Rafael Belchior
% Thanks to Prof. RUI SANTOS CRUZ for providing the template
%
\usepackage{helvet}
\renewcommand{\familydefault}{\sfdefault}
\usepackage{a4wide}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SELECT ONE OF THE FOLLOWING PACKAGES FOR THE LANGUAGE
\usepackage[english]{babel}
% \usepackage[portuges]{babel}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{subfig}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{cite}
\usepackage[absolute]{textpos}
\usepackage{tabularx}
\usepackage{tabulary}
\usepackage{fancyhdr}
\usepackage[table]{xcolor}
\pagestyle{fancy}
\headsep=50pt
\setlength{\headheight}{50pt}
\usepackage{listings}
\usepackage{minted}
\definecolor{LightGray}{rgb}{0.95, 0.95, 0.95}
\definecolor{darkblue}{rgb}{0.0,0.0,0.6}
\definecolor{editorOcher}{rgb}{1, 0.5, 0}

% Clever Referencing of document parts
\usepackage{cleveref}

\lstdefinestyle{commandline} {%
language={[WinXP]command.com},
breaklines=true,
%aboveskip=\baselineskip,
belowskip=\baselineskip,
showstringspaces=false,
backgroundcolor=\color{LightGray},
basicstyle=\small\color{black}\ttfamily,
showstringspaces=false,
keywordstyle=\color{cyan}\bfseries,
stringstyle=\color{cyan}\ttfamily,
commentstyle=\color{green}\itshape,
moredelim=[s][\color{blue}\bfseries]{C:}{\>}
}

\lstdefinestyle{Bash} {%
language=bash,
breaklines=true,
belowskip=\baselineskip,
backgroundcolor=\color{LightGray},
showstringspaces=false,
keywordstyle=\color{black}\bfseries,
basicstyle=\small\color{black}\ttfamily,
stringstyle=\color{editorOcher}\ttfamily,
commentstyle=\color{cyan}\itshape,
otherkeywords={xcode-select, mkdir,rm},
moredelim=[s][\color{red}]{~$},
literate={~} {$\sim$}{1}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLEASE FILL THE ADEQUATE DATA IN THE TABLE REPLACING
% THE VALUES EXEMPLIFIED
\lhead{}
{\renewcommand{\arraystretch}{1.1}
\fancyhead[C]{\begin{tabularx}{1.0\textwidth}{|l|X|l|l|}
\hline
% In the following line change Course Name: PPIII, PPB
\textbf{EB 20/21} & \textbf{Enterprise Blockchain Technologies} & \textbf{Number:}  &  1 \\
\hline
% In the following line insert your Name and IST ID
\multicolumn{2}{|l|}{Module I - Introduction} & \textbf{Issue Date:}  &  14 Sept 2020 \\
\hline
% In the following line insert the Activity CODE and Title (abridged)
%\textbf{WP n.} (99) & (Subject) & \textbf{Group:} & (99) \\
\multicolumn{2}{|l|}{Background: Distributed Systems} & \textbf{Due Date:} &  21 Sept 2020\\
\hline
\end{tabularx}}
\rhead{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DO NOT CHANGE THIS BLOCK
\begin{document}
\textblockorigin{-34pt}{-12pt}
\begin{textblock*}{10cm}(2cm,1cm)
\includegraphics[width=6cm]{hyperledger.png}
\end{textblock*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%,sdist2017

\section*{Preliminary Notes}
This class recalls background on distributed systems and crpytography, essencial building blocks for studying blockchain technology. This laboratory is based on several sources \cite{dscd,sdist2017,Wong2014,Ousterhout_presentation,Verissimo2001,correia2019byzantine,raft_paper}. We recommend students to conduct additional research on other sources \cite{mit,princeton}.

%The reference section should be viewed as a ``additional readings reference'' - if you would like more information about the topic.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR TEXT STARTS HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Distributed Systems}
\label{sec:ds}
A Distributed System is a system comprised of software and hardware components, that are connected on a network, coordinating their actions via message exchange.

\begin{figure}[h!]
    \centering
    \includegraphics[scale=0.25]{figures/itsystem.png}
    \caption{Enterprise IT system \cite{Wong2014}}
    \label{fig:internet}
\end{figure}

Nowadays, most modern enterprise information systems are distributed systems, popularized due to the growth of the Internet.

Servers communicate with other servers using protocols, such as the HTTP protocol, based on the TCP/IP protocol. Communication can also happen locally (i.e., between computer processes, using interfaces such as sockets). Typically, such protocols are independent of the system architecture.

In fact, there are different architectures for organizing distributed systems. One of the simplest is the client-server architecture, where servers can expose endpoints to clients.
In this architecture, servers keep resources and serve those resources according to requests made by clients. Client-server architectures have a popularized programming model: RPC, and the more recent gRPC. Such protocols allow to structure distributed system programming, based on client calls on server's code. While this is an efficient scheme, it constitutes a centralized paradigm, that suffers from issues such as lack of scalability, and constitutes a single point of failure. Other architectures offer different tradeoffs: for instance, the peer-to-peer architecture is decentralized, which allows for a higher degree of openness, and fault tolerance. This makes the whole system more robust to faults.

\subsection{Faults, Errors, and Failures}
A fault is an event that alters the normal behaviour of a system, which typically causes errors. Errors comprise the transition from a correct state of the system to an incorrect state. Errors can provoke failures, if the system deviates from its specification. The protection against faults allows participants to remain consistent: this way it is easier to maintain a correct state, that can be shared with other participants. Faults can be defined as crash faults or byzantine faults. Crash faults are faults where a node becomes unresponsive. This can be due to either physical problems, like a power outage or a faulty component, or due to networking problem, like network partitions which lead to the node being unreachable. Byzantine faults occur when a node acts arbitrarily and strays from the specification. This can include fake information, different responses to different nodes or even a completely random answer.

\subsection{Fault Tolerant Models}
There are different ways to deal with each type of faults.
Crash Fault Tolerant (CFT) models can usually withstand up to $\frac{N}{2}$ crashes, with N being the total number of nodes. This model assumes a quorum of $\frac{n}{2}+1$ nodes which has to agree on the global state. This means that as long as there is a \textbf{majority of nodes operational}, a consistent answer can be returned.

However, if one of the nodes is malicious, the system may enter a faulty (or inconsistent) state with the node giving fake information. As such, for Byzantine Fault Tolerant (BFT) models, we need at least $3f+1$ correct nodes to support $f$ malicious nodes. This way, we can guarantee that in any given quorum, we have at least a \textbf{majority of correct nodes} answers.

%figure client server vs distributed

But a question may arise: how to make sure that, in decentralized networks, the participants agree on a single state? How can we assure agreement under uncertain conditions, and assuming that participants can fail (e.g., connection problems)? What if there are malicious participants? This problem is also referred as the consensus problem, a well known problem studied the early stages of distributed systems.

\subsection{Consensus}
Consensus is a fundamental problem in distributed systems, in the area of fault tolerance. Consensus involves multiple servers (or nodes) agreeing on a global state (set of values). %The problem of consensus has been studied since the late 1970s, .
The problem of reaching consensus was popularized by Pease, Shostak, and Lamport's Byzantine generals story, where the generals had to agree on a common plan to attack a city, where the messengers that exchanged their messages where unreliable \cite{byzantine_generals}.

Distributed systems ought to be \emph{Byzantine fault-tolerant}, as there may be malicious nodes on the network \cite{correia2019byzantine}.  A consensus algorithm is used to create agreement on a global ledger state in the presence of crash faults or Byzantine faults. Consensus algorithms depend on the assumptions made about the environment and the system that the algorithm is running. For instance, what are the faults that can happen (crash vs byzantine), how are messages passed, synchronous (communication and processing delays are bounded) vs asynchronous (communication and processing delays are not bounded, i.e., may take an arbitrary amount of time to be completed) environment.

\subsection{State Machine Replication}
A general approach to building fault-tolerant systems is state machine replication, that is studied in the context of consensus. A state machine can be described by a set of possible initial states and next relation. This relation determines the possible steps, given the current step.

In a distributed system, each participating node has a state machine and a log. A log is a collection of log entries, which typically contain timestamped information. The state machine is a device (computer) that stores and manipulates state (e.g., values written in an hash table, the value of a variable, an array with commands to be executed).

The idea of state machine replication is that clients can interact with a distributed system that agrees upon a global state, while it appears that it is interacting with a single state machine. Even if one of the server fails, the system still responds to the requests. Each node takes commands from the logs, triggering the execution of a function. The consensus ensures that the commands received are the same while ensuring a common ordering for the commands received. This ensures that each state machine yields the same series of results, achieving the same series of states, at a common order. Each state machine executes the same function calls, in the same order, yielding the same output. For this reason, state machines need to be deterministic (otherwise, the same input could yield different outputs, causing inconsistencies in the system).

 \begin{figure}[h!]
    \centering
    \includegraphics[scale=0.25]{figures/rsm.png}
    \caption{Replicated state machines \cite{replicated}}
    \label{fig:rms}
\end{figure}

Figure \ref{fig:rms} shows an example of a replicated state machine. A group of machines executes a deterministic function based on a replicated log they receive (input). The idea is that logs are kept identical so the execution on different machines is performed using the same commands, and on the same order. When the first machine receives an update request (set variable Z to 6), it first replicates that command into all the logs, passing that command to the other machines on the cluster. All machines add the newest command to the end of their logs. After that, the command can be executed by each state machine. This triggers the original state machine to return the answer of the execution to the client. If one of the servers (that comprises a state machine) fails, the other servers can still process the request, as far as a majority of servers are operational. There are consensus algorithms that apply either Crash Fault Tolerant or Byzantine Fault Tolerant models. These systems allow to develop and deliver resilient, dependable systems that can satisfy today's needs in terms of availability and scalability.

With this in mind, consensus algorithms then aim to replicate commands into the replicated logs, allowing us to assure desirable properties about a system, like validity, agreement and termination.

A system is valid if the value decided by the system was a value proposed by a node. This way, we avoid cases where a node can inject a decision that no one has seen. Moreover, a system is in agreement if no two correct processes decide differently (e.g., output different values for the same input). Termination is a liveness property, that states ``every correct process eventually decide``. Those are all desirable properties for fault-tolerant systems.

In this class, you will work in greater detail the RAFT algorithm: a well-adopted consensus algorithm, which focuses on understandability, and has many practical applications, inclusively on the blockchain area.

\subsection{The RAFT consensus algorithm}

RAFT is a crash fault tolerant consensus algorithm that is designed to be easy to understand.  %Typically, crash fault tolerant leader-based algorithms have a linear message complexity, O(n), as opposed to Byzantine leader-based consensus, which incurrs in a quadratic message complexity (On2).

It belongs to a class of consensus algorithms called leader-based (or coordinator-based or primary-backup). In leader-based consensus, there is a leader that calculates a decision on the state and communicates it to the followers; if the leader is faulty, a new leader is elected.
RAFT considers several participants:
\begin{itemize}
    \item Leader: the node that receives client requests, manages the log replication, and sends log update requests to other nodes. While having a leader that commands the global state is sufficient fine in some scenarios, a malicious leader can constitute a problem, as eventual inconsistencies are corrected (overwritten) by the leader.
    \item  Followers: constituting all nodes except the leader;
    \item Candidates: followers requesting to be the leader.
\end{itemize}

As explained before, a distributed system built with RAFT has several servers running an instance of RAFT and a persistent log to store the commands received from clients. The replication of logs and consequent command execution is given by three phases:
\begin{enumerate}
    \item Leader Election
    \item Log Replication
    \item Safety
\end{enumerate}

Figure \ref{fig:state_empty} shows the global state, which corresponds to the state of each node (server), across the different terms (time periods in which the algorithm executes). Terms are used to detect obsolete information: there is at most one leader per term, and followers do not accept information coming from a previous term than the one they record. This way, the algorithm guarantees that each node has the latest information.

\begin{figure}
    \centering
    \subfloat[Initial state of the shared log of operations]{    \label{fig:state_empty}{\includegraphics[width=7cm]{figures/RAFT_state_empty.png} }}%
    \qquad
    \subfloat[Leader S2 received a command from a client. It wishes to updated state]{    \label{fig:RAFT5}{\includegraphics[width=7cm]{figures/RAFT_state.png} }}%
    \caption{States at different terms for each server}%
    \label{fig:state_empty}%
\end{figure}

\subsubsection{Leader Election}
Time in RAFT is divided into terms. Each terms has a unique identifier, a monotonically increasing number, which is set to 0 at server startup time. RAFT uses randomized timers to elect the leader, illustrated by Figure  \ref{fig:le_1}. In order to know if the leader is alive, RAFT uses an heartbeat system. At the beginning of the consensus process, nodes are awaiting for heartbeats from the leader (Figure \ref{fig:RAFT1}). When a follower receives a heartbeat, their timeouts are reset. As there is no leader, no heartbeats are received within the election timeout, and eventually a timeout is triggered. The node that triggered the timeout sends a proposal to become the new leader to every other node (Figure \ref{fig:RAFT2}), and eventually becomes the new leader.

If a node does not receive a heartbeat from a leader within a certain time frame (for example, in case the leader crashes), there is another timeout that initiates the beginning of a new term and a new leader election process starts.

\begin{figure}
    \centering
    \subfloat[There is no leader, each node awaits instructions]{    \label{fig:RAFT1}{\includegraphics[width=7cm]{figures/RAFT1.png} }}%
    \qquad
    \subfloat[After the timeout, the node sends a proposal to become the new leader]{    \label{fig:RAFT2}{\includegraphics[width=7cm]{figures/RAFT2.png} }}%
    \caption{Leader election process}%
    \label{fig:le_1}%
\end{figure}

The terms also allows for liveness: nodes reject requests coming with an invalid term number (older than the current term). Nodes ask to be a leader in a  first-come-first-served basis. As the majority rule applies, only one node can win an election and become the leader. In case of a draw, another round of election proceeds. Note that at the beginning of the execution, the log is empty (Figure \ref{fig:state_empty}).

\subsubsection{Log Replication}
The idea of log replication is the leader to accept requests from clients, which translate to commands that the server has to execute. Those commands are appended to a log. After that, the leader replicates its logs to other nodes on the network, so they can execute the received commands as well, in order to keep the states of all nodes consistent. Part of this procedure is illustrated by Figure \ref{fig:su_2}.

Let us now suppose that a leader receives a request from a client (Figure \ref{fig:RAFT_state}). The log is updated, but it is still not safe to execute the received commands (as the log entry is considered uncommitted). This happens because followers are not aware of the request (illustrated by the dashed lines around an update). The leader broadcasts the most recent log entry to their followers, and receives a message from each one acknowledging that they received the updated log entry (Figure \ref{fig:RAFT6}). After the leader confirms the state, after receiving the acknowledgments from the followers, the followers can append the command to their logs, for later executing the command.

\begin{figure}[h!]
    \centering
    \subfloat[Followers acknowledge the latest log entry]{    \label{fig:RAFT_state}{\includegraphics[width=7cm]{figures/RAFT5.png} }}%
    \qquad
    \subfloat[Global status is committed, that is, the received command from the client is present in the logs of the majority of nodes, so each node can now execute it]{    \label{fig:RAFT6}{\includegraphics[width=7cm]{figures/RAFT6.png} }}%
    \caption{State updates}%
    \label{fig:su_2}%
\end{figure}

Imagine that now, node S2 crashes due to a hardware fail. This causes that S2 cannot send heartbeats back to each follower. We now assume that S1 timeouts - it becomes a candidate and tries to become a leader (Figure \ref{fig:crash}). When a follower becomes a candidate, it increments the current term, votes for itself, and send requests to other followers, to become the new leader. Eventually, S1 wins the election and becomes the new leader (Figure \ref{fig:new_leader}).

\begin{figure}[h!]
    \centering
    \subfloat[Server S2 crashes. A follower triggers a timeout, becoming a candidate. The candidate S1 tries to become the leader]{    \label{fig:crash}{\includegraphics[width=7cm]{figures/s2_crash.png} }}%
    \qquad
    \subfloat[Candidate S1 becomes the leader]{    \label{fig:new_leader}{\includegraphics[width=7cm]{figures/new_leader.PNG} }}%
    \caption{New leader election process}%
    \label{fig:leader_election}%
\end{figure}

\subsubsection{Safety}
Logs need to be kept consistent, and, furthermore, that only nodes with the latest log version can become leaders.
If a crash occurs, log inconsistencies may happen. In other words, the information stored by each node may not be the same at the same point. To solve this problem, the log of the current leader is always considered the correct version. The leader repairs inconsistencies, if applicable. For that, there cannot be missing log entries, and log entries need to be ordered.

RAFT's approach to maintain consistency is by matching logs and solving conflicts: each time a leader sends a log entry to a follower, it includes the term number and log index of the preceding entry. If there is a match, the follower accepts the new log entry; otherwise, the follower rejects the log entry. In case a log entry is rejected, the leader tries again, sending the current log entry, the preceding one, and the one further before., until a match happens. When a match happens, the follower's log is updated, accordingly to the leader.

The explained mechanisms guarantee the log matching property: if a given entry is committed, all previous entries are committed and, if two logs share the same log index and term, then those logs are identical (represent the same log entry). Lastly, the desirable safety tackles leader completeness: once a log entry is committed, the leader will have that log entry stored. This assures a committed entry is never lost. Nodes with incomplete logs are not chosen in the election process. Essentially, a node denies the vote if their log is more complete than the candidate's log.

\section{Hands on RAFT}
Now, that we illustrated the RAFT algorithm, let's consolidate that knowledge. It is recommended for students to read the additional support materials available: \cite{raft_paper,raft_homepage}. You may use the simulator provided at RAFT's homepage to aid you in your responses\footnote{https://RAFT.github.io/}.

\subsubsection*{Exercise 1:} To consolidate your knowledge, watch the guided visualization of the RAFT algorithm \cite{raft_viz}.
Elaborate on how the RAFT algorithm can contribute to increase the availability of the system, from a client perspective.

\subsubsection*{Exercise 2:} When does a new term start?

\subsubsection*{Exercise 3:} How does candidates act if two of them attempt to become leaders at exactly the same time?

\subsubsection*{Exercise 4:} If a follower receives a message has a different term than the one recorded, does it accept the message?

\subsubsection*{Exercise 5:} Suppose that a leader receives a command from a client, and notifies all followers. Upon sending the second transaction, one of the followers crashes. How does RAFT guarantee that the follower is updated?

\subsubsection*{Exercise 6:} Consider the following RAFT cluster, composed by a leader and four followers. The leader has sent eight messages, corresponding to eight commands. Those eight commands translate into eight log entries. The current term is term three.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.6]{figures/exercise.png}
    \caption{RAFT cluster composed of 5 machines}
    \label{fig:RAFT_cluster}
\end{figure}

What are the commands, represented by each log index, that are safe to be executed by each node?

\subsubsection*{Exercise 7:} Consider the following RAFT cluster, composed by a leader and four followers, in which there are inconsistencies on the logs across nodes. Assume that node S4 crashed on term 1. Node S5 became the leader on term 2, but was only successful to replicate logs 4 and 5, before crashing. Node S3 is the leader on term 3.

\begin{figure}[H]
    \centering
    \includegraphics[scale=0.6]{figures/inconsistent_logs.png}
    \caption{RAFT cluster composed of 5 servers, showing inconsistencies regarding the log}
    \label{fig:inconsistent_logs}
\end{figure}

What will be the final state of the system, assuming no new commands are received, and that S5 recovers from crash?

\subsubsection*{Exercise 8:} During the leader election process, a mechanism that ensure consistency across logs is the rejection of candidates that have outdated logs. What are the checks the follower does against the candidate's proposal? In which occasion does it reject voting on the candidate?

\subsubsection*{Exercise 9:} Consider the following RAFT cluster, composed by seven nodes.

\begin{figure}[H]
    \centering
    \includegraphics[scale=0.6]{figures/discussion.png}
    \caption{RAFT cluster composed of 7 nodes}
    \label{fig:inconsistent_logs}
\end{figure}

Given that the leader of term eight has ten log entries, and interacts with six followers, how did the system get into this state?

\subsubsection*{Exercise 10:} How does the system re-organize itself if a node leaves the system?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR TEXT ENDS HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{IEEEtran}
\bibliography{lab1.bib}

\end{document}                             % The required last line

max_stars_repo_name: RafaelAPB/university-course

max_issues_repo_name: RafaelAPB/university-course

max_forks_repo_name: RafaelAPB/university-course


text:
\documentclass[a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage{parskip}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{biblatex}
\addbibresource{sources.bib}

\title{Techniques in Functional Programming\\
Term Project\\
Fold right and fold left}
\author{Benjamin Nørgaard \#201209884}
\date{April 2015}

\begin{document}

\maketitle
\tableofcontents
\clearpage
\section{Introduction}
This project explores the functions \texttt{fold\_right} and
\texttt{fold\_left}. In order to understand what the fold functions do and how
they work, we will define a few functions using folds. We will then see if we
can define one fold function in term of the other. We will find that there are
some cases where it does not make any difference whether we use
\texttt{fold\_right} or \texttt{fold\_left}, and we will try to formalise when
this applies. In the end we will draw parallels to primitive iteration and
primitive recursion.

\subsection{About this report}
This report was written to document a project written using Coq Proof Assistant.
The report is meant to be read alongside the actual Coq source file. There will
not be be code in the report, but I will instead refer to and talk about the
code written in the source file.

\subsection{About the source file}
The source file is long, but it is structured more or less like the sections in
this report, so it should be easy enough to follow along in the code.

We use a few conventions in the source code. Most noticeably is the meaning of
the indentation: Proofs are indented by two spaces per subgoal. This makes it
easier to for example see the cases of an induction proof visually in the
indentation.

The code has quite a few comments, with the purpose of getting the reader to
faster understand the intention behind the steps.

\section{Fold right and fold left}
A fold is a function which traverses a list data structure and replaces the
list constructors with some other given function. This means that the given
function, like the list constructor should take two arguments. The first should
be of the same type as the type that list holds. The other should be the same
type as the return type of the function.  The last type requirement is due to
the result of the function being fed to the next application of the function.
One should also pick a value which will be used as the initial value of the
second argument.

\texttt{fold\_left} means that the first application of the function will be
made on the left most element in the list. \texttt{fold\_right} is the opposite
where the first application of the function is made on the last element.
Figure~\ref{fig:folds} illustrates this very well.

\begin{figure}[h]
\includegraphics[width=0.5\textwidth]{fold_left}
\includegraphics[width=0.5\textwidth]{fold_right}
\caption{Illustrations of fold left and fold right.~\cite{foldl}\cite{foldr}}
\label{fig:folds}
\end{figure}

We will now look into how we could implement this in Coq.

\subsection{Fold right}
We implement \texttt{fold\_right} by matching on the input list. If the list is
empty, then we return the \texttt{nil\_case}. If the list is a head and a tail,
we apply the \texttt{cons\_case} to the head and call \texttt{fold\_right}
recursively on the tail.

\subsection{Fold left}
We implement \texttt{fold\_left} by matching input list. If the list is empty,
we return the \texttt{nil\_case}. If the list is a head and a tail, we call
\texttt{fold\_left} recursively, where the \texttt{nil\_case} is replaced by
an application of the \texttt{cons\_case} on the head and the former
\texttt{nil\_case}, and the list on which we now do the fold is the tail of the
list.

\section{A few functions}
In this section we will specify, define and implement a few different functions
using \texttt{fold\_right} and \texttt{fold\_left}. This will act as a stepping
stone, while we are learning to better understand \texttt{fold\_right} and
\texttt{fold\_left}.

One should note that all of the following functions could be implemented without
using fold. But when we have the fold mechanism it turns out to be very easy to
implement these functions.

\subsection{List Identity}
We start with the simple example of the \texttt{list\_identity} function. The
\linebreak\texttt{list\_identity} function takes a type and a list of that
type, and returns the list it is given.

We define the \texttt{list\_identity} with \texttt{fold\_right}. The idea is
that we let the \texttt{cons\_case} function be the \texttt{cons} function, and
the \texttt{nil\_case} be \texttt{nil}. In this way any input list would simply
be taken apart and be put back together by the fold.

We prove that it fits the specification with a simple structural induction
proof.

\subsection{List append}
In order to define a \texttt{reverse\_list} function, we will first need an
implementation of \texttt{list\_append}. We will also need
\texttt{list\_append} when we look when we look at the \texttt{map\_append}
function. It happens that we can implement \texttt{list\_append} with fold as
well. The \texttt{list\_append} function takes a type and two lists of that
type, and returns a list where one has been appended on the other.

We implement this almost like the \texttt{list\_identity} function, but this
time we let the \texttt{nil\_case} be the list that we want to append to the
other. This means that the first list is taken apart and constructed with the
\texttt{cons} function, when we reach \texttt{nil}, we replace it by the
\texttt{nil\_case} which is the second list. In turn this gives us exactly what
we want.

Later on in our proofs about \texttt{reverse\_list}, we will need the property
that \texttt{list\_append} is associative, in order to complete the proof. We
prove this property using structural induction on the first list.

\subsection{Reverse list}
With a specification and an implementation of \texttt{list\_append}, we are now
in position to define \texttt{reverse\_list} using fold. The
\texttt{reverse\_list} function takes a type and a list of that type and
returns the reversed list.

This function is actually almost the same as the \texttt{list\_identity}
function. The only difference is that we use \texttt{fold\_left} instead of
\texttt{fold\_right}. This means that we take apart the list and then we
construct it backwards. This in turn gives us the reversed list.

The reason why we needed \texttt{list\_append} is because the specification of
\linebreak \texttt{reverse\_list} relies on \texttt{list\_append}. In fact we
don't need an actual implementation of \texttt{list\_append} at least not in
the proofs about \texttt{reverse\_list}. What we \emph{do} need is the
specification of \texttt{list\_append}. When we have the specification, we can
just make an assumption that we have an implementation, and the use this
assumption to do the proof.

When we prove that this fold implementation of \texttt{reverse\_list} fits the
specification of \texttt{reverse\_list}, we need a master lemma in the
\texttt{cons\_case}. The master lemma is the following property about
\texttt{fold\_left}: having a non-empty list as \texttt{nil\_case} is the
equivalent to letting the \texttt{nil\_case} be \texttt{nil} and then append
the list to the result of the fold. This is proven by structural induction and
with the help of the property that \texttt{list\_append} is associative.

We will see this pattern of using a master lemma in the \texttt{cons\_case} a
few more times in this project.

\subsection{Map}
In order to further explore how folds work, we decided to go a bit beyond the
project and look into \texttt{map} and a variation of it called
\texttt{map\_append}. We start by looking at the ordinary \texttt{map}
function.

\texttt{map} takes two types, a function from the first to the second type and
a list of the first type.  It then returns a list where the function has been
applied to each element of the former list.

We implemented this using \texttt{fold\_right}. We let the \texttt{cons\_case}
of the fold apply the function to the current element, and we use \texttt{cons}
to construct the new list.

\subsection{Map append}
\texttt{map\_append} is a variation of \texttt{map} but where the function it
is given always returns a list. When we use an ordinary \texttt{map}, and the
function that we give it returns a list, then the result of the \texttt{map}
will be a list of lists.  Sometimes we would rather have that we just got a
flattened list back. This is exactly what \texttt{map\_append} achieves.

The implementation looks a lot like the \texttt{map} implementation, but the
key difference is that we in the \texttt{cons\_case} have replaced the list
constructor with our implementation of \texttt{list\_append}. This means that
instead of constructing a list by using \texttt{cons} on lists, we construct a
list by appending lists together. In turn this yields a flattened list instead
of a list of lists.

\section{Back and forth}
In this section we will look at a connection between \texttt{fold\_right} and
\texttt{fold\_left}. We will in turn be able to define each in term of the
other. We will accomplish this with an accumulator.

\subsection{Fold left from fold right}
We implement this by letting the fold have a \texttt{nil\_case} function which
simply returns what it is given. Lets call this function $n$.

\begin{center}
  $n(a) = a$
\end{center}

And the \texttt{cons\_case} function which we will denote with a $c$, looks
like this:

\begin{center}
  $c(x, h, a) = h (cons\_case(x, a))$
\end{center}

Please note that these are the functions we pass to \texttt{fold\_right} as
\texttt{nil\_case} and \texttt{cons\_case}, and the $cons\_case$ mentioned here
is the one we would normally pass to \texttt{fold\_left}. We know from the
specification of folds that they return the same type as their
\texttt{nil\_case}, which means that when we give the above two to
\texttt{fold\_right}, the result will be a function taking one argument just
like the \texttt{nil\_case} function. This turns out to be a perfect fit for
the function $h$.

In order to better explain this, consider the list $cons(1,cons(2, nil))$.
\begin{align*}
        &fold\_left(nil\_case, cons\_case, cons(1,cons(2, nil)))
  \\ &= fold\_right(n, c, cons(1,cons(2, nil)))(nil\_case)
  \\ &= (a \rightarrow c(1, fold\_right(n, c, cons(2,nil), a)))(nil\_case)
  \\ &= c(1, fold\_right(n, c, cons(2,nil), nil\_case))
  \\ &= fold\_right(n, c, cons(2,nil))(cons\_case(1,nil\_case))
  \\ &= fold\_right(n, c, nil)(cons\_case(2,cons\_case(1,nil\_case)))
  \\ &= n(cons\_case(2,cons\_case(1,nil\_case)))
  \\ &= cons\_case(2,cons\_case(1,nil\_case))
\end{align*}

In the final line we see that we get exactly what we would expect an ordinary
\texttt{fold\_left} would do.

Proving that this implementation of \texttt{fold\_left} fits the specification,
and thereby always \emph{acts} as a \texttt{fold\_left} is done by a simple
structural induction proof on the input list.

\subsection{Fold right from fold left}
The implementation of \texttt{fold\_right} is very similar to that of the
\texttt{fold\_left} implementation which we just looked at. $n$ and $c$ will
remain the same, and if we do the same example as the one above we get:
\begin{align*}
        &fold\_right(nil\_case, cons\_case, cons(1,cons(2, nil)))
  \\ &= fold\_left(n, c, cons(1,cons(2, nil)))(nil\_case)
  \\ &= fold\_left((a \rightarrow c(1, n, a)), c, cons(2,nil))(nil\_case)
  \\ &= fold\_left((a \rightarrow n (cons\_case(1, a))), c, cons(2,nil))(nil\_case)
  \\ &= fold\_left((a \rightarrow cons\_case(1, a)), c, cons(2,nil))(nil\_case)
  \\ &= fold\_left((a \rightarrow (a \rightarrow
  cons\_case(1,a))(cons\_case(2,a))), c, nil)(nil\_case)
  \\ &= (a \rightarrow (a \rightarrow
  cons\_case(1,a))(cons\_case(2,a)))(nil\_case)
  \\ &= (a \rightarrow cons\_case(1,a))(cons\_case(2, nil\_case))
  \\ &= cons\_case(1,cons\_case(2, nil\_case))
\end{align*}

Again we see that the result is what we would expect from a
\texttt{fold\_right}.

In order to prover that this is the case for all lists and not just the one in
the example, we do a proof by structural induction on the list. Where the
\texttt{cons\_case} uses a master lemma.

\section{When cons case is an associative and commutative function}
When we wrote the code for the project, we noticed that some of the unit tests
for \texttt{fold\_right} and \texttt{fold\_left} were identical, and still
passed for both of the functions, while others only passed for one of them.

This happens when the \texttt{cons\_case} is an associative and commutative
function. In order for the \texttt{cons\_case} function to be associative and
commutative, it must take two arguments of the same type. This means that if we
can prove that the \texttt{cons\_case} function is associative and commutative,
then it should make no difference whether we choose to use \texttt{fold\_right}
or \texttt{fold\_left} --- The result should be the same.

When we did the project we started by showing that the property held for
\texttt{plus} as the \texttt{cons\_case} function. We did this as stepping
stone before writing a more generalised version of the property. One should
note that we choose to use \texttt{ring} as the tactic that applies the
associativity, commutativity and neutral element rules. We could have written
it out using the actual lemmas as well. When writing the generalised version,
one sees that the actual proofs are more or less the same thing.

\subsection{Structure of the proof}
The proof consists of a proof by structural induction on the input list.

In the \texttt{cons\_case} we will need a master lemma, which shows that after
applying the \texttt{cons\_case} for each of the folds once,
\texttt{fold\_right} can then be replaced by \texttt{fold\_left}. The master
lemma is also proven by structural induction. It is reasonable to this nested
induction because we will need three values in order to apply the associativity
rules.

The \texttt{cons\_case} of the master lemma is where we actually use the
assumption that the \texttt{cons\_case} function is associative and
commutative, and that is really the core of the proof.

\section{Primitive iteration and primitive recursion over polymorphic lists}
In the course for which this is a term project, we visited primitive recursion
and primitive iteration. In this section we will compare the two with
\texttt{fold\_left} and \texttt{fold\_right}.

The versions of primitive iteration and recursion that we looked at during the
course, worked on natural numbers. In order to compare with \texttt{fold\_left}
and right, we will need to specify and implement new versions of primitive
iteration and primitive recursion. First we will specify
\texttt{p\_i\_over\_polymorphic\_lists} (primitive iteration over polymorphic
lists), and compare this to \texttt{fold\_left} and \texttt{fold\_right}.
Afterwards we will specify \texttt{p\_r\_over\_polymorphic\_lists} (primitive
recursion over polymorphic lists) and do the same thing for that function.

\subsection{Primitive iteration over lists}
When we implement \texttt{p\_i\_over\_polymorphic\_lists}, we must define two
things.
\begin{enumerate}
  \item What to do when the list is empty.
  \item What to do when the list consists of a head and a tail
\end{enumerate}

In the first case, we return what we want to happen when we reach \texttt{nil}.
In the second case, we return the result of applying a \texttt{cons\_case}
function with the head of the list as the first argument and in the second
argument we call \texttt{p\_i\_over\_polymorphic\_lists}.

\subsubsection{Fold right from primitive iteration}
If we at this point stop and look at the specification of \texttt{fold\_right},
we will see that the specifications are in fact the same thing. If we have an
implementation of \texttt{p\_i\_over\_polymorphic\_lists}, we can easily get an
implementation of \texttt{fold\_right} because
\texttt{p\_i\_over\_polymorphic\_lists} also satisfies the specification of
\texttt{fold\_right}.

In fact this means that \texttt{p\_i\_over\_polymorphic\_lists} is the same thing as \texttt{fold\_right}.

\subsubsection{Fold left from primitive iteration}
We know that \texttt{p\_i\_over\_polymorphic\_lists} satisfies the
specification of \linebreak\texttt{fold\_right} and earlier in the project we
found a way of using a \texttt{fold\_right} to get a \texttt{fold\_left}. If we
combine this knowledge, we will have an implementation of \texttt{fold\_left}
from \texttt{p\_i\_over\_polymorphic\_lists}.

\subsection{Primitive recursion over lists}
\texttt{p\_r\_over\_polymorphic\_lists} is almost like
\texttt{p\_i\_over\_polymorphic\_lists}, only with a slight change to the
\texttt{cons\_case} function. The \texttt{cons\_case} function now takes one
more argument, which is a snapshot of the current list.

We can implement \texttt{fold\_right} and \texttt{fold\_left} in the exact same
way as with \texttt{p\_i\_over\_polymorphic\_lists}, if we just ignore that the
snapshot exists.

\subsubsection{Usage of the snapshot}
This project was focused on folds, and we can easily implement folds without
using the snapshots. But in order to show some example usage of the snapshot in
\texttt{p\_r\_over\_polymorphic\_lists}, we implemented a function which maps a
list to the list of it's suffixes.

The last two unit tests for \texttt{p\_r\_over\_polymorphic\_lists} show two
functions which make use of the snapshots. The first one returns a list of all
the snapshots, the second returns a flattened list of snapshots.

We took a step further and defined three versions of the
\texttt{list\_of\_suffixes} function (the first one of the two mentioned
above), where one of these use the snapshot like the one in the unit test. The
other two are alternative implementations. We then prove these to be
equivalent.

This is done by structural induction on the input list. When proving the
equivalency with the one that uses \texttt{p\_i\_over\_polymorphic\_lists}, we
need a master lemma. Apart from that the proof is straight forward.

\section{Conclusion}
In this report we looked at \texttt{fold\_right} and \texttt{fold\_left}. We
saw that they share the same expressive power, by showing that one could be
defined in terms of the other. Afterwards we found out that when the
\texttt{cons\_case} function that we pass to a fold function is associative and
commutative, then we get the same result no matter if we use
\texttt{fold\_left} or \texttt{fold\_right}.

During this we saw a recurring proof pattern where we in the induction case of a
structural induction proof, needed a master lemma in order to complete the
proof.

In the end we saw that \texttt{p\_i\_over\_polymorphic\_lists} is the exact
same thing as \texttt{fold\_right}. Using the connection between
\texttt{fold\_right} and \texttt{fold\_left}, we were able to define \texttt{fold\_left} in terms
of \texttt{p\_i\_over\_polymorphic\_lists}. We also saw how the snapshot in \texttt{p\_r\_over\_polymorphic\_lists} was
irrelevant in relation to folding.

\printbibliography[heading=bibnumbered]{}

\end{document}

ext: tex

max_stars_repo_name: blacksails/tfp-term-fold

max_issues_repo_name: blacksails/tfp-term-fold

max_forks_repo_name: blacksails/tfp-term-fold


text:
% 2012-12-16-handy,greg-cone-to-parallel.tex

\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{amssymb}
\usepackage[dvips,text={7.2in,9.2in}]{geometry}

\newcommand{\infig}[2][] {\includegraphics[#1]{fig/#2.pdf}}

\title{3D CT Cone Beam Reconstruction Via Fan to Parallel Rebinning}
\author{Gregory Handy}
\date{2012-12-16}

\begin{document}
\maketitle

\section{Overview}

This semester (Fall 2012), I spent time reading and implementing the paper by
Grass et al.\ titled ``3D cone-beam CT reconstruction for circular
trajectories''. This paper presented a new reconstruction algorithm that
modifies the well-known Feldkamp algorithm.  Their new method can be broken
into two major steps.  The first step rebins the standard 3D cone beam
geometry into a parallel-beam geometry along the source trajectory.  It is
possible to perform a filtered backprojection of the data at this step, and
this method is known as P-FDK (parallel-FDK). Despite the algorithmic
differences with FDK, this method does not differ in image quality/artifacts.
The second step interpolates the projections onto a rectangular virtual
detector plane and then proceeds with the filtered backprojection.  This
method is known as T-FDK (tent-FDK) and shows improved image quality.

In this paper I will provide a mathematical formalism of the algorithm, which
is not found in the Grass et al.\ paper, along with directions on how to use
the corresponding code.  The formulas will be consistent with the variables
used in the book chapters written by Fessler, and will reference sections out
of these chapters, code found in the irt toolbox, as well as new code.  The P-
FDK and T-FDK methods have been completely implemented, but both should be
tested on additional phantoms and real datasets.

\section{P-FDK}

Assume that we are given cone-beam projection data $p^f(s,t,\beta)$, where $f$
denotes the cone/fan-beam geometry, $s$ is the column position on the
detector, $t$ is the row position on the detector, and $\beta$ is the angular
position of the x-ray source.  Rebinning to parallel fixes $t$, and rebins in
the $s$ direction.  For each fixed $t$, we are in the geometry presented in
Section 3.9.

\subsection{Fan to Parallel Rebinning}

In Section 3.9.1, fan to parallel rebinning is defined assuming $r_{\mathrm{off}}=0$. We will also assume that we are using third-generation X-ray CT systems ($D_\mathrm{fs}=0$). This simplifies expression
(3.9.11) to
\begin{eqnarray*}
p^p(r,\phi) = p^f(D_{sd}\arcsin(r/D_{s0}),\phi-\arcsin(r/D_{s0})),
\end{eqnarray*}
where $\phi$ denotes the angle between the ray and the y-axis, $r$ is the distance from the center of rotation to the ray, $D_{s0}$ is the distance from the source to the center of rotation, and $D_{sd}$ is
the distance from the source to the detector. This geometry can be seen in Figure \ref{fig:above}. Extending this to our 3D case yields,
\begin{eqnarray*}
p^p(r,t,\phi) = p^f(D_{sd}\arcsin(r/D_{s0}),t,\phi-\arcsin(r/D_{s0})).
\end{eqnarray*}
The code implementing this step simply loops over the $t$-coordinate and makes
repeated calls to the function rebin\_fan2par.  It should be noted that this
function requires a 2D projection matrix and a 2D CT geometry.  Both of these
can be created using the 3D CT geometry and 3D projection matrix provided.
The user has control over the resulting parallel beam geometry (can change
orbit, ds, etc.). These changes may offer a slight improvement in image
quality.  When decreasing ds to make a finer grid, the only requirement is
that ds$\cdot$ns remains constant.

After rebinning, one must define a new 3D CT geometry, and make any necessary
changes (new ds, etc.). When doing so, the user should use the function
ct\_geom\_par.  This is essentially the same as ct\_geom, but allows for one
additional variable, called `rebinned', which should now be set to 1.  In the
future, the ct\_geom code should be rewritten to allow for a new type, called
par\_rebinned, but for now, the modified code checks if cg.rebinned == 1.
This small change allows for our 3D CT geometry to use variables such as
D$_\mathrm{sd}$, which can only be set to `inf' for the parallel type
currently.

\subsection{New Geometry}

In order to complete the reconstruction algorithm, the new geometry must be
understood.  Figure \ref{fig:above} offers an above view of the rebinnned
projection data while Figure \ref{fig:side} offers a side view.  The distances
in these figures will be referenced in the following sections.

It is important to note that the total distance from the source to the
detector (D$_\mathrm{sd}$) is kept constant.  As a result, the `new' detector
has an inverted

\begin{figure}[tbh]
\centering
\infig[width=0.98\linewidth]{officialAbove}
\caption{New geometry after rebinning. View is from above, and the x-y axis is provided. Note that the
total distance D$_\mathrm{sd}$ is preserved for the rebinned data, which causes the new detector
to curve outwards. Point of interest is the black dot.}
\label{fig:above}
\end{figure}

\begin{figure}[tbh]
\centering
\infig[scale=0.6]{officialSide}
\caption{The side view of the projection ray through the point of interest for the new geometry. The v-z axis is provided. The
total horizontal distance from the source to the detector is D$_\mathrm{sd}$.}
\label{fig:side}
\end{figure}
\clearpage
\noindent arc. Figure \ref{fig:side} is given relative to the v-z
axis, where $v(x,y,\phi) = y_\phi$, which can be found in Equation (3.9.16), and is reproduced below:
\begin{eqnarray*}
v(x,y,\phi) = (-x \sin(\phi) + y \cos(\phi)).
\end{eqnarray*}

\subsection{Filtered Backprojection}

The pre-weighting and filtering step proceed as normal with
\begin{eqnarray*}
w^p(t) &=& \frac{D_\mathrm{sd}}{\sqrt{D_\mathrm{sd}^2 + t^2}} \\
\hat{p}^p(r,t,\phi) &=& \Big(w^p(t) p^p(r,t,\phi)\Big) * g^p(r).
\end{eqnarray*}
The cosine weighting function can be derived from Figure \ref{fig:side}, and
$g^p(r)$ is the usual ramp filter used in the parallel-beam case (Equation
(3.4.14)).

The backprojection now follows with
\begin{eqnarray*}
f_p(x,y,z) = \int_0^{2\pi} \hat{p}^p(r(x,y,\phi),t(x,y,z),\phi) d\phi,
\end{eqnarray*}
where
\begin{eqnarray*}
r(x,y,\phi) &=& x \cos(\phi) + y \sin(\phi) \\
t(x,y,z,\phi) &=& \frac{D_\mathrm{sd}z}{\sqrt{D_\mathrm{s0}^2-r(x,y,\phi)^2}+v(x,y,\phi)}. \
\end{eqnarray*}
The formula for $r(x,y,\phi)$ follows directly from the parallel beam geometry. In order to derive
$t(x,y,z,\phi)$, note that by the similar triangles seen in Figure \ref{fig:side}, we have
\begin{eqnarray*}
\frac{t}{D_\mathrm{sd}} &=& \frac{z}{d'+v}, \
\end{eqnarray*}
where
\begin{eqnarray*}
d'&=& \sqrt{D_\mathrm{s0}^2-r^2}, \
\end{eqnarray*}
can be derived from Figure \ref{fig:above}.
It is interesting to note that additional weighting is not required
during the backprojection step.
Figure \ref{fig:comp} shows the reconstruction results using
the P-FDK and FDK algorithms.

\begin{figure}
\centering
\infig[scale=0.6]{comparison}
\caption{True and reconstruct images of a phantom for the first and middle slices, along with a side
view of the phantom.}
\label{fig:comp}
\end{figure}

\section{T-FDK}
Instead of filtering and backprojecting after rebinning to the parallel
geometry, the T-FDK method interpolates the projection data to form a
Cartesian grid for the projection data on a virtual detector found at the axis
of rotation.  Before interpolation, the grid at this virtual detector can be
seen in Figure \ref{grid}. Instead of simply interpolating (in the vertical
direction) a rectangle that is contained within this curved grid, we are
interested in extrapolating in order to create a rectangle that contains the
curved grid.

\begin{figure}[tbh]
\centering
\mbox{
\subfigure{\infig[width=2in]{rectGrid}
\quad
\subfigure{\infig[width=2in]{curved}}
}}
\caption{The left figure represents the projection grid on the detector, while the right figure represents the
same grid on the virtual detector located at the center of rotation.}
\label{grid}
\end{figure}

The function gridding.m performs this interpolation/extrapolation.  It takes
in the CT geometry and the projection matrix, and returns the new projection
matrix, along with coordinates for the new grid.  This function uses interp1
and currently performs linear interpolation, but this method can be changed.
We will refer to the new projection after interpolation/extrapolation in the
$t$-direction as $p^{pr}(r,t',\phi)$.

After performing the interpolation/extrapolation, the algorithm continues with pre-weighting
and filtering, according to the equations
\begin{eqnarray*}
w^{pr}(r,t') &=& \frac{\sqrt{D_\mathrm{s0}^2-r^2}}{\sqrt{D_\mathrm{s0}^2 -r^2+ t'^2}} \\
\hat{p}^{pr}(r,t',\phi) &=& \Big(w^{pr}(r,t') p^{pr}(r,t',\phi)\Big) * g^p(r).
\end{eqnarray*}
The new cosine weighting function, $w^{pr}(r,t')$, can be derived directly from Figure \ref{fig:side}.
Finally, the backprojection occurs with
\begin{eqnarray*}
f_{pr}(x,y,z) = \int_0^{2\pi} \hat{p}^{pr}(r(x,y,\phi),t'(x,y,z),\phi) d\phi,
\end{eqnarray*}
where $r(x,y,\phi)$ is the same as before, and
\begin{eqnarray*}
t'(x,y,z,\phi) &=& \frac{z\sqrt{D_\mathrm{s0}^2-r(x,y,\phi)^2}}{\sqrt{D_\mathrm{s0}^2-r(x,y,\phi)^2}+v(x,y,\phi)}. \
\end{eqnarray*}
Again, this follows from the similar triangle identity
\begin{eqnarray*}
\frac{t'}{d'} &=& \frac{z}{d'+v}. \
\end{eqnarray*}

\section{Conclusion/Future Work}

Preliminary results show similar images for all three algorithms (FDK, T-FDK,
and P-FDK). Future work should be done testing the methods on additional
phantoms and real CT datasets.  Additional work can be done on the
interpolation/extrapolation step.  Specifically, examining the effects of
different interpolation methods, and adjusting the dt and nt values of the new
grid.  Currently, the new grid keeps nt constant.  Similar to how our
rebinning steps can be used to decrease ds, the interpolation step can be used
to decrease dt.  Lastly, parallel\_feldkamp\_example.m provides an example
that runs the code found in the folder rebinnedFunctions, which can be placed
in the folder irt/fbp.

\section{References}

\cite{grass:00:3cb}
\cite[Section 3.9]{fessler::ir}
\cite{turbell:01:cbr}

\bibliographystyle{unsrt}
\bibliography{master}

\end{document}

max_stars_repo_name: tianrluo/mirt

max_issues_repo_name: tianrluo/mirt

max_forks_repo_name: tianrluo/mirt


text:
\chapter{Maxwell's Equations for an waveguide}
This chapter will dedicate to some fundamental theory of Maxwell's equations. More detailed
discussions can be found in Ref.~\cite{Snyder1983Optical}. Some detailed solution of Maxwell's equation
applied to cylindrical structures can be found in Ref.~\cite{Wait1959}.

For non-magnetic materials which normally constitute an optical waveguide with $ \mu=\mu_0 $, the
spatial dependence of the electrical field $ \boldsymbol{\mathcal{E}}(\br) $ and the magnetic field $
\boldsymbol{\mathcal{H}}(\br) $ of an optical waveguide is determined by Maxwell's equations:
\begin{align}
\nabla\times \boldsymbol{\mathcal{E}} &= i(\mu_0/\varepsilon_0)^{1/2} k \boldsymbol{\mathcal{H}}, &
\nabla\times \boldsymbol{\mathcal{H}} &= \boldsymbol{\mathcal{J}}-i(\varepsilon_0/\mu_0)^{1/2}kn^2
\boldsymbol{\mathcal{E}}, \label{EHtimesMKS}\\
\nabla\cdot (n^2 \boldsymbol{\mathcal{E}}) &= \rho/\varepsilon_0, & \nabla\cdot
\boldsymbol{\mathcal{H}} &=0, \label{EHdotsMKS}
\end{align}
where $ \boldsymbol{\mathcal{J}} $ and $ \rho $ are the current density and charge density, $
\varepsilon=n^2 \varepsilon_0 $ is the dielectric constant of the waveguide as a function of space, and $
k_0=2\pi/\lambda=\omega_0/c $.  We usually assume an implicit chromatic time dependence factor $ \exp(-i\omega_0 t) $ in the
full field vectors. These equations above are written in MKS units.

Sometimes, to avoid the dielectric constant $ \varepsilon_0 $ and magnetic permeability $ \mu_0 $ of
free space, we can rewrite the Maxwell's equations in Gaussian-cgs units, which yield
\begin{align}
\nabla \times \boldsymbol{\mathcal{E}} &= -\frac{1}{c} \pt{\boldsymbol{\mathcal{B}}}=i
k_0\boldsymbol{\mathcal{B}}, & \nabla \times
\boldsymbol{\mathcal{H}} &= \frac{4\pi}{c} \boldsymbol{\mathcal{J}} +
\frac{1}{c}\pt{\boldsymbol{\mathcal{D}}}= \frac{4\pi}{c}\boldsymbol{\mathcal{J}}
-ik_0\boldsymbol{\mathcal{D}}, \label{eq:maxwelldiv}\\
\nabla \cdot \boldsymbol{\mathcal{D}} &= 4\pi \rho, & \nabla \cdot \boldsymbol{\mathcal{B}} &=0,\label{eq:maxwellgrad}
\end{align}
where $ \boldsymbol{\mathcal{D}} = \varepsilon \boldsymbol{\mathcal{E}}=\varepsilon_r \boldsymbol{\mathcal{E}} $ and $
\boldsymbol{\mathcal{B}}= \mu\boldsymbol{\mathcal{H}}\approx \boldsymbol{\mathcal{H}} $ for non-magnetic materials. By default, we will use Gauss units in this piece of work, and only consider non-magnetic materials. Notice that, we make the dielectric constant equal to the relative dielectric constant ($ \varepsilon=\varepsilon $) in the Gaussian-cgs units convention.

Equations~\eqref{eq:maxwelldiv} can lead to the wave equation for a chromatic electric field in Gauss units as
\begin{align}
-\nabla \times (\nabla \times \boldsymbol{\mathcal{E}})+\varepsilon k_0^2\boldsymbol{\mathcal{E}}=-ik_0\frac{4\pi}{c} \boldsymbol{\mathcal{J}}. \label{eq:waveeqGaussU}
\end{align}
Similarly, the corresponding wave equation in SI units can be given by
\begin{align}
-\nabla \times (\nabla \times \boldsymbol{\mathcal{E}})+ n^2 k_0^2\boldsymbol{\mathcal{E}}=-ik_0\sqrt{\frac{\mu_0}{\varepsilon_0}} \boldsymbol{\mathcal{J}}.
\end{align}

We can further simplify the equation above by using
\begin{align}
\nabla\times (\nabla\times \bmc{E}) &= -\nabla(\nabla\cdot \bmc{E})+ \nabla^2\bmc{E}.
\end{align}
In the case that there is no net charge and current in the space, the wave equations above in both Gauss and SI units can be simplified to be
\begin{align}\label{eq:freespacewaveeq}
(\nabla^2 + n^2k_0^2) \boldsymbol{\mathcal{E}}=0.
\end{align}

\section{Fields of translationally invariant waveguides}
We define the axis of the waveguide is along $ z $-direction, and the refractive index profile of the
waveguide is independent of $ z $, i.e. $ n=n(\br\!_\perp) $, which means the waveguide is
translationally invariant. The fields of the waveguide can then be rewritten in a separable form
\begin{align}
\bmc{E}(\br) &= \bmc{E}(\br\!_\perp)\exp(i\beta z), & \bmc{H}(\br) &= \bmc{H}(\br\!_\perp)\exp(i\beta
z),
\end{align}
where $ \beta $ is the propagation constant and $ \br\!_\perp $ is the position vector in the transverse
plane perpendicular to the $ z $-axis. We can further decompose these fields into longitudinal and
transverse components, parallel and orthogonal to the waveguide axis, respectively, and denoted by
subscripts $ z $ and $ \perp $. That is
\begin{align}\label{EHtranslong}
\bmc{E}(\br) &= \left( \bmc{E}\!_\perp + \mathcal{E}_z \hat{e}_z \right)\exp(i\beta z), & \bmc{H}(\br) &=
\left( \bmc{H}\!_\perp + \mathcal{H}_z \hat{e}_z\right) \exp(i\beta z)
\end{align}

Notice that, in literature, many authors use the sign parameter $ f=\pm 1 $ in front of $ \beta $ to indicate the propagating direction of the fields. Here, we define forward- and backward-propagating waves using the sign of the product of $ \beta z $. This is convenient to differ the contour integration paths when we decompose the radiation and guided modes for $ z=z-z'>0 $ or $ z=z-z'<0 $ cases.

\section{Relationships between field components}\label{MWE:components}
By substituting Equ.~\ref{EHtranslong} into the source-free Maxwell equations
(Equs.~\ref{EHtimesMKS} and~\ref{EHdotsMKS} with $ \bmc{J}=0,\, \rho=0 $), and comparing
longitudinal and transverse components, we obtain the relationships between fields components as
below
\begin{subequations}
\label{EHtz0}
\begin{align}
\bmc{E}\!_\perp &= -\left(\frac{\mu_0}{\varepsilon_0} \right)^{1/2} \frac{1}{k_0n^2} \hat{e}_z \times
\left(\beta \bmc{H}\!_\perp + i\nabla\!_\perp \mathcal{H}_z \right),\\
\bmc{H}\!_\perp &= \left(\frac{\varepsilon_0}{\mu_0} \right)^{1/2} \frac{1}{k_0} \hat{e}_z \times \left(\beta
\bmc{E}\!_\perp +i\nabla\!_\perp \mathcal{E}_z \right),\\
\mathcal{E}_z &= i\left(\frac{\mu_0}{\varepsilon_0} \right)^{1/2} \frac{1}{k_0n^2} \hat{e}_z \cdot
\nabla\!_\perp \times \bmc{H}\!_\perp = \frac{i}{\beta} \left( \nabla\!_\perp \cdot \bmc{E}\!_\perp +
(\bmc{E}\!_\perp \cdot \nabla\!_\perp) \ln n^2 \right),\\
\mathcal{H}_z &= -i \left(\frac{\varepsilon_0}{\mu_0} \right)^{1/2} \frac{1}{k_0}\hat{e}_z \cdot
\nabla\!_\perp \times \bmc{E}\!_\perp = \frac{i}{\beta} \nabla\!_\perp \cdot \bmc{H}\!_\perp.
\end{align}
\end{subequations}

Now we consider the waveguide of cylindrical fiber case. Due to symmetry, we tend to use the
cylindrical coordinate system, which gives
\begin{align}
\nabla\!_\perp = \hat{r}\!_\perp \pp{}{r\!_\perp} + \hat{\phi} \frac{1}{r\!_\perp} \pp{}{\phi}.
\end{align}
Therefore, Equ.~\ref{EHtz0} implies that the transverse components can be expressed in terms of the
longitudinal components by
\begin{subequations}
\begin{align}
\mathcal{E}_{r\!_\perp} &= \frac{i}{\kappa_i^2} \left[ \beta \pp{\mathcal{E}_z }{r\!_\perp} +
\left(\frac{\mu_0}{\varepsilon_0} \right)^{1/2}  \frac{k_0}{r\!_\perp}
\pp{\mathcal{H}_z}{\phi} \right], \\
\mathcal{E}_\phi &= \frac{i}{\kappa_i^2} \left[\frac{\beta}{r\!_\perp}
\pp{\mathcal{E}_z}{\phi} -\left(\frac{\mu_0}{\varepsilon_0} \right) ^{1/2} k_0
\pp{\mathcal{H}_z}{r\!_\perp}
\right],\\
\mathcal{H}_{r\!_\perp} &= \frac{i}{\kappa_i^2} \left[ \beta \pp{\mathcal{H}_z}{r\!_\perp}
-\left(\frac{\varepsilon_0}{\mu_0} \right)^{1/2} \frac{k_0n^2}{r\!_\perp}\pp{\mathcal{E}_z}{\phi} \right],\\
\mathcal{H}_\phi &= \frac{i}{\kappa_i^2} \left[ \frac{\beta}{r\!_\perp} \pp{\mathcal{H}_z}{\phi} +
\left(\frac{\varepsilon_0}{\mu_0} \right)^{1/2} k_0n^2 \pp{\mathcal{E}_z}{r\!_\perp} \right],
\end{align}
\end{subequations}
with $ \kappa_i^2= k_0^2n_i^2-\beta^2=k_0^2 \varepsilon_i -\beta^2  $ and $ n=n(\br\!_\perp)=n_i$ for corresponding region $i$.

Correspondingly, the component relationships in Gaussian-cgs units are
\begin{subequations}\label{EHzgauss}
\begin{align}
\mathcal{E}_{r\!_\perp} &= \frac{i}{\kappa_i^2} \left[ \beta \pp{\mathcal{E}_z }{r\!_\perp} +
  \frac{k_0}{r\!_\perp}
\pp{\mathcal{H}_z}{\phi} \right]= \frac{i\beta}{\kappa_i^2} \pp{\mathcal{E}_z }{r\!_\perp} -
  \frac{k_0m}{r\!_\perp \kappa_i^2} {\mathcal{H}_z}, \\
\mathcal{E}_\phi &= \frac{i}{\kappa_i^2} \left[\frac{\beta}{r\!_\perp}
\pp{\mathcal{E}_z}{\phi} - k_0
\pp{\mathcal{H}_z}{r\!_\perp}
\right] = -\frac{\beta m}{r\!_\perp \kappa_i^2}
{\mathcal{E}_z} - \frac{ik_0}{\kappa_i^2}
\pp{\mathcal{H}_z}{r\!_\perp},\\
\mathcal{H}_{r\!_\perp} &= \frac{i}{\kappa_i^2} \left[ \beta \pp{\mathcal{H}_z}{r\!_\perp}
- \frac{k_0n^2}{r\!_\perp}\pp{\mathcal{E}_z}{\phi} \right]= \frac{i\beta}{\kappa_i^2} \pp{\mathcal{H}_z}{r\!_\perp}
+ \frac{k_0n^2m}{r\!_\perp \kappa_i^2} {\mathcal{E}_z},\\
\mathcal{H}_\phi &= \frac{i}{\kappa_i^2} \left[ \frac{\beta}{r\!_\perp} \pp{\mathcal{H}_z}{\phi} +
 k_0n^2 \pp{\mathcal{E}_z}{r\!_\perp} \right] = -\frac{\beta m}{r\!_\perp \kappa_i^2} {\mathcal{H}_z} +
  \frac{ik_0n^2}{\kappa_i^2} \pp{\mathcal{E}_z}{r\!_\perp}.
\end{align}
\end{subequations}

%\textcolor{red}{Vector and Scalar operators in various coordinate systems...}


text:
\begin{quote}
	I am motivated to present this architecture, which is linked to antiquity and doubtless to other cultures, because it is an elegant and lively witness to what I have tried to define as an outside-time category, \textit{algebra}, or structure of music, as opposed to its other two categories, in-time and temporal. \im \parencite[192]{Xen92:For}
\end{quote}

\begin{quote}
	With [the relational] model any formatted data base is viewed as a collection of time-varying relations of assorted degrees\dots this collection is called a \textit{relational algebra}\dots a query language could be directly based on it\dots The primary purpose of [relational] algebra is to provide a collection of operations on relations of all degrees\dots suitable for selecting data from a relational data base. \im \parencite[1-5]{Codd72relationalcompleteness}
\end{quote}

\paragraph{Querying the Sieves}
If we consider pitches as an outside-time (relational) database, one way of understanding Xenakis' sieve theory is as a query method, for which E.F. Codd's \nameref{model:relational} model would fit perfectly. The nature of this consideration stems from the application of algebra as a programmable selection mechanism or simply, filters. Both concepts (sieves and relational algebra) have a common link which, not surprisingly, is the \gls{ibm-7090} mainframe computer.\footnote{Among other things, the \gls{ibm-7090} computer was used in the computation of the first 100,000 digits of $\pi$ \parencite{picalc}, Roger Shepard's computation of the homonymous `shepard' tone \parencite{shepard}, Alexander Hurwitz's computation of the 19th and 20th mersenne prime numbers,\footnote{\url{https://www.mersenne.org/primes/}} and Peter Sellers' plot-twisting moment in Stanley Kubrick's ``Dr. Strangelove or: How I Learned to Stop Worrying and Love the Bomb:'' \url{https://en.wikipedia.org/wiki/Dr._Strangelove}} While Xenakis' experiments were carried out on the \gls{ibm-7090} mainframe computer located at IBM-France in Paris, Codd himself worked at the IBM Research Laboratory in San Jose, California. Furthermore, this same computer was used by Hiller and Baker in their realization of \gls{musicomp}, a pioneering language for algorithmic composition \parencite[44]{Ari05:Ano}. Most important, the \gls{ibm-7090} used the programming language FORTRAN IV, as can be seen by the printed FORTRAN routines for Xenakis' 1962 work \textit{Atrées (ST/10-3 060962)} \parencite[145]{Xen92:For}.\footnote{Interestingly, given that \citeauthor{arizaSieves} finds Xenakis' sieves code unusable \parencite[1]{arizaSieves}, chances are that the printed code for the ST/10-3 composition is likewise useless.} Xenakis's work on sieves came a few years after his experiments on the \gls{ibm-7090}, and his sieves program was written in Basic and then in C. However, the experience with FORTRAN IV at the \gls{ibm-7090} serves nonetheless as a common ancestor to both Xenakis and Codd. For example, Xenakis' transcriptions in early \gls{cac} systems were performed with tables of outputted computer data. Further, \citeauthor{Ari05:Ano} \parencite{Ari05:Ano} writes how ``the early systems of Hiller, Xenakis, and Koenig all required manual transcription of computer output into Western notation. The computer output of these early systems was in the form of \textit{alpha-numeric data tables}: each row represents an event, each column represents an event parameter value'' \im \parencite[94]{Ari05:Ano}. In this sense, performance in \gls{cac} meant interpreting results out of a database.\footnote{For further reference on the early uses of computers in \gls{cac}, I refer the reader to \citeauthor{Ari05:Ano}'s PhD thesis \parencite{Ari05:Ano}.}

% This should be a footnote as it is not directly related to your argument byt only as a reference for speed. or am I reading something wrong?

\paragraph{Sound Synthesis Parenthesis}
(Before continuing, a sound synthesis parenthesis must be opened. While Xenakis praised the speed at which the \gls{ibm-7090} could perform computations, Max Mathews \parencite{Mat63:The}, then director of the Behavioral Research Laboratory at Bell Telephone Laboratories, wrote:

\begin{quote}
	A high-speed machine such as the \gls{ibm-7090}, using the programs described later in this article, can compute \textit{only about 5000 numbers per second} when generating a reasonably complex sound. However, the numbers can be temporarily stored on one of the computer's digital magnetic tapes, and this tape can subsequently be replayed at rates up to 30,000 numbers per second (each number being a 12-bit binary number). \im \parencite[553]{Mat63:The}
\end{quote}

Mathews' concern for speed was grounded on the need to achieve sound synthesis, which meant fast computations of the sample theorem. Initially, the first synthesized sound was obtained in 1957, with the (assembly-code written) MUSIC 1 program with the IBM-704 (a predecessor of the \gls{ibm-7090}). Later, when Bell Labs obtained the IBM-7094 ---which ``was a very, very effective machine'' \parencite[16]{Roa80:Int}---, and in combination with the (then) widely available FORTRAN compiler, Mathews could develop the MUSIC I program, into MUSIC V, which became the first portable computer music language designed for computer music synthesis.\footnote{As an example, I would refer the reader to James Tenney's work from 1962 ``Five Stochastic Studies,'' which can be found on a YouTube account on his name: \url{https://www.youtube.com/channel/UCEzSaoPnxCJVzXxA9obuRWg/videos}. Roads, while interviewing Matthews recalls this piece to be named ``Noise Studies'' \parencite[18]{Roa80:Int}, which fades out the reference to Xenakis' music.} I return to this discussion in \ref{improv})

\paragraph{Algebraic Abstractions for Freedom}
Xenakis' and Codd's papers came out around the same time: Xenakis' English publication of \textit{Towards a Metamusic} was in 1970, Codd's papers were published in 1970 and 1972. Sieve theory was aimed at providing a plethora of computable sets (or relations) of pitches, according to different temperings of the smallest displacement unit and the selected value for the modulo operator. Contemporaneously, Codd's relational algebra was meant to be the internal structure of a query language for selecting elements based on their relations. Both of these can be considered algebraic abstractions of a selection process. In the case of Xenakis, the abstraction was outside-time, as the composer could make a snapshot, or a tomography of the pitch space in order to analyze it and extrapolate structural relations. In Codd's case, the abstraction was spatial: the query language would be separated from the database, allowing databasers to perform queries in the `frontend' without worrying about internal data structures, memory allocation, and so on, since these operations would occur in the background of the `backend.' By black-boxing and hardware-specific programming, both methods freed the human operator to devise any kind of algebraic queries. Thus, the operator worked on a higher level of abstraction that enabled a less problematic kind of envisioning. Xenakis wrote convinced: ``freed from tedious calculations, the composer is able to devote himself to the general problems that the new musical form poses, and to explore the nooks and crannies of this form while \textit{modifying the values of the input data}\dots'' \im \parencite[144]{Xen92:For}.

\paragraph{A Cosmic Vessel and an Armchair}
Therefore, the composer delegates to the computer the minutiae of arduous iterative computations: precisely what the computer is better at than the human. As a result, in Xenakis' view, and in resonance with programmer Charles Bachman's claim in \textit{The Programmer as Navigator} \parencite{Bachman:1973:PN:355611.362534}, the composer became a pilot:

\begin{quote}
	With the aid of electronic computers the composer becomes a sort of pilot: he presses the buttons, introduces coordinates, and supervises the controls of \textit{a cosmic vessel sailing in the space of sound}, across sonic constellations and galaxies that he could formerly glimpse only as a distant dream. \textit{Now he can explore them at his ease, seated in an armchair} \im \parencite[144]{Xen92:For}
\end{quote}

Codd's and Xenakis' propositions were abstractions deeply rooted in and contextualized against a backdrop of their own fields. Xenakis wrote against the current state of Western Music with its ``degradation of outside-time structures'', the ``followers of information theory'' and the ``intuitionists.''\footnote{Xenakis defined both groups as poles on an extreme. On the first extreme (``followers of information theory''), and on the other (``intuitionists''), Xenakis identified two groups: ``graphists'' and those involved with ``happenings'' \parencite[180]{Xen92:For}. We can only speculate that Xenakis was referring to two American music currents during the 1950s: Hiller and Isaacson can be linked to the information theorists, and the intuitionists can be represented by Cage and the New York School composers.} Codd wrote against the previously developed hierarchical and network database models. Most important, these tools and their development had the human operator's considerations in mind. The composer, like the databaser, would engage in a rudimentary and limited, but still present, feedback process at the \textit{input} level. That is to say, unless rewriting the code, which consisted in a very long and economically expensive process combining punch cards and magnetic tapes, the composer and the databaser could change the input several times, achieving different outputs in a matter of hours.\footnote{As a reference, the computation of the first 100,000 values of $\pi$ took about eight and a half hours \parencite{picalc}.} For example, queries made on the relational model would appear on screen at a very fast rate, thus enabling better tuning of the input in relation to a wanted output. Likewise, the composer could modify the input values to highly complex calculations that would otherwise take a long time, or be error prone. The limitation, of course, is the level of intervention with the code, which the overall circuitry would thus complicate; criticism on account of this shortcoming of the circuit would thus be rendered anachronistic, but recalling these limitations places composition and databasing in perspective.


text:
% Created 2016-05-03 tis 13:36
% Intended LaTeX compiler: pdflatex
\documentclass{scrartcl}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{grffile}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{textcomp}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\usepackage{khpreamble}
\author{Kjartan Halvorsen}
\date{Due Friday 2016-04-22}
\title{Computerized control - homework 4}
\hypersetup{
 pdfauthor={Kjartan Halvorsen},
 pdftitle={Computerized control - homework 4},
 pdfkeywords={},
 pdfsubject={},
 pdfcreator={Emacs 24.5.1 (Org mode 8.3.4)},
 pdflang={English}}
\begin{document}

\maketitle

\section*{Active suspension}
\label{sec:orgheadline5}
A so-called \emph{quarter model} for the suspension of a car is shown in figure \ref{fig:quarter}. The system consists of two masses, one is the mass of 1/4 of the car, the other is the much smaller suspension mass. The two masses are connected by two passive elements: a spring and a damper, as well as an active element: a linear force actuator. The suspension mass is connected with the ground via a spring, representing the tyre (sometimes there is also a damper included in the model of the tyre).
\begin{figure}
\begin{center}
\includegraphics[width=0.35\linewidth]{../notes/figures/active-suspension-nodamper}
\caption{Active suspension model. The distances $z_1$, $z_2$ and $w$ are displacements from a static equilibrium. The displacements are measured with respect to a stationary frame of reference.}
\label{fig:quarter}
\end{center}
\end{figure}

A slightly more complex model is discussed \href{http://ctms.engin.umich.edu/CTMS/index.php?example=Suspension&section=ControlStateSpace}{here}. That model also includes a damping in the tyre. You are welcome to study the example. The parameter values are the same as for this homework assignment.

\subsection*{Determine the state space model}
\label{sec:orgheadline1}

\begin{enumerate}
\item Draw free-body diagrams for each of the two bodies, and set up the equation of motion of each of the two masses. Note that the variables \(z_1\), \(z_2\) and \(w\) are displacements from a static equilibrium, measured with respect to a stationary frame of reference. This means that the force acting on mass \(m_1\) from the spring between the two masses is \(F_s = m_1g -k_1(z_1-z_2)\), and the force acting on the suspension mass from the spring (tyre) between the suspension mass and ground is \(F_t= (m_1+m_2)g - k_2(z_2 - w)\). Moreover, the force in the damper is proportional with the relative velocity between the car and the suspension mass: \(F_d = -b_1(\dot{z}_1 - \dot{z}_2)\).
\item Show that the following state-space model of the system is obtained with the choice
\[ x = \bbm x_1\\x_2\\x_3\\x_4 \ebm = \bbm z_1 - z_2\\ \dot{z}_1 - \dot{z}_2 \\ z_2\\ \dot{z}_2 \ebm \]
\begin{align*}
\dot{x} &= \bbm 0 & 1 & 0 & 0\\
-\Big(\frac{1}{m_1} + \frac{1}{m_2}\Big)k_1 & -\Big(\frac{1}{m_1} + \frac{1}{m_2}\Big)b_1 & \frac{k_2}{m_2} & 0\\
0 & 0 & 0 & 1\\
              \frac{k_1}{m_2} & \frac{b_1}{m_2} & - \frac{k_2}{m_2} & 0\ebm
              \bbm x_1\\x_2\\x_3\\x_4 \ebm + \bbm 0\\-\frac{k_2}{m_2}\\0\\\frac{k_2}{m_2} \ebm w + \bbm 0\\\frac{1}{m_1} + \frac{1}{m_2}\\0\\-\frac{1}{m_2} \ebm u \\
y &= \bbm 1 & 0 & 0 & 0 \ebm x.
\end{align*}
Note that the model has two inputs: The force from the active suspension \(u\), and the disturbance from the road \(w\).
\end{enumerate}
\subsection*{Simulate step-responses}
\label{sec:orgheadline2}
Implement the model in matlab. Use the following values for the parameters
\begin{center}
\begin{tabular}{l}
\(m1 = 2500\)\\
\(m2 = 320\)\\
\(k1 = 80000\)\\
\(k2 = 500000\)\\
\(b1 = 350\)\\
\end{tabular}
\end{center}
Note that the model has two input signals: The disturbance from the road, \(w\), and the actuator input \(u\). Simulate a step-response from each of the two input signals. Below is some code to help you
\begin{verbatim}
m1 = 2500;
m2 = 320;
k1 = 80000;
k2 = 500000;
b1 = 350;

A =   % Your code here
Bu = [0
     1/m1+1/m2
     0
     -1/m2 ];
Bw = [ 0
      -k2/m2
      0
      k2/m2];
C=[1 0 0 0];
D=[0];
sys_uy=ss(A,Bu,C,D);
sys_wy=ss(A,Bw,C,D);

% Step responses
figure(1)
clf
step(sys_uy, sys_wy)

% Simulate step response for first half second only
T = linspace(0, 0.5, 800);
figure(2)
clf
subplot(121)
stepplot(sys_uy, T);
title('Response from u')
xlabel('Time (seconds)')
ylabel('y (meters)')
subplot(122)
stepplot(sys_wy, T);
title('Response from w')
xlabel('Time (seconds)')
ylabel('y (meters)')
\end{verbatim}

\subsection*{Sample the system}
\label{sec:orgheadline3}
\begin{enumerate}
\item The system is quite oscillative. Choose a sampling period \(h\), such that \(\omega_0h = 0.2\), where \(\omega_0\) is the period of the oscillations in radians.
\item Sample the system (numerically, using \texttt{c2d} in matlab)
\end{enumerate}

\subsection*{Observability and reachability}
\label{sec:orgheadline4}
Determine observability and reachability of the sampled system. There are matlab commands to do this for you (\texttt{obsv} and \texttt{ctrb}), but do it also by forming the observability and controllability matrices explicitly. Recall
\[ W_C = \bbm \Gamma & \Phi\Gamma & \cdots & \Phi^{n-1}\Gamma \ebm \]
and
\[ W_O = \bbm C\\C\Phi\\\vdots\\C\Phi^{n-1} \ebm. \]

\section*{Solutions}
\label{sec:orgheadline10}

\subsection*{State space model}
\label{sec:orgheadline6}
\begin{enumerate}
\item There are four forces acting on the car body
\begin{center}
\begin{tikzpicture}
\node (M1) [draw, align=center, minimum width=3.5cm,minimum height=2cm] {Body mass};
\draw[->, thick] (M1.center) -- node[left, near end] {$m_1g$} ++(0, -3cm);
\draw[<-, thick] (M1.south) ++(-1cm, 0) -- node[left, near end] {$F_s$} ++ (0, -3cm);
\draw[<-, thick] (M1.south) ++(1cm, 0) -- node[left, near end] {$F_d$} ++ (0, -3cm);
\draw[<-, thick] (M1.south) ++(0.2cm, 0) -- node[right, near end] {$u$} ++ (0, -1.76cm);
\end{tikzpicture}
\end{center}
The equation of motion in the vertical direction becomes
\begin{equation}
 m_1 \ddot{z}_1 = \sum_i F_i = -m_1g + \underbrace{m_1g -k_1(z_1-z_2)}_{F_s} + \underbrace{\big(-b_1(\dot{z}_1 - \dot{z}_2) \big) }_{F_d} + u  = -k_1(z_1 - z_2) - b_1(\dot{z}_1 - \dot{z}_2) + u.
 \label{eq:eom1}
 \end{equation}
On the suspension mass, there are five forces acting:
\begin{center}
\begin{tikzpicture}
\node (M2) [draw, align=center, minimum width=3.5cm,minimum height=2cm] {Suspension mass};
\draw[->, thick] (M2.center) -- node[left, near end] {$m_2g$} ++(0, -3cm);
\draw[<-, thick] (M2.north) ++(-1cm, 0) -- node[left, near end] {$F_s$} ++ (0, 3cm);
\draw[<-, thick] (M2.north) ++(1cm, 0) -- node[left, near end] {$F_d$} ++ (0, 3cm);
\draw[<-, thick] (M2.north) ++(0cm, 0) -- node[right, near end] {$u$} ++ (0, 1.76cm);
\draw[<-, thick] (M2.south) ++(0.2cm, 0) -- node[right, near end] {$F_t$} ++ (0, -3.76cm);
\end{tikzpicture}
\end{center}
The equation of motion in the vertical direction becomes
\begin{equation}
\begin{split}
m_2\ddot{z}_2 &= -m2_g - \underbrace{\big(m_1g -k_1(z_1-z_2)\big)}_{F_s}
 - \underbrace{\big(- b_1(\dot{z}_1-\dot{z}_2)\big)}_{F_d}
 - \underbrace{\big( (m_1 + m_2)g - k_2(z_2 - w) \big)}_{F_t} - u\\
	       &= k_1(z_1 - z_2) + b_1(z_1 - z_2) - k_2(z_2-w) - u
 \end{split}
 \end{equation}

\item With the state vector
\[ x = \bbm x_1\\x_2\\x_3\\x_4 \ebm = \bbm z_1 - z_2\\ \dot{z}_1 - \dot{z}_2 \\ z_2\\ \dot{z}_2 \ebm \]
we can write the four first-order differential equations of the state space model
 \begin{align*}
  \dot{x}_1 &= \dot{z}_1 - \dot{z}_2 = x_2\\
  \dot{x}_2 &= \ddot{z}_1 - \ddot{z}_2 = \frac{1}{m_1}\big( -k_1(z_1 - z_2) - b_1(\dot{z}_1 - \dot{z}_2) + u \Big) - \frac{1}{m_2} \Big(k_1(z_1 - z_2) + b_1(z_1 - z_2) - k_2(z_2-w) - u\Big)  \\
            &= -\Big(\frac{1}{m_1} + \frac{1}{m_2}\Big)k_1x_1 - \Big(\frac{1}{m_1} + \frac{1}{m_2}\Big)b_1x_2 + \frac{k_2}{m_2}x_3 - \frac{k_2}{m_2}w + \Big(\frac{1}{m_1} + \frac{1}{m_2}\Big) u\\
  \dot{x}_3 &= \dot{z}_2 = x_4\\
  \dot{x}_4 &= \ddot{z}_2 = \frac{1}{m_2} \Big( k_1(z_1 - z_2) + b_1(z_1 - z_2) - k_2(z_2-w) - u \Big)\\
            &= \frac{k_1}{m_2} x_1 + \frac{b_1}{m_2} x_2 - \frac{k_2}{m_2}x_3 + \frac{k_2}{m_2}w - \frac{1}{m_2}u.
\end{align*}
In state-space form, this can be written
\begin{align*}
\dot{x} &= \bbm 0 & 1 & 0 & 0\\
               -\Big(\frac{1}{m_1} + \frac{1}{m_2}\Big)k_1 & \Big(\frac{1}{m_1} + \frac{1}{m_2}\Big)b_1 & \frac{k_2}{m_2} & 0\\
               0 & 0 & 0 & 1\\
               \frac{k_1}{m_2} & \frac{b_1}{m_2} & - \frac{k_2}{m_2} & 0\ebm
               \bbm x_1\\x_2\\x_3\\x_4 \ebm + \bbm 0\\-\frac{k_2}{m_2}\\0\\\frac{k_2}{m_2} \ebm w + \bbm 0\\\frac{1}{m_1} + \frac{1}{m_2}\\0\\-\frac{1}{m_2} \ebm u \\
 y &= \bbm 1 & 0 & 0 & 0 \ebm x.
\end{align*}
\end{enumerate}

\subsection*{step-response}
\label{sec:orgheadline7}
The figure below shows the step-responses of both the continuous- and the sampled systems.
\begin{center}
\includegraphics[width=0.6\linewidth]{active-susp-plant-response}
\end{center}

\subsection*{Sampling the system}
\label{sec:orgheadline8}
\begin{enumerate}
\item The fast oscillations in the beginning of the step-response have period \unit{0.13}{\second},  so the frequency is \(\omega_0 = \frac{2\pi}{0.13}\approx \unit{48}{\radians\per\second}\). This gives
\[ h = \frac{0.2}{\omega_0} \approx \unit{0.0042}{\second}. \]
\item Sampling the system is straightforward
\begin{verbatim}
sys_uy_d = c2d(sys_uy, h); % The system from input u to y
sys_wy_d = c2d(sys_wy, h); % From input w to y
\end{verbatim}
\end{enumerate}

\subsection*{Observability and reachability}
\label{sec:orgheadline9}
Obtain the \(\Phi\), \(\Gamma\), \(C\) and \(D\) matrices from the discretized system. Then form \(W_o\) and \(W_c\) and calculate the determinant
\begin{verbatim}
% Check observability and reachability
[Ad, Bd, Cd, Dd] = ssdata(sys_uy_d);
Wo = [Cd;Cd*Ad;Cd*Ad*Ad; Cd*Ad*Ad*Ad]
det(Wo)
Wc = [Bd  Ad*Bd  Ad*Ad*Bd Ad*Ad*Ad*Bd]
det(Wc)
\end{verbatim}
The result is
\[ \det W_o = 1.26 \cdot 10^{-8}\]
and
\[ \det W_c = 5.92 \cdot 10^{-30}.\]

So, the system is observable, but not reachable.
\end{document}