text:
module MarkdownLiteral

import HypertextLiteral, CommonMark

"""
```julia
@markdown(""\"
# MarkdownLiteral.jl

The macro `@markdown` lets you write [Markdown](https://www.markdownguide.org/getting-started/) inside Pluto notebooks. *Here is an example:*
""\")
```

You can also use the macro to write HTML!

```julia
@markdown(""\"
<p>
	The macro <code>@markdown</code> lets you write <a href="https://developer.mozilla.org/docs/Web/HTML">HTML</a> inside Pluto notebooks.
	<em>Here is an example:</em>
</p>
""\")
```

## Interpolation

You can unlock superpowers by combining `@markdown` with **interpolation** (using `\$`). For our example, let's create some data:

```julia
films = [
	(title="Frances Ha", director="Noah Baumbach", year=2012)
	(title="Portrait de la jeune fille en feu", director="Céline Sciamma", year=2019)
	(title="De noorderlingen", director="Alex van Warmerdam", year=1992)
]
```
Now, we can use *interpolation* to display our data:
```julia
@markdown(""\"
My films:
\$([
	"- **\$(f.title)** (\$(f.year)) by _\$(f.director)_\n"
	for f in films
])
""\")
```

This gives us:

> My films:
> - **Frances Ha** (2012) by _Noah Baumbach_
> - **Portrait de la jeune fille en feu** (2019) by _Céline Sciamma_
> - **De noorderlingen** (1992) by _Alex van Warmerdam_

Alternatively, you could write this using HTML instead of Markdown (*with the same macro!*):

```julia
@markdown(""\"
<p>My films:</p>
<ul>
\$([
	@markdown("<li>
		<b>\$(f.title)</b> (\$(f.year)) by <em>\$(f.director)</em>
	</li>")
	for f in films
])
</ul>
""\")
```

## Advanced interpolation

Because interpolation is powered by [HypertextLiteral.jl](https://github.com/MechanicalRabbit/HypertextLiteral.jl), you can use advanced features:
- Interpolated attributes are automatically escaped
- You can use a `NamedTuple` or `Dict` for the CSS `style` attribute
- Interpolating Julia objects into a `<script>` will automatically convert to JavaScript code(!)

For
```julia
logs = [
	(text="Info", urgent=false),
	(text="Alert", urgent=true),
	(text="Update", urgent=false),
]
```

```julia
@markdown("\$((
	@markdown("<div style=\$((
        font_weight=900,
		padding=".5em",
		background=log.urgent ? "pink" : "lightblue",
	))>\$(log.text)</div>")
	for log in logs
))")
```
Result:

![](https://user-images.githubusercontent.com/6933510/146623300-316e5a17-2daf-43ed-b70c-6c33278faf32.png)
"""
macro markdown(expr)
    cm_parser = CommonMark.Parser()
    CommonMark.enable!(cm_parser, [
        CommonMark.AdmonitionRule(),
        CommonMark.AttributeRule(),
        CommonMark.AutoIdentifierRule(),
        CommonMark.CitationRule(),
        CommonMark.FootnoteRule(),
        CommonMark.MathRule(),
        CommonMark.RawContentRule(),
        CommonMark.TableRule(),
        CommonMark.TypographyRule(),
    ])
    quote
        result = $(esc(Expr(:macrocall, getfield(HypertextLiteral, Symbol("@htl")), __source__, expr)))
        htl_output = repr(MIME"text/html"(), result)

        $(cm_parser)(htl_output)
    end
end

# Some aliases, it's up to you which one to import.

var"@markdownliteral" = var"@markdown"
var"@mdx" = var"@markdown"

end


text:
using Gurobi

using Compat
using Compat.Test, Compat.SparseArrays, Compat.Random

@testset "C API" begin
    include("c_wrapper.jl")
end

const mpb_tests = [
    "lp_01a",
    "lp_01b",
    "lp_02",
    "lp_03",
    "lp_04",
    "mip_01",
    "qp_01",
    "qp_02",
    "qcqp_01",
    "mathprog",
    "test_grb_attrs",
    "env",
    "range_constraints",
    "test_get_strarray",
    "large_coefficients",
    "multiobj",
    "test_read"
]

@testset "MathProgBase Tests" begin
    for file in mpb_tests
        evalfile(joinpath("MathProgBase", "$(file).jl"))
    end
end

@testset "MathOptInterface Tests" begin
    include("MOIWrapper.jl")
end

ext: jl


text:
"""
Near field solution for a crack in a rectilinear anisotropic elastic medium.
See:
G. C. Sih, P. C. Paris and G. R. Irwin, Int. J. Frac. Mech. 1, 189 (1965)

Ported from `matscipy.fracture_mechanics.crack` module by James Kermode
"""
module Fracture

using LinearAlgebra
using Polynomials: Polynomial, roots
using Einsum
using JuLIP
using JuLIPMaterials.CLE

export RectilinearAnisotropicCrack, PlaneStrain, PlaneStress
export displacements, stresses, deformation_gradient, k1g, u_CLE

abstract type Crack end

abstract type StressState end
struct PlaneStress <: StressState end
struct PlaneStrain <: StressState end

abstract type CoordinateSystem end

struct Cartesian <: CoordinateSystem
    x::AbstractVector
    y::AbstractVector
end

struct Cylindrical <: CoordinateSystem
    r::AbstractVector
    θ::AbstractVector
end

Cylindrical(cart::Cartesian) = Cylindrical(sqrt.(cart.x.^2 + cart.y.^2), atan.(cart.y, cart.x))

Cartesian(cyl::Cylindrical) = Cartesian(cyl.r * cos.(cyl.θ), cyl.r * sin.(cyl.θ))

struct RectilinearAnisotropicCrack{T <: Number} <: Crack
    a22::T
    μ1::Complex{T}
    μ2::Complex{T}
    p1::Complex{T}
    p2::Complex{T}
    q1::Complex{T}
    q2::Complex{T}
    inv_μ1_μ2::Complex{T}
    μ1_p2::Complex{T}
    μ2_p1::Complex{T}
    μ1_q2::Complex{T}
    μ2_q1::Complex{T}
end

function RectilinearAnisotropicCrack(a22::T, μ1::Complex{T}, μ2::Complex{T}, p1::Complex{T}, p2::Complex{T},
                                     q1::Complex{T}, q2::Complex{T}) where {T}
    inv_μ1_μ2 = 1/(μ1 - μ2)
    μ1_p2 = μ1 * p2
    μ2_p1 = μ2 * p1
    μ1_q2 = μ1 * q2
    μ2_q1 = μ2 * q1
    return RectilinearAnisotropicCrack(a22, μ1, μ2, p1, p2, q1, q2, inv_μ1_μ2,
                                       μ1_p2, μ2_p1, μ1_q2, μ2_q1)
end

function RectilinearAnisotropicCrack(a11::T, a22::T, a12::T, a16::T, a26::T, a66::T) where {T}
    p = Polynomial(reverse([ a11, -2*a16, 2*a12 + a66, -2*a26, a22 ] )) # NB: opposite order to np.poly1d()
    μ1, μ1s, μ2, μ2s = roots(p)
    (μ1 != conj(μ1s) ||  μ2 != conj(μ2s)) && error("Roots not in pairs.")

    p1 = a11 * μ1^2 + a12 - a16 * μ1
    p2 = a11 * μ2^2 + a12 - a16 * μ2

    q1 = a12 * μ1 + a22 / μ1 - a26
    q2 = a12 * μ2 + a22 / μ2 - a26

    return RectilinearAnisotropicCrack(a22, μ1, μ2, p1, p2, q1, q2)
end

RectilinearAnisotropicCrack(::PlaneStress, S::AbstractMatrix) = RectilinearAnisotropicCrack(S[1, 1], S[2, 2], S[1, 2], S[1, 6], S[2, 6], S[6, 6])

function RectilinearAnisotropicCrack(::PlaneStrain, S::AbstractMatrix)
    b11, b22, b33, b12, b13, b23, b16, b26, b36, b66 = (S[1, 1], S[2, 2], S[3, 3], S[1, 2], S[1, 3],
                                                        S[2, 3], S[1, 6], S[2, 6], S[3, 6], S[6, 6])
    a11 = b11 - (b13 * b13) / b33
    a22 = b22 - (b23 * b23) / b33
    a12 = b12 - (b13 * b23) / b33
    a16 = b16 - (b13 * b36) / b33
    a26 = b26 - (b23 * b36) / b33
    a66 = b66
    return RectilinearAnisotropicCrack(a11, a22, a12, a16, a26, a66)
end

function RectilinearAnisotropicCrack(ss::StressState, C::AbstractMatrix, crack_surface::AbstractVector, crack_front::AbstractVector)
    A = rotation_matrix(y=crack_surface, z=crack_front)
    Crot = rotate_moduli(C, A)
    S = inv(Crot)
    return RectilinearAnisotropicCrack(ss, S)
end

RectilinearAnisotropicCrack(ss::StressState,
                            C11, C12, C44,
                            crack_surface::AbstractVector, crack_front::AbstractVector) =
                                RectilinearAnisotropicCrack(ss, voigt_moduli(C11, C12, C44), crack_surface, crack_front)

"""
Displacement field in mode I fracture.

Parameters
----------
r : array
    Distances from the crack tip.
theta : array
    Angles with respect to the plane of the crack.

Returns
-------
u : array
    Displacements parallel to the plane of the crack.
v : array
    Displacements normal to the plane of the crack.
"""
function displacements(crack::RectilinearAnisotropicCrack, cyl::Cylindrical)
    h1 = sqrt.(2.0 * cyl.r / π)
    h2 = sqrt.( cos.(cyl.θ) .+ crack.μ2 * sin.(cyl.θ) )
    h3 = sqrt.( cos.(cyl.θ) .+ crack.μ1 * sin.(cyl.θ) )

    u = h1 .* real.( crack.inv_μ1_μ2 * ( crack.μ1_p2 * h2 - crack.μ2_p1 * h3 ) )
    v = h1 .* real.( crack.inv_μ1_μ2 * ( crack.μ1_q2 * h2 - crack.μ2_q1 * h3 ) )

    return u, v
end

"""
Deformation gradient tensor in mode I fracture.

Parameters
----------
r : array_like
    Distances from the crack tip.
theta : array_like
    Angles with respect to the plane of the crack.

Returns a Jacobian ```[du_dx du_dy
                       dv_dx dv_dz]```

where

* `du_dx` - Derivatives of displacements parallel to the plane within the plane.
* `du_dy`- Derivatives of displacements parallel to the plane perpendicular to the plane.
* `dv_dx`- Derivatives of displacements normal to the plane of the crack within the plane.
* `dv_dy` - Derivatives of displacements normal to the plane of the crack perpendicular to the plane.
"""
function deformation_gradient(crack::RectilinearAnisotropicCrack, cyl::Cylindrical)
    f = 1 ./ sqrt.(2 * π * cyl.r)

    h1 = (crack.μ1 * crack.μ2) * crack.inv_μ1_μ2
    h2 = sqrt.( cos.(cyl.θ) + crack.μ2 * sin.(cyl.θ) )
    h3 = sqrt.( cos.(cyl.θ) + crack.μ1 * sin.(cyl.θ) )

    du_dx = f .* real.( crack.inv_μ1_μ2 * ( crack.μ1_p2 ./ h2 - crack.μ2_p1 ./ h3 ) )
    du_dy = f .* real.( h1 * ( crack.p2 ./ h2 - crack.p1 ./ h3 ) )

    dv_dx = f .* real.( crack.inv_μ1_μ2 * ( crack.μ1_q2 ./ h2 - crack.μ2_q1 ./ h3 ) )
    dv_dy = f .* real.( h1 * ( crack.q2 ./ h2 - crack.q1 ./ h3 ) )

    return [[du_dx] [du_dy]
            [dv_dx] [dv_dy]]
end

"""
Stress field in mode I fracture.

Parameters
----------
r : array
    Distances from the crack tip.
theta : array
    Angles with respect to the plane of the crack.

Returns
-------
sig_x : array
    Diagonal component of stress tensor parallel to the plane of the
    crack.
sig_y : array
    Diagonal component of stress tensor normal to the plane of the
    crack.
sig_xy : array
    Off-diagonal component of the stress tensor.
"""
function stresses(crack::RectilinearAnisotropicCrack, cyl::Cylindrical)
    f = 1 ./ sqrt.(2.0 * π * cyl.r)

    h1 = (crack.μ1 * crack.μ2) * crack.inv_μ1_μ2
    h2 = sqrt.( cos.(cyl.θ) + crack.μ2 * sin.(cyl.θ) )
    h3 = sqrt.( cos.(cyl.θ) + crack.μ1 * sin.(cyl.θ) )

    sig_x  = f * real.(h1 * (crack.μ2 ./ h2 - crack.μ1 ./ h3))
    sig_y  = f * real.(crack.inv_μ1_μ2 * (crack.μ1 ./ h2 - crack.μ2 ./ h3))
    sig_xy = f * real.(h1 * (1 ./ h3 - 1 ./ h2))

    return sig_x, sig_y, sig_xy
end

"""
K1G, Griffith critical stress intensity in mode I fracture
"""
function k1g(crack::RectilinearAnisotropicCrack, surface_energy)
    return sqrt(abs(4 * surface_energy /
                     imag(crack.a22 * ((crack.μ1 + crack.μ2) / (crack.μ1  * crack.μ2 )))))
end

# Cartesian coordinate convenience wrappers

displacements(crack::RectilinearAnisotropicCrack, cart::Cartesian) = displacements(crack, Cylindrical(cart))
displacements(crack::RectilinearAnisotropicCrack, x::AbstractVector, y::AbstractVector) = displacements(crack, Cartesian(x, y))

deformation_gradient(crack::RectilinearAnisotropicCrack, cart::Cartesian) = deformation_gradient(crack, Cylindrical(cart))
deformation_gradient(crack::RectilinearAnisotropicCrack, x::AbstractVector, y::AbstractVector) = deformation_gradient(crack, Cartesian(x, y))

stresses(crack::RectilinearAnisotropicCrack, cart::Cartesian) = stresses(crack, Cylindrical(cart))
stresses(crack::RectilinearAnisotropicCrack, x::AbstractVector, y::AbstractVector) = stresses(crack, Cartesian(x, y))

function u_CLE(crack::RectilinearAnisotropicCrack, x::AbstractVector, y::AbstractVector)
    function u(k, α)
        ux, uy = displacements(crack, x .- α, y)
        return k * [ux uy]'
    end

    function du(k, α)
        D = deformation_gradient(crack, x .- α, y)
        return  -k * [D[1,1] D[2,1]]'
    end

    return u, du
end

function u_CLE(crack::RectilinearAnisotropicCrack, crystal::AbstractAtoms{T}, x0::T, y0::T) where {T}
    X = positions(crystal) |> mat
    return u_CLE(crack, X[1, :] .- x0, X[2, :] .- y0)
end

end


text:
using Printf: @sprintf

const superscriptNum = Dict(['0'=>'⁰', '1'=>'¹', '2'=>'²', '3'=>'³', '4'=>'⁴', '5'=>'⁵',
                             '6'=>'⁶', '7'=>'⁷', '8'=>'⁸', '9'=>'⁹'])
const subscriptNum   = Dict(['0'=>'₀', '1'=>'₁', '2'=>'₂', '3'=>'₃', '4'=>'₄', '5'=>'₅',
                             '6'=>'₆', '7'=>'₇', '8'=>'₈', '9'=>'₉'])
const superscriptSym = Dict(['+'=>'⁺', '-'=>'⁻', '('=>'⁽', ')'=>'⁾', '!'=>'ꜝ'])

"""

    checkFname(Fname::String; showWarning::Bool=true) -> String

Check if there is a file with the same name in the current directory. If so, will add an
`"_N"` at the end of the file name `String`. `showWarning` determines whether prints out
the WARNING info when there is a file with the same name.
"""
function checkFname(Fname::String; showWarning::Bool=true)
    FnameN = Fname
    while isfile(FnameN) == true
        i=0
        contains(FnameN, ".") ? i=findlast(".", FnameN)|>last : i=(FnameN|>length)+1
        FnameN = FnameN[1:i-1]*"_N"*FnameN[i:end]
    end
    FnamePrint, FnameNPrint = map([Fname, FnameN]) do f
        contains(f, "/") ? f[(findlast("/", f) |> last)+1 : end] : f
    end
    FnameN != Fname && showWarning && (@warn """The file expected to create already exists.
                                                Adding another suffix.
                                                Old name: $(FnamePrint)
                                                New name: $(FnameNPrint)""")
    FnameN
end

function advancedParse(content::AbstractString,
                       ParseFunc::F=adaptiveParse) where {F<:Function}
    res = ParseFunc(content)
    res === nothing && (res = content)
    res
end

function adaptiveParse(content::AbstractString)
    res = tryparse(Int, content)
    res === nothing && (res = tryparse(Float64, content))
    res === nothing && (res = tryparse(Complex{Float64}, content))
    res
end

function numToSups(num::Int)
    str = string(num)
    [superscriptNum[i] for i in str] |> prod
end

function numToSubs(num::Int)
    str = string(num)
    [subscriptNum[i] for i in str] |> prod
end

function alignNum(x::Number, lpadN::Int=8, rpadN::Int=21; roundDigits::Int=-1)
    if roundDigits < 0
        str = x |> string
    else
        format = "%.$(roundDigits)f"
        str = :(@sprintf($format, $x)) |> eval
    end
    body = split(str, '.')
    if length(body) == 2
        head, tail = body
        tail = "."*tail
        rpadN += 1
    else
        head = body[]
        tail = ""
    end
    lpad(head, lpadN) * rpad(tail, rpadN)
end

function alignNumSign(c::Real; roundDigits::Int=-1)
    if c < 0
        alignNum(c, 0, 0; roundDigits)
    else
        " "*alignNum(c, 0, 0; roundDigits)
    end
end


text:
cd(@__DIR__)
using Revise, DifferentialEvolutionMCMC, Random, Distributions

Random.seed!(50514)

function prior_loglike(μ, σ)
    LL = 0.0
    LL += logpdf(Normal(0, 1), μ)
    LL += logpdf(truncated(Cauchy(0, 1), 0, Inf), σ)
    return LL
end

function sample_prior()
    μ = rand(Normal(0, 1))
    σ = rand(truncated(Cauchy(0, 1), 0, Inf))
    return [μ,σ]
end

bounds = ((-Inf,Inf),(0.0,Inf))

data = rand(Normal(0.0, 1.0), 50)

function loglike(data, θ...)
    μ, σ = θ
    return sum(logpdf.(Normal(μ, σ), data))
end

names = (:μ,:σ)

model = DEModel(;
    sample_prior,
    prior_loglike,
    loglike,
    data,
    names
)

de = DE(;sample_prior, bounds, burnin=1000, Np=6)
n_iter = 2000
chains = sample(model, de, MCMCThreads(), n_iter, progress=true)


text: Dict{String, Any}("delete" => Dict{String, Any}("value" => "This string has a \x7f delete control code.", "type" => "string"), "formfeed" => Dict{String, Any}("value" => "This string has a \f form feed character.", "type" => "string"), "notunicode2" => Dict{String, Any}("value" => "This string does not have a unicode \\u escape.", "type" => "string"), "unitseparator" => Dict{String, Any}("value" => "This string has a \x1f unit separator control code.", "type" => "string"), "backslash" => Dict{String, Any}("value" => "This string has a \\ backslash character.", "type" => "string"), "notunicode3" => Dict{String, Any}("value" => "This string does not have a unicode \\u0075 escape.", "type" => "string"), "notunicode4" => Dict{String, Any}("value" => "This string does not have a unicode \\u escape.", "type" => "string"), "tab" => Dict{String, Any}("value" => "This string has a \t tab character.", "type" => "string"), "carriage" => Dict{String, Any}("value" => "This string has a \r carriage return character.", "type" => "string"), "quote" => Dict{String, Any}("value" => "This string has a \" quote character.", "type" => "string"), "newline" => Dict{String, Any}("value" => "This string has a \n new line character.", "type" => "string"), "notunicode1" => Dict{String, Any}("value" => "This string does not have a unicode \\u escape.", "type" => "string"), "backspace" => Dict{String, Any}("value" => "This string has a \b backspace character.", "type" => "string"))


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Generates a Linear Programming model for solving an L1 SVM binary
# classification problem from the LIBSVM repository.

# Sample usage:
# $ julia --project=. generate_l1_svm_lp.jl --input_filename=../data/duke.tr
#     --output_filename=/tmp/duke.mps.gz --regularizer_weight=1.5

# The LP formulation for L1 SVM is taken from equation (5) in:
# Ji Zhu, Saharon Rosset, Trevor Hastie, and Rob Tibshirani. 2003. 1-norm
# support vector machines. In Proceedings of the 16th International Conference
# on Neural Information Processing Systems (NIPS'03). MIT Press, Cambridge, MA,
# USA, 49–56.
# https://papers.nips.cc/paper/2003/file/49d4b2faeb4b7b9e745775793141e2b2-Paper.pdf

import ArgParse
import JuMP
import SparseArrays
import LinearAlgebra

const SparseMatrixCSC = SparseArrays.SparseMatrixCSC
const sparse = SparseArrays.sparse
const norm = LinearAlgebra.norm
const Diagonal = LinearAlgebra.Diagonal
const nzrange = SparseArrays.nzrange

mutable struct SvmTrainingData
  feature_matrix::SparseMatrixCSC{Float64,Int64}
  labels::Vector{Float64}
end

"""
Populates 'model' with the linear program model for solving the L1 SVM
problem with training set (X, y). The model is given by
  min  sum w[i] + regularizer_weight * sum z[i]
  s.t. w >= 1 - diag(y) * [X 1] * beta
       w >= 0
       z >= beta
       z >= -beta
"""
function populate_libsvm_model(
  model::JuMP.Model,
  data::SvmTrainingData,
  regularizer_weight::Float64,
)
  n, d = size(data.feature_matrix)
  println("Generating a model with $n datapoints and $(d - 1) features.")
  JuMP.@variable(model, beta[i = 1:d])
  JuMP.@variable(model, w[i = 1:n], lower_bound = 0.0)
  JuMP.@variable(model, z[i = 1:d])
  JuMP.@objective(model, Min, sum(w) + regularizer_weight * sum(z))
  JuMP.@constraint(model, z .>= beta)
  JuMP.@constraint(model, z .>= -beta)
  JuMP.@constraint(
    model,
    w .>= (1.0 .- data.labels .* (data.feature_matrix * beta))
  )
  return model
end

"""
Defines parses and args.

# Returns
A dictionary with the values of the command-line arguments.
"""
function parse_command_line()
  arg_parse = ArgParse.ArgParseSettings()

  ArgParse.@add_arg_table! arg_parse begin
    "--input_filename"
    help = "The .train file containing the problem data."
    arg_type = String
    required = true

    "--output_filename"
    help = "Filename for the output .mps (or .mps.gz) model."
    arg_type = String
    required = true

    "--regularizer_weight"
    help = "Weight of the L1 regularizer."
    arg_type = Float64
    required = true
  end

  return ArgParse.parse_args(arg_parse)
end

"""
Loads a LIBSVM file and returns the loaded SvmTrainingData struct.
"""
function load_libsvm_file(file_name::String)
  open(file_name, "r") do io
    labels = Vector{Float64}()
    row_indices = Vector{Int64}()
    col_indices = Vector{Int64}()
    matrix_values = Vector{Float64}()

    row_index = 0
    found_label_one = false
    for line in eachline(io)
      row_index += 1
      split_line = split(line)

      label = parse(Float64, split_line[1])
      # This ensures that labels are 1 or -1. Different datasets use {-1, 1}, {0, 1}, and {1, 2}.
      if label == 1.0
        found_label_one = true
      else
        label = -1.0
      end
      push!(labels, label)

      for i in 2:length(split_line)
        push!(row_indices, row_index)
        matrix_coef = split(split_line[i], ":")
        push!(col_indices, parse(Int64, matrix_coef[1]))
        push!(matrix_values, parse(Float64, matrix_coef[2]))
      end
    end
    @assert found_label_one
    feature_matrix = sparse(row_indices, col_indices, matrix_values)
    return SvmTrainingData(feature_matrix, labels)
  end
end

function normalize_columns(feature_matrix::SparseMatrixCSC{Float64,Int64})
  norm_of_columns = vec(sqrt.(sum(t -> t^2, feature_matrix, dims = 1)))
  norm_of_columns[iszero.(norm_of_columns)] .= 1.0
  return feature_matrix * Diagonal(1.0 ./ norm_of_columns)
end

function remove_empty_columns(feature_matrix::SparseMatrixCSC{Float64,Int64})
  keep_cols = Vector{Int64}()
  for j in 1:size(feature_matrix, 2)
    if length(nzrange(feature_matrix, j)) > 0
      push!(keep_cols, j)
    end
  end
  return feature_matrix[:, keep_cols]
end

function add_intercept(feature_matrix::SparseMatrixCSC{Float64,Int64})
  return [sparse(ones(size(feature_matrix, 1))) feature_matrix]
end

function preprocess_training_data(result::SvmTrainingData)
  result.feature_matrix = remove_empty_columns(result.feature_matrix)
  result.feature_matrix = add_intercept(result.feature_matrix)
  result.feature_matrix = normalize_columns(result.feature_matrix)
  return result
end

function main()
  parsed_args = parse_command_line()

  filename = parsed_args["output_filename"]
  model = JuMP.Model()
  regularizer_weight = parsed_args["regularizer_weight"]
  input_filename = parsed_args["input_filename"]
  data = load_libsvm_file(input_filename)
  data = preprocess_training_data(data)
  populate_libsvm_model(model, data, regularizer_weight)
  JuMP.write_to_file(model, parsed_args["output_filename"])
end

main()


text:
# export names
# do typealiases have to be exported?  I don't think so

typealias comm_type MPI.CComm

export PETSC_NULL, PETSC_IGNORE, PETSC_DECIDE, PETSC_DETERMINE, PETSCDEFAULT, PETSC_COMM_SELF

export PetscInt, PetscScalar, PetscBool, PetscErrorCode, PetscDataType, PetscReal

export PetscIntNullArray
export PETSC_INSERT_VALUES, PETSC_ADD_VALUES, PETSC_COPY_VALUESA
export PETSC_TRUE, PETSC_FALSE
export NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY, NORM_MAX

export PETSC_MAT_FLUSH_ASSEMBLY, PETSC_MAT_FINAL_ASSEMBLY

export KSPRICHARDSON, KSPCHEBYSHEV, KSPCG, KSPGROPPCG, KSPIPECG, KSPCGNE, KSPNASH, KSPSTCG, KSPGLTR, KSPFCG, KSPGMRES, KSPFGMRES, KSPLGMRES, KSPDGMRES, KSPPGMRES, KSPTCQMR, SKPBCGS, KSPIBCGS, KSPFBCGS, KSPFBCGSR, KSPBCGSL, KSPCGS, KSPTFQMR, KSPCR, KSPPIPECR, KSPSLQR, SKPPREONLY, KSPQCQ, KSPBICG, KSPMINRES, KSPSYMMLQ, KSPLCD, KSPPYTHON, KSPCR

export KSP_NORM_DEFAULT, KSP_NORM_NONE, KSP_NORM_PRECONDITIONED, KSP_NORM_UNPRECONGDITIONED, KSP_NORM_NATURAL

export KSP_NORM_MAX

export VECSEQ
export VECMPI
export VECSTANDARD
export VECSHARED
export VECSEQCUSP
export VECMPICUSP
export VECCUSP
export VECSEQVIENNACL
export VECMPIVIENNACL
export VECVIENNACL
export VECNEST
export VECSEQPTHREAD
export VECMPIPTHREAD
export VECPTHREAD

# decide which version of Petsc to use
# if no PETSC_DIR or PETSC_ARCH defined, use the one (hopefully) built
# during package installation
if !haskey(ENV, "PETSC_DIR") && !haskey(ENV, "PETSC_ARCH")
  file_path = joinpath(Pkg.dir("PETSc"), "deps/petsc_evars")
  args = open(readdlm, file_path)
  println("args = ", args)
  ENV["PETSC_DIR"] = args[1]
  ENV["PETSC_ARCH"] = args[2]
end

global const PETSC_DIR = ENV["PETSC_DIR"]
global const PETSC_ARCH = ENV["PETSC_ARCH"]

#PETSC_DIR = readall(`echo $PETSC_DIR`)
#PETSC_ARCH = readall(`echo $PETSC_ARCH`)

#PETSC_DIR = "/home/jared/build/petsc-3.6.0"
#PETSC_ARCH = "arch-linux2-c-debug"

#PETSC_DIR = getenv("PETSC_DIR");
#PETSC_ARCH = getenv("PETSC_ARCH");
#=
if (length(PETSC_DIR) == 0)
  disp("Must have environmental variable PETSC_DIR set")
end
if (length(PETSC_ARCH) == 0)
  disp("Must have environmental variable PETSC_ARCH set")
end
=#

global const libpetsclocation = string(PETSC_DIR, "/", PETSC_ARCH, "/lib/", "libpetsc")
global const petsc = libpetsclocation # for compatability with auto generated wrappers
global const libpetsc = Libdl.dlopen(libpetsclocation)

# definitions of Petsc types
# a way to automatically generate these would be preferable
# define all the data types that are known a priori
typealias Petsc64bitInt Int64
#typealias PetscBLASInt Int32
#typealias PetscScalar Cint
typealias PetscBool UInt32
typealias PetscDataType Cint  # C enums are Int32

typealias PetscLogDouble Cdouble
typealias PetscErrorCode Cint
#=
const PETSC_PRECISION_SINGLE = (UInt32)(4)
const PETSC_PRECISION_DOUBLE = (UInt32)(8)
=#

const PETSC_FALSE = (UInt32)(0)
const PETSC_TRUE = (UInt32)(1)

const PETSC_INT = (Int32)(0)
const PETSC_DOUBLE = (Int32)(1)
const PETSC_COMPLEX = (Int32)(2)
const PETSC_LONG = (Int32)(3)
const PETSC_SHORT = (Int32)(4)
const PETSC_FLOAT = (Int32)(5)
const PETSC_CHAR = (Int32)(6)
const PETSC_BIT_LOGICAL = (Int32)(7)
const PETSC_ENUM = (Int32)(8)
const PETSC_BOOL = (Int32)(9)
const PETSC___FLOAT128 = (Int32)(10)
const PETSC_OBJECT = (Int32)(11)
const PETSC_FUNCTION = (Int32)(12)
const PETSC_STRING = (Int32)(12)

# these functions are used for figuring out the sizes of the datatypes
# use the ones in PETSc.jl for writing programs
function PetscDataTypeFromString_(name::AbstractString)
    ptype = Array(Cint, 1)
    found = Array(PetscBool, 1)
    ccall((:PetscDataTypeFromString,petsc),PetscErrorCode,(Cstring,Ptr{PetscDataType},Ptr{PetscBool}), name, ptype, found)

    return ptype[1], convert(Bool, found[1])
end

function PetscDataTypeGetSize_(dtype::PetscDataType)
    datasize = Array(Csize_t, 1)
    ccall((:PetscDataTypeGetSize,petsc),PetscErrorCode,(PetscDataType,Ptr{Csize_t}), dtype, datasize)

    return datasize[1]
end

# define types that depend on the options Petsc was compiled with
(real_type_enum, found_real) = PetscDataTypeFromString_("Real")
(scalar_type_enum, found_scalar) = PetscDataTypeFromString_("Scalar")
int_size = PetscDataTypeGetSize_(PETSC_INT)

# confirm a value was found
@assert(found_real)
@assert(found_scalar)

# figure out what the values mean
if real_type_enum == PETSC_DOUBLE
  precision = "double"
elseif real_type_enum == PETSC_FLOAT
  precision = "single"
else
  println("unknown type of Real")
  println("real_type_enum = ", real_type_enum)
end

if scalar_type_enum == real_type_enum
  scalar_type = "real"
elseif scalar_type_enum == PETSC_COMPLEX
  scalar_type = "complex"
else
  println("unknown type of Scalar")
  println("scalar_type_enum = ", scalar_type_enum)
end

# figure out types

if scalar_type == "real"
  # single or double precision real
  if precision == "single"
    scalar_dtype = Float32
    real_dtype = Float32
  else
    scalar_dtype = Float64
    real_dtype = Float64
  end
else  # scalar_type = complex
  if precision == "single"
    scalar_dtype = Complex64
    real_dtype = Float32
  else
    scalar_dtype = Complex128
    real_dtype = Float64
  end
end

if int_size == 4
   int_dtype = Int32
elseif int_size == 8
   int_dtype = Int64
else
  println("unknown Int size")
  println("int_size = ", int_size)
end

println("PetscScalar type = ", scalar_dtype)
println("PetscReal type = ", real_dtype)
println("PetscInt type = ", int_dtype)

typealias PetscScalar scalar_dtype
typealias PetscReal real_dtype
typealias PetscInt int_dtype

# some useful type unions
PetscInt_arr_or_null = Union{AbstractArray{PetscInt}, Ptr{Void}}

# used to C_NULL in place of an array
const PetscIntNullArray = pointer_to_array(Ptr{PetscInt}(0), 0)

#=
const PETSC_PI = pi
const PETSC_MAX_INT = 2147483647
const PETSC_MIN_INT = -PETSC_MAX_INT - 1
const PETSC_MAX_REAL = PETSC_MAX_INT  # made up
const PETSC_INFINITY = PETSC_MAX_REAL / 4.0
const PETSC_NINFINITY = -PETSC_INFINITY
#const PassiveReal = Cint  # what is this?
=#

# some useful constants
#const PassiveScalar = PetscScalar
#const MPIU_MATSCALAR = MPIU_SCALAR
#const MPIU_2INT = Cint
global const PETSC_NULL = C_NULL
global const PETSC_IGNORE = C_NULL
global const PETSC_DECIDE = convert(Int32, -1)
global const PETSC_DETERMINE = PETSC_DECIDE
global const PETSC_DEFAULT = convert(Int32, -2)
global const PETSC_COMM_SELF = MPI.COMM_SELF

typealias MPI_Comm MPI.Comm  # use MPI package communicator type
# typealias to size of MPI communicator value
# this is not defined by the MPI C standard, so it might be
# 32 bits (MPICH), or possibly 64 bits
typealias comm_type MPI.CComm

global const PETSC_INSERT_VALUES = convert(Int32, 1);
global const PETSC_ADD_VALUES    = convert(Int32, 2);
global const PETSC_COPY_VALUES   = convert(Int32, 0);

typealias NormType Int32
global const NORM_1         = convert(Int32, 0);
global const NORM_2         = convert(Int32, 1);
global const NORM_FROBENIUS = convert(Int32 ,2);
global const NORM_INFINITY  = convert(Int32 ,3);
global const NORM_MAX       = NORM_INFINITY;

global const  PETSC_MAT_FLUSH_ASSEMBLY = convert(Int32, 1)
global const  PETSC_MAT_FINAL_ASSEMBLY = convert(Int32, 0)

# vector formats
global const VECSEQ = "seq"
global const VECMPI = "mpi"
global const VECSTANDARD = "standard"
global const VECSHARED = "shared"
global const VECSEQCUSP = "seqcusp"
global const VECMPICUSP = "mpicusp"
global const VECCUSP = "cusp"
global const VECSEQVIENNACL = "seqviennacl"
global const VECMPIVIENNACL = "mpiviennacl"
global const VECVIENNACL = "viennacl"
global const VECNEST = "nest"
global const VECSEQPTHREAD = "seqpthread"
global const VECMPIPTHREAD = "mpipthread"
global const VECPTHREAD = "pthread"

export PetscMatStructure
export DIFFERENT_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, SAME_NONZERO_PATTERN
# begin enum MatStructure
typealias PetscMatStructure Int32
global const DIFFERENT_NONZERO_PATTERN = (Int32)(0)
global const SUBSET_NONZERO_PATTERN = (Int32)(1)
global const SAME_NONZERO_PATTERN = (Int32)(2)
# end enum MatStructure

export MATSEQAIJ, MATMPIAIJ, MATMPIBAIJ
# matrix types
global const MATSAME = "same"
global const MATMAIJ = "maij"
global const MATSEQMAIJ = "seqmaij"
global const MATMPIMAIJ = "mpimaij"
global const MATIS = "is"
global const MATAIJ = "aij"
global const MATSEQAIJ = "seqaij"
global const MATSEQAIJPTHREAD = "seqaijpthread"
global const MATAIJPTHREAD = "aijpthread"
global const MATMPIAIJ = "mpiaij"
global const MATAIJCRL = "aijcrl"
global const MATSEQAIJCRL = "seqaijcrl"
global const MATMPIAIJCRL = "mpiaijcrl"
global const MATAIJCUSP = "aijcusp"
global const MATSEQAIJCUSP = "seqaijcusp"
global const MATMPIAIJCUSP = "mpiaijcusp"
global const MATAIJCUSPARSE = "aijcusparse"
global const MATSEQAIJCUSPARSE = "seqaijcusparse"
global const MATMPIAIJCUSPARSE = "mpiaijcusparse"
global const MATAIJVIENNACL = "aijviennacl"
global const MATSEQAIJVIENNACL = "seqaijviennacl"
global const MATMPIAIJVIENNACL = "mpiaijviennacl"
global const MATAIJPERM = "aijperm"
global const MATSEQAIJPERM = "seqaijperm"
global const MATMPIAIJPERM = "mpiaijperm"
global const MATSHELL = "shell"
global const MATDENSE = "dense"
global const MATSEQDENSE = "seqdense"
global const MATMPIDENSE = "mpidense"
global const MATELEMENTAL = "elemental"
global const MATBAIJ = "baij"
global const MATSEQBAIJ = "seqbaij"
global const MATMPIBAIJ = "mpibaij"
global const MATMPIADJ = "mpiadj"
global const MATSBAIJ = "sbaij"
global const MATSEQSBAIJ = "seqsbaij"
global const MATMPISBAIJ = "mpisbaij"
global const MATSEQBSTRM = "seqbstrm"
global const MATMPIBSTRM = "mpibstrm"
global const MATBSTRM = "bstrm"
global const MATSEQSBSTRM = "seqsbstrm"
global const MATMPISBSTRM = "mpisbstrm"
global const MATSBSTRM = "sbstrm"
global const MATDAAD = "daad"
global const MATMFFD = "mffd"
global const MATNORMAL = "normal"
global const MATLRC = "lrc"
global const MATSCATTER = "scatter"
global const MATBLOCKMAT = "blockmat"
global const MATCOMPOSITE = "composite"
global const MATFFT = "fft"
global const MATFFTW = "fftw"
global const MATSEQCUFFT = "seqcufft"
global const MATTRANSPOSEMAT = "transpose"
global const MATSCHURCOMPLEMENT = "schurcomplement"
global const MATPYTHON = "python"
global const MATHYPRESTRUCT = "hyprestruct"
global const MATHYPRESSTRUCT = "hypresstruct"
global const MATSUBMATRIX = "submatrix"
global const MATLOCALREF = "localref"
global const MATNEST = "nest"

export MAT_LOCAL, MAT_GLOBAL_MAX, MAT_GLOBAL_SUM
typealias MatInfoType Int32
const MAT_LOCAL = (Int32)(1)
const MAT_GLOBAL_MAX = (Int32)(2)
const MAT_GLOBAL_SUM = (Int32)(3)
#

export MAT_INITIAL_MATRIX, MAT_REUSE_MATRIX, MAT_IGNORE_MATRIX
typealias MatReuse Int32
const MAT_INITIAL_MATRIX = (Int32)(0)
const MAT_REUSE_MATRIX = (Int32)(1)
const MAT_IGNORE_MATRIX = (Int32)(2)
#

export KSPType
# types of KSP solvers
typealias KSPType ASCIIString
global const KSPRICHARDSON = "richardson"
global const KSPCHEBYSHEV = "chebyshev"
global const KSPCG       =  "cg"
global const KSPGROPPCG  =  "groppcg"
global const KSPPIPECG   =  "pipecg"
global const   KSPCGNE   =    "cgne"
global const   KSPNASH   =    "nash"
global const   KSPSTCG   =    "stcg"
global const   KSPGLTR   =    "gltr"
global const KSPFCG      =  "fcg"
global const KSPGMRES    =  "gmres"
global const   KSPFGMRES  =   "fgmres"
global const   KSPLGMRES  =   "lgmres"
global const   KSPDGMRES  =   "dgmres"
global const   KSPPGMRES  =   "pgmres"
global const KSPTCQMR     = "tcqmr"
global const KSPBCGS      = "bcgs"
global const   KSPIBCGS   =   "ibcgs"
global const   KSPFBCGS   =   "fbcgs"
global const   KSPFBCGSR  =   "fbcgsr"
global const   KSPBCGSL   =   "bcgsl"
global const KSPCGS       = "cgs"
global const KSPTFQMR     = "tfqmr"
global const KSPCR        = "cr"
global const KSPPIPECR    = "pipecr"
global const KSPLSQR      = "lsqr"
global const KSPPREONLY   = "preonly"
global const KSPQCG       = "qcg"
global const KSPBICG      = "bicg"
global const KSPMINRES    = "minres"
global const KSPSYMMLQ    = "symmlq"
global const KSPLCD       = "lcd"
global const KSPPYTHON    = "python"
global const KSPGCR       = "gcr"

typealias KSPNormType Cint
global const KSP_NORM_DEFAULT = (Int32)(-1)
global const KSP_NORM_NONE = (Int32)(0)
global const KSP_NORM_PRECONDITIONED = (Int32)(1)
global const KSP_NORM_UNPRECONDITIONED = (Int32)(2)
global const KSP_NORM_NATURAL = (Int32)(3)
# end enum KSPNormType

export KSPConvergedReason

typealias KSPConvergedReason Cint
global const KSP_CONVERGED_RTOL_NORMAL = (Int32)(1)
global const KSP_CONVERGED_ATOL_NORMAL = (Int32)(9)
global const KSP_CONVERGED_RTOL = (Int32)(2)
global const KSP_CONVERGED_ATOL = (Int32)(3)
global const KSP_CONVERGED_ITS = (Int32)(4)
global const KSP_CONVERGED_CG_NEG_CURVE = (Int32)(5)
global const KSP_CONVERGED_CG_CONSTRAINED = (Int32)(6)
global const KSP_CONVERGED_STEP_LENGTH = (Int32)(7)
global const KSP_CONVERGED_HAPPY_BREAKDOWN = (Int32)(8)
global const KSP_DIVERGED_NULL = (Int32)(-2)
global const KSP_DIVERGED_ITS = (Int32)(-3)
global const KSP_DIVERGED_DTOL = (Int32)(-4)
global const KSP_DIVERGED_BREAKDOWN = (Int32)(-5)
global const KSP_DIVERGED_BREAKDOWN_BICG = (Int32)(-6)
global const KSP_DIVERGED_NONSYMMETRIC = (Int32)(-7)
global const KSP_DIVERGED_INDEFINITE_PC = (Int32)(-8)
global const KSP_DIVERGED_NANORINF = (Int32)(-9)
global const KSP_DIVERGED_INDEFINITE_MAT = (Int32)(-10)
global const KSP_DIVERGED_PCSETUP_FAILED = (Int32)(-11)
global const KSP_CONVERGED_ITERATING = (Int32)(0)
#
# map values to string for printing
export KSPConvergedReasonDict
global const KSPConvergedReasonDict = Dict{KSPConvergedReason, ASCIIString}(
 KSP_CONVERGED_RTOL_NORMAL => "Converged: RTol normal",
 KSP_CONVERGED_ATOL_NORMAL => "Converged: ATol normal",
 KSP_CONVERGED_RTOL => "Converged: RTol",
 KSP_CONVERGED_ATOL => "converged: ATol",
 KSP_CONVERGED_ITS => "Converged: ITS",
 KSP_CONVERGED_CG_NEG_CURVE => "Converged: CG Negative Curvature",
 KSP_CONVERGED_CG_CONSTRAINED => "Converged: CG Constrained",
 KSP_CONVERGED_STEP_LENGTH => "Converged: step length",
 KSP_CONVERGED_HAPPY_BREAKDOWN => "Converged: happy breakdown",
 KSP_DIVERGED_NULL => "Diverged: Null",
 KSP_DIVERGED_ITS => "Diverged: iteration max",
 KSP_DIVERGED_DTOL => "Diverged divergence tolerance",
 KSP_DIVERGED_BREAKDOWN => "Diverged: Krylov breakdown",
 KSP_DIVERGED_BREAKDOWN_BICG => "Diverged: BICG breakdown",
 KSP_DIVERGED_NONSYMMETRIC => "Diverged: non symmetric",
 KSP_DIVERGED_INDEFINITE_PC => "Diverged: indefinitate PC",
 KSP_DIVERGED_NANORINF => "Diverged nan or inf",
 KSP_DIVERGED_INDEFINITE_MAT => "Diverged: indefinate matrix",
 KSP_DIVERGED_PCSETUP_FAILED => "Diverged: PC setup failed",
 KSP_CONVERGED_ITERATING => "Still running: please be patient",
 )
#

global const KSP_NORM_MAX = KSP_NORM_NATURAL + 1

typealias PCType ASCIIString

global const PCNONE = "none"
global const PCJACOBI = "jacobi"
global const PCSOR = "sor"
global const PCLU = "lu"
global const PCSHELL = "shell"
global const PCBJACOBI = "bjacobi"
global const PCMG = "mg"
global const PCEISENSTAT = "eisenstat"
global const PCILU = "ilu"
global const PCICC = "icc"
global const PCASM = "asm"
global const PCGASM = "gasm"
global const PCKSP = "ksp"
global const PCCOMPOSITE = "composite"
global const PCREDUNDANT = "redundant"
global const PCSPAI = "spai"
global const PCNN = "nn"
global const PCCHOLESKY = "cholesky"
global const PCPBJACOBI = "pbjacobi"
global const PCMAT = "mat"
global const PCHYPRE = "hypre"
global const PCPARMS = "parms"
global const PCFIELDSPLIT = "fieldsplit"
global const PCTFS = "tfs"
global const PCML = "ml"
global const PCGALERKIN = "galerkin"
global const PCEXOTIC = "exotic"
global const PCCP = "cp"
global const PCBFBT = "bfbt"
global const PCLSC = "lsc"
global const PCPYTHON = "python"
global const PCPFMG = "pfmg"
global const PCSYSPFMG = "syspfmg"
global const PCREDISTRIBUTE = "redistribute"
global const PCSVD = "svd"
global const PCGAMG = "gamg"
global const PCSACUSP = "sacusp"
global const PCSACUSPPOLY = "sacusppoly"
global const PCBICGSTABCUSP = "bicgstabcusp"
global const PCAINVCUSP = "ainvcusp"
global const PCBDDC = "bddc"
global const PCKACZMARZ = "kaczmarz"

# begin enum PCSide
typealias PCSide Cint
global const PC_SIDE_DEFAULT = (Int32)(-1)
global const PC_LEFT = (Int32)(0)
global const PC_RIGHT = (Int32)(1)
global const PC_SYMMETRIC = (Int32)(2)
# end enum PCSide

typealias PCJacobiType Int32
global const PC_JACOBI_DIAGONAL = (Int32)(0)
global const PC_JACOBI_ROWMAX = (Int32)(1)
global const PC_JACOBI_ROWSUM = (Int32)(2)
#

typealias PetscErrorType Int32
global const PETSC_ERROR_INITIAL = (Int32)(0)
global const PETSC_ERROR_REPEAT = (Int32)(1)
global const PETSC_ERROR_IN_CXX = (Int32)(2)
#

# begin enum MatOption
typealias MatOption Cint
const MAT_OPTION_MIN = (Int32)(-5)
const MAT_NEW_NONZERO_LOCATION_ERR = (Int32)(-4)
const MAT_UNUSED_NONZERO_LOCATION_ERR = (Int32)(-3)
const MAT_NEW_NONZERO_ALLOCATION_ERR = (Int32)(-2)
const MAT_ROW_ORIENTED = (Int32)(-1)
const MAT_SYMMETRIC = (Int32)(1)
const MAT_STRUCTURALLY_SYMMETRIC = (Int32)(2)
const MAT_NEW_DIAGONALS = (Int32)(3)
const MAT_IGNORE_OFF_PROC_ENTRIES = (Int32)(4)
const MAT_USE_HASH_TABLE = (Int32)(5)
const MAT_KEEP_NONZERO_PATTERN = (Int32)(6)
const MAT_IGNORE_ZERO_ENTRIES = (Int32)(7)
const MAT_USE_INODES = (Int32)(8)
const MAT_HERMITIAN = (Int32)(9)
const MAT_SYMMETRY_ETERNAL = (Int32)(10)
const MAT_DUMMY = (Int32)(11)
const MAT_IGNORE_LOWER_TRIANGULAR = (Int32)(12)
const MAT_ERROR_LOWER_TRIANGULAR = (Int32)(13)
const MAT_GETROW_UPPERTRIANGULAR = (Int32)(14)
const MAT_SPD = (Int32)(15)
const MAT_NO_OFF_PROC_ZERO_ROWS = (Int32)(16)
const MAT_NO_OFF_PROC_ENTRIES = (Int32)(17)
const MAT_NEW_NONZERO_LOCATIONS = (Int32)(18)
const MAT_OPTION_MAX = (Int32)(19)
# end enum MatOption

# begin enum MatOperation
typealias MatOperation UInt32
const MATOP_SET_VALUES = (UInt32)(0)
const MATOP_GET_ROW = (UInt32)(1)
const MATOP_RESTORE_ROW = (UInt32)(2)
const MATOP_MULT = (UInt32)(3)
const MATOP_MULT_ADD = (UInt32)(4)
const MATOP_MULT_TRANSPOSE = (UInt32)(5)
const MATOP_MULT_TRANSPOSE_ADD = (UInt32)(6)
const MATOP_SOLVE = (UInt32)(7)
const MATOP_SOLVE_ADD = (UInt32)(8)
const MATOP_SOLVE_TRANSPOSE = (UInt32)(9)
const MATOP_SOLVE_TRANSPOSE_ADD = (UInt32)(10)
const MATOP_LUFACTOR = (UInt32)(11)
const MATOP_CHOLESKYFACTOR = (UInt32)(12)
const MATOP_SOR = (UInt32)(13)
const MATOP_TRANSPOSE = (UInt32)(14)
const MATOP_GETINFO = (UInt32)(15)
const MATOP_EQUAL = (UInt32)(16)
const MATOP_GET_DIAGONAL = (UInt32)(17)
const MATOP_DIAGONAL_SCALE = (UInt32)(18)
const MATOP_NORM = (UInt32)(19)
const MATOP_ASSEMBLY_BEGIN = (UInt32)(20)
const MATOP_ASSEMBLY_END = (UInt32)(21)
const MATOP_SET_OPTION = (UInt32)(22)
const MATOP_ZERO_ENTRIES = (UInt32)(23)
const MATOP_ZERO_ROWS = (UInt32)(24)
const MATOP_LUFACTOR_SYMBOLIC = (UInt32)(25)
const MATOP_LUFACTOR_NUMERIC = (UInt32)(26)
const MATOP_CHOLESKY_FACTOR_SYMBOLIC = (UInt32)(27)
const MATOP_CHOLESKY_FACTOR_NUMERIC = (UInt32)(28)
const MATOP_SETUP_PREALLOCATION = (UInt32)(29)
const MATOP_ILUFACTOR_SYMBOLIC = (UInt32)(30)
const MATOP_ICCFACTOR_SYMBOLIC = (UInt32)(31)
const MATOP_DUPLICATE = (UInt32)(34)
const MATOP_FORWARD_SOLVE = (UInt32)(35)
const MATOP_BACKWARD_SOLVE = (UInt32)(36)
const MATOP_ILUFACTOR = (UInt32)(37)
const MATOP_ICCFACTOR = (UInt32)(38)
const MATOP_AXPY = (UInt32)(39)
const MATOP_GET_SUBMATRICES = (UInt32)(40)
const MATOP_INCREASE_OVERLAP = (UInt32)(41)
const MATOP_GET_VALUES = (UInt32)(42)
const MATOP_COPY = (UInt32)(43)
const MATOP_GET_ROW_MAX = (UInt32)(44)
const MATOP_SCALE = (UInt32)(45)
const MATOP_SHIFT = (UInt32)(46)
const MATOP_DIAGONAL_SET = (UInt32)(47)
const MATOP_ZERO_ROWS_COLUMNS = (UInt32)(48)
const MATOP_SET_RANDOM = (UInt32)(49)
const MATOP_GET_ROW_IJ = (UInt32)(50)
const MATOP_RESTORE_ROW_IJ = (UInt32)(51)
const MATOP_GET_COLUMN_IJ = (UInt32)(52)
const MATOP_RESTORE_COLUMN_IJ = (UInt32)(53)
const MATOP_FDCOLORING_CREATE = (UInt32)(54)
const MATOP_COLORING_PATCH = (UInt32)(55)
const MATOP_SET_UNFACTORED = (UInt32)(56)
const MATOP_PERMUTE = (UInt32)(57)
const MATOP_SET_VALUES_BLOCKED = (UInt32)(58)
const MATOP_GET_SUBMATRIX = (UInt32)(59)
const MATOP_DESTROY = (UInt32)(60)
const MATOP_VIEW = (UInt32)(61)
const MATOP_CONVERT_FROM = (UInt32)(62)
const MATOP_MATMAT_MULT = (UInt32)(63)
const MATOP_MATMAT_MULT_SYMBOLIC = (UInt32)(64)
const MATOP_MATMAT_MULT_NUMERIC = (UInt32)(65)
const MATOP_SET_LOCAL_TO_GLOBAL_MAP = (UInt32)(66)
const MATOP_SET_VALUES_LOCAL = (UInt32)(67)
const MATOP_ZERO_ROWS_LOCAL = (UInt32)(68)
const MATOP_GET_ROW_MAX_ABS = (UInt32)(69)
const MATOP_GET_ROW_MIN_ABS = (UInt32)(70)
const MATOP_CONVERT = (UInt32)(71)
const MATOP_SET_COLORING = (UInt32)(72)
const MATOP_SET_VALUES_ADIFOR = (UInt32)(74)
const MATOP_FD_COLORING_APPLY = (UInt32)(75)
const MATOP_SET_FROM_OPTIONS = (UInt32)(76)
const MATOP_MULT_CONSTRAINED = (UInt32)(77)
const MATOP_MULT_TRANSPOSE_CONSTRAIN = (UInt32)(78)
const MATOP_FIND_ZERO_DIAGONALS = (UInt32)(79)
const MATOP_MULT_MULTIPLE = (UInt32)(80)
const MATOP_SOLVE_MULTIPLE = (UInt32)(81)
const MATOP_GET_INERTIA = (UInt32)(82)
const MATOP_LOAD = (UInt32)(83)
const MATOP_IS_SYMMETRIC = (UInt32)(84)
const MATOP_IS_HERMITIAN = (UInt32)(85)
const MATOP_IS_STRUCTURALLY_SYMMETRIC = (UInt32)(86)
const MATOP_SET_VALUES_BLOCKEDLOCAL = (UInt32)(87)
const MATOP_GET_VECS = (UInt32)(88)
const MATOP_MAT_MULT = (UInt32)(89)
const MATOP_MAT_MULT_SYMBOLIC = (UInt32)(90)
const MATOP_MAT_MULT_NUMERIC = (UInt32)(91)
const MATOP_PTAP = (UInt32)(92)
const MATOP_PTAP_SYMBOLIC = (UInt32)(93)
const MATOP_PTAP_NUMERIC = (UInt32)(94)
const MATOP_MAT_TRANSPOSE_MULT = (UInt32)(95)
const MATOP_MAT_TRANSPOSE_MULT_SYMBO = (UInt32)(96)
const MATOP_MAT_TRANSPOSE_MULT_NUMER = (UInt32)(97)
const MATOP_CONJUGATE = (UInt32)(102)
const MATOP_SET_VALUES_ROW = (UInt32)(104)
const MATOP_REAL_PART = (UInt32)(105)
const MATOP_IMAGINARY_PART = (UInt32)(106)
const MATOP_GET_ROW_UPPER_TRIANGULAR = (UInt32)(107)
const MATOP_RESTORE_ROW_UPPER_TRIANG = (UInt32)(108)
const MATOP_MAT_SOLVE = (UInt32)(109)
const MATOP_GET_REDUNDANT_MATRIX = (UInt32)(110)
const MATOP_GET_ROW_MIN = (UInt32)(111)
const MATOP_GET_COLUMN_VECTOR = (UInt32)(112)
const MATOP_MISSING_DIAGONAL = (UInt32)(113)
const MATOP_GET_SEQ_NONZERO_STRUCTUR = (UInt32)(114)
const MATOP_CREATE = (UInt32)(115)
const MATOP_GET_GHOSTS = (UInt32)(116)
const MATOP_GET_LOCAL_SUB_MATRIX = (UInt32)(117)
const MATOP_RESTORE_LOCALSUB_MATRIX = (UInt32)(118)
const MATOP_MULT_DIAGONAL_BLOCK = (UInt32)(119)
const MATOP_HERMITIAN_TRANSPOSE = (UInt32)(120)
const MATOP_MULT_HERMITIAN_TRANSPOSE = (UInt32)(121)
const MATOP_MULT_HERMITIAN_TRANS_ADD = (UInt32)(122)
const MATOP_GET_MULTI_PROC_BLOCK = (UInt32)(123)
const MATOP_FIND_NONZERO_ROWS = (UInt32)(124)
const MATOP_GET_COLUMN_NORMS = (UInt32)(125)
const MATOP_INVERT_BLOCK_DIAGONAL = (UInt32)(126)
const MATOP_GET_SUB_MATRICES_PARALLE = (UInt32)(128)
const MATOP_SET_VALUES_BATCH = (UInt32)(129)
const MATOP_TRANSPOSE_MAT_MULT = (UInt32)(130)
const MATOP_TRANSPOSE_MAT_MULT_SYMBO = (UInt32)(131)
const MATOP_TRANSPOSE_MAT_MULT_NUMER = (UInt32)(132)
const MATOP_TRANSPOSE_COLORING_CREAT = (UInt32)(133)
const MATOP_TRANS_COLORING_APPLY_SPT = (UInt32)(134)
const MATOP_TRANS_COLORING_APPLY_DEN = (UInt32)(135)
const MATOP_RART = (UInt32)(136)
const MATOP_RART_SYMBOLIC = (UInt32)(137)
const MATOP_RART_NUMERIC = (UInt32)(138)
const MATOP_SET_BLOCK_SIZES = (UInt32)(139)
const MATOP_AYPX = (UInt32)(140)
const MATOP_RESIDUAL = (UInt32)(141)
const MATOP_FDCOLORING_SETUP = (UInt32)(142)
const MATOP_MPICONCATENATESEQ = (UInt32)(144)
# end enum MatOperation


text:
# save: frequencies, positions, coverage, qualityThreshold, removeambiguous, algoritm,
#       method, newton regularization (if applicable)

function savedictJLD(fileprefix, dict)
	save("$fileprefix.jld", dict, compress=true)
end

function savedict(fileprefix, dict, outFormat=:JLD)
	if outFormat==:JLD
		savedictJLD(fileprefix, dict)
	else
		error("Unknown format: $outFormat")
	end
end

function computecodonfrequencies(samplePath::String, sampleName::String,
                                 outFolder::String;
                                 strands=:both, mappingQualityThreshold=30,
                                 baseQualityThreshold=30, removeAmbiguous=true,
                                 method=:Newton, newtonRegularization=1e-6,
                                 maxIter=10000,
                                 outFormat=:JLD)
	log = IOBuffer()

	println(log, "Computing codon frequencies: ", sampleName)
	startTime = time()

	bamFile = BamFile(samplePath)
	freqs,positions,coverage = mlcodonfreqs(bamFile,strands=strands,
	                                        mappingQualityThreshold=mappingQualityThreshold,
	                                        baseQualityThreshold=baseQualityThreshold,
	                                        log=log,removeAmbiguous=removeAmbiguous,method=method,
	                                        newtonRegularization=newtonRegularization,
	                                        maxIter=maxIter)

	fileprefix = joinpath(outFolder,sampleName)
	d = Dict{String,Any}("codonFreqs"=>freqs,
	                          "positions"=>positions,
	                          "coverage"=>coverage,
	                          "segmentInfo"=>sequences(bamFile),
	                          "strands"=>string(strands),
	                          "mappingQualityThreshold" => mappingQualityThreshold,
	                          "baseQualityThreshold" => baseQualityThreshold,
	                          "removeAmbiguousBases" => removeAmbiguous,
	                          "algorithm" => "Maximum Likelihood",
	                          "solver" => string(method))
	method==:Newton && (d["newtonRegularization"] = newtonRegularization)

	typeof(outFormat) <: AbstractArray || (outFormat = [outFormat])
	for of in outFormat
		savedict(fileprefix,d,of)
	end

	duration = time()-startTime
	println(log, "Done in ", duration, "s.")
	String(take!(log))
end

function computecodonfrequencies(syn, samplePaths::Vector{String},
                                 sampleNames::Vector{String},
                                 outFolder::String;
                                 log=STDOUT, bamDir="bam",
                                 kwargs...)

	# Threading modelled after the pmap implementation in http://docs.julialang.org/en/release-0.4/manual/parallel-computing/ (NB: not the same as actual pmap())
	# --------------------------------------------------------------------------
	procList = procs() # find the available processes
	n = length(samplePaths)
	i = 1
	# function to produce the next work item from the queue.
	# in this case it's just an index.
	nextidx() = (idx=i; i+=1; idx)
	@sync begin
		for p in procList
			if p != myid() || length(procList)==1
				@async begin
					while true
						idx = nextidx()
						if idx > n
							break
						end

						sPath = samplePaths[idx]
						sName = sampleNames[idx]

						println("Downloading sample $sName")
						sLocal = localpath(syn,sPath,downloadLocation=bamDir,ifcollision="overwrite.local")

						println("Computing codon frequencies for sample $sName")

						# compute codon frequencies in in worker thread
						logStr = fetch(@spawnat p computecodonfrequencies(sLocal,sName,outFolder;kwargs...))
						print(log, logStr); flush(log)
						println("Finished computing codon frequencies for sample $sName")
					end
				end
			end
		end
	end
	# --------------------------------------------------------------------------

	nothing
end


text:
using DataFrames
using CSV
using Lathe.preprocess:TrainTestSplit
using GLM
using StatsBase
using StatsModels
using ROC
using Plots

data = DataFrame(CSV.File("pacijenti1000a.csv"))

dataTrain, dataTest = TrainTestSplit(data, .80)

fm = @formula(bolest ~ visina + tezina + dbp + sbp)

logisticRegressor = glm(fm, dataTrain, Binomial(), ProbitLink())

dataPredictTest = predict(logisticRegressor, dataTest)

dataPredictTestClass = zeros(length(dataPredictTest))

for i in 1:length(dataPredictBoja)
    if dataPredictTest[i] < 0.5
        dataPredictTestClass[i] = 0
    else
        dataPredictTestClass[i] = 1
    end
end

brojObolelihOdDijabetesaPredict = sum(dataPredictTestClass)
brojObolelihOdDijabetesa = sum(dataTest.bolest)

brojZdravihPredict = length(dataPredictTestClass) - brojObolelihOdDijabetesaPredict
brojZdravih = length(data.bolest) - brojObolelihOdDijabetesa

println("Predvidjen broj obolelih: $brojObolelihOdDijabetesaPredict")
println("Broj osoba koje imaju dijabetes: $brojObolelihOdDijabetesa")

println("Predvidjen broj zdravih: $brojZdravihPredict")
println("Broj zdravih osoba: $brojZdravih")

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
struct Transition
    utc_datetime::DateTime  # Instant where new zone applies
    zone::FixedTimeZone
end

Base.isless(x::Transition,y::Transition) = isless(x.utc_datetime,y.utc_datetime)

"""
    VariableTimeZone

A `TimeZone` with an offset that changes over time.
"""
struct VariableTimeZone <: TimeZone
    name::String
    transitions::Vector{Transition}
    cutoff::Union{DateTime,Nothing}

    function VariableTimeZone(name::AbstractString, transitions::Vector{Transition}, cutoff::Union{DateTime,Nothing}=nothing)
        new(name, transitions, cutoff)
    end
end

name(tz::VariableTimeZone) = tz.name

function rename(tz::VariableTimeZone, name::AbstractString)
    VariableTimeZone(name, tz.transitions, tz.cutoff)
end

function Base.:(==)(a::VariableTimeZone, b::VariableTimeZone)
    a.name == b.name && a.transitions == b.transitions
end

function Base.isequal(a::VariableTimeZone, b::VariableTimeZone)
    return (
        isequal(a.name, b.name) &&
        isequal(a.transitions, b.transitions) &&
        isequal(a.cutoff, b.cutoff)
    )
end

function Base.hash(tz::VariableTimeZone, h::UInt)
    h = hash(tz.name, h)
    h = hash(tz.transitions, h)
    h = hash(tz.cutoff, h)
    return h
end


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{String, Any}(
"physics" => "Euler",
"run_type" => 5,
"jac_method" => 2,
"jac_type" => 2,
"order" => 1,
"dimensions" => 2,
"Q_transpose" => true, # Weak differentiate option in eval volume integrals
"IC_name" => "ICIsentropicVortex",
#"numBC" => 2,
"numBC" => 1,
#"BC1" => [ 7, 13],
"BC1" => [0],
"BC1_name" => "isentropicVortexBC", # "allOnesBC",
"delta_t" => 0.001, # default 0.005
"t_max" => 500000.000,
"smb_name" => "SRCMESHES/Test1el.smb",
#"dmg_name" => "SRCMESHES/vortex.dmg",
"dmg_name" => ".null",
"res_abstol" => 1e-10,
"res_reltol" => 1e-9,
"step_tol" => 1e-10,
"itermax" => 30,
"writeq" => false,
"step_tol" => 1e-10,
"itermax" => 30,
"use_edgestab" => false,
"edgestab_gamma" => -0.01,
"use_GLS" => true,
"use_filter" => false,
#"use_res_filter" => true,
#"filter_name" => "raisedCosineFilter",
"use_dissipation" => false,
"dissipation_name" => "damp1",
"dissipation_const" => 12.00,
"writeq" => false,
# "perturb_ic" => true,
#"perturb_mag" => 0.001,
#"write_sparsity" => true,
#"write_jac" => true,
"write_edge_vertnums" => false,
"write_face_vertnums" => false,
"write_qic" => false,
"writeboundary" => false,
"write_res" => false,
"print_cond" => false,
#"write_counts" => true,
"write_vis" => true,
"solve" => false,
)


text:
###############################################################################
#
# This script contains Julia code for performing uncertainty analysis of camera
# calbiration and georeferencing, as described in Schweitzer & Cowen, 2021; and
# Schweitzer & Cowen, 2022
#
# # References
# Schweitzer, S. A., & Cowen, E. A. (2021). "Instantaneous river-wide water surface
# velocity field measurements at centimeter scales using infrared quantitative
# image velocimetry". Water Resources Research, 57,
# e2020WR029279. https://doi.org/10.1029/2020WR029279
#
# Schweitzer, S. A., & Cowen, E. A. (2022, exptected) "A Method for Analysis of
# Spatial Uncertainty in Image Based Surface Velocimetry". Frontiers in Water,
# Water and Hydrocomplexity.
#
# # License
# This software is released under the MIT license, however, the authors request to
# be credited if you use the software, preferably by citing the above-mentioned
# Schweitzer & Cowen papers.
#
# MIT License
#
# Copyright (c) 2022 Seth A. Schweitzer
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
###############################################################################

###############################################################################
using Statistics
#
import ExtrinsicCameraCalibration as Calib
import ImageRectification as Rect
import GeoreferencingUncertaintySimulations as GeorefSim
#
dropdim3(x) = dropdims(x;dims=3)
dim3quantiles(data; q=0.95) = mapslices(d->quantile(d, q), data; dims=3) |> d -> dropdims(d; dims=(findall(size(d).==1)[]))
###############################################################################

###############################################################################
#
### Define camera position and orientation :
#
## Camera optical center position:
x_c = 0.0
y_c = 0.0
z_c = 25.0
## Camera orientation:
# Generic setup is 45deg tilt for oblique view, 0.0deg tilt for UAV.
φ = deg2rad(0.0) # φ azimuth. clockwise from the north (y-axis). This is meaningless in a camera-centric coordinate system
τ = deg2rad(45.0) # τ tilt. angle away from the vertical of the camera centeral ray (0 is nadir view).
σ = deg2rad(0.0) # σ roll. rotation of the camera around the centeral ray. (0.0 is no roll, rotation is from the positive x axis)
###############################################################################
# Define actual water surface elevation:
trueWSE = 0.0
# Number of GCPs to use (may later add restrictions on their distribution):
nGCPs = 4 # 6 # 8
GCPzrange = (trueWSE, trueWSE+3.0) # we'll randomly generate a value in this range for each GCP
#
GCPdistribution =  "corners" # "left_and_right" # "corners", "perimeter", "perimeter_ideal", "anywhere", "left_and_right", "top_and_bottom"
###############################################################################
SimulationParametersInfo = join(["φ=$(round(rad2deg(φ),digits=1))°",
                                 "τ=$(round(rad2deg(τ),digits=1))°",
                                 "σ=$(round(rad2deg(σ),digits=1))°",
                                 "z_c=$(z_c)",
                                 "$(nGCPs) GCPs",
                                 GCPdistribution], ", ")
###############################################################################

###############################################################################
#
### GCP perturbation parameters:
#
# "NoiseFactors" represent the standard deviation of the distribution from which
# perturbations will be randomly drawn:
xyzNoiseFactor = [0.03,0.03,0.03]
uvNoiseFactor  = [0.5,0.5]
wseNoiseFactor =  0.03
#
φNoiseFactor = deg2rad(0.8)
τNoiseFactor = deg2rad(0.1)
σNoiseFactor = deg2rad(0.1)
x_cNoiseFactor =  0.02           # horizontal position (RTK)
y_cNoiseFactor =  0.02           # horizontal position (RTK)
z_cNoiseFactor = 0.05           # meters
#
## Do we want to simulate radial distortion?
AddRadialDistortion = true     # false
#
nSimulations = 700
###############################################################################

###############################################################################
#
## Define camera intrinsic and extrinsic parameters, and DLT coefficients:
# Here we use the parameters of the camera used in Schweitzer & Cowen, 2021 WRR.
#
CameraPixelResolution = [1344 784]
LensFocalLength_m = 17e-3
#
ιtrue = Calib.IntrinsicCalibrationParameters(mean([1 CameraPixelResolution[1]]), # u₀
                                             mean([1 CameraPixelResolution[2]]), # v₀
                                             1.0,   # λᵤ
                                             1.0,   # λᵥ
                                             14e-6, # nominal pixel pitch
                                             [1, 3, 5], # radial_distortion_poly_degrees
                                             [0.020522291524825344, -3.156396095419575e-8, -8.859823756718999e-14]) # radial_distortion_poly_coeffs
#
χtrue = Calib.ExtrinsicCalibrationParameters(
    φ, τ, σ,
    -LensFocalLength_m/ιtrue.PixelPitch,    # f; distance b/w FPA and optical center of camera. Focal length/pixel pitch. Note that this may be a negative number (since the FPA is behind the focal point)
    x_c, y_c, z_c,
)
#
trueDLT = Calib.CalculateDLTcoefficients(χtrue, ιtrue)
###############################################################################

###############################################################################
# Work on only a subset of pixels in the image, to improve performance:
row_skip_px = 6                 # 2, 4, 8, etc.
col_skip_px = 6                 # 2, 4, 8, etc.
###############################################################################

###############################################################################
# Set up the pixel grid over which we will run the calibration sensitivity analysis:
u = [u for v in 1:row_skip_px:CameraPixelResolution[2], u in 1:col_skip_px:CameraPixelResolution[1]];
v = [v for v in 1:row_skip_px:CameraPixelResolution[2], u in 1:col_skip_px:CameraPixelResolution[1]];
###############################################################################

###############################################################################
#
### Generate the simulated GCP:
#
# random uvz coordinates:
trueGCPuv, trueGCPz = GeorefSim.GenerateGCPuvCoordinatesForSimulation(
    nGCPs,
    u,v,
    GCPzrange;
    EdgeWidth = 50, # how far in from the edges can GCPs be?
    IdealEdge = 10, # distance from edge for fixed, "ideal" GCPs
    GCPdistribution = GCPdistribution) # how to distribute the GCPs throughout the frame. defined above
# -----------------------------------------------------------------------------
# Calculate GCP x,y coordinates from the uvz coordinates:
trueGCPx, trueGCPy = Rect.FindXYFromUVZ(trueGCPuv[:,1], trueGCPuv[:,2], trueGCPz, trueDLT)
trueGCPxyz = hcat(trueGCPx, trueGCPy, trueGCPz)
###############################################################################

###############################################################################
#
# Calculate the "reference grid" of true pixel footprints:
#
χcalculated = Calib.HollandEtAlExtrinsicCalibration(trueGCPuv, trueGCPxyz, χtrue, ιtrue; MaxIterations = 50)
DLTcalculated = Calib.CalculateDLTcoefficients(χcalculated, ιtrue)
XCalculated, YCalculated = Rect.FindXYFromUVZ(u,v, trueWSE, DLTcalculated)
###############################################################################

###############################################################################
#
# We can now define the complete set of simulation input parameters:
#
GCPSimulationInputs = (; trueGCPxyz, trueGCPuv, trueWSE, u, v, χcalculated, ιtrue)
ExtrinsicParamsSimulationInputs = (; χtrue, ιtrue, trueWSE, u, v)
#
# And the set of uncertainty parameters to use for the simulations:
GCPSimulationUncertaintyParams = (; xyzNoiseFactor,
                                  uvNoiseFactor,
                                  wseNoiseFactor,
                                  AddRadialDistortion,
                                  nSimulations, # how many simulation scenarios to run
                                  MaxIterations = 50, # max iterations of calibration to run before giving up
                                  )
#
ExtrinsicParamsSimulationUncertaintyParams = (; φNoiseFactor, τNoiseFactor, σNoiseFactor,
                                              x_cNoiseFactor,
                                              y_cNoiseFactor,
                                              z_cNoiseFactor,
                                              wseNoiseFactor, #
                                              AddRadialDistortion,
                                              nSimulations, # how many simulation scenarios to run
                                              MaxIterations=50) # max iterations of calibration to run before giving up
###############################################################################

###############################################################################
#
# Apply perturbations, and calculate results!
#
###############################################################################

###############################################################################
# Perturb xyz and uv coordinates of GCPs, and WSE:
ObliqueView_Outcomes_xyz_uv_wse = GeorefSim.GeorefGCPSimulationWrapper(
    GCPSimulationInputs...;
    GCPSimulationUncertaintyParams...,
    nSimulations = nSimulations*2) #
###############################################################################

###############################################################################
# Perturb only xyz coordinates of GCPs:
ObliqueView_Outcomes_xyz = GeorefSim.GeorefGCPSimulationWrapper(
    GCPSimulationInputs...;
    GCPSimulationUncertaintyParams...,
    uvNoiseFactor = 0.0,
    wseNoiseFactor = 0.0)
###############################################################################

###############################################################################
# Perturb only uv coordinates of GCPs:
ObliqueView_Outcomes_uv = GeorefSim.GeorefGCPSimulationWrapper(
    GCPSimulationInputs...;
    GCPSimulationUncertaintyParams...,
    xyzNoiseFactor = 0.0,
    wseNoiseFactor = 0.0)
###############################################################################

###############################################################################
# Perturb only WSE:
ObliqueView_Outcomes_wse = GeorefSim.GeorefGCPSimulationWrapper(
    GCPSimulationInputs...;
    GCPSimulationUncertaintyParams...,
    xyzNoiseFactor = 0.0,
    uvNoiseFactor = 0.0)
###############################################################################

###############################################################################
# Perturb Direct Georeferencing parameters:
ObliqueView_Outcomes_Extrinsic = GeorefSim.GeorefExtrinsicParamsSimulationWrapper(
    ExtrinsicParamsSimulationInputs...;
    ExtrinsicParamsSimulationUncertaintyParams...,
    nSimulations = nSimulations*2,
)
###############################################################################

###############################################################################
###############################################################################
#
# Gradients of uncertainty from here!!
#
###############################################################################
###############################################################################

###############################################################################
""" Calculate gradients in 1 or 2 dimensions, using a central differencing scheme. """
function central_diff(data_array::AbstractMatrix;
                      dims=nothing,
                      grid_spacing_1=1.0,
                      grid_spacing_2=1.0)
    if any(dims .== 1)
        d1    = diff(data_array, dims=1)/2
        _g1   = cat(d1[[1],:], d1; dims=1) # put the first row back in to recover the original shape
        _g1 .+= cat(d1, d1[[end],:]; dims=1)
        _g1 ./= grid_spacing_1
    end
    if any(dims .== 2)
        d2    = diff(permutedims(data_array), dims=1)/2
        _g2   = cat(d2[[1],:], d2; dims=1) # put the first row back in to recover the original shape
        _g2 .+= cat(d2, d2[[end],:]; dims=1)
        _g2 = permutedims(_g2)
        _g2 ./= grid_spacing_2
    end
    if dims == 1
        return _g1
    elseif dims == 2
        return _g2
    else
        return _g1 .+ _g2 .* im # return complex number, where _g1 is real part, _g2 is imaginary
    end
end
###############################################################################

###############################################################################
"""
Calculate the 2D gradient of difference between reference and simulated
coordinates (i.e., the georeferencing error), at each pixel footprint
location. Results are complex numbers at each node, in which the real part
represents dim1 (y direction) and the imaginary part represents dim2 (x
direction).
"""
function Calculate2DCoordinateGradient(SimulationOutcomes)
    # Calculate the (irregular) "grid" spacing along each coordinate:
    grid_spac1 = central_diff(SimulationOutcomes.Calculated_y[:,:]; dims=1) # distance between y coordinate of grid-adjacent points (y is dim 1)
    grid_spac2 = central_diff(SimulationOutcomes.Calculated_x[:,:]; dims=2) # distance between x coordinate of grid-adjacent points (x is dim 2)
    #
    # Calculate the gradients of uncertainty along each dimension, for each simulation:
    # y coordinate (varies along dim 1!)
    _y_coordinate_gradient = mapslices(y->central_diff(y; dims=(1,2),
                                                          grid_spacing_1=grid_spac1,
                                                          grid_spacing_2=grid_spac2),
                                       SimulationOutcomes.Modified_y .- SimulationOutcomes.Calculated_y;
                                       dims=(1,2))
    # Note: this is 2D, despite only considering the y-direction, because the
    # reference points are not laid out on a regular grid, and hence "grid-adjacent"
    # nodes have different coordinates in two dimensions. We calculate the slope
    # with which the difference between (1D) simulated and reference "y" coordinate
    # changes along each of the two directions (x,y), at each location on the
    # (reference) grid.
    #
    # Similar to above, but for the "x" coordinate (varies along dim 2!):
    _x_coordinate_gradient = mapslices(x->central_diff(x; dims=(1,2),
                                                          grid_spacing_1=grid_spac1,
                                                          grid_spacing_2=grid_spac2),
                                       SimulationOutcomes.Modified_x .- SimulationOutcomes.Calculated_x;
                                       dims=(1,2))
    #
    # return an array of complex numbers, in which the real component is the
    # total gradient in the y direction, and the the imaginary component is the
    # gradient in the x direction:
    return (; _x_coordinate_gradient, _y_coordinate_gradient)
end
###############################################################################

###############################################################################
#
# Calculate the gradients of uncertainty resulting from spatial uncertainty:
#
ObliqueView_UncertGradients_xyz_uv_wse = Calculate2DCoordinateGradient(ObliqueView_Outcomes_xyz_uv_wse)
ObliqueView_UncertGradients_xyz = Calculate2DCoordinateGradient(ObliqueView_Outcomes_xyz)
ObliqueView_UncertGradients_uv = Calculate2DCoordinateGradient(ObliqueView_Outcomes_uv)
ObliqueView_UncertGradients_wse = Calculate2DCoordinateGradient(ObliqueView_Outcomes_wse)
ObliqueView_UncertGradients_Extrinsic = Calculate2DCoordinateGradient(ObliqueView_Outcomes_Extrinsic)
###############################################################################

max_stars_repo_name: saschweitzer/camera-calibration

max_issues_repo_name: saschweitzer/camera-calibration

max_forks_repo_name: saschweitzer/camera-calibration


text:
"Wrapper around JavaSparkContext"
mutable struct SparkContext
    jsc::JJavaSparkContext
    master::AbstractString
    appname::AbstractString
    tempdir::AbstractString
end

"""
Params:
 * master - address of application master. Currently only local and standalone modes
            are supported. Default is 'local'
 * appname - name of application
"""
function SparkContext(;master::AbstractString="local",
                      deploymode::AbstractString="client",
                      appname::AbstractString="Julia App on Spark",
                      conf::SparkConf=SparkConf(SPARK_DEFAULT_PROPS))
    setmaster(conf, master)
    setappname(conf, appname)
    setdeploy(conf, deploymode)
    jsc = JJavaSparkContext((JSparkConf,), conf.jconf)
    sc = SparkContext(jsc, master, appname, "")
    add_jar(sc, joinpath(dirname(@__FILE__), "..", "jvm", "sparkjl", "target", "sparkjl-0.1.jar"))
    return sc
end

function SparkContext(jsc::JJavaSparkContext)
    appname = jcall(jsc, "appName", JString, ())
    master = jcall(jsc, "master", JString, ())
    sc = SparkContext(jsc, master, appname, "")
    return sc
end

function get_temp_dir(sc::SparkContext)
    if sc.tempdir == ""
        sc.tempdir = mktempdir()
        x = deepcopy(sc.tempdir)
        atexit(()->rm(x, recursive=true))
    end
    return sc.tempdir
end

function Base.show(io::IO, sc::SparkContext)
    print(io, "SparkContext($(sc.master),$(sc.appname))")
end

"Close SparkContext"
function close(sc::SparkContext)
    jcall(sc.jsc, "close", Nothing, ())
end

"Add JAR file to SparkContext. Classes from this JAR will then be available to all tasks"
function add_jar(sc::SparkContext, path::AbstractString)
    jrdd = jcall(sc.jsc, "addJar", Nothing, (JString,), path)
end

"Add file to SparkContext. This file will be downloaded to each executor's work directory"
function add_file(sc::SparkContext, path::AbstractString)
    # isfile(path) || error("File `$path` doesn't exist")
    jrdd = jcall(sc.jsc, "addFile", Nothing, (JString,), path)
end

"Create RDD from a text file"
function text_file(sc::SparkContext, path::AbstractString)
    jrdd = jcall(sc.jsc, "textFile", JJavaRDD, (JString,), path)
    return JavaRDD(jrdd)
end

function default_parallelism(sc::SparkContext)
    scala_sc = jcall(sc.jsc, "sc", JSparkContext, ())
    return jcall(scala_sc, "defaultParallelism", jint, ())
end

function parallelize(sc::SparkContext, coll; n_split::Int=-1)
    n_split = n_split == -1 ? default_parallelism(sc) : n_split
    tmp_path, f = mktemp()
    dump_stream(f, coll)
    close(f)
    jrdd = jcall(JJuliaRDD, "readRDDFromFile", JJavaRDD,
                 (JJavaSparkContext, JString, jint),
                 sc.jsc, tmp_path, n_split)
    rm(tmp_path)
    JavaRDD(jrdd)
end

function set_log_level(sc::SparkContext, level::AbstractString)
    jcall(sc.jsc, "setLogLevel", Nothing, (JString,), level)
end