text:
### A Pluto.jl notebook ###
# v0.11.2

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7548fe76-b328-11ea-1dff-b16fc38e0cc5
#O comando md seguido de três aspas duplas inicia o modo Markdown, para fechá-lo basta usar novamente três aspas duplas.
md"""

# Julia 101: Básico interativo com Pluto.jl

Vamos começar falando de algumas coisas que só farão sentido depois de um tempo.

Como estamos usando o [Pluto.jl](https://github.com/fonsp/Pluto.jl), então cada célula deveria ser *atômica*, ou seja, conter somente uma atribuição ou operação. No entanto, em alguns casos, podemos fazer várias operações em uma mesma célula por meio de blocos de [*begin*](https://docs.julialang.org/en/v1/base/base/#begin) com `begin ... end` ou de [*let*](https://docs.julialang.org/en/v1/manual/variables-and-scoping/#Let-Blocks-1) com `let ... end`. Veremos mais detalhes sobre esses blocos mais adiante!

Para importar um pacote, usamos o comando `using` ou `import`. Em ambos os casos, estamos carregando e permitindo que os *nomes* contidos no pacote sejam visíveis e acessíveis a partir do nosso *espaço de nomes*. A [diferença](https://docs.julialang.org/en/v1/manual/modules/#modules-1) é principalmente quanto à disponibilidade dos *nomes* e acesso aos *tipos* contidos no pacote.

Ao usar `using`, traremos os *nomes* do pacote para o nosso *espaço de nomes*. Por outro lado, ao usarmos `import`, apenas permitiremos o acesso ao *espaço de nomes* do pacote por meio de prefixos. Ou seja, se o pacote `foo` contém a *função* `bar()`, então podemos chamar a *função* como `foo.bar()` por qualquer um dos métodos e como `bar()` somente se importamos `foo` por meio do `using`.

Portanto, o `using` é mais cômodo, no entanto pode trazer conflitos e não permite que uma *função* ou um *tipo* de um pacote sejam trivialmente alterados ou expandidos. No entanto, em geral não é necessário se preocupar, dado que qualquer conflito será notificado pelo compilador!

Dito isso, podemos importar alguns pacotes e começar com a introdução.
"""

# ╔═╡ 11bce208-d445-11ea-277d-d1f081a3f39f
begin
	md"""

	### Introdução: O que é o Julia?

	A linguagem [Julia](https://julialang.org/) é composta de 4 elementos centrais: **nomes**, **tipos**, **variáveis** e **valores**. Podemos apresentar a relação destes elementos em uma única frase: **Uma *variável* é um *nome* associado a um *valor* de um *tipo*.** A principal força do Julia vem do seu sistema de *tipos*.

	Afinal, o que é um *tipo*?

	Bem, *tipos* são formatos para armazenar informação de modo significativo. Para a máquina, não seria possível decidir que *tipo* de *valor* irá assumir uma sequência de *bits* antes de tê-la percorrido por inteiro. *Tipos* resolvem isso reservando alguns *bits* para indicar à máquina quais são os resultados possíveis que uma sequência de *bits* produzirá. Os **tipos primitivos** mais comuns no Julia são `Bool`, `Char`, `String`, `Int64` e `Float64`. Outros *tipos primitivos* podem ser vistos [aqui](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Integers-and-Floating-Point-Numbers-1).

	O sistema de tipagem do Julia não contém somente *tipos primitivos*. Na realidade, há uma [árvore de tipos](https://docs.julialang.org/en/v1/manual/types/index.html) em que os *tipos folhas*, como os *tipos primitivos*, são ditos **tipos concretos**. Todos os demais *tipos* são ditos **tipos abstratos**. Cada *tipo* possui um único **supertipo** e somente *tipos* abstratos possuem **subtipos**. A única exceção à regra é o *tipo* `Any` que não possui um *supertipo*, sendo o *tipo* de maior hierarquia na linguagem. Como decorrência desta estrutura, dois *valores* de *tipos* distintos sempre podem encontrar o menor *supertipo* comum dentro da árvore. Além disso, **os *valores* são sempre de um *tipo concreto*.**

	Apesar do compilador ser capaz de inferir os *tipos* dos *valores* de cada variável, como em qualquer outra linguagem de tipagem dinâmica, o Julia opera de maneira mais eficiente quando utilizamos uma tipagem mais forte. Podemos fazer isso utilizando basicamente três instrumentos: a [função](https://docs.julialang.org/en/v1/base/base/#Core.typeof) `typeof()`, o [operador](https://docs.julialang.org/en/v1/manual/types/#Type-Declarations-1) `::` e o [construtor](https://docs.julialang.org/en/v1/base/base/#struct) `struct()`.

	A linguagem também utiliza um sistema de *despacho múltiplo*, na qual uma mesma *função* pode aplicar diferentes procedimentos a depender do *tipo* dos *valores* de entrada (*inputs*). O compilador é capaz de otimizar o código usando procedimentos mais específicos para os *tipos* de *valores* usados, levando a um ganho significativo de desempenho e melhora nas práticas de programação.
	"""
end

# ╔═╡ 4e8cda56-d37f-11ea-2044-f97c6be953de
begin
	md"""
	Agora que vimos um pouco sobre *nomes*, *tipos*, *variáveis* e *valores*, podemos explorar o último elemento básico: **funções**. O Julia é uma *linguagem funcional*, ou seja, a *função* é unidade básica de organização do código e a melhora do desempenho está intimamente relacionada com o quanto conseguimos discretizar o código em *funções* e evitar declarações de *variáveis* globais.

	Como vimos, pela característica de *despacho múltiplo*, uma *função* pode aceitar *valores* de entrada de diferentes *tipos* e usar um bloco de código específico para cada um. Esses blocos são chamados de *métodos* da *função*. Podemos ver os *métodos* por meio da [função](https://docs.julialang.org/en/v1/base/base/#Base.methods) `methods()`. Além disso, é possível definir e adicionar *métodos* às *funções* já existentes, bastando usar a sintaxe usual para a construção de *funções*:

	* Forma abreviada

	```
	function nomedafunção(x::<tipo>, y::<tipo> = <valor padrão>) = ...
	```

	* Forma completa

	```
	function nomedafunção(x::<tipo>, y::<tipo> = <valor padrão>)
		<comando_1>
		...
		<comando_n>
	end
	```

	Um novo *método* será adicionado à *função* caso não exista algum *método* para os *valores* dos *tipos* especificados. Além disso, *funções* podem ser *frutíferas* ou *infrutívera*, a depender se apresentarem ou não um *valor* de retorno por meio do [operador](https://docs.julialang.org/en/v1/base/base/#return) `return`.

	Vamos experimentar criar uma *função* com alguns *métodos*.
	"""
end

# ╔═╡ d1257c48-d500-11ea-2bcc-8925c4e468f6
md"""
Vejamos agora algumas operações lógicas e de comparação.

Note que para os *valores* do *tipo* `Bool`, Julia usa as palavras-chave `true` e `false`. Estas palavras são reservadas pelo Julia.

* Os operadores lógicos são `||`, `&&` e `!`, em ordem crescente de hierarquia.

* Os operadores de comparação são `>`, `≥` (`\geq <Tab>`) ou `>=`, `<`, `≤` (`\leq <Tab>`) ou `<=`, e `==`. Como exceção, existe o operador `≡` (`\equiv <Tab>`) ou `===` que irá checar se dois *valores* são iguais e possuem o mesmo *tipo*.

* Existem dois operadores `&` e `|` que são *bitwise AND* e *bitwise OR*. Estes operadores funcionam **somente** entre pares de números inteiros expressos em dígitos binários.

Se uma sentença conter operações aritméticas, lógicas e de comparação, então a ordem de avaliação será: aritmética, comparação e então lógica.

A regra de resolução de conflitos de ordem de operações segue o posto acima, com preferência para avaliar primeiro o que está dentro de parênteses e, ocorrendo empate de hierarquia de operações, realizando primeiro a operação mais à esquerda e indo em ordem até a operação mais à direita.
"""

# ╔═╡ 43ad0664-d3a3-11ea-2166-879e9399ead3
#Formas de se criar o mesmo vetor 3 × 1
[3, 1, 1]

# ╔═╡ 4bc35814-d386-11ea-22e1-e1e19f5a34d8
#Formas de se criar o mesmo vetor 3 × 1
[3; 1 ;1]

# ╔═╡ 37601e7e-d393-11ea-3abe-ef3da3e63de2
#Note no entanto que esse é um vetor 1 × 3
[3 1 1]

# ╔═╡ 490a365c-d3ee-11ea-3ac7-65cdebf1ac8d
#Criando um arranjo de arranjos, em que cada arranjo é uma linha de um vetor-coluna
[[1, 2, 3], [1 2 3]]

# ╔═╡ 3cd0dd4a-d407-11ea-328a-1dd685dc6ea2
md"""
Arranjos podem ter índices arbitrários ou que não representem o número e ordem dos seus elementos. Mais detalhes [aqui](https://docs.julialang.org/en/v1/devdocs/offset-arrays/#man-custom-indices-1).
"""

# ╔═╡ 4b7a5266-d3fd-11ea-3493-2d952e1fb5ba
md"""

Uma peculiaridade do Pluto.jl é que a atualização automática só ocorre quando alguma variável *imutável* é alterada. Logo, se `x = 2` e alteramos o programa para que `x = 3`, o Pluto.jl atualizará todas as células que contenham referência à variável x.

No entanto, os elementos do arranjos são *mutáveis*, então se mudarmos algum elemento do arranjo não teremos atualização automática. Isso só ocorrerá se mudarmos a própria variável do arranjo, ou seja:

* Não irá ativar a reatividade: `x = [1, 2, 3] ; x[1] = 2 #[2, 2, 3]`.

* Irá ativar a reatividade: `x = [1, 2, 3], x = [2, 2, 3] #[2, 2, 3]`.

Lembrem-se, no entanto, que **não** é possível múltiplas atribuições de valor à mesma variável em células diferentes!
"""

# ╔═╡ d6f2fe18-d3e6-11ea-17f6-3baada805cce
md"""
* Operações comuns entre **vetores** usando o `LinearAlgebra`, supondo `x` e `y` vetores com dimensões apropriadas e `c` um escalar:

1) [Adição](https://docs.julialang.org/en/v1/base/math/#Base.:+): `x + y`, ou `+(x, y)`;

2) [Subtração](https://docs.julialang.org/en/v1/base/math/#Base.:--Tuple{Any,Any}): `x - y`, ou `-(x, y)`;

3) [Multiplicação](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.lmul!) por escalar: `c * x`, ou `*(c, x)`, ou `x .* c`, ou `lmul!(c, x)`;

4) Adição por escalar: `x + c * ones(length(x))` ou `x .+ c`;

5) [Transposição](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#Base.transpose): `x'`, ou `transpose(x)`;

6) [Produto interno](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.dot) ($x^{T} y$): `dot(x, y)`, ou `x . y` (⋅ feito com `\cdot <Tab>`);

7) [Produto](https://docs.julialang.org/en/v1/base/math/#Base.:*-Tuple{Any,Vararg{Any,N}%20where%20N}) matricial ($x y$): `x * y`;

8) [Produto cruzado](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cross) ($x \times y$): `cross(x, y)`, ou `x × y` (× feito com `\times <Tab>`, somente com vetores 3 × 1);

9) [Norma](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.norm): `norm(x)`;

10) Distância: `norm(x - y)`;

11) [Normalização](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.normalize) $\left(\frac{x}{||x||}\right)$: `normalize(x)`;

12) [Somatório](https://docs.julialang.org/en/v1/base/collections/#Base.sum) $\left(\sum^{n}_{i = 1} x_{i}\right)$: `sum(x)`;

13) [Média](https://docs.julialang.org/en/v1/stdlib/Statistics/#Statistics.mean): `mean(x)`, no entanto, neste caso, a função estão no [pacote](https://docs.julialang.org/en/v1/stdlib/Statistics/) básico `Statistics`.

* Operações comuns entre **matrizes** usando o `LinearAlgebra`, supondo `A` e `B` matrizes com dimensões apropriadas, `x` vetor com dimensões apropriadas e `c` um escalar:

1) [Adição](https://docs.julialang.org/en/v1/base/math/#Base.:+): `A + B`, ou `+(A, B)`;

2) [Subtração](https://docs.julialang.org/en/v1/base/math/#Base.:--Tuple{Any,Any}): `A - B`, ou `-(A, B)`;

3) [Transposição](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#Base.transpose): `A'`, ou `transpose(A)`;

4) [Inversão](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#Base.inv-Tuple{AbstractArray{T,2}%20where%20T}): `inv(A)`;

5) [Pseudo-inversa de Moore-Penrose](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.pinv): `pinv(A)`;

6) [Multiplicação](https://docs.julialang.org/en/v1/base/math/#Base.:*-Tuple{Any,Vararg{Any,N}%20where%20N}) por escalar: `c * A`, ou `*(c, A)`;

7) [Multiplicação](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.mul!) por vetor ($Ax$): `A * x`, ou `*(A, x)`, ou `mul!(Y, A, x)`;

8) [Multiplicação](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.mul!) por matriz ($AB$): `A * B`, ou `*(A, B)`, ou `mul!(Y, A, B)`;

9) [Exponenciação](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.:^-Tuple{AbstractArray{T,2}%20where%20T,Number}) ($A^{c}$): `A^c`;

10) [Determinante](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.det) $\left(|A|\right)$: `det(A)`;

11) [Log-determinante](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.logdet) $\left(ln(|A|)\right)$: `logdet(A)`;

12) [Log-determinante absoluto](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.logabsdet) $\left(ln(||A||)\right)$: `logabsdet(A)`;

13) Soma da coluna $\left(\sum_{m \in M} a_{m,n}\right)$: `sum(A, dims = 1)`, atente que o arranjo resultante é um arranjo-linha;

14) Soma da linha $\left(\sum_{n \in N} a_{m,n}\right)$: `sum(A, dims = 2)`, atente que o arranjo resultante é um arranjo-coluna;

15) [Traço](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.tr) $\left(tr(A)\right)$: `tr(A)`;

16) [Espaço Nulo](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.nullspace): `nullspace(A)`.

* Fatorações e decomposições comuns entre **matrizes** usando o `LinearAlgebra`:

1. A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.eigen) `eigen()` executa a [decomposição por autovalores](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix). Os seus atributos são `values` e `vectors`, contendo os autovalores e autovetores da matriz. Podemos acessar os atributos por meio da notação `.`: `D = eigen(A) ; D.values # [λ₁, ..., λₙ]`.

2. A [função]() `svd()` executa a [decomposição por valores singulares](https://pt.wikipedia.org/wiki/Decomposi%C3%A7%C3%A3o_em_valores_singulares). Os seus atributos são `U`, `S`, `V` e `Vt`, onde `U` é uma matriz unitária real ou complexa, `S` é um vetor com números reais não-negativos e `Vt` é uma matriz unitária real ou complexa, de modo que `A = U * Diagonal(S) * Vt`.

3. A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.lu) `lu()` executa a [decomposição LU](https://pt.wikipedia.org/wiki/Decomposi%C3%A7%C3%A3o_LU). Os seus atributos são `L`, `U`, `p` e `P`, onde `L` é a matriz triangular inferior, `U` é a matriz triangular superior, `p` é o vetor de permutação à direita e `P` é a matriz de permutação à direita.

4. A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.cholesky) `cholesky()` executa a [fatoração de Cholesky](https://pt.wikipedia.org/wiki/Fatora%C3%A7%C3%A3o_de_Cholesky). Os seus atributos são `L` e `U`, onde `U` é a transposta de `L`.

5. A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.qr) `qr()` executa a [decomposição QR](https://en.wikipedia.org/wiki/QR_decomposition). Os seus atributos são `Q` e `R`.

6. A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.schur) `shur()` executa a [decomposição de Shur](https://en.wikipedia.org/wiki/Schur_decomposition). Os seus atributos são `T`, `Z` e `values`, onde `T` é o fator (quase) triagunlar de Shur, `Z` é a matriz ortogonal dos vetores de Shur e `values` é o vetor dos autovalores.

Como o Julia utiliza dispacho múltiplo, é interessante checar a documentação para verificar o que as funções retornam e seus parâmetros. Por exemplo, é especialmente importante checar sobre a função `pinv()` pelo seu papel no método dos mínimos quadrados.

"""

# ╔═╡ 9e6b5a0e-d454-11ea-3cd8-69b687166340
#A tralha ("hashtag") permite escrever comentário no corpo do código
#Carregando os pacotes
begin
	using LinearAlgebra
	using Statistics
	using PlutoUI
	import Plots
end

# ╔═╡ 5e55de14-d452-11ea-20bd-114c97ba4df2
#Verificando o tipo do valor 1
typeof(1)

# ╔═╡ 4589e708-d453-11ea-15f4-7f4efb68b159
#Verificando o supertipo de Int64
supertype(typeof(1))

# ╔═╡ 0a5f13f0-d454-11ea-2fb9-f3b15278c9ee
#Verificando o supertipo de Signed
supertype(supertype(typeof(1)))

# ╔═╡ 244d9106-d454-11ea-2672-6d3da0d2fff6
#Verificando o supertipo de Integer
supertype(supertype(supertype(typeof(1))))

# ╔═╡ ab8b9182-d453-11ea-1ddb-89717ea927fb
#Verificando o tamanho em bytes do Int64
Base.summarysize(1)

# ╔═╡ 3a4d9dee-d453-11ea-1024-9d8b156212c7
#Verificado o tipo do valor 2.5
typeof(1 + 1.5)

# ╔═╡ 7d0cf292-d453-11ea-34d8-4bd8ff8b925c
#Verificando o supertipo de Float64
supertype(typeof(1 + 1.5))

# ╔═╡ 18ca704c-d454-11ea-1900-bf83468082d0
#Verificando o supertipo de AbstractFloat
supertype(supertype(typeof(1 + 1.5)))

# ╔═╡ cd2754ca-d453-11ea-0c1e-0507fd4a323c
#Verificando o tamanho em bytes do Float64
Base.summarysize(1 + 1.5)

# ╔═╡ c34df1f0-d455-11ea-09e3-c9982355f070
#Verificando o tipo do nothing
typeof(nothing)

# ╔═╡ 499979dc-d456-11ea-1c5c-27988e76b40a
#Verificando o tipo do nothing
supertype(typeof(nothing))

# ╔═╡ 3c66adf2-d456-11ea-2503-89147b7afded
#Verificando o tipo de um arranjo [1, 2, 3]
typeof([1, 2, 3])

# ╔═╡ cc640950-d455-11ea-1700-79b264b0bb2c
#Verificando o supertipo de Array{Int64, 1}
supertype(typeof([1, 2, 3]))

# ╔═╡ 8a721040-d456-11ea-3625-63fcaaf109d2
#Verificando o supertipo de DenseArray{Int64, 1}
supertype(supertype(typeof([1, 2, 3])))

# ╔═╡ 98e5f79a-d456-11ea-036f-6148262b7d68
#Verificando o tamanho em bytes da Array{Int64, 1}
Base.summarysize([1, 2, 3])

# ╔═╡ ba514512-d378-11ea-2baa-9d3ed72cec53
begin

	a = 10.0

	md"""

	### 1. Operações básicas

	O *nome* de uma *variável* pode conter qualquer caracter [Unicode](https://en.wikipedia.org/wiki/List_of_Unicode_characters), no entanto não pode começar com um dígito. Por convenção, *variáveis* e *funções* começam com letras minúsculas, enquanto *tipos* começam com letras maiúsculas. Os seguintes *nomes* são reservados e não podem ser usados como *nome* de variáveis: `abstract type`, `baremodule`, `begin`, `break`, `catch`, `const`, `continue`, `do`, `else`, `elseif`, `end`, `export`, `false`, `finally`, `for`, `function`, `global`, `if`, `import`, `importall`, `in`, `let`, `local`, `macro`, `module`, `mutable struct`, `primitive type`, `quote`, `return`, `true`, `try`, `using`, `struct`, `where`, `while`. Não é necessário decorar essa lista, dado que, na maior parte dos IDEs, esses *nomes* serão coloridos de modo a diferenciá-los dos demais.

	Um caracter Unicode que não tenha um comando usual do teclado pode ser inserido por meio de [abreviações](https://docs.julialang.org/en/v1/manual/unicode-input/#) semelhantes a do LaTeX seguidas de `Tab`. Por exemplo: a letra grega gama maiúsculo, `Γ`, é obtida digitando `\Gamma`, seguido da tecla `Tab`. Mais interessante do que isso, a tartaruga 🐢 é obtida digitando `\:turtle:`, seguido da tecla `Tab`.

	A atribuição do *valor* à *variável* se dá por meio do operador `=`, onde o *nome* da *variável* estará sempre à esquerda da igualdade. Abaixo, na célula de código, vamos atribuir à *variável* de *nome* `a` o *valor* `10.0` do *tipo* `Float64`.

	> a = 10.0
	"""
end

# ╔═╡ e1125f26-d37b-11ea-1eff-816c70048c72
#Verificando o valor da variável a
a

# ╔═╡ 4aa9ef08-d390-11ea-2d2d-a1538323a05d
#Verificando o tipo da variável a
typeof(a)

# ╔═╡ 3edc3ff6-d461-11ea-0bbf-8788f190f750
#Criando uma função com vários métodos
begin
	function função_adicionar(x::Number, y::Number)
		return(x+y)
	end

	function função_adicionar(x::Number, y::String)
		dic = Dict("um"=> 1, "dois" => 2, "três" => 3, "quatro" => 4, "cinco" =>5, "seis"=> 6, "sete" => 7, "oito" => 8, "nove" => 9)
		return(x + dic[lowercase(y)])
	end

	function função_adicionar(x::String, y::Number)
		dic = Dict("um"=> 1, "dois" => 2, "três" => 3, "quatro" => 4, "cinco" =>5, "seis"=> 6, "sete" => 7, "oito" => 8, "nove" => 9)
		return(y + dic[lowercase(x)])
	end

	function função_adicionar(x::String, y::String)
		dic = Dict("um" => 1, "dois" => 2, "três" => 3, "quatro" => 4, "cinco" =>5,"seis"=> 6, "sete" => 7, "oito" => 8, "nove" => 9)
		return(dic[lowercase(x)] + dic[lowercase(y)])
	end

	function função_adicionar(x::Char, y::String)
		dic = Dict("um"=> 1, "dois" => 2, "três" => 3, "quatro" => 4, "cinco" =>5, "seis"=> 6, "sete" => 7, "oito" => 8, "nove" => 9)
		return(parse(Int64,x) + dic[lowercase(y)])
	end

	function função_adicionar(x::String, y::Char)
		dic = Dict("um"=> 1, "dois" => 2, "três" => 3, "quatro" => 4, "cinco" =>5, "seis"=> 6, "sete" => 7, "oito" => 8, "nove" => 9)
		return(parse(Int64,y) + dic[lowercase(x)])
	end

	function função_adicionar(x::Number, y::Char)
		return(x + parse(Int64,y))
	end

	function função_adicionar(x::Char, y::Number)
		return(y + parse(Int64,x))
	end

	function função_adicionar(x::Char, y::Char)
		return(parse(Int64,x) + parse(Int64,y))
	end
end

# ╔═╡ 4b2e3e5e-d465-11ea-37b0-87f8a3d719b6
#Vendo os métodos da função_adicionar()
methods(função_adicionar)

# ╔═╡ 4f8ce150-d466-11ea-3d3b-612a86a20dd5
#Usando os vários métodos da função
função_adicionar('5',5) + função_adicionar('5',"cinco") + função_adicionar("cinco",5)

# ╔═╡ 7283abe8-d467-11ea-3097-3b063f8a3ca3
begin
	md"""
	Vejamos algumas [operações e funções](https://docs.julialang.org/en/v1/manual/mathematical-operations/) básicas com escalares. Para isso vamos usar uma nova variável.

	**Selecione o *valor* de 🐢 na barra de correr abaixo**:

	$(@bind 🐢 Slider(-25.0:0.01:25.0))

	Observe que em todas as células abaixo os *valores* de 🐢 são atualizados, com estado inicial igual a `-25`. Além disso, vemos que 🐢 é do *tipo* `Int64` no seu estado inicial. No entanto, ao alterar o *valor*, o *tipo* pode mudar para `Float64`.

	Se levarmos a barra ao outro extremo, no *valor* igual a `25`, veremos que ocorrerá uma [promoção](https://docs.julialang.org/en/v1/manual/conversion-and-promotion/) automática em que o *tipo* do *valor* de 🐢 será convertido para `Int64`. Isso não é o usual. Geralmente, somente ocorrerá o contrário, dado que qualquer número inteiro tem uma representação em ponto flutuante, mas nem todo ponto flutuante terá uma representação em número inteiro.

	Logo abaixo temos o gráfico de uma função

	\begin{equation}
	f(x) = sen\left( \frac{x}{b} \right) \text{ ,}
	\end{equation}

	onde $$x \in [0,30]$$ e $$b$$ como a nossa tartaruga.
	"""
end

# ╔═╡ 71e90f4a-d467-11ea-14ea-8ddad3ed52b3
#Verificando o valor da variável b
md"""

**Valor de 🐢**: $🐢

**Tipo da variável 🐢**: $(typeof(🐢))
"""

# ╔═╡ 94563de6-d381-11ea-08f4-fbfeaf4cc4d5
#Vamos brincar um pouco com uma visualização
#As funções trigonométricas como sin(), cos() e tan() usam como entrada radianos!
begin
	println("Vamos chamar o $🐢")
	f(x, b = 🐢) = sin(x / b)
	Plots.plot(f, 0:0.01:30, legend = false)
end

# ╔═╡ 25960f18-d380-11ea-1a4b-f9ba52963806
#Soma
a + 🐢

# ╔═╡ 2c61bc66-d380-11ea-1a6e-5dc74f8e39c6
#Subtração
a - 🐢

# ╔═╡ 328fa512-d380-11ea-0509-edff01e88423
#Multiplicação
a * 🐢

# ╔═╡ 36e87f44-d380-11ea-0632-6fafa86b2ba2
#Divisão com a no numerador e b no denominador
a / 🐢

# ╔═╡ 3fa8d246-d380-11ea-2195-bd3862da65f0
#Divisão com b no numerador e a no denominador
a \ 🐢

# ╔═╡ 5a3e856a-d380-11ea-25e5-f9bcb219d70e
#Potenciação
a^🐢

# ╔═╡ 6b1f3104-d380-11ea-3f61-b3aa1408c4d9
#Divisão inteira de a por b: div ou ÷ (\div <Tab>)
div(a, 🐢)

# ╔═╡ 0ff9a65a-d381-11ea-2eb5-c3d6dc0e7cf1
#Mesmo acima
a ÷ 🐢

# ╔═╡ 9fc3a96c-d380-11ea-2d22-c3b0f3800736
#Resto da divisão inteira de a por b: mod, rem ou %
mod(a, 🐢)

# ╔═╡ c2a6e6f6-d380-11ea-2848-e765432972dd
#Mesmo acima
a % 🐢

# ╔═╡ d91246be-d381-11ea-22d3-bb8b8ce382c8
#Representação como racional de a / b
rationalize(a / 🐢)

# ╔═╡ 982d5204-d3e4-11ea-3feb-6d17061562d5
#Checar uma igualdade com ==
#Neste caso, verifica se referem-se aos mesmos valores
a == 🐢

# ╔═╡ f2d99db2-d409-11ea-363d-1d58daadc79a
#Checar uma equivalência com ≡ (\equiv <Tab>) ou ===
a ≡ 10

# ╔═╡ f898ebd8-d409-11ea-39a8-c5bf1392e99b
#Checar uma igualdade com ==
a == 10

# ╔═╡ 73dca420-d3e4-11ea-2437-e5e894e94587
#Checar uma desigualdade com ≠ (\ne <Tab>)
a ≠ 🐢

# ╔═╡ 3c9b3894-d382-11ea-32cc-1f075f8fff54
#Comparação se a é maior ou igual (\geq <Tab>) ao b
a ≥ 🐢

# ╔═╡ 0af09224-d389-11ea-115e-4ba278ee876a
#Comparação dupla usando menor ou igual (\leq <Tab>)
0 ≤ a ≤ 🐢

# ╔═╡ 573f4b9c-d38a-11ea-137e-cfbf238c3d14
begin

	c = Vector([10.0; 15.0; 20.0])

	md"""

	### 2. Operações vetoriais

	O [construtor](https://docs.julialang.org/en/v1/base/arrays/#Core.Array) `Array()`, junto ao operador relacionado `[]`, criará um *arranjo-coluna*. O delimitador padrão entre os elementos de arranjos unidimensionais é a vírgula (`,`), no entanto, em arranjos com mais dimensões, o delimitador para concatenar elementos na horizontal é o espaço em branco, enquanto que para concatenar elementos na vertical é o ponto-e-vírgula (`;`). Para criar arranjos $n$-dimensionais, pode ser necessário usar o parâmetro de dimensão no comando `Array{<Tipo>, <n>}([<dados>], (<nᵒ de elementos d₁>, ..., <nᵒ de elementos dₙ>))`.

	O arranjo é uma estrutura de dados *mutável* e *indexável* que pode ser fatiado pela notação de colchete: `x = [1, 2, 3, 4]; x[2:end] #[2,3,4]`. Note que `end` é um argumento de índice que indica o último elemento de um objeto. Além disso, é possível criar peridiciocidade no fatiamento por meio de um terceiro parâmetro: `x[<começo>:<passo>:<fim>]`. Também podemos usar o [operador](https://docs.julialang.org/en/v1/base/math/#Base.::) `:` para criar uma amplitude (*range*).

	Vetores também podem ser concatenados verticalmente com a [função](https://docs.julialang.org/en/v1/base/arrays/#Base.vcat) `vcat(<vetor1>, <vetor2>)` ou horizontalmente com a [função](https://docs.julialang.org/en/v1/base/arrays/#Base.hcat) `hcat(<vetor1>, <vetor2>)`. As funções exigem que os vetores tenham dimensões compatíveis ou retornarão erro.

	Caso os elementos sejam iteráveis ou uma amplitude (*range*), então o delimitador usado mudará o resultado final de maneira drástica.

	Espaçoem branco implicará que cada coluna será ocupada pelos elementos da amplitude, ou seja, a primeira coluna será dos termos da primeira amplitude, a segunda coluna será dos termos da segunda amplitude *etc*.;

	1. Vírgula(`,`) implicará que os próprios objetos de amplitude serão os elementos em um vetor-coluna;

	2. Ponto-e-vírgula(`;`) implicará que haverá uma única coluna com os elementos da amplitude, serão colocados um abaixo do outro seguindo a ordem.

	Uma mistura destes separadores pode ser usada para gerar algum efeito desejado. Por exemplo, `m6 = [1:3 4:6; 7:9 10:12]` gerará uma matriz com primeira coluna composta pelo vetor `[1, 2, 3, 7, 8, 9]` e a segunda coluna composta pelo vetor `[4, 5, 6, 10, 11, 12]`.

	Julia é uma linguagem *column-major*, ou seja, as estruturas de dados $n$-dimensionais são acessadas "verticalmente", seguindo as colunas. Assim, se `m1` é uma matriz $M$ 2 × 2, então `m1[1]` é o elemento $m_{1,1}$, `m1[2]` é o elemento $m_{2,1}$, `m1[3]` é o elemento $m_{1,2}$ e `m1[4]` é o elemento $m_{2,2}$ da matriz.

	Abaixo, na célula de código, vamos atribuir à variável `c` o arranjo 3 × 1 com elementos `[10.0; 10.0; 10.0]` do tipo `Array{Float64, 1}`. De maneira geral, em operações matemáticas, o Julia tratará arranjos como vetores. Para declarar explicitamente um arranjo é do *tipo* vetor, podemos usar o [construtor](https://docs.julialang.org/en/v1/base/arrays/#Base.Vector) `Vector()`. No caso de arranjos bidimensionais, podemos fazer o mesmo para o *tipo* matriz por meio do [construtor](https://docs.julialang.org/en/v1/base/arrays/#Base.Matrix) `Matrix()`.
	"""
end

# ╔═╡ dea4dd2c-d3e2-11ea-00ff-8d164680707f
#Verificando o valor da variável c
c

# ╔═╡ ac7996ba-d3ff-11ea-1c9e-754c7b9af8af
#Verificando o tipo da variável c
typeof(c)

# ╔═╡ abb5b982-d408-11ea-070a-3d2163c19406
#Note que se criarmos um vetor d com os mesmos valores, teremos outro objeto
d = Vector([10.0, 15.0, 20.0])

# ╔═╡ c01277a6-d408-11ea-0dd5-317c0ee3cf6e
#Verificando se c e d são equivalentes com ≡ (\equiv <Tab>) ou ===
c ≡ d

# ╔═╡ d71ce17a-d408-11ea-3e1d-39bc64921448
#Verificando se c e d são idênticos com ==
c == d

# ╔═╡ 4852ed98-d40a-11ea-326b-5bc4cb41ef1c
#Note que agora somente estamos dando um outro nome ao mesmo objeto
e = c

# ╔═╡ 6b204b2a-d40a-11ea-12c2-6bb7734bee65
#Verificando se c e e são equivalentes com ≡ (\equiv <Tab>) ou ===
c ≡ e

# ╔═╡ 70a6b630-d40a-11ea-2598-4752fbd46e78
#Verificando se c e e são idênticos com ==
c == e

# ╔═╡ e25ec324-d38a-11ea-24ce-f7113692e2db
#Formas de se criar o mesmo vetor 3 × 1
Array([3, 1, 1])

# ╔═╡ 9779777c-b329-11ea-3b3b-47ed2b4cf59c
#Criando uma matriz identidade A com dimensões 3x3
A = Matrix([2 0 0; 0 1 0; 0 0 1])

# ╔═╡ 4290aa8c-d3a4-11ea-1177-ff69a45735b0
#Verificando o tipo de A
typeof(A)

# ╔═╡ 698ac79a-d3a8-11ea-0efb-9da5e7268b34
#Primeira coluna de A
A[:, 1]

# ╔═╡ c8a7b4b8-d3a8-11ea-1cd7-c7dfa1be6c49
#A primeira linha de A
A[1, :]

# ╔═╡ 2e1bad0e-d3a9-11ea-3b66-2f70af73242f
#Criando arranjo bimensional com os elementos de iteradores como colunas
[2:4 7:9 12:14]

# ╔═╡ 33d37e94-d3ee-11ea-2ae1-958b308a8113
#Criando um arranjo com iteradores como elementos
[2:4, 3:5, 4:6]

# ╔═╡ 345fbf30-d3ee-11ea-2a55-e1ad67b68838
#Criando um arranjo com os elementos de iteradores como uma única coluna
[2:4; 3:5; 4:6]

# ╔═╡ 5213eae6-d3e2-11ea-04bd-b1b48fe6f533
begin
	B = Matrix{Float64}(undef, 2, 2)

	md"""
	O pacote `LinearAlgebra` do Julia básico contém alguns construtores, operadores e funções especializados para Álgebra Linear, como, por exemplo, `I` para gerar uma matriz identidade de dimensão arbitrária. A documentação do pacote pode ser encontrada [aqui](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) e é recomendado que qualquer pessoa com interesse em computação matemática se familiarize com este pacote.

	Algumas funções especialmente importantes para a manipulação de arranjos são:

	* A [função](https://docs.julialang.org/en/v1/base/arrays/#Base.ndims) `ndims()` retorna o número de dimensões de um arranjo n-dimensional.

	* A [função](https://docs.julialang.org/en/v1/base/arrays/#Base.length-Tuple{AbstractArray}) `length()` retorna o número de elementos de um arranjo.

	* A [função](https://docs.julialang.org/en/v1/base/arrays/#Base.size) `size()` retorna o comprimento de cada dimensão de um arranjo.

	* A [função](https://docs.julialang.org/en/v1/base/arrays/#Base.zeros) `zeros()` retorna um arranjo com todas as entradas nulas com a dimensão desejada.

	* A [função](https://docs.julialang.org/en/v1/base/arrays/#Base.ones) `ones()` retorna um arranjo com todas as entradas iguais a 1 com a dimensão desejada.

	* A [função](https://docs.julialang.org/en/v1/base/arrays/#Base.reshape) `reshape(<objeto>, (<dimensões>))` permite alterar as dimensões de um arranjo. Deste modo, se `v2 = [1 2 3]` é um arranjo 1x3, então podemos transformá-lo em um arranjo 3x1 com o comando `reshape(v2, (3, 1))`.

	* O [construtor](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.Diagonal) `Diagonal()` construirá uma matriz diagonal n × n se for passado um arranjo n × 1 como argumento.

	* A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.diag) `diag()` extrairá a diagonal de uma matriz, retornando-a como um arranjo.

	É possível pré-alocar valores arbitrários em um arranjo para que a estrutura tenha as dimensões correta por meio de `undef`. Por exemplo, podemos criar uma matriz 2 × 2 de elementos do tipo `Float64` com `Matrix{Float64}(undef, 2, 2)`. As [constantes](https://docs.julialang.org/en/v1/base/arrays/#Core.undef) `undef` irão utilizar valores que estejam disponíveis na memória ou marcados para serem removidos pelo *garbage collector* como *placeholders*.

	Além disso, há um [operador](https://docs.julialang.org/en/v1/manual/mathematical-operations/#man-dot-operators-1) *dot* (`.`) para funções escalares que aplicam a função a cada elemento de um arranjo. Por exemplo: `cos.(A)` aplica a função `cos()` a todos os elementos da matriz $A$.
	"""
end

# ╔═╡ 3f99ff12-d3eb-11ea-11f3-5b5df7549571
#Verificando o valor da variável B
B

# ╔═╡ b62ed4c0-d3fc-11ea-2370-2fd698ef96ae
#Vamos atribuir alguns valores para os elementos de B
begin
	B[1] = 1
	B[2] = rand(0:10, 1)[1]
	B[:,2] = rand(0:10, 2)
	B
end

# ╔═╡ 6e286974-d411-11ea-3fe9-f3c3ae0dc280
#Vejamos o comprimento de B
length(B)

# ╔═╡ 71b3d632-d411-11ea-0e22-0d78174b087d
#Vejamos o tamanho de B
size(B)

# ╔═╡ 744cdbbe-d411-11ea-2c98-8fc6abc6f4ee
##Vejamos o número de dimensões de B
ndims(B)

# ╔═╡ ac7bfdf4-d40b-11ea-1fa1-3beb6ee6ca54
#Extraindo a diagonal de B
diag(B)

# ╔═╡ eaa83194-d3eb-11ea-0478-79370e80b20b
#Subtraindo 1 dos elementos da matriz B
B .- 1

# ╔═╡ aa47de40-d410-11ea-09a5-15ab5bf763f0
#Subtraindo 1 somente dos elementos fora da diagonal de B
B - ones(size(B)) + I

# ╔═╡ eeba65e8-d401-11ea-0a37-d32e0ac61cc9
#Multiplicando cada um dos elementos de B por 2
B .* 2

# ╔═╡ f7ba275a-d401-11ea-1b89-fdfd7e1c48bc
#Note que neste caso é o mesmo que simplesmente multiplicar a matriz por um escalar
B * 2

# ╔═╡ 7b85d7ba-d3eb-11ea-1dd3-3b5f8a105916
#Aplicando sin() aos elementos matriz [1 0; 0 1] multiplicados por 3π/2
sin.(3π/2 .* B)

# ╔═╡ ea4d1a52-d3f5-11ea-01b3-0f58262e6d9e
#Note que é igual ao seno da matriz somente porque a matriz é quadrada!
sin(3π/2 .* B)

# ╔═╡ 20703f8a-d3f6-11ea-03aa-f7adfb9e3654
#Aplicando sin() aos elementos matriz 3 × 2 [1 0; 0 1; 2 2] multiplicados por 3π/2
sin.(3π/2 .* [1 0; 0 1; 2 2])

# ╔═╡ 0daccb54-d3ed-11ea-109f-3d6315ab66b5
#Elevando cada elemento da matriz ao quadrado
(B .- 1).^ 2

# ╔═╡ 21fcc64a-d3ed-11ea-15dc-d35654f03f99
#Note que usualmente é diferente da matriz ao quadrado!
(B .- 1)^2

# ╔═╡ 75a9c00c-d3f6-11ea-1205-19aad70fa1d9
begin
	md"""

	Selecione abaixo os valores para atribuirmos à nossa matriz $M$. Note que inicialmente a matriz será uma matriz de elementos `Int64`, no entanto se optarmos por colocar algum número com casa decimal, imediatamente será convertida em uma matriz de elementos `Float64`. Os valores permitidos estão no intervalo $[-100, 100]$.

	**Primeira linha**: $(@bind m1 html"<input type ='number', min = '-100', max = '100', value = '1', size = '1'>")
	$(@bind m4 html"<input type ='number', min = '-100', max = '100', value = '0', size = '1'>")
	$(@bind m7 html"<input type ='number', min = '-100', max = '100', value = '0', size = '1'>")

	**Segunda linha**: $(@bind m2 html"<input type ='number', min = '-100', max = '100', value = '0', size = '1'>")
	$(@bind m5 html"<input type ='number', min = '-100', max = '100', value = '1', size = '1'>")
	$(@bind m8 html"<input type ='number', min = '-100', max = '100', value = '0', size = '1'>")

	**Terceira linha**: $(@bind m3 html"<input type ='number', min = '-100', max = '100', value = '0', size = '1'>")
	$(@bind m6 html"<input type ='number', min = '-100', max = '100', value = '0', size = '1'>")
	$(@bind m9 html"<input type ='number', min = '-100', max = '100', value = '1', size = '1'>")
	"""
end

# ╔═╡ efaf99ba-d3f8-11ea-2eab-9d3cae4f8697
#Criando uma matriz 3 × 3
begin
	M = [m1 m4 m7; m2 m5 m8; m3 m6 m9]

	md"""

	Vamos criar uma matriz interativa e realizar alguns testes sobre suas propriedades. A nossa primeira matriz será 3 × 3. As funções que usaremos serão:

	* A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.issymmetric) `issymmetric()` para testar se uma matriz é simétrica.

	* A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.isposdef) `isposdef()` para testar se uma matriz é positiva definida. Esse teste é realizado tentando aplicar a fatoração Cholesky e checando o resultado com `issuccess()`.

	* As [funções](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.istril) `istril()` e `istriu()` para testar se uma matriz é triangular inferior ou superior, respectivamente.

	* A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.isdiag) `isdiag()` para testar se uma matriz é diagonal.

	* A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.ishermitian) `ishermitian()` para testar se uma matriz é Hermitiana.

	* A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.checksquare) `checksquare()` para testar se uma matriz é quadrada, retornando a dimensão comum.

	* A [função](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.rank) `rank()` retorna o posto da matriz.

	Também podemos verificar se uma variável é um dos elementos de um arranjo por meio do [operador](https://docs.julialang.org/en/v1/base/collections/#Base.in) `∈` (`\in <Tab>`) ou `in`. Caso a ordem esteja inversa, podemos usar o operador inverso `∋` (`\ni <Tab>`). Deste modo, tomanto a variável como `a` e o arranjo como `c`: `a ∈ c` é equivalente a `c ∋ a`. A negação desta relação é dada pelo [operador](https://docs.julialang.org/en/v1/base/collections/#Base.:%E2%88%89) `∉` (`\notin <Tab>`), com inverso `∌` (`\nni <Tab>`).

	Após os testes, veremos algumas operações e decomposições de matrizes.
	"""
end

# ╔═╡ 393945e8-d3f9-11ea-1003-5f6dfe66e9e7
#Visualizando a matriz M
M

# ╔═╡ d672b758-d3fa-11ea-0741-25003cc6bb63
#Testando se M é simétrica
issymmetric(M)

# ╔═╡ e094a62e-d3fa-11ea-2f92-71c3a1bd739b
#Testando se M é positiva definida
isposdef(M)

# ╔═╡ e9db294e-d3fa-11ea-2d1e-8fb2c74f307c
#Testando se M é triangular superior
istriu(M)

# ╔═╡ f41fac98-d3fa-11ea-1438-9125dfcc7ceb
#Testando se M é diagonal
isdiag(M)

# ╔═╡ fd9b0132-d3fa-11ea-0dd7-591bbcedb9eb
#Testando se M é Hermitiana
ishermitian(M)

# ╔═╡ cad69f28-d3f9-11ea-3001-b9ae34a6d095
#Checando se a matriz é quadrada e retornando a dimensão comum
LinearAlgebra.checksquare(M)

# ╔═╡ 5a71ae82-d401-11ea-250a-f9833e649895
#Checando se a matriz é quadrada e retornando a dimensão comum
LinearAlgebra.checksquare(B)

# ╔═╡ 47fb5076-d40d-11ea-2add-e315e0284c6d
#Checando o posto da matriz M
rank(M)

# ╔═╡ 72de5950-d417-11ea-0dd7-09aed198885b
#Checando se a matriz M contém o número 1
1 ∈ M

# ╔═╡ 7cf8b9da-d417-11ea-0194-25f516dcac23
#Checando se a matriz M não contém o número 1
1 ∉ M

# ╔═╡ 867adb1e-d417-11ea-3478-f98e26ecfb77
#Note que é o mesmo que usar o operador de negação ! na expressão inteira
!(1 ∈ M)

# ╔═╡ 13e69962-d401-11ea-2876-19822449cec5
#Lembrando do valor de c
c

# ╔═╡ 1c63b770-d401-11ea-350b-6182e0a78a13
#Lembrando do valor de c
A

# ╔═╡ 86c24c9a-d400-11ea-2291-3f47ae05b0d7
#Realizando o produto A * x
A * c

# ╔═╡ 8aa49656-d400-11ea-2805-19d2558de524
#Criando um vetor x_t com dimensões 1x3, transposto x com a função transpose()
c_t = transpose(c)

# ╔═╡ 8e79a7f8-d400-11ea-2e6d-772aceae30c6
#Mesma coisa que utilizar o operador ' para transpor um arranjo, vetor ou matriz
c'

# ╔═╡ 963e8e68-d400-11ea-04f0-59a9f4991015
#Verificando o tipo de x_t
typeof(c_t)

# ╔═╡ a3ce9d7a-d400-11ea-2f6c-a9d649b816b6
#Transposta de A
M_t = M'

# ╔═╡ a89e8a40-d400-11ea-377d-75be788e91cf
#Realizando o produto x_t * A
c_t * M

# ╔═╡ abd800a6-d400-11ea-219f-99d3f4fd501f
#Vamos verificar se é o caso de (Ac)' == c_t A_t
(M * c)' == c_t * M_t

# ╔═╡ c92d4c7c-d3a8-11ea-30a8-fb9f07f3bf83
#A constante I utiliza matriz identidade de tamanho arbitrário
c_t * I

# ╔═╡ 2e753018-d3a9-11ea-3ef5-054db5875158
#O operador ⋅ (\cdot ,<Tab>) realiza o produto interno
c_t ⋅ c

# ╔═╡ 2e60c100-d3a9-11ea-2a0d-e182a55603ca
#Dependendo do contexto, o mesmo será feito com o operador *
c_t * c

# ╔═╡ 2e4b0f2c-d3a9-11ea-131f-6f0ed0af20c1
#Nem sempre será o caso! A multiplicação de um vetor 3 × 1 por um 1 × 3 gera uma matriz
c * c_t

# ╔═╡ 2e33f544-d3a9-11ea-21d4-317b43764db9
#E o produto escalar destes mesmo vetores ainda gera o resultado de xt * x
c ⋅ c_t

# ╔═╡ 22b97952-d401-11ea-1971-57b730a767fd
#Inversa do produto (M_t M)
inv((M_t * M))

# ╔═╡ 2dff201c-d3a9-11ea-29fa-d3ffa50117a3
#Usando a fatoração Cholesky se a matriz M for positiva definida
begin
	if isposdef(M)
		cholesky(M)
	end
end

# ╔═╡ 2dc1df34-d3a9-11ea-1bd4-dd0b5a9b8fb7
#Usando a decomposição espectral. Note que os autovalores retornam como vetor 3 × 1
eigen(M)

# ╔═╡ Cell order:
# ╟─7548fe76-b328-11ea-1dff-b16fc38e0cc5
# ╠═9e6b5a0e-d454-11ea-3cd8-69b687166340
# ╟─11bce208-d445-11ea-277d-d1f081a3f39f
# ╠═5e55de14-d452-11ea-20bd-114c97ba4df2
# ╠═4589e708-d453-11ea-15f4-7f4efb68b159
# ╠═0a5f13f0-d454-11ea-2fb9-f3b15278c9ee
# ╠═244d9106-d454-11ea-2672-6d3da0d2fff6
# ╠═ab8b9182-d453-11ea-1ddb-89717ea927fb
# ╠═3a4d9dee-d453-11ea-1024-9d8b156212c7
# ╠═7d0cf292-d453-11ea-34d8-4bd8ff8b925c
# ╠═18ca704c-d454-11ea-1900-bf83468082d0
# ╠═cd2754ca-d453-11ea-0c1e-0507fd4a323c
# ╠═c34df1f0-d455-11ea-09e3-c9982355f070
# ╠═499979dc-d456-11ea-1c5c-27988e76b40a
# ╠═3c66adf2-d456-11ea-2503-89147b7afded
# ╠═cc640950-d455-11ea-1700-79b264b0bb2c
# ╠═8a721040-d456-11ea-3625-63fcaaf109d2
# ╠═98e5f79a-d456-11ea-036f-6148262b7d68
# ╟─ba514512-d378-11ea-2baa-9d3ed72cec53
# ╠═e1125f26-d37b-11ea-1eff-816c70048c72
# ╠═4aa9ef08-d390-11ea-2d2d-a1538323a05d
# ╟─4e8cda56-d37f-11ea-2044-f97c6be953de
# ╠═3edc3ff6-d461-11ea-0bbf-8788f190f750
# ╠═4b2e3e5e-d465-11ea-37b0-87f8a3d719b6
# ╠═4f8ce150-d466-11ea-3d3b-612a86a20dd5
# ╟─71e90f4a-d467-11ea-14ea-8ddad3ed52b3
# ╟─7283abe8-d467-11ea-3097-3b063f8a3ca3
# ╟─94563de6-d381-11ea-08f4-fbfeaf4cc4d5
# ╠═25960f18-d380-11ea-1a4b-f9ba52963806
# ╠═2c61bc66-d380-11ea-1a6e-5dc74f8e39c6
# ╠═328fa512-d380-11ea-0509-edff01e88423
# ╠═36e87f44-d380-11ea-0632-6fafa86b2ba2
# ╠═3fa8d246-d380-11ea-2195-bd3862da65f0
# ╠═5a3e856a-d380-11ea-25e5-f9bcb219d70e
# ╠═6b1f3104-d380-11ea-3f61-b3aa1408c4d9
# ╠═0ff9a65a-d381-11ea-2eb5-c3d6dc0e7cf1
# ╠═9fc3a96c-d380-11ea-2d22-c3b0f3800736
# ╠═c2a6e6f6-d380-11ea-2848-e765432972dd
# ╠═d91246be-d381-11ea-22d3-bb8b8ce382c8
# ╟─d1257c48-d500-11ea-2bcc-8925c4e468f6
# ╠═982d5204-d3e4-11ea-3feb-6d17061562d5
# ╠═f2d99db2-d409-11ea-363d-1d58daadc79a
# ╠═f898ebd8-d409-11ea-39a8-c5bf1392e99b
# ╠═73dca420-d3e4-11ea-2437-e5e894e94587
# ╠═3c9b3894-d382-11ea-32cc-1f075f8fff54
# ╠═0af09224-d389-11ea-115e-4ba278ee876a
# ╟─573f4b9c-d38a-11ea-137e-cfbf238c3d14
# ╠═dea4dd2c-d3e2-11ea-00ff-8d164680707f
# ╠═ac7996ba-d3ff-11ea-1c9e-754c7b9af8af
# ╠═abb5b982-d408-11ea-070a-3d2163c19406
# ╠═c01277a6-d408-11ea-0dd5-317c0ee3cf6e
# ╠═d71ce17a-d408-11ea-3e1d-39bc64921448
# ╠═4852ed98-d40a-11ea-326b-5bc4cb41ef1c
# ╠═6b204b2a-d40a-11ea-12c2-6bb7734bee65
# ╠═70a6b630-d40a-11ea-2598-4752fbd46e78
# ╠═e25ec324-d38a-11ea-24ce-f7113692e2db
# ╠═43ad0664-d3a3-11ea-2166-879e9399ead3
# ╠═4bc35814-d386-11ea-22e1-e1e19f5a34d8
# ╠═37601e7e-d393-11ea-3abe-ef3da3e63de2
# ╠═9779777c-b329-11ea-3b3b-47ed2b4cf59c
# ╠═4290aa8c-d3a4-11ea-1177-ff69a45735b0
# ╠═698ac79a-d3a8-11ea-0efb-9da5e7268b34
# ╠═c8a7b4b8-d3a8-11ea-1cd7-c7dfa1be6c49
# ╠═2e1bad0e-d3a9-11ea-3b66-2f70af73242f
# ╠═33d37e94-d3ee-11ea-2ae1-958b308a8113
# ╠═345fbf30-d3ee-11ea-2a55-e1ad67b68838
# ╠═490a365c-d3ee-11ea-3ac7-65cdebf1ac8d
# ╟─3cd0dd4a-d407-11ea-328a-1dd685dc6ea2
# ╟─5213eae6-d3e2-11ea-04bd-b1b48fe6f533
# ╠═3f99ff12-d3eb-11ea-11f3-5b5df7549571
# ╟─4b7a5266-d3fd-11ea-3493-2d952e1fb5ba
# ╠═b62ed4c0-d3fc-11ea-2370-2fd698ef96ae
# ╠═6e286974-d411-11ea-3fe9-f3c3ae0dc280
# ╠═71b3d632-d411-11ea-0e22-0d78174b087d
# ╠═744cdbbe-d411-11ea-2c98-8fc6abc6f4ee
# ╠═ac7bfdf4-d40b-11ea-1fa1-3beb6ee6ca54
# ╠═eaa83194-d3eb-11ea-0478-79370e80b20b
# ╠═aa47de40-d410-11ea-09a5-15ab5bf763f0
# ╠═eeba65e8-d401-11ea-0a37-d32e0ac61cc9
# ╠═f7ba275a-d401-11ea-1b89-fdfd7e1c48bc
# ╠═7b85d7ba-d3eb-11ea-1dd3-3b5f8a105916
# ╠═ea4d1a52-d3f5-11ea-01b3-0f58262e6d9e
# ╠═20703f8a-d3f6-11ea-03aa-f7adfb9e3654
# ╠═0daccb54-d3ed-11ea-109f-3d6315ab66b5
# ╠═21fcc64a-d3ed-11ea-15dc-d35654f03f99
# ╟─efaf99ba-d3f8-11ea-2eab-9d3cae4f8697
# ╠═393945e8-d3f9-11ea-1003-5f6dfe66e9e7
# ╟─75a9c00c-d3f6-11ea-1205-19aad70fa1d9
# ╠═d672b758-d3fa-11ea-0741-25003cc6bb63
# ╠═e094a62e-d3fa-11ea-2f92-71c3a1bd739b
# ╠═e9db294e-d3fa-11ea-2d1e-8fb2c74f307c
# ╠═f41fac98-d3fa-11ea-1438-9125dfcc7ceb
# ╠═fd9b0132-d3fa-11ea-0dd7-591bbcedb9eb
# ╠═cad69f28-d3f9-11ea-3001-b9ae34a6d095
# ╠═5a71ae82-d401-11ea-250a-f9833e649895
# ╠═47fb5076-d40d-11ea-2add-e315e0284c6d
# ╠═72de5950-d417-11ea-0dd7-09aed198885b
# ╠═7cf8b9da-d417-11ea-0194-25f516dcac23
# ╠═867adb1e-d417-11ea-3478-f98e26ecfb77
# ╟─d6f2fe18-d3e6-11ea-17f6-3baada805cce
# ╠═13e69962-d401-11ea-2876-19822449cec5
# ╠═1c63b770-d401-11ea-350b-6182e0a78a13
# ╠═86c24c9a-d400-11ea-2291-3f47ae05b0d7
# ╠═8aa49656-d400-11ea-2805-19d2558de524
# ╠═8e79a7f8-d400-11ea-2e6d-772aceae30c6
# ╠═963e8e68-d400-11ea-04f0-59a9f4991015
# ╠═a3ce9d7a-d400-11ea-2f6c-a9d649b816b6
# ╠═a89e8a40-d400-11ea-377d-75be788e91cf
# ╠═abd800a6-d400-11ea-219f-99d3f4fd501f
# ╠═c92d4c7c-d3a8-11ea-30a8-fb9f07f3bf83
# ╠═2e753018-d3a9-11ea-3ef5-054db5875158
# ╠═2e60c100-d3a9-11ea-2a0d-e182a55603ca
# ╠═2e4b0f2c-d3a9-11ea-131f-6f0ed0af20c1
# ╠═2e33f544-d3a9-11ea-21d4-317b43764db9
# ╠═22b97952-d401-11ea-1971-57b730a767fd
# ╠═2dff201c-d3a9-11ea-29fa-d3ffa50117a3
# ╠═2dc1df34-d3a9-11ea-1bd4-dd0b5a9b8fb7

max_stars_repo_name: fausto-mpj/EduPluto

max_issues_repo_name: fausto-mpj/EduPluto

max_forks_repo_name: fausto-mpj/EduPluto


text:
function allreduce(sendbuff, recvbuff, count::Int, datatype::Int, op::Int, comm, stream::CuStream)
    @nccl(:ncclAllReduce, (Ptr{Cvoid},Ptr{Cvoid},Csize_t,Cint,Cint,Ptr{Cvoid},Ptr{Cvoid}),
        sendbuffm, recvbuff, count, datatype, op, comm, stream)
end

function bcast(buff, count::Int, datatype::Int, root::Int, comm::Comm, stream::CuStream)
    @nccl(:ncclBcast, (Ptr{Cvoid},Csize_t,Cint,Cint,Ptr{Cvoid},Ptr{Cvoid}),
        buff, count, datatype, root, comm, stream)
end

function Base.reduce(sendbuff, recvbuff, count::Int, datatype::Int, op::Int, root::Int, comm::Comm, stream::CuStream)
    @nccl(:ncclReduce, (Ptr{Cvoid},Ptr{Cvoid},Csize_t,Cint,Cint,Cint,Ptr{Cvoid},Ptr{Cvoid}),
        sendbuff, recvbuff, count, datatype, op, root, comm, stream)
end

ext: jl


text:
function numbersYield(n::Integer)
  @sync begin
    for i in 1:3
      @async begin
        for a in 1:n
          print(i)
          yield()
        end
      end
    end
  end
end

function produceNumber(number::Integer, n::Integer)
  for i = 1:n
    produce(number)
  end
end

function consumeNumber(p::Task)
  print(consume(p))
  print(" ")
end

function numbersProdCon(n::Integer)
  p1 = @task produceNumber(1,n)
  p3 = @task produceNumber(3,n)
  p2 = @task produceNumber(2,n)
  for i = 1:n
    consumeNumber(p1)
    consumeNumber(p2)
    consumeNumber(p3)
  end
end

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
"""
* `sxp` : how to create the inversion vector for sa?
   * =:positive  only positive lags
   * =:all all lags
"""
mutable struct IBD{T}
	om::ObsModel{T}
	optm::OptimModel{T}
	gx::X{T}
	sx::X{T}
	sxp::Sxparam
	sx_fix_zero_lag_flag::Bool
	attrib_inv::Symbol
	verbose::Bool
	err::DataFrames.DataFrame
	opG::LinearMaps.LinearMap{T}
	opS::LinearMaps.LinearMap{T}
	pbandpassG::P_bandpass{T}
	pbandpassS::P_bandpass{T}
	# if eigenvalue decomposition of toeplitz(s) is necessary
	toepS::Array{T,2}
end

"""
`gprecon` : a preconditioner applied to each Greens functions [ntg]
"""
function IBD(ntg, nt, nr, nts;
	       fft_threads=true,
	       fftwflag=FFTW.PATIENT,
	       sxp=Sxparam(1,:positive),
	       sx_fix_zero_lag_flag=true,
	       dobs=nothing,
	       gobs=nothing,
	       sobs=nothing,
	       verbose=false,
	       fmin=0.15,
	       fmax=0.45,
	       size_check=true,
	       )

	# determine type of IBD
	if(!(dobs===nothing))
		T=eltype(dobs)
	else
		T1=eltype(sobs)
		T2=eltype(gobs)
		(T1≠T2) ? error("type difference") : (T=T1)
	end

	if(size_check)
		if(ntg+nts-1 ≠ nt)
			error("invalid sizes for convolutional model")
		end
	end

	# use maximum threads for fft only for patient
	fft_threads &&  (FFTW.set_num_threads(Sys.CPU_THREADS))

	# store observed data
	om=ObsModel(ntg, nt, nr, nts, T, d=dobs, g=gobs, s=sobs)

	# create interferometric Optim
	optm=OptimModel(2*ntg-1, 2*nt-1,
		 binomial(nr, 2)+nr, 2*nts-1, T, fftwflag=fftwflag,
	slags=[nts-1, nts-1],
	dlags=[nt-1, nt-1],
	glags=[ntg-1, ntg-1],
		 )

	# inversion variables allocation
	gx=X(length(optm.cal.g),T)

	if(sxp.attrib==:positive)
		if(sx_fix_zero_lag_flag)
			sx=X(nts-1,T) # not including zero lag
		else
			sx=X(nts,T)
		end
	elseif(sxp.attrib==:all)
		if(sx_fix_zero_lag_flag)
			sx=X(length(optm.cal.s)-1,T) # not including zero lag
		else
			sx=X(length(optm.cal.s),T)
		end
	else
		error("invalid sxp")
	end

	err=DataFrame(g=[], g_nodecon=[], s=[], d=[], whiteness_obs=[], whiteness_cal=[])

	# dummy
	opG=LinearMap{T}(x->0.0, y->0.0,1,1, ismutating=true)
	opS=LinearMap{T}(x->0.0, y->0.0,1,1, ismutating=true)

	pbandpassG=P_bandpass(T, fmin=fmin, fmax=fmax, nt=optm.ntg)
	pbandpassS=P_bandpass(T, fmin=fmin, fmax=fmax, nt=optm.nts)

	pa=IBD{T}(om, optm, gx, sx, sxp, sx_fix_zero_lag_flag,
	:g, verbose, err, opG, opS, pbandpassG, pbandpassS,
	zeros(nts,nts))

	# update operators
	pa.opG=create_operator(pa, G())
	pa.opS=create_operator(pa, S())

	#=
	if(sx_fix_zero_lag_flag)
		if(sxp==:positive)
			# adjust sprecon
			pa.sx.precon[1]=0.0 # do not update zero lag
			pa.sx.preconI[1]=0.0 # do not update zero lag
		elseif(sxp==:all)
			pa.sx.precon[nts]=0.0 # do not update zero lag
			pa.sx.preconI[nts]=0.0 # do not update zero lag
		else
			error("invalid sxp")
		end
	end
	=#

	gobs=hcat(Conv.xcorr(pa.om.g)...)
	sobs=hcat(Conv.xcorr(pa.om.s)...)

	# obs.g <-- gobs
	replace!(pa.optm, gobs, :obs, :g )
	# obs.s <-- sobs
	replace!(pa.optm, sobs, :obs, :s )
	# obs.d <-- dobs
	dobs=hcat(Conv.xcorr(pa.om.d)...) # do a cross-correlation

	copyto!(pa.optm.obs.d, dobs) # overwrites the forward modelling done in previous steps

	# normalize the observed data to 1.0
	scal=maximum(pa.optm.obs.d)

	rmul!(pa.optm.obs.d, inv(scal))
	rmul!(pa.optm.obs.g, inv(scal)) # to be consistent, mutiple gobs, if necessary

	initialize!(pa)

	return pa

end

function update_prepare!(pa::IBD, ::S)
	# fixing zero lag of source is equivalent to fitting d-g, for sxp.n==1
	# note that it is necessary while using IterativeSolvers.jl
	if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
		for ir in 1:pa.optm.nr
			for it in 1:pa.optm.ntg
				pa.optm.obs.d[pa.om.nt-pa.om.ntg+it,ir] -= pa.optm.cal.g[it,ir]
			end
		end
	end
end
function update_prepare!(pa::IBD, ::G)
	if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
		pa.optm.cal.s[pa.om.nts]=1.0
	end
end
function update_finalize!(pa::IBD, ::S)
	# add g back to d, after optimizing s
	if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
		for ir in 1:pa.optm.nr
			for it in 1:pa.optm.ntg
				pa.optm.obs.d[pa.om.nt-pa.om.ntg+it,ir] += pa.optm.cal.g[it,ir]
			end
		end
	end

	if(FILT_FLAG)
		# put sa[0]=1.0, before filtering
		if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
			pa.optm.cal.s[pa.om.nts]=1.0
		end
		bp=create_operator(pa.pbandpassS)
		copyto!(pa.optm.ds,pa.optm.cal.s)
		mul!(pa.optm.sfilt, bp, pa.optm.cal.s)
		copyto!(pa.optm.cal.s,pa.optm.sfilt)
		#
		if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
			pa.optm.cal.s[pa.om.nts]=0.0
		end
	end

	if(SDP_FLAG)
		if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
			pa.optm.cal.s[pa.om.nts]=1.0
		end
		# view the positive lags
		sp=view(pa.optm.cal.s, pa.om.nts:2*pa.om.nts-1)
		# create a Toeplitz matrix
		toeplitz!(pa.toepS,sp)
		eigS=eigen(pa.toepS)
		v=eigS.values; vv=eigS.vectors;
		if(any(v.<0.0))
			println("SDP PROJECTED\t", v)
		end
		v=Diagonal(broadcast(x->max(x,0),v))
		copyto!(pa.toepS,vv*v*vv')

		# view postive and zero lags of first column
		sp=view(pa.toepS,:,1)
		for i in 1:pa.om.nts-1
			# put same in positive lags (normalize with zero lag)
			pa.optm.cal.s[pa.om.nts+i]=sp[i+1]*inv(sp[1])
			# put same in negative lags (normalize with zero lag)
			pa.optm.cal.s[pa.om.nts-i]=sp[i+1]*inv(sp[1])
		end

		if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
			pa.optm.cal.s[pa.om.nts]=0.0
		end
	end

	if(STF_FLAG)
		for i in eachindex(pa.optm.cal.s)
			if(pa.optm.cal.s[i] <0.0)
				pa.optm.cal.s[i]=0.0
			end
		end
	end
end
function update_finalize!(pa::IBD, ::G)
	if(isequal(pa.sxp.n,1) && pa.sx_fix_zero_lag_flag)
		pa.optm.cal.s[pa.om.nts]=0.0
	end

	if(FILT_FLAG)
		bp=create_operator(pa.pbandpassG)
		copyto!(pa.optm.dg,pa.optm.cal.g)
		mul!(pa.optm.gfilt, bp, pa.optm.cal.g)
		copyto!(pa.optm.cal.g,pa.optm.gfilt)

		# put back for autocorr, if focusing
		if(minimum(pa.gx.preconI) == 0.0)
			irr=1  # auto correlation index
			for ir in 1:pa.om.nr
				gcopyv=view(pa.optm.dg,:,irr)
				gv=view(pa.optm.cal.g, :, irr)
				copyto!(gv, gcopyv)
				irr+=pa.om.nr-(ir-1)
			end
		end
	end
end

function model_to_x!(x, pa::IBD, ::S)
	if(pa.sxp.attrib==:all)
		if(pa.sx_fix_zero_lag_flag)
			nts=pa.om.nts
			for i in 1:nts-1
				x[i]=^(pa.optm.cal.s[i],inv(pa.sxp.n))*pa.sx.precon[i]
			end
			for i in nts+1:2*nts-1
				x[i-1]=^(pa.optm.cal.s[i],inv(pa.sxp.n))*pa.sx.precon[i-1]
			end
		else
			for i in eachindex(x)
				x[i]=^(pa.optm.cal.s[i],inv(pa.sxp.n))*pa.sx.precon[i]
			end
		end
	elseif(pa.sxp.attrib==:positive)
		if(pa.sx_fix_zero_lag_flag)
			for i in eachindex(x)
				x[i]=^(pa.optm.cal.s[pa.om.nts+i],inv(pa.sxp.n))*pa.sx.precon[i] # just take positive lags
			end
		else
			for i in eachindex(x)
				x[i]=^(pa.optm.cal.s[i+pa.om.nts-1],inv(pa.sxp.n))*pa.sx.precon[i] # just take positive lags
			end
		end
	else
		error("invalid sxp")
	end
	return nothing
end
function model_to_x!(x, pa::IBD, ::G)
	for i in eachindex(x)
		x[i]=pa.optm.cal.g[i]*pa.gx.precon[i] 		# multiply by gprecon
	end
	return nothing
end

function x_to_model!(x, pa::IBD, ::S)
	if(pa.sxp.attrib==:all)
		if(pa.sx_fix_zero_lag_flag)
			for i in 1:pa.om.nts-1
				pa.optm.cal.s[i]=^(x[i],pa.sxp.n)*pa.sx.preconI[i]
			end
			for i in pa.om.nts+1:2*pa.om.nts-1
				pa.optm.cal.s[i]=^(x[i-1],pa.sxp.n)*pa.sx.preconI[i-1]
			end
			if(isequal(pa.sxp.n,1))
				pa.optm.cal.s[pa.om.nts]=0.0
			else
				pa.optm.cal.s[pa.om.nts]=1.0
			end

		else
			for i in eachindex(pa.optm.cal.s)
				pa.optm.cal.s[i]=^(x[i],pa.sxp.n)*pa.sx.preconI[i]
			end
		end
	elseif(pa.sxp.attrib==:positive)
		if(pa.sx_fix_zero_lag_flag)
			for i in 1:pa.om.nts-1
				# put same in positive lags
				pa.optm.cal.s[pa.om.nts+i]=^(x[i],pa.sxp.n)*pa.sx.preconI[i]
				# put same in negative lags
				pa.optm.cal.s[pa.om.nts-i]=^(x[i],pa.sxp.n)*pa.sx.preconI[i]
			end
			if(isequal(pa.sxp.n,1))
				pa.optm.cal.s[pa.om.nts]=0.0
			else
				pa.optm.cal.s[pa.om.nts]=1.0
			end
		else
			for i in 1:pa.om.nts-1
				# put same in positive lags
				pa.optm.cal.s[pa.om.nts+i]=^(x[i+1],pa.sxp.n)*pa.sx.preconI[i+1]
				# put same in negative lags
				pa.optm.cal.s[pa.om.nts-i]=^(x[i+1],pa.sxp.n)*pa.sx.preconI[i+1]
			end
			pa.optm.cal.s[pa.om.nts]=^(x[1],pa.sxp.n)*pa.sx.preconI[1]
		end
	else
		error("invalid sxp")
	end
	return pa
end

function x_to_model!(x, pa::IBD, ::G)
	for i in eachindex(pa.optm.cal.g)
		pa.optm.cal.g[i]=x[i]*pa.gx.preconI[i]
	end
	return pa
end

"""
Zero out elements of g_ii, depending on nlags
"""
function focus!(pa::IBD, nlags=1)
	nr=pa.om.nr
	ntg=pa.om.ntg
	irr=1  # auto correlation index
	for ir in 1:nr
		for i in nlags:ntg-1
			# zero out Greens functions at non zero lags
			pa.optm.cal.g[ntg+i,irr]=0.0
			pa.optm.cal.g[ntg-i,irr]=0.0
		end
		irr+=nr-(ir-1)
	end

end

"""
Add a preconditioner to update only a few lags of the g_ii
Total lags that the precon is non-zero is given by 2*nlags-1
"""
function add_focusing_gprecon!(pa::IBD, nlags=1; k=1.0)

	nr=pa.om.nr
	ntg=pa.om.ntg
	gprecon=ones(pa.optm.ntg, pa.optm.nr).*k;

	irr=1  # auto correlation index
	for ir in 1:nr
		for i in nlags:ntg-1
			gprecon[ntg+i,irr]=0.0    # put zero at +ve lags
			gprecon[ntg-i,irr]=0.0    # put zero at -ve lags
		end
		gprecon[ntg,irr]=1.0    # put zero at -ve lags
		irr+=nr-(ir-1)
	end
	add_gprecon!(pa, gprecon)
	return pa
end

"""
Add a preconditioner to update only g_ii without focusing
"""
function add_autocorr_gprecon!(pa::IBD)

	nr=pa.om.nr
	ntg=pa.om.ntg
	gprecon=zeros(pa.optm.ntg, pa.optm.nr);

	irr=1  # auto correlation index
	for ir in 1:nr
		gprecon[:,irr]=1.0    # put one
		irr+=nr-(ir-1)
	end
	add_gprecon!(pa, gprecon)
	return pa
end

"""
Add focusing in pa.
Zero out pa.optm.cal.g for autocorrelations and non zero lags
"""
function add_focusing_gweights!(pa::IBD)

	nr=pa.om.nr
	ntg=pa.om.ntg
	gweights=zeros(pa.optm.ntg, pa.optm.nr);

	irr=1  # auto correlation index
	for ir in 1:nr
		for i in 1:ntg-1
			gweights[ntg+i,irr]=abs((i)/(ntg-1))
			gweights[ntg-i,irr]=abs((i)/(ntg-1))
		end
		irr+=nr-(ir-1)
	end
	add_gweights!(pa, gweights)
	return pa
end

"""
compute the amount of focusing in cross-correlated g
"""
function whiteness_focusing(pa::IBD, attrib=:obs)
	g=getfield(pa.optm, attrib).g
	nr=pa.om.nr
	ntg=pa.om.ntg

	irr=1  # auto correlation index
	J = 0.0
	for ir in 1:nr
		gg=view(g, :, irr)
		JJ=0.0
		#fact=inv(gg[ntg]*gg[ntg])
		for i in eachindex(gg)
			JJ += abs2(gg[i]) * abs2((ntg-i)/(ntg-1))
		end
		J += JJ
		irr+=nr-(ir-1)
	end
	J = J * inv(nr) # scale with number of receivers
	return J
end

function ibd!(pa::IBD, io=stdout)

	remove_gprecon!(pa, including_zeros=true)  # remove precon
	remove_gweights!(pa, including_zeros=true)
	#update_func_grad!(pa,goptim=[:ls,], gαvec=[1.]);
	update_all!(pa, io, max_reroundtrips=1, max_roundtrips=100000, roundtrip_tol=1e-8)

	err!(pa)
end

"""
Focused Interferometric Blind Deconvolution
"""
function fit!(pa::IBD, io=stdout; verbose=true, α=[Inf, 0.0], tol=[1e-10,1e-4])

	if(io===nothing)
		logfilename=joinpath(pwd(),string("XFIBD",Dates.now(),".log"))
		io=open(logfilename, "a+")
	end

	updates = (x,y) -> update_all!(pa, io,
				    max_reroundtrips=1, max_roundtrips=y,
				    roundtrip_tol=x, verbose=true)

	add_focusing_gweights!(pa) # adds weights
	for (iα, αv) in enumerate(α)
		write(io, string("## (",iα,"/",length(α),") homotopy for alpha=", αv,"\n"))
		# if set α=∞
		if(αv==Inf)
			focus!(pa)
			add_focusing_gprecon!(pa) # adds precon
			#update_func_grad!(pa,goptim=[:ls,], gαvec=[1.]);
			updates(tol[iα], 10000000)
			#add_autocorr_gprecon!(pa) # adds precon
			#updates(tol[iα])
			remove_gprecon!(pa, including_zeros=true)  # remove precon
		elseif(αv==0.0)
			remove_gprecon!(pa, including_zeros=true)  # remove precon
			remove_gweights!(pa, including_zeros=true)
			#update_func_grad!(pa,goptim=[:ls,], gαvec=[1.]);
			updates(tol[iα], 500)
		else
			error("not functional for arbitary alpha")
			remove_gprecon!(pa, including_zeros=true)  # remove precon
			add_focusing_gweights!(pa) # adds weights
			#update_func_grad!(pa,goptim=[:ls, :weights], gαvec=[1., αv]);
			updates(tol[iα])
		end
	end
	remove_gweights!(pa, including_zeros=true)
	err!(pa)
end

function initialize!(pa::IBD, all_zero=false)
	for i in eachindex(pa.optm.cal.s)
		pa.optm.cal.s[i]= ^(abs.(randn()),(pa.sxp.n-1)) # zero for sxp.n==1, otherwise randn
	end
#	pa.optm.cal.s[pa.om.nts]=1.0 # initialize zero lag to one
	if(all_zero)
		pa.optm.cal.g[:]=0.0
	else
		for i in eachindex(pa.optm.cal.g)
			x=(pa.gx.precon[i]≠0.0) ? randn() : 0.0
			pa.optm.cal.g[i]=x
		end
	end
end

function F!(pa::IBD,	x::AbstractVector, attrib::S)
	compute=(x!=pa.sx.last_x)
	if(compute)
		x_to_model!(x, pa, attrib) # modify pa.optm.cal.s
		copyto!(pa.sx.last_x, x)
		Conv.mod!(pa.optm.cal, Conv.D()) # modify pa.optm.cal.d
		return pa
	end
end

function F!(pa::IBD,	x::AbstractVector, ::G)
	compute=(x!=pa.gx.last_x)
	if(compute)
		x_to_model!(x, pa, G()) # modify pa.optm.cal.g
		copyto!(pa.gx.last_x, x)
		Conv.mod!(pa.optm.cal, Conv.D()) # modify pa.optm.cal.d
		return pa
	end
end

"""
Apply Fadj to dcal and return storage
"""
function Fadj!(pa::IBD, storage, x, dcal, ::S)
	fill!(storage, 0.)
	Conv.mod!(pa.optm.cal, Conv.S(), d=dcal, s=pa.optm.ds)
	if(pa.sxp.attrib == :positive)
		# stacking over +ve and -ve lags
		if(pa.sx_fix_zero_lag_flag)
			for j in 2:pa.om.nts
				storage[j-1] += pa.optm.ds[pa.om.nts-j+1] # -ve lags
				storage[j-1] += pa.optm.ds[pa.om.nts+j-1] # +ve lags
			end
		else
			for j in 2:pa.om.nts
				storage[j] += pa.optm.ds[pa.om.nts-j+1] # -ve lags
				storage[j] += pa.optm.ds[pa.om.nts+j-1] # +ve lags
			end
			storage[1]=pa.optm.ds[pa.om.nts]
		end
	elseif(pa.sxp.attrib == :all)
		if(pa.sx_fix_zero_lag_flag)
			nts=pa.om.nts
			for i in 1:nts-1
				storage[i] = pa.optm.ds[i] #
			end
			for i in nts+1:2*nts-1
				storage[i-1] = pa.optm.ds[i] #
			end
		else
			for i in eachindex(storage)
				storage[i] = pa.optm.ds[i] #
			end
		end
	else
		error("invalid attrib" )
	end
	# apply precon
	for i in eachindex(storage)
		if(iszero(pa.sx.precon[i]))
			storage[i]=0.0
		else
			storage[i] = storage[i]*pa.sxp.n*^(x[i],pa.sxp.n-1)*pa.sx.preconI[i]
		end
	end
end

function Fadj!(pa::IBD, storage, x, dcal, ::G)
	Conv.mod!(pa.optm.cal, Conv.G(), g=pa.optm.dg, d=dcal)
	for i in eachindex(storage)
		if(iszero(pa.gx.precon[i]))
			storage[i]=0.0
		else
			storage[i]=pa.optm.dg[i]/pa.gx.precon[i]
		end
	end
	return storage
end

#=

"""
Force real and non negative spectrum for g
"""
function weak_autocorr_constraint_g!(g, pac)
	Conv.pad!(g, pac.gpad, pac.glags[1], pac.glags[2], pac.np2)
	A_mul_B!(pac.gfreq, pac.gfftp, pac.gpad)
	ntot=1
	nr=size(pac.gfreq,2)
	for ir in 1:nr
		for i in size(pac.gfreq,1)
			pac.gfreq[i,ntot]=complex(real(pac.gfreq[i,ntot]),0.0)
			if(real(pac.gfreq[i,ntot]) < 0.0)
				pac.gfreq[i,ntot]=complex(0.0,0.0)
			end
		end
		ntot += (nr-ir+1)
	end
	A_mul_B!(pac.gpad, pac.gifftp, pac.gfreq)
	Conv.truncate!(g, pac.gpad, pac.glags[1], pac.glags[2], pac.np2)
end

"""
Force real and non negative spectrum for s
"""
function weak_autocorr_constraint_s!(s, pac)
	Conv.pad!(s, pac.spad, pac.slags[1], pac.slags[2], pac.np2)
	A_mul_B!(pac.sfreq, pac.sfftp, pac.spad)
	for i in eachindex(pac.sfreq)
		pac.sfreq[i]=complex(real(pac.sfreq[i]),0.0)
		if(real(pac.sfreq[i]) < 0.0)
			pac.sfreq[i]=complex(0.0,0.0)
		end
	end
	A_mul_B!(pac.spad, pac.sifftp, pac.sfreq)
	Conv.truncate!(s, pac.spad, pac.slags[1], pac.slags[2], pac.np2)
end
=#

"""
The cross-correlated Green's functions and the auto-correlated source signature have to be reconstructed exactly, except for a scaling factor.
"""
function err!(pa::IBD, io=stdout; cal=pa.optm.cal)
	(fg, b)	= Misfits.error_after_scaling(cal.g, pa.optm.obs.g) # error upto a scalar b
	(fs, b)	= Misfits.error_after_scaling(cal.s, pa.optm.obs.s) # error upto a scalar b
	f = Misfits.error_squared_euclidean!(nothing, cal.d, pa.optm.obs.d, nothing, norm_flag=true)

	xg_nodecon=hcat(Conv.xcorr(pa.om.d,Conv.P_xcorr(pa.om.nt, pa.om.nr, cglags=[pa.om.ntg-1, pa.om.ntg-1]))...)
	xgobs=hcat(Conv.xcorr(pa.om.g)...) # compute xcorr with reference g
	fg_nodecon = Misfits.error_squared_euclidean!(nothing, xg_nodecon, xgobs, nothing, norm_flag=true)

	push!(pa.err[!,:whiteness_obs], whiteness_focusing(pa, :obs))
	push!(pa.err[!,:whiteness_cal], whiteness_focusing(pa, :cal))

	push!(pa.err[!,:s],fs)
	push!(pa.err[!,:d],f)
	push!(pa.err[!,:g],fg)
	push!(pa.err[!,:g_nodecon],fg_nodecon)
	if(!(io===nothing))
		write(io,"Interferometric Blind Decon Errors\t\n")
		write(io,"==================\n")
		if(io==stdout)
			display(pa.err)
		else
			write(io, string(pa.err))
		end
		write(io, "\n")
	end
end


text:
abstract type Boundary end

abstract type SimpleBoundary <: Boundary end
# Assumption of SimpleBoundary is that it has fields
# ni::Complex{Float64} - Internal index of refraction
# ne::Complex{Float64} - External index of refraction

nInt(s::SimpleBoundary) = s.ni
nExt(s::SimpleBoundary) = s.ne

struct Point2D
  x::Float64
  y::Float64
end

import Base:-
import Base:+
import Base.*
import LinearAlgebra.dot
import LinearAlgebra.norm
convert(::Type{Array{Float64,1}}, p::Point2D) = [p.x, p.y]
function convert(::Type{Point2D}, a::Array{Float64,1})
  if length(a) != 2
    error("Only 2 element arrays can be converted to points")
  end
  Point2D(a[1],a[2])
end
function -(a::Point2D, b::Point2D)
  Point2D(convert(Array{Float64,1}, a)-convert(Array{Float64,1}, b))
end
function +(a::Point2D, b::Point2D)
  Point2D(convert(Array{Float64,1}, a)+convert(Array{Float64,1}, b))
end
function +(a::Point2D, b::Array{Float64,1})
  Point2D(convert(Array{Float64,1}, a)+b)
end
function *(a::Number, p::BoundaryIntegralMethod.Point2D)
  Point2D(a*p.x, a*p.y)
end
Point2D(a::Array{Float64}) = convert(Point2D, a)
function dot(a::Point2D, b::Point2D)
  dot(convert(Array{Float64,1}, a), convert(Array{Float64,1}, b))
end
function dot(a::Array{Float64,1}, b::Point2D)
  dot(a, convert(Array{Float64,1}, b))
end
function norm(a::Point2D)
  sqrt(dot(a, a))
end

"""
Return boundary point as array.
"""
function boundaryPointVector(b::Boundary, t::Float64)
  convert(Array{Float64, 1}, boundaryPoint(b, t))
end

"""
Return boundary normal vector as array
"""
function normVector(b::Boundary, t::Float64)
  xp = boundaryPointDeriv(b, t)
  [xp[2], -xp[1]]
end

"""
Return unit scaled boundary normal vector as array
"""
function unitNormVector(b::Boundary, t::Float64)
  n = normVector(b, t)
  n/norm(n)
end

# This doesn't work with some boundaries like criss cross, probably due to discontinuities
# in boundary derivative
# function arcLen(b::Boundary, t::Float64, steps::Int)
#   function ds(t′::Float64)
#     (dxdt, dydt) = boundaryPointDeriv(b, t′)
#     sqrt(dxdt^2 + dydt^2)
#   end
#   quadgk(ds, t-pi/steps, t+pi/steps)[1]
# end

function arcLen(b::Boundary, t::Float64, steps::Int)
  p1 = boundaryPoint(b, t)
  nextT = t+2π/steps
  if nextT > 2π
    nextT = nextT - 2π
  end
  p2 = boundaryPoint(b, nextT)
  norm(p2-p1)
end

function curvature(b::Boundary, t::Float64)
  xp = boundaryPointDeriv(b, t)
  xpp = boundaryPointSecondDeriv(b, t)
  (xp[1]xpp[2]-xp[2]xpp[1])/(xp[1]^2+xp[2]^2)^(3/2)
end

function nIdxRefr(b::Boundary, p::Point2D)
  isInside(b, p) ? nInt(b) : nExt(b)
end

struct Circle <: SimpleBoundary
  # Center of circle
  center::Point2D
  # Radius
  r::Float64
  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}
end

function Circle(r, ni, ne)
  Circle(Point2D(0,0), r, ni, ne)
end

function boundaryPoint(c::Circle, t::Float64)
  Point2D(c.center.x+cos(t)*c.r, c.center.y+sin(t)*c.r)
end

function boundaryPointDeriv(c::Circle, t::Float64)
  c.r*[-sin(t), cos(t)]
end

function boundaryPointSecondDeriv(c::Circle, t::Float64)
  -c.r*[cos(t), sin(t)]
end

# Inverse of radius of curvature: ρ = 1/κ
function curvature(c::Circle, t::Float64)
  1/c.r
end

function arcLen(c::Circle, t::Float64, steps::Int)
  c.r*2π/steps
end

function isInside(c::Circle, p::Point2D)
  p = p-c.center
  norm(p) < c.r
end

function boundaryField(c::Circle, N::Int, m::Int, k::Complex{Float64})
  xbar = [exp(im*m*j*π/N)*besselj(m, c.ni*k*c.r) for j in 0:2N-1]
  append!(xbar, [exp(im*m*j*π/N)*1/2*c.ni*k*(besselj(m-1,c.ni*k*c.r)-besselj(m+1,c.ni*k*c.r)) for j in 0:2N-1])
  xbar
end

struct DiscretizedBoundary <: Boundary
  underlying::Boundary
  N::Int64
  boundaryPointVectors::Array{Array{Float64}}
  boundaryPointDerivs::Array{Array{Float64}}
  normVectors::Array{Array{Float64}}
  absDiffMatrix::Array{Float64,2}
end

function DiscretizedBoundary(b::Boundary, N::Int64)
  boundaryPointVectors = [boundaryPointVector(b, j*π/N) for j in 0:(2N-1)]
  boundaryPointDerivs = [boundaryPointDeriv(b, j*π/N) for j in 0:(2N-1)]
  normVectors = [normVector(b, j*π/N) for j in 0:(2N-1)]
  absDiffMatrix = [norm(boundaryPoint(b, j*π/N)-boundaryPoint(b, l*π/N)) for j in 0:(2N-1), l in 0:(2N-1)]
  DiscretizedBoundary(b, N, boundaryPointVectors, boundaryPointDerivs, normVectors, absDiffMatrix)
end

function tToIdx(d::DiscretizedBoundary, t::Float64)
  Int64(round(d.N*t/π))+1
end

function boundaryPoint(d::DiscretizedBoundary, t::Float64)
  i = tToIdx(d, t)
  x = d.boundaryPointVectors[i]
  Point2D(x[1], x[2])
end

function boundaryPointVector(d::DiscretizedBoundary, t::Float64)
  i = tToIdx(d, t)
  d.boundaryPointVectors[i]
end

function boundaryPointDeriv(d::DiscretizedBoundary, t::Float64)
  i = tToIdx(d, t)
  d.boundaryPointDerivs[i]
end

function normVector(d::DiscretizedBoundary, t::Float64)
  i = tToIdx(d, t)
  d.normVectors[i]
end

function nInt(d::DiscretizedBoundary)
  nInt(d.underlying)
end

function nExt(d::DiscretizedBoundary)
  nExt(d.underlying)
end

function absDiff(d::DiscretizedBoundary, t::Float64, τ::Float64)
  i = tToIdx(d, t)
  j = tToIdx(d, τ)
  d.absDiffMatrix[i,j]
end

function isInside(b::DiscretizedBoundary, p::Point2D)
  try
    polygonWindingNumber(b.boundaryPointVectors, p) != 0
  catch e
    # Consider points on edge or vertex to be outside
    false
  end
end

function polygonWindingNumber(boundaryVecs::Array{Array{Float64}}, p::Point2D)
  # See: http://www.sciencedirect.com/science/article/pii/S0925772101000128
  # "The point in polygon problem for arbitrary polygons"
  # An implementation of Hormann-Agathos (2001) Point in Polygon algorithm
  r = BoundaryIntegralMethod.convert(Array{Float64, 1}, p)
  ω = 0
  @inline detq(q1,q2) = (q1[1]-r[1])*(q2[2]-r[2])-(q2[1]-r[1])*(q1[2]-r[2])
  @inline crossing(q1, q2) = (q1[2] < r[2]) != (q2[2] < r[2])
  @inline right_crossing(q1, q2) = (detq(q1,q2) > 0) == (q2[2] > q1[2])
  @inline modify_ω(q1, q2) = ω = ω + 2*(q2[2] > q1[2]) - 1

  if r == boundaryVecs[1]
    throw(VertexException())
  end
  for idx in 1:length(boundaryVecs)
    q1 = boundaryVecs[idx]
    q2 = idx == length(boundaryVecs) ? boundaryVecs[1] : boundaryVecs[idx+1]
    if q2[2] == r[2]
      if q2[1] == r[1]
          throw(VertexException())
      elseif (q1[2] == r[2]) && ((q2[1] > r[1]) == (q1[1] < r[1]))
          throw(EdgeException())
      end
    end
    if crossing(q1, q2)
      if q1[1] >= r[1]
        if q2[1] > r[1]
          modify_ω(q1, q2)
        elseif right_crossing(q1, q2)
          modify_ω(q1, q2)
        end
      elseif q2[1] > r[1] && right_crossing(q1, q2)
        modify_ω(q1, q2)
      end
    end
  end
  return ω
end

struct HeiderConcaveBoundary <: SimpleBoundary
  # Center of structure
  center::Point2D
  # Parameters for (x(t), y(t)) = (a cos t + b cos 2t - c, d sin t)
  # Heider paper uses a=1, b=.65, c=-.65, d=1.5
  a::Float64
  b::Float64
  c::Float64
  d::Float64
  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}
end

function HeiderConcaveBoundary(ni::Complex{Float64}, ne::Complex{Float64})
  HeiderConcaveBoundary(Point2D(0,0), 1, .65, -.65, 1.5, ni, ne)
end

function shapeFunc(b::HeiderConcaveBoundary, t::Float64)
  Point2D(b.a*cos(t) + b.b*cos(2t) - b.c, b.d*sin(t))
end

function boundaryPoint(b::HeiderConcaveBoundary, t::Float64)
  b.center+shapeFunc(b, t)
end

function boundaryPointDeriv(b::HeiderConcaveBoundary, t::Float64)
  [-b.a*sin(t)-2b.b*sin(2t), b.d*cos(t)]
end

function isInside(b::HeiderConcaveBoundary, p::Point2D)
  # figure out what t the angle of the poing would be at
  p = p - b.center
  ϕ = atan(p.y, p.x)
  # Range is -π < ϕ < π, place in range of t
  t = ϕ > 0 ? ϕ : 2π+ϕ
  p = p-c.center
  norm(p) < norm(boundaryPoint(b, t))
end

struct BowtieBoundary <: SimpleBoundary
  # Center of structure
  center::Point2D

  # Parameters for shape in polar coords r(ϕ) = r0(1+ϵ cos(2ϕ))
  r0::Float64
  ϵ::Float64

  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}
end

function BowtieBoundary(r0::Float64, ϵ::Float64, ni::Complex{Float64}, ne::Complex{Float64})
  BowtieBoundary(Point2D(0,0), r0, ϵ, ni, ne)
end

function radius(b::BowtieBoundary, t::Float64)
  b.r0*(1+b.ϵ*cos(2t))
end

function boundaryPoint(b::BowtieBoundary, t::Float64)
  r = radius(b, t)
  b.center+Point2D([r*cos(t), r*sin(t)])
end

function boundaryPointDeriv(b::BowtieBoundary, t::Float64)
  [-(1+2*b.ϵ+3*b.ϵ*cos(2t))*sin(t),
  (1-2*b.ϵ+3*b.ϵ*cos(2t))*cos(t)]
end

function isInside(b::BowtieBoundary, p::Point2D)
  # figure out what t the angle of the Point2D would be at
  p = p - b.center
  ϕ = atan(p.y,p.x)
  # Range is -π < ϕ < π, place in range of t
  t = ϕ > 0 ? ϕ : 2π+ϕ
  norm(p) < norm(boundaryPoint(b, t)-b.center)
end

# Another quadrupole boundary with normalized area, used in publications by J. Noeckel
struct NoeckelQuadrupoleBoundary <: SimpleBoundary
  # Center of structure
  center::Point2D

  # Parameters for shape in polar coords r(ϕ) = (1 + ϵ Cos[ϕ]) (1+ϵ^2 /2 )^(-1/2)
  ϵ::Float64

  # Normally this is 2 for the standard quadrupole boundary but it can be increased to
  # create more oscillations in the boundary
  cosFactor::Int

  # Scaling applied to x-axis lobes only
  xScale::Float64

  #Phase for angle argument, allows rotation of the shape
  phase::Float64

  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}

  # Scale up shape by this factor
  scale::Float64
end

function NoeckelQuadrupoleBoundary(ϵ::Float64, ni::Complex{Float64}, ne::Complex{Float64})
  NoeckelQuadrupoleBoundary(Point2D(0,0), ϵ, 2, 1.0, 1.0, 0., ni, ne, 1.0)
end

function NoeckelQuadrupoleBoundary(ϵ::Float64, ni::Complex{Float64}, ne::Complex{Float64}, scale::Float64)
  NoeckelQuadrupoleBoundary(Point2D(0,0), ϵ, 2, 1.0, 0., ni, ne, scale)
end

function NoeckelQuadrupoleBoundary(ϵ::Float64, cosFactor::Int, phase::Float64, ni::Complex{Float64}, ne::Complex{Float64})
  NoeckelQuadrupoleBoundary(Point2D(0,0), ϵ, cosFactor, 1.0, phase, ni, ne, 1.0)
end

function NoeckelQuadrupoleBoundary(center::Point2D, ϵ::Float64, cosFactor::Int, phase::Float64, ni::Complex{Float64}, ne::Complex{Float64})
  NoeckelQuadrupoleBoundary(center, ϵ, cosFactor, 1.0, phase, ni, ne, 1.0)
end

function radius(b::NoeckelQuadrupoleBoundary, t::Float64)
  (1 + b.ϵ*cos(b.cosFactor*t+b.phase))*(1 + b.ϵ^2/2)^(-1/2)
end

function boundaryPoint(b::NoeckelQuadrupoleBoundary, t::Float64)
  r = radius(b, t)*b.scale
  b.center+Point2D([b.xScale*r*cos(t), r*sin(t)])
end

function boundaryPointDeriv(b::NoeckelQuadrupoleBoundary, t::Float64)
  sq2 = sqrt(2)

  sint = sin(t)
  cost = cos(t)

  b.scale*[
    -b.xScale*sq2*(sint + b.ϵ*cos(b.cosFactor*t + b.phase)*sint + b.cosFactor*b.ϵ*cost*sin(b.cosFactor*t + b.phase)),
    sq2*(cost + b.ϵ*cos(b.cosFactor*t + b.phase)*cost - b.cosFactor*b.ϵ*sint*sin(b.cosFactor*t + b.phase))
  ]/sqrt(2+b.ϵ^2)
end

# Another shape based on octopole, like a 4-leaf clover
struct CloverBoundary <: SimpleBoundary
  # Center of structure
  center::Point2D

  # Deformation parameter
  ϵ::Float64

  # scaling applied to x lobes only
  xScale::Float64

  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}

  # Scale up shape by this factor
  scale::Float64
end

function CloverBoundary(ϵ::Float64, ni::Complex{Float64}, ne::Complex{Float64})
  CloverBoundary(Point2D(0,0), ϵ, 9.0/8.0, ni, ne, 1.0)
end

function radius(b::CloverBoundary, t::Float64)
  1 + 2*b.ϵ*atan(3*(1+cos(4t)))
end

function boundaryPoint(b::CloverBoundary, t::Float64)
  r = radius(b, t)*b.scale
  b.center+Point2D([b.xScale*r*cos(t), r*sin(t)])
end

function boundaryPointDeriv(b::CloverBoundary, t::Float64)

  sint = sin(t)
  cost = cos(t)
  term1 = (1+2*b.ϵ*atan(6*cos(2t)^2))
  term2 = -24*b.ϵ*sin(4t)/(1+36*cos(2t)^4)

  b.scale*[
    -b.xScale*term1*sint-b.xScale*term2*cost,
    term1*cost-term2*sint
  ]
end

# Another shape based on octopole, like a 4-leaf clover
struct SkewedOctupole <: SimpleBoundary
  # Center of structure
  center::Point2D

  # Deformation amount
  ϵ::Float64

  # Deformation applied to one corner
  η::Float64

  # t0 sets placement of deformation
  t0::Float64

  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}

  # Scale up shape by this factor
  scale::Float64
end

function SkewedOctupole(ϵ::Float64, η::Float64, t0::Float64, ni::Complex{Float64}, ne::Complex{Float64})
  SkewedOctupole(Point2D(0,0), ϵ, η, t0, ni, ne, 1.0)
end

function radius(b::SkewedOctupole, t::Float64)
  1 + b.ϵ*(1-b.η*((1+cos(t-b.t0))/2)^4)*cos(4t)
end

function boundaryPoint(b::SkewedOctupole, t::Float64)
  r = radius(b, t)*b.scale
  b.center+Point2D([r*cos(t), r*sin(t)])
end

function boundaryPointDeriv(b::SkewedOctupole, t::Float64)

  sintt0 = sin(t-b.t0)
  costt0 = cos(t-b.t0)
  cos4t = cos(4t)
  sin4t = sin(4t)
  term1 = 1/4*b.ϵ*b.η*cos4t*sintt0*(costt0+1)^3 -
    4b.ϵ*sin4t*(1-1/16*b.η*(costt0+1)^4)
  term2 = b.ϵ*cos4t*(1-1/16*b.η*(costt0+1)^4)+1

  b.scale*[
    cos(t)*term1-sin(t)*term2
    cos(t)*term2+sin(t)*term1
  ]
end

# Needed for Wiersig
# function boundaryPointSecondDeriv(b::NoeckelQuadrupoleBoundary, t::Float64)
#   sq2 = sqrt(2)
#
#   sint = sin(t)
#   cost = cos(t)
#
#   -b.scale*sq2/sqrt(2+b.ϵ^2)*[
#     cost + (1+b.cosFactor^2)*b.ϵ*cos(b.cosFactor*t + b.phase)*cost - 2*b.cosFactor*b.ϵ*sint*sin(b.cosFactor*t + b.phase),
#     sint + (1+b.cosFactor^2)*b.ϵ*cos(b.cosFactor*t + b.phase)*sint + 2*b.cosFactor*b.ϵ*cost*sin(b.cosFactor*t + b.phase)
#   ]
#
# end

# Another quadrupole boundary with normalized area, used in publications by J. Noeckel
struct ReuleauxBoundary <: SimpleBoundary
  # Center of structure
  center::Point2D
  # Scale radius
  scale::Float64
  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}
  # (n, a_n) tuple coefficients defining the parameterized boundary shape
  coefficientTuples::Array{Tuple{Int, Complex},1}
end

function B1ReuleauxBoundary(ni::Complex{Float64}, ne::Complex{Float64}; scale::Float64=1.0/12)
  # B1 coefficients from 1539-3755/2014/90(2)/022903(15)
  b1Coeffs = [12., 0., 0., 3.0im/2, 0., 3.]
  rest = b1Coeffs[2:end]
  coeffTuples = vcat((0, b1Coeffs[1]), [(-n, conj(an)) for (n, an) in enumerate(rest)], [(n,an) for (n, an) in enumerate(rest)])
  ReuleauxBoundary(Point2D(0.,0.), scale, ni, ne, coeffTuples)
end

# The complex parameter describing the shape, where (x, y) = (re(z), im(z))
function z(r::ReuleauxBoundary, t::Float64)
  # Check for an != 0 is so we don't divide by 0 for a_1
  -im*r.coefficientTuples[1][2]-im*sum([an/(n+1)*(exp(im*(n+1)*t)-1) for (n, an) in r.coefficientTuples if an != 0])
end

# dz/dt
function z′(r::ReuleauxBoundary, t::Float64)
  # Check for an != 0 is so we don't divide by 0 for a_1
  sum([an*exp(im*(n+1)*t) for (n, an) in r.coefficientTuples if an != 0])
end

function boundaryPoint(r::ReuleauxBoundary, t::Float64)
  zPt = z(r, t)
  r.center+r.scale*Point2D([real(zPt), imag(zPt)])
end

function boundaryPointDeriv(r::ReuleauxBoundary, t::Float64)
  z′Pt = z′(r,t)
  r.scale*[real(z′Pt), imag(z′Pt)]
end

struct CrissCrossBoundary <: SimpleBoundary
  # Center of structure
  center::Point2D

  # Deformation parameter
  ϵ::Float64

  # Internal index of refraction
  ni::Complex{Float64}
  # External index of refraction
  ne::Complex{Float64}

  # Scale up shape by this factor
  scale::Float64

  # Length of right arm
  rLen::Float64

  # Length of left arm
  lLen::Float64

  tFlat::Float64

  # Shift parameter makes arm thinner but leaves the lobeSpeed
  shift::Float64
  CrissCrossBoundary(center, ϵ, ni, ne, scale, rLen, lLen) = new(center, ϵ, ni, ne, scale, rLen, lLen, tFlat(ϵ), 0.0)
  CrissCrossBoundary(center, ϵ, ni, ne, scale, rLen, lLen, shift) = new(center, ϵ, ni, ne, scale, rLen, lLen, tFlat(ϵ), shift)
  CrissCrossBoundary(ϵ, ni, ne) = new(Point2D(0,0), ϵ, ni, ne, 1.0, 2.0, 1.7, tFlat(ϵ), 0.0)
end

function rVec(b::CrissCrossBoundary, t::Float64)
  [cos(t), sin(t)]*(1+b.ϵ*cos(4t))
end

# Parameter at which shape becomes horizontal
function tFlat(ϵ::Float64)
  numeratorFactor = (sqrt(2)*sqrt(ϵ*(-5+13*ϵ)))/ϵ
  atan(sqrt(6+numeratorFactor)/(2*sqrt(5)),
    sqrt(14-numeratorFactor)/(2*sqrt(5)))
end

# Define squircle shape
function squircle(a::Float64, b::Float64, δ::Float64, sign1::Int,
    sign2::Int, t::Float64)
  [sign1*a*abs((sign1*cos(t)))^(1/2)+δ, sign2*b*abs(sin(sign2*t))^(1/2)]
end

function rEllipse(a::Float64, b::Float64, t::Float64)
  a*b/sqrt(b^2*cos(t)^2+a^2*sin(t)^2)
end

function circleMap(a::Float64, b::Float64, t::Float64)
  π/2*(a-rEllipse(a, b, t))/(a-b)
end

function radius(b::CrissCrossBoundary, t::Float64)
  (x0, y0) = rVec(b, b.tFlat)
  if b.tFlat <= t <= (π-b.tFlat)
    [0.0, -b.shift/2]+rVec(b, t)
  elseif (π+b.tFlat) <= t <= (2π-b.tFlat)
    [0.0, b.shift/2]+rVec(b, t)
  elseif 0 <= t <= b.tFlat
    squircle(b.rLen, y0-b.shift/2, x0, 1, 1, circleMap(b.rLen, y0, t/b.tFlat*π/2))
  elseif (π-b.tFlat) <= t <= π
    squircle(b.lLen, y0-b.shift/2, -x0, -1, 1, π-circleMap(b.lLen, y0, (π-t)/b.tFlat*π/2))
  elseif π <= t <= π+b.tFlat
    squircle(b.lLen, y0-b.shift/2, -x0, -1, -1, π+circleMap(b.lLen, y0, (t-π)/b.tFlat*π/2))
  else
    squircle(b.rLen, y0-b.shift/2, x0, 1, -1, 2π-circleMap(b.rLen, y0, (2π-t)/b.tFlat*π/2))
  end
end

function boundaryPoint(b::CrissCrossBoundary, t::Float64)
  r = radius(b, t)*b.scale
  b.center+Point2D([r[1], r[2]])
end

function lobeTangent(b::CrissCrossBoundary, t::Float64)
  denom = sqrt(4+34b.ϵ^2+8b.ϵ*cos(4t)-30b.ϵ^2*cos(8t))
  [(-2*sin(t)-3*b.ϵ*sin(3t)-5*b.ϵ*sin(5t))/denom,
    (2*cos(t)-3*b.ϵ*cos(3t)+5*b.ϵ*cos(5t))/denom]
end

function squircleTangent(a::Float64, b::Float64, δ::Float64, sign1::Int,
    sign2::Int, u::Float64)
  # Avoid issue with negative 0 sometimes being passed in which is outside domain
  u = abs(u)
  try
    if sign1 == 1 && sign2 == 1
      [-a/sqrt(abs(a^2+b^2*cot(u)^3)), b/sqrt(abs(b^2+a^2*tan(u)^3))]
    elseif sign1 == -1 && sign2 == 1
      [-a/sqrt(abs(a^2-b^2*cot(u)^3)), -b/sqrt(abs(b^2-a^2*tan(u)^3))]
    elseif sign1 == -1 && sign2 == -1
      [a/sqrt(abs(a^2+b^2*cot(u)^3)), -b/sqrt(abs(b^2+a^2*tan(u)^3))]
    else
      [a/sqrt(abs(a^2-b^2*cot(u)^3)), b/sqrt(abs(b^2-a^2*tan(u)^3))]
    end
  catch e
    @error "Domain error? $a, $b, $δ, $sign1, $sign2, $u\n$e" exception=e
  end
end

# Normalized tangent, needs speed factor to work for Heider
function normTangent(b::CrissCrossBoundary, t::Float64)
  (x0, y0) = rVec(b, b.tFlat)
  if b.tFlat <= t < (π-b.tFlat) || (π+b.tFlat) <= t < (2π-b.tFlat)
    lobeTangent(b, t)
  elseif 0 <= t < b.tFlat
    squircleTangent(b.rLen, y0-b.shift/2, x0, 1, 1, circleMap(b.rLen, y0, t/b.tFlat*π/2))
  elseif (π-b.tFlat) <= t < π
    squircleTangent(b.lLen, y0-b.shift/2, -x0, -1, 1, π-circleMap(b.lLen, y0, (π-t)/b.tFlat*π/2))
  elseif π <= t < π+b.tFlat
    squircleTangent(b.lLen, y0-b.shift/2, -x0, -1, -1, π+circleMap(b.lLen, y0, (t-π)/b.tFlat*π/2))
  else
    squircleTangent(b.rLen, y0-b.shift/2, x0, 1, -1, 2π-circleMap(b.rLen, y0, (2π-t)/b.tFlat*π/2))
  end
end

function squircleSpeed(a::Float64, b::Float64, t::Float64, tFlat::Float64)
  prefactor = 1/tFlat*π/2
  trigFactor = π*(a - a*b/sqrt(b^2*cos(t)^2+a^2*sin(t)^2))/(2*(a-b))
  spd = if 0 < t < π/2 && !isapprox(0, trigFactor) && !isapprox(π/2, trigFactor)
    spdSquared = (a+b)^2*π^2*cos(t)^2*sin(t)^2*
      (a^2*b^4*cos(trigFactor)*cot(trigFactor)+a^4*b^2*sin(trigFactor))*tan(trigFactor) /
      (16*(b^2*cos(t)^2 + a^2*sin(t)^2)^3)
      sqrt(spdSquared)
  elseif t <= 0 || !isapprox(0, trigFactor)
    sqrt(1/4*a*(a+b)*π)
  else
    sqrt(1/4*b*(a+b)*π)
  end
  prefactor*spd
end

function lobeSpeed(b::CrissCrossBoundary, t::Float64)
  sqrt((2+17*b.ϵ^2+4*b.ϵ*cos(4t)-15*b.ϵ^2*cos(8t))/2)
end

function boundaryPointDeriv(b::CrissCrossBoundary, t::Float64)
  (x0, y0) = rVec(b, b.tFlat)
  speed = if b.tFlat <= t <= (π-b.tFlat) || (π+b.tFlat) <= t <= (2π-b.tFlat)
    lobeSpeed(b, t)
  elseif 0 <= t <= b.tFlat
    squircleSpeed(b.rLen, y0-b.shift/2, t/b.tFlat*π/2, b.tFlat)
  elseif (π-b.tFlat) <= t <= π
    squircleSpeed(b.lLen, y0-b.shift/2, (π-t)/b.tFlat*π/2, b.tFlat)
  elseif π <= t <= π+b.tFlat
    squircleSpeed(b.lLen, y0-b.shift/2, (t-π)/b.tFlat*π/2, b.tFlat)
  else
    squircleSpeed(b.rLen, y0-b.shift/2, (2π-t)/b.tFlat*π/2, b.tFlat)
  end
  speed*normTangent(b, t)
end


text:
# Solvers
using Ipopt
using GLPK
using SCS
using HiGHS
using Cbc

ipopt_optimizer =
    JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6, "print_level" => 0)
fast_ipopt_optimizer = JuMP.optimizer_with_attributes(
    Ipopt.Optimizer,
    "print_level" => 0,
    "max_cpu_time" => 5.0,
)
# use default print_level = 5 # set to 0 to disable
GLPK_optimizer =
    JuMP.optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)
scs_solver = JuMP.optimizer_with_attributes(
    SCS.Optimizer,
    "max_iters" => 100000,
    "eps" => 1e-4,
    "verbose" => 0,
)

HiGHS_optimizer = JuMP.optimizer_with_attributes(
    HiGHS.Optimizer,
    "time_limit" => 100.0,
    "log_to_console" => false,
)

cbc_optimizer = JuMP.optimizer_with_attributes(Cbc.Optimizer)


text:
using ProgressBars,Plots, Statistics, JLD, LaTeXStrings
function Boundary(i,L)
    if i>L
        return i-L
    elseif i<1
        return L+i
    else
        return i
    end
end
function height_cal(arr)
    h=[]
    for col in eachcol(arr)
        height=0
        for element in col
            if element!=0
                height+=1
            end
        end
        push!(h, height)
    end
    return maximum(h)
end

function RandomWalk(Network_, L, FirstPos, t_limit, DirectionList)
    CurrentPos=FirstPos
    while true
        temp=CurrentPos .+ rand(DirectionList)
        temp= [temp[1],Boundary(temp[2],L)]
        if Network_[temp...]<=0
            CurrentPos=temp
            if CurrentPos[1]==1
                Network_[CurrentPos...]-=1
                return Network_, CurrentPos
            elseif CurrentPos[1]> t_limit
                return Network_, 0
            end
            for neighbor in DirectionList
                if neighbor[2]+CurrentPos[2]==0
                    if Network_[neighbor[1]+CurrentPos[1], L] > 0
                        Network_[CurrentPos...] =-1
                        return Network_, CurrentPos
                    end
                elseif neighbor[2]+CurrentPos[2]==L+1
                    if Network_[neighbor[1]+CurrentPos[1], 1] > 0
                        Network_[CurrentPos...] =-1
                        return Network_, CurrentPos
                    end
                elseif Network_[(CurrentPos .+ neighbor)...] > 0
                    Network_[CurrentPos...]=-1
                    return Network_, CurrentPos
                end
            end
        end
    end
end
direction_list=[[0,1],[0,-1],[1,0],[-1,0]]
L=20
Network_=zeros(Int,L,L)
N=800
color=1
seperation=5
for particle in ProgressBar(1:N)
    b_limit=height_cal(Network_) + 3
    t_limit= b_limit+seperation
    if t_limit<L
        Network_,FinalPos=RandomWalk(Network_, L, [b_limit, rand(1:L)], t_limit, direction_list)
        if FinalPos==0
            continue
        elseif Network_[FinalPos...]==-3
            Network_[FinalPos...]=color
        end
        if particle%(10*L*color)==0
            color+=1
        end
    end
end
Network_
Network_[findall(x->x<0, Network_)].=0
heatmap(Network_, dpi=400)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.6\\Figs\\4.6_strange.png")


text:
function jet_mixing(settings::PyObject, data::PyObject, ac::PyObject, n_t::Int64, M_0, c_0, theta, V_j_star, rho_j_star, A_j_star, Tt_j_star)

    # Initialize solution
    T = eltype(V_j_star)
    msap_jet_mixing = zeros(T, (n_t, settings.N_f))

    r_s_star = 0.3048/sqrt(settings.A_e)
    jet_delta = 0.

    # Calculate jet mixing
    log10Vja0 = log10.(V_j_star)

    # Calculate density exponent (omega)
    # Source: Zorumski report 1982 part 2. Chapter 8.4 Table II
    # Source: Hoch - Studies of the influence of density on jet noise: extend the
    array_1 = range(-0.45, 0.6, step=0.05)
    array_2 = [-1.0, -0.9, -0.76, -0.58, -0.41, -0.22, 0.0, 0.22, 0.5, 0.77, 1.07, 1.39, 1.74, 1.95, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
    if any(x->x<-0.45, log10Vja0) || any(x->x>0.6, log10Vja0)
        throw(DomainError(x," log10(V_jet/c_0) contains value outside the interpolation range of the jet density exponent, omega."))
    end
    f_omega = PCHIPInterpolation.Interpolator(array_1, array_2)
    omega = reshape(f_omega.(log10Vja0), (n_t, 1))

    # Calculate power deviation factor (P)
    # Source: Zorumski report 1982 part 2. Chapter 8.4 Table III
    array_1 = range(-0.45, 0.4, step=0.05)
    array_2 = [-0.13, -0.13, -0.13, -0.13, -0.13, -0.13, -0.12, -0.1, -0.05, 0.0, 0.1, 0.21, 0.32, 0.41, 0.43, 0.41,0.31, 0.14]
    if any(x->x<-0.4, log10Vja0) || any(x->x>0.4, log10Vja0)
        throw(DomainError(x," log10(V_jet/c_0) contains value outside the interpolation range of the power deviation factor, P."))
    end
    f_log10P = PCHIPInterpolation.Interpolator(array_1, array_2)
    log10P = f_log10P.(log10Vja0)
    P_function = 10 .^log10P

    # Calculate acoustic power (Pi_star)
    # Source: Zorumski report 1982 part 2. Chapter 8.4 Equation 3
    K = 6.67e-5
    Pi_star = reshape(K * rho_j_star .^omega .* V_j_star.^8 .* P_function, (n_t,1))

    # Calculate directivity function (D)
    # Source: Zorumski report 1982 part 2. Chapter 8.4 Table IV
    f_log10D = LinearInterpolation((data.jet_D_velocity, data.jet_D_angles), data.jet_D)
    log10D = f_log10D.(theta, log10Vja0)
    D_function = reshape(10 .^log10D, (n_t,1))

    # Calculate Strouhal frequency adjustment factor (xi)
    # Source: Zorumski report 1982 part 2. Chapter 8.4 Table V
    # Note: added extra lines in the table for V_j_star = 0, to not have set_index
    f_xi = LinearInterpolation((data.jet_xi_velocity, data.jet_xi_angles), data.jet_xi)
    xi = reshape(f_xi.(V_j_star, theta), (n_t, 1))

    # Calculate Strouhal number (St)
    D_j_star = sqrt.(4 * A_j_star / π)  # Jet diamater [-] (rel. to sqrt(settings.A_e))
    f_star = reshape(data.f, (1, settings.N_f)) .* sqrt(settings.A_e) ./ reshape(c_0, (n_t,1))
    St = (f_star .* reshape(D_j_star, (n_t,1))) ./ reshape(xi .* (V_j_star .- M_0), (n_t,1))
    log10St = log10.(St)

    # Source: Zorumski report 1982 part 2. Chapter 8.4 Table VI
    # Note: extended the data table temperature range ([1, 2, 2.5, 3, 3.5]) with linearly extrapolated values ([0, 1, 2, 2.5, 3, 3.5, 4, 5, 6, 7]) to avoid set_index for backward_diff
    f_log10F = LinearInterpolation((data.jet_F_angles, data.jet_F_temperature, data.jet_F_velocity, data.jet_F_strouhal), data.jet_F)
    #mlog10F_a_lg = f_log10F.(theta.*ones(1, settings.N_f), 3.5*ones(n_t, settings.N_f), log10Vja0*ones(1, settings.N_f), log10St)
    #mlog10F_b_lg = f_log10F.(theta.*ones(1, settings.N_f), 3.4*ones(n_t, settings.N_f), log10Vja0*ones(1, settings.N_f), log10St)
    #mlog10F_a_sm = f_log10F.(theta.*ones(1, settings.N_f), 1.1*ones(n_t, settings.N_f), log10Vja0*ones(1, settings.N_f), log10St)
    #mlog10F_b_sm = f_log10F.(theta.*ones(1, settings.N_f), 1.0*ones(n_t, settings.N_f), log10Vja0*ones(1, settings.N_f), log10St)

    # Linear extrapolation from the data table
    mlog10F = f_log10F.(theta.*ones(1, settings.N_f), Tt_j_star.*ones(1, settings.N_f), log10Vja0 .*ones(1, settings.N_f), log10St)
    #mlog10F[findall(Tt_j_star.*ones(1, settings.N_f) .> 3.5)] = ((mlog10F_a_lg .- mlog10F_b_lg) / (0.1) .* (Tt_j_star .- 3.5) .+ mlog10F_a_lg)[findall(Tt_j_star.*ones(1, settings.N_f) .> 3.5)]
    #mlog10F[findall(Tt_j_star.*ones(1, settings.N_f) .< 1)] = ((mlog10F_a_sm .- mlog10F_b_sm) / (0.1) .* (Tt_j_star .- 1.0) .+ mlog10F_b_sm)[findall(Tt_j_star.*ones(1, settings.N_f) .< 1)]
    F_function = 10 .^(-mlog10F / 10)

    # Calculate forward velocity index (m_theta)
    # Source: Zorumski report 1982 part 2. Chapter 8.4 Table VII
    array_1 = range(0, 180, step=10)
    array_2 = [3, 1.65, 1.1, 0.5, 0.2, 0, 0, 0.1, 0.4, 1, 1.9, 3, 4.7, 7, 8.5, 8.5, 8.5, 8.5, 8.5]
    f_m_theta = PCHIPInterpolation.Interpolator(array_1, array_2)
    m_theta = f_m_theta.(theta)

    # Calculate mean-square acoustic pressure (msap)
    # Source: Zorumski report 1982 part 2. Chapter 8.4 Equation 8
    msap_j = Pi_star .* reshape(A_j_star,(n_t,1)) / (4 * π * r_s_star^2) .* D_function .* F_function ./ reshape((1 .- M_0 .* cos.(π / 180. * (theta .- jet_delta))), (n_t, 1)) .* reshape(((V_j_star .- M_0) ./ V_j_star) .^ m_theta, (n_t,1))

    # Multiply with number of engines
    msap_j = msap_j * ac.n_eng

    # Normalize msap by reference pressure
    msap_jet_mixing = msap_j/settings.p_ref^2

    return msap_jet_mixing
end

# Function
function jet_shock(settings::PyObject, data::PyObject, ac::PyObject, n_t::Int64, M_0, c_0, theta, V_j_star, M_j, A_j_star, Tt_j_star)

    # Get elements of jet Mach number vector larger than 1
    T = eltype(V_j_star)
    msap_jet_shock = zeros(T, (n_t, settings.N_f))

    for i in range(1, n_t, step=1)

        # Only shock noise if jet Mach number is larger than 1
        if M_j[i] > 1

            r_s_star = 0.3048/sqrt(settings.A_e)
            jet_delta = 0.

            # Calculate beta function
            # Source: Zorumski report 1982 part 2. Chapter 8.5 Equation 4
            beta = (M_j[i].^2 .- 1).^0.5

            # Calculate eta (exponent of the pressure ratio parameter)
            # Source: Zorumski report 1982 part 2. Chapter 8.5 Equation 5
            if beta > 1
                if Tt_j_star[i] < 1.1
                    eta = 1.
                elseif Tt_j_star[i] >= 1.1
                    eta = 2.
                end
            else
                eta = 4.
            end

            # Calculate f_star
            # Source: Zorumski report 1982 part 2. Chapter 8.5 page 8-5-1 (symbols)
            f_star = reshape(data.f, (1, settings.N_f)) .* sqrt(settings.A_e) ./ c_0[i]

            # Calculate sigma parameter
            # Source: Zorumski report 1982 part 2. Chapter 8.5 Equation 3
            sigma = 7.80 * beta .* (1 .- M_0[i] .* cos.(π / 180 .* theta[i])) .* sqrt.(A_j_star[i]) .* f_star
            log10sigma = log10.(sigma)

            # Calculate correlation coefficient spectrum (C-function)
            # Source: Zorumski report 1982 part 2. Chapter 8.5 Table II
            # Note: extended array for log10sigma > 2 and log10sigma < -0.7
            array_1_c = range(-3.5, 3.5, step=0.1)
            array_2_c = [0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.703,0.71,0.714,0.719,0.724,0.729,0.735,0.74,0.74,0.74,0.735,0.714,0.681,0.635,0.579,0.52,0.46,0.4,0.345,0.29,0.235,0.195,0.15,0.1,0.06,0.03,0.015,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
            f_C = PCHIPInterpolation.Interpolator(array_1_c, array_2_c)
            C = f_C.(log10sigma)

            # Calculate W function
            # Source: Zorumski report 1982 part 2. Chapter 8.5 Equation 6-7
            b = 0.23077
            W = 0.
            for k in range(1, settings.N_shock-1, step=1)
                sum_inner = 0.
                for m in range(0,settings.N_shock - (k+1), step=1)
                    # Calculate q_km
                    q_km = 1.70 * k ./ V_j_star[i] .* (1 .- 0.06 * (m + (k + 1) / 2)) .* (1 .+ 0.7 .* V_j_star[i] .* cos.(π / 180 * theta[i]))

                    # Calculate inner sum (note: the factor b in the denominator below the sine should not be there: to get same graph as Figure 4)
                    sum_inner = sum_inner .+ sin.((b * sigma .* q_km / 2)) ./ (sigma .* q_km) .* cos.(sigma .* q_km)
                end

                W = W .+ sum_inner .* C.^(k.^2)
            end
            W = (4. ./ (settings.N_shock .* b)) .* W

            # Calculate the H function
            # Source: Zorumski report 1982 part 2. Chapter 8.5 Table III (+ linear extrapolation in logspace for log10sigma < 0; as given in SAEARP876)
            array_1_H = range(-2.5, 3.5, step=0.1)
            array_2_H = [-12.19,-11.81,-11.43,-11.05,-10.67,-10.29,-9.91,-9.53,-9.15,-8.77,-8.39,-8.01,-7.63,-7.25,-6.87,-6.49,-6.11,-5.73,-5.35,-4.97,-4.59,-4.21,-3.83,-3.45,-3.07,-2.69,-2.31,-1.94,-1.59,-1.33,-1.1,-0.94,-0.88,-0.91,-0.99,-1.09,-1.17,-1.3,-1.42,-1.55,-1.67,-1.81,-1.92,-2.06,-2.18,-2.3,-2.42,-2.54,-2.66,-2.78,-2.9,-3.02,-3.14,-3.26,-3.38,-3.5,-3.62,-3.74,-3.86,-3.98,-4.1]
            f_H = PCHIPInterpolation.Interpolator(array_1_H, array_2_H)
            log10H = f_H.(log10sigma)
            if Tt_j_star[i] < 1.1
                log10H = log10H - 0.2
            end

            # Calculate mean-square acoustic pressure (msap)
            # Source: Zorumski report 1982 part 2. Chapter 8.5 Equation 1
            msap_j = 1.92e-3 * A_j_star[i] / (4 * π * r_s_star^2) .* (1 .+ W) ./ (1 .- M_0[i] .* cos.(π / 180. * (theta[i] .- jet_delta))).^4 .* beta.^eta .* (10 .^log10H)

            # Multiply with number of engines
            msap_j = msap_j * ac.n_eng

            # Normalize msap by reference pressure
            msap_jet_shock[i,:] = msap_j/settings.p_ref^2
        end
    end

    return msap_jet_shock
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
"""
Stochastic Scheme for Uncertainty Quantification
Copyright (c) 2020-2021 Tianbai Xiao
"""

module Langevin

using Reexport
@reexport using PolyChaos
@reexport using KitBase
using KitBase.OffsetArrays
using KitBase.Plots
using KitBase.JLD2

export AbstractUQ,
       UQ1D,
       UQ2D,
       ran_chaos,
       chaos_ran,
       lambda_tchaos,
       t_lambdachaos
export uq_moments_conserve,
       uq_maxwellian,
       uq_prim_conserve,
       uq_conserve_prim,
       uq_conserve_prim!,
       uq_prim_conserve!,
       uq_sound_speed,
       uq_vhs_collision_time,
       uq_aap_hs_collision_time,
       uq_aap_hs_prim

include("struct.jl")
include("uq.jl")
include("filter.jl")
include("kinetic.jl")
include("initialize.jl")
include("flux_ctr.jl")
include("flux_sol.jl")
include("out.jl")
include("solver.jl")
include("step.jl")

end


text:
"""
Abstract type for interfaces for Nextion objects
"""
abstract type AbstractINextion end

NexID(interface::AbstractINextion) = interface._nid

include("IViewable.jl")
include("IBooleanValued.jl")
include("INumericalValued.jl")
include("IStringValued.jl")
include("IColourable.jl")
include("IFontStyleable.jl")
include("IPicturable.jl")
# include("ICallback.jl")
include("ITouchable.jl")
include("IWidthable.jl")
include("IHeightable.jl")

function abilities(obj::AbstractNexObject)
    ab = Symbol[]
    for field in fieldnames(typeof(obj))
        if typeof(getfield(obj, field)) <: AbstractINextion
            push!(ab, field)
        end
    end
    ab
end


text:
@testset "Label components" begin
    A = [true  true  false true;
         true  false true  true]
    lbltarget = [1 1 0 2;
                 1 0 2 2]
    lbltarget1 = [1 2 0 4;
                  1 0 3 4]
    @test label_components(A) == lbltarget
    @test label_components(A, [1]) == lbltarget1
    connectivity = [false true  false;
                    true  false true;
                    false true  false]
    @test label_components(A, connectivity) == lbltarget
    connectivity = trues(3,3)
    lbltarget2 = [1 1 0 1;
                  1 0 1 1]
    @test label_components(A, connectivity) == lbltarget2
    @test component_boxes(lbltarget) == Vector{Tuple}[[(1,2),(2,3)],[(1,1),(2,2)],[(1,3),(2,4)]]
    @test component_lengths(lbltarget) == [2,3,3]
    @test component_indices(lbltarget) == Array{Int64}[[4,5],[1,2,3],[6,7,8]]
    @test component_subscripts(lbltarget) == Array{Tuple}[[(2,2),(1,3)],[(1,1),(2,1),(1,2)],[(2,3),(1,4),(2,4)]]
    @test @inferred(component_centroids(lbltarget)) == Tuple[(1.5,2.5),(4/3,4/3),(5/3,11/3)]
end


text:
load_exch_met_map() = ldat(HumanGEM,
    "exch_met_map", ".jls";
    verbose = false
)