# julia

### Site Tools

arraysmatrix

 Arrays Introduction Matrices+

snippet.juliarepl
`julia> pkgchk.( [ "julia" => v"1.0.3" ] );`

# Matrix (2-Dimensional Arrays)

## Creating A Matrix from Values

A 2-D array can be created by specifying all elements in blank separated columns and semicolon (not comma!) separated rows:

snippet.juliarepl
```julia> x= [ 1 2 3 4 ; 5 6 7 8 ]
2×4 Array{Int64,2}:
1  2  3  4
5  6  7  8```
• A newline can be used instead of the semi-colon Watch the spaces: `x= [ 1 -2 ]` is a 1×2 array (matrix). In contrast, `x= [ 1-2 ]` is a 1-element vector.

A matrix is fundamentally not a vector of vectors. It is a rectangular object. Compare to a vector of vectors:

snippet.juliarepl
```julia> x= [ 1, 2, 3 ] ; y= [ 10, 11, 12 ];

julia> [ x, y ]
2-element Array{Array{Int64,1},1}:
[1, 2, 3]
[10, 11, 12]```

### Cat (hcat/vcat/hvcat) and Bracket Alternatives

Vectors are created with comma-separated lists. Matrices are created with space and semicolon separated lists.

The parenthesis notation above are “just” shortcuts.

snippet.juliarepl
```julia> [ 1 2 3 ] == hcat(1,2,3)
true

julia> [ 1;2;3 ] == vcat(1,2,3)
true

julia> hvcat( (3,3,3), 1,2,3,4,50,6,7,8,9 )
3×3 Array{Int64,2}:
1  2  3
4 50  6
7  8  9

julia> x= collect( 1:5 );  y= collect(11:15); hcat( x, y )	## place multiple vectors into matrix
5×2 Array{Int64,2}:
1  11
2  12
3  13
4  14
5  15```

## Printing A Matrix

snippet.juliarepl
```julia> x= [1 2 3 ; 4 5 6];

julia> println(x)
[1 2 3; 4 5 6]

julia> display(x)
2×3 Array{Int64,2}:
1  2  3
4  5  6```

## Initializing a Matrix with Same Values

`fill()` can initialize arbitrarily-dimensional arrays:

snippet.juliarepl
```julia> fill( NaN , (3,2) )	# or use fill( NaN, 3, 2 )
3×2 Array{Float64,2}:
NaN  NaN
NaN  NaN
NaN  NaN```

There are type-specialized versions for `zeros(Float64, 3, 3)` and `ones(Bool, 3, 3)`.

You could also create uninitialized matrices. Even if you immediately assign to your matrix afterwards, such uninitialized variables are rarely a good idea, and the cost of initialization with nada is trivial. Chances are that sooner or later, you may forget. Please always initialize your variables!

### Diagonal Matrices

The identity matrix has special support:

snippet.juliarepl
```julia> using LinearAlgebra

julia> Matrix( 2.0I, 3, 3 )		## An 'I' is following '2.0' creates a particular diagonal matrix
3×3 Array{Float64,2}:
2.0  0.0  0.0
0.0  2.0  0.0
0.0  0.0  2.0

julia> Matrix( I, 3, 3 )		## I is shortcut for 1.0I
3×3 Array{Float64,2}:
1.0  0.0  0.0
0.0  1.0  0.0
0.0  0.0  1.0

julia> Matrix( 10.0I, 3, 3 ) + I	## I is an identity matrix that reshapes itself to whatever is needed
3×3 Array{Float64,2}:
11.0   0.0   0.0
0.0  11.0   0.0
0.0   0.0  11.0```

For arbitrary diagonal matrices, you can use

snippet.juliarepl
```julia> using LinearAlgebra

julia> diagm( 0 => 1:4 )
4×4 Array{Int64,2}:
1  0  0  0
0  2  0  0
0  0  3  0
0  0  0  4```

## Pushing a Vector onto an Array

Remember: A matrix is not like a vector of vectors.

### Pushing Vectors onto Vectors

snippet.juliarepl
```julia> mtype= Vector{ Vector{Float64} }
Array{Array{Float64,1},1}

julia> M= mtype(undef,0)	## oddly, 'undef' is now required and nothing else is accepted
0-element Array{Array{Float64,1},1}

julia> M= mtype(undef, 2)	## an object of type M with two undef items
2-element Array{Array{Float64,1},1}:
#undef
#undef

julia> push!( M, [2.0, 3.0] )	## add a vector to the vector (stack)
3-element Array{Array{Float64,1},1}:
#undef
#undef
[2.0, 3.0]```

### Side-Append Column Vectors to Matrix

snippet.juliarepl
```julia> Matrix{Float64}(undef,2,0)			## Method 1: direct empty 2×0 array: 2 rows, 0 columns
2×0 Array{Float64,2}

julia> M= reshape( Vector{Float64}([]), 2, 0 )       ## Method 2: via reshape of Vector
2×0 Array{Float64,2}

julia> [M [ 1.0, 2.0 ]]                              ## space means hcat. hcat= cat(2,...) , i.e., dimension 2
2×1 Array{Float64,2}:
1.0
2.0

julia> [ [M [1.0, 2.0]]  [2.0, 3.0]]                            ## two 2 vectors
2×2 Array{Float64,2}:
1.0  2.0
2.0  3.0

julia> M= hcat(M, [ 3.0, 4.0 ] ); M=hcat(M, [ 5.0, 6.0 ] )      ## two 2 vectors
2×2 Array{Float64,2}:
3.0  5.0
4.0  6.0```

### Bottom-Append Row Vectors to Matrix

snippet.juliarepl
```julia> M= Matrix{Float64}(undef,0,2)                         ## empty 0×2 array: 2 columns, 0 rows
0×2 Array{Float64,2}

julia> [M;[1.0,2.0]']                                  ## ; means vcat. vcat= cat(1,...) , i.e., dimension 1
1×2 Array{Float64,2}:
1.0  2.0

julia> vcat( vcat( vcat(M, [1.0,2.0]'), [3.0,4.0]' ), [5.0,6.0]' )  ## push row vector (transposed column vecs)
3×2 Array{Float64,2}:
1.0  2.0
3.0  4.0
5.0  6.0```

## Multidimensional Iterators ("Comprehensions")

A sophisticated “comprehensions” versions of “for” can initialize multidimensional arrays with iterators:

snippet.juliarepl
```julia> [r * c for r in 1:5, c in 1:5]
5×5 Array{Int64,2}:
1   2   3   4   5
2   4   6   8  10
3   6   9  12  15
4   8  12  16  20
5  10  15  20  25```

Comprehensions are well explained in https://en.wikibooks.org/wiki/Introducing_Julia/Arrays_and_tuples#Creating_arrays_using_range_objects

“Generator Expressions” can omit the surrounding brackets, as in

snippet.juliarepl
```julia> sum(1/n^2 for n=1:10)
1.5497677311665408```

## Matrix Function Over Grid From Two Vectors (for 3D or Contourplots)

Julia makes it easy to create a rectangular grid with values. Here is how the grid will look like:

snippet.juliarepl
```julia> x= 1:5; y= 3:9; g(x,y)= (x,y);		## show off the grid that we will generate

julia> g.(x,y')			        	## note element-wise use with transpose
5×7 Array{Tuple{Int64,Int64},2}:
(1, 3)  (1, 4)  (1, 5)  (1, 6)  (1, 7)  (1, 8)  (1, 9)
(2, 3)  (2, 4)  (2, 5)  (2, 6)  (2, 7)  (2, 8)  (2, 9)
(3, 3)  (3, 4)  (3, 5)  (3, 6)  (3, 7)  (3, 8)  (3, 9)
(4, 3)  (4, 4)  (4, 5)  (4, 6)  (4, 7)  (4, 8)  (4, 9)
(5, 3)  (5, 4)  (5, 5)  (5, 6)  (5, 7)  (5, 8)  (5, 9)```

The prime on the y defines the transpose (or adjoint). Instead of a `g()` function that returns tuples to illustrate what will be used in each array cell, we can use a simple function:

snippet.juliarepl
```julia> x= 1:5; y= 3:9; f(x,y)= x^2*sqrt(y);	## apply a useful function

julia> f.(x,y')
5×7 Array{Float64,2}:
1.73205   2.0   2.23607   2.44949   2.64575   2.82843   3.0
6.9282    8.0   8.94427   9.79796  10.583    11.3137   12.0
15.5885   18.0  20.1246   22.0454   23.8118   25.4558   27.0
27.7128   32.0  35.7771   39.1918   42.332    45.2548   48.0
43.3013   50.0  55.9017   61.2372   66.1438   70.7107   75.0```

## Number of Dimensions, Elements, and Dimensions

snippet.juliarepl
```julia> M= [r * c for r in 1:3, c in 1:6]	## a sample two-dimensional comprehension
3×6 Array{Int64,2}:
1  2  3   4   5   6
2  4  6   8  10  12
3  6  9  12  15  18

julia> ndims(M)					## array as two dimensions
2

julia> size(M)				## number of rows (aka nrows, 3) and number of columns (aka ncols, 6)
(3, 6)

julia> size(M,2)			## the second dimension
6```

## Linear Indexing (Folded Matrix as Vector)

Julia can interpret matrices as long vectors, placed in column-row order. This is called linear indexing:

snippet.juliarepl
```julia> M= [r * c for r in 1:3, c in 1:6];	## the same two-dimensional matrix

julia>

julia> length(M)				## the total number of elements
18

julia> ( M, M, M )			## 5=(2,2); 8=(2,3); 11=(2,4)
( 4, 6, 8 )

julia> eachindex(M)				## An iterator (OneTo). (could be collect()ed)
Base.OneTo(18)

julia> eachindex(M) == 1:length(M)
true

julia> ( stride(M,1) , stride(M,2) )	## num indexes to advance in vector M to advance in matrix M direction 1 or 2
(1, 3)```

### Mixing Linear and Non-Linear Access

snippet.juliarepl
```julia> M= [ 1 2 3 4 ; 5 6 7 8 ]
2×4 Array{Int64,2}:
1  2  3  4
5  6  7  8

julia> M[:]                ## matrix linear access as vector
8-element Array{Int64,1}:
1
5
2
6
3
7
4
8

julia> M[:,1]              		## column (same as x[ [1,2] ])
2-element Array{Int64,1}:
1
5

julia> M[1,:]              		## row access
4-element Array{Int64,1}:
1
2
3
4

julia> M[2,3]              		## cell access
7

julia> M[ 2, [3,4] ]       		## arbitrary column and row access
2-element Array{Int64,1}:
7
8

julia> M[ :, end–1]        		## `end` is often easier
2-element Array{Int64,1}:
3
7

julia> M[ :, [false,true,false,true] ]  ## via Boolean vector
2×2 Array{Int64,2}:
2  4
6  8```

## REPL-Pastable Expression

snippet.juliarepl
```julia> M= [r * c for r in 1:3, c in 1:6];        ## the two-dimensional matrix

julia> @show(M);
M = [1 2 3 4 5 6; 2 4 6 8 10 12; 3 6 9 12 15 18]```

## Converting an X-Dimensional Array into a Different (e.g., Y-Dimensional) Array

To convert an existing 1-D array into a 2-D array, or reshape an existing array, use the `reshape` function:

snippet.juliarepl
```julia> reshape( 1:20 , (4, 5) )              ## sequence into 4×5 matrix
4×5 reshape(::UnitRange{Int64}, 4, 5) with eltype Int64:
1  5   9  13  17
2  6  10  14  18
3  7  11  15  19
4  8  12  16  20```

Usually, the type is ok, but you could also `convert( Matrix{Int64}, M )`. Reshape is an Alias: If you change a value in the reshaped array, the original array will still be affected, too. If this is not desired, work with a `copy(v)`.

snippet.juliarepl
```julia> M= [1,2,3,4,5,6];	## the original matrix M

julia> b= reshape(M, 2, 3)	## b is only an alias to M
2×3 Array{Int64,2}:
1  3  5
2  4  6

julia> b= 99;		## assigning to b changes M, too!

julia> M			## as you can see here
6-element Array{Int64,1}:
99
2
3
4
5
6```

## Creating Multidimensional Grids

Use the `product` function from `Iterators`:

snippet.juliarepl
```julia> using Base.Iterators

julia> arr= collect( product( [1,2,3], [4,5] ) )
3×2 Array{Tuple{Int64,Int64},2}:
(1, 4)  (1, 5)
(2, 4)  (2, 5)
(3, 4)  (3, 5)

julia> reshape( ans, (6,1) )
6×1 Array{Tuple{Int64,Int64},2}:
(1, 4)
(2, 4)
(3, 4)
(1, 5)
(2, 5)
(3, 5)```

Tuples were explained in Arrays and Tuples. To untupelize (convert) the resulting object into a matrix, use the reinterpret function (this is a bit tricky due to memory representation formats):

snippet.juliarepl
```julia> using Base.Iterators

julia> vt= collect( product( [1,2,3], [4,5] ) )  ## 6-element Array{Tuple{Int64,Int64},1}: vector of 2-tuples
3×2 Array{Tuple{Int64,Int64},2}:
(1, 4)  (1, 5)
(2, 4)  (2, 5)
(3, 4)  (3, 5)

julia> reshape( collect(flatten(vt)), (6,2) )
6×2 Array{Int64,2}:
1  1
4  5
2  2
4  5
3  3
4  5```
• It may not be necessary to materialize the grid with `collect`, which can save memory.

## Learning the Dimension of Higher-Dimensional Arrays

In Julia, rows come first, columns come second. However, this also works for more than two dimensions.

snippet.juliarepl
```julia> z= zeros(2,3,4)
2×3×4 Array{Float64,3}:
[:, :, 1] =
0.0  0.0  0.0
0.0  0.0  0.0

[:, :, 2] =
0.0  0.0  0.0
0.0  0.0  0.0

[:, :, 3] =
0.0  0.0  0.0
0.0  0.0  0.0

[:, :, 4] =
0.0  0.0  0.0
0.0  0.0  0.0

julia> size( z )
(2, 3, 4)```

The second argument to size selects the element from the results vector. Thus, `size( zeros(2,3,4), 2 ) == size( zeros(2,3,4) )`, which here is 3.

snippet.juliarepl
```julia> ncol( d ) = size( d, 2 )
ncol (generic function with 1 method)

julia> ncol( zeros(2,3) )
3```

## Mapping Over Matrices (Applying Operations f(X,Y) -> Z)

### For Loops

The standard code is:

snippet.juliarepl
```julia> f(x,y) = sqrt( x^2 + y^2 );

julia> xset= 10:11; yset= 12:13; i= 1;

julia> m= fill( NaN, length(xset)*length(yset), 3 )

julia> i
1

julia> for x in xset
for y in yset
global i
m[i,:]=	[ x, y, f(x,y) ]
i+=1
end#for y#
end#for x##

julia> m
4×3 Array{Float64,2}:
10.0  12.0  15.6205
10.0  13.0  16.4012
11.0  12.0  16.2788
11.0  13.0  17.0294```

Possible speed enhancements:

• use `Matrix{Float64}(undef, length(xset)*length(yset), 3);` to work with unitialized matrix
• place `@inbounds` before the `for` statement tells Julia not to check the bounds. However, this risks crashing Julia and probably works only inside functions.
• use `@simd @inbounds` before the `for` statement to ask for modestly advanced Intel CPU vector features.

### Compact Alternatives

snippet.julia
```using Base.Iterators

reduce(vcat, [[x, y, f(x, y)]' for (y, x) in Iterators.product(yset, xset)])```

or

snippet.julia
`hcat([ [z(x,y) for x in xset for y in yset] for z in ((x,y)->x, (x,y)->y, f) ]...)`

or

snippet.julia
```using StaticArrays

vv = [ @SVector [x, y, f(x,y)]  for x in xset for y in yset ]

reinterpret(Float64, vv, (3, length(vv) ))```

or

snippet.julia
```using Base.Iterators
reduce(hcat, [[x, y, f(x, y)] for (y, x) in Iterators.product(yset, xset)]) ## 8.531 μs```

dot_postfix_functions operate on scalars inside arrays. To operate on arrays inside arrays (e.g., vectors inside matrices),

snippet.juliarepl
```julia> broadcast( / , [1 2 3 4 ; 5 6 7 8], [10,20] )	## shorthand for [1 2 3 4 ; 5 6 7 8 ] ./ [10,20]
2×4 Array{Float64,2}:
0.1   0.2  0.3   0.4
0.25  0.3  0.35  0.4

julia> broadcast( / , [1 2 3 4 ; 5 6 7 8], [10,20] )	==  [1 2 3 4 ; 5 6 7 8 ] ./ [10,20]
true

julia> broadcast( / , [10,20], [1 2 3 4 ; 5 6 7 8] )	## shorthand for [10,20] ./ [1 2 3 4 ; 5 6 7 8 ]
2×4 Array{Float64,2}:
10.0  5.0      3.33333  2.5
4.0  3.33333  2.85714  2.5

julia> broadcast( / , [10,20], [1 2 3 4 ; 5 6 7 8] )	== [10,20] ./ [1 2 3 4 ; 5 6 7 8 ]
true```

The dot function versions are faster than the broadcasts, because the dot functions know how to “fuse” operations. That is, they do not create temporary intermediate matrices—they fuse all the necessary operations into each destination cell.

## Removing Sets of Rows Based on Condition (e.g., containing at least one NaN)

snippet.juliarepl
```julia> M= collect(reshape( 1.0:18.0, 6, 3 ));  M[2,3]=NaN;  M[4,2]=NaN;  ## rows 4 and 6 are bad

julia> M
6×3 Array{Float64,2}:
1.0    7.0   13.0
2.0    8.0  NaN
3.0    9.0   15.0
4.0  NaN     16.0
5.0   11.0   17.0
6.0   12.0   18.0

julia> goodobs= map( i ->(!any(isnan, M[i,:])), 1:size(M,1) )   ## Method 1
6-element Array{Bool,1}:
true
false
true
false
true
true

julia> M[ goodobs , : ]
4×3 Array{Float64,2}:
1.0   7.0  13.0
3.0   9.0  15.0
5.0  11.0  17.0
6.0  12.0  18.0

julia> goodobs= filter( i->(all(!isnan, M[i,:])), 1:size(M,1) )		## Method 2
4-element Array{Int64,1}:
1
3
5
6

julia> M[ goodobs , : ]
4×3 Array{Float64,2}:
1.0   7.0  13.0
3.0   9.0  15.0
5.0  11.0  17.0
6.0  12.0  18.0```

## Matrix Overall, Columnwise, and Rowwise Statistics

Some statistics work over matrices. For example,

snippet.juliarepl
```julia> using Statistics

julia> M= collect(reshape( 1.0:15.0, 3, 5 ))
3×5 Array{Float64,2}:
1.0  4.0  7.0  10.0  13.0
2.0  5.0  8.0  11.0  14.0
3.0  6.0  9.0  12.0  15.0

julia> ( minimum(M), mean(M), maximum(M) )   ## note: not min() and max(), which operate over their args instead
(1.0, 8.0, 15.0)

julia> mean( M; dims=1 )              ## column-wise row means
1×5 Array{Float64,2}:
2.0  5.0  8.0  11.0  14.0

julia> std( M; dims=2 )                  ## row-wise column standard deviations
3×1 Array{Float64,2}:
4.743416490252569
4.743416490252569
4.743416490252569```

## Matrix Math (Linear Algebra)

• Not shown: `inv(A)`, `trace(A)`, `det(A)`, `eigvecs(A)`, `eigvals(A)`, `factorize(A)`. Many of these, esp factorize, can be clever in exploting aspects of its input matrix.

For matrix multiplication, use the standard `*`. For element-by-element multiplication, use `.*`. Powering matrices also works as expected (either matrix with `^`, or element-by-element with `.^`).

snippet.juliarepl
```julia> using LinearAlgebra

julia> ones(2,2)^4          ## matrix multiplication, repeated
2×2 Array{Float64,2}:
8.0  8.0
8.0  8.0

julia> ones(2,2).^4         ## element multiplicated, repeated
2×2 Array{Float64,2}:
1.0  1.0
1.0  1.0```

### Sequential Matrix Multiplications

Consider multiplying four matrices with dimensions `(1x10) * (10x10) * (10x10) * (10x1)`. If you calculate `(1x10) * (10x10)` and `(10x10) * (10x1)` first, you need less storage and fewer calculations than if you calculate `(10x10) * (10x10)` first. Although mathematically identical, the order in which you do calculations matters.

Fortunately, there are two packages that automatically find the speediest way to do sequential matrix multiplications:

# Backmatter

• Julia has built-in support for sparse arrays, `SparseMatrix`. Use `dropzeros`. The implied non-entered value is zero.
• Julia has built-in support for a variety of special matrices, such as Triangular or Symmetric ones.
• Julia can easily deal with more than 2-dimensions in arrays (tensors).

### Page Tools 