text:
using Gen

include("../gp_structure/shared.jl")

@gen function covariance_prior()
    node_type = @trace(categorical(node_dist), :type)

    if node_type == CONSTANT
        param = @trace(uniform_continuous(0, 1), :param)
        node = Constant(param)

    # linear kernel
    elseif node_type == LINEAR
        param = @trace(uniform_continuous(0, 1), :param)
        node = Linear(param)

    # squared exponential kernel
    elseif node_type == SQUARED_EXP
        length_scale= @trace(uniform_continuous(0, 1), :length_scale)
        node = SquaredExponential(length_scale)

    # periodic kernel
    elseif node_type == PERIODIC
        scale = @trace(uniform_continuous(0, 1), :scale)
        period = @trace(uniform_continuous(0, 1), :period)
        node = Periodic(scale, period)

    # plus combinator
    elseif node_type == PLUS
        left = @trace(covariance_prior(), :left)
        right = @trace(covariance_prior(), :right)
        node = Plus(left, right)

    # times combinator
    elseif node_type == TIMES
        left = @trace(covariance_prior(), :left)
        right = @trace(covariance_prior(), :right)
        node = Times(left, right)

    # unknown node type
    else
        error("Unknown node type: $node_type")
    end

    return node
end

@gen function model(xs::Vector{Float64})
    n = length(xs)

    # sample covariance function
    covariance_fn::Node = @trace(covariance_prior(), :tree)

    # sample diagonal noise
    noise = @trace(gamma(1, 1), :noise) + 0.01

    # compute covariance matrix
    cov_matrix = compute_cov_matrix_vectorized(covariance_fn, noise, xs)

    # sample from multivariate normal
    @trace(mvnormal(zeros(n), cov_matrix), :ys)

    return covariance_fn
end

##############
# noise move #
##############

@gen function noise_proposal(prev_trace)
    @trace(gamma(1, 1), :noise)
end

noise_move(trace) = metropolis_hastings(trace, noise_proposal, ())[1]

########################
# replace subtree move #
########################

@gen function pick_random_node_path(node::Node, path::Vector{Symbol})
    if isa(node, LeafNode)
        @trace(bernoulli(1), :done)
        path
    elseif @trace(bernoulli(0.5), :done)
        path
    elseif @trace(bernoulli(0.5), :recurse_left)
        push!(path, :left)
        @trace(pick_random_node_path(node.left, path), :left)
    else
        push!(path, :right)
        @trace(pick_random_node_path(node.right, path), :right)
    end
end
@gen function subtree_proposal(prev_trace)
    prev_subtree_node::Node = get_retval(prev_trace)
    (path::Vector{Symbol}) = @trace(pick_random_node_path(prev_subtree_node, Symbol[]), :choose_subtree_root)
    new_subtree_node::Node = @trace(covariance_prior(), :subtree) # mixed discrete / continuous
    (path, new_subtree_node)
end

@transform subtree_involution (model_in, aux_in) to (model_out, aux_out) begin

    (path::Vector{Symbol}, new_subtree_node) = @read(aux_in[], :discrete)

    # populate backward assignment with choice of root
    @copy(aux_in[:choose_subtree_root], aux_out[:choose_subtree_root])

    # swap subtrees
    model_subtree_addr = isempty(path) ? :tree : (:tree => foldr(=>, path))
    @copy(aux_in[:subtree], model_out[model_subtree_addr])
    @copy(model_in[model_subtree_addr], aux_out[:subtree])
end

is_involution!(subtree_involution)

replace_subtree_move(trace) = metropolis_hastings(
    trace, subtree_proposal, (), subtree_involution; check=true)[1]

#####################
# inference program #
#####################

function inference(xs::Vector{Float64}, ys::Vector{Float64}, num_iters::Int, callback)

    # observed data
    constraints = choicemap()
    constraints[:ys] = ys

    # generate initial trace consistent with observed data
    (trace, _) = generate(model, (xs,), constraints)

    # do MCMC
    local covariance_fn::Node
    local noise::Float64
    for iter=1:num_iters

        covariance_fn = get_retval(trace)
        noise = trace[:noise]
        #callback(covariance_fn, noise)

        # do MH move on the subtree
        trace = replace_subtree_move(trace)

        # do MH move on the top-level white noise
        trace = noise_move(trace)
    end
    (covariance_fn, noise)
end

function experiment()

    # load and rescale the airline dataset
    (xs, ys) = get_airline_dataset()
    xs_train = xs[1:100]
    ys_train = ys[1:100]
    xs_test = xs[101:end]
    ys_test = ys[101:end]

    # set seed
    Random.seed!(0)

    # print MSE and predictive log likelihood on test data
    callback = (covariance_fn, noise) -> begin
        pred_ll = predictive_ll(covariance_fn, noise, xs_train, ys_train, xs_test, ys_test)
        mse = compute_mse(covariance_fn, noise, xs_train, ys_train, xs_test, ys_test)
        println("mse: $mse, predictive log likelihood: $pred_ll")
    end

    # do inference, time it
    @time (covariance_fn, noise) = inference(xs_train, ys_train, 1000, callback)

    # sample predictions
    pred_ys = predict_ys(covariance_fn, noise, xs_train, ys_train, xs_test)
end

@time experiment()
@time experiment()

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
# used in shortest path calculations
# has_distances{T}(distmx::AEdgeMap{T}) =
#     issparse(distmx)? (nnz(distmx) > 0) : !isempty(distmx)

"""
    eccentricity(g, v, distmx=weights(g))

Calculates the eccentricity[ies] of a vertex `v`,
An optional matrix of edge distances may be supplied.

The eccentricity of a vertex is the maximum shortest-path distance between it
and all other vertices in the graph.
"""
function eccentricity(
    g::AGraphOrDiGraph,
    v::Int,
    distmx::AEdgeMap=weights(g)
)
    e = maximum(dijkstra_shortest_paths(g,v,distmx).dists)
    e == typemax(valtype(distmx)) && error("Infinite path length detected")

    return e
end

"""

    eccentricities(g, distmx=weights(g))
    eccentricities(g, vs, distmx=weights(g))

Returns `[eccentricity(g,v,distmx) for v in vs]`. When `vs` it is not supplied,
considers all node in the graph.

See also [`eccentricity`](@ref).

Note: the eccentricity vector returned by `eccentricity` may be eventually used as input
in some eccentricity related measures ([`periphery`](@ref), [`center`](@ref)).
"""
function eccentricities(
    g::AGraphOrDiGraph,
    vs::AbstractVector,
    distmx::AEdgeMap=weights(g)
)
    [eccentricity(g,v,distmx) for v in vs]
end

eccentricities(
    g::AGraphOrDiGraph,
    distmx::AEdgeMap=weights(g)
) = eccentricities(g, 1:nv(g), distmx)

"""
    diameter(g, distmx=weights(g))

Returns the maximum distance between any two vertices in `g`.
Distances  between two adjacent nodes are given by `distmx`.

See also [`eccentricities`](@ref), [`radius`](@ref).
"""
diameter(g::AGraphOrDiGraph, distmx::AEdgeMap = weights(g)) =
    maximum(eccentricities(g, distmx))

"""
    radius(g, distmx=weights(g))

Returns the minimum distance between any two vertices in `g`.
Distances  between two adjacent nodes are given by `distmx`.

See [`eccentricities`](@ref), [`diameter`](@ref).
"""
radius(g::AGraphOrDiGraph, distmx=weights(g)) =
    minimum(eccentricities(g, distmx))

"""
    periphery(g, distmx=weights(g))
    periphery(all_ecc)

Returns the set of all vertices whose eccentricity is equal to the graph's
diameter (that is, the set of vertices with the largest eccentricity).

Eventually a vector `all_ecc` contain the eccentricity of each node
can be passed as argument.

See [`eccentricities`](@ref).
"""
function periphery(all_e::Vector)
    diam = maximum(all_e)
    return filter((x)->all_e[x] == diam, 1:length(all_e))
end

periphery(g::AGraphOrDiGraph, distmx::AEdgeMap=weights(g)) =
    periphery(eccentricities(g, distmx))

"""
    center(g, distmx=weights(g))
    center(all_ecc)

Returns the set of all vertices whose eccentricity is equal to the graph's
radius (that is, the set of vertices with the smallest eccentricity).

Eventually a vector `all_ecc` contain the eccentricity of each node
can be passed as argument.

See [`eccentricities`](@ref).
"""
function center(all_e::Vector)
    rad = minimum(all_e)
    return filter((x)->all_e[x] == rad, 1:length(all_e))
end

center(g::AGraphOrDiGraph, distmx::AEdgeMap = weights(g)) =
    center(eccentricities(g, distmx))

ext: jl


text:
const PM_BUSTYPES = Dict{PSY.BusTypes, Int}(
    PSY.BusTypes.ISOLATED => 4,
    PSY.BusTypes.PQ => 1,
    PSY.BusTypes.PV => 2,
    PSY.BusTypes.REF => 3,
    PSY.BusTypes.SLACK => 3,
)

const ACCEPTED_PM_DATA_KWARGS = [:name, :start_time, :time_periods, :period]

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
#   Unit "simulations.jl" of the PosDefManifoldML Package for Julia language
#   v 0.4.5 - last update June 26 2021
#
#   MIT License
#   Copyright (c) 2019-2021,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   Main module of the PosDefManifoldML package.
#
#   This pachwge works in conjunction with the PosDefManifold package
#   https://github.com/Marco-Congedo/PosDefManifold.jl

 __precompile__()

module PosDefManifoldML

using LinearAlgebra, Base.Threads, Random, StatsBase, Statistics, PosDefManifold
using Dates:now
using GLMNet:GLMNet, glmnet, glmnetcv, GLMNetPath, GLMNetCrossValidation
using Distributions:Distributions, Binomial
using LIBSVM: svmpredict, svmtrain, SVC, NuSVC, OneClassSVM, NuSVR, EpsilonSVR,
      LinearSVC, Linearsolver, Kernel

# Special instructions and variables
BLAS.set_num_threads(Sys.CPU_THREADS)

# constants #

const titleFont     = "\x1b[32m"
const separatorFont = "\x1b[92m"
const defaultFont   = "\x1b[0m"
const greyFont      = "\x1b[90m"
const 📌            = titleFont*"PosDefManifoldML"*defaultFont
const dice = ("⚀", "⚁", "⚂", "⚃", "⚄", "⚅")

# shortcut to LIBSVM kernels enum type
const Linear 		= Kernel.KERNEL(0)
const Polynomial 	= Kernel.KERNEL(1)
const RadialBasis 	= Kernel.KERNEL(2)
const Sigmoid 		= Kernel.KERNEL(3)
const Precomputed 	= Kernel.KERNEL(4)

# types #
abstract type MLmodel end # all machine learning models
abstract type PDmodel<:MLmodel end # PD manifold models
abstract type TSmodel<:MLmodel end # tangent space models

IntVector=Vector{Int}

import Base:show
#import GLMNet.predict
#import Distributions.fit
# import LIBSVM.predict
import StatsBase:fit, predict

export

    # from this module
    MLmodel,
    PDmodel,
    TSmodel,
    IntVector,
	Linear,
	Polynomial,
	RadialBasis,
	Sigmoid,
	Precomputed,

    # from mdm.jl
    MDMmodel,
    MDM,
    fit,
    predict,
    getMean,
    getDistances,

    # from enlr.jl
    ENLRmodel,
    ENLR,

	# from libSVM.jl
	SVMmodel,
	SVM,
	SVC,
	NuSVC,
	OneClassSVM,
	NuSVR,
	EpsilonSVR,
	LinearSVC,
	Linearsolver,
	Kernel,

    # from cv.jl
    CVacc,
    cvAcc,
    cvSetup,

    # from tools.jl
    tsMap,
    tsWeights,
    gen2ClassData,
	confusionMat,
	predictAcc,
    predictErr,
	rescale!

include("mdm.jl")
include("enlr.jl")
include("cv.jl")
include("tools.jl")
include("svm.jl")

println("\n⭐ "," Welcome to the ", 📌," package", " ⭐\n")
@info " "
println(" Your Machine `",gethostname(),"` (",Sys.MACHINE, ")")
println(" runs on kernel ",Sys.KERNEL," with word size ",Sys.WORD_SIZE,".")
println(" CPU  Threads: ", Sys.CPU_THREADS)
println(" Base.Threads: ", "$(Threads.nthreads())")
println(" BLAS Threads: ", "$(Sys.CPU_THREADS)", "\n")

end # module


text:
module TestSugar

using Atomix: @atomic, @atomicreplace, @atomicswap
using Test

function test_get()
    A = [42]
    @test (@atomic A[1]) === 42
    @test (@atomic A[end]) === 42
    @test (@atomic :monotonic A[begin]) === 42
    order = :monotonic
    @test (@atomic order A[begin]) === 42
end

function test_get_2d()
    A = view([11 12; 21 22], 1:2, 1:2)
    @test IndexStyle(A) isa IndexCartesian
    @test (@atomic A[1]) === 11
    @test (@atomic A[2]) === 21
    @test (@atomic A[end]) === 22
    @test (@atomic A[2, 1]) === 21
    @test (@atomic A[end, 1]) === 21
    @test (@atomic A[1, end]) === 12
end

function test_set()
    A = [42, 43]
    @atomic A[1] = 123
    @atomic A[end] = 124
    @test A[1] === 123
    @test A[end] === 124
end

function test_inc()
    A = [1, 1]
    @test (@atomic A[1] += 123) === 124
    @test A[1] === 124
    @test (@atomic A[end] += 123) === 124
    @test A[end] === 124
end

function test_swap()
    A = [1, 1]
    @test (@atomicswap A[1] = 123) === 1
    @test A[1] === 123
    @test (@atomicswap A[end] = 456) === 1
    @test A[end] === 456
end

function test_cas()
    A = [1, 1]
    @test (@atomicreplace A[1] 1 => 123) == (old = 1, success = true)
    @test A[1] === 123
    @test (@atomicreplace A[1] 1 => 456) == (old = 123, success = false)
    @test A[1] === 123
    @test (@atomicreplace A[end] 1 => 789) == (old = 1, success = true)
    @test A[end] === 789
    update = 789 => 123
    @test (@atomicreplace A[end] update) == (old = 789, success = true)
    @test A[end] === 123
end

end  # module


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

using Markdown
using InteractiveUtils

# ╔═╡ db039f9e-8787-11eb-3037-6566e5cd35e1
begin
  using Pkg; Pkg.activate("MLJ_env", shared=true)
  using PlutoUI

    using CategoricalArrays
    import MLJLinearModels
    import DataFrames
    import CSV
    import DecisionTree
    using MLJ
    import MLJLinearModels
    import MultivariateStats
    import MLJFlux
    using Plots
	using Test
end

# ╔═╡ 922be298-878e-11eb-16a8-35f47ace05e1
begin
	using Clustering
	using MLJClusteringInterface
end

# ╔═╡ 303866da-8791-11eb-0872-b908cba30c10
using Random

# ╔═╡ b0de5b8e-8787-11eb-1a63-e94d308ecf79
md"""
## A workshop introducing the machine learning toolbox MLJ - Part3 - Transformers and Pipelines

  - source:[MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
  - based on [Machine Learning in Julia using MLJ, JuliaCon2020](https://github.com/ablaom/MachineLearningInJulia2020)
"""

# ╔═╡ dae6f0a6-8787-11eb-18be-5bee665ff211
md"""
##### General resources

  - [List of methods introduced in this tutorial](methods.md)
  - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
  - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
  - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
  - [Data Science Tutorials in Julia](https://alan-turing-institute.github.io/DataScienceTutorials.jl/)
"""

# ╔═╡ dacb2c5e-8787-11eb-0f18-03004883918e
html"""

"""

# ╔═╡ daadfbca-8787-11eb-1e8b-df76827065f3
md"""
##### TOC
  - [Transformers](#transformers)
  - [More transformers](#more-transformers)
  - [Pipelines](#pipelines)
  - [Training of composite models is "smart"](#training)
  - [Incorporating target transformations](#target-transfo)
  - [Resources or this part](#resources)
  - [Exercises for Part 3](#exercises)
"""

# ╔═╡ 336759dc-8788-11eb-2f1f-47edaaa7eeb6
html"""

  back to TOC

"""

# ╔═╡ 3468ea44-8788-11eb-0043-c777b00c423c
md"""
#### Transformers

Unsupervised models, which receive no target `y` during training,
always have a `transform` operation. They sometimes also support an
`inverse_transform` operation, with obvious meaning, and sometimes
support a `predict` operation (see the clustering example discussed
[here](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/#Transformers-that-also-predict-1)).
Otherwise, they are handled much like supervised models.

Here's a simple standardization example:
"""

# ╔═╡ 66ad95f4-8788-11eb-0d86-21ccf44a07ad
begin
	xᵣ = rand(100);
	@show (μ=mean(xᵣ), σ=std(xᵣ));
end

# ╔═╡ 6692ca9e-8788-11eb-0e22-23255f139aa5
begin
	modelᵣ = Standardizer() # a built-in model

	machᵣ = machine(modelᵣ, xᵣ)
	fit!(machᵣ)
	x̂ᵣ = transform(machᵣ, xᵣ);
	@show (μ=mean(x̂ᵣ,), σ=std(x̂ᵣ));
end

# ╔═╡ 667ae6fe-8788-11eb-37e0-937840f15342
md"""
This particular model has an `inverse_transform`:
"""

# ╔═╡ 665e3cac-8788-11eb-28ab-e1a34482956c
inverse_transform(machᵣ, x̂ᵣ) ≈ xᵣ

# ╔═╡ 662d8c94-8788-11eb-1cc2-ebadc5980bc4
md"""
### Re-encoding the King County House data as continuous

For further illustrations of transformers, let's re-encode *all* of the
King County House input features (see **[Ex
3/Part 1]**) into a set of `Continuous`
features. We do this with the `ContinuousEncoder` model, which, by
default, will:

  - one-hot encode all `Multiclass` features
  - coerce all `OrderedFactor` features to `Continuous` ones
  - coerce all `Count` features to `Continuous` ones (there aren't any)
  - drop any remaining non-Continuous features (none of these either)

First, we reload the data and fix the scitypes (Exercise 3):
"""

# ╔═╡ 66140f1a-8788-11eb-1b90-138a63cdc086
begin
	const DIR="."
	hfile = CSV.File(joinpath(DIR, "data", "house.csv"));

	house = DataFrames.DataFrame(hfile)
	coerce!(house, autotype(hfile))
	coerce!(house, Count => Continuous, :zipcode => Multiclass);
	schema(house)
end

# ╔═╡ da835c8e-8788-11eb-3905-0b961fb28506
yₕ, Xₕ = unpack(house, ==(:price), name -> true, rng=123);

# ╔═╡ da5ab4fa-8788-11eb-1aa3-c1a921e861c8
md"""
Instantiate the unsupervised model (transformer):
"""

# ╔═╡ da285618-8788-11eb-0635-8dcc0c9b5968
encoder = ContinuousEncoder() # a built-in model; no need to @load it

# ╔═╡ fe5fd3bc-8788-11eb-095a-29f563116783
md"""
Bind the model to the data and fit! $(html"
")
Then, Transform and inspect the result:
"""

# ╔═╡ fe3c4352-8788-11eb-02ab-a39f52eed4b8
begin
	machₕ = machine(encoder, Xₕ) |> fit!;
	Xcont = transform(machₕ, Xₕ);
	schema(Xcont)
end

# ╔═╡ 45ad2292-8789-11eb-1b70-e7df1c2198b4
html"""

  back to TOC

"""

# ╔═╡ 2688f3aa-8789-11eb-2ccc-57a834adb9d2
md"""
#### More transformers

Here's how to list all of MLJ's unsupervised models:
"""

# ╔═╡ fe1da1ce-8788-11eb-0235-47ac875ab7fb
models(m -> !m.is_supervised)

# ╔═╡ fe055a04-8788-11eb-2b38-2b23cbae597c
md"""
Some commonly used ones are built-in (do not require `@load`ing):

model type                  | does what?
----------------------------|----------------------------------------------
ContinuousEncoder | transform input table to a table of `Continuous` features (see above)
FeatureSelector | retain or dump selected features
FillImputer | impute missing values
OneHotEncoder | one-hot encoder `Multiclass` (and optionally `OrderedFactor`) features
Standardizer | standardize (whiten) a vector or all `Continuous` features of a table
UnivariateBoxCoxTransformer | apply a learned Box-Cox transformation to a vector
UnivariateDiscretizer | discretize a `Continuous` vector, and hence render its elscitypw `OrderedFactor`

In addition to "dynamic" transformers (ones that learn something
from the data and must be `fit!`) users can wrap ordinary functions
as transformers, and such *static* transformers can depend on
parameters, like the dynamic ones. See
[here](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/#Static-transformers-1)
for how to define your own static transformers.
"""

# ╔═╡ fdeaedae-8788-11eb-3dab-f16f800d2b77
html"""

  back to TOC

"""

# ╔═╡ da147a30-8788-11eb-01d0-ed61ac48ddee
md"""
#### Pipelines
"""

# ╔═╡ 65fa258c-8788-11eb-00c6-9f0b8f477d5d
length(schema(Xcont).names)

# ╔═╡ 013a548a-878a-11eb-160a-1da458bd2324
md"""
Let's suppose that additionally we'd like to reduce the dimension of our data.  A model that will do this is `PCA` from `MultivariateStats`:
"""

# ╔═╡ 01205ec2-878a-11eb-2fad-f9b213f4560f
# Pkg.add("MLJMultivariateStatsInterface")

# ╔═╡ 0103cd2a-878a-11eb-34df-8fb2e0261c30
reducer = @load PCA

# ╔═╡ 00ec055a-878a-11eb-3dac-1d0f6e10ea8c
md"""
Now, rather than simply repeating the work-flow above, applying the new
transformation to `Xcont`, we can combine both the encoding and the
dimension-reducing models into a single model, known as a
*pipeline*.

While MLJ offers a powerful interface for composing
models in a variety of ways, we'll stick to these simplest class of
composite models for now. The easiest way to construct them is using
the `@pipeline` macro:
"""

# ╔═╡ 00d28b20-878a-11eb-108d-99e7981ae462
pipe = @pipeline encoder reducer

# ╔═╡ 266cd818-878a-11eb-34d6-c173735c9604
md"""
Notice that `pipe` is an *instance* of an automatically generated type (called `Pipeline<some digits>`).

The new model behaves like any other transformer:
"""

# ╔═╡ 2651f002-878a-11eb-2f30-71384dae63d3
begin
	machₚ = machine(pipe, Xₕ)
	fit!(machₚ)
	Xsmallₚ = transform(machₚ, Xₕ)
	schema(Xsmallₚ)
end

# ╔═╡ 2637ca10-878a-11eb-07be-af6d5a769213
md"""
Want to combine this pre-processing with ridge regression?
"""

# ╔═╡ 67bde8de-878a-11eb-19d9-4176c3871746
@load RidgeRegressor pkg=MLJLinearModels

# ╔═╡ 67a2c568-878a-11eb-2c52-1f3503de2742
rregr = MLJLinearModels.RidgeRegressor()

# ╔═╡ 6788dc0e-878a-11eb-1552-0d7d020eaa2c
pipe₂ = @pipeline encoder reducer rregr

# ╔═╡ 77a77936-878a-11eb-255c-116e90f96738
md"""
Now our pipeline is a supervised model, instead of a transformer, whose performance we can evaluate:
"""

# ╔═╡ 77885420-878a-11eb-08a8-43f03f192651
begin
	mach₂ = machine(pipe₂, Xₕ, yₕ)
	evaluate!(mach₂, measure=mae, resampling=Holdout()) # CV(nfolds=6) is default
end

# ╔═╡ abebfdde-878a-11eb-09ad-b1440c37497c
html"""

  back to TOC

"""

# ╔═╡ 775acdfc-878a-11eb-0007-5519985c1642
md"""
#### Training of composite models is "smart"

Now notice what happens if we train on all the data, then change a regressor hyper-parameter and retrain:
"""

# ╔═╡ 6776361a-878a-11eb-3fca-b721f6de3c4b
fit!(mach₂);

# ╔═╡ 261a54a8-878a-11eb-0b82-a55b8f646bb2
begin
	pipe₂.ridge_regressor.lambda = 0.1
	fit!(mach₂)
end

# ╔═╡ 235d9a16-878a-11eb-219f-7750704b1e0b
md"""
Second time only the ridge regressor is retrained!

Mutate a hyper-parameter of the `PCA` model and every model except the `ContinuousEncoder` (which comes before it will be retrained):
"""

# ╔═╡ 235c6dbc-878a-11eb-22b2-1f7f9f4fec5b
begin
	pipe₂.pca.pratio = 0.9999
	fit!(mach₂)
end

# ╔═╡ fa99cd76-878a-11eb-00f5-170c7825df20
md"""
##### Inspecting composite models

The dot syntax used above to change the values of *nested* hyper-parameters is also useful when inspecting the learned parameters and report generated when training a composite model:
"""

# ╔═╡ fa8042de-878a-11eb-10db-35eb7404e05b
fitted_params(mach₂).ridge_regressor

# ╔═╡ fa613254-878a-11eb-0d6a-1f07472f012c
report(mach₂).pca

# ╔═╡ fa4524ec-878a-11eb-2754-cd5ca0dcee5d
html"""

  back to TOC

"""

# ╔═╡ fa29314c-878a-11eb-2fce-e95e75060931
md"""
##### Incorporating target transformations

Next, suppose that instead of using the raw `:price` as the training target, we want to use the log-price (a common practice in dealing with house price data). However, suppose that we still want to report final *predictions* on the original linear scale (and use these for evaluation purposes). Then we supply appropriate functions to key-word arguments `target` and `inverse`.

First we'll overload `log` and `exp` for broadcasting:
"""

# ╔═╡ 235c47bc-878a-11eb-3511-7b85a574207a
begin
	Base.log(v::AbstractArray) = log.(v)
	Base.exp(v::AbstractArray) = exp.(v)
end

# ╔═╡ 5282e4a0-878b-11eb-36df-73cdf38aaa7e
md"""
Now for the new pipeline:
"""

# ╔═╡ 52629a88-878b-11eb-3cf0-77bc93decb35
begin
	pipe₃ = @pipeline encoder reducer rregr target=log inverse=exp
	mach₃ = machine(pipe₃, Xₕ, yₕ)
	evaluate!(mach₃, measure=mae);
end

# ╔═╡ 52499fb0-878b-11eb-191c-7f89628d3a85
md"""
MLJ will also allow you to insert *learned* target transformations. For example, we might want to apply `Standardizer()` to the target, to standardize it, or `UnivariateBoxCoxTransformer()` to make it look Gaussian. Then instead of specifying a *function* for `target`, we specify an unsupervised *model* (or model type). One does not specify `inverse` because only models implementing `inverse_transform` are allowed.

Let's see which of these two options results in a better outcome:
"""

# ╔═╡ 5230d4dc-878b-11eb-2f83-e160aa434041
begin
	box = UnivariateBoxCoxTransformer(n=20)
	stand = Standardizer()
end

# ╔═╡ 52178a98-878b-11eb-191d-a3617b11098c
begin
	pipe₄ = @pipeline encoder reducer rregr target=box
	mach₄ = machine(pipe₄, Xₕ, yₕ)
	evaluate!(mach₄, measure=mae);
end

# ╔═╡ 51ff4a0a-878b-11eb-1663-3175d3023d3d
begin
	pipe₄.target = stand
	evaluate!(mach₄, measure=mae);
end

# ╔═╡ c419e15e-878b-11eb-12ae-f321152db3b4
html"""

  back to TOC

"""

# ╔═╡ c18d4b80-878b-11eb-03df-57554a304eac
md"""
#### Resources for Part 3

- From the MLJ manual:
    - [Transformers and other unsupervised models](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/)
    - [Linear pipelines](https://alan-turing-institute.github.io/MLJ.jl/dev/composing_models/#Linear-pipelines-1)
- From Data Science Tutorials:
    - [Composing models](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/composing-models/)
"""

# ╔═╡ 4ae60c30-8791-11eb-3547-051e3ec0269e
html"""

  back to TOC

"""

# ╔═╡ c178cd36-878b-11eb-019a-ab7f9371e0f0
md"""
#### Exercises for Part 3

##### Exercise 1

Consider again the Horse Colic classification problem considered in Exercise 6, but with all features, `Finite` and `Infinite`:
"""

# ╔═╡ 117c134a-878c-11eb-0062-910219a6efb9
begin
  h_file = CSV.File(joinpath(DIR, "data", "horse.csv"));
  horse = DataFrames.DataFrame(h_file); # convert to data frame without copying columns
end

# ╔═╡ 22e61c98-878c-11eb-07f6-075b60e2979a
begin
  coerce!(horse, autotype(horse));
  coerce!(horse, Count => Continuous);
  coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);
  schema(horse)
end

# ╔═╡ c15af1f6-878b-11eb-207b-ff38182ab074
begin
	y₁, X₁ = unpack(horse, ==(:outcome), name -> true);
	schema(X₁)
end

# ╔═╡ 51e39daa-878b-11eb-2ac7-ff55af65453d
md"""
(a) Define a pipeline that:
- uses `Standardizer` to ensure that features that are already
  continuous are centered at zero and have unit variance
- re-encodes the full set of features as `Continuous`, using
  `ContinuousEncoder`
- uses the `KMeans` clustering model from `Clustering.jl`
  to reduce the dimension of the feature space to `k=10`.
- trains a `EvoTreeClassifier` (a gradient tree boosting
  algorithm in `EvoTrees.jl`) on the reduced data, using
  `nrounds=50` and default values for the other
   hyper-parameters

(b) Evaluate the pipeline on all data, using 6-fold cross-validation
and `cross_entropy` loss.

(c) Plot a learning curve which examines the effect on this loss
as the tree booster parameter `max_depth` varies from 2 to 10.
"""

# ╔═╡ 69ab3d78-878e-11eb-3b38-bf4f92b0f82a
# Pkg.add("Clustering")
# Pkg.add("MLJClusteringInterface")

# ╔═╡ 98847588-878e-11eb-1394-75d948d73c5c
@load KMeans pkg=Clustering

# ╔═╡ ab66fdd8-878e-11eb-1187-d11889a8895b
@load EvoTreeClassifier

# ╔═╡ e39ff3da-878e-11eb-18c7-1babe2dfe6e1
# (a) part
pipe₁ = @pipeline(Standardizer,
	ContinuousEncoder,
	KMeans(k=10),
	EvoTrees.EvoTreeClassifier(nrounds=50))

# ╔═╡ 661d2bd2-878c-11eb-0565-d9d420533cfc
# (b) part
begin
	mach₁ = machine(pipe₁, X₁, y₁)
	evaluate!(mach₁, resampling=CV(nfolds=6), measure=cross_entropy)
end

# ╔═╡ 6606b832-878c-11eb-3830-a3d32f3b07bd
# (c) part
begin
	r₁ = range(pipe₁, :(evo_tree_classifier.max_depth), lower=1, upper=10)

	curve₁ = learning_curve(mach₁,
		range=r₁,
		resampling=CV(nfolds=6),
		measure=cross_entropy)
end

# ╔═╡ 65ee2404-878c-11eb-028b-2d6b037a55f4
plot(curve₁.parameter_values, curve₁.measurements,
	xlabel="max_depth",
	ylabel="CV estimate of cross entropy",
	legend=false)

# ╔═╡ 65d52fe4-878c-11eb-3ba8-3d7f2375c137
md"""
(*from solutions to exercises*) Here's a second curve using a different random seed for the booster:
"""

# ╔═╡ 65bbf126-878c-11eb-2d17-151e58fe91f1
begin
	pipe₁.evo_tree_classifier.rng = Random.MersenneTwister(123)

	curve₂ = learning_curve(mach₁,
		range=r₁,
		resampling=CV(nfolds=6),
		measure=cross_entropy)

	plot!(curve₂.parameter_values, curve₂.measurements)
end

# ╔═╡ Cell order:
# ╟─b0de5b8e-8787-11eb-1a63-e94d308ecf79
# ╠═db039f9e-8787-11eb-3037-6566e5cd35e1
# ╟─dae6f0a6-8787-11eb-18be-5bee665ff211
# ╟─dacb2c5e-8787-11eb-0f18-03004883918e
# ╟─daadfbca-8787-11eb-1e8b-df76827065f3
# ╟─336759dc-8788-11eb-2f1f-47edaaa7eeb6
# ╟─3468ea44-8788-11eb-0043-c777b00c423c
# ╠═66ad95f4-8788-11eb-0d86-21ccf44a07ad
# ╠═6692ca9e-8788-11eb-0e22-23255f139aa5
# ╟─667ae6fe-8788-11eb-37e0-937840f15342
# ╠═665e3cac-8788-11eb-28ab-e1a34482956c
# ╟─662d8c94-8788-11eb-1cc2-ebadc5980bc4
# ╠═66140f1a-8788-11eb-1b90-138a63cdc086
# ╠═da835c8e-8788-11eb-3905-0b961fb28506
# ╟─da5ab4fa-8788-11eb-1aa3-c1a921e861c8
# ╠═da285618-8788-11eb-0635-8dcc0c9b5968
# ╟─fe5fd3bc-8788-11eb-095a-29f563116783
# ╠═fe3c4352-8788-11eb-02ab-a39f52eed4b8
# ╟─45ad2292-8789-11eb-1b70-e7df1c2198b4
# ╟─2688f3aa-8789-11eb-2ccc-57a834adb9d2
# ╠═fe1da1ce-8788-11eb-0235-47ac875ab7fb
# ╟─fe055a04-8788-11eb-2b38-2b23cbae597c
# ╟─fdeaedae-8788-11eb-3dab-f16f800d2b77
# ╟─da147a30-8788-11eb-01d0-ed61ac48ddee
# ╠═65fa258c-8788-11eb-00c6-9f0b8f477d5d
# ╟─013a548a-878a-11eb-160a-1da458bd2324
# ╠═01205ec2-878a-11eb-2fad-f9b213f4560f
# ╠═0103cd2a-878a-11eb-34df-8fb2e0261c30
# ╟─00ec055a-878a-11eb-3dac-1d0f6e10ea8c
# ╠═00d28b20-878a-11eb-108d-99e7981ae462
# ╟─266cd818-878a-11eb-34d6-c173735c9604
# ╠═2651f002-878a-11eb-2f30-71384dae63d3
# ╟─2637ca10-878a-11eb-07be-af6d5a769213
# ╠═67bde8de-878a-11eb-19d9-4176c3871746
# ╠═67a2c568-878a-11eb-2c52-1f3503de2742
# ╠═6788dc0e-878a-11eb-1552-0d7d020eaa2c
# ╟─77a77936-878a-11eb-255c-116e90f96738
# ╠═77885420-878a-11eb-08a8-43f03f192651
# ╟─abebfdde-878a-11eb-09ad-b1440c37497c
# ╟─775acdfc-878a-11eb-0007-5519985c1642
# ╠═6776361a-878a-11eb-3fca-b721f6de3c4b
# ╠═261a54a8-878a-11eb-0b82-a55b8f646bb2
# ╟─235d9a16-878a-11eb-219f-7750704b1e0b
# ╠═235c6dbc-878a-11eb-22b2-1f7f9f4fec5b
# ╟─fa99cd76-878a-11eb-00f5-170c7825df20
# ╠═fa8042de-878a-11eb-10db-35eb7404e05b
# ╠═fa613254-878a-11eb-0d6a-1f07472f012c
# ╟─fa4524ec-878a-11eb-2754-cd5ca0dcee5d
# ╟─fa29314c-878a-11eb-2fce-e95e75060931
# ╠═235c47bc-878a-11eb-3511-7b85a574207a
# ╟─5282e4a0-878b-11eb-36df-73cdf38aaa7e
# ╠═52629a88-878b-11eb-3cf0-77bc93decb35
# ╟─52499fb0-878b-11eb-191c-7f89628d3a85
# ╠═5230d4dc-878b-11eb-2f83-e160aa434041
# ╠═52178a98-878b-11eb-191d-a3617b11098c
# ╠═51ff4a0a-878b-11eb-1663-3175d3023d3d
# ╠═c419e15e-878b-11eb-12ae-f321152db3b4
# ╟─c18d4b80-878b-11eb-03df-57554a304eac
# ╟─4ae60c30-8791-11eb-3547-051e3ec0269e
# ╟─c178cd36-878b-11eb-019a-ab7f9371e0f0
# ╠═117c134a-878c-11eb-0062-910219a6efb9
# ╠═22e61c98-878c-11eb-07f6-075b60e2979a
# ╠═c15af1f6-878b-11eb-207b-ff38182ab074
# ╟─51e39daa-878b-11eb-2ac7-ff55af65453d
# ╠═69ab3d78-878e-11eb-3b38-bf4f92b0f82a
# ╠═922be298-878e-11eb-16a8-35f47ace05e1
# ╠═98847588-878e-11eb-1394-75d948d73c5c
# ╠═ab66fdd8-878e-11eb-1187-d11889a8895b
# ╠═e39ff3da-878e-11eb-18c7-1babe2dfe6e1
# ╠═661d2bd2-878c-11eb-0565-d9d420533cfc
# ╠═6606b832-878c-11eb-3830-a3d32f3b07bd
# ╠═65ee2404-878c-11eb-028b-2d6b037a55f4
# ╟─65d52fe4-878c-11eb-3ba8-3d7f2375c137
# ╠═303866da-8791-11eb-0872-b908cba30c10
# ╠═65bbf126-878c-11eb-2d17-151e58fe91f1

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
import TrajectoryOptimization
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics/urdf")
urdf_doublependulum = joinpath(urdf_folder, "doublependulum.urdf")

acrobot_model = Model(urdf_doublependulum,[0.;1.]) # underactuated, only control for second "elbow" joint

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
mutable struct Portfolio
    id::UUID
    owner::String
    families::Vector{Family}
    Portfolio(o::String, f::Vector{Family}) = new(UUIDs.uuid4(), o, f)
end

function Base.show(io::IO, p::Portfolio)
    print(io, "Portfolio of $(owner(p)) with $(size(p)) families" )
end

owner(p::Portfolio) = p.owner
Base.size(p::Portfolio) = length(p.families)
families(p::Portfolio) = p.families

function diversity(p::Portfolio, levelfun)
    fams = families(p)
	diversity(fams, levelfun)
end

# struct PortfolioSummary
#     size::Int
#     # Think of useful summaries
# end


text:
using DelimitedFiles, Plots
pyplot()

labels = ["Planet b","c","d","e","f","g","h"]
nbg = readdlm("nbg_tts.txt", ',')[:,1:3000]
reb = readdlm("reb_tts.txt", ',')[:,1:3000]
fast = readdlm("ttvfast_tts.txt", ',')[:,1:3000]

reb_diff = nbg .- reb
plot()
for i in 1:7
    ind = findfirst(isequal(0.0),reb_diff[i,:])
    ind = ind === nothing ? ind = 3000 : ind = ind - 5
    plot!(nbg[i,1:ind], reb_diff[i,1:ind] .* 86400, label=labels[i])
end
ylabel!("Difference in transit times [sec]")
xlabel!("Time [Days]")
savefig("nbg_vs_reb_transit_times.pdf")

fast_diff = nbg .- fast
plot()
for i in 1:7
    ind = findfirst(isequal(0.0),fast_diff[i,:])
    ind = ind === nothing ? ind = 3000 : ind = ind - 5
    plot!(nbg[i,1:ind], fast_diff[i,1:ind] .* 86400, label=labels[i])
end
ylabel!("Difference in transit times [sec]")
xlabel!("Time [Days]")
savefig("nbg_vs_ttvfast_transit_times.pdf")

max_forks_repo_name: langfzac/nbg-papers


text:
import ThreeJS

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
    w = Signal(1.0)
    h = Signal(1.0)
    d = Signal(1.0)
    rx = Signal(0.0)
    ry = Signal(0.0)
    rz = Signal(0.0)
    o = Signal(1.0)
    vbox(
        "ThreeJS Example",
        vskip(2em),
        hbox(
            vbox(
                "Sides",
                hbox("width",slider(1.0:5.0) >>> w),
                hbox("height",slider(1.0:5.0) >>> h),
                hbox("depth",slider(1.0:5.0) >>> d),
            ),
            hskip(2em),
            vbox(
                "Rotations",
                hbox("x", slider(0.0:5.0:360.0) >>> rx),
                hbox("y", slider(0.0:5.0:360.0) >>> ry),
                hbox("z", slider(0.0:5.0:360.0) >>> rz),
            ),
            hbox("opacity",slider(0.0:0.1:1.0, value = 1.0) >>> o),
        ),
        vskip(2em),
        map(w, h, d, rx, ry, rz, o) do w, h, d, rx, ry, rz, o
        ThreeJS.outerdiv() <<
            (ThreeJS.initscene() <<
                [
                    ThreeJS.mesh(0.0, 0.0, 0.0, rx = rx, ry = ry, rz = rz) <<
                    [
                        ThreeJS.box(w, h, d), ThreeJS.material(Dict(:kind=>"lambert",:color=>"red", :transparent=>true, :opacity=>o))
                    ],
                    ThreeJS.pointlight(10.0, 10.0, 10.0),
                    ThreeJS.camera(0.0, 0.0, 10.0)
                ]
            )
        end
    ) |> pad(2em)
end


text:
using TestPackage
using Test

@testset "TestPackage.jl" begin
    # integration parameters
    a = 0
    b = 1
    n = 1000

    # Error tolerance for testing
    # total_error < interval_width / (n_samples * 0.1)
    total_error = (b - a) / (n * 0.1)

    f(x) = x^2
    F(x) = x^3/3
    res = F(b) - F(a) # actual result

    lr = left_riemann(f, a, b, n)
    rr = right_riemann(f, a, b, n)
    tr = trapezoid(f, a, b, n)

    @test abs(lr - res) < total_error
    @test abs(rr - res) < total_error
    @test abs(tr - res) < total_error

end


text:
using Logging, LeafAreaIndex, ParallelDataTransfer, CoordinateTransformations
using PyCall
# TODO use MicroLogging on julia 0.6 or Base.Logging on julia 1.0
using Images #also imports FileIO for reading jpg
using Statistics
import StatsBase, JLD2, FileIO

CAMERALENSES = "CameraLenses.jld2"
if !isfile(CAMERALENSES)
    @warn "file with previous CameraLens calibrations not found, will create empty one called $CAMERALENSES"
    close(JLD2.jldopen(CAMERALENSES,"w"))
end

@everywhere begin
    # Lg = Logging

    abstract type LAIresultInfo; end
    "Convenience type to hold results from LAI calculation."
    struct LAIresult <: LAIresultInfo
        imagepath::AbstractString
        LAI::Float64
        LAIe::Float64
        thresh::Float64
        clump::Float64
        overexposure::Float64
        csv_gapfraction::AbstractString
        csv_histogram::AbstractString
        csv_exif::AbstractString
        csv_stats::AbstractString
        jpgpath::AbstractString
        binpath::AbstractString
        ALIA::Float64
    end
    struct NoLAIresult <: LAIresultInfo
        exception::Exception
    end

    using LeafAreaIndex
    using Logging, Dates
    function getLAI(imagepath::AbstractString, cl::LeafAreaIndex.CameraLens,
                    slp::Union{LeafAreaIndex.SlopeParams, Missing})
        # This function gets executed in parallel, so need to set up new logger
        # on each processor.
        id = myid() # ID of current processor for logging file
        #@show "before logger"
        #baselog, logext = splitext(mainlogfile)
        #locallogfile = baselog * string(myid()) * logext
        #writecsv(locallogfile, "") #clear logfile
        #println("csv written")
        logger_local_io = open(joinpath("logs", "locallog$(id).log"), "a")
        logger_local = SimpleLogger(logger_local_io, Logging.Debug)
        with_logger(logger_local) do

        try
            #@show ("start getLAI on $imagepath")
			@debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - start processing on $imagepath with slope params $slp"
            println("start processing on $imagepath with slope params $slp")
            img = readrawjpg(imagepath, slp)
            #@show "image read"
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - image read"
            polim = LeafAreaIndex.PolarImage(img, cl, slp)
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - PolarImage created"
            thresh = LeafAreaIndex.threshold(polim)
			#@show "thresh calculated"
			@debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - threshold calculated"
            csv_gf = csv_gapfraction(polim, thresh)
			#@show "csv_gapfraction calculated"
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - created csv_gapfraction"
			csv_hist = csv_histogram(polim.img)
			@debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - created csv_histogram"
            csv_ex = csv_exif(imagepath)
			#@debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - created csv_exif"
            csv_st = csv_stats(polim)
			@debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - created csv_stats"
            jpgfn, binfn = write_bin_jpg(polim, thresh, imagepath)
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - wrote bin and jpg"
            LAIe, ALIA = LeafAreaIndex.inverse(polim, thresh)
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - effective LAI: $LAIe"
            clump = LeafAreaIndex.langxiang(polim, thresh, 0, pi/2)
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - clumping: $clump"
            LAI = LAIe / clump
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - LAI: $LAI"
            overexposure = sum(img .== 1) / (pi * cl.fθρ(pi/2)^2)
			@debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - overexposure calculated"
			res = LAIresult(imagepath, LAI, LAIe, thresh, clump, overexposure, csv_gf, csv_hist, csv_ex, csv_st, jpgfn, binfn, ALIA)
			@debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - LAIresult created"
            flush(logger_local_io)
            close(logger_local_io)
            return res
        catch lai_err
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - error: $lai_err"
            flush(logger_local_io)
            close(logger_local_io)
            return NoLAIresult(lai_err)
        end
        end #with_logger

        #close(logger_local_io)
    end

    function readrawjpg(imp::AbstractString, slp::Union{LeafAreaIndex.SlopeParams, Missing})
        #i = myid()# ID of current processor
        #@debug "$i start reading $imp")
        @assert imp != nothing
        @assert isfile(imp)

        ext = lowercase(splitext(imp)[end])

        if ext in LeafAreaIndex.RAW_EXT
            #imgblue = LeafAreaIndex.rawblueread(imp)
            img = LeafAreaIndex.rawblueread(imp)
			imgblue = Images.blue.(img)
        elseif ext in [".jpg",".jpeg", ".tiff"]
            img = FileIO.load(imp)
            imgblue = Images.blue.(img)
            gamma_decode!(imgblue)
        else
            @error("image has unknown extension at $imp")
        end
        @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - image read"

        #@show "check overexposure"
        if sum(imgblue .== 1) > 0.005 * length(imgblue)
            @warn("Image overexposed: $imp")
            #warn(setlog, "$i Image overexposed: $imp")
        end

        #rotate if in portrait mode
        if size(imgblue,1) > size(imgblue,2)
            ismissing(slp) || @error("image with slope in portrait mode, don't know which way to turn: $imp")
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - will rotate image in portrait mode" imp
            imgblue = rotate90(imgblue) #default clockwise, could influence result due to lens center
        end
        return imgblue
    end

    # Rotate sometimes because currently LeafAreaIndex expects landscape *in memory*.
    "Rotates an image (or in general an `AbstractMatrix`) 90 degrees."
    function rotate90(img; clockwise=true)
        transf = recenter(RotMatrix(ifelse(clockwise, 1, -1)*pi/2), center(img))
        img = warp(img, transf)
        #fix for images.jl #717
        return parent(img)
    end

    "Gamma decode a gray image taken from single channel in sRGB colorspace."
    function gamma_decode!(A::AbstractMatrix)
        @fastmath for i in eachindex(A)
            # See https://en.wikipedia.org/wiki/SRGB
            A[i] = A[i] <= 0.04045 ? A[i]/12.92 : ((A[i]+0.055)/1.055)^2.4
        end
    end

using Images
    function write_bin_jpg(polim::LeafAreaIndex.PolarImage, thresh, imgfilepath)
        jpgfilepath, binfilepath = make_bin_jpg_paths(imgfilepath)
        left, right, down, up = cropbox(polim)
        image = polim.img[down:up, left:right]
        Images.save(jpgfilepath, image)
        imgage_gray = Images.Gray.(image .> thresh)
        Images.save(binfilepath, imgage_gray)
        return jpgfilepath, binfilepath
    end
    function cropbox(polim::LeafAreaIndex.PolarImage)
        ci = polim.cl.params.lenscenter[1]
        cj = polim.cl.params.lenscenter[2]
        radius = floor(Int, polim.cl.fθρ(pi/2))
        # FIXME use CameraLensParams
        left, right = cj - radius, cj + radius
        down, up    = ci - radius, ci + radius
        # prevent out of bounds
        left, right  = max(1, left), min(size(polim.cl)[2], right)
        down, up     = max(1, down), min(size(polim.cl)[1], up)
        return left, right, down, up
    end
    function make_bin_jpg_paths(imgfilepath)
        imgdir, imgfile = splitdir(imgfilepath)
        imgbase, imgext = splitext(imgfile)
        bindir = joinpath(imgdir, "bin")
        isdir(bindir) || mkpath(bindir)
        test =Dates.format(Dates.now(), "_yyyymmdd-HHhMM")
        binfilepath = joinpath(bindir, imgbase*"$(test)_bin.png")
        jpgdir = joinpath(imgdir, "jpg")
        isdir(jpgdir) || mkpath(jpgdir)
        jpgfilepath = joinpath(jpgdir, imgbase*"$(test)_jpg.jpg")
        return jpgfilepath, binfilepath
    end

    function csv_gapfraction(polim::LeafAreaIndex.PolarImage, thresh)
        Nrings = LeafAreaIndex.Nrings_def(polim)
        θmax = LeafAreaIndex.maxviewangle(polim)
        θedges, θmid, K = LeafAreaIndex.contactfreqs(polim, 0.0, θmax, Nrings, thresh)
        T = exp.(-K./cos.(θmid))
        csv = "view_angle, gapfraction\n "
        for i in 1:length(T)
            csv *= "$(θmid[i]), $(T[i])\n "
        end
        return csv
    end

    function csv_histogram(img, bins=256)
        hist_range = -1 / (bins-1) : 1 / (bins-1) : 1
        log10nz(x) = x == 0 ? 0 : log10(x)
        counts = log10nz.(LeafAreaIndex.fasthist(reshape(img, length(img)), hist_range))
        csv = "intensity, count_log10\n "
        for i in 1:length(counts)
            csv *= "$(i/bins), $(counts[i])\n "
        end
        return csv
    end

    #PyCall.py"""
    #py = PyCall.pyimport("py")
    #py"""
    #import exifread

    #def exif(path):
    #    tags = {}
    #    with open(path, 'rb') as f:
    #        tags = exifread.process_file(f, details=False)
    #    kinds = ['EXIF ExposureTime', 'EXIF FNumber', 'EXIF ISOSpeedRatings', 'Image Orientation',
    #            'Image DateTime', "Image Make", "Image Model", "EXIF FocalLength" ]
        # for kind in kinds:
        #     print(str(tags.get(kind)))
    #    res = {k:str(v) for (k,v) in tags.items() if k in kinds}
    #    return res
    #"""
    function csv_exif(imgfilepath)
    #    tags = py"exif"(imgfilepath)
    #    csv = "key, value\n "
    csv = "todo"
    #    for (k,v) in tags
    #        csv *= "$k, $v\n "
    #    end
        return csv
    end
    function csv_stats(polim::LeafAreaIndex.PolarImage)
        left, right, down, up = cropbox(polim)
        image = polim.img[down:up, left:right]
        len = length(image)
        image = reshape(image, len)
        sort!(image)
        csv = "percentile, intensity\n "
        for p in [0.95, 0.98, 0.99, 0.999]
            csv *= "$(p), $(float(image[floor(Int, len*p)])) \n "
        end
        csv *= "max, $(float(image[end])) \n "
        return csv
    end
end

function processcenterfile(dfcenter, height, width, logfile)
    writecsv(logfile, "") #clear logfile
    # setlog = Logger("setlog")
    # Logging.configure(setlog, filename=logfile, level=DEBUG)

    # debug(setlog, "Start calibrate center ")
    calres = calibrate_center(dfcenter, height, width)
    # debug(setlog, "calibration result: $calres")
    return (calres)
end

function processprojfile(dfproj, height, width, logfile)
    # writecsv(logfile, "") #clear logfile
    # setlog = Logger("setlog")
    # Logging.configure(setlog, filename=logfile, level=DEBUG)

    # debug(setlog, "Start calibrate projection ")
    calres = calibrate_projfun(dfproj, height, width)
    # debug(setlog, "calibration result: $calres")
    return (calres)
end

function processimages(images, lensparams, logfile, datafile)
    N = nrow(images)

    ## LOGGING
    # Create specific logger per set with debug info
    #writecsv(logfile, "") #clear logfile
    logger_set_io = open(logfile, "w+")
    logger_set = SimpleLogger(logger_set_io, Logging.Debug)
    with_logger(logger_set) do

    println("Start `processimages` with lens parameters $lensparams" )
    # and slope parameters $slopeparams")
    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - Start `processimages` with lens parameters $lensparams"# and slope parameters $slopeparams"
    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - received $N image paths"
println(images)
    # create result dictionary
    result = Dict{String, Any}("success" => false)

    slope = images[1,2]
    slopeaspect = images[1,3]
    if (slope == nothing || slope == zero(slope) )
        myslopeparams = missing
    else
        myslopeparams= SlopeParams(slope/180*pi, slopeaspect/180*pi)
    end

    # load first image for image size, required for calibration
    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - load first image for image size from $(images[1, :path])"
    imgsize = size(readrawjpg(images[1, 1], myslopeparams))

    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - calibrate CameraLens or load previous calibration"
    mycamlens = load_or_create_CameraLens(imgsize, lensparams, logfile)

    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - parallel process getLAI"
    #needed for anon functions in CameraLens
    sendto(procs(), lensparams=lensparams, mycamlens=mycamlens, myslopeparams=myslopeparams)
    @everywhere lensx, lensy, lensa, lensb, lensρ = lensparams
    #remotecall_fetch(2, println, mycamlens)

    resultset = []
    for x in eachrow(images)
        path = x.path
        slope = x.slope
        slopeaspect = x.slopeAspect
        if (slope == nothing || slope == zero(slope) )
            myslopeparams = missing
        else
            myslopeparams= SlopeParams(slope/180*pi, slopeaspect/180*pi)
        end
        res = getLAI(path, mycamlens, myslopeparams)
        push!(resultset,res)
    end
    @debug  "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - parallel process done"

    # Create datafile with calculated values
    datalog = open(datafile, "w")
    truncate(datalog, 0)
    close(datalog)
    datalog = open(datafile, "a+")
    write(datalog, "Filename, LAI, LAIe, Threshold_RC, Clumping_LX, Overexposure\n")
    witherror = false
    result["csv_gapfraction"] = Dict{String, String}()
	result["csv_histogram"] = Dict{String, String}()
	result["csv_exif"] = Dict{String, String}()
    result["csv_stats"] = Dict{String, String}()
    result["jpgpath"] = Dict{String, String}()
    result["binpath"] = Dict{String, String}()
    result["LAIs"] = Dict{String, Float64}()
    result["LAIe"] = Dict{String, Float64}()
    result["threshold"] = Dict{String, Float64}()
    result["clumping"] = Dict{String, Float64}()
    result["overexposure"] = Dict{String, Float64}()
    result["ALIA"] = Dict{String, Float64}()
    for lai in resultset
        if !isa(lai, LAIresult)
            witherror = true
            @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - found error in LAIresult $lai"
            continue
        end
        #overexp_str = @sprintf("%.7f", lai.overexposure)
        overexp_str = lai.overexposure
        write(datalog, "$(basename(lai.imagepath)), $(lai.LAI), $(lai.LAIe), $(lai.thresh), $(lai.clump), $(overexp_str)\n")
        result["csv_gapfraction"][lai.imagepath] = lai.csv_gapfraction
        result["csv_histogram"][lai.imagepath] = lai.csv_histogram
        result["csv_exif"][lai.imagepath] = lai.csv_exif
        result["csv_stats"][lai.imagepath] = lai.csv_stats
        result["jpgpath"][lai.imagepath] = lai.jpgpath
        result["binpath"][lai.imagepath] = lai.binpath
        result["LAIs"][lai.imagepath] = lai.LAI
        result["LAIe"][lai.imagepath] = lai.LAIe
        result["threshold"][lai.imagepath] = lai.thresh
        result["clumping"][lai.imagepath] = lai.clump
        result["overexposure"][lai.imagepath] = overexp_str
        result["ALIA"][lai.imagepath] = lai.ALIA
    end
    close(datalog)
    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - closed $datafile"
    flush(logger_set_io)
    witherror && (return result)

    LAIs = Float64[r.LAI for r in resultset]
    result["LAI"] = median(LAIs)
    result["LAIsd"] = StatsBase.mad(LAIs,normalize=false)
    result["success"] = true
    #end #with_logger

    result
end
end

function load_or_create_CameraLens(imgsize, lensparams, setlog)
    @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - start load_or_create_CameraLens"
    @assert isfile(CAMERALENSES)

    lenshash = string(hash( (imgsize,lensparams) ))  #create unique cameralens identifier
    past_hashes = JLD2.jldopen(CAMERALENSES, "r") do file
         keys(file)
    end

    if lenshash in past_hashes
        @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - previous calibration found for hash $lenshash"
        mycamlens = FileIO.load(CAMERALENSES, lenshash)
    else
        lensx, lensy, lensa, lensb, lensρ = lensparams
        # Generic functions can't serialize, so need anonymous function to save
        projfθρ = θ -> (lensa*θ + lensb*θ^2) * lensρ
        invprojfρθ = ρ ->(-lensa + sqrt(lensa^2+4lensb*ρ/lensρ)) / 2lensb
        lensb == zero(lensb) && (invprojfρθ = ρ -> ρ / (lensρ * lensa))
        @assert projfθρ(pi/2) > 2
        @assert projfθρ(pi/2) < maximum(imgsize)

        # Fix likely lens coordinates mistake.
        if lensx > lensy # rowcoord > colcoord
            lensx, lensy = lensy, lensx
            @warn("lensx > lensy, probably a mistake, values have been swapped.")
            @warn(setlog,"lensx > lensy, probably a mistake, values have been swapped.")
        end

        @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - calibrate new mycamlens"
        #mycamlens = CameraLens(imgsize...,lensx,lensy,projfθρ,invprojfρθ)
        mycamlens = CameraLens(imgsize,(lensx,lensy), lensρ, [lensa, lensb])
        @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - calibrated new mycamlens, now save to file"
        JLD2.jldopen(CAMERALENSES, "r+") do file #"r+" to append writing data
            file[lenshash] = mycamlens
        end
        @debug "$(Dates.format(Dates.now(), "dd u yyyy HH:MM:SS")) - new mycamlens saved to file: $lenshash"
    end
    mycamlens
end

max_stars_repo_name: ETC-UA/LAIscript

max_issues_repo_name: ETC-UA/LAIscript

max_forks_repo_name: ETC-UA/LAIscript


text:
export Native
module Native
using ..Mocha
using Compat

const libname = "libmochaext.so"
const libpath = abspath(joinpath(dirname(@__FILE__), "..", "deps", libname))
try
  println("Loading native extension $libname...")
  const global library = Libdl.dlopen(libpath)
catch y
  println("ERROR: Could not load native extension at $libpath.")
  println("To use native extension, run deps/build.jl to compile the native code.")
  throw(y)
end
println("Native extension loaded")

end # module Native


text:
"""
Provide an API for the interface defined in StatsBase.

cov(ce::CovarianceEstimator, X::AbstractMatrix, [w::AbstractWeights]; mean=nothing, dims::Int=1)

Compute the covariance matrix of the matrix `X` along dimension `dims`
using estimator `ce`. A weighting vector `w` can be specified.

The keyword argument `mean` can be:
* `nothing` (default) in which case the mean is estimated and subtracted
  from the data `X`,
* a precalculated mean in which case it is subtracted from the data `X`.
  Assuming `size(X)` is `(N,M)`, `mean` can either be:
  * when `dims=1`, an `AbstractMatrix` of size `(1,M)`,
  * when `dims=2`, an `AbstractVector` of length `N` or an `AbstractMatrix`
    of size `(N,1)`.
"""

# StatsBase.cov(ce::RobustVariance, X::Vector, args...; kwargs...) =
#         StatsBase.cov(ce, reshape(X, (length(X),1)), args...; kwargs...)

# function StatsBase.cov(ce::RobustVariance, X::AbstractMatrix; mean=nothing, dims::Int=1, corrected::Bool=false)
#     dims ∈ (1, 2) || throw(ArgumentError("Argument dims can only be 1 or 2 (given: $dims)"))
#     Z = dims == 1 ? X : X'
#     if mean === nothing
#         covariance(Z, ce, Matrix, scale=size(X, dims)-corrected, demean=true)
#     else
#         @assert length(mean) == size(X, dims == 1 ? 2 : 1)
#         covariance(Z .- mean, ce, Matrix, scale=size(Z, dims)-corrected, demean=false)
#     end
# end

# function StatsBase.cov(ce::RobustVariance, X::AbstractMatrix, cache::AbstractCache; mean=nothing, dims::Int=1, corrected::Bool=false)
#     if mean === nothing
#         covariance(X, ce, cache, Matrix, SVD, scale = size(X, dims)-corrected, demean = true)
#     else
#         @assert length(mean) == size(X, dims == 1 ? 2 : 1)
#         covariance(X .- mean, ce, Matrix, SVD, scale = size(X, dims)-corrected, demean = false)
#     end
# end


text:
macro _remove_fn(objname, model_type, idx_type, args...)
    body = last(args)
    typeof(body) == Expr || throw(DomainError(body, "missing function body"))
    plural = :plural in args
    plural_s = plural ? "s" : ""
    inplace = :inplace in args
    fname = Symbol(:remove_, objname, plural_s, inplace ? "!" : "")
    idx_var = Symbol(
        objname,
        idx_type == :Int ? "_idx" :
        idx_type == :String ? "_id" :
        throw(DomainError(idx_type, "unsupported index type for _remove_fn macro")),
        plural_s,
    )

    if plural
        idx_type = AbstractVector{eval(idx_type)}
    end

    docstring = """
        $fname(model::$model_type, $idx_var::$idx_type)

    Remove $objname$plural_s from the model of type `$model_type`
    $(inplace ? "in-place" : "and return the modified model").
    """

    return :(@doc $docstring $fname(model::$model_type, $idx_var::$idx_type) = $body)
end


text:
module Formulas

export @~,
       Term,
       Terms,
       Formula

# package code goes here
include("formula.jl")
include("terms.jl")

end # module


text:
"""
**FUNÇÃO VOLUME CILINDRO**
Sintaxe:
volcilindro(raio, altura_cilindro)
"""
function volcilindro(raio, altura)
    area_base = areaBC(raio)
    println("Volume do Cilindro")
    return  area_base.*altura
end


text:
#### Restrict/prolong operators ####

# These are operators for solving PDEs via multigrid. These are also
# very fast routines for changing the size of an array by twofold,
# using linear interpolation for enlarging, and properly anti-aliasing
# for shrinking.

typealias LapackScalar Union(Float32, Float64, Complex64, Complex128)
# prolong for BLAS types
function prolong{T<:LapackScalar}(A::Array{T}, dim::Integer, len::Integer)
    sz = [size(A)...]
    sz[dim] = prolong_size(sz[dim], len)
    P = zeros(eltype(A), to_tuple(sz))
    sA = [strides(A)...]
    skipA = sA[dim]
    sP = [strides(P)...]
    skipP = sP[dim]
    n = size(A, dim)
    sz[dim] = 1
    if isodd(len)
        # Copy on-grid, and interpolate at half-grid points
        for indices = Counter(sz)
            startA = coords2lin(indices, sA)
            startP = coords2lin(indices, sP)
            # handle the on-grid points
            copy!(P, range(startP, 2*skipP, n), A, range(startA, skipA, n))
            # handle the off-grid points (linear interpolation)
            rP = range(startP+skipP, 2*skipP, n-1)
            Base.LinAlg.BLAS.axpy!(0.5, A, range(startA, skipA, n-1), P, rP)
            Base.LinAlg.BLAS.axpy!(0.5, A, range(startA+skipA, skipA, n-1), P, rP)
        end
    else
        # Interpolate at 1/4 and 3/4 points
        for indices = Counter(sz)
            startA = coords2lin(indices, sA)
            startP = coords2lin(indices, sP)
            rA1 = range(startA, skipA, n-1)
            rA2 = range(startA+skipA, skipA, n-1)
            rP1 = range(startP, 2*skipP, n-1)
            rP2 = range(startP+skipP, 2*skipP, n-1)
            Base.LinAlg.BLAS.axpy!(0.75, A, rA1, P, rP1)
            Base.LinAlg.BLAS.axpy!(0.25, A, rA2, P, rP1)
            Base.LinAlg.BLAS.axpy!(0.25, A, rA1, P, rP2)
            Base.LinAlg.BLAS.axpy!(0.75, A, rA2, P, rP2)
        end
    end
    return P
end
# prolong for non-BLAS types (including tiled arrays)
function prolong{T}(A::Array{T}, dim::Integer, len::Integer)
    sz = [size(A)...]
    sz[dim] = prolong_size(sz[dim], len)
    sA = [strides(A)...]
    skipA = sA[dim]
    P = Array(typeof(0.5*A[1]+0.5*A[1+skipA]), to_tuple(sz))
    sP = [strides(P)...]
    skipP = sP[dim]
    n = size(A, dim)
    sz[dim] = 1   # prolonged dimension is handled directly
    if isodd(len)
        for indices = Counter(sz)
            iA = coords2lin(indices, sA)
            iP = coords2lin(indices, sP)
            for i = 1:n-1
                P[iP] = A[iA]
                P[iP+skipP] = lincomb(0.5, A[iA], 0.5, A[iA+skipA])
                iA += skipA
                iP += 2*skipP
            end
            P[iP] = A[iA]
        end
    else
        for indices = Counter(sz)
            iA = coords2lin(indices, sA)
            iP = coords2lin(indices, sP)
            for i = 1:n-1
                P[iP] = lincomb(0.75, A[iA], 0.25, A[iA+skipA])
                P[iP+skipP] = lincomb(0.25, A[iA], 0.75, A[iA+skipA])
                iA += skipA
                iP += 2*skipP
            end
        end
    end
    return P
end

function prolong(A::Array, sz::Array{Int})
    if ndims(A) != size(sz,1)
        error("Array of sizes must have as many rows as dimensions of A")
    end
    Ap = A
    for j = 1:size(sz, 2)
        for i = 1:size(sz, 1)
            if sz[i,j] > size(Ap, i)
                Ap = prolong(Ap, i, sz[i,j])
            end
        end
    end
    return Ap
end

# restrict for BLAS types
function restrict{T<:LapackScalar}(A::Array{T}, dim::Integer, scale::Real)
    sz = [size(A)...]
    sz[dim] = restrict_size(sz[dim])
    R = zeros(eltype(A), to_tuple(sz))
    sA = [strides(A)...]
    skipA = sA[dim]
    sR = [strides(R)...]
    skipR = sR[dim]
    n = sz[dim]
    sz[dim] = 1
    if isodd(size(A, dim))
        for indices = Counter(sz)
            startA = coords2lin(indices, sA)
            startR = coords2lin(indices, sR)
            Base.LinAlg.BLAS.axpy!(scale, A, range(startA, 2*skipA, n), R, range(startR, skipR, n))
            rA = range(startA+skipA, 2*skipA, n-1)
            rR = range(startR, skipR, n-1)
            Base.LinAlg.BLAS.axpy!(scale/2, A, rA, R, rR)
            rR = range(startR+skipR, skipR, n-1)
            Base.LinAlg.BLAS.axpy!(scale/2, A, rA, R, rR)
        end
    else
        for indices = Counter(sz)
            startA = coords2lin(indices, sA)
            startR = coords2lin(indices, sR)
            rA = range(startA, 2*skipA, n-1)
            rR = range(startR, skipR, n-1)
            Base.LinAlg.BLAS.axpy!(0.75*scale, A, rA, R, rR)
            rA = range(startA+skipA, 2*skipA, n-1)
            Base.LinAlg.BLAS.axpy!(0.25*scale, A, rA, R, rR)
            rR = range(startR+skipR, skipR, n-1)
            Base.LinAlg.BLAS.axpy!(0.75*scale, A, rA, R, rR)
            rA = range(startA, 2*skipA, n-1)
            Base.LinAlg.BLAS.axpy!(0.25*scale, A, rA, R, rR)
        end
    end
    return R
end
# restrict for non-BLAS types
function restrict{T}(A::Array{T}, dim::Integer, scale::Real)
    sz = [size(A)...]
    sz[dim] = restrict_size(sz[dim])
    sA = [strides(A)...]
    skipA = sA[dim]
    R = Array(typeof(0.5*A[1]+0.5*A[1+skipA]), to_tuple(sz))
    sR = [strides(R)...]
    skipR = sR[dim]
    n = sz[dim]
    sz[dim] = 1
    if isodd(size(A, dim))
        for indices = Counter(sz)
            iA = coords2lin(indices, sA)
            iR = coords2lin(indices, sR)
            R[iR] = lincomb(scale, A[iA], scale/2, A[iA+skipA])
            for i = 2:n-1
                iA += 2*skipA
                iR += skipR
                R[iR] = lincomb(scale/2, A[iA-skipA], scale, A[iA], scale/2, A[iA+skipA])
            end
            R[iR+skipR] = lincomb(scale/2, A[iA+skipA], scale, A[iA+2*skipA])
        end
    else
        for indices = Counter(sz)
            iA = coords2lin(indices, sA)
            iR = coords2lin(indices, sR)
            R[iR] = lincomb(0.75*scale, A[iA], 0.25*scale, A[iA+skipA])
            iR += skipR
            for i = 2:n-1
                R[iR] = lincomb(0.25*scale, A[iA], 0.75*scale, A[iA+skipA], 0.75*scale, A[iA+2*skipA], 0.25*scale, A[iA+3*skipA])
                iR += skipR
                iA += 2*skipA
            end
            R[iR] = lincomb(0.25*scale, A[iA], 0.75*scale, A[iA+skipA])
        end
    end
    return R
end
restrict{T}(A::Array{T}, dim::Integer) = restrict(A, dim, one(T))
function restrict(A::Array, flag::Union(Array{Bool},BitArray), scale::Real)
    func = (A, dim) -> restrict(A, dim, scale)
    return mapdim(func, A, flag)
end
restrict{T}(A::Array{T}, flag::Union(Array{Bool},BitArray)) = restrict(A, flag, one(eltype(T)))

# "Balanced" restriction and prolongation

# If you take an array of ones, restrict it, and then prolong it
# (i.e., compute P*R*A, where P is the prolongation operator and R is
# the restriction operator), it smooths the array. However, by
# comparison to all other points in the array the edges will tend
# towards zero. The "balanced" versions of restrict/prolong ensure
# that each point of the array is used equally, i.e., that sum(P*R, 1)
# is constant, and hence the edges will be preserved.
function restrictb{T}(A::Array{T}, dim::Integer, scale::Real)
    a::T
    b::T
    Ar = restrict(A, dim, scale)
    sz = [size(Ar)...]
    sz[dim] = 1
    s = [strides(A)...]
    sr = [strides(Ar)...]
    offset_end = (size(A,dim)-1)*s[dim]
    offset_endr = (size(Ar,dim)-1)*sr[dim]
    Astep = s[dim]
    if isodd(size(A,dim))
        b = scale/sqrt(3)
        a = 2*b
        for indices = Counter(sz)
            index = coords2lin(indices, s)
            indexr = coords2lin(indices, sr)
            Ar[indexr] = a*A[index] + b*A[index+Astep]
            Ar[indexr+offset_endr] = a*A[index+offset_end] + b*A[index+offset_end-Astep]
        end
    else
        b = scale/(2*sqrt(2))
        a = 3*b
        for indices = Counter(sz)
            index = coords2lin(indices, s)
            indexr = coords2lin(indices, sr)
            Ar[indexr] = a*A[index] + b*A[index+Astep]
            Ar[indexr+offset_endr] = a*A[index+offset_end] + b*A[index+offset_end-Astep]
        end
    end
    return Ar
end
restrictb{T}(A::Array{T}, dim::Integer) = restrictb(A, dim, one(T))
function restrictb(A::Array, flag::Union(Array{Bool},BitArray), scale::Real)
    func = (A, dim) -> restrictb(A, dim, scale)
    return mapdim(func, A, flag)
end
restrictb{T}(A::Array{T}, flag::Union(Array{Bool},BitArray)) = restrictb(A, flag, one(eltype(T)))
function prolongb{T}(A::Array{T}, dim::Integer, len::Integer)
    a::T
    b::T
    Ap = prolong(A, dim, len)
    sz = [size(Ap)...]
    sz[dim] = 1
    s = [strides(A)...]
    sp = [strides(Ap)...]
    offset_end = (size(A,dim)-1)*s[dim]
    offset_endp = (size(Ap,dim)-1)*sp[dim]
    Astep = s[dim]
    Apstep = sp[dim]
    if isodd(len)
        b = 1/sqrt(3)
        a = 2*b
        for indices = Counter(sz)
            index = coords2lin(indices, s)
            indexp = coords2lin(indices, sp)
            Ap[indexp] = a*A[index]
            Ap[indexp+Apstep] = b*A[index] + A[index+Astep]/2
            Ap[indexp+offset_endp] = a*A[index+offset_end]
            Ap[indexp+offset_endp-Apstep] = b*A[index+offset_end] + A[index+offset_end-Astep]/2
        end
    else
        b = 1/(2*sqrt(2))
        a = 3*b
        for indices = Counter(sz)
            index = coords2lin(indices, s)
            indexp = coords2lin(indices, sp)
            Ap[indexp] = a*A[index] + A[index+Astep]/4
            Ap[indexp+Apstep] = b*A[index] + (3/4)*A[index+Astep]
            Ap[indexp+offset_endp] = a*A[index+offset_end] + A[index+offset_end-Astep]/4
            Ap[indexp+offset_endp-Apstep] = b*A[index+offset_end] + (3/4)*A[index+offset_end-Astep]
        end
    end
    return Ap
end
function prolongb(A::Array, sz::Array{Int})
    if ndims(A) != size(sz,1)
        error("Array of sizes must have as many rows as dimensions of A")
    end
    Ap = A
    for j = 1:size(sz, 2)
        for i = 1:size(sz, 1)
            if sz[i,j] > size(Ap, i)
                Ap = prolongb(Ap, i, sz[i,j])
            end
        end
    end
    return Ap
end

# restrict, then "repair" the edges by linear extrapolation
# Do _not_ use this if you're doing gradient restriction/prolongation,
# as this breaks the chain rule. Just use ordinary restrict/prolong or
# their balanced "b" variants.
function restrict_extrap(A::Array, dim::Integer, scale::Real)
    Ar = restrict(A, dim, scale)
    sz = [size(Ar)...]
    sz[dim] = 1
    s = [strides(A)...]
    sr = [strides(Ar)...]
    offset_end = (size(A,dim)-1)*s[dim]
    offset_endr = (size(Ar,dim)-1)*sr[dim]
    if isodd(size(A,dim))
        for indices = Counter(sz)
            index = coords2lin(indices, s)
            indexr = coords2lin(indices, sr)
            Ar[indexr] = (2*scale)*A[index]
            Ar[indexr+offset_endr] = (2*scale)*A[index+offset_end]
        end
    else
        for indices = Counter(sz)
            index = coords2lin(indices, s)
            indexr = coords2lin(indices, sr)
            Ar[indexr] = (3*scale)*A[index] - scale*A[index+s[dim]]
            Ar[indexr+offset_endr] = (3*scale)*A[index+offset_end] - scale*A[index+offset_end-s[dim]]
        end
    end
    return Ar
end
restrict_extrap{T}(A::Array{T}, dim::Integer) = restrict_extrap(A, dim, one(T))
function restrict_extrap(A::Array, flag::Union(Array{Bool},BitArray), scale::Real)
    func = (A, dim) -> restrict_extrap(A, dim, scale)
    return mapdim(func, A, flag)
end
restrict_extrap{T}(A::Array{T}, flag::Union(Array{Bool},BitArray)) = restrict_extrap(A, flag, one(T))

function prolong_size(sz::Int, len::Int)
    if isodd(len)
        newsz = 2*sz-1
    else
        newsz = 2*(sz-1)
    end
    if newsz != len
        error("Cannot prolong a dimension of length ", sz, " to ", len)
    end
    return newsz
end

restrict_size(len::Int) = isodd(len) ? div(len+1,2) : div(len,2)+1
function restrict_size(szin::Union(Dims, Vector{Int}), flag::Union(Array{Bool},BitArray))
    if length(szin) != size(flag,1)
        error("Boolean flag must have as many rows as dimensions of A")
    end
    sz = Array(Int, size(flag))
    # First round of restriction comes from original size
    for i = 1:length(szin)
        sz[i] = flag[i] ? restrict_size(szin[i]) : szin[i]
    end
    # Remaining rounds of restriction come from the previous size
    for j = 2:size(flag, 2)
        for i = 1:size(flag, 1)
            sz[i,j] = flag[i,j] ? restrict_size(sz[i,j-1]) : sz[i,j-1]
        end
    end
    return sz
end
function prolong_size(szin::Union(Dims, Vector{Int}), flag::Union(Array{Bool},BitArray))
    sz = restrict_size(szin, flag)
    return [sz[:,end-1:-1:1] [szin...]]
end


text:
import Random
using Test: @test, @testset, @inferred
using Statistics: mean, var
using Random: bitrand, rand!, randn!
using RestrictedBoltzmannMachines: RBM, Binary, visible, hidden, weights, free_energy
using RestrictedBoltzmannMachines: Gaussian, ReLU, dReLU, pReLU, xReLU
using RestrictedBoltzmannMachines: sample_v_from_v, sample_h_from_h, mean_from_inputs, var_from_inputs
using RestrictedBoltzmannMachines: rescale_hidden!, rescale_activations!

Random.seed!(23)

@testset "rescale_activations!" begin
    N = (7,2)
    layers = (
        Gaussian(randn(N...), randn(N...)),
        ReLU(randn(N...), randn(N...)),
        dReLU(randn(N...), randn(N...), randn(N...), randn(N...)),
        pReLU(randn(N...), randn(N...), randn(N...), 2rand(N...) .- 1),
        xReLU(randn(N...), randn(N...), randn(N...), randn(N...))
    )
    λ = rand(N...)
    for layer in layers
        μ = mean_from_inputs(layer)
        ν = var_from_inputs(layer)
        rescale_activations!(layer, λ)
        @test mean_from_inputs(layer) ≈ μ ./ λ
        @test var_from_inputs(layer) ≈ ν ./ λ.^2
    end
end

@testset "rescale_hidden!" begin
    rbm = RBM(Binary(2), ReLU(1), randn(2,1))
    randn!(rbm.visible.θ)
    randn!(rbm.hidden.θ)
    rand!(rbm.hidden.γ)
    rbm.hidden.γ .+= 0.5

    v = sample_v_from_v(rbm, bitrand(size(visible(rbm))..., 20000); steps=100)
    h = sample_h_from_h(rbm, rand(size(hidden(rbm))..., 20000); steps=100)
    ave_v = mean(v; dims=2)
    ave_h = mean(h; dims=2)
    var_v = var(v; dims=2)
    var_h = var(h; dims=2)
    F = free_energy(rbm, v)

    λ = [1 + rand()]
    rescale_hidden!(rbm, λ)

    @test free_energy(rbm, v) ≈ F .+ sum(log, λ)

    v = sample_v_from_v(rbm, bitrand(size(visible(rbm))..., 20000); steps=100)
    h = sample_h_from_h(rbm, rand(size(hidden(rbm))..., 20000); steps=100)

    @test mean(v; dims=2) ≈ ave_v rtol=0.1
    @test mean(h; dims=2) ≈ ave_h ./ λ rtol=0.1
    @test var(v; dims=2) ≈ var_v rtol=0.1
    @test var(h; dims=2) ≈ var_h ./ λ.^2 rtol=0.1
end


text:
using Clang

HEADER_BASE = "/usr/include" # using system-wide installation of SCIP
all_headers = readdir(joinpath(HEADER_BASE, "scip"))
scip_headers = vcat(
    filter(h -> startswith(h, "type_"), all_headers),
    filter(h -> startswith(h, "pub_"), all_headers),
    filter(h -> startswith(h, "scip_"), all_headers),
    "scipdefplugins.h",
    filter(h -> startswith(h, "cons_"), all_headers),
    "intervalarith.h",                                 # for nlpi/pub_expr
)
lpi_headers = ["type_lpi.h"]
nlpi_headers = ["type_expr.h", "type_nlpi.h", "type_exprinterpret.h", "pub_expr.h"]

headers = vcat(
    [joinpath(HEADER_BASE, "scip", h) for h in scip_headers],
    [joinpath(HEADER_BASE, "lpi", h) for h in lpi_headers],
    [joinpath(HEADER_BASE, "nlpi", h) for h in nlpi_headers],
)

context = Clang.init(
    headers=headers,
    common_file="commons.jl",
    output_dir="../src/wrapper",
    clang_includes=vcat(HEADER_BASE, LLVM_INCLUDE),
    clang_args = ["-I", HEADER_BASE],
    clang_diagnostics=true,
    header_wrapped=(header, cursorname) -> header == cursorname,
    header_library=header_name -> "libscip"
)
Clang.run(context)


text:
"""
A context (`Ctx`) instance encapsulates the index location and provides a useful way to split indexing across multiple indices. It also holds a lock to handle it being called across tasks.
"""
struct Ctx
    store::String       # directory where to store indices
    resolver::Any       # function to determine index filename (as subpath of store) if needed
    lock::ReentrantLock # lock (`withenv` not task safe, could have used per index lock if it were)

    function Ctx(; store::String=joinpath(homedir(), ".googlecodesearchjl"), resolver=(ctx,inpath)->joinpath(ctx.store,"index"))
        mkpath(store)
        new(store, resolver, ReentrantLock())
    end
end

show(io::IO, ctx::Ctx) = print(io, "GoogleCodeSearch.Ctx(store=\"" * ctx.store * "\")")

function readcmd_with_index(ctx::Ctx, cmd::Cmd, idxpath::String)
    oc = lock(ctx.lock) do
        withenv("CSEARCHINDEX"=>idxpath) do
            OutputCollector(cmd)
        end
    end
    success = wait(oc)
    success, oc.stdout_linestream.lines, oc.stderr_linestream.lines
end

"""
Returns a list of index files that are being used in this context.
"""
indices(ctx::Ctx) = map(x->joinpath(ctx.store,x), readdir(ctx.store))

"""
Clears all indices. Fresh index will be created on next call to `index`.
"""
function clear_indices(ctx::Ctx)
    for file in indices(ctx)
        rm(file; force=true)
    end
    nothing
end

"""
Returns a list of paths that have been indexed.
"""
function paths_indexed(ctx::Ctx)
    paths = Set{String}()
    cmd = Cmd([cindex, "-list"])
    for idxpath in indices(ctx)
        success, out, err = readcmd_with_index(ctx, cmd, idxpath)
        if success
            for (t,m) in out
                push!(paths, strip(m))
            end
        else
            error("error reading index $idxpath")
        end
    end
    paths
end

"""
Index paths by calling the index method. While indexing, ensure paths are sorted such that paths appearing later are not substrings of those earlier. Otherwise, the earlier indexed entries are erased (strange behavior of `cindex`).
"""
function index(ctx::Ctx, path::String)
    cmd = Cmd([cindex, path])
    success, out, err = readcmd_with_index(ctx, cmd, ctx.resolver(ctx,path))
    success
end

function index(ctx::Ctx, paths::Vector{String})
    idxpaths = Dict{String,Vector{String}}()
    for path in paths
        idxpath = ctx.resolver(ctx,path)
        paths = get!(()->String[], idxpaths, idxpath)
        push!(paths, path)
    end
    results = Bool[]
    for (idxpath, paths) in idxpaths
        cmd = Cmd(append!([cindex], paths))
        success, out, err = readcmd_with_index(ctx, cmd, idxpath)
        push!(results, success)
    end
    results
end

"""
Search by calling the search method with a regex pattern to search for. Optionally pass the following parameters:
- `ignorecase`: boolean, whether to ignore case during search (default false)
- `pathfilter`: a regular expression string to restrict search only to matching paths

The search method returns a vector of named tuples, each describing a match.
- `file`: path that matched
- `line`: line number therein that matched
- `text`: text that matched
"""
function search(ctx::Ctx, pattern::String; ignorecase::Bool=false, pathfilter::Union{Nothing,String}=nothing, maxresults::Int=20)
    cmdparts = [csearch]
    if pathfilter !== nothing
        push!(cmdparts, "-f")
        push!(cmdparts, pathfilter)
    end
    ignorecase && push!(cmdparts, "-i")
    push!(cmdparts, "-n")
    push!(cmdparts, pattern)
    cmd = Cmd(cmdparts)
    results = Vector{NamedTuple{(:file,:line,:text),Tuple{String,Int,String}}}()
    for idx in readdir(ctx.store)
        idxpath = joinpath(ctx.store, idx)
        success, out, err = readcmd_with_index(ctx, cmd, idxpath)
        if success
            for (_,s) in out
                s = strip(s)
                ( isempty(s) || !startswith(s, "/")) && continue
                parts = split(s, ':'; limit=3)
                (length(parts) != 3) && continue
                (length(results) > maxresults) && (return results)
                try
                    push!(results, (file=String(parts[1]), line=parse(Int,parts[2]), text=String(parts[3])))
                catch ex
                    @info "Exception: ",ex
                end
            end
        end
    end
    results
end